Floquet and quantum circuit models

Floquet models, or more generally quantum circuit models, are two common systems in the study of quantum many body physics. They offer a minimal environment for the study of dynamics and chaos, while also being closely related quantum information and computing contexts. In this tutorial we will discuss both classes of models and how to work with them numerically.

Kicked spin $\frac{1}{2}$ Ising model

A prototypical example of a Floquet model is the kicked Ising spin model. The time dependent Hamiltonian is written as,

$$\hat{H}_{KI}(t) = \hat{H}_I + \delta_p(t) \hat{H}_X,$$

where $\hat{H}_I = \sum_{m=0}^{L-1} J Z_m Z_{m+1} + h_m Z_m$ and $\hat{H}_X = b \sum_{m=0}^{L-1} X_m$. The $J$ and $b$ constants are model dependent while the $h_m$ are position dependent field strengths typically drawn randomly from a distribution. The time dependent $\delta_p(t) = \sum_{m =-\infty}^{\infty}\delta(t-m)$ is the periodic delta function. The Floquet operator generated by these dynamics reads,

$$U_F = T \exp \left[ T \int_0^1 \hat{H}_{KI}(t) dt \right] = e^{-i \hat{H}_X} e^{-i \hat{H}_I}.$$

Quantum states can therefore be evolved in discrete steps as,

$$|\psi(t)\rangle = U_F^t |\psi\rangle,$$

where $t$ is taken to be an integer. $U_F$ can be broken down into a circuit like structure. To see this first note that every operator in both exponentials commutes with one another so we may write,

$$e^{-i \hat{H}_X} = \prod_{m=0}^{L-1} e^{-i b X_m},$$

and likewise we split $\hat{H}_I$

$$e^{-i \hat{H}_I} = \prod_{m=0}^{L-1} e^{-i \left(J Z_m Z_{m+1} + h_m Z_m\right)}.$$

This model has therefore been effectively decomposed into local unitary (operators)[https://docs.kotharcomputing.com/docs/learning-center/Concepts/Many-Body/Operators]. We can represent $U_F$ in aleph with an operator_prod

var L = 28
var b = 1.0
var J = 1.0
var h = random([L])
var circ_Hx = operator_prod()
var circ_HI = operator_prod()
for(m : [0..L])
{
    // we use rotation in X as a shorthand for our term
    // taking care to cancel the 1/2 factor
    circ_Hx *= RotX(2*b,m)
    // we likewise extract the diagonal of our Ising Z terms
    // and store it in an operator diagonal, which make the code faster
    var op_HI = 1i*(J*Z(m)*Z((m+1)%L)+h[m]*Z(m))
    var diag = exp(matrix(op_HI)).diagonal()
    circ_HI *= operator_diagonal(diag,[m,m+1])
    
}
// combining everything is easy
var circ_KI = circ_Hx * circ_HI

tip

Most local operators benefit from specialized in place kernels when acting on state vectors. To use in place kernels, update vectors with >> instead of *. However be cautious, this modifies the vector in place, meaning the vector itself is always modified.

This is of course just an example, Floquet physics and likewise more broadly quantum circuit models can take many forms and typically aren't this simple.

Local unitary operators

Operators that typically appear in the dynamics of Floquet or quantum circuit models are local, meaning they act on one, two or at most a few degrees of freedom at once. An example of this might be a layer of two-local, nearest neighbor gates coupling even to odd degrees of freedom. This would be one layer of a brickwork circuit model. Lets build this with Haar random gates for simplicity. We will take PBC.

var L = 30
var circ_haar = operator_prod()
// only select even m
for(m: [0..L][::2])
{
    // we generate a new Haar random unitary for each coupling
    var mat = haar_matrix(4)
    circ_haar *= operator_dense(mat,[m,(m+1)%L])
}

tip

Aleph supports a comprehensive list of possible operators, which need not be unitary. To see a comprehensive list checkout this concept

Global unitary operators

Sometimes Floquet evolutions cannot be broken down into local quantum circuit models without doing an approximation like a Trotter decomposition or Strang splitting. Unless you are putting this circuit on a quantum computer, these approximations are unnecessary and build up error quickly. Instead one can use operator_function similar to this tutorial. Suppose I have an evolution that looks like,

$$U = exp(- i J \sum_{m=0}^{L-1} X_m Z_{m+1} X_{m+2})$$

Note, each term in the sum does not commute with it's neighbor and we cannot straightforwardly break this unitary up into local quantum gates acting on three qubits. Unitary matrices like this can be represented easily by combining operator_sum with an operator_function. We just need to choose the size of our Krylov basis, and instead of a time increment we specify $J$.

var L = 30
var J = 1.0;
var opexpress = operator_sum()
for(m: [0..L])
{
    opexpress += X(m)*Z((m+1)%L)*X((m+2)%L)
}

var Nkrylov = 10;
var expiop = fun[J](x) { return exp(-1i * J * x) }
var U = operator_function(opexpress, expiop, ["krylov_dimension" : Nkrylov]);

tip

operator_function, like all operator expressions, can be combined into products. So you can take advantage of the exponential accuracy of operator_function in conjunction with decomposed quantum circuit expressions in the same operator_prod.

Worked example

We can put everything together that we have looked at so far in this tutorial and create a strange, but instructive example. Let us begin our quantum system in the all zero state $|\psi\rangle = |0\rangle$ and we will do a Floquet evolution,

$$U_F = U_1 U_0$$

and we take $|\psi(t)\rangle = U_F^t |0\rangle$ with $t$ an integer. We will use the unitary operators from the previous sections to define $U_1$ and $U_0$,

$$U_0 = \prod_{m=0}^{L/2-1} V_{2m,2m+1}$$

where the matrix $V$ is a Haar random unitary. We next take,

$$U_1 = exp(- i J \sum_{m=0}^{L-1} X_m Z_{m+1} X_{m+2})$$

The following code time evolves $|\psi\rangle$ with $U_F$,

// system size
var L = 18
// coupling coefficient
var J = 1.0;
// max time
var tmax = 10;

// first lets build the operator function
var opexpress = operator_sum()
for(m: [0..L])
{
    opexpress += X(m)*Z((m+1)%L)*X((m+2)%L)
}
// operator function construction
var Nkrylov = 6;
var expiop = fun[J](x) { return exp(-1i * J * x) }
var U1 = operator_function(opexpress, expiop, ["krylov_dimension" : Nkrylov]);

// next we do the Haar layer
var U0 = operator_prod()
// only select even m
for(m: [0..L][::2])
{
    // we generate a new Haar random unitary for each coupling
    var mat = haar_matrix(4)
    U0 *= operator_dense(mat,[m,(m+1)%L])
}

// we put both together now
var U_F = U1*U0

// |0> initial state
var psi = zeros([2**L],as_matrix)
psi[0] = 1
// acting in place will speed up the U0 layer. 
for(t : [0..tmax])
{
    U_F >> psi
}

What's next

• Take our example and track the entanglement entropy of our state $|\psi\rangle$. This can be accomplished fetching the reduced density matrix of a subsystem with this function and then calculating the entropy with this function.
• Design a new Floquet or quantum circuit model. Can you put onsite random unitary gates?
• How large of a system can you simulate with this approach?
• Implement a trotter decomposition using what you've learned in this tutorial. How much more accurate is our evolution compared to a Trotter decomposed version? Which is ultimately faster to simulate?