Chapter 2 dynamite: massively parallel numerics for many-body quantum spin systems

Related published works:

• Exponentially slow heating [arxiv] [journal]

• Long-range prethermal phases [arxiv] [journal]

Closed quantum systems subjected to a drive tend to heat up to infinite temperature, because the drive inputs energy to the system but no energy is ever dissipated. At first sight, this means that most driven closed quantum systems are not particularly interesting—the infinite temperature state corresponds with a density matrix that is the identity, and is thus featureless. One way of avoiding this infinite-temperature fate is by setting up the system to exhibit many-body localization (MBL), which prevents the system from thermalizing, even at infinite time (MBL is explored further in Section 2.3.2). However, another way to observe interesting physics is to look at the system’s behavior pre-thermal regime: the period of time before it has had the chance to fully heat up. Amazingly, analytic studies have shown that in certain Floquet systems—those in which the drive is periodic—this pre-thermal regime can exist for a time scale at least exponentially long in the frequency of the drive. [MKS16, KMS16, ADH+17a, EBN17, ADH+17b] The hope is that such a long prethermal regime provides a sufficient window of time to observe interesting phenomena. Here, we numerically study the time evolution of a set of related Hamiltonians, and observe that long-range order induced by an effective Hamiltonian in the prethermal regime can lead to a novel out-of-equilibrium phase of matter called a discrete time crystal—in which a discrete time translation symmetry is spontaneously broken.

We begin by considering a 1D spin chain with open boundary conditions, with a long- or short-range $ZZ$ interaction (either decaying as a power law, or nearest-neighbor). In addition we apply a nearest-neighbor $XX$ interaction and a uniform, static magnetic field $\overrightarrow{h}$. This interaction has the form:

where the angle brackets on the second term indicates it is only a nearest-neighbor interaction. We may consider the nearest-neighbor interaction to correspond to the case of $\alpha =\infty$. In dynamite, this Hamiltonian itself can be implemented in just a few lines of code, as shown in Code Listing 2.1.

For the Floquet drive, after every period $T$ of time evolution the system will undergo a global $\pi$-pulse which we denote as $X$, rotating all spins by ${180}^{\circ }$ around the $x$-axis. (Note that this is equivalent to flipping the direction of the magnetic field $\overrightarrow{h}$ across the $x$ axis every time $T$).

There are a few quantities we would like to track during the evolution of this system. First, we’d like to measure the energy of the effective Hamiltonian $D^{\ast }$ under which the system evolves (in a frame that flips along with the $\pi$-pulses $X$). The exact form of $D^{\ast }$ cannot be expressed concisely, but we can approximate it as the Hamiltonian of Eq. 2.1 averaged over a full cycle of length $2T$, which we denote as $D$:

In addition to the expectation value of $D$, we compute the magnetization of each spin in the $Z$ direction as the expectation values of $S_i^z$ for all $i$, as well as the entanglement entropy. With that, we can implement the Floquet evolution of our system, as shown in Code Listing 2.2. Note how $D$ and the global pi-pulse are extremely easy to implement in dynamite via its support for arbitrary operator arithmetic. Furthermore, both operators are very fast to compute because the arithmetic is implemented symbolically via Pauli operations, rather than matrix-matrix products.

Using dynamite to evolve this system yields the results plotted in Figure 2.2. Here, we have used a system size of $L=22$ spins, with the parameters $\alpha \in \{1.13,\infty \}$ (for long- and short-range interactions respectively), $J_z=1$, $J_x=0.75$, and $\overrightarrow{h}=(0.21,0.17,0.13)$.¹¹ 1 For the time evolution in the figure we have built the Hamiltonian with closed rather than open boundary conditions, for which we define the long range-interaction strength by replacing $|i-j{|}^{-\alpha }$ with ${\left[\right(L/\pi )\sin |i-j|\pi /L]}^{-\alpha }$, which matches the periodicity of the boundary conditions. The frequency $\omega$ is the inverse of the time $T$ that is allowed to elapse before the spin flip operator $X$ is applied. In all cases the evolution starts from a product state; the precise product state is labeled in Figure 2.2 and is varied in the long-range case to explore both low-energy (“cold”) and high-energy (“hot”) initial states.

Begin by observing the top row of plots of Figure 2.2, which show the expectation value of the approximate effective Hamiltonian $D$ throughout the time evolution. In every case (short-range interactions, and both temperatures of long-range interactions) we observe a clear signature of a pre-thermal regime corresponding to the effective Hamiltonian approximated by $D$. The expectation value of $D$ is conserved for the period in which the system’s evolution is well-approximated by it, and only after a sufficiently long time does the system begin to heat up and eventually reach a thermal state in which $⟨D⟩=0$. As expected, we observe that the length of the prethermal regime depends exponentially on the frequency $\omega$. The same pattern is shown in the entanglement entropy, in the middle row of plots, with the addition of the fact that we see the initial effective thermalization of the product state to the thermal state with respect to $D$! In particular, the initial product state has no entanglement, and after a time ${\tau }_{pre}$ (which does not depend on $\omega$) the system reaches an entanglement entropy corresponding to a thermal state of the same temperature as the initial state (with respect to $D$). This is made particularly clear by comparing to the black dotted line, which plots the limit of $\omega =\infty$ (implemented by just evolving under $D$ directly, which is exactly equal to $D^{\ast }$ in this limit). The “kink” corresponding to where the system has thermalized with respect to $D$ but has not yet left the prethermal regime is called the “prethermal plateau.” Finally, in the bottom row of plots, we observe the time-crystalline behavior that was alluded to at the beginning of this section, in the form of period-doubling: the magnetization only returns to its original value every other period of the drive, corresponding to a spontaneously broken time translation symmetry. Note in particular that in the short-range and hot long-range cases, this time crystalline phase is not stable: it disappears by the time ${\tau }_{pre}$ is reached, because in these cases the effective Hamiltonian $D^{\ast }$ is not capable of supporting long-range magnetic order (the time crystal “melts”). On the other hand, in the cold long-range case, $D^{\ast }$ supports stable long-range ferromagnetic order, and the time crystalline behavior survives all the way until the system leaves the prethermal regime and begins to thermalize to its final infinite temperature state—a length of time that is exponentially dependent on the frequency. For more details and corresponding analytic analysis, see the published manuscripts corresponding to this study. [MKE+19, MEK+20]

One of the main motivations for pushing numerics to the largest possible system sizes is to study emergent phenomena which simply cannot exist in just a few spins. A prototypical example is many-body localization (MBL): a surprising phenomenon in which certain many-body systems with sufficiently strong disorder fail to thermalize. In particular we will study the MBL transition, in which the system moves from thermalizing to localized as the disorder strength is increased. Characterization of this transition has remained elusive: finite size effects seem hard to avoid, and tensor network methods break down due to extensive entanglement in states near the transition. Thus, iterative methods like the Krylov subspace methods used by dynamite have proved to be one of the best tools for its study. [PML+18, LLA15] We refer readers interested in learning more about the physics of MBL to one of the excellent review papers on the subject. [AAB+19]

Here, we show how to implement a model of nearest-neighbor Heisenberg interactions on a 1D chain with open boundary conditions, with disorder implemented as random $Z$ fields on each site:

where $\overrightarrow{S}=(S^x,S^y,S^z)$, the subscripts indicate the index of the spin in the chain, and the angle brackets indicate that the indices run over nearest neighbors. The values of $h_i$ are drawn from a uniform distribution $[-h,h]$ where $h$ is a parameter that controls the strength of the disorder. The implementation in dynamite is presented in Code Listing 2.3. This Hamiltonian conserves total magnetization $\sum S^z$, which we use to reduce the size of the Hilbert space for our computations.

In this example we identify the MBL transition via the half-chain entanglement entropy of eigenstates, $S_{L/2}$, which is simply the bipartite von Neumann entropy when half the spins are traced out. The MBL transition should correspond to a transition from volume law to area law entanglement. In the full example included with the dynamite source, we also examine an eigenvalue statistic called the “adjacent gap ratio.”

The key feature that makes MBL so interesting is that the transition from volume to area law entanglement does not only occur in the ground state, but in excited states as well. This presents a great use case for dynamite’s “target” eigenvalue solver, which finds the $k$ eigenvalues (and eigenvectors) closest to a target energy, where $k$ is user configurable. Code Listing 2.4 implements the following computation: 1) choose an energy in the middle of the spectrum, 2) solve for 32 eigenpairs near this point, and then 3) compute the entanglement entropy for all of those eigenpairs.

An extensive numerical study of the MBL problem is presented in Chapter 3.

• Related published work: [arxiv] [journal]

In spirit the Sachdev-Ye-Kitaev (SYK) model represents the opposite of the localization in the previous example: it is expected to scramble information highly efficiently via chaotic many-body dynamics. [KIT15, MS16] Indeed, it exhibits maximal chaos: the Lyapunov exponent, which characterizes the rate of chaos, saturates its upper bound of $2\pi T$, where $T$ is the temperature of the system. [MSS16] Its physics can be connected to the dynamics of quantum information in black holes, providing a testbed for exotic phenomena such as scrambling-based teleportation. [GJW17, MSY17, BGL+23, SKG+22]

The SYK model gives us a chance to look at how quantum systems other than spins can be explored with dynamite, by transforming them onto a spin system. The SYK model we’ll use consists of Majoranas interacting in 0D, with random couplings. Specifically it consists of every possible 4-body interaction among N Majoranas, with each term having a random coupling strength:

where $J_{ijkl}$ are randomly chosen from a Gaussian distribution with variance 1.

To map the Majoranas onto the spin systems that are natively supported in dynamite, we can use the following transformation. For the Majorana with index $i$, let $q=\lfloor i/2\rfloor$. Then

where the first Pauli is ${\sigma }^x$ if $i$ is even and ${\sigma }^y$ if it’s odd. In words, the Majorana consists of a ${\sigma }^x$ or ${\sigma }^y$ with a string of ${\sigma }^z$ extending to the edge of the spin chain. Note that we get two Majoranas for each spin! This is straightforward to implement in dynamite, but is actually already built in in the dynamite.extras module so we don’t have to do it ourselves. Using that, an implementation of this Hamiltonian is shown in Code Listing 2.5. Note that this Hamiltonian has enough terms that building it can take an appreciable amount of time; a more efficient but less transparent version of the build_hamiltonian function can be found in the full SYK example, distributed with the source code of dynamite. The full example also takes advantage of a $Z_2$ symmetry in the Hamiltonian to speed up the computations.

We can study the fast scrambling behavior of the SYK model by examining out of time order correlators (OTOCs). In particular, we will measure to what extent two local operators $V\left(0\right)$ and $W\left(t\right)$ anticommute for various times $t$, where the anticommutator at time $t=0$ is zero:

It is helpful to reduce this to the following equivalent expression

which is the formulation of $C\left(t\right)$ that we will use here.

Defining $O\left(t\right)=W\left(t\right)V\left(0\right)W\left(t\right)V\left(0\right)$, we are specifically interested in this operator’s expectation value with respect to thermal states, of various (inverse) temperatures $\beta$. Now, dynamite’s speed comes from the fact that it works with pure states, rather than mixed states—so the obvious plan to just compute $Tr\left[O\right(t\left)e^{-\beta H}\right]$ is out of the question. Instead, we can take advantage of an idea called “quantum typicality” to get an estimate of the expectation value more efficiently. [BG09, SGB14] Quantum typicality says that $Tr\left[O\right(t\left)e^{-\beta H}\right]$ is well approximated by the expectation value of $e^{-\beta H/2}O\left(t\right)e^{-\beta H/2}$ with respect to random states (where we have split the thermal operator in half across the trace to get a more symmetric result). For simplicity we can set $W={\chi }_0$ and $V={\chi }_1$. With that, for a uniformly random state $|{\psi }_r⟩$ we will compute (writing things out in full):

In Code Listing 2.6 we demonstrate how the expectation value of this operator can be computed in dynamite.

In Figure 2.3, we compute $F\left(t\right)$ which is closely related to $C\left(t\right)$: it is the “regularized” OTOC in which the thermal density matrix is spread equally across the operator, see the supplementary information of the published manuscript for details [KYK+21]. Note that there are a lot of independent parameters: the time $t$, temperature $\beta$, and system size $N$ can all be varied, and to achieve good statistical significance requires then averaging over both the disorder in $H$ and the random states used in the quantum typicality computation. Figure 2.3(a) shows the dependence of $F\left(t\right)$ on $t$ for various system sizes, at $\beta J=10$. (The plotted value is $\tilde{F}\left(t\right)=F\left(t\right)/F\left(0\right)$). There are two main observations to be made from this figure. First, the scrambling behavior is clear: $F\left(t\right)$ decays toward zero with a certain time scale, determined at early times by the Lyanpunov exponent $\lambda$ as $1-\tilde{F}\left(t\right)\sim e^{\lambda t}/N$. Second, there is clear flow of $\lambda$ with system size. A finite-size rescaling procedure can be used to extrapolate $\lambda$ to the limit of infinite system size, at which point it can be plotted as a function of the temperature $\beta J$ (Figure 2.3(b)). The observed dependence of $\lambda$ tracks closely with the result of a semiclassical solution (the Schwinger-Dyson equations, labeled “SD” in the figure), and at large $\beta$ (low temperature), the theoretically predicted fast-scrambling limit of $\lambda \le 2\pi /\beta$. For more details, including of the finite-size rescaling procedure and the prediction from Schwinger-Dyson equations, see the published manuscript presenting these results. [KYK+21]

The results summarized in Figure 2.3 and elucidated in [KYK+21] provided the first numerical evidence supporting the claim that the SYK model exhibits fast scrambling behavior at low temperature. Previous efforts at smaller system sizes were unable to observe this behavior due to finite size effects; this study was made possible by dynamite’s speed at very large system sizes which enabled analysis of Hamiltonians on up to 60 Majoranas. Due to the four-body all-to-all connectivity of the model (Eq. 2.4), the Hamiltonian has an extremely large number of terms in the SYK Hamiltonian (487,635 for a system size of 60 Majoranas!), explicitly storing the (sparse!) matrix in memory is infeasible for large problems: for 60 Majoranas it would require over 350 terabytes of RAM! Thus dynamite’s matrix-free (“shell”) mode, which only stores and manipulates the operator symbolically, was crucial for this study.

After these results were published, two studies followed in which I was not directly involved, but warrant mention for the curious reader. The many-body teleportation of quantum states via scrambling, a phenomenon that emphasizes the duality between quantum mechanical models like SYK and the dynamics of gravity models, in particular black holes and wormholes, was numerically observed with dynamite. [SKG+22] In another study, dynamite was used to show that sparse SYK, in which many of the $J_{ijkl}$ are zero, can reproduce many of the features of the regular SYK model, but with a much lower computational cost to the simulation. [CMP21]

The two main computationally intensive tasks that dynamite performs are time evolution and eigensolving. At first glance, for an operator $H$, the computational task of quantum time evolution is to compute the unitary $U=\exp (-iHt)$, and for eigensolving it is to compute a matrix $A$ and diagonal matrix $\Lambda$ such that $AH=\Lambda H$. Both $U$ and $(A,\Lambda )$ can be computed in a straightforward manner using standard linear algebra subroutines (for example, SciPy’s linalg.expm and linalg.eigh routines); however, as we explain below, doing so is exponentially more expensive than is necessary for most physics studies. For an quantum operator on a Hilbert space of dimension $N$, matrices $U$ and $A$ both of dimension $N\times N$, consisting of $4^L$ complex numbers for a generic system of $L$ spins-1/2 particles. The power of the Krylov subspace algorithms underlying dynamite stems from the fact that they avoid ever explicitly computing $U$ or $A$, only requiring $O(2^L\cdot poly(L\left)\right)$ memory to operate and having a corresponding improvement in runtime. Here we give a brief overview of how they work; for details please see the SLEPc Users Manual and Technical Reports. [RCD+23, HRT+07] Note that if the dynamite user wants to tune the solver by passing any of the “command line” options to SLEPc as described in those documents, it can be achieved via the slepc_args argument to dynamite.config.initialize().

An initial concern seems to be that the matrix $H$ itself is of dimension $N\times N$, which ostensibly is an obstacle to achieving the efficiency just promised. However, $H$ has a crucial feature that $U$ and $A$ do not: it is sparse. This means that most of the matrix elements are zero, and the cost of storing and working with the matrix can be reduced dramatically by only storing the non-zero elements. In fact, the memory usage can be reduced even further by not storing any of the matrix elements at all, as we discuss below in Section 2.4.2. With that, we move on to discussing how we can use this sparse $H$ to achieve performant time evolution and eigensolving.

For time evolution, the key is that we are often more interested in the action of the matrix $U$ on a state vector than we are in the matrix itself. For an initial state $|{\psi }_0⟩$, we desire to directly compute the time evolved vector $|{\psi }_t⟩=\exp (-iHt)|{\psi }_0⟩$ without needing to compute the unitary as an intermediate step. The Taylor expansion of the matrix exponential suggests that for reasonably small $||Ht||$, $|{\psi }_t⟩$ will be dominated by a linear combination of the vectors $H^i|{\psi }_0⟩$ for $i\in [0,m]$, where $m$ is a small integer, say 30. These vectors are easy to compute by repeatedly multiplication with $H$. This intuition forms the core idea behind Krylov subspace methods: construct the subspace $K_m=span\left(\right\{|{\psi }_0⟩,H|{\psi }_0⟩,H^2|{\psi }_0⟩,\cdots ,H^{m-1}|{\psi }_0⟩\left\}\right)$, and solve the linear algebra problem on this much smaller subspace. For time evolution, this means projecting $H$ onto $K_n$ to yield a matrix $\tilde{H}$ of dimension only $m\times m$, computing the matrix exponential $\tilde{U}=\exp (-i\tilde{H}t)$ (a trivial computation for a small matrix like $\tilde{H}$), and then computing $|{\psi }_t⟩\approx \tilde{U}|{\psi }_0⟩$. Projected back into the full Hilbert space, for small $||Ht||$, $\tilde{U}$ should be a good approximation of $U$, and thus the result should be a good approximation of $|{\psi }_t⟩$.

The case in which $||Ht||$ is not sufficiently small presents the opportunity to discuss another powerful technique crucial to these algorithms, which is called restarting. Instead of computing a single Krylov subspace and attempting to perform the whole computation in it, the subspace is iteratively updated throughout the computation until convergence is achieved. The default restarting scheme for time evolution in dynamite is quite simple: the desired evolution of time $t$ is broken down into a series of shorter evolutions of time $\Delta t$, which is sufficiently small that $\tilde{U}$ is a good approximation of $U$. After each evolution of $\Delta t$, the Krylov subspace is constructed anew, such that it tracks the evolution of the state vector through the Hilbert space. This is the algorithm used by the Expokit package, which is reproduced in SLEPc; [SID98, RCD+23] there is also a Krylov-Schur solver with a less intuitive restarting scheme available in SLEPc but we find it to be less stable.

For eigensolving, we take advantage of the fact that computing the matrix $A$ of all eigenvectors, which is enormous, may be overkill. Physics studies are often interested in only a few particular eigenstates, such as the ground state and first few excited states. Matrix-vector products are again useful here, with the simplest algorithm being the so-called “power method.” Consider the decomposition of a random vector in the basis of (unknown) eigenstates of $H$: it will have the form $|{\psi }_r⟩={\sum }_i{c}_i|{\varphi }_i⟩$, where each $|{\varphi }_i⟩$ is an eigenstate. When $H$ is multiplied by this vector several times, the result is $H^m|{\psi }_r⟩={\sum }_i{c}_i{\lambda }_i^m|{\varphi }_i⟩$—the eigenvectors with eigenvalues of largest absolute value have been enhanced exponentially; for sufficiently large $m$ the result will converge to the most extremal eigenvector. Once again, this process can be enhanced through the use of a Krylov subspace: instead of discarding the vectors $|{\psi }_r⟩$ through $H^{m-1}|{\psi }_r⟩$, they are used to construct a Krylov subspace. The matrix $H$ is projected into this subspace and the full spectrum is computed of the much smaller $\tilde{H}$. As in the case of time evolution, a single iteration of this algorithm is unlikely to converge within a reasonable error tolerance, again calling for a restarting scheme. A great deal of effort has gone into designing optimal restarting schemes for Krylov eigensolving; they are out of the scope of this work but we refer interested readers to the SLEPc technical papers for a discussion of the solvers underlying dynamite [HRT+07].

Perhaps surprisingly, since it solves for extremal eigenvalues, this technique can also be used to compute sets of neighboring excited states anywhere in the spectrum. This is achieved via a spectral transform: rather than performing the solve on the matrix $H$ itself, it is performed on a function of $H$ which moves the desired eigenvalues to the outside of the spectrum. In dynamite, the default spectral transform for finding interior eigenpairs is shift-and-invert, which, for a target energy $\sigma$, performs the transformation $H\to (H-\sigma {)}^{-1}$ such that the eigenvalues closest to $\sigma$ lie on the outside of the transformed matrix. This technique has been used with great success in the study of many-body localization in particular, as discussed in Sec. 2.3.2.

Finally, we note an important feature of all of these techniques: they almost exclusively rely on the matrix-vector product as the definition of the matrix, rather than ever needing to access individual matrix elements themselves. This fact is crucial for the ability to use “matrix-free” implementations with these algorithms, as discussing in Sec. 2.4.2 below.

Most operators that one encounters in the study of many-body spin physics are sparse: an overwhelming fraction of the matrix elements are zero. It is for this reason that dynamite uses the sparse numerical linear algebra libraries PETSc and SLEPc for its solvers. By default, the solvers use the compressed sparse row (CSR) storage format for matrices, in which only nonzero matrix elements (and the necessary indices) are stored in memory. For an operator of dimension $2^L\times 2^L$, this reduces the cost of storage (and matrix-vector multiplication) from $O\left(4^L\right)$ to $O\left(2^{L}k\right)$, where $k$ is the number of nonzero elements per row—and $k$ is polynomial in $L$ for virtually all physically relevant operators. Of course, $O\left(2^{L}k\right)$ still has the potential to be quite large, and frequently is in practice. For the case that even the sparse representation of an operator uses too much memory, dynamite implements a custom representation of operators which we call the “XZC” representation, that leverages the tensor product structure of spin chain operators to store them using as little as $O\left(k\right)$ memory.²² 2 The XZC representation is currently denoted “MSC” in various parts of the dynamite source code. The intuition is that the operator is represented symbolically as a decomposition into strings of Pauli matrices; instead of being stored, the elements of the large sparse matrix are computed “on-the-fly” as they are needed. In the field of numerical linear algebra, this type of technique is known as “matrix-free” computing; for historical reasons it is referred to as “shell matrices” in dynamite. It dramatically reduces the memory usage, enabling computations that wouldn’t be possible otherwise due to memory limitations—especially on GPUs, where VRAM is limited. Usually, this comes with a trade-off: the need to compute matrix elements on-the-fly adds to the computational cost, so the user pays for the reduction in memory usage with an increase in the runtime. However, CSR matrix operations are frequently limited by memory bandwidth rather than CPU usage—the speed is limited by how quickly the CPU can retrieve matrix elements from main memory. By reducing the amount of data transferred across the memory bus, it’s actually possible for matrix-free methods to be faster, while simultaneously decreasing the memory usage! Whether this is the case in practice depends on how fast the matrix elements can be generated in real time—and much effort in the development of dynamite has been directed at optimizing this process to an extreme. In Section 2.5, we compare the performance of shell and non-shell computations for various problems. Here we describe how it works, and how it interfaces with the solvers.

The fundamental idea behind the XZC representation is that any operator on a set of $L$ spins can be decomposed as a linear combination of products of single-site Pauli operators. Explicitly, for an operator on $L$ spins, which has $m$ terms in this decomposition, we may define a pair of $m\times L$ binary matrices $M^{\left(x\right)}$ and $M^{\left(z\right)}$, and a length-$m$ complex-valued vector $C$, to represent an operator $H$ as

where ${\sigma }_j^x$ and ${\sigma }_j^z$ are the Pauli $X$ and $Z$ operators on site $j$. In words, each row of $M^{\left(x\right)}$ and $M^{\left(z\right)}$, and corresponding element of $C$, corresponds to one term of the operator. In each row, the $1$’s in each binary matrix represents the locations in which a ${\sigma }^x$ or ${\sigma }^z$ Pauli should be present (with $0$’s corresponding to the identity Pauli on that site). ${\sigma }^y$ operators are represented via the product ${\sigma }_x{\sigma }_z$, when both $M_{ij}^{\left(x\right)}$ and $M_{ij}^{\left(z\right)}$ are equal to 1.

Because $L$ in practice never exceeds $64$, $M^{\left(x\right)}$ and $M^{\left(z\right)}$ are stored as vectors of integers.³³ 3 dynamite uses 32-bit integers by default, but can be configured to use 64-bit integers if $L$ is to exceed 31. The bits of each integer represent the spin indices for which a ${\sigma }_x$ and/or ${\sigma }_z$ operator should be included for that term. For the $C_i$, it is straightforward to see that each value must be either entirely real or entirely imaginary if the operator is Hermitian; this allows $C$ to be stored as a vector of real values (64-bit floating point numbers). All together, an operator with $m$ terms is stored as two length-$m$ vectors of integers plus one length-$m$ vector of real numbers, yielding an overall storage cost of $O\left(m\right)$. The quantity $k$ mentioned in the previous paragraph corresponds to the number of unique rows of $M^{\left(x\right)}$; in most situations $m$ is proportional to $k$.

We now describe how computations are implemented on this format—specifically matrix-vector multiplication, which is by far the most important matrix operation for the iterative solvers described in Section 2.4.1. By default, dynamite uses product states in the $Z$ basis as a basis for its representation of the Hilbert space. For the full Hilbert space of length $2^L$, dynamite takes basis state of index $\alpha$ to correspond to the product state $|\alpha ⟩={⨂}_j|{\alpha }_j⟩$, where ${\alpha }_j$ is the $j^{th}$ bit of the binary representation of the integer $\alpha$ (and with $|0⟩$ corresponding to an “up” spin, with ${\sigma }_z$ eigenvalue $+1$, and $|1⟩$ corresponding to a “down” spin, with eigenvalue $-1$). Then, an arbitrary state vector $|\psi ⟩$ is stored as an array $\overrightarrow{x}$ of $2^L$ complex numbers $x_{\alpha }$, such that $|\psi ⟩={\sum }_{\alpha }{x}_{\alpha }|\alpha ⟩$. This yields a very straightforward implementation of multiplication of operators by states. The operator-state multiplication can be decomposed as as the sum over multiplications by each of the $H_i$ of Eq. 2.14:

It is straightforward to see that the product of a single term $H_i$ from Eq. 2.14 with a product state $|\alpha ⟩$ yields another product state $|\beta ⟩$ multiplied by a coefficient:

where $M_{i,\ast }^{\left(x\right)}$ is the $i^{th}$ row of $M$. In words, the result product state is $|\beta ⟩$, where $\beta$ is the bitwise XOR of the integer $\alpha$ and the $i^{th}$ row of $M^{\left(x\right)}$, multiplied by the coefficient $C_i$, with a sign flip if an odd number of bits are set in the bitwise AND of the $i^{th}$ row of $M^{\left(z\right)}$ and $\alpha$.

This is very promising from an implementation perspective. As described above, the rows of $M^{\left(x\right)}$ and $M^{\left(z\right)}$, as well as the product states $\alpha$, will be represented as integers, this operation can be performed extremely quickly using bitwise operations. Computing $\beta =M_{i,\ast }^{\left(x\right)}⊻\alpha$ costs a single bitwise XOR. Computing the sign flip costs just a bitwise AND, a popcount operation to count the number of bits set⁴⁴ 4 popcount is a native machine instruction in most modern CPUs., and another bitwise AND to get the least significant bit of the result (the parity of the number of bits set). Finally flipping the sign of $C_i$ costs a single XOR instruction, and then a fused multiply-add (one machine instruction in modern CPUs) can be used to multiply $x_{\alpha }$ by the sign-flipped $C_i$ and add it into the appropriate element $y_{\beta }$ of the output vector. Overall, computing the action of each term of $H_i$ on each basis state costs two bitwise XORs, two bitwise ANDs, a popcount, and a fused multiply-add—not very expensive at all!

With that, we have a straightforward algorithm for computing a matrix-vector multiplication using the XZC representation of operators, given as pseudocode in Code Listing 2.7. The code listing includes one small addition that has not yet been discussed, which is the conversion of indices to states and vice versa. As previously mentioned, when using the full Hilbert space this operation corresponds to the identity—index $\alpha$ corresponds to product state $|\alpha ⟩$. However when the Hilbert space dimension is reduced into a symmetry subspace, this may not be the case. The functions state_to_idx and idx_to_state represent the conversion between indices and integers representing the corresponding product states.