A solvable model for graph state decoherence dynamics

We present an exactly solvable toy model for the continuous dissipative dynamics of permutation-invariant graph states of $N$ qubits. Such states are locally equivalent to an $N$-qubit Greenberger-Horne-Zeilinger (GHZ) state, a fundamental resource in many quantum information processing setups. We focus on the time evolution of the state governed by a Lindblad master equation with the three standard single-qubit jump operators, the Hamiltonian part being set to zero. Deriving analytic expressions for the expectation values of observables expanded in the Pauli basis at all times, we analyze the nontrivial intermediate-time dynamics. Using a numerical solver based on matrix product operators, we simulate the time evolution for systems with up to 64 qubits and verify a numerically exact agreement with the analytical results. We find that the evolution of the operator space entanglement entropy of a bipartition of the system manifests a plateau whose duration increases logarithmically with the number of qubits, whereas all Pauli-operator products have expectation values decaying at most in constant time.


Introduction
Graph states were introduced by Briegel and Raussendorf in 2001 [1] as special entangled states of N qubits.These states with multipartite entanglement play an important role in quantum information theory because they can be employed as a resource in a measurementbased quantum computation [2,3], and they can be used in error correction codes [4] and for quantum communications [5].In particular, permutation-invariant graph states [6], which are locally equivalent to an N -qubit Greenberger-Horne-Zeilinger (GHZ) state, are the subject of extensive research [7][8][9][10] and their creation and characterization serve as one of a few standard benchmarks of quantum computation hardware [11][12][13][14][15][16][17][18].The use of graph states for information processing in current quantum devices will inevitably have to face uncontrolled decoherence processes and some aspects of graph state entanglement under the presence of decoherence have already been investigated [19][20][21][22].Most of these previous works focused on discrete evolutions of the density matrix via completely positive maps (noisy channels).
In this work, we introduce and discuss a model of a graph state realized with qubits (or spin-one-half particles), which decoheres continuously in time as described by a Lindblad master equation for the density operator [23,24].We account for the three most prevalent local jump operators (dissipators); Two jump terms describe the incoherent transitions from |0〉 to |1〉 (and vice versa), and the third one is the so-called dephasing term.The initial state is a pure (graph) state, and it evolves into a mixed state under the action of the dissipation.Although it is a many-body problem, the structure of the model is simple enough that the expectation values of any observable can here be computed exactly by solving the equations of motion for the expectation values of products of Pauli matrices.
We complement our analytic treatment with the use of a numerical Lindblad solver [25,26], which is internally based on the C++ ITensor library [27] (see also [28] for a review on available numerical methods for this type of problem).In the solver, the state of the systema many-body density matrix ρ -is stored in the form a matrix-product operator (MPO).Since ρ is in general a matrix of size 2 N × 2 N , a brute-force numerical simulation of the Lindblad dynamics generally becomes very demanding beyond a dozen of qubits.Taking advantage of the fact that the states produced along the time evolution are only mildly correlated in the present model, the MPO approach allows us to reach a very high accuracy with modest computing resources (i.e.low MPO bond dimension) even with as many as N = 64 qubits.Other numerical approaches would also be efficient in the context of the current setup [6,29].
The presented model can be viewed as a toy model illustrating the basic mechanisms at play and allowing one to gain an understanding of the dominant dynamical behavior in similar setups.It may also be used as a starting point for more realistic studies with different graph states and some Hamiltonian terms competing with the dissipators.With our numerical approach, we are able to calculate global quantities that are not immediately accessible analytically, and observe an interesting scaling with size of bipartite correlations in the system.
In the next section, we present the model.In Sec. 3, we show how to compute exactly the evolution of the expectation value of any product of Pauli matrices in this model.From this result, one can then obtain the expectation value of any observable as a function of time.In Sec. 4, we then compare these analytical results with numerical simulations based on a matrix product operator (MPO) representation of the density matrix.Sec. 5 provides some concluding remarks.

Notations and definition of the model
We consider a system composed of N qubits (with basis states |0〉 = | ↑〉 and |1〉 = | ↓〉), indexed with i = 1 . . .N .We denote the usual Pauli operators by σ x , σ y and σ z , or alternatively by X , Y and Z.Additionally, we define σ ± = 1 2 (σ x ± iσ y ).For completeness, we start by recalling the definition of a graph state.A graph state is an entangled state that can be produced by the symmetrical two-qubit gate "controlled-Z" (CZ), which is defined by its matrix in the basis {|00|〉, |10|〉, |01|〉, |11|〉}.Given an undirected graph G(V, E) where V = 1 . . .N represents the qubits and E ⊂ V × V is the set of edges, the corresponding graph state |g〉 is defined by where |+〉 = 1 2 (|0〉 + |1〉).It is well known [1], that the graph state |g〉 is characterized by its stabilizers through In the present study, we focus on the case where the system is invariant under any permutation of the N qubits.We start at t = 0 from the only non-trivial fully symmetrical graph state, that is the graph state associated to the complete graph.A complete graph is a graph where all possible edges are present, so that each vertex is linked to the N − 1 other vertices.Our initial state (at t = 0) is thus given by As a side remark, we note that, thanks to the transformation rules of graph states under the action of local Clifford (LC) gates [20], the complete graph is LC-equivalent to the star graph. 1  In turn, the star graph state can be transformed into the N -qubit Greenberger-Horne-Zeilinger (GHZ) state [30] by application of Hadamard gates to all qubits except the center of the star.The complete graph is thus LC-equivalent to the GHZ state.We consider a time evolution generated by a Lindblad equation (see for example [31]) where the Hamiltonian part is set to zero.The state of the system is described by its density matrix ρ whose time evolution is given by where D is the dissipator, a linear superoperator acting on ρ.We consider three possible terms 1 The star graph has a central vertex connected to all the other vertices.

Closed-form observable dynamics
In this section, we derive some exact results for the time evolution of a large class of observables.

Observable dynamics
To compute the evolution of the mean value of a given (time-independent) observable O, we start from the fact that 〈O〉 = Tr(ρO) and we use Eq. 6 to get First let us consider D 0 where the superoperator Λ i 0 is given by Similarly, we have with It is clear that if O does not operate on qubit i then Λ i * [O] = 0.Moreover, for an operator O i acting on qubit i only, Λ i * [O i ] commutes with operators acting on the other qubits.The nontrivial action of the Λ i * can then be summarized by the following relations: where the qubit index of the Pauli operators are identical in the l.h.s and r.h.s and have been omitted for brevity.

Observables at t = 0
To lighten the notations, we will write X i instead of σ x i and likewise for Y and Z.As our system is invariant under any permutation of the qubits, specific indices are irrelevant.More generally, we note that which is independent of the actual order of the operators or the specific indices as long as they are all different.When indices are necessary, we will add them, for example To compute the expectation value of a product of Pauli operators at t = 0, that is on the complete graph state |g〉, we start with two remarks.First the stabilizer S = X Z N −1 leaves |g〉 unchanged (see Eqs. 3 and 4) so that Second, since Z commutes with C Z and since C Z 2 = 1, we have for n > 0 We start by computing observables of the form 〈X n Z l 〉.To do this, we introduce the stabilizer at one of the indices of the X .So for n > 0 To be real, this last expression above must be 0 when n is even.We can do the same for 〈Y m Z l 〉, which gives for m > 0 And again this must be 0 if m is odd.Substituting Eq. 19 in Eq. 20, we can conclude that for even m which is 0 if l > 0 and 1 otherwise.We can also finish the computation for X for odd n which is 1 if n + l = N , and zero otherwise.The last product we have not yet computed is the general one 〈X n Y m Z l 〉 with n > 0 and m > 0.
This can be nonzero only if n + m is odd.But if n = 1 there are no Y left and this is zero according to Eq. 22.And if n > 1, the right-hand side can be nonzero only if n + m − 1 is odd (thus n+ m is even) which we have already excluded.Thus, To summarize, only the following products of Pauli operators have nonzero mean values in the complete graph state: and

Solution of the equations of motion
We start with an example to show how the equation of motion leads to some differential equations.Here we consider 〈X Z〉 in the case where only g 0 is not zero.
Now the general formula for Globally, we obtain for l > 0 where Here, α is the rate of dephasing (decoherence) associated to X and Y (with its inverse equal to the characteristic T 2 timescale), β is the decay parameter associated to Z -the inverse of the relaxation time T 1 , γ is the global drive towards the steady state, and γ/β determines the thermal steady state population (value of 〈Z〉) reached in the limit of large time. 3In the case where γ = 0 (that is g 0 = g 1 ), all observables have a simple exponential decay.These equations can be solved by recursion starting at l = 0 using the initial conditions of the previous section.Indeed, at l = 0 (that is no Z), we get and all the others products without Z give zero because they start at 0 and stay there.Now we can increase l and get Note that if γ = 0 or γ = β = 0 then 〈Y 2n Z l 〉 = 0 for l > 0. Finally for 〈X 2n+1 Z N −2n−1 〉, the equations are directly solved and we obtain with all the others being identically zero.
It is interesting to note that the stabilizers that characterize the initial graph state decay very rapidly, i.e. with a timescale inversely proportional to the size of the system.This reflects some relative fragility of the graph state correlations in the present dissipative context, and may be related to the extensive number of neighbors in the initial complete graph.

Reduced density matrices
The two-qubit density matrix can be obtained from the two-point correlations computed previously.Writing z ± = 1 ± 〈Z〉 = 1 ± γ β 1 − e −β t (see Eq. 31 with n = 0 and l = 1) and y 2 = 〈Y Y 〉 = e −2αt (see Eq. 30 with n = 1), this matrix reads (for N > 2 only) Likewise for 3 qubits and N > 3, we have More generally, reduced density matrices for larger subsystem can be obtained by noting that each matrix element is the expectation value of a product of N operators which are of the type (Z i +1)/2 (if the matrix element connects two states where Z i = 1), (1− Z i )/2 (matrix element between two states where Z i = −1), or σ ± i (matrix element between two states with opposite Z i ).
Since ρ 2 can be decomposed as it is easy to see that it is separable (at all times).We also checked that ρ 3 and ρ 4 (not shown) show no negativity, independent of the time t.
It is in fact a well known property of GHZ states [20] that the reduced density matrices are separable, even though the system has some global multipartite entanglement.In the present model, the reduced density matrices are thus all separable at t = 0.In addition, since the dynamics is only due to single qubit disspators the reduced density matrices must remain separable at all times (unless one considers the system globally (all qubits)).

Computational results
The possibility to simulate quantum open systems by representing the density matrix as a (onedimensional) tensor network was originally discussed in Refs.[32,33].In the present work, we use a representation of ρ as an MPO, closely related to what was done by Prosen and Žnidarič in [34].In our implementation, the density matrix ρ is encoded in a vectorized form (denoted by |ρ〉〉) as an ITensor matrix-product state (MPS) over an Hilbert space with dimension 4 at each site. 4More details on this representation can be found in App. A. Concerning the implementation of the time evolution, the first step amounts to expressing the Lindbladian superoperator L as a (super-) MPO that can act on a vectorized |ρ〉〉.Next, one constructs a (super-) MPO approximation of the exponential exp (τL) associated to a small time step τ.This approximation is based on an extension of the W II scheme proposed in Ref. [35], and it allows for long-range interactions.In the scheme we use (see appendix A of Ref. [36]), the Trotter error scales as O(τ 5 ).All the simulations used a time step τ = 0.004.The truncation in the MPO was carried out with a maximum discarded weight of 10 −16 .
As a pure state, the initial state can be written exactly has a matrix-product state with bond dimension equal to 2, and as a density matrix ρ, it can be represented exactly by an MPO of bond dimension equal to 4. It turns out that the local dissipation terms of the model do not increase the required bond dimension.Note however that due to small numerical errors the bond dimension was sometimes observed to be above 4 in the simulations (but never exceeding 15).

Dissipative dynamics of observables
In this section, we present the dynamics of Pauli observables calculated from the analytic expressions of Eqs.30-32, together with numerical simulation results from the MPO solver.App.A gives more details on the numerical simulations.To explore the parameter space, we studied five representative cases varying the values of g 0 , g 1 and g 2 .The values used are shown in Tab. 1, together with some characterization of the environments that would produce such parameter values.Table 1: Different sets of physical parameters used in the simulations.The dissipation parameters g 0 , g 1 and g 2 are defined in Eqs.7-9.The parameters α, β and γ are linear combinations of the g i defined by Eq. 29.First, we look at Eq. 30 in the left panel of Fig. 1.The numerical results are in perfect agreement with the theory.Comparing the numerical data with the exact expression shows that the absolute value of the error is always below 10 −9 .In the following figures, the maximum error also never exceeds this value.Since the expectation value of a product of an odd number of Y vanishes, 〈Y Y 〉 is a connected correlation, and it shows a decay with a characteristic timescale ∼ 1/α.Now we turn to Eq.31, first in a simple case for the single qubit observable 〈Z〉.The result is displayed in the right panel of Fig. 1 (this corresponds to n = 0 and l = 1).The cases 2 and 5 are not shown since they have γ = 0 and thus 〈Z〉 = 0. 〈Z〉 displays a relaxation toward the steady state value 〈Z〉 t→∞ = γ β = g 0 −g 1 g 0 +g 1 .When g 1 = 0 this is simply a relaxation toward the |0〉 state.We also note that the effects of the correlations in the system are not visible in this observable, in the sense that the exact same behavior would be observed independently of the initial state provided that 〈Z〉 = 0 on all qubits at time 0.
In the more complex case where n ̸ = 0, we cannot scale all the cases on one curve, so we chose to show the dependence in the number of qubits for one case.In the left panel of Fig. 2, we show the time evolution of 〈Y Y Z〉 (that is n = 1 and l = 1) for case 1 and different number of qubits.Again the simulations are in perfect agreement with the theory.This observable has a non-monotonous time evolution for a simple reason: from Eq. 31, we see that this threepoint observable factorizes into 〈Y Y Z〉 = 〈Y Y 〉〈Z〉, that is a product of a decreasing function  Finally, we also checked the expectation value of the stabilizer of the complete graph state, that is Eq.32 with n = 1.The results are displayed in the right panel of Fig. 2, they are again in perfect agreement with the theory.The initial value is 1, as it should be since the initial state is an eigenstate of the stabilizer for the eigenvalue 1.We then observe an exponential decay with a characteristic timescale given by (α + (N − 1)β) −1 and which decreases with the number of qubits.This size-dependence of the decay rate can be viewed as a consequence of the fact that this specific observable involves all the qubits and reflects some global correlations in the system.

Operator space entanglement entropy
The method to compute the von Neumann entanglement entropy associated to a given bipartition of a graph state is explained in Ref. [4].In the case of the complete graph state, the result is S vN = ln 2 independent of the the bipartition (for two non-empty subsystems).This result is also easy to obtain using the fact that the complete graph state is LC-equivalent to the GHZ state.
For a mixed state ρ, it is interesting to consider the operator space entanglement entropy (OSEE) [37], a quantity that naturally arises in simulations of the density matrix dynamics represented using MPO.In a closely related context, entanglement entropies associated to (unitary) operators were also discussed in Ref. [38].The OSEE has been the subject of many theoretical and numerical studies in the context of many-body problems.For instance, Ref. [39] considered the OSEE in thermal density matrices of spin chain models, Ref. [40] studied the OSEE in one dimensional many-body systems using (conformal) field theory techniques and [41] studied the long-time dynamics of OSEE in a spin chain subject to a dephasing dissipation (equivalent to the g 2 term in the present study).Another recent work [42] compared two different approaches for the dynamics of a dissipative quantum spin chain, Lindblad master equation versus quantum trajectories, and provided a comparison between the OSEE in the Lindblad approach with the entanglement entropy in the trajectory formalism.
The OSEE can be defined by considering the vectorization |ρ〉〉 of ρ, which is a pure state in an enlarged Hilbert of dimension 4 per site (spanned by the 3 Pauli matrices plus the identity matrix).The OSEE associated to a given bipartition into two subsystems A and B is by definition the von Neumann entanglement entropy OSEE A|B = S vN,|ρ〉〉 associated to this partition of the vectorized pure state |ρ〉〉 (more details in App.A).The OSEE quantifies the total amount of correlations between the two subsystems.We note that the OSEE alone does not indicate whether the correlations are mostly classical or quantum.
For a pure state |g〉, the associated density matrix is ρ = |g〉〈g| and the OSEE of ρ for a given bipartition is by construction twice the von Neumann entropy associated to the same bipartition in |g〉.So, in our model, the OSEE at time t = 0 is 2 ln 2 for any nontrivial bipartition of the N qubits.At long times, the system reaches a product state (if g 1 = 0 it simply corresponds to all the qubits in state |0〉) which is a state with bond dimension equal to 1 and vanishing OSEE.At intermediate finite times t > 0, the OSEE must therefore decrease from its initial value and eventually converge to 0 for any bipartition of the system.It is, however, no longer independent of the bipartition and in the rest of the paper, we will focus of the bipartition in two subsystems of equal size (N /2).
The dynamics of the OSEE for this bipartition in two equal halves of system is shown in Fig. 3.The interesting feature is the appearance of a very clear plateau at OSEE = ln 2, whose duration grows with the number of qubits.To explore this behavior, we determined numerically the scaling laws for both the time at which the plateau begins and the time at which it ends.We observe a 1/N behavior at early times and ln N behavior for the time at the end of the OSEE plateau.The corresponding plots are shown in Fig. 4. At early times, the behavior is similar for cases 3, 4 and 5 with a 1/N behavior (left panel Fig. 4) and a plateau at OSEE = ln 2. Case 2 is different with a plateau that starts at t = 0, see Fig. 5.
This behavior arises only in the case with β = 0, which corresponds to an absence of flip terms (pure dephasing).At late times, the ln N behavior is universal, with different values of the prefactor δ (up to some finite size effects the OSEE depends only on t − δ ln N at late times).In case 2, the data collapse is particularly striking (see right panel of Fig. 5).
Comparing the values of δ to the parameters in each case, we see that for all our cases the value of δ is compatible with δ = 1/(2α).It is remarkable that the OSEE stays essentially Figure 4: Time evolution of the OSEE between the two halves of the system for different values of N in case 1 (same data as Fig. 3).Left: time rescaled by the system size N .The collapse of the curves at early times shows that the initial drop of the OSEE takes place over a timescale proportional to 1/N .Right: same data with time shifted by δ ln N , here δ = 1.In this panel the collapse of the curves at the end of the plateau illustrates the fact that the duration of the plateau is proportional to δ ln N .constant for a duration that increases with the number of qubits, although the correlations decrease exponentially at best in constant time.A large enough system can be in a regime where 1/β ≪ t and t ≲ δ ln(N ).In case 1 (g 0 only), the first condition implies that 〈Z〉 is arbitrary close to 1, while the second condition puts the system in the OSEE plateau, with OSEE ≃ ln (2).In other words, the system can have an arbitrary low density of qubits in the |1〉 state and, still, some sizeable correlation between the two halves of the system.
We checked that such plateau is absent if the initial state is a graph state with a lower connectivity.As an example, Fig. 6 shows the OSEE in a case where the initial state is a graph state constructed from a periodic one-dimensional lattice with N sites (a ring).For a subsystem A of the form in such a ring graph state [4], hence the value OSEE = 4 ln(2) at time t = 0.The correlations are plausibly only short-ranged (with a finite correlation length) in that case, so that the correlations between the two halves of the bipartition essentially come from the qubits close to the boundaries/cutting between the two halves.Thus, when the system size becomes significantly larger than the correlation length, the OSEE becomes independent of N .The OSEE then drops to zero over a characteristic timescale which is independent of the system size, contrary to the cases where the initial state is a complete graph.Finally, we also looked at what happens with a simple nonzero Hamiltonian.As an example we considered an Ising interaction.Namely, we chose H = Z a Z b where a and b refer to two fixed qubits.When measuring the OSEE, one subsystem contained qubit a and the other qubit b.With these parameters, the maximum bond dimension was 16.The results for case 1 are shown on Fig. 7.These figures are to be compared with Fig. 3 and 4 (right panel) where H was zero.There are now four regimes: as before, at early times, a fast initial drop of the OSEE over a timescale proportional to 1/N (scaling not shown), at constant time (around t ∼ 0.8) a peak, then a quasi plateau up to time δ ln N and finally a decay to 0 in constant time (δ = 1 in this case, as with H = 0).The peak is associated to a local OSEE maximum at time t ∼ 0.8.We interpret this peak as a consequence of some new correlations created by the Ising term.These correlations are then washed out by the dissipation.We then have a regime of slow decay of the OSEE.In this regime the slope of the OSEE appears to scale as 1/ ln N at large N .As far as the OSEE is concerned, in this regime H plays here the role of a perturbation to the plateau observed when H = 0, and the effect of this perturbation decreases when N becomes large.Due to the observed finite-size scaling we expect a plateau to asymptotically appear in the large N limit.However, to clearly separate the plateau from the peak one would need to go to much larger sizes, which is out of reach at the moment.We obtain similar results in the other cases (data not shown).

Conclusion
In this paper, we have introduced an exactly solvable toy model for the decoherence of a graph state.We have considered the time evolution of the complete graph state under a Lindblad equation with general single-qubit dissipators and no Hamiltonian.Exploiting the permutation invariance of the model has allowed obtaining simple closed formulae for the expectation values of any product of Pauli operators at any time.The method can be extended to other types of graph states, although the number of equations will grow for less symmetric situa- tions.The availability of analytic solutions is valuable as a guiding tool in understanding the dissipative mechanisms acting in numerical studies of complex setups.
We have compared the theoretical results with a numerical solver that was pushed up to 64 qubits.The results are in perfect agreement with the theory, showing that the MPO approach is adapted for simulating this type of problem.
A peculiar long-lasting plateau has been identified in the OSEE between the two halves of the system, pointing to the presence of some nonlocal long-lived correlations in this setup.Despite the dissipative processes acting everywhere on the qubits, the correlations have been observed to survive for a time that increases with the system size.The t ∼ ln(N ) scaling of the OSEE decay survives with a simple nonzero Hamiltonian.In a future study, we consider it valuable to compute exactly the OSEE for t > 0 in this model.[43].Generally speaking, for a given precision the higher the von Neumann entropy accross a certain cut, the larger the bond dimension of the MPS along that cut.See [44] for a precise discussion concerning the Schmidt spectrum, Rényi entropies, and the "simulability" using MPS.

) 3
In the long time limit and for n = 0 and l = 1 the Eq. 31 below gives 〈Z〉 t≫β = γ β which is the thermal equilibrium value if we introduce a local Hamiltonian H = ∆ i Z i and a temperature T satisfying tanh(∆/(k B T )) = γ β .

Figure 1 :Figure 2 :
Figure 1: Left: Time evolution of 〈Y Y 〉 in the different parameter cases with N = 64 qubits.The rescaled time α • t in the horizontal axis allows the collapse of the curves associated to different sets of parameters.The line corresponds to Eq. 30 in the case n = 1.For this quantity, the maximum deviation between the numerical result and the exact value is less than 10 −11 .Right: same for 〈Z〉.For this quantity, the relevant rescaling of the time is β • t.The line corresponds to Eq. 31 in the case n = 0 and l = 1.Cases 2 and 5 are not shown as they have γ = 0 and thus 〈Z〉 = 0.The maximum deviation is here less than 10 −9 .

Figure 3 :
Figure3: Time evolution of the OSEE between the two halves of the system for different number of qubits in case 1.A plateau in the OSEE value at intermediate times t, whose duration grows (logarithmically) with the number of qubits is clear.

Figure 5 :
Figure 5: Time evolution of the OSEE between the two halves of the system for different number of qubits in case 2. Left: data plotted as a function of time.Right: data plotted as a function of time shifted by δ ln N , here δ = 0.25.All the curves are essentially on top of one another and are indistinguishable at the scale of the figure.

Figure 6 :
Figure6: Time evolution of the OSEE between the two halves of the system, for a simulation of a ring graph state (initialized at t = 0), with the dissipator D 0 .The data for N = 16, 32 and 64 are almost on top of each other, showing rapid convergence to a limiting curve in the large N limit.Contrary to the cases where the initial state is a complete graph, these curve do not exhibit any plateau.In this simulation the MPO bond dimension is 16 or below.

Figure 7 :
Figure7: Left: Time evolution of the OSEE between the two halves of the system for different values of N in case 1 with a nonzero Hamiltonian (a ZZ interaction between two qubits one in each half of the system).Right: same data with time shifted by δ ln N , here δ = 1.

5
If {|a〉} and {|b〉} represent respectively orthonormal bases of the Hilbert spaces of the regions A and B, we can decompose |ψ〉 as |ψ〉 = a,b ψ a,b |a〉 ⊗ |b〉.This can in turn be rewritten |ψ〉 = k λ k |φ A k 〉 ⊗ |φ B k 〉 where the λ k can be taken real and positive and are the singular values of the matrix ψ a,b and the vectors |φ A k 〉 and |φB k 〉 satisfy 〈φ A k |φ A k ′ 〉 = δ k,k ′ = 〈φ B k |φ B k ′ 〉.The von Neuman entropy associated to this bipartition has a simple expression in terms of the Schmidt singular values: S vN = − k p k ln(p k ) with p k = (λ k ) 2 and k p k = 1 = 〈ψ|ψ〉.A central idea in the construction of an MPS approximation of |ψ〉 is to truncate this Schmidt spectrum by removing the singular values below a certain threshold