Operations

In this section we will see how to use the functions defined in this package to create operators acting on mixed states.

Left- and right-multiplication

For each vectorised site type, each operator from the “original”, non-vectorised site type can be used to multiply a state (in the sense of this package, that is a mixed state or an operator) on the left on the right. In other words, if there exists an operator "A" for the "T" type, then the operators "A⋅" and "⋅A" (the dot is a \cdot) are automatically defined for the "vT" site type, which perform left- and right-multiplication by "A", respectively.

For example, let's play around with a vFermion site. We create a site index and define on it the empty state \(\proj{0}\), then apply the creation operator on the left.

julia> s = siteind("vFermion");

julia> emp = state(s, "Emp");

julia> Adagemp = apply(op("Adag⋅", s), emp);

We have now obtained \(\adj{a}\proj{0}\), and we can check for example that its inner product with \(\proj1\) is zero, that is, \(\tr(\proj1 \adj{a}\proj{0}) = \braket{0}{1} \bra1 \adj{a}\ket0 = 0\):

julia> occ = state(s, "Occ");

julia> scalar(occ * Adagemp)
0.0 + 0.0im

If we multiply \(\adj{a}\proj{0}\) by \(a\) on the right too, then we get \(\tr(\proj1 \adj{a}\proj{0}a) = \bra0 a \ket1 \bra1 \adj{a}\ket0 = 1\):

julia> AdagempA = apply(op("⋅A", s), Adagemp);

julia> scalar(occ * AdagempA)
0.9999999999999998 + 0.0im

Left- and right-multiplication operators can be composed with the usual ITensor rules: for example, they can be multiplied or added together. In fact, we could have defined AdagempA in one shot by applying the "Adag⋅ * ⋅A" operator:

julia> AdagempA = apply(op("Adag⋅ * ⋅A", s), emp);

Commutators

The time evolution of mixed states is often described by the von Neumann equation, or the GKSL equation, both of which involve the commutator of the state and an operator. For this reason, this package provides a useful feature with which we can create the \(\rho \mapsto -\iu [A,\rho]\) “superoperator”: this is the gkslcommutator function.

Given a list of operator names and integers, in an alternating fashion like "A1, n1, A2, n2, ...", the gkslcommutator function returns an OpSum object that represents the \(-\iu[A,\blank]\) map where \(A\) is the product of A1 on site n1, A2 on site n2 and so on, in a similar syntax as OpSum itself. For example,

julia> gkslcommutator("A", 1)
sum(
  0.0 - 1.0im A⋅(1,)
  0.0 + 1.0im ⋅A(1,)
)

julia> gkslcommutator("A", 1, "B", 3)
sum(
  0.0 - 1.0im A⋅(1,) B⋅(3,)
  0.0 + 1.0im ⋅A(1,) ⋅B(3,)
)

These OpSum objects can then be turned into MPOs with the usual syntax MPO(opsum, sites).

If we also provide an array of site indices, we can use the gkslcommutator_itensor function, that directly gives us an ITensor object.

julia> sites = siteinds("vS=1/2", 3);

julia> gkslcommutator_itensor(sites, "Sx", 2)
ITensor ord=2 (dim=4|id=652|"Site,n=2,vS=1/2")' (dim=4|id=652|"Site,n=2,vS=1/2")
NDTensors.Dense{ComplexF64, Vector{ComplexF64}}

Example: GKSL equation

Consider the Hamiltonian operator

\[H = \sum_{n=1}^{N} \omega_n\phantomadj \adj{a_n} a_n\phantomadj\]

for a harmonic oscillator. Let's choose 8 as the dimension of the Hilbert space, \(N=3\) and all \(\omega\sb{n}\)s equal to one, for simplicity. The \(-\iu[H,\blank]\) map in the von Neumann equation \(\dot\rho\sb{t} = -\iu[H,\rho\sb{t}]\) can be defined as follows.

julia> N = 3; s = siteinds("vBoson", N; dim=8); ω = ones(N);

julia> h = OpSum();

julia> for n in 1:N
           h += ω[n] * gkslcommutator("N", n)
       end

julia> h
sum(
  0.0 - 1.0im N⋅(1,)
  0.0 + 1.0im ⋅N(1,)
  0.0 - 1.0im N⋅(2,)
  0.0 + 1.0im ⋅N(2,)
  0.0 - 1.0im N⋅(3,)
  0.0 + 1.0im ⋅N(3,)
)

We could also add to the Hamiltonian an exchange interaction

\[\sum_{n=1}^{N} \lambda_n\phantomadj (\adj{a_n} a_{n+1}\phantomadj + \adj{a_{n+1}} a_n\phantomadj)\]

with the following commands (we set \(\lambda\sb{n} = \frac12\)).

julia> λ = fill(0.5, N);

julia> for n in 1:N-1
           h += λ[n] * gkslcommutator("Adag", n, "A", n+1)
           h += λ[n] * gkslcommutator("A", n, "Adag", n+1)
       end

julia> h
sum(
  0.0 - 1.0im N⋅(1,)
  0.0 + 1.0im ⋅N(1,)
  0.0 - 1.0im N⋅(2,)
  0.0 + 1.0im ⋅N(2,)
  0.0 - 1.0im N⋅(3,)
  0.0 + 1.0im ⋅N(3,)
  0.0 - 0.5im Adag⋅(1,) A⋅(2,)
  0.0 + 0.5im ⋅Adag(1,) ⋅A(2,)
  0.0 - 0.5im A⋅(1,) Adag⋅(2,)
  0.0 + 0.5im ⋅A(1,) ⋅Adag(2,)
  0.0 - 0.5im Adag⋅(2,) A⋅(3,)
  0.0 + 0.5im ⋅Adag(2,) ⋅A(3,)
  0.0 - 0.5im A⋅(2,) Adag⋅(3,)
  0.0 + 0.5im ⋅A(2,) ⋅Adag(3,)
)

Adjoint map

Another common operation on mixed states is \(\rho \mapsto X \rho \adj{X}\), where \(X\) is some operator. For example, in the definition of the dissipation terms in the GKSL equation (see the relative example) we have \(a\sb{n}\phantomadj \rho \adj{a\sb{n}}\), which can be simply written as apply(op("A⋅ * ⋅Adag", s, n), ρ), where s is the list of Index object that ρ is defined on. This is straightforward because the A and Adag operators are already defined, but not all operators have their adjoint available in the library.

The adjointmap_itensor^[1] function provides a convenient way of creating the tensor representing the above map. Define, for example, the following tensors t1 and t2:

julia> s = siteinds("Boson", 4; dim=4);

julia> s_vec = siteinds("vBoson", 2; dim=4);

julia> ρ_vec = state("0", s_vec[1]) * state("2", s_vec[2]);

julia> t1 = random_itensor(s[1], s[1]');

julia> t2 = random_itensor(s[1], s[2], s[2]', s[1]');

The adjointmap_itensor takes the tensor as the first argument; then, we need to provide the original sites on which the tensor acts (the orig_sites keyword argument, mandatory) and the sites of the vectorised site type which the returned tensor will act on (the vec_sites keyword argument, mandatory). This is necessary since the original tensor and the result of adjointmap_itensor are defined on different sets of indices which in principle might have nothing to do with each other. Note that orig_sites and vec_sites are always vectors of indices, even when there is only one index (per argument).

The function automatically computes the adjoint of the given operator and creates the tensor corresponding to the vectorisation of the \(\rho \mapsto X \rho \adj{X}\) map, that can then be applied to MPSs or other ITensors of the vectorised site type.

julia> t1_vec = adjointmap_itensor(t1; orig_sites=[s[1]], vec_sites=[s_vec[1]]);

julia> apply(t1_vec, ρ_vec);

julia> t2_vec = adjointmap_itensor(t2; orig_sites=s[1:2], vec_sites=s_vec[1:2]);

julia> apply(t2_vec, ρ_vec);

If the operator is already defined with an OpName for the non-vectorised site type, then there is a simpler syntax, similar to the op function: for example, we can write

julia> s = siteinds("vQubit", 6);

julia> ρ = MPS(s, "0");

julia> apply(adjointmap_itensor("CNOT", s[1], s[2]), ρ);

julia> apply(adjointmap_itensor("Rx", s, 5; θ=pi/3), ρ);

In these cases we skip the creation of the non-vectorised tensor altogether. We also did not need to define the adjoint of "CNOT" or "Rx": the adjointmap_itensor already takes care of computing it internally, without having to define a new operator.