Manual

• Manual
• • How to read this manual
  • Sites and Systems
  • States
  • • Prepared states
  • Conserved quantities
  • Limits
  • Operators
  • • Fermions
  • Algorithms
  • • Gates
    • Ground states
    • Time evolution
    • Noise
    • Steady states
    • Thermal states
  • Measurements
  • Accuracy

How to read this manual

This manual presents the objects of TMS one after the other, sites, states, operators, then the algorithms and the measurements, with small examples run directly from Julia, and ends with how to check the accuracy of the results. The pages that follow build on it: High Level Interface shows how to write a whole simulation as a list of phases, which writes its results to files and can be stopped and resumed, the way most programs are written; Conserving a quantity tells how the conservation of a quantity is declared, what it forbids and what it saves; Threads and performance is for when a computation gets large; and Extending TMS is for sites, phases, algorithms and representations of one's own. The reference pages give every function in detail.

The examples assume TMS has been imported:

using TensorMixedStates

Sites and Systems

The first step in using TMS is the definition of your quantum system. In TMS, a system is composed of a finite number of sites numbered from 1 (1, 2, ..., N). These sites may be all identical or not.

There are eight different predefined types of site: Qubit, Qudit, Fermion, Boson, Spin, Electron, Tj and Qboson.

To use each of these sites and the corresponding predefined operators you need first to import the corresponding module. For example to use qubits, you need to write

using .Qubits

Note the "." before the name and the "s" at the end.

To define a site, call the corresponding constructor, for example

s = Qubit()

Qubit()

Four site constructors need an argument: Qudit(dim) and Boson(dim) for the dimension of the local Hilbert space, Spin(s) for the spin, and Qboson(q, dim) for the deformation parameter and the dimension. For example

using .Bosons, .Spins

s = Boson(4)

Boson(4)

s = Spin(3/2)

Spin(1.5)

You can now define a quantum system by declaring the sites it contains:

system1 = System(10, Qubit())

gives you a system with 10 qubits. Systems may have sites of different types, in which case System is given an array of sites

using .Fermions

system2 = System([Qubit(), Boson(4), Fermion()])

gives you a three site system.

States

A state is either a wave function, its pure representation, written Pure, or a density matrix, its mixed representation, written Mixed, which open systems need. Both are stored as a matrix product state (MPS): a chain of tensors, one per site, joined by links whose size, the bond dimension, measures how much entanglement, or for a density matrix how many correlations, the state can hold. A density matrix is stored as the vector of its entries, so each site has dimension $d^2$ instead of $d$, its norm is the Hilbert-Schmidt norm $\sqrt{\mathrm{tr}\,\rho^\dagger\rho}$ rather than its trace, and it costs more than a wave function of the same system.

To create a state, we call the State constructor

state1 = State{Pure}(system1, "Up")

returns a pure up state in a 10 qubit system. Predefined local states are designated by their name. Here "Up" is a predefined state of site Qubit.

Sites need not be all in the same local state, in which case we give State an array of local states

state2 = State{Mixed}(system2, ["+", "2", "Occ"])

Here we choose a mixed representation.

States may be added or multiplied by a number (they need to be based on the same system). For example

ghz = (State{Pure}(system1, "Up") + State{Pure}(system1, "Dn")) / sqrt(2)

We can transform a pure representation into a mixed representation by

mixedstate = mix(state1)

For mixed states there is a local mixed state "FullyMixed" which corresponds to a density matrix proportional to the identity matrix (that is the infinite temperature state).

If you need a local state which is not predefined, it is possible to pass its vector (or matrix for mixed states) directly. For example, we could also define state1 by

state1 = State{Pure}(system1, [1., 0.])

Prepared states

States that are not products of local states, but have an exact MPS of small bond dimension, are built by functions of their own, which take the local states they are made of:

myneel = superposition(System(4, Qubit()), [1 => ["Up", "Dn", "Up", "Dn"], 1 => ["Dn", "Up", "Dn", "Up"]])
myghz = ghz_state(System(6, Qubit()), "Up", "Dn")          # (|↑↑…↑⟩ + |↓↓…↓⟩)/√2
mydicke = dicke_state(System(6, Qubit()), 2, "Up", "Dn")   # two sites down, in every way
mydimers = dimer_state(System(6, Spin(1/2)), [(1, 2), (3, 4), (5, 6)], "1/2", "-1/2")
myaklt = aklt_state(System(6, Spin(1)))
mysector = fully_mixed(System(6, Fermion()), N => 3)       # 3 fermions at infinite temperature
measure(mydicke, Z)

6-element Vector{Float64}:
 0.3333333333333333
 0.3333333333333333
 0.33333333333333337
 0.33333333333333337
 0.3333333333333333
 0.3333333333333333

superposition superposes product states with the coefficients given, here the two Néel states, as mixture mixes them with weights in a mixed representation, w_state is the Dicke state of one site, dimer_state puts a singlet on each pair of sites given, aklt_state is the ground state of the AKLT chain of spins one, and fully_mixed the fully mixed state of a sector, here of 3 fermions. The Slater determinants and the Fermi seas of free fermions are built by slater_state and fermi_sea, see Fermions, and the thermal states by thermalize, see Thermal states. Each is a State, which CreateState takes as its state, see Prepared states for the details.

Conserved quantities

A site can be told that a quantity is conserved, for instance the number of fermions:

mysystem = System(10, Fermion(conserve = N))
mystate = State{Pure}(mysystem, [isodd(i) ? "Occ" : "Emp" for i in 1:10])

The tensors are then block sparse, which makes a large computation smaller and faster, and in exchange the state stays in the sector it was built in: a local state such as [1, 1] / √2, which superposes two numbers of particles, is refused, and so is an operator of no definite charge acting on the state, which can still be measured. For a mixed state, a quantity may also be conserved strongly, Fermion(conserve = strong(N)). All of this is described in Conserving a quantity.

Limits

TMS uses Matrix Product States to represent quantum states internally. It is important to control the parameters of this approximation, in particular the maximum bond dimension and the cutoff, the total weight of the singular values a truncation of the tensors may drop. To achieve this, many functions accept a Limits object as keyword argument containing those parameters. It is built thus

lim = Limits(cutoff = 1e-10, maxdim = 50)

Limits(cutoff = 1.0e-10, maxdim = 50)

Any of the arguments may be omitted: maxdim and mindim then impose nothing, and the default cutoff, 1e-26, only discards what rounding leaves, the singular values below about 1e-13 of the norm, while cutoff = 0 discards nothing. The truncation is a choice to be checked, see Checking the accuracy.

A third parameter, mindim, sets the bond dimension the truncation is not allowed to go below, as in Limits(cutoff = 1e-10, maxdim = 50, mindim = 4). maxdim takes precedence when the two conflict.

To apply the constraints on a state, one uses

newstate = truncate(ghz; limits = lim)

Many functions accept such an argument. For example, when adding states instead of

state3 = (state1 + ghz) / 2

one can write

state3 = +(state1, ghz; limits = lim) / 2

Operators

In TMS, there are two kinds of operators: generic operators and indexed operators. For example,

X

X

represents the $\sigma_x$ Pauli operator for qubits. This is a generic operator, it is not applied to a specific site.

X(3)

X(3)

represents the $\sigma_x$ Pauli operator applied to the system site number 3. This is an indexed operator.

Note that all predefined operator names start with a capital letter, so it is better to keep your own identifiers lowercase to prevent name collisions with them.

Lowercase is not a safe harbour by itself, though. using TensorMixedStates also brings in more than eighty lowercase names, among them output, measure, trace, dim, matrix, tensor, apply, norm and sample. Assigning to one of them at the top level of your program shadows the function for the rest of the file, and if you happen to have used it before assigning to it, Julia 1.10 and 1.11 refuse the assignment outright with cannot assign a value to imported variable. There is no such issue inside a function, where state = ... is an ordinary local variable. This manual prefixes its own variables with my, as in mystate and myop, which is one way of staying clear. Another is to take from a site module only the names one uses, using .Qubits: Qubit, X, Z, which leaves N free, or to write them qualified after import .Qubits as Q, as Q.X(3).

The operator system is very rich and flexible. For example, if you want to use this Hamiltonian

\[H = \sum_{i=1}^{n-1} \sigma_x(i) \sigma_x(i+1)\]

you will simply write

n = 10
h = sum(X(i)X(i+1) for i in 1:n-1)

X(1)*X(2)+X(2)*X(3)+X(3)*X(4)+X(4)*X(5)+X(5)*X(6)+X(6)*X(7)+X(7)*X(8)+X(8)*X(9)+X(9)*X(10)

Many operations are defined on generic operators:

• addition, multiplication and power by a number
• tensor product: X⊗X is a two site operator (⊗ is usually obtained by typing \otimes in your editor, just in case, one can also write tensor(X, X))
• dag represents the adjoint operator, for example C is the c operator for fermions and dag(C) is $c^\dagger$; C' is the same, and both take placed operators too, h + h'.
• Dissipator represents a Lindblad dissipator, for example Dissipator(Sp) is the dissipator whose jump operator, Sp, the $S^+$ operator, flips a qubit toward up
• Gate represents an operator to be applied as a gate on a mixed state. It is useful to define noisy gate operators, for example 0.9Gate(Id) + 0.1Gate(X) is a noisy gate operator that will apply a $\sigma_x$ gate 10 percent of the time.
• Proj represents an operator that projects on the given state, for example Proj("Up") projects qubits on the up state, or on an eigenspace, Proj(Sz => 0) projecting a spin 1 on the state of zero $S_z$.
• Level numbers the basis states of any site, $\mathrm{diag}(0, 1, \dots)$: its eigenbasis is the basis of the site, and Level => n the basis state n, counted from 0.
• the functions exp and sqrt: for example sqrt(Swap)
• controlled makes controlled gates, which do nothing when their first site is in its first basis state, Level => 0, as "Up" for a qubit: controlled(X) is the CX of the module Circuits, see Gates

For example one can define the Rxy 2-site operator by

Rxy(t) = exp(-im * t * (X⊗X + Y⊗Y) / 4)

Rxy (generic function with 1 method)

If this is not enough to define your favorite operator you can create new ones with named, from an expression, a matrix or a function of the sites, see Operators and measurements of one's own:

myop = named([1 1 ; 1 -1] / √2, "MyOp")

MyOp

Finally from generic operators, we define indexed operators by simply applying them to the corresponding sites

Rxy(0.2)(2, 5)

exp(-0.05im*X⊗X-0.05im*Y⊗Y)(2,5)

myop(3)

MyOp(3)

Fermions

A fermionic operator is written as it is: C(i) destroys a fermion on site i and dag(C)(j) creates one on site j, the Jordan-Wigner strings being inserted for you wherever the operator is used, in expectation values, in MPOs and in gates. A hopping term is thus simply

hopping = sum(dag(C)(i)C(i+1) + dag(C)(i+1)C(i) for i in 1:n-1)

The strings follow the order of the sites: C(j) is the matrix of C on site j with F on every site before it. A product of two fermionic operators placed in increasing order of the sites, A(i) * B(j) with i < j, is thus A * F on site i, F on every site in between and B on site j, and in the other order it takes the sign of the swap. A matrix given for an operator of several sites, see Operators and measurements of one's own, is read in the basis of the sites as it is, with no string.

A function of an even fermionic operator of several sites, as the exponential of a hopping term, exp(-0.1im * (dag(C) ⊗ C + dag(dag(C) ⊗ C)))(2, 5), is applied as a gate with its strings, on sites in any order and apart. A function of an odd operator, as exp(C), mixes the two parities: it is split into its even part and its odd part, which takes the strings.

The ground state of free fermions, the Fermi sea of a quadratic hamiltonian, is built at once by fermi_sea, here with three fermions, and any Slater determinant by slater_state, see Prepared states:

sea = fermi_sea(System(n, Fermion()), -hopping, 3)
measure(sea, N)

10-element Vector{Float64}:
 0.1714223902218396
 0.3817216995052644
 0.3351266180805347
 0.2687201853988544
 0.3430091067935101
 0.3430091067935089
 0.2687201853988528
 0.33512661808053745
 0.38172169950526236
 0.17142239022183706

Algorithms

We can now work with states and operators.

Gates

The gates of quantum circuits are in the module Circuits: H, S, T, Phase, the rotations, Rx(θ) being $e^{-i\theta X/2}$ on a qubit as in OpenQASM and Qiskit and Rxx(θ) $e^{-i\theta X \otimes X/2}$, and the controlled gates CX or CNOT, CY, CZ, and CCX or Toffoli:

using .Circuits

Gates are applied with apply. In a product of gates the rightmost acts first: here the Hadamard gate on the first qubit, then the controlled not, which turns two up qubits into a Bell pair

mybell = apply(CX(1, 2) * H(1), State{Pure}(System(2, Qubit()), "Up"))
measure(mybell, [Z(1)Z(2), X(1)X(2)])

2-element Vector{Pair{String, Any}}:
 "Z(1)*Z(2)" => 1.0
 "X(1)*X(2)" => 1.0

Applying all the gates in a single call is much more efficient than one by one, and the keyword argument limits constrains the truncations made along the way.

Ground states

The ground state of a hamiltonian is computed by ground_state, from a starting state, here a random state of bond dimension 8, by the method algo, Dmrg. It returns the energy and the ground state, which can then be measured. For the Ising chain in a transverse field

mysystem = System(10, Qubit())
hamiltonian = -sum(Z(i)Z(i + 1) for i in 1:9) - sum(X(i) for i in 1:10)
energy, ground = ground_state(hamiltonian, random_state(Pure, mysystem, 8);
                              algo = Dmrg(nsweeps = 10), limits = Limits(maxdim = 20))
energy

-12.381489999654798

measure(ground, X)

10-element Vector{Float64}:
 0.8512118670042681
 0.7324519047454634
 0.7021738068131279
 0.6900669649442622
 0.6853707301414419
 0.6853707301414402
 0.6900669649442619
 0.7021738068131284
 0.7324519047454645
 0.8512118670042708

nsweeps is the number of sweeps, passes of dmrg along the chain and back, Dmrg(tol = …) stops the search earlier once the energy changes less than tol between two sweeps, and limits may give one value per sweep, as in Limits(maxdim = [10, 20, 50]).

Given eigenstates found before, orthogonal_to, the search gives the lowest eigenstate orthogonal to them, an excited state: here the first excited state, and the gap

energy1, excited = ground_state(hamiltonian, random_state(Pure, mysystem, 8);
                                algo = Dmrg(nsweeps = 20), limits = Limits(maxdim = 20),
                                orthogonal_to = [ground])
energy1 - energy

0.29892037434574625

A penalty, weight times the projectors on the states given, is added to the hamiltonian for that. It must exceed the gap between the state searched and each of them, and its default, 1.1 times a bound of the width of the spectrum, always does but slows the convergence, which a weight just above the gap speeds up, see ground_state.

Time evolution

Time evolution is done by evolve, which takes an evolver, the time to evolve for and the state, and the method as algo, Tdvp() by default, or ApproxW. The evolver follows the convention

evolver = -im * hamiltonian + dissipators

dissipators being accepted on a mixed state only. On a pure state, it integrates $\frac{d}{dt}|\psi\rangle = -iH|\psi\rangle$, and on a mixed state the Lindblad equation

\[\frac{d\rho}{dt} = -i[H, \rho] + \sum_k \left(L_k \rho L_k^\dagger - \tfrac{1}{2}\{L_k^\dagger L_k, \rho\}\right)\]

each term Dissipator(L) giving one jump operator $L_k$, with $\hbar = 1$. A rate $\gamma$ is written γ * Dissipator(L), or equally Dissipator(sqrt(γ) * L), since Dissipator(c * L) is abs2(c) * Dissipator(L). Under $H = \sum_i \sigma_z^i$, qubits starting in "+" precess, with $\langle \sigma_x \rangle = \cos 2t$, which is $\cos 1 \approx 0.5403$ at $t = 0.5$:

myevolved = evolve(-im * sum(Z(i) for i in 1:4), 0.5, State{Pure}(System(4, Qubit()), "+");
                   nsteps = 10)
measure(myevolved, X)

4-element Vector{Float64}:
 0.5403023058681404
 0.54030230586814
 0.5403023058681382
 0.5403023058681393

nsteps is the number of steps, each of length t / nsteps, and limits constrains the state as before. With no hamiltonian and Dissipator(Sm) on each site, the qubits of a mixed state go from up to down, with $\langle \sigma_z \rangle = 2e^{-t} - 1$, which is $-0.2642$ at $t = 1$:

myrho = State{Mixed}(System(4, Qubit()), "Up")
mydecayed = evolve(sum(Dissipator(Sm)(i) for i in 1:4), 1.0, myrho; nsteps = 10)
measure(mydecayed, Z)

4-element Vector{Float64}:
 -0.26424111765711483
 -0.2642411176571154
 -0.26424111765711555
 -0.26424111765711494

ApproxW approximates the exponential of the evolver over a step by MPOs, built from the approximations WI or WII of Zaletel et al. It is given the order of its approximation, from 1 to 4, and w, 1 or 2 for WI or WII, for which order = 4, w = 2 is usually a good choice:

mydecayed = evolve(sum(Dissipator(Sm)(i) for i in 1:4), 1.0, myrho;
                   algo = ApproxW(order = 4, w = 2), nsteps = 10)
measure(mydecayed, Z)

4-element Vector{Float64}:
 -0.2642411176571155
 -0.26424111765711683
 -0.26424111765711644
 -0.2642411176571165

An evolver may also depend on time, see Time dependent evolvers.

Noise

Noise acts on a mixed state, either as a channel applied between gates, or as a term of an evolver. Three kinds are ready for sites of any type, each as a channel …_gate(p, …), taking place with probability p, and as a generator …_dissipator(γ, …), of rate γ:

• reset_channel(p, state) resets a site to state, see Reset; with n, reset_channel(p, n, state) resets n sites together, and reset_channel(p, [s1, s2]) each site to its state;
• depolarizing_channel(p, n) is the reset to "FullyMixed": on a qubit, the channel applying each Pauli operator with probability p/3 is depolarizing_channel(4p/3);
• dephasing_channel(p, A) erases the coherences between the eigenspaces of A, or in the basis of the site without A, see Dephase.

Evolving under a generator for a time t is the channel of probability $1 - e^{-\gamma t}$. Dephasing qubits in "+" at rate 1 gives $\langle \sigma_x \rangle = e^{-t}$, which is $0.6065$ at $t = 0.5$, both ways:

myplus = State{Mixed}(System(2, Qubit()), "+")
mydephased = evolve(sum(dephasing_dissipator(1.0)(i) for i in 1:2), 0.5, myplus; nsteps = 10)
measure(mydephased, X)

2-element Vector{Float64}:
 0.6065306597126334
 0.6065306597126329

measure(apply(prod(dephasing_channel(1 - exp(-0.5))(i) for i in 1:2), myplus), X)

2-element Vector{Float64}:
 0.6065306597126334
 0.6065306597126334

Sites reset together are not sites reset each on its own: depolarizing_channel(p, 2) depolarizes both sites with probability p, where depolarizing_channel(p) ⊗ depolarizing_channel(p) depolarizes each with probability p. And dephasing_dissipator(γ, A) damps every coherence between two eigenspaces of A at the rate γ, where Dissipator(A) damps it at $(\lambda_a - \lambda_b)^2/2$, faster for distant eigenvalues; they agree on a qubit, Dissipator(Z) being dephasing_dissipator(2, Z).

The reset, and so the depolarization, changes the charges and is refused on sites conserving something strongly, where a dephasing in the eigenbasis of an operator commuting with the charge, as N, is accepted. Any other channel is a sum of gates, Gate(K) being $\rho \mapsto K \rho K^\dagger$, or the evolution under a Lindbladian L of a few sites for a time t, exp(t * L).

On a qubit, "Up" is $|0\rangle$ and "Dn" $|1\rangle$, the excitation N counts: the decay of an excitation, from "Dn" to "Up", is Dissipator(Sp), and Dissipator(Sm), which takes "Up" to "Dn" in the evolution above, excites. The decay is given for sites of any type by amplitude_damping_channel and amplitude_damping_dissipator, of the excitations an operator destroys, Sp by default, and the relaxation of times $T_1$ and $T_2$ of a qubit, its decay and its dephasing, by thermal_relaxation_channel and thermal_relaxation_dissipator of the module Qubits.

Steady states

steady_state looks for the state $\rho$ with $L\rho = 0$, $L$ being the evolver, starting from a mixed state, by running dmrg on $L^\dagger L$. It returns two things: the residual $\|L\rho\|^2$, for $\rho$ of unit Hilbert-Schmidt norm, which is zero at a true steady state and tells how well the search converged, and the steady state itself, normalized to trace one. A Lindbladian may have several steady states, see Checking the accuracy. Qubits jumping from up to down at rate 1 and from down to up at rate 0.5 settle at $\langle \sigma_z \rangle = -1/3$:

value, mysteady = steady_state(sum(Dissipator(Sm)(i) + 0.5Dissipator(Sp)(i) for i in 1:4),
                               myrho; algo = Dmrg(nsweeps = 10))
measure(mysteady, Z)

4-element Vector{Float64}:
 -0.33333333333333437
 -0.3333333333333349
 -0.3333333333333338
 -0.3333333333333334

Thermal states

thermalize takes a mixed state $\rho$ to $e^{-\beta H/2} \rho \, e^{-\beta H/2}$, normalized, by an evolution in imaginary time, by Tdvp unless its algo says otherwise. From "FullyMixed", the state at infinite temperature, it gives the thermal state $e^{-\beta H}/Z$. It returns the logarithm of the trace the state would have without being normalized, here $\log Z - 4 \log 2$, and the state. Qubits in the field $H = -\sum_i \sigma_z^i$ have $\langle \sigma_z \rangle = \tanh \beta$, which is $0.4621$ at $\beta = 0.5$, and $\log Z - 4 \log 2 = 4 \log \cosh \beta \approx 0.4805$:

mylog, mythermal = thermalize(-sum(Z(i) for i in 1:4), 0.5,
                              State{Mixed}(System(4, Qubit()), "FullyMixed"); nsteps = 10)
mylog

0.48045802783310904

measure(mythermal, Z)

4-element Vector{Float64}:
 0.46211715726000985
 0.46211715726000985
 0.4621171572600098
 0.46211715726000835

For more details, see the Algorithms page of the reference or the inline help.

Measurements

Once we have created a state, we may want to measure it. Take for example a three qubit state

mystate = State{Pure}(System(3, Qubit()), ["+", "Up", "+"])

result = measure(mystate, X(1)X(3))

1.0

will give $\langle \psi | \sigma_x^1 \sigma_x^3 | \psi \rangle$

result = measure(mystate, X)

3-element Vector{Float64}:
 1.0
 0.0
 1.0

will give the array of the $\langle \psi | \sigma_x^i | \psi \rangle$

result = measure(mystate, (X, Y))

3×3 Matrix{ComplexF64}:
 0.0+0.0im          0.0+0.0im  0.0+0.0im
 0.0+0.0im          0.0+1.0im  0.0+0.0im
 0.0-4.26642e-17im  0.0+0.0im  0.0+0.0im

will give the matrix of the $\langle \psi | \sigma_x^i \sigma_y^j | \psi \rangle$, complex since its diagonal is $\langle \sigma_x \sigma_y \rangle = i \langle \sigma_z \rangle$

A fermionic operator is measured as any other, its Jordan-Wigner strings laid for it, in a correlation as (dag(C), C) too. A single one, as C(3), has a value only on a state that superposes parities, and vanishes on any other.

We can also measure properties of the state as a whole, with measurements such as Trace, Purity, EntanglementEntropy(cut) or Fidelity(ref): the Measurements page has the table of them all.

A few sites can be looked at together: ReducedDensityMatrix(positions) measures their density matrix, VonNeumannEntropy(positions) their entropy and LogNegativity(a, b) the entanglement between two parts of them, which on a mixed state the entropy no longer tells. Both are $\log 2$ for one qubit of a Bell pair:

mybell = apply(CX(1, 2) * H(1), State{Pure}(System(3, Qubit()), "Up"))
measure(mybell, [VonNeumannEntropy([1]), LogNegativity([1], [2])])

2-element Vector{Pair{String, Any}}:
 "VonNeumannEntropy(1)" => 0.6931471805599453
   "LogNegativity(1;2)" => 0.6931471805599453

We can also ask for several measurements at the same time

results = measure(mystate, [X, X(2)Z(3), (X, Y), Trace, MemoryUsage])

5-element Vector{Pair{String, Any}}:
           "X" => [1.0, 0.0, 1.0]
   "X(2)*Z(3)" => 0.0
      "(X, Y)" => ComplexF64[0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 1.0im 0.0 + 0.0im; 0.0 - 4.266421588589643e-17im 0.0 + 0.0im 0.0 + 0.0im]
       "Trace" => 0.9999999999999997 + 0.0im
 "MemoryUsage" => 5488

A measurement can also be drawn, as in an experiment: probabilities(state, op) gives the results of measuring an operator placed on sites with their probabilities, sample(state, op) draws one, and collapse(state, op) draws one and gives it with its probability and the state projected onto it, (result, probability, state). Any Hermitian operator of a few sites can be measured, and on any number of sites a string of Pauli operators or the number of particles of a region. The parity of two qubits in "+" is ±1 with probability 1/2, and either result leaves an entangled pair:

mytwo = State{Pure}(System(2, Qubit()), "+")
probabilities(mytwo, Z(1) * Z(2))

2-element Vector{Pair{Float64, Float64}}:
 -1.0 => 0.5
  1.0 => 0.5

myparity, myprob, mypair = collapse(mytwo, Z(1) * Z(2))
myparity, measure(mypair, X(1)X(2))

(-1.0, 1.0000000000000002)

Given the result, collapse(state, op => λ) returns the same, the state projected onto that result, which post-selects it:

_, myprob, mypair = collapse(mytwo, Z(1) * Z(2) => -1)
myprob, measure(mypair, [Z(1)Z(2), X(1)X(2)])

(0.5, Pair{String, Any}["Z(1)*Z(2)" => -1.0000000000000002, "X(1)*X(2)" => 1.0000000000000002])

For more details see the reference or the inline help.

Accuracy

A result is trusted once it no longer moves when the truncation, the time step and the sweeps are improved: how to choose them and how to check what they leave is the subject of Choosing and checking the parameters.