Operators & Dual Unitary Chaos

In this tutorial we use the quantum toolkit, and specifically the operator framework to efficiently validate dual unitary chaos. For further reading you can read the general analytical setup and proof here. We are going to use the numerical setup here. This problem is interesting for a few reasons. First it showcases the power of the operator framework. We will use features that can build quantum circuits as well as construct a custom kernel for a global, diagonal operator that is not directly supported in the quantum toolkit. Second it is physically interesting because it features a space-time swap. We will be computing things on the time lattice and computing thermodynamic limit properties of our system in real space.

The articles aren't necessary to understand this tutorial.

We wish to understand the spectrum of the following operator expression,

$$\mathcal{T} = (\mathbb U_o \otimes \mathbb U_o^* ) \mathcal{O}^\dagger (\mathbb U_e \otimes \mathbb U_e^* ) \mathcal{O}$$

where our system sits on a time lattice of size $4t$ coupling a forward and backward lattice of size $2t$ each. In the expression,

$$\mathbb U_o \otimes \mathbb U_o^*$$

we have a tensor product between the forward and backward lattice where $U_o$ couples odd to even nearest neighbour sites. We take both the forward and backward time lattices to be periodic: if we label our forward lattice with coordinates $j\in \{0, \ldots, 2t-1 \}$ then $j=2t \to 0$.

$$\mathbb U_o = \prod_{m=0}^{t} U_{2m-1,2m}$$

where our gate $U$ is the space-time swapped version of the real-space gate,

$$V_{m,n} = \exp \left[ i\frac{\pi}{4} \left( X_m X_n + Y_m Y_n + J Z_m Z_n \right) \right]$$

you can see a quick explanation for space time swapping a gate here Similar to the odd to even layer of quantum gates, $\mathbb U_e$ couples even to odd nearest neighbour time lattice coordinates.

The final detail of our expression to understand is the operator $\mathcal{O}$. This operator is split into two acting on odd and even sites respectively,

$$\mathcal{O} = \mathcal{O}_o \mathcal{O}_e,$$

and we further break down $\mathcal{O}_s$ into $X$ and $Z$ components. Taking the odd operator as an example,

$$\mathcal{O}_o = \mathcal{O}_o^{(x)}\mathcal{O}_o^{(z)}\mathrm{.}$$

Finally, $\mathcal{O}_o^{(x)}$ acts as a non-expanding operator in the $x$-basis that either leaves a product state untouched (scales by 1) or completely suppresses a product state (scales the coefficient by 0). For this operator, you need to check the magnetization in the x direction for the odd sites in the forward lattice, and likewise for the odd sites in the backward lattice. If the magnetization is identical between the two lattices on odd sites then you scale by unity, and if the magnetization is different you set the coefficient to zero. The other operators are the same, for example $\mathcal{O}_e^{(z)}$ you test the $z$ basis magnetization instead and calculate it only for even lattice sites.

This completes our setup for the problem, let's start coding in aleph!

Investigating the spectrum

Space time swapping

The first step to solve this problem is to construct our real space gates $V_{m,n}$ and space time swap them. This is accomplished by reorganizing the entries of the matrix $V$ that represents our operator.

def st_swapped_du(Js)
{
    // this line spawns the matrix V we want to swap.
    var U = exp(1i*pi/4*matrix(X(0)*X(1) + Y(0)*Y(1)+Js*Z(0)*Z(1)));

    var U_swapped  = zeros([4,4]);
    //here we handle the swapping by collecting our "in" qubits and "out" qubits
    // with a Qbit with 4 sites, in qubits are treated as in the first two
    // out qubits are the second two degrees of freedom. 
    // U_swapped is then filled from the spatial evolution perspective
    for(psi : qbit_range(4))
    {
        var psi_pair = psi.split(2,2);
        U_swapped[psi[2]+2*psi[0],psi[3]+2*psi[1]] = U[psi_pair[0],psi_pair[1]]
    }
    return U_swapped
}

Constructing $\mathbb U_o$ and $\mathbb U_e$

$\mathbb U_o$ and $\mathbb U_e$ can now be constructed from the resulting space-time swapped gate. To do this we can initialize operator products and fill them.

var tmax = 4
// Ug is the matrix that represents the spaced-time swapped V
// assume J has been declared elsewhere
var Ug = st_swapped_du(J)
// Ug sits on the forward lattice, Ugc sits on the backwards lattice.
var Ugc = Ug.conjugated()

// We fill 4 operator products. 
// Uodd for example sits on the forward lattice and couples odd to even sites
// Uoddc is the same but sits on the backward lattice and is constructed with Ugc

var Uodd = operator_prod();
var Ueven = operator_prod();
var Uoddc = operator_prod();
var Uevenc = operator_prod();

// we have to be careful and ensure we obey the PBC
var s1; var s2;
for(var m = 0; m<2*tmax-1;m+=2)
{
    s1 = m;
    s2 = (m+1)%(2*tmax)
    Ueven *= operator_dense(Ug,[s1,s2])

    // the forward lattice is treated as the sites 0 to 2t-1
    // the backward lattice is treated as sites 2t to 4t-1
    s1+=2*tmax;
    s2+=2*tmax;
    Uevenc *= operator_dense(Ugc,[s1,s2])
}
for(var m = 1; m<2*tmax;m+=2)
{
    s1 = m;
    s2 = (m+1)%(2*tmax)
    Uodd *= operator_dense(Ug,[s1,s2])

    s1+=2*tmax;
    s2+=2*tmax;
    Uoddc*= operator_dense(Ugc,[s1,s2])
}

Setting up $\mathcal{O}_s$

Next we want to setup the operators which represent the action of $\mathcal{O}_s$ ($s=o,e$) onto states. Remember that the operators we need to build are diagonal in either the $x$ or $z$ basis, and so we can always rotate the $z$ basis into the $x$ basis and vice-versa. Let's focus on $\mathcal{O}_s^{(z)}$, and then build the required transformation in the $x$ basis via a rotation.

We can build this operator in two ways.

operator_diagonal

The first way is to use operator_diagonal. Here we need to build a 1 dimensional vector that acts as the diagonal of a matrix. To do this is simple, we cycle through every possible product state and test the subsystem magnetizations we care about.

def Oz(listkeep, listdiscard, integer L)
{
    var dim = 2**L
    var diag = zeros([dim]);
    var q = Qbit(L);
    var t = L/4;
    for(var m = 0; m<dim; ++m)
    {
        q.set_representation(m);
        // listkeep picks out either even or odd sites 
        var qt = q.split(listkeep,listdiscard)
        // since we initialized to all 0s to start we just need to test
        // which diagonal entries should instead be a 1. 
        var qs = qt[0].split(t, t)
        if (qs[0].count() == qs[1].count())
        {
            diag[m] = 1;
        }
    }

    return diag;
}

Now that we have the diagonal of the operator inside a vector we can construct the operator_diagonal easily,

// assume L = 4t was declared above 
var evenlist = [0..L][::2]
var oddlist = [0..L][1::2]

var diageven = Oz(evenlist,oddlist,L);
var diagodd = Oz(oddlist,evenlist,L);

var Oze = operator_diagonal(diageven,[0..L]);
var Ozo = operator_diagonal(diagodd,[0..L]);

operator_generic

operator_diagonal might force you to pay a storage cost that you don't want, it is after all storing a vector of size $2^{4t}$. Instead we could just write a custom kernel using operator_generic. This looks of quite similar to constructing operator_diagonal. operator_generic does not store anything but the function you build. The drawbacks of this are clear, you are exchanging the memory overhead issues of operator_diagonal for more computation. At the extremes of system sizes it is often preferable to choose more computation to save memory, but operator_diagonal will always be faster for diagonal operators. To construct an operator_generic we need to write a kernel that accepts an array of coefficients and integers representing product states. Below we show the kernel for even site magnetizations or $\mathcal{O}_e^{(z)}$.

var evenlist = [0..L][::2]
var oddlist = [0..L][1::2]
var f0 = fun[L,evenlist,oddlist](x, y)
{
    // x is the coefficient corresponding to state y
    var q = Qbit(L)
    // y is represented with an integer
    q.set_representation(y[0])
    var qt = q.split(evenlist,oddlist)
    var qs = qt[0].split(L/4, L/4)
    // if the magnetizations don't match we set x to 0
    if (qs[0].count() != qs[1].count())
    {
        x[0] = 0.0
    }
}

This then allows us to construct operator_generic objects to represent our operator.

// the second argument is a list of lattice sites
// here the operator is global so we leave this empty
var Oze = operator_generic(f0, [], as_complex);
var Ozo = operator_generic(f1, [], as_complex);

Rotating basis

We can easily switch basis by exploiting the RotY operator and rotating by $\frac{\pi}{2}$ from the $z$ basis to the $x$ basis. Constructing the operator_prod that does this is simple. Don't forget to build the expression that rotates back to the $z$ basis.

var Uy = operator_prod()
var Uydag = Uy;
for(var m = 0; m<L; ++m)
{
    Uy*= RotY(pi/2,m);
    Uydag*= RotY(-pi/2,m);
}

Putting it all together

Since we are at the dual unitary point of this problem, the maximal eigenvalue should be $1$, in fact, there should be exactly $t$ unit eigenvalues. In this tutorial, we will write a simple power method to project out the leading eigenvector, but first we need to assemble the full operator and build a vector to act on.

// random guess vector
var psi = random([2**L]);
// normalize it
psi.normalize();
// full operator expression for T
var Utotal = (Uodd*Uoddc) * (Oze*Ozo) * (Uydag * Oze * Ozo * Uy) * (Ueven*Uevenc) *  (Uydag*Ozo*Oze*Uy) * (Ozo*Oze);

// we use >> to perform all operations in place, avoiding copies, just the vector is modified. 
for(var m = 0; m<8; ++m)
{
    Utotal >> psi;
    print(psi.norm())
    psi.normalize();
}
// when the state's norm isn't changing away from unity anymore we have reached the eigenvector we want

What's next?

• Put everything together into a workable file. You'll need to define tmax, L = 4*tmax and pick a J for your gate
• Compare the speed of operator_diagonal and operator_generic? Can you make the operator_generic implementation faster? What about operator_diagonal?
• Here we did a simple power method, to improve this which sparse eigensolver should you use? Is T / Utotal hermitian?
• What is the memory overhead for this numerical experiment?

Documentation Contributors

Jonathon Riddell