Scalable characterization of localizable entanglement in noisy topological quantum codes

Topological quantum error correcting codes have emerged as leading candidates towards the goal of achieving large-scale fault-tolerant quantum computers. However, quantifying entanglement in these systems of large size in the presence of noise is a challenging task. In this paper, we provide two different prescriptions to characterize noisy stabilizer states, including the surface and the color codes, in terms of localizable entanglement over a subset of qubits. In one approach, we exploit appropriately constructed entanglement witness operators to estimate a witness-based lower bound of localizable entanglement, which is directly accessible in experiments. In the other recipe, we use graph states that are local unitary equivalent to the stabilizer state to determine a computable measurement-based lower bound of localizable entanglement. If used experimentally, this translates to a lower bound of localizable entanglement obtained from single-qubit measurements in specific bases to be performed on the qubits outside the subsystem of interest. Towards computing these lower bounds, we discuss in detail the methodology of obtaining a local unitary equivalent graph state from a stabilizer state, which includes a new and scalable geometric recipe as well as an algebraic method that applies to general stabilizer states of arbitrary size. Moreover, as a crucial step of the latter recipe, we develop a scalable graph-transformation algorithm that creates a link between two specific nodes in a graph using a sequence of local complementation operations. We develop open-source Python packages for these transformations, and illustrate the methodology by applying it to a noisy topological color code, and study how the witness and measurement-based lower bounds of localizable entanglement varies with the distance between the chosen qubits.


Introduction
The advancements in the science of quantum information and computation [1] has put entanglement [2] in a position of extreme importance for a number of quantum protocols including quantum teleportation [2][3][4], quantum dense coding [5][6][7], quantum cryptography [8,9], and measurement-based quantum computation [10][11][12]. Besides, entanglement has also been identified as the key ingredient in several other seemingly unrelated problems, such as the study of topological [13][14][15] and non-topological [16][17][18][19] quantum phases and corresponding quantum phase transitions in many-body systems, understanding transport properties in biological systems such as light-harvesting complexes [20][21][22][23], investigating the role of radical-pair mechanism in the navigability of animals in weak magnetic fields [24], and aspects of AdS/CFT correspondence in models of quantum gravity [25][26][27][28]. Tremendous technological development has made the laboratory realization of entangled states, in both biparty and multiparty scenario, possible by using trapped ions [29][30][31], photonic systems [32][33][34], nuclear magnetic resonance (NMR) molecules [35], superconducting qubits [36,37], cold atoms [38][39][40][41], solid-state systems [42], and set-ups involving hybrids of these systems [43,44]. This highlights the potential of realizing different quantum protocols that uses entanglement as resource in these systems. The theoretical aspect of this line of study involves quantum information processing using many-body systems realizable in the above substrates, which has particularly brought the importance of the study of quantum many-body systems in a language consistent with quantum information theory in focus. Being a key resource in quantum protocols, entanglement has been the natural choice as the characterising feature of quantum many-body systems for this purpose.
In the last decade, the extraordinary potential of quantum computation [45,46] in addressing otherwise intractable problems, such as simulating large quantum systems or decrypting codes efficiently, has led to major efforts towards fabricating scalable systems, with a long-term goal of fault-tolerance, using the available many-body substrates such as trapped ions [47] and superconducting qubits [48,49]. There already exist noisy intermediate-scale quantum (NISQ) devices [50][51][52] made of 50-100 qubits, which are being viewed and investigated as platforms to potentially achieve 'quantum supremacy' [53], and as possible candidate systems to host logical qubits as building blocks of envisioned large-scale quantum computers [54]. To date, the figures of merit for the usefulness of a quantum state, prepared in these systems, in a given quantum protocol are the different quantum correlations, such as entanglement, that serves as resource in that protocol. However, characterizing such systems using entanglement proves to be difficult mainly due to exceeding resource requirements in traditional techniques such as complete quantum state or process tomography, even for moderately sized systems of a few qubits [55,56]. Another obstacle, from the computational point of view, is the presence of noise, which requires entanglement to be computed for a mixed state of a large system-a long-standing problem of quantum information theory [2,57].
In the vision of building large-scale fault-tolerant quantum computers, topological quantum error correcting codes [58][59][60], such as, eg., the surface codes [61,62] and the color codes [63,64] (see figure 1 and appendix A for a description), are being considered as leading candidate systems. These systems are realized by arranging qubits on lattices of specific geometry, and are robust against external perturbation [65][66][67][68], qubit loss [69][70][71], as well as computational errors [72]. Attempts have recently been made to implement such systems in the laboratory with, for example, trapped ions [73][74][75], and superconducting qubits [76,77]. In order to host a single logical qubit with low or moderate code distance, and to perform error correction protocols taking into account errors on multiple physical qubits, one needs to deal with systems with a large number of physical qubits in the presence of noise. This makes the characterization of these systems using entanglement measures difficult. Therefore, quantifying bipartite as well as multipartite entanglement in subsystems of a large-scale quantum many-body system like the topological quantum codes in the presence of noise has been an active field of research in recent times [78][79][80].
There exists a plethora of avenues towards quantifying entanglement in subsystems of a quantum many-body system, which mainly follow two approaches-(1) the partial trace-based [2], and (2) the local measurement-based [81][82][83][84] approach. In the former, an entanglement measure is computed either between two subsystems denoted by Ω and Ω of the full quantum many-body system, or over a subsystem Ω of a quantum many-body system by using the reduced state ρ Ω = Tr Ω (ρ) of Ω. The reduced state is obtained by performing a partial trace operation on the rest of the system Ω, where Ω ∪ Ω represents the entire system, and Ω ∩ Ω = ∅. This protocol is particularly useful in scenarios where the state of the entire system is pure, leading to concepts like topological entanglement entropy [85][86][87][88] between the subsystems Ω and Ω. This is also effective in situations where the partial trace operation results in a reduced density matrix ρ Ω that faithfully provides a non-zero value of a chosen bipartite or multipartite entanglement measure thereby quantifying entanglement over Ω, when the subsystems constituting the region Ω are entangled accordingly. Apart from leading to the study of topological entanglement entropy as a function of the system parameters in topological quantum codes [89,90], the partial trace-based prescription has also resulted in, for example, investigations of the behaviour of entanglement over a collection of subsystems. These include pair of nearest-neighbor spins, in a large set of quantum many-body Hamiltonians, such as the transverse field Ising and the XY models, and the XXZ model in the presence of a magnetic field [16][17][18][19].
However, topological entanglement entropy fails to provide a faithful characterization of bipartite entanglement in a topological code when the state of the system is mixed-for example, when there is noise in the system [2,78,79]. Also, there exist logical states in topological error correcting codes for which the reduced state ρ Ω is classical [73], thereby providing zero entanglement as quantified by a chosen entanglement measure despite the existence of entanglement over Ω. In these situations, the appropriate protocol is the local measurement-based approach [82][83][84], which relies on performing local projection measurements on the subsystems constituting Ω in specific measurement bases in order to create entangled post-measurement states over the region Ω. This leads to non-zero average entanglement in the region Ω after the measurement, thereby appropriately quantifying the entanglement present in the subsystem Ω. The maximum average post-measurement entanglement, maximized over all possible local projection measurements on the subsystems in Ω, is referred to as the localizable entanglement [82][83][84], which is defined for both pure as well as mixed states. Apart from being a good quantifier of local entanglement in stabilizer states [54,[91][92][93], including topological quantum codes, with or without noise, localizable entanglement is crucial also in other scenarios. For example, it has been used in conceptualizing the correlation length in certain quantum many-body systems [82,83,94,95], for characterizing local entanglement in cluster-Ising [96,97] and cluster-XY models [98], and in protocols including measurement-based quantum computation [54,91,92] and entanglement percolation in quantum network [99].
Although the definition of localizable entanglement [82][83][84] is straightforward to understand, the maximization involved in the definition often makes the computation of the quantity difficult if the size of the quantum many-body system is too large. Also, in the case of mixed states, which may originate in situations where there is noise in the system-a practical scenario being topological quantum error correcting codes, the computation of the value of localizable entanglement requires calculation of a chosen entanglement measure for the states in the post-measurement ensemble over the region Ω, which may turn out to be a challenging task. These obstacles highlight the need of a recipe to determine the localizable entanglement over a region of a large-scale topological quantum error correcting code in the presence of noise, which is the aim of this paper.
In this paper, we calculate, if not the actual value, an effective lower bound of localizable entanglement in an efficient and experiment-friendly way via the use of entanglement witness operators. We illustrate our prescriptions for a region of two qubits in a color code under local uncorrelated Pauli noise. However, the recipe is fully applicable to the more general case of arbitrary stabilizer states, and it has the potential to be generalized to the case of arbitrary regions in arbitrary topological quantum error correcting stabilizer codes. The major new results reported in this paper are as follows (see also figure 2).
(a) We propose two prescriptions for computing non-trivial lower bounds of localizable entanglement over a group of qubits in noisy topological quantum codes, and establish a connection between the two seemingly different prescriptions. The complete prescriptions are given in section 2. 1. One of the prescriptions is a method based on entanglement witness operators [100][101][102][103][104][105][106][107], which exploits appropriate construction of a local witness operator [105,107] from the topological code via a projection operator-based approach and the fact that a witness-based lower bound of localizable entanglement is computable from local witness operators [93]. We propose a specific construction of the local witness operator in the case of topological color codes, and explain the method in detail in section 2.1. 2. The second approach is graph-based. It exploits graph states [91,92] obtained via local unitary transformation from the stabilizer states representing topological codes, and subsequently . Schematic representation of the objectives of this paper. This paper focuses on computing the localizable entanglement over a specific subsystem Ω of a noisy topological quantum error correcting code, via single-qubit projection measurements on the qubits in the rest of the system, designated by Ω. Two prescriptions (section 2), one witness-based (section 2.1) and the other graph-based (section 2.2), are proposed. The latter approach requires a transformation of the stabilizer state to a graph state, where a geometric recipe for building the graph from the stabilizer structure is discussed in section 3, and an algebraic method to compute the adjacency matrices of the graphs is given in appendix E. The graph-based method further requires a fail-safe recipe to create a link between any two qubits if they do not already share a link. A graph-transformation algorithm is developed and is discussed in appendix G for this purpose. The corresponding methodologies are applied to a topological color code on a square-hexagonal lattice (section 4), and the dependence of the bound on the distance between two chosen qubits is investigated.
transforming the graphs such that the chosen region Ω becomes connected in the transformed graph. The salient features of this method are given in section 2.2.
(b) The graph-based method requires obtaining the possible graph states from an aritrary stabilizer state in a topological quantum error correcting code. A geometric prescription for the transformation of the surface code to the local unitary equivalent graph state has been proposed in reference [108]. In this paper, we discuss the geometric recipe for transforming a general stabilizer state of a topological color code to a graph state (section 3). We also present, in connection to the geometric recipe, an algebraic methodology that transforms an arbitrary stabilizer state to a graph state, which is convenient for obtaining the adjacency matrix of the graph from the structures of the stabilizer operators. This part is fairly technical, and to make the main text of the paper comprehensive and appropriate for a broad readership, we discuss the technicalities of this conversion in appendix E. The method has also been given the form of an open source Python package. (c) Moreover, for the graph-based approach, we develop a fail-safe adaptive algorithm for connecting a region of two qubits in arbitrary connected graphs by local complementation operations, and implement the algorithm in the form of an open-source Python package. This involves technical details on graphs and their transformations under specific operations. Similar to the case of the algebraic approach for obtaining the adjacency matrices of graphs, the technical details of the methodology are discussed in appendix G.
The detailed calculations of the decomposition of local entanglement witnesses, and a few examples and the necessary information on the topological color codes, graph states, and noise models have been included in the appendices A-L. To demonstrate how the newly developed methodology can be used for topological quantum codes, we apply them to the specific example of the two-dimensional color code on the square hexagonal lattice, and discuss the results regarding the noise as well as distance dependence of localizable entanglement between a pair of qubits. The results are discussed in detail in section 4. Section 5 contains the concluding remarks.

Bounds of localizable entanglement
Localizable entanglement (LE) [82][83][84]93] over a region Ω composed of the N − m qubits in an N-qubit state ρ is the maximum average entanglement that can be localized over Ω, by performing local projection measurements on the m qubits in Ω-the region outside the chosen subsystem Ω. Without any loss in generality, we assign the first m qubits in the set of N qubits labelled as 1, 2, . . . , m, m + 1, . . . , N to the region Ω and the rest in Ω, with Ω ∪ Ω representing the complete N-partite system and Ω ∩ Ω = ∅, (1 m N − 2). Considering only rank-1 projection measurements on the qubits in Ω, LE over the region Ω is given by where ρ k Ω = Tr Ω p −1 k M k ρM k is the normalized post-measurement state over the region Ω, p k = Tr[M k ρM k ] is the probability of obtaining the measurement outcome k with 2 m −1 k=0 p k = 1, and M k = [⊗ i∈Ω |k i k i |] ⊗ I Ω forms the complete set M of rank-1 projection measurements over Ω. Here |k i ∈ {|0 i , |1 i } are two arbitrary and mutually orthogonal single-qubit states in the Bloch sphere, with , |1 } being the computational basis, and {θ i , φ i } are the real azimuthal and polar angles of the Bloch sphere (0 θ i π, 0 φ i < 2π). The index k here can be identified as the multi-index k ≡ k 1 k 2 . . . k m , with k i = 0, 1. The optimization in equation (1), therefore, reduces to an optimization over 2m real parameters, and is difficult to perform when m is a large number. However, even in cases where m is small, analytical determination of LE is possible only in the cases of a very limited number of multiqubit quantum states with special properties.
In order to extract useful information about the properties of LE in situations where computing the exact value of LE proves difficult, a possible approach is to determine computable lower bounds of LE, which can provide insight of the behaviour of the actual quantity. In this spirit, one can perform the optimization over a restricted subset of the complete set of local projection measurements, for example, by allowing only local Pauli measurements over the qubits, thereby computing the restricted LE (RLE), E R Ω (ρ), of ρ [93]. However, computation of RLE is also non-trivial in the case of large N, where one has to consider a total of 3 m measurement settings, corresponding to three possible local Pauli measurements on each qubit, which is a large number when m is large. A further lower bound can be found by considering a specific Pauli measurement setting where each qubit i in Ω is measured in the basis of a specific Pauli operator σ i , where σ i = X, Y, Z. The average entanglement in the region Ω, corresponding to the chosen measurement setting P, is given by where the superscript P represents an appropriately chosen specific Pauli measurement setting, and E Ω (ρ) E R Ω (ρ) E P Ω (ρ) by the definition of LE. The challenge now is to choose an appropriate Pauli measurement setting P, which would provide a non-trivial lower bound of LE. In [93], two avenues for constructing such lower bounds of the LE have been discussed-(1) an experimentally accessible entanglement witness-based lower bound (WLB) by using entanglement witness operators appropriate for the post-measurement states on the region Ω, and (2) a graph-based approach to determine an appropriate Pauli measurement set-up P, which provides a non-trivial value of E P Ω (ρ), denoted by the measurement-based lower bound (MLB). We consider the state ρ to be a mixed one in general, originated from, for example, a stabilizer state ρ 0 due to application of noise ρ 0 → ρ = Λ(ρ 0 ), Λ(.) representing the noise model. In the following subsections, we present the underlying intricacies of computing the WLB and the MLB for LE in the case of arbitrary stabilizer states under noise. We use topological color codes [63,64] (see also appendix A for definitions) for demonstration.

Witness-based lower bound
We begin our discussion with the WLB of LE, which employs an entanglement witness operator [100][101][102][103][104][105][106][107], W, that indicates non-zero entanglement content in a quantum state ρ via a negative expectation value, i.e., Tr(Wρ) < 0. A lower bound of the entanglement content in a quantum state ρ can be determined using the set of expectation values of appropriately chosen entanglement witness operators as the solution of the optimization problem [93,[109][110][111][112][113] E min (ω) = inf E(ρ), subject to Tr (ρW) = ω, ρ 0, and Tr (ρ) = 1, where E is the chosen entanglement measure. This can be used in equation (2) to write The stabilizers S z and S x are obtained by multiplying respectively the z-and x-type stabilizers corresponding to the plaquettes on two adjacent paths of plaquettes, connecting the two qubits a and b. In the present example, the two stabilizers S j=1 and S j=2 for constructing the local witness operator W Ω are constituted of the x-and z-type plaquette stabilizers as (a) ( being the plaquette index. The distance between a and b is the length of the path constituted of the lattice links common to the two adjacent paths (shown by double continuous lines).
where ω k is the expectation value of an appropriately chosen witness operator W k Ω for the post-measurement state ρ k Ω over the region Ω such that E min (ω k ) E(ρ k Ω ), and E W Ω (ρ) is the WLB. Note that to obtain ω k , one has to (1) perform local Pauli measurements on the qubits in the region Ω, (2) choose appropriate witness operators W k Ω for the post-measurement states ρ k Ω subject to the specific measurement outcome k for the chosen Pauli measurement setting, and then (3) measure the expectation values ω k = Tr W k Ω ρ k Ω in the post-measurements states. However, this protocol can prove difficult to carry out in experiments when N and m are large numbers. An alternative approach would be to look for a witness operator W Ω , called the local witness operator, such that the expectation value of W Ω , when determined with respect to the state ρ, detects entanglement in the region Ω [105,107], and a functional relation between Tr(W Ω ρ) and E W Ω (ρ) exists. The challenge, however, in this approach is choosing an appropriate form of W Ω which can be connected to the specific measurement outcomes of the chosen Pauli measurement setting P, such that the witnesses W k Ω corresponding to different values of k can be constructed out of W Ω .
For this purpose, we construct a local witness operator in terms of stabilizers describing the stabilizer state. A local witness operator W Ω that can detect genuine multiparty entanglement in the region Ω in an arbitrary stabilizer state |ψ S , or a state ρ S close to it can be chosen to be of the form [107] where S = {S j } is a subset of the complete set of stabilizers defining the state | ψ S , given by S j = ⊗ i τ u i,j , u i,j = 0, 1, 2, 3, i is the qubit-index, j indicates which stabilizer the qubit belongs to, and τ u i,j = I, X, Y, and Z, for u i,j = 0, 1, 2, and 3 respectively (this is the same definitions of stabilizers as given in appendix A, with a new variable u i,j introduced in order to conveniently represent the Pauli matrices, which will be clear in subsequent discussions). One can write the stabilizers S j constructing W Ω by distinguishing the supporting qubits according to whether they belong to the region Ω or Ω, as where S Ω j = ⊗ l∈Ω τ Ω u l,j is the part of the stabilizer S j with support on Ω. For W Ω to detect entanglement in the region Ω, the stabilizers S j in equation (5) have to be such that [107] (a) τ Ω u i,j , τ Ω u i ,j = 0 ∀ qubits i, i , and ∀ stabilizer pairs j, j , i.e., the supports of the stabilizers involved in constructing W Ω must commute outside the region Ω, and (b) the set {S Ω j } obtained from the subset S of stabilizers S j of | ψ S is a complete set of stabilizer generators of a genuinely multipartite entangled state |ψ Ω over Ω.
See figures 3(a) and (b) for demonstrations.
The constructed local witness operator can be decomposed as a sum of the products of Pauli projection operators on Ω and the witness operators W k Ω , as (see appendix B for the derivation) with W k . v m are multi-indices with v i = 1, 2, 3, such that k represents the outcome of the projection operation, τ 1 = X, τ 2 = Y, τ 3 = Z, I Ω is the identity operator in the Hilbert space of the qubits in Ω, and η j = ±1. For a specific set of Pauli operators chosen for the qubits in Ω (i.e., for a specific value of v), the expectation value of W Ω with respect to the state ρ S is where p k = Tr P Ω (k,v) ⊗ I Ω ρ S is the probability of obtaining the measurement outcome k if the projection operator P Ω (k,v) is applied on the qubits in Ω, and to obtain this, we have used with ρ k Ω = Tr Ω P Ω (k,v) ρ S . Note that we have now been able to write the expectation value of the witness operator as the sum of the product of a number of expectation values of witness operators confined to the region Ω and the corresponding probabilities of obtaining those witness operators in Ω via a projection operation on qubits in Ω (equation (8)). This has direct resemblance with the definition of localizable entanglement, and therefore opens up a pathway to define a lower bound of LE according to equation (4). However, choice of an appropriate entanglement measure E for the state ρ k Ω still remains a challenge. The required characteristics of the chosen entanglement measure would be as follows.
(a) In order to exploit equation (8), the functional relation E min (ω) should be such that (b) The forms of the witness operator W k Ω suggest that the chosen measure E should capture the genuine multiparty entanglement in ρ k Ω . The first requirement indicates E min (ω) to be a linear function of ω, while the second requirement demands a computable genuine multiparty entanglement measure for mixed multiparty states. From here onward, we focus on regions Ω ≡ ab of size 2 (as in the example in figure 3), constituted of qubits, say, a and b, and choose negativity as the entanglement measure (see appendix C for definition) for mixed states. It has been shown in [93] that for a two-qubit state ρ ab and the corresponding expectation value ω of witness operator W having the form given in equation (5), the lower bound of entanglement, E min , as measured by negativity [116,117] and as a function of ω, can be obtained as E min(ω) = −2ω, which is linear in ω. Also, negativity captures the genuine multiparty entanglement over the two-qubit region ab. Therefore, from equations (4) and (8), which can be computed from the original state by measuring only the expectation value of the constructed local witness operator W Ω .

Measurement-based lower bound: graph-based method
Here we present a complete description of the methodology for determining the MLB of LE between any two qubits a and b in a large stabilizer state under noise. More specifically, we prescribe a logical choice for the specific Pauli measurements over the qubits in Ω which guarantees a non-trivial MLB, by exploiting our previously reported results [93] on graph states (a primer on the graph states and the graph transformations required for our purpose is provided in appendix D). In order to make the paper comprehensive for a general audience, we postpone discussions on the technical details of the different steps involved for the appendices.

Graph-based method: underlying mechanism
Let us denote the stabilizer state of N physical qubits, |ψ S , under noise represented by the map Λ(.): Any stabilizer state can be transformed to a graph state |ψ G [91] corresponding to a connected bicolorable graph G by local unitary transformations, such that where U S→G = ⊗ i U i , U i being either a single-qubit Clifford unitary operation, or the identity operator in the qubit Hilbert space (a detailed discussion on a geometric approach towards this transformation is given in section 3, and the algebraic details related to this methodology can be available at appendix E). The graph is bicolorable since the qubits situated on the nodes can be divided into two disjoint sets-a set of control qubits and a set of target qubits, where inter-set links are present in the graph, but intra-set links are not allowed. In the case of a graph state, we have provided a prescription for determining a non-trivial MLB of LE over a region Ω in [93] as long as Ω is connected, which, for a two-qubit region, implies the existence of a link between the two qubits constituting the region. However, in the present case, the underlying connected graph G may or may not contain the link (a, b) corresponding to the two-qubit region Ω ≡ ab.
In such scenario, a set of local complementation (LC) operations on a number of strategically chosen qubits situated on a selected simple path L ab connecting the nodes a and b may result in a graph transformation G → G , where the link (a, b) exists in G [93]. A sequence of LC operations on a graph G is equivalent to a local unitary operation U G→G = ⊗ i U i such that and U i is either a single-qubit Clifford unitary operation or the identity operator corresponding to the qubit i situated on the chosen simple path L ab . Therefore, the complete transition of the stabilizer state to the graph state |ψ G can be described by a unitary operation (see appendix F for a demonstration with a 7-qubit color code) Let us proceed with the assumption that the graph state |ψ G corresponding to the graph G containing the link (a, b) has been created from |ψ S via local unitary transformation U S→G = U G→G U S→G . In presence of the noise represented by the map Λ(.), where ρ = Λ (|ψ G ψ G |), and Λ(.) → Λ (.) is the transformation of the noise due to the local unitary operation. Note that the LE and the RLE for ρ S are the same as respectively the LE and the RLE of ρ due to the local Clifford unitary connection between ρ S and ρ . But computation of the LE and the RLE over the region ab in ρ still remains difficult in the case of large N which results in large m = N − 2, and one has to look for an appropriate Pauli measurement set-up in Ω which can provide a computable non-trivial MLB for ρ . The local unitary connection between ρ S and ρ can then be exploited to connect the value of the MLB with a specific measurement setup in the case of ρ S . We have shown in [93] that in situations where noise in the system is low, and the link between qubits a and b exists in the graph-state representation, an appropriate choice of measurement basis for a non-trivial MLB E P ab (ρ ) of E ab (ρ ) is local Z measurements on all the qubits except qubits a and b, which is an optimal basis for the LE over ab in |ψ G . The value of E P ab (ρ ) = E P ab (ρ S ) represents the value of the MLB corresponding to a specific Pauli measurement setup P for the original state ρ S , where P can be obtained by transforming the Z measurements on the qubits in Ω ≡ ab according to the over-all local unitary operation U S→G .
The graph-based algorithm is summarized in figure 4 (a pseudo-code for the algorithm can be found in appendix H). Evidently, the graph-based algorithm has two parts-transforming the stabilizer state to the graph states, for which the graph adjacency matrix is to be determined, and transforming the graph via local complementation operations to another graph in which a link exists between the chosen qubits. The first part and its different aspects have been discussed in section 3 and appendix E, and the algorithm has been transformed into a Python open-source package, namely, StabGraph [114], which generates the adjacency matrix of a graph corresponding to a graph state which is connected to a given stabilizer state via local unitary operators. For the second part of the graph-based method, the key challenge is to ensure the A schematic representation of the graph-based method to obtain non-zero measurement-based lower bounds of localizable entanglement, as presented in section 2.2. The input state is the stabilizer state of the topological quantum error correcting code under specific noise, and a pair of chosen qubits on which entanglement is to be localized. The next step is to determine the set of bicolorable graphs from the stabilizer structure by using Stabgraph (see section 3 and appendix E for the details on the procedure, and reference [114] for the package). For each graph in the set, if the chosen qubits are connected, the methodology developed in reference [93] is to be used to compute the lower bound of localizable entanglement. On the other hand, if the qubits are not connected, the adaptive local complementation algorithm is to be used in the form of the ALCPack on the graph so that a link between the two qubits is generated via a graph transformation (see appendix G for details on the algorithm, and reference [115] for the package). Once the link is created, the methodology developed in reference [93] can be used to compute the lower bound of localizable entanglement. The maximum of all of these values of lower bounds is to be chosen to obtain the best measurement-based lower bound of localizable entanglement on the two chosen qubits. certainty of creating a link between two chosen qubits by an optimal sequence of LC operations on a set of qubits in the graph. Towards this goal, we have developed an adaptive LC (ALC) algorithm, which, along with the corresponding graph transformations, is discussed in detail in appendix G. The crux of the algorithm depends on adapting itself according to the change in the graph after each local complementation operation on individual qubits, and subsequently choosing the qubit for the next local complementation operation according to the updated information. This algorithm has been made available as a package named ALCPack [115], which creates a link between any two given nodes in a simple, connected, and undirected graph via the adaptive local complementation method.

Numerical aspects of the graph-based approach
Note here that there exists a set of bicolorable graphs {G} that can be obtained from a specific stabilizer state |ψ S by appropriately varying the local unitary operation U S→G . Moreover, for each such bicolorable graph G where a and b are not connected, one can choose a set of simple paths so that LC operations on each of these paths would provide a graph G with the link (a, b). Let us denote the complete set of all possible bicolorable graphs G obtained from |ψ S by S G , having cardinality N G . In situations where the link (a, b) / ∈ G, let us denote the set of all possible simple paths L ab connecting the qubits a and b be S L ab , having cardinality N L ab . In the case of an arbitrary stabilizer state with large N, N G is usually a large number 6 . Also, for a given graph and a given pair of nodes {a, b}, determination of all possible paths 6 For a color code with N p plaquettes, there can be N p control qubits among N qubits, and the number N G varies as N N p .
between a and b can be difficult when N as well as the number of links in the graph is large. Therefore, for a large system S described by a stabilizer state |ψ S , it is difficult to obtain the optimal value of E P ab (ρ ). However, since each graph G results in a non-zero value of MLB, one can choose only a subset S G of S G with cardinality n G N G , and a subset S L ab of S L ab with cardinality n L ab N L ab for each G ∈ S G . This leads to a set of n G n L ab graphs G , for each of which a value of E P ab (ρ ) can be obtained. The bound max E P ab (ρ ) can be tightened by increasing the value of n G n L ab according to the available numerical resources. Note that the analytical computation of E P ab (ρ ) depends also on the noise model Λ. As shown in [93], in the case of local uncorrelated Pauli noise, E P ab (ρ ) is analytically computable for arbitrary system size as long as the neighborhood of the region Ω ≡ ab in G and the noise present in this region is fully known.
A word on the dependence of the run-time of the graph-based method, in particular, the adaptive local-complementation algorithm and the algorithm used to get the graph from the stabilizer structure, on the system size N is in order here. Since the ALC algorithm takes into account the transformed graphs at each of its steps, it is difficult to determine an exact dependence of the run-time of the algorithm on system size. However, one can determine a bound on how the run-time scales with N, by determining the maximum number of link operations during the ALC algorithm, which is N 3 (see appendix G for an explicit derivation). On the other hand, the Stabgraph algorithm uses Gauss elimination technique, and scales as ∼ N 3 . These indicate an overall polynomial scaling of the graph-based method with system size.
We point out here that one can also control the transformation S → G in such a way that the bicolorable graph G directly contains the link (a, b) (see section 3 for a discussion on how to ensure the creation of the link in G). Note that in the modified algorithm, which we refer to as the modified graph-based method, the optimization of E P ab (ρ ) is performed over a set of bicolorable graphs the corresponding states of which are connected to each other via local unitary operations. On the other hand, the former algorithm additionally uses graph states outside the set of bicolorable graphs. Therefore, the maximum value of E P ab (ρ ) obtained from the former algorithm is higher than the same obtained from the latter. See appendix H for a pseudo code of the modified graph-based method.

Graphs from topological codes: geometric approach
An arbitrary stabilizer state |ψ S describing a system S can be shown to be connected to a graph state |ψ G defined on a bicolorable graph G via a local unitary transformation [91]. While this transformation has an established algebraic formulation, it has also been shown [108] in the case of Kitaev's toric code that the graph G can also be constructed from the structure of the stabilizers via a geometric recipe. In this section, we introduce the geometric recipe for the topological color codes (see appendix A for the details on the topological color codes) and explain the underlying idea which emerges from the preparation protocols of the logical states of the code.
To discuss how a graph state can be obtained from a logical state of a topological color code via a geometric construction, we point out that the procedure for the creation of the logical states of a code on a lattice involves (1) initializing the qubits to either |0 or |+ so that they collectively form a product state, and then (2) creating plaquette-wise GHZ state-type [118] entanglement [73] via controlled entangling gates, which naturally labels one of the qubits as control (c), and the rest of the qubits as target (t) (see figure 5 and appendix I for an example with four qubits). The state on each of the plaquettes, however, can be further transformed to a graph state corresponding to a simple, connected, and undirected star-shapped graph with the control qubit c as the central qubit, and the target qubits t as the peripheral qubits, by applying local unitary transformations in the form of Hadamard operations on the target qubits. In terms of the stabilizer operators, the plaquette stabilizers corresponding to the plaquette are transformed to the graph-state generators via application of Hadamard operations on the target qubits. As an example, following this prescription, the |0 L state of the 7-qubit color code can be created by (1) choosing the qubits 1, 5, and 7, located at the corners of the triangular code, as control qubits controlling the rest of the qubits in their respective plaquettes, (2) initializing the seven qubits to the state |+ 1 0 2 0 3 0 4 + 5 0 6 + 7 , and (3) applying controlled phase gate to the (c, t i ) pairs. Subsequently, the local unitary connected graph state, which is obtained by applying Hadamard operations on all the target qubits in the 7-qubit lattice, corresponds to a graph that is obtained by creating the three star graphs with central qubits c = 1, 5, 7, and their respective target qubits (see figure 6(a)).
Note that in the case of topological color codes, the term plaquette does not always represent the original plaquettes of the color code lattice. It can be used in a broader sense, since a topological code can also be defined in terms of the products of its original plaquette stabilizers, defining the resultant plaquettes corresponding to the stabilizer operators obtained by multiplying two or more than two of the original Figure 5. Quantum circuit creating GHZ-type entanglement. The quantum circuit that takes the state |+000 to the four-qubit GHZ state |ψ P = 1 √ 2 (|0000 + |1111 ) (denoted by operation O 1 ), and then further takes it to the four-qubit graph state . The implication of applying this circuit to four qubits arranged as a plaquette of a topological color code has also been shown at the top, where the graph G P corresponding to the state ψ G P is created on the four qubits, starting from four qubits in the product state |+000 and via the single-plaquette GHZ state |ψ P . See appendix I for a detailed description. . . , 7} and initialized at |+ 1 0 2 0 3 0 4 + 5 0 6 + 7 , leads to the |0 L state of the 7-qubit color code, where qubits 1, 5, and 7 are chosen as control for the four groups of qubits respectively. Application of O 2 on the same groups of qubits, i.e., application of Hadamard operation on the target qubits 2, 3, 4, 6 in the next step leads to a graph state corresponding to a graph obtained by creating three star graphs on the three groups of qubits, {1, 2, 3, 4}, {5, 2, 3, 6}, and {7, 3, 4, 6}, where the control qubits 1, 5, and 7 are used as central qubits respectively. (c) If one assumes the qubits 3, 5, and 7 to be the control qubits, controlling the target quits in the plaquettes P 1 , P 1 P 2 , and P 1 P 3 respectively, then the resulting graph is obtained by combining the start graphs on the groups of qubits {1, 2, 3, 4}, {1, 4, 5, 6}, and {1, 2, 6, 7}, using respectively qubits 3, 5, and 7 as central qubits.
plaquette stabilizer operators of the original lattice. Therefore, applying the above prescription for creating GHZ state-type entanglement needs to be suitably generalized for larger codes. Nevertheless, the above recipe can be applied to all plaquettes in a topological code in its logical state, where the challenge is identifying the appropriate set of control qubits at correct positions of the code, and determining the correct sets of target qubits that are controlled by each of these control qubits. The graph corresponding to a graph state that is local unitary equivalent to the logical state of the topological code can then simply be obtained by creating all the star-shaped graphs that involve a control qubit and all its target qubits, with the control qubit as the central qubit. For example, the seven-qubit color code shown in figure 6(a) can also be expressed in terms of the plaquette stabilizer operators corresponding to the re-combined plaquettes P 1 , P 1 P 2 , and P 1 P 3 . This can be understood from the example shown in figure 6(b), where the qubit 3 controls the plaquette P 1 , while the qubit 5 controls the plaquette P 1 P 2 constituted of qubits 1, 2, 5, 6, 3, 4. Although there are two control qubits, qubits 5 and 3, among these, the target qubits {1, 6, 4} in plaquette P 1 P 2 are controlled by only the control qubit 5, while the target qubit 2 is controlled by none of the control qubits from the plaquette P 1 P 2 . The plaquette P 1 P 3 , on the other hand, is controlled by the qubit 7. Note that the x-type stabilizers associated to the plaquettes {P 1 , P 1 P 2 , P 3 } are an equivalent subset of the x-type generators associated to the plaquettes {P 1 , P 2 , P 3 }, which ensures the validity of the geometric recipe. It is also clear from the above discussion that the number of chosen control qubits has to be equal to the number of plaquettes in the topological color code.
In view of all these aspects, for an arbitrary topological color code with an appropriate choice of the set of control qubits and the set of target qubits that are controlled by them, one has to ensure the following conditions: (a) a qubit, once chosen as a target (control) qubit, can not be chosen as a control (control or target) qubit any more, and (b) the support of all the plaquettes and the logical operators, as well as the products of them, must contain at least one control and one target qubits .
Note here that the condition 1 ensures that the transformation from the logical state of a TCC to graph state ensures only (c, t i )-type links, thereby dividing the qubits into two mutually disjoint sets of control and target qubits, and ensuring that the resulting graph is bicolorable, i.e., (c, c) and (t i , t j ) links are not present. Note also that in the modified graph-based method discussed in section 2.2 one needs to ensure creation of a link between two given qubits. This can be done by choosing a path of adjacent plaquettes so that the given qubit-pair is contained by the large plaquette constituted by the plaquettes on the path, and then using one of the given qubit-pair as control qubit to create the graph by creating the (c, t i )-type links (see figure 7 for a demonstration).
The success of the geometric method depends explicitly on the correct choice of a set of control qubits, and the determination of the sets of target qubits that are controlled by the control qubits, ensuring the conditions (i) and (ii). While this is possible irrespective of the size of the code, the choice of a correct set of control qubits may prove difficult in the case of larger codes. This leads us to an algebraic approach for creating the graph, exploiting the binary picture of the code. Note that an algebraic treatment of the same problem was considered in [91], although the connection between the geometric recipe and the algebraic method was absent. We revisit the treatment, and provide a slightly modified version of the algebraic calculation in appendix E for a general stabilizer state. The geometric recipe for the color code presented above is inherent in the algebraic approach for determining the adjacency matrix of the graph from the stabilizer structure of the color code.

Application: color code on a square hexagonal lattice
In this section, we apply the methodology developed through sections 2 and 3 to the case of a color code hosted in a square hexagonanl lattice with open boundary condition on a plane, where the 'square' indicates the shape of the lattice (our methodology works irrespective of whether the lattice is square or triangular; for an example of the 'triangular' lattice, see the 7-qubit code shown in figure 4(b)). The number of qubits, N, in the code is represented by the distance, D, of the code, where, for the square hexagonal lattice, N increases quadratically with increasing D (see appendix A). An example of the square hexagonal lattice hosting a color code with D = 4 and containing two logical qubits is given in figure 8(a), where N = 18. We will be computing the lower bound of LE over a qubit pair {a, b} in the bulk, which is constructed from the lattice by removing D/4 qubits in the direction towards the center from each boundary (see figure 8(b)). Note that the choice of the control and the target qubits depends explicitly on the choice of the stabilizer state. In the present case, we consider the color code to be in |+ L for all our discussions, and the local unitary connected graph state is obtained by applying Hadamard operations over the control qubits. This choice is justified as the logical Pauli-eigenstates in 2D color codes are connected by local unitary operators due to the transversality of logical Clifford gate operations [63,64].

Computation of witness-based lower bound
Here, we explicitly compute the expectation value ω of the local witness operator W ab of the form given in equation (5), constructed for a pair {a, b} of qubits in the bulk. As discussed in section 2.1 and demonstrated in figure 3, we build these witnesses from two stabilizers denoted by S x and S z , one being x-type and the other z-type, obtained by multiplying respectively the x-and z-type stabilizers corresponding to the plaquettes on two adjacent paths of plaquettes connecting the qubits a and b, so that the necessary conditions for the construction of W ab (see section 2.1) are satisfied. The path constituted of the common lattice-links between the adjacent paths of plaquettes provide the path connecting the two chosen qubits, and the length of this path is the distance d between the chosen qubits. In all our discussions, we compute the distance between the chosen qubits with respect to the square hexagonal color code lattice.
Expanding the form of W ab from equation (5) using S x and S z , the WLB −2ω is calculated as where ω x(z) = Tr ρ S S x(z) , and ω xz = Tr [ρ S S x S z ]. For demonstration, we consider the state ρ S originating from the application of single-qubit uncorrelated Pauli noise channels [119,120] to the stabilizer state |ψ S (see appendix K for a brief description), including the bit-flip (BF), phase-flip (PF), bit-phase-flip (BPF), and depolarizing (DP) noise. For the BF (PF) noise, one can show that ω x = 1(ω z = 1) and where q is the strength of the noise (0 q 1) which we asssume to be the same for all qubits, and n x (n z ) is the number of qubits in the support R x (R z ) of the stabilizer S x (S z ). See appendix L for the calculation in the case of PF noise. The calculation in the case of the DP channel is similar to the same in the case of the PF channel, where the expectation values are ω x = (1 − q) n x , ω z = (1 − q) n z , and ω xz = (1 − q) n x +n z −2 , where n x + n z − 2 is the number of qubits in the support of S x S z . Using these, equation (17) becomes Note from the design of the local witness operator (see figure 3) that the types of the plaquette stabilizer operators corresponding to each of the two adjacent paths of plaquettes-one above and the other below the path made of lattice-links connecting qubits a and b-are different from each other, one being z-type while the other x-type. One obtains a valid local witness operator even when x-and z-types of the stabilizers above and below the path connecting qubits a and b, and contributing to W ab , are interchanged. Note also that the values of n x and n z depends on the distance between the chosen qubits, implying that the dependence of WLB on d is decided by how n x and n z grow with increasing d. The exact dependence of n x and n z on d depends explicitly on the construction of W ab , and the layout of the path made of lattice-links connecting the qubits a and b on the lattice. In general, the number of plaquettes involved in the local witness operators constructed in this way, and therefore the support of the stabilizers grows as ∼ a + bd, with some constants a and b, implying an exponential dependence of ω on d. As an example, let us consider the layout of the path of lattice-links connecting the qubits a and b in the bulk as shown in figure 3(b), for which

Computation of measurement-based lower bound
We now investigate the variation of MLB with the distance between the chosen qubits in the bulk of the square hexagonal lattice of code-distance D. Before moving on to the analysis of the numerical results, a word on the notion of the existence of local and non-local links in the graph, and its relation with the topological properties of the system is in order here. In the color code lattice, the links are all local links since they connect qubits belonging to a specific plaquette. This notion of locality comes from the fact that the plaquette stabilizer operators are intrinsically local in the sense that they operate on the qubits belonging to the same plaquette. Therefore, a graph constructed following the lattice of the color code (where one considers lattice sites as nodes in the graph where the qubits are situated, and introduces links in the graph according to the links in the color-code lattice-see figure 8(c) for an example) contains only local links. However, the graph obtained from the code using the methodology presented in section 3 may contain a number of non-local links connecting a control and a target qubit belonging to two different and distant plaquettes (see figure 8(d)), which is in contrast with the characteristics of a local graph. It has been shown that small and simple setups diminishes the effect of these non-local links, and a critical size is to be achieved in order to observe the effect of the topological properties in terms of the existence of the non-local links in the graph [108]. However, one can also take a different perspective, and ask whether a differentiation can be made in terms of entanglement. Our graph-based algorithms are appropriate for such investigations.
We use the developed packages StabGraph and ALCPack to apply the graph-based algorithms H1 and H2 to determine the MLB for LE over qubit-pairs situated in the bulk of the lattice of the square hexagonal code. In the case of the graph-based algorithm H1, for a pair of chosen qubits separated from each other by distance d, we set n G = 1, and optimize the MLB of LE over a set of graphs G generated via LC operations on the qubits situated on randomly chosen paths L ab connecting the qubits a and b, where the size of the set of graphs G equals to the number of random paths n L ab . We first test how the ALC algorithm scales with the system size by looking at the average time t taken by the ALC algorithm to create a link between the chosen qubits a and b. We vary the code-distance D of a square hexagonal code as D = 12, 16, 20, 24, such that the number of qubits in the system are N = 194, 354, 562, 818 respectively. In figure 9, we present the variation of t as a function of N for qubit pairs with different distances d = 2, 4, 6, 8, 10, where the average value t is determined over a sample of size n L ab = 10 4 for each value of d. We separately consider a local graph that follows the structure of the square hexagonal lattice, and a non-local graph obtained from the TCC defined on the square hexagonal lattice by using the methodology described in section 3. From figure 9, it is evident that t increases with N for a fixed d, and increases with d for a fixed N when the graph is local, which is in contrast with the variation of t with N in the case of a non-local graph obtained from the square hexagonal code. In the latter case, t increases only negligibly with d for a fixed value of N. This can be understood from the fact that there exists considerable number of non-local links in the case of  the local unitary equivalent graph obtained from the TCC, which results in comparable lengths of the paths connecting the chosen qubits in the graph irrespective of the actual distance d between the qubits. This leads to a similar number of required LC operations, which results in slowly increasing values of t with d for a fixed N. On the other hand, in the local graph, the typical length of a path L ab connecting a and b increases with increasing distance d between qubits a and b in the TCC lattice, subsequently increasing the required number of LC operations, and therefore the average value of t . Also, the variation of t with increasing N for a fixed d clearly validates the polynomial scaling of the ALC algorithm as discussed in section 2.2. Given the above discussion, it is interesting to investigate whether the average number of LC operation, n LC , required to create a link between two chosen qubits in the bulk varies with the distance between the qubits, when the system-size is fixed. Figure 10 depicts the variation of n LC with d in the case of the local and non-local graphs corresponding to the square hexagonal color code lattice of D = 20, where the averaging has been done over a sample size of n L ab = 10 4 for each value of d. In the case of the local graph, n LC rapidly increases with increasing d, while for the non-local graph, the increasing trend of n LC slows down considerably when d increases. These results are in agreement with the variations of t against d for a fixed value of N.
Since both U S→G and U G→G are constituted of local Clifford unitary operators (see section 2.2), the transformed noise Λ is also local uncorrelated Pauli noise similar to Λ, although the individual bases in which the noise processes take place corresponding to Λ on each qubit may differ from that in Λ. As shown in [93], for uncorrelated single-qubit Pauli noise applied to a graph G in which a link between the two chosen qubits is present, a high value of MLB is favourable when the size 'n' of the neighborhood of the qubit-pair {a, b}, for which the noise does not commute with the Z-measurement, is low. In figures 11(a) and (b), we plot the variation of the minimum value of n, represented by n min , as a function of d for the local as well as non-local graphs corresponding to the square hexagonal code with D = 20, where the minimization of n is achieved over a sample size of n L ab = 10 4 for all values of d, in the case of (a) the DP and (b) the PF noise. Note that in the former case, the value of n equals the number of qubits in the full neighborhood, while in the latter, n is the number of qubits in the neighborhood with BF or BPF noise. This implies a higher value of n in the former case than the latter, which is clearly demonstrated also in the values of n min in the figure 11. In both cases of the local and the non-local graphs corresponding to the TCC, and for both types of noise, the value of n min increases monotonically with d.
Next, we plot the natural logarithm of the maximum value of MLB, denoted by max E P ab (ρ ) = max E P ab (ρ S ) and computed following the methodology developed in [93] as a function of d. The plots are shown in figure 11(c) (for the DP noise) and figure 11(d) (for the PF noise), where we choose the noise strength q = 10 −2 , and the maximization is achieved over the same set of n L ab graphs as in the cases of n LC and n min . We fit the variation of the value of the natural logarithm of max E P ab (ρ ) with d using the equation ln max E P ab (ρ ) = a + bd, such that the MLB decays exponentially with d according to the equation max E P ab (ρ ) = ae bd with a = lna. Here, a and b are fitting parameter, which are expected to be functions of the noise strength q. See figure 11 for the values of the fitting parameters a and b, obtained in the example.
So far, we have considered a variant of the graph-based algorithm by setting n G = 1, n L ab = 10 4 . We now apply the modified graph-based algorithm in order to investigate the features of MLB where a link between a and b is created every time that one transforms the stabilizer state into a graph state. Therefore, the optimization this time is over a large number n G of bicolorable graphs G, obtained directly from S (see section 2.2). The data obtained for n min and ln max E P ab (ρ ) as functions of d by using graph-based algorithm H2 are presented in figure 12. The numerical results presented in this section are illustrations of the applicability of the witness-and graph-based methodologies developed and discussed in section 2 for determining non-trivial lower bounds of the localizable entanglement over bulk qubit-pairs in the case of arbitrary TCC lattices, such as the square hexagonal lattice with open boundary condition. The lower bounds decrease with increasing d, and the trend agrees with our understanding of the growth of the size of the neighborhood with noise having bases that do not commute with Z-measurements around the link connecting the qubits in the chosen pair, with the increasing distance between them. However, in order to infer the exact dependence of the MLB on d, one could consider the optimality of the algorithms for obtaining the maximum value of the MLB. This issue could be thoroughly investigated, which is beyond the scope of the present work.

Conclusions and outlook
In this paper, we provide two specific pathways to estimate the lower bound of localizable entanglement over a subset of qubits in a large system of noisy topological codes, including the surface and the color codes. In one approach, we use an appropriately constructed local entanglement witness operator constituted of the stabilizer generators of the code, and estimate a lower bound of localizable entanglement over the chosen region of qubits via the expectation value of the witness operator. We also propose a specific construction of the the local witness operator for this purpose. On the other hand, we also propose a methodology for estimating a lower bound of localizable entanglement corresponding to a specific projection measurement setup on the qubits outside the specified subset of qubits in a noisy stabilizer state. This uses the fact that an arbitrary stabilizer state can be connected to a graph state via local unitary transformations. We discuss in detail a scalable geometric recipe for determining the graph underlying the local unitary connected graph state from the stabilizer state describing a color code. We also present an algebraic methodology to determine the adjacency matrix of the graph corresponding to a graph state obtained via local unitary transformation from a stabilizer state of arbitrary size. Moreover, we develop an algorithm which creates a link between any two chosen nodes in a simple, connected, and undirected graph by using a sequence of local complementation. We also develop appropriate numerical packages [114,115] for these purposes. We predict that the runtime of the graph-based algorithm scales polynomially with the system size, which is supported by our numerical findings corresponding to a topological color code on a square hexagonal lattice.
We also determine the witness-and measurement-based lower bounds of localizable entanglement in the case of qubit pairs situated in the bulk of a topological color code described on a square hexagonal lattice. We explicitly compute the expectation value of the local entanglement witness operator constructed according to our prescription in the case of a stabilizer state of the system under local uncorrelated Pauli noise. Our calculations show that the bound obtained from the proposed construction of the local witness operator exponentially decreases with increasing support of the witness operator, and therefore with increasing distance between the chosen qubits. In the case of the measurement-based method, along with computing the bounds of localizable entanglement over a qubit-pair in the bulk of a topological color code via the graph-based methods I and II, we also determine the bounds of localizable entanglement corresponding to a qubit pair in the case of a local graph that follows the square hexagonal lattice.
Note that the numerical data presented in this paper corresponding to the measurement-based lower bound is obtained as a proof of the functionality of the graph-based method and the modified graph-based method developed in this paper. It could be interesting to study the optimality of the algorithm in order to optimize the bound. The algorithms proposed in this paper could also be generalized for regions beyond two qubits leading to the notion of multiparty localized entanglement [84] in topological quantum codes, which could allow to reveal long-range multiparty quantum correlations. have been performed on the Swansea SUNBIRD system. The Swansea SUNBIRD system is part of the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via the Welsh Government. We thank Davide Vodola for letting us use the numerical tools developed by him to construct the square hexagonal lattice. We also thank Ciarán Ryan-Anderson for useful discussions.

Appendix A. Topological color codes
We use topological color codes (TCC) [63,64] as the testing ground for our results, and the defining features of a TCC are briefly discussed in this appendix. A TCC model is constructed on a two-dimensional (2D), three-colorable, and trivalent lattice, where each vertex of the lattice contains a physical qubit, and the lattice can be embedded on a compact surface having arbitrary topology of genus g (for example, a torus with g = 1). The three-colorability of the lattice implies that the faces of the lattice, also known as the plaquettes, can be painted with three different colors, where neighboring plaquettes always have different colors. The trivalency of the lattice implies that each vertex is connected to three links. The lattice can also be characterized by coloring the three links connected to each vertex with the same set of colors as the colors of the plaquettes, such that each link has a color different than both colors of the plaquettes sharing the link. There is a number of such regular lattices available, such as the hexagonal (honeycomb) lattice, the square-octagonal lattice, and the square-hexagon-dodecahedron lattice. The number of logical qubits in a TCC is given by k = 4 − 2χ, where χ is the Euler characteristic of the surface, thereby ensuring that k depends on the topology of the surface. The methodology developed in this paper are aimed for arbitrary stabilizer states, and therefore applies to arbitrary TCC. For the purpose of testing our prescriptions, in this paper, we shall focus on the 2D honeycomb lattice with open boundary conditions. To keep the figures uncluttered, we shall use the three-colorability of only the plaquettes for demonstration [see figure 1 for a schematic of a honeycomb lattice, where the three colors corresponding to the TCC are red (R), green (G), and blue (B)].
We now set up the terminology for the stabilizer description of a TCC, which requires the definition of two key concepts, (1) the stabilizer subspace, and (2) the logical operators. The stabilizer subspace of the TCC is determined by the stabilizer group of operators, which is generated by a set of plaquette operators denoted by S α p . There are two types of plaquette operators, corresponding to α = x and α = z for each plaquette P p , called the x-type and the z-type operators, given by where X and Z are respectively the x and z components of Pauli matrices. Each of the plaquette operators squares to the identity, i.e., S α p 2 = I ∀ p, α = x, z, and they mutually commute, [S α p , S α p ] = 0, since all plaquettes have an even number of vertices in a trivalent three-colorable lattice, and since adjacent plaquettes share even number of vertices. The stabilizer subspace H S in the full Hilbert space H of the TCC is given by Note that the plaquette stabilizer operators are intrinsically local in the sense that they operate on the qubits belonging to the same plaquette. This notion of locality is crucial for the discussions presented in this paper. In all our considerations, we define the physical distance between any two lattice points by the length of the shortest path constituted of the lattice links and connecting the two lattice points. A TCC hosts a total of 2k independent logical generators L (q) α , where q = 1, 2, . . . , k is the number of logical qubit, and α = x, z indicates the local Pauli matrices corresponding to the lattice sites that constitute the logical operators. They are defined on homologically non-trivial (i.e., non-contractible) strings across the lattice, and they commute with all stabilizers: The x-and z-type logical operators define the computational basis {|0 L , |1 L } for the logical qubits, such that where the subscript 'L' denotes the logical states. The eigenbasis of L (q) x and L (q) y in terms of {|0 L , |1 L } are |± L = (|0 L ± |1 L )/ √ 2 and |±i L = (|0 L ± i|1 L )/ √ 2 respectively. In this paper, we focus on a 2D color code defined on a square-hexagonal lattice, as described in section 4. The number of physical qubits, N, in a 2D topological color code depends on the code distance D. In the case of the topological color code defined on the square-hexagonal lattice with code distance D, N varies with D as where D = 4l(l = 1, 2, 3, 4, . . .).

Appendix B. Decomposition of local witness operators
In this appendix, we present the detailed calculation for decomposing a local entanglement witness operator of the form in equation (5) into the form given in equation (7). In order to decompose W Ω in terms of local projection operators on Ω and witness operators W k Ω on Ω, notice that the commutation property in (i) splits Ω into two regions, (1) Ω 1 constituted of qubits for which u i,j = 0 ∀j, and (2) Ω 2 consisting of qubits such that there exists at least one stabilizer S j ∈ S in which u i,j = 0. Therefore, for qubits i ∈ Ω 1 , ⊗ i∈Ω 1 τ Ω 1 u i,j = I Ω 1 , I Ω 1 being the identity operator in the Hilbert space of the qubits in Ω 1 . On the other hand, in the case of a qubit i in Ω 2 , if u i,j = 0 for more than one stabilizers S j ∈ S, then the values of u i,j are identical for all j for which u i,j = 0. This assigns a specific Pauli operator τ Ω 2 v i (v i = 1, 2, or 3) for the qubits i ∈ Ω 2 . The corresponding projection operators can be written as identity operator in the Hilbert space of the qubit i ∈ Ω 2 in stabilizer S j , and k i (= 0, 1) can be interpreted as the outcome of the projection measurement via P Ω 2 (k i ,v i ) . The constructions of the local witness operator has Figure B1. The supports of the witness operator in figure 3(b). The two-qubit region is marked by the white nodes, while the black, turquoise, and yellow nodes mark the region Ω. The black nodes represent the sub-region Ω 1 , while the turquoise and yellow nodes stand for the sub-region Ω 2 . The two different colors, turquoise and yellow on a qubit i ∈ Ω 2 represent τ been demonstrated in figures 3(a) and (b). The supports of the stabilizers S z and S x , denoted by R z and R x respectively, are the sets of nodes corresponding to which the Pauli operator contributing to the stabilizer is not an identity. For example, the sizes of R z and R x , denoted respectively by n z and n x , are n z = 8, n x = 6, for the figure (a), while for (b), n x = n z = 8. See figure B1 for an illustration of the supports is the case depicted in figure 3(b).
Before breaking the mathematics any further, let us consider the effect of the application of a projection operator P Ω 2 (k i ,v i ) on each of the qubits in Ω 2 , which results in Therefore, application of a projection operator of the form P Ω 2 (k i ,v i ) on each qubit i ∈ Ω 2 yields with η j = i∈Ω 2 η i,j , and . v m are multi-indices 7 , and we denote the size of Ω 2 by m ( m). From equation (B3), it is clear that the application of local projection operations in the basis of Pauli operators fixed by the witness W Ω outside the region Ω results in local stabilizers S Ω j . These stabilizers correspond to the genuine multiparty entangled state |ψ Ω -the same state that is obtained over the region Ω by performing projection operation P Ω 2 (k,v) on the stabilizer state |ψ [see condition (ii)]. We now point out that the completeness of the basis formed by the eigenstates of the Pauli matrices in the Hilbert space H i of a qubit allows one to write I i as This can be extended to the identity operator in the Hilbert space H = ⊗ n i=1 H i of the full system, which, when multiplied with W Ω from the left and expanded in terms of the projectors P Ω 2 (k,v) yields where I Ω is the identity operator in the Hilbert space of the qubits in Ω. Note here that given the condition (ii) for the construction of witness operators, W k Ω detects genuine multiparty entanglement in the region Ω. Note also that one can use equation (B5) for the identity operators in the Hilbert space of the qubits in Ω 1 as well, by choosing any one of the three Pauli operators, X, Y, or Z to define the projectors, since the form of W k Ω remains independent of this choice. Each set of Pauli operators defining the projectors over the region Ω 1 leads to a specific Pauli measurement setup over the qubits in Ω. For each of these setups, equation (B6) takes the form with {m 2 , m 3 , . . . , m n−1 } being the nodes that are visited while traveling from the source a to the target b along L ab , where none of the nodes is repeated in the sequence, and the link (m i , m i+1 ) ∈ G, i = 1, . . . , n − 1. The number of links traversed while going from a to b is the length l = n − 1 of the path. We denote a simple path between a and b having length l as L (l) . To keep our notations uncluttered, we discard the subscript 'ab' and the superscript '(l)', and denote a simple path of the form given in equation (D3) by L, unless we need to distinguish between two paths of different source and/or target nodes, or of different length. Shortest paths. There can be more than one simple paths of different or same lengths between two nodes a and b in a graph. The shortest path is the simple path between a and b having the minimal length l = l min , where the minimization is taken over all possible simple paths between a and b. There can be more than one shortest paths between a specific pair of nodes.
where the LC operation is performed on the node m 1 first, and then according to the sequence {m 1 , m 2 , . . . , m n }. Note that this specific sequence is important, as LC operations, in general, do not commute. Stabilizer formalism. The stabilizer description of a graph state |ψ G corresponding to a graph G uses N generators g i of the form which forms a subset of the Pauli group. The generators {g i } mutually commute, thereby sharing a common eigenbasis, and the graph state |ψ G represents the common eigenstate with eigenvalue +1. The rest of the eigenstates of the generators can be represented as

Local complementation as local unitary operation.
An LC operation on a graph G with respect to the node i is equivalent to a local unitary transformation of the corresponding graph state |ψ G , where the unitary operator is given by A sequence of LC operations on a number of chosen nodes in a graph G, denoted by equation (D4), is given by where U i is as in equation (D6).

Appendix E. Algebraic approach
In this appendix, we put the graphical recipe discussed in section 3 in a mathematical footing, and show how the adjacency matrix of the graph can be obtained from the structure of the stabilizers, revisiting the results obtained in [91]. We use the binary picture [91,92] for the description of the stabilizer state of the N-qubit system, which is a 2N × N binary matrix composed of two N × N blocks corresponding to the Zand X-type stabilizers, given by Here each column of A, marked by the index j, represents one stabilizer S j , in which the Pauli matrix corresponding to the qubit i is represented by where Z i,j , X i,j = 0, 1 (see appendix J for details). Since the stabilizers, or more precisely, the generators are independent, the matrix A is of rank(A) = N. Assuming that the rank of X is n, the stabilizer state can be represented by a new set of stabilizers obtained by re-combining the original stabilizers through Gaussian elimination among the columns in such a way that X ij = 0 for j = n + 1 . . . , N, where we arrange the columns as j = 1, 2, 3, . . . , n, n + 1, . . . , N without any loss in generality. We point out here that all the operations in the binary picture are performed modulo 2. Denoting 'left' by l and 'right' by r, the stabilizer state now becomes where X l is an N × n matrix of full rank. Note here that in the case of a CSS code [119], Z l can be considered to be 0, ensuring that the left block of columns in A corresponds to only x-type stabilizers, while the stabilizers represented by the right block of columns are only z-type.
Since the matrix X l has rank n, the same number of linearly independent rows can be chosen from it, which will then form n × n invertible matrix X l,c corresponding to the n control qubits, the subscript c denoting control. As a consequence, X l,c is an invertible matrix. The rest of the rows in X l stand for the target qubits, denoted by the subscript t. We always label the rows as i = 1, 2, . . . , n, n + 1, . . . , N from the top, and consider the first n rows to be corresponding to the control qubits, while the rest N − n rows are for target qubits. Therefore the matrix A takes the form Our approach towards extracting the adjacency matrix Γ underlying a local unitary connected graph state from a stabilizer state represented by A can be summarized via the following equation: where R is an invertible binary matrix, and Q is constituted of local Clifford unitary operations, having the form Q = U Z H target . Here, U Z is a local π 2 rotation w.r.t. the z-axis on a subset of the control qubits that we are going to specify later, and H target represents Hadamard operations on all the target qubits. In the subsequent discussion, we explicitly calculate the L.H.S of equation (E4), and demonstrate the extraction of Γ.
The first step is to apply Hadamard operations on all the target qubits, which, in the binary picture, implies the interchange of the elements above and below the horizontal line in equation (E2), i.e., It can be proved that the lower block of H target A is invertible, and its inverse, R, can be determined. To do this, we take a different approach that the one used in [91]. We apply to H target A a Clifford operation represented by U that makes the upper block vanish and leaves the lower block unchanged. The unitary that we construct for the purpose of the proof is represented by where explicit forms of B, C will be given in equations (E15) and (E16) respectively. The matrix UH target A is given by ⎛ From the form of the matrix C in equation (E15), the non-zero diagonal terms of the upper block vanish, while the off-diagonal terms vanish due to equation (E11). Hence the proof of the invertibility of the lower block of H target A . Note that Clifford operations are represented by full-rank matrices, and therefore U has preserved the rank N of H target A . Consequently, the matrix UH target A , which contains only the lower block, is of full rank, and therefore, the lower block is invertible.
In order to determine R, we note that the lower block of H target A has the form of a lower triangular matrix, and the diagonal terms X l,c and Z r,t of the lower block of H target A are also invertible due to the invertibility of the lower block of H target A . Therefore, R is another lower triangular matrix by blocks having the form where V must satisfy that Z l,t X −1 l,c + Z r,t V = 0, leading to V = Z −1 r,t Z l,t X −1 l,c . The lower block of H target A , when multiplied by R from the right, becomes the identity, while the upper block becomes We now show that the upper block of H target A R is symmetric, for which we exploit the fact that every stabilizer state, say, A , due to the communitativity between all pairs of stabilizers, has to satisfy A T DA = 0, where D is a 2N × 2N binary matrix with zeros in the two N × N diagonal blocks, and N × N identity matrices in the off-diagonal blocks [92] (also see appendix J). This leads to Multiplying equation (E11) from the left hand side by X T l,c −1 , and from the right hand side by Z −1 Next, we multiply equation (E10) from the left by X T l,c −1 and from the right by X −1 l,c to obtain Use of equation (E12) leads to the modified form of the upper block of H target A R as where with C being symmetric, given equation (E13). Therefore, in view of our goal to explicitly calculate Γ from equation (E4), we have now obtained For the above matrix Γ to be the adjacency matrix Γ of a graph, the diagonal elements of C must vanish, which is achieved by multiplying H target A R by U Z from the left, given by (see appendix J) where U Z represents the single-qubit π 2 rotations with respect to the z axis on the control qubits for which C ii = 1. This leads to the adjacency matrix Γ of the form (E19) Based on the above discussion, one can develop an algorithm for obtaining the adjacency matrix Γ corresponding to the graph G representing the graph state that is local unitary equivalent to the stabilizer state, by utilizing the stabilizer structure of the state only. This algorithm has been realized in the form of a Python open-source package called StabGraph [114], which generates the adjacency matrix corresponding to the graph underlying the graph state which is connected via local unitary operations with the stabilizer state (see appendix H for a pseudo code of the algorithm). Our recipe using the Gaussian elimination technique scales as ∼ N 3 with the system size, N. Note here that in the case of CSS codes [119], Z l,c and Z l,t in A can be set to zero by stabilizer recombination, implying C = 0, thereby ensuring that the adjacency matrix (E20) has vanishing diagonal blocks. This also implies that links between a pair of control qubits and a pair of target qubits are prohibited, thereby ensuring that the resulting graph is bicolorable, although in the general stabilizer state, the existence of (c, c)-type links is allowed. Moreover, transformation of the stabilizer state to the graph state does not require the application of the local unitary operation U Z , if C = 0.

Appendix F. Graph-based algorithm on a 7-qubit color code
The graph-based method for computing the measurement-based lower bound of LE has been discussed in section 2.2. In this appendix, for illustration, we apply the graph-based method described in section 2.2 on the smallest 2D color code, constructed over a triangular lattice of seven qubits arranged in three adjoined plaquettes, where one physical qubit is placed at each vertex. The stabilizer operators corresponding to the 7-qubit color code are constituted of 4-qubit X-and Z-type operators for each plaquette, given by Let us assume that the two-qubit region is given by the qubits a ≡ 1 and b ≡ 5 (see figure G1) D). Therefore, the unitary operator U S→G consists of the local unitary operators U 2 = exp(−iπX 2 /4), U 3 = U 4 = U 6 = I, and U 7 = H 7 in the region Ω ≡ ab. Consequently, E P ab (ρ ) as an MLB with P ≡ Z measurement on all qubits ∈ ab in ρ is equivalent to Y measurement on qubit 2, X measurement on qubit 7, and Z measurements on qubits 3, 4, and 6, denoted altogether by P in ρ S , so that E P ab (ρ ) = E P ab (ρ S ).

Appendix G. Graph transformations via adaptive local complementation
As mentioned in section 2.2, for the successful implementation of the graph-based protocol, one needs to develop a graph trannsformation algorithm that creates a link between any two chosen nodes in a simple, connected, undirected graph via a sequence of local complementation operations over a set of nodes in the graph. We have developed an adaptive local complementation technique that updates itself at every step of the graph transformation via local complementation over one qubit, and chooses the next node for local complementation using the updated information. Here, we discuss the technical details of the adaptive LC technique developed for the graph transformation G → G used in section 2.2. The prescription works for a generic simple, connected, and undirected graph G, irrespective of whether it is bicolorable or not. We shall denote a link connecting two nodes i and j by (i, j). We shall also represent a simple path in the graph, connecting the two nodes a and b and having length l, by L (l) ab , which, without any loss in generality, is given by (see appendix D) (G1) Figure G1. Transformation of a 7-qubit topological color code to a graph G by using the graph-based algorithm described in section 2.2, where the nodes 1 and 5 are connected by a link in G . A non-zero value of E P 15 (ρ ) in G requires local projection measurement in the Z-basis at all other qubits. This corresponds to X-and Y measurement on qubits 7 and 2, respectively, and Z-measurement everywhere else. See appendix F for details.
Here, {2, 3, . . . , n − 1} are the nodes that are visited while traveling from a to b along L (l) ab , such that (i, i + 1) ∈ G, i = 1, . . . , n − 1, and none of the nodes is repeated in the sequence. The number of links traversed while going from a to b is the length l = n − 1 of the path. To keep the notation uncluttered, we will drop the subscript ab and the superscript (l) from the notation of the simple path unless it is explicitly required for comparison. For the purpose of this discussion, one needs a specific classification of the simple paths connecting two nodes, and a definition of the distance between them, as follows. Category of simple paths. We first divide the set of all possible simple paths between two given nodes a and b into two categories, C 1 and C 2 . A simple path L of the form given in equation (G1) belongs to the category C 1 if a qubit on L is not connected to another qubit on L via a direct link, unless the second qubit is right before or after the first qubit in the sequence L. Formally put, for two qubits i and j on L, iff (i, j) / ∈ G ∀ i, j ∈ L, i < j n, j = i + 1, then L ∈ C 1 . Otherwise, L ∈ C 2 (see figure G2 for an example). Note that C 1 ∩ C 2 = ∅, while C 1 ∪ C 2 constitutes the set of all possible simple paths between a and b. Distance between two nodes. There can be a number of ways in which the distance between a pair of nodes in a graph can be quantified. However, the distance between two nodes i and j that is specific to a given path L of the form in equation (G1) is unique. It is measured by the number of links that one has to traverse while going from i to j along L, and can be represented by d L (i,j) = j − i, where we have assumed j > i without any loss of generality. Note that for a path L, d L (1,n) = l = n − 1. All shortest paths between two given nodes a and b in a simple, connected, and undirected graph belong to C 1 . This can be proved by assuming that there exists a shortest path (see appendix D for a definition) L between two nodes a ≡ 1 and b ≡ n that belongs to C 2 . However, by definition of C 2 , there exists at least one link (i, j) such that i, j ∈ L, i < j n, and j = i + 1. This implies that there exists a simple path between a ≡ 1 and b ≡ n, given by of length l = l − (j − i) + 1, which is < l since j > i + 1 by the definition of C 2 . Therefore, L is not a shortest path. However, note that not all paths in C 1 are shortest paths (see figure G2 for examples).
We are now in a position to prove the following theorem, which will be crucial for our purpose. Proof. Let us consider the path L ∈ C 1 , as given in equation (G1), between the nodes a and b in a simple, connected, and undirected graph G. Let us denote the neighborhood of the node 2 as N 2 . The definition of C 1 implies that only 1, 3 ∈ L belongs to N 2 , while {4, . . . , n} / ∈ N 2 . Let us assume that the LC operation (see appendix D) on node 2 transforms the graph as G → G 2 = τ 2 (G). The definition of LC operation implies creation of the link (1, 3) thereby introducing 1 in the neighborhood N 3 of 3, while keeping the neighborhoods of the rest of the nodes {4, 5, . . . , n} unchanged. Also, the LC operation τ 2 does not affect the links {(i, i + 1), i = 3, 4, . . . , n − 1}. Therefore, in G 2 , there exists a simple path of length l = l − 1, given by [a ≡ 1, 3, 4, . . . , n ≡ b] ∈ C 1 , between the pair of nodes a and b. One can continue performing a total of n − 2 LC operations successively on the nodes 2, 3, . . . , n − 1 in the same order as they are in the sequence L, where during each individual LC operations on the node i, i = 2, 3, . . . , n − 1, the above arguments apply, and a link is created between the node 1 and the node i + 1. The last LC operation on the node n − 1 creates a link between the nodes 1 and n. Hence the proof. Figure G2. Examples of C 1 and C 2 paths. A simple, connected, and undirected graph of 8 nodes. For demonstration, we consider the path L = [a ≡ 1, 2, 3, 4, 5, 6, 7 ≡ b], which is a C 2 path due to the existence of the links (1, 3), (1,5), (2,4), (3,7), and (4,6). From these links and the links constructing the path L, three C 1 paths of different lengths l = 2, 3, 4 can be distilled, given by L Theorem 1 suggests that the creation of a link between a and b in a connected graph is ensured if one chooses only C 1 paths between a and b and performs LC operations on the qubits in those paths. However, there usually exists a large number of C 2 paths in a connected graph, and there is no intuition whatsoever for choosing a C 1 path in order to obtain a tight MLB of LE. The challenge therefore is to develop an effective algorithm that can guarantee the creation of the link (a, b) irrespective of whether the chosen path belongs to C 1 or C 2 , when LC operations only on the nodes belonging to the chosen path are performed. Towards this goal, note that by virtue of the definition of C 2 , for a path L ∈ C 2 , one can distil a number of C 1 paths of different lengths from the links constructing L and the links connecting the nodes in the L, but not contributing to L (see figure G2 for examples). These C 1 paths have support on the nodes ∈ L only, and therefore the condition of performing LC operations only on the nodes belonging to L is satisfied. Theorem 1 ensures that each of these C 1 paths can be chosen to create (a, b) with certainty. Therefore, the challenge is to find a C 1 path from the C 2 path, which, in turn, leads to the creation of the link (a, b) according to Theorem 1. It is easily understood how this is achieved via the C 2 → C 1 algorithm represented by the pseudocode in appendix H.4.
Note here that to minimize effort, one has to use the distilled C 1 path having the shortest length, which may or may not be the shortest path between a and b. However, there is no intuitive reason behind using the shortest path, or the shortest distilled path to create the link, since the value of the MLB computed via the graph-based method in section 2.2 depends explicitly on the structure of the graph, and how the size of the neighborhood of the qubit-pair {a, b} in the transformed graph depends on the length of the chosen path between a and b in the original graph. In a graph with high connectivity, there may exist other paths having length close to the shortest length which may provide a tighter value of the MLB compared to the same obtained from the shortest path (we elaborate more on this in section 4). The C 2 → C 1 algorithm focuses on finding any one of the distilled C 1 paths that creates the link (a, b). For example, in figure G2, the C 2 → C 1 algorithm distils the C 1 path L (3) 2 of length 3, while the shortest distilled C 1 path is L (2) 2 of length 2, which is also one of the shortest paths between the qubits 1 and 7.
We now combine all these insights, and develop an LC-based graph-transformation algorithm that creates a link (a, b) by performing LC operations on a subset of the set of nodes situated on a chosen path, where the subset of nodes is chosen by adapting the chosen path according to its category. We call this algorithm the adaptive LC (ALC) algorithm, which has been represented by the pseudocode in appendix H.5. The LC operation on a simple, connected, and undirected graph w.r.t. a chosen node, the distillation of a C 1 path from a C 2 path via the C 2 → C 1 algorithm, and the creation of a link between any two given nodes in a graph by using ALC algorithm have been realized in the form of the Python open-source package ALCPack [115].
A word on the dependence of the run-time of the adaptive LC algorithm on the system size N is in order here. Since the algorithm takes into account the transformed graphs at each of its steps, it is difficult to determine an exact dependence of the run-time of the algorithm on system size. However, one can determine a bound on how the run-time scales with N. It is easy to see that the maximum size of the neighborhood N i of a node i can be N − 1, and it can host at most N−1 2 links if G is connected and undirected. Therefore, the maximum number of links that can be created or deleted during an LC operation on a single node is N−1 2 . Since there can be at most N − 2 nodes on a path L between a and b, the total number of link operations during the ALC algorithm is (N − 2) N−1 2 N 3 , which indicates a polynomial scaling with system size.
connected, and undirected star-shaped graph G P with the control qubit c as the central qubit via local unitary transformations H t 1 H t 2 H t 3 on the target qubits [92], where H represents a Hadamard operator (see figure 5). Note here that the convention of creating the graph state |ψ G P involves initializing all the qubits in |+ , and subsequently applying a controlled phase gate U z (c,t i ) on all (c, t i ) pairs of nodes [92], where the controlled phase gate, U z (c,t i ) = 1 2 (I c + Z c )I t i + (I c − Z c )Z t i , ( I 3 ) is connected to U x (c,t i ) via a Hadamard operation on the target qubit: U z (c,t i ) = H t i U x (c,t i ) H t i . In terms of the stabilizer operators, the plaquette stabilizer S x P = X c X t 1 X t 2 X t 3 is transformed to the graph-state generator X c Z t 1 Z t 2 Z t 3 via application of Hadamard operations on the target qubits. or equivalently, rank(A) = N. Moreover, an invertible recombination of the stabilizers operators leads to the same stabilizer operator. In the binary picture this recombination is represented by an invertible N × N binary matrix R that multiplies A from the right. Finally, the commutation of the stabilizer operators defining a stabilizer state is guaranteed by the relation A T DA = 0, where D is the 2N × 2N matrix: (J7) Clifford operation on stabilizer state. Clifford operations performed on multiple qubits of the stabilizer state are represented by 2N × 2N binary matrices Q that multiply A from the left. Clifford operations can not map a stabilizer operator or a recombination of them into the identity operator, implying that multiplication of Q to AR for all R can not result in vanishing a column. Consequently, Q is a full-rank matrix. For local operations like H or U Z in equation (J2), Q is the tensor product of the single-qubit matrices. For example, represents a π/2z-rotation applied on every qubit i for which Λ ii = 1, Λ being an N × N diagonal matrix.