LREI: A fast numerical solver for quantum Landau-Lifshitz equations

TL;DR

The paper introduces LREI, a low-rank numerical solver that significantly accelerates and reduces memory usage for simulating quantum spin dynamics governed by Landau-Lifshitz equations, enabling larger system simulations.

Contribution

It develops a low-rank, efficient numerical scheme for quantum Landau-Lifshitz equations that handles large spin systems by exploiting matrix structure and invariant subspaces.

Findings

Reduces computational complexity from cubic to quadratic in system size.

Decreases memory requirements from quadratic to linear in system size.

Enables simulation of systems with over a million density matrix elements.

Abstract

We develop LREI (Low-Rank Eigenmode Integration), a memory- and time-efficient scheme for solving quantum Landau-Lifshitz (q-LL) and quantum Landau-Lifshitz-Gilbert (q-LLG) equations, which govern spin dynamics in open quantum systems. Although system size grows exponentially with the number of spins, our approach exploits the low-rank structure of the density matrix and the sparsity of Hamiltonians to avoid full matrix computations. By representing density matrices via low-rank factors and applying Krylov subspace methods for partial eigendecompositions, we reduce the per-step complexity of Runge-Kutta and Adams-Bashforth schemes from $\mathcal{O}(N^3)$ to $\mathcal{O}(r^2N)$, where $N = 2^n$ is the Hilbert space dimension for $n$ spins and $r \ll N$ the effective rank. Similarly, memory costs shrink from $\mathcal{O}(N^2)$ to $\mathcal{O}(rN)$, since no full $N\times N$ matrices are formed. A key advance is handling the invariant subspace of zero eigenvalues. By using Householder reflectors built for the dominant eigenspace, we perform the solution entirely without large matrices. For example, a time step of a twenty-spin system, with density matrix size over one million, now takes only seconds on a standard laptop. Both Runge-Kutta and Adams-Bashforth methods are reformulated to preserve physical properties of the density matrix throughout evolution. This low-rank algorithm enables simulations of much larger spin systems, which were previously infeasible, providing a powerful tool for comparing q-LL and q-LLG dynamics, testing each model validity, and probing how quantum features such as correlations and entanglement evolve across different regimes of system size and damping.

Equations

$$\begin{array}{rl}& \dot{\rho }=\frac{\text{i}}{\hslash }[\rho ,H]-\frac{\kappa }{\hslash }[\rho ,[\rho ,H]]\\ & \rho (0)={\rho }_0\end{array}$$

$$\begin{array}{rl}& \dot{\rho }=\frac{\text{i}}{\hslash }[\rho ,H]+\text{i}\kappa [\rho ,\dot{\rho }]\\ & \rho (0)={\rho }_0\end{array}$$

$$\frac{d}{dt}{\lambda }_k(\rho (t))=0$$

$$\frac{d}{dt}\text{Tr}({\rho }^m)=0,m=1,2,\dots .$$

$$\rho (t)=\frac{\exp \left(-\frac{\text{i}}{\hslash }\tilde{H}t\right){\rho }_0\exp \left(\frac{\text{i}}{\hslash }\tilde{H}t\right)}{\mathrm{Tr}\left(\exp \left(-\frac{\text{i}}{\hslash }\tilde{H}t\right){\rho }_0\exp \left(\frac{\text{i}}{\hslash }\tilde{H}t\right)\right)},$$

$$\rho =V\Lambda V^{\ast }$$

$$X:=V^{\ast }\dot{\rho }V,$$

$$(I+S)X-XS=D,$$

$$D=\frac{\text{i}}{\hslash }{V}^{\ast }[\rho ,H]V$$

$$\begin{array}{rl}& \dot{\rho }=VX{V}^{\ast },\\ & \rho (0)={\rho }_0,\end{array}$$

$${\mathbf{x}}_{\ell }=\text{cwd}({\mathbf{d}}_{\ell },\mathbf{e}+\mathbf{s}-s_{\ell }\mathbf{e}),\ell =1,2,\dots ,N,$$

$$\mathrm{rank}({\rho }_0)=r$$

$$\mathrm{rank}(\rho (t))=r,t\ge 0.$$

$$\rho =\left[\begin{array}{cc}\tilde{V}& \hat{V}\end{array}\right]\left[\begin{array}{cc}\tilde{\Lambda }& \\ & 0\end{array}\right]\left[\begin{array}{c}{\tilde{V}}^{\ast }\\ {\hat{V}}^{\ast }\end{array}\right]$$

$$\rho =\tilde{V}\tilde{\Lambda }{\tilde{V}}^{\ast },$$

$$X=\text{cwd}(D,E+L-L^T)$$

$$L=\left[\begin{array}{cc}{L}_{11}& L_{12}\\ 0& 0\end{array}\right]$$

$$E+L-L^T=E+L+L^{\ast }=\left[\begin{array}{cc}{E}_{11}+L_{11}+L_{11}^{\ast }& E_{12}+L_{12}\\ (E_{12}+L_{12}{)}^{\ast }& E_{22}\end{array}\right].$$

$$D$$

$$=\frac{\text{i}}{\hslash }\left[\begin{array}{cc}\tilde{\Lambda }{\tilde{V}}^{\ast }H\tilde{V}-{\tilde{V}}^{\ast }H\tilde{V}\tilde{\Lambda }& \tilde{\Lambda }{\tilde{V}}^{\ast }H\hat{V}\\ -{\hat{V}}^{\ast }H\tilde{V}\tilde{\Lambda }& 0\end{array}\right]$$

$$=:\left[\begin{array}{cc}{D}_{11}& D_{12}\\ D_{12}^{\ast }& 0\end{array}\right]$$

$$X=\left[\begin{array}{cc}{X}_{11}& X_{12}\\ X_{12}^{\ast }& 0\end{array}\right]$$

$$\begin{array}{rl}& X_{11}=\text{cwd}(D_{11},E_{11}+L_{11}+L_{11}^{\ast }),D_{11}=\frac{\text{i}}{\hslash }(\tilde{\Lambda }{\tilde{V}}^{\ast }H\tilde{V}-{\tilde{V}}^{\ast }H\tilde{V}\tilde{\Lambda }),\\ & X_{12}=\text{cwd}(D_{12},E_{12}+L_{12}),D_{12}=\frac{\text{i}}{\hslash }\tilde{\Lambda }{\tilde{V}}^{\ast }H\hat{V}.\end{array}$$

$$\tilde{V}=QR=\left[\begin{array}{cc}\tilde{Q}& \hat{Q}\end{array}\right]\left[\begin{array}{c}\tilde{R}\\ 0\end{array}\right],$$

$$P_k=I-\frac{2}{{\mathbf{u}}_k^{\ast }{\mathbf{u}}_k}{\mathbf{u}}_k{\mathbf{u}}_k^{\ast },k=1,\dots ,r,$$

$$Q=V=P_1{P}_2\cdots P_r.$$

$$P=I-\beta \mathbf{u}{\mathbf{u}}^{\ast },\beta =\frac{2}{{\mathbf{u}}^{\ast }\mathbf{u}},$$

$$AP=A(I-\beta \mathbf{u}{\mathbf{u}}^{\ast })=A-\beta \mathbf{w}{\mathbf{u}}^{\ast },\mathbf{w}=A\mathbf{u},$$

$$PA=(I-\beta \mathbf{u}{\mathbf{u}}^{\ast })A=A-\beta \mathbf{u}{\mathbf{w}}^{\ast },\mathbf{w}=A^{\ast }\mathbf{u},$$

$$AV=A[\tilde{V}\hat{V}]=[A\tilde{V}A\hat{V}],$$

$$\hat{V}A=[\tilde{V}\hat{V}]\left[\begin{array}{c}0\\ A\end{array}\right]=:V{A}_0.$$

Full text

Keywords: Quantum Landau-Lifshitz equation, Quantum Landau-Lifshitz-Gilbert equation, Low-rank representations, Krylov subspace methods, Eigenvalue decomposition, Runge-Kutta methods, Adams-Bashforth methods.

Mathematics Subject Classification (2020): 65L06, 81Q05, 65F55.

1 Introduction

While the Landau–Lifshitz (LL) and Landau–Lifshitz–Gilbert (LLG) equations have long been central tools for modeling the microscopic dynamics of magnetic systems and materials, it is well understood that the underlying processes at the atomic scale are inherently quantum mechanical and governed by quantum dynamics. Crucially, intrinsically quantum features such as many-body correlations and entanglement, which are at the heart of modern quantum science and technology [1, 6, 9, 10] lie beyond the descriptive power of classical equations. Thus, to obtain a more accurate and fundamental understanding of the dynamics of many-body spin systems, it is necessary to study quantum-mechanical counterparts of the LL and LLG equations.

A quantum analogue of the LL equation was introduced in [19] based on a phenomenological open-system framework. The motivation was to derive spin dynamics from a more fundamental quantum perspective than that offered by the classical LL equation. The q-LL equation reads as

[TABLE]

where $H$ is the Hamiltonian, $\rho\in\mathbb{C}^{N\times N}$ is the density operator111In line with standard conventions in numerical linear algebra, we denote matrices by capital letters. The only exceptions are the variable $\rho$ (with or without subscripts), and the Pauli matrices $\sigma_{x}$, $\sigma_{y}$, and $\sigma_{z}$ where we retain the standard notation from physics texts., $\dot{\rho}=\frac{d\rho}{dt}$, $\hslash$ is the Planck constant, $\kappa$ is the dimensionless damping rate, and $\rho_{0}$ is the initial condition. Here, the commutator, $[A,B]$, of two operators $A$ and $B$, is defined to be $[A,B]=AB-BA$.

Building on this concept, in the recent paper [12] a quantum analog of the classical Landau–Lifshitz–Gilbert (LLG) equation was presented as

[TABLE]

The q-LL model (1.1) has a simpler damping term compared with the q-LLG equation (1.2), since the time derivative $\dot{\rho}$ does not appear on its right-hand side. In fact, the q-LL equation is an approximation of the full model (1.2), obtained by substituting $\dot{\rho}$ on the right-hand side of (1.2) with the whole right-hand side and truncating after the first term.

These models reveal several nontrivial quantum phenomena in systems composed of spin-$\frac{1}{2}$ particle pairs, including the emergence of spinless local states in antiferromagnetically coupled particles and emergence of long-lived highly entangled states, despite the presence of dissipation.

However, the studies in [19, 12] are limited to either pure state initial conditions or single- and two-spin systems, while quantum effects become increasingly intricate in general cases of initial conditions and/or many-spin settings. Although exact solutions exist for some special cases, a robust numerical solver is needed for studying the models in their general forms to further investigate their theoretical and practical aspects.

Two of the essential properties of the solution $\rho(t)$ of either equation (1.1) or equation (1.2) are that it is Hermitian and possesses the conservation of the spectrum, i.e.,

[TABLE]

where $\lambda_{k}\equiv\lambda_{k}(\rho(t))$, $k=1,\ldots,N$ are eigenvalues of $\rho(t)$. This means that the spectrum is independent of time; $\lambda_{k}(\rho(t))=\lambda_{k}(\rho(0))$ for all $t\geq 0$. As a consequence, the traces of all powers of $\rho$ are also conserved;

[TABLE]

The eigenvalues of $\rho(t)$ ($\equiv$ eigenvalues of $\rho_{0}$) are nonnegative, hence $\rho(t)$ remains positive semi-definite throughout the time integration.

Exact solutions of (1.1) and (1.2) exist only for pure (rank-one) initial states $\rho_{0}$. For the q-LL equation, the solution is

[TABLE]

where $\tilde{H}=\left(1-\text{i}\kappa\right)H$. The exact solution of the q-LLG equation in the pure-state case coincides with (1.4), up to a rescaling of time $t\mapsto t/(1+\kappa^{2})$. We use these closed-form solutions to assess the accuracy of our proposed algorithm in such special cases. However, the exact solutions involve computing matrix exponentials of the form $\exp(A)$ for complex matrices $A$, which in turn requires robust and efficient numerical solvers.

Any numerical method for solving the q-LL and q-LLG equations in their general form should ideally be robust, highly accurate, preserve the physical properties of the models, and scale to large spin systems. Numerical simulation becomes particularly challenging in many-body settings due to the exponential growth of the density matrix $\rho$, which is a full complex matrix of size $N\times N$ with $N=2^{n}$, where $n$ is the number of particles. This exponential scaling, combined with the inherent nonlinearity of the equations, poses severe computational challenges for simulating large quantum systems. For related numerical approaches to simulating large quantum systems, we refer the reader to [5, 18].

In [3], we proposed a numerical method for solving the q-LLG equation (1.2), which can be readily adapted to the q-LL equation (1.1). The method efficiently handles the complexities of quantum spin systems while preserving the key physical properties of the q-LLG dynamics. While it represents the first numerical solution of this equation, the approach relies on eigendecomposition of the evolving density matrix at each time step. As a result, the method becomes computationally expensive for large-scale systems due to the exponential growth of the matrix sizes with increasing numbers of spins.

In this paper we exploit the observation that in many practical scenarios the matrices and operators involved are initially of low rank and retain this property throughout the time evolution. By combining low-rank matrix representations with Krylov subspace methods and structure-preserving time integration, we develop an efficient computational method which reduces the memory usage and the time complexity, significantly. This enables the simulation and study of much larger quantum spin networks.

The remainder of the paper is organized as follows. In Section 2, we briefly review the method of [3]. Section 3 introduces our new method for solving the q-LLG equation (1.2), and Section 4 demonstrates how this method can be adapted to the q-LL equation (1.1). In fact, we first address the more complex q-LLG equation and then tailor the approach to the simpler q-LL model. Interestingly, from a computational perspective, both equations exhibit essentially the same complexity. Section 5 presents numerical experiments, and Section 6 concludes with a summary and directions for future research.

2 Available numerical solution

In [3], a conservative modification of explicit Runge-Kutta methods was developed for solving the q-LLG equation (1.2). The approach is briefly outlined below.

We starting with the eigendecomposition

[TABLE]

where $V\in\mathbb{C}^{N\times N}$ is a unitary matrix whose columns are the eigenvectors of $\rho$, and $\Lambda\in\mathbb{R}^{N\times N}$ is a diagonal matrix containing the eigenvalues222Here $V^{*}$ denotes the Hermitian conjugate (adjoint) of $V$. This convention differs from that in standard quantum information texts such as [13], where the star symbol refers only to the element-wise complex conjugate of a matrix.. The change of variables

[TABLE]

transforms the original equation into a Sylvester-type equation of the form

[TABLE]

where

[TABLE]

depends nonlinearly on $\rho$ and $S=-\text{i}\kappa\Lambda$. The matrix $X$ is computed from (2.3) and used to evolve $\rho$ via

[TABLE]

instead of the original form (1.2). Notably, the right-hand side $VXV^{*}$ is now independent of $\dot{\rho}$ and depends only on $\rho$ through its spectral decomposition.

The solution to the Sylvester-type equation (2.3) is given column-wise as

[TABLE]

where $\bm{x}_{\ell}$ and $\bm{d}_{\ell}$ denote the $\ell$-th columns of $X$ and $D$, respectively, $\bm{s}$ is the diagonal of $S$; and $\bm{e}=[1,1,\ldots,1]^{T}$. The operator $\text{cwd}(\bm{x},\bm{y})$ denotes the componentwise division of vectors $\bm{x}$ and $\bm{y}$. The solution is well-defined because the matrices $I+S$ and $-S$ have no eigenvalues in common, which, in turn, follows from the fact that $\rho$, being Hermitian, has no eigenvalues equal to ${-\text{i}}/{2\kappa}$.

Explicit Runge-Kutta (RK) methods were then employed to integrate equation (2.5) in time, as described in [3]. To preserve the spectrum of the density matrix during the time integration, the approach in [3] projects all the intermediate RK states back onto the manifold of matrices that share the same spectrum as the initial density matrix $\rho_{0}$ by replacing the evolving eigenvalue matrix $\Lambda(t)$ with the fixed initial eigenvalue matrix $\Lambda(0)$ at each stage of the RK method. We refer to such techniques as the eigenmode integration (EI) methods for brevity. For instance, when combined with RK schemes, we refer to them as EI-RK methods.

2.1 Computational bottleneck

The dominant computational cost of EI methods comes from the repeated eigendecompositions of the full density matrix $\rho\in\mathbb{C}^{N\times N}$ at each time step. Specifically, for an $s$-stage Runge-Kutta method, $s$ eigendecompositions are required per time step. Given that $N=2^{n}$, where $n$ is the number of spins, and each decomposition costs $\mathcal{O}(N^{3})$, the total cost at each time step grows exponentially with $n$, more precisely $\mathcal{O}(sN^{3})$. As a result, both the time and memory complexities of the algorithm increase significantly as the system size grows. This is a major computational bottleneck for large spin systems. The aim of this paper is to benefit from low-rank structure of the density matrix $\rho$ and the sparsity of other matrices involved to reduce the time complexity to $\mathcal{O}(sr^{2}N)$, and the space complexity to $\mathcal{O}(srN)$, where $r\ll N$ is the rank of $\rho$.

3 Low-rank eigenmode integration (LREI) scheme for q-LLG

In typical physical systems, such as spin chains, lattices, and quantum many-body models, the Hamiltonian $H$ is a sparse matrix due to local interactions and is generally full-rank in the absence of special symmetries or fine-tuning. Moreover, in simulations and experiments—particularly in quantum dynamics—the initial density matrix $\rho_{0}$ is often low-rank, either inherently or as a result of renormalization or truncation, especially at low temperatures or when representing only a small subset of relevant eigenstates. We assume that

[TABLE]

with $r\ll N$. Owing to the structure-preserving nature of the q-LLG evolution, the rank of $\rho(t)$ remains unchanged and satisfies

[TABLE]

We benefit from the sparsity of $H$ and mostly the low-rank structure of $\rho(t)$ to develop an efficient algorithm, which we refer to as the low-rank eigenmode integration (LREI) method. In this section, we first develop the LREI method for solving the q-LLG equation and then adapt it to the q-LL equation in Section 4.

3.1 Partial eigendecomposition

The eigendecomposition (2.1) can be rewritten as

[TABLE]

where ${\widetilde{V}}\in\mathbb{C}^{N\times r}$ contains orthonormal columns corresponding to the $r$ positive eigenvalues of $\rho$, and ${\widetilde{\Lambda}}\in\mathbb{R}^{r\times r}$ is a diagonal matrix with these eigenvalues sorted in a decreasing order on its diagonal. In the quantum computing language, ${\widetilde{V}}$ is called the support of operator $\rho$. The partial eigendecomposition

[TABLE]

can be computed efficiently using Krylov subspace iterative methods, which are well-suited for large-scale Hermitian matrices; see Section 3.7 for more details.

However, the contribution of the complementary subspace, associated with the matrix ${\widehat{V}}\in\mathbb{C}^{N\times(N-r)}$ also plays a significant role in the numerical procedure described above. In Section 3.4, by relying only on partial decomposition (3.1), we will effectively handle the contribution of big matrix ${\widehat{V}}$ without explicitly forming and storing it.

3.2 Space complexity

To manage the space (memory) complexity of the algorithm, we avoid storing or manipulating any full matrices of size $N\times N$ or $N\times(N-r)$, since $N=2^{n}$ grows exponentially with the number of spins $n$. This applies in particular to matrices such as $\rho$, ${\widehat{V}}$, $D$, and $X$.

According to the underlying theory, the spectrum $\Lambda$ of $\rho$ remains invariant throughout the time integration. To enforce this property in our numerical scheme, we adopt the same projection method proposed in [3]. As a result, instead of explicitly working with $\rho$, we store and evolve only its factor ${\widetilde{V}}$ of size $N\times r$ throughout the computation. Furthermore, in the remainder of this section, we explain in detail how to completely avoid manipulating the large matrices ${\widehat{V}}$, $D$, and $X$.

3.3 Efficient computation of $X$

Equations (2.6) can be written in a matrix form as

[TABLE]

where $D$ is defined in (2.4), $E$ denotes the $N\times N$ matrix of all ones, and $L=-\text{i}\kappa[\bm{\lambda}\;\bm{\lambda}\;\ldots\;\bm{\lambda}]\in\mathbb{C}^{N\times N}$, where $\bm{\lambda}=[\lambda_{1},\ldots,\lambda_{r},0,\ldots,0]^{T}$ is the vector of all eigenvalues of $\rho$. Here, the operator $\text{cwd}(A,B)$ denotes the componentwise division of matrices $A$ and $B$. We also note that according to the structure of $L$ we simply have $-L^{T}=L^{*}$. On the other hand, we can split the matrix $L$ into four blocks as

[TABLE]

where $L_{11}$ is an $r\times r$ matrix and $L_{12}$ is of size $r\times(N-r)$. If we split the all-ones matrix $E$ in the same way, we can write the denominator as

[TABLE]

In addition, from (2.4) the numerator $D$ has the representation

[TABLE]

where the top-left block $D_{11}$ is of size $r\times r$ and $D_{12}$ is of size $r\times(N-r)$. This shows that the major part of $D$ (the fourth block) is zero. In the computation above we used the facts that ${\widetilde{V}}^{*}{\widetilde{V}}=I$, ${\widetilde{V}}^{*}{\widehat{V}}=\bm{0}$, ${\widehat{V}}^{*}{\widetilde{V}}=\bm{0}$, and $H=H^{*}$.

Consequently, the solution $X$ can be similarly split into four blocks as

[TABLE]

where

[TABLE]

The computation of $X_{11}$ is relatively inexpensive due to the sparsity of the large matrix $H$, the small size of ${\widetilde{\Lambda}}$, $L_{11}$ and $E_{11}$ (all $r\times r$), and the tall-skinny structure of ${\widetilde{V}}$.

The main computational challenge lies in computing $X_{12}$, as it involves the dense and high-dimensional matrix ${\widehat{V}}$ through the computation of $D_{12}$. To address this, we avoid explicitly constructing ${\widehat{V}}$ and instead employ an efficient technique for computing $D_{12}$ implicitly, as detailed below.

3.4 Efficient computation of products $A{\widehat{V}}$ and ${\widehat{V}}\!A$

In this section, we describe how to efficiently compute matrix products of the form $A{\widehat{V}}$, which are required for evaluating $D_{12}$, and ${\widehat{V}}\!A$ which will later be needed for assembling the right-hand side of the ODE (2.5). Recall that the matrices ${\widetilde{V}}$ and ${\widehat{V}}$ contain eigenvectors corresponding to the positive and zero eigenvalues, respectively, of the positive semi-definite matrix $\rho$. As such, the columns of ${\widehat{V}}$ form the orthonormal complement of the columns of ${\widetilde{V}}$. One common approach to computing ${\widehat{V}}$, without having it at hand, is via the full QR factorization of the $N\times r$ matrix ${\widetilde{V}}$:

[TABLE]

where $Q$ is a unitary matrix with $\widetilde{Q}\in\mathbb{C}^{N\times r}$, $\widehat{Q}\in\mathbb{C}^{N\times(N-r)}$, and $\widetilde{R}\in\mathbb{C}^{r\times r}$ is upper triangular. By normalizing the diagonal entries of $\widetilde{R}$ to be real positive numbers, we can ensure that $\widetilde{Q}={\widetilde{V}}$ (up to a column permutation), and hence $\widehat{Q}$ provides an orthonormal basis for the orthogonal complement, i.e., ${\widehat{V}}=\widehat{Q}$. The full unitary matrix $Q=[\widetilde{Q}\;\;\widehat{Q}]$ therefore coincides with the complete eigenvector matrix $V$.

The QR factorization of ${\widetilde{V}}$ via Householder algorithm involves $r$ steps, each using a Householder transformation

[TABLE]

where each $P_{k}$ is Hermitian and unitary [7, 16]. These matrices are applied sequentially to introduce zeros below the diagonal in the $k$-th column of ${\widetilde{V}}$, eventually producing the upper triangular matrix $R$. The full matrix $Q=V$ is then expressed as the product of Householder reflectors, i.e.,

[TABLE]

To compute matrix products such as $AV$ or $V\!A$ without explicitly forming the large matrix $V$, we can exploit the special structure of the Householder matrices. Instead of constructing each $P_{k}$ in full, we work directly with the Householder vectors $\bm{u}_{k}$ that define them via (3.3) to have an efficient and memory-saving implementations of such matrix-matrix products.

More precisely, instead of performing a full QR factorization, we compute only the $r$ Householder vectors $\bm{u}_{k}$, for $k=1,2,\ldots,r$. This requires approximately $2Nr^{2}$ floating-point operations (flops), which is computationally cheap when $r$ is small. To compute matrix products of the form $AV=AP_{1}P_{2}\cdots P_{r}$, one would naively perform $r$ successive matrix-matrix multiplications of the form $AP$. For a general $m\times N$ matrix $A$ and a dense $N\times N$ matrix $P$, this would cost approximately $2mN^{2}$ flops. However, by using the structure of the Householder matrix

[TABLE]

we can compute the product efficiently without forming $P$ explicitly. We use the identity

[TABLE]

which requires only $4mN$ flops. A similar strategy applies for computing $V\!A=P_{1}P_{2}\cdots P_{r}A$ for a matrix $A\in\mathbb{C}^{N\times m}$. Using the same formulation,

[TABLE]

also costs for $4mN$ flops. In our algorithm, $m=r$, and since this product must be repeated $r$ times, the total cost of computing $AV$ or $V\!A$ becomes $4Nr^{2}$. However, we are only interested in the portion involving ${\widehat{V}}$, the orthogonal complement of ${\widetilde{V}}$. For computing $A{\widehat{V}}$, observe that

[TABLE]

so we compute $AV$ and extract only the last $N-r$ columns corresponding to $A{\widehat{V}}$.

To compute ${\widehat{V}}\!A$, where $A\in\mathbb{C}^{(N-r)\times r}$, we augment $A$ by placing a zero block of size $r\times r$ on top and obtain

[TABLE]

Then the product can be computed using the same Householder-based approach discussed above.

Consequently, the total cost of computing the products $A{\widehat{V}}$ for $A\in\mathbb{C}^{r\times N}$, and ${\widehat{V}}\!A$ for $A\in\mathbb{C}^{r\times(N-r)}$, is approximately $6Nr^{2}$. Since $r$, the rank of $\rho$, is fixed and typically small, the overall time complexity becomes $\mathcal{O}(N)$, which is optimal in terms of $N$, the size of the system. Importantly, throughout this process, we never form the large matrix ${\widehat{V}}$ or the full Householder matrices $P_{k}$ explicitly. Instead, we operate directly on the small portion ${\widetilde{V}}$ and the compact Householder vectors $\bm{u}_{k}$, to reduce both memory and computational costs, significantly. The space complexity is clearly of order $rN$.

3.5 Computation of the ODE right-hand side $VXV^{*}$

All ODE solvers used to integrate (2.5) require efficient evaluation of the right-hand side matrix $VXV^{*}$ at each iteration. Although this matrix is of size $N\times N$, we avoid forming it explicitly by working only with its low-rank ingredients. Suppose that, at a given time step, we are provided with the partial factor ${\widetilde{V}}$ and positive part of the spectrum, ${\widetilde{\Lambda}}$, of $\rho$. As described in the previous subsections, this information suffices to compute the submatrices $X_{11}\in\mathbb{C}^{r\times r}$ and $X_{12}\in\mathbb{C}^{r\times(N-r)}$ of the matrix $X$, with time and memory complexities of $\mathcal{O}(r^{2}N)$ and $\mathcal{O}(rN)$, respectively. Now, employing the block representations, we can write

[TABLE]

where

[TABLE]

The computation of $W$ involves the orthogonal complement matrix ${\widehat{V}}$, but as shown in the previous subsection, we avoid constructing ${\widehat{V}}$ explicitly. Instead, we compute $W$ efficiently using the Householder-based approach with a cost of $\mathcal{O}(Nr^{2})$. The first term ${\widetilde{V}}X_{11}$ of $Z$ is inexpensive to compute since $X_{11}$ is a small $r\times r$ matrix.

Thus, rather than forming the full $N\times N$ matrix $VXV^{*}$, we return the pair of matrices $Z\in\mathbb{C}^{N\times r}$ and $W\in\mathbb{C}^{N\times r}$, and express the time derivative as

[TABLE]

We indeed write the right-hand side of the ODE in terms of its low-rank components ${\widetilde{V}}$, $Z$, and $W$. We will see in Subsections 3.7 and 3.8, how this specific representation can be efficiently used in numerical ODE solvers.

3.6 Memory-efficient representation of $\rho(t)$

Since the spectrum of $\rho(t)$ is preserved over time, we have

[TABLE]

where ${\widetilde{\Lambda}}(0)=:{\widetilde{\Lambda}}_{0}$ denotes the fixed diagonal matrix of positive eigenvalues of $\rho(0)$. This allows us to avoid storing or evolving the full $N\times N$ matrix $\rho(t)$. Instead, we restrict the evolution to the low-rank factor ${\widetilde{V}}(t)$, which yields substantial dimensionality reduction without loss of any information. In what follows, we show how physical observables such as energy, magnetization, and entanglement measures can be computed efficiently using only ${\widetilde{V}}(t)$, without forming $\rho(t)$ explicitly.

For the energy and magnetization, we compute the corresponding expectation values

[TABLE]

where $A$ represents either the Hamiltonian $H$ (for energy), or the magnetization operators $M_{v}$, for $v=x,y,z$ (for magnetization). The magnetization operator along direction $v$ is defined by

[TABLE]

with $S_{i}^{v}$ denoting the spin operator at site $i$ in the $v$-direction,

[TABLE]

where $I$ is the identity matrix, and $\sigma_{x},\sigma_{y},\sigma_{z}\in\mathrm{SU}(2)$ are the Pauli matrices [14].

Although both $H$ and $M_{v}$ are large $N\times N$ matrices, their sparse structure allows them to be stored in memory-efficient formats. Using the factorized form of $\rho(t)$, the expectation value becomes

[TABLE]

where the last equality follows from the cyclic property of the trace. The parentheses are deliberately placed to enforce an optimized order of matrix multiplication. We first multiply the sparse matrix $A$ by the tall-skinny matrix ${\widetilde{V}}$, yielding an $N\times r$ matrix. This result is then left-multiplied by ${\widetilde{V}}^{*}$ to produce a small $r\times r$ matrix, which is finally multiplied by the diagonal matrix ${\widetilde{\Lambda}}_{0}$ of size $r\times r$. This ordering minimizes both memory usage and computational cost by replacing the full trace computation with one of complexity $\mathcal{O}(Nr^{2})$, where $r\ll N$.

To compute the entanglement between two spins located at sites $k$ and $l$, the measurement is no longer linear with respect to the density operator and cannot be directly expressed in the simple form of an expectation value as in (3.6). In fact, two-spin entanglement can be quantified by measures such as concurrence [20], negativity [17], or nonlocality [8], all of which involve non-linear operations on the reduced two-spin density matrix $R_{kl}(t)\in\mathbb{C}^{4\times 4}$ of the global state $\rho(t)$. This reduced density matrix is obtained via the partial trace

[TABLE]

where the trace is taken over all degrees of freedom (i.e., basis states) corresponding to spins not in $\{\mathbf{S}_{k},\mathbf{S}_{l}\}$, with each $\mathbf{S}_{i}=[S_{i}^{x},S_{i}^{y},S_{i}^{z}]$ denoting the spin operator vector at site $i$ [14]333A boldface capital letter, such as $\mathbf{S}$, denotes a horizontal concatenation of three matrices.. The computation of the 16 entries of the matrix $R_{kl}$ reduces to evaluating expressions of the form $\text{Tr}(A\rho(t))$ for certain sparse observables $A$, similar to the energy and magnetization computations described earlier. Subsequent entanglement measures (such as concurrence, negativity, and locality) are then computed directly from the small $4\times 4$ matrix $R_{kl}(t)$.

3.7 Explicit one-step schemes

Explicit one-step ODE solvers start with the initial density matrix $\rho_{0}$ and compute an approximate solution $\rho_{k+1}\approx\rho(t_{k+1})$ at each time level $t_{k+1}=(k+1)h$ using a recurrence of the form

[TABLE]

where $h=t_{k+1}-t_{k}$ denotes the time step and $\psi$ is a function of $\rho_{k}$ and $h$ depending on the right-hand side of the ODE. The iteration continues until the final time $t_{F}$ is reached. For example, in Euler’s method $\psi(\rho;h)=Z{\widetilde{V}}^{*}+{\widetilde{V}}W^{*}=:f(\rho)$ and for a Runge-Kutta method of order two (RK2) $\psi(\rho;h)=\frac{1}{2}[f(\rho)+f(\rho+hf(\rho))]$ which can be reformulated as a two-stage scheme

[TABLE]

where the function $f(K)$ is given by

[TABLE]

with ${\widetilde{V}}$ being the partial eigenvector matrix of $K$ corresponding to its nonzero eigenvalues. The matrices $Z$ and $W$ are computed as described previously, by evaluating the right-hand side of the evolution equation with $\rho$ replaced by $K$.

In general, an explicit $s$-stage RK method has the form

[TABLE]

for some real coefficients $a_{j,\ell}$ and $b_{\ell}$. For example, in the above RK2 scheme, $a_{2,1}=1$ and $b_{1}=b_{2}=\frac{1}{2}$. We refer the reader to [4] for details and the list of coefficients for higher order RK methods. At each stage, we compute $f(K_{j})=Z_{j}{\widetilde{V}}_{j}^{*}+{\widetilde{V}}_{j}W_{j}^{*}$ and update the stage to obtain the new matrix $K_{j+1}$ via

[TABLE]

where

[TABLE]

and finally

[TABLE]

where the new matrices $A_{\ell}$ and $B_{\ell}$ are defined by

[TABLE]

We first point out that this scheme preserves the spectrum of $\rho_{k}$ for all $k=0,1,\ldots$, since the initial eigenvalue matrix ${\widetilde{\Lambda}}_{0}$ is consistently used in place of ${\widetilde{\Lambda}}_{1}$ in the second row of equations (3.8) and (3.9). This guarantees that the eigenvalues of $\rho_{k}$ remain unchanged throughout the evolution.

We also do not store the full matrices $K_{j}$ and $\rho_{k+1}$ as dense $N\times N$ arrays. We also do not store the full matrices $K_{j}$ and $\rho_{k+1}$ as dense $N\times N$ arrays. Instead, we use them in their low-rank representations (not approximation!) as in equations (3.8) and (3.9), involving matrices $A_{\ell}$, $B_{\ell}$, and ${\widetilde{V}}_{\ell}$, each of size $N\times r$. Note that $\rho_{k+1}$ computed at time step $k$ serves as $K_{1}$ for the next time step $k+1$.

Given the sparse initial density matrix $\rho_{0}$ and its rank $r=\operatorname{rank}(\rho_{0})$, we compute its partial eigendecomposition

[TABLE]

and start the RK iterations. Each iteration $k$ contains $s$ stages. At each stage $j+1$, we must compute the partial eigendecomposition of the matrix $K_{j+1}$. This is performed efficiently using iterative Krylov subspace methods, more specifically the Implicitly Restarted Lanczos Method (IRLM) [11, 15].

The idea is to build a length-$m$ Lanczos decomposition

[TABLE]

where $m$ is the dimension of the Krylov subspace

[TABLE]

to project onto, $U_{m}$ is an $N\times m$ matrix whose columns form an orthogonal basis for the Krylov subspace, $e_{m}$ denotes the $m$-th column of the $m\times m$ identity matrix, $\bm{f}_{m}$ is the residual vector which is orthogonal to the columns of $U_{m}$, and $T_{m}=U_{m}^{\ast}K_{j+1}U_{m}$ is tridiagonal. The eigenvalues of the small $m\times m$ matrix $T_{m}$, known as the Ritz values, approximate those of $K_{j+1}$, and its eigenvectors when multiplied by $U_{m}$, yield the corresponding approximate eigenvectors of $K_{j+1}$.

Two points merit highlighting in practice. First, to obtain $r$ eigenvalues of $K_{j+1}$ to high-accuracy, it is enough to project onto a Krylov subspace with dimension as small as $m=2r+1$, which is far less than $N$. Second, the method does not necessarily require the full matrix $K_{j+1}$ to be explicitly formed [11, 15]. Instead, these methods only require access to a function that computes matrix-vector products $K_{j+1}\bm{x}$ for arbitrary vectors $\bm{x}\in\mathbb{C}^{N}$. This allows us to supply the Krylov solver with a function $g:\mathbb{C}^{N}\to\mathbb{C}^{N}$, defined by

[TABLE]

where all matrix-vector multiplications are performed using only the low-rank components. The parentheses in the expression for $g(\bm{x})$ are deliberately placed to avoid forming full $N\times N$ matrices, instead, allowing only matrix-vector products involving matrices of size $N\times r$. There are $m$ Lanczos steps, each involving matrix-vector products whose cost is $\mathcal{O}(Nr)$, summing to $\mathcal{O}(mNr)=\mathcal{O}(Nr^{2})$ operations which is similar to the cost of orthonormalization. The resulting $m\times m$ tridiagonal eigenvalue problem requires $\mathcal{O}(m^{2})=\mathcal{O}(r^{2})$ operations for the Ritz values which is negligible for $r\ll N$ and forming $r$ Ritz vectors adds $\mathcal{O}(Nr^{2})$. Overall, the time complexity of each execution of the Lanczos method remains $\mathcal{O}(Nr^{2})$ and its storage is $\mathcal{O}(Nm)=\mathcal{O}(Nr)$.

Similarly, we avoid forming the full matrix $\rho_{k+1}$ and instead work directly with its low-rank representation (3.9), just as we did for the intermediate stages $K_{j+1}$. From this representation, we compute its partial eigenvector matrix ${\widetilde{V}}_{1}$, which serves as the representative of $K_{1}$ matrix for the next time step. This iteration is repeated until the final time $t_{F}$ is reached.

We refer to this class of methods as LREI-RK schemes to highlight their foundation in low-rank eigendecomposition combined with Runge-Kutta time integration.

3.8 Explicit multi-step schemes

Runge–Kutta methods belong to the class of single-step schemes, where at each time step $k+1$, the solution $\rho_{k+1}$ depends solely on the previous value $\rho_{k}$. However, each RK step involves $s$ intermediate stages, each requiring one partial eigendecomposition.

In contrast, $m$-step methods compute $\rho_{k+1}$ based on $m$ previous steps $\rho_{k},\rho_{k-1},\ldots,\rho_{k-m+1}$. Among these, the Adams–Bashforth (AB) methods form a class of explicit schemes with the general form

[TABLE]

given initial conditions $\rho_{0},\rho_{1},\ldots,\rho_{m-1}$. For example, the coefficients of the 2-step AB method are $b_{1}=\frac{3}{2}$, $b_{2}=-\frac{3}{2}$. Refer to [2] for more details and the list of higher order AB methods. AB methods give an accuracy of order $\mathcal{O}(h^{m})$, but they typically have relatively small regions of absolute stability. To initialize the method, the starting values $\rho_{1},\ldots,\rho_{m-1}$ can be computed using an RK method of global order $m-1$ (local order $m$).

The main advantage of an AB method of order $m$ over an RK method of the same order is a more computational efficiency in the sense that each AB step requires only one partial eigendecomposition, whereas an RK method of the same order requires $m$ such decompositions per step. However, RK methods have larger stability regions and usually give more accurate results for a fixed time step $h$.

For an efficient implementation, similar to (3.9) we can write $\rho_{k+1}$ in terms of its low-rank components as

[TABLE]

where ${\widetilde{V}}_{\ell}$, $Z_{\ell}$, $W_{\ell}$ are the low-rank components derived from previous steps $\ell=k-m+1,\ldots,k$. The same efficient strategy, outlined earlier for RK methods, is applied for storage, manipulation, and eigendecomposition in the AB schemes. These schemes are refereed to as LREI-AB methods.

4 Solving q-LL equation

At first glance, one might expect solving the q-LL equation to be computationally cheaper than solving the q-LLG equation. However, if we manage to have numerical schemes that respect the physical invariants, and/or we aim to simulate for large spin systems (e.g., more than 13) by employing the low-rank representations, solving the q-LL equation exhibits roughly the same time and space complexities as the q-LLG equation. More precisely, in the LREI process for the q-LL equation, we do not need to solve the Sylvester equation. The only difference then lies in the computation of the matrix blocks of $X$, which are now given by

[TABLE]

A comparison of (3.2) and (4.1) shows that the computational cost of computing $X_{11}$ and $X_{12}$ shares the same dominant matrix-matrix operations resulting in a complexity of the same order. In particular, solving the Sylvester equation in the q-LLG case does not introduce any additional cost for the LREI algorithm. All other components of the algorithms (including the implementation in RK and AB settings) are exactly the same for both the q-LL and q-LLG equations.

5 Numerical results

The code was implemented in MATLAB and executed on a MacBook Pro with an Apple M1 Max chip and 64 GB of RAM. As a Krylov subspace method for computing partial eigendecompositions we use the eigs function of MATLAB which internally uses ARPACK, an iterative Arnoldi/Lanczos method [11]. This function can also work with function handles such as $g$ in (3.10).

The spin system we consider here is governed by the spin Hamiltonian

[TABLE]

where $J\in\mathbb{R}$ is the isotropic Heisenberg exchange interaction, $\bm{d}_{ij}\in\mathbb{R}^{3}$ ($\bm{d}_{ij}=-\bm{d}_{ji}$) represents the Dzyaloshinskii–Moriya interaction (DMI), and $\bm{b}\in\mathbb{R}^{3}$ is a uniform external magnetic field. The gyromagnetic constant is given by $\mu=-\frac{\mu_{B}g}{\hslash}$, where $\mu_{B}=5.8\times 10^{-2}\,\text{meV/T}$ is the Bohr magneton and $g=2$ is the Landé $g$-factor. The spin operators $\mathbf{S}_{i}=[S_{i}^{x},S_{i}^{y},S_{i}^{z}]$ are spin-$\frac{1}{2}$ operators defined using the Pauli matrices as described in (3.5). An exact solution to the system is known when the initial density matrix $\rho_{0}$ has rank one; see [3]. In all examples, the fixed dissipation parameter $\kappa=0.5$ will be used.

5.1 Accuracy test

First we reproduce the results of [3] using the LREI methods to verify their accuracy for a rank-one initial density matrix for which the exact solutions are available. We consider an antiferromagnetically ordered pure state $\rho_{0}=|\text{AF}_{1}\rangle\langle\text{AF}_{1}|$ as such a rank-one density matrix444 The Dirac notation $|\psi\rangle$ denotes a column vector and $\langle\psi|$ its conjugate transpose (row vector). Consequently, the outer product $|\psi\rangle\langle\psi|$ is a rank-1 matrix. In particular, the state $|\text{AF}_{1}\rangle$ corresponds to a vector with a single nonzero entry at the index $i=({0101}\ldots)_{2}+1,$ where the binary string of length $n$ alternates between [math] and $1$. .

Following the setup in [3], we consider a triangular lattice of 16 periodic spins (with 9 effective degrees of freedom) with Hamiltonian parameters $J=1$ meV, $\bm{d}_{ij}=-\bm{d}_{ji}=\|\bm{d}\|(0,0,1)$, $\|\bm{d}\|_{2}=0.4$ meV, and $\bm{b}=\|\bm{b}\|_{2}(1,0,0)$, where $\|\bm{b}\|_{2}=1$ T is the strength of magnetic field vector. Given the pure state initial condition $\rho_{0}=|\text{AF}_{1}\rangle\langle\text{AF}_{1}|$, Figure 1 presents the error plots for LREI-RK and LREI-AB methods in terms of the time step size $h$ for both q-LL and q-LLG equations.

As seen in the plots, RK methods achieve their expected convergence orders $m=1,2,3,4$ when $h\lesssim 0.05$. The AB methods also display their theoretical orders $m=2,3,4$, but only for smaller step sizes, here $h\lesssim 0.0125$. This delay is attributed to the relatively smaller stability regions of high-order AB schemes. For sufficiently small $h$, both RK$m$ and AB$m$ methods attain their theoretical convergence orders $\mathcal{O}(h^{m})$. However, for the same order, RK methods tend to produce more accurate results than their AB counterparts, albeit at the cost of requiring multiple partial eigendecompositions per time step.

We also emphasize that the construction of LREI methods guarantees that both RK and AB schemes introduced in Sections 3.7 and 3.8 are conservative, in the sense that the spectrum of $\rho_{k}$ is preserved at each time step. As a result, physical properties such as the non-negativity of $\rho_{k}$ and the trace of all powers $\text{Tr}(\rho_{k}^{\ell}),\,\ell=1,2,3,\ldots$ are preserved throughout the time evolution.

5.2 Complexity test

We now compare the CPU time required by the original method presented in [3] (the EI method) and LREI method developed in this work. To evaluate performance, both algorithms are applied to spin clusters of different sizes, with spin number $n$ ranging from $2$ to $24$. In all cases, we use triangular lattice configurations (either periodic or non-periodic) and the Hamiltonian defined in (5.1), with parameters identical to those used in the previous accuracy tests.

We focus on the RK4 and AB4, and perform a comparison based solely on execution times, since both low-rank and full-rank algorithms yield approximately the same accuracy in cases where the later one is computationally feasible. In Figure 2 we report the average CPU time required to complete ‘one’ time step of the time evolution of q-LLG dynamics including computation of the density operator and the physical quantities such as energy, magnetization, and entanglement. For the new algorithm we report the results for initial conditions with different ranks ranging from $r=1$ to $4$ while the original EI method is tested only for an initial data of rank $1$, as its performance is approximately independent of matrix rank. To construct initial conditions with different ranks, we use weighted combinations of rank-1 state $|\text{AF}_{1}\rangle\langle\text{AF}_{1}|$ (as defined above), the antiferromagnetic state $|\text{AF}_{2}\rangle\langle\text{AF}_{2}|$ with the index $i=({1010...})_{2}+1$, the GHZ-State (Greenberger–Horne–Zeilinger) $|\text{GHZ}\rangle\langle\text{GHZ}|$, and the W-state $|\text{W}\rangle\langle\text{W}|$. See [3] for definitions.

Figure 2 illustrates that while both algorithms execute in a fraction of a second for small systems (up to 5 spins), the low-rank approach demonstrates a clear and substantial computational advantage as the number of spins increases. For instance, with $n=13$ spins, the LREI algorithm completes in less than one second, whereas the standard method requires approximately $10^{4}$ seconds. Moreover, the EI algorithm fails to execute for systems with more than $13$ effective spins due to memory-bound and computational constraints, while the LREI method reliably scales up to $n=24$ effective spins ($35$ periodic spins) on our personal laptop. Another notable observation is that increasing the matrix rank has only a modest and manageable impact on the computational complexity of the LREI algorithm. In practice, we observed linear growth in $r$, whereas the theoretical bound predicts a quadratic dependence.

Figure 3 compares the runtimes per time step of the LREI-RK4 and LREI-AB4 algorithms, tested on an initial density matrix of rank $r=3$. The LREI-AB4 method achieves an average speedup of about $2\times$ for $n>10$. This is because LREI-RK4 requires four partial eigendecompositions per step, whereas LREI-AB4 needs only one. However, additional shared matrix–vector and matrix–matrix operations reduce the theoretical speedup from $4\times$ to $2\times$. This does not imply that LREI-AB schemes are preferable. As discussed earlier, LREI-RK methods generally achieve higher accuracy for a given step size $h$. Therefore, to reach comparable accuracy, the LREI-AB method typically requires a smaller step size, which diminishes its apparent efficiency advantage.

5.3 q-LL vs. q-LLG dynamics

As noted in the introduction, for a pure initial state $\rho_{0}$, the solutions of the q-LL and q-LLG equations coincide up to a time rescaling $t\mapsto t/(1+\kappa^{2})$. However, no theoretical results are currently available for the case of mixed (non-pure) states. In this section, we compare numerical results for the q-LL and q-LLG equations in a system of nine effective spins, considering both pure and mixed initial states.

We compare the dynamics of magnetization and entanglement (concurrence) under the q-LL and q-LLG models with damping parameter fixed at $\kappa=0.5$. For the pure state case, as we observe from Figure 4, the system is initialized in the antiferromagnetic state $\rho_{0}=|\text{AF}_{2}\rangle\langle\text{AF}_{2}|$, and the two dynamics coincide up to a simple rescaling of time, $t\mapsto t/(1+\kappa^{2})$. In the mixed state case, the initial density matrix is given by a convex combination $\rho_{0}=0.75|\text{AF}_{2}\rangle\langle\text{AF}_{2}|+0.25|\text{GHZ}\rangle\langle\text{GHZ}|$. Figure 5 shows that in this case the q-LL and q-LLG dynamics again exhibit qualitative similarities but no longer coincide under the same time rescaling. This highlights the distinct behavior of the two models beyond the pure-state regime.

Moreover, as expected since the q-LL model is only an approximation of the full q-LLG dynamics, our numerical experiments (not shown here) indicate that the solutions of the two models converge for smaller values of $\kappa$, while they diverge significantly for larger $\kappa$.

5.4 A full-rank state

We close this section with a remark about a well-known high-rank quantum state and how it can be efficiently handled using the LREI method.

The Werner state is a well-known class of mixed quantum states sometimes used as an initial condition in simulations. It is formed as a convex combination of a maximally entangled state and the maximally mixed state, and can be expressed as

[TABLE]

where $\hat{\rho}_{0}$ is a low-rank maximally entangled state, $I$ is the identity matrix, and $I/2^{n}$ represents the maximally mixed state of $n$ qubits. The parameter $p$ controls the balance between the low-rank entangled part and the mixed noise contribution which allows the entanglement to be tuned continuously. This makes the Werner state an ideal test case for studying how the dynamics preserve or degrade quantum correlations. Although $\rho_{W}$ is not low-rank (because the identity term $I/2^{n}$ is full-rank), it differs from a low-rank density matrix $\hat{\rho}_{0}$ only by a renormalization. Therefore, the LREI method can still be applied to solve the q-LL and q-LLG equations with $\rho_{W}$ as an initial condition. In fact, if $\rho(t)$ is the solution of either (1.1) or (1.2) with initial condition $\rho_{W}$, we can write

[TABLE]

where $\hat{\rho}(t)$ is the solution of modified versions of (1.1) and (1.2) where the damping parameter $\kappa$ is replaced by $(1-p)\kappa$, and with the low-rank initial state $\hat{\rho}_{0}$. We then apply the LREI method to approximate $\hat{\rho}(t)$ and finally recover $\rho(t)$ via (5.2). The computation of physical quantities involving expressions of the form $\text{Tr}(A\rho(t))$ for certain sparse observables $A$, simplifies to $p/2^{n}\text{Tr}(A)+(1-p)\text{Tr}(A\hat{\rho}(t))$. The first term is cheap to evaluate due to the sparsity of $A$, while the second term can be efficiently computed using the method described in Subsection 3.6.

6 Summary and future studies

We presented a low-rank numerical algorithm, called low-rank eigenmode integration (LREI) method for solving quantum Landau–Lifshitz (q-LL) and quantum Landau–Lifshitz–Gilbert (q-LLG) equations using Runge-Kutta and Adams-Bashfoth ODE solvers. The LREI method resulted in a significant improvement over the previous full-rank approach. By exploiting the low-rank structure of the density matrix and the sparsity of the operators involved, the proposed method significantly reduces both memory usage and computational cost, making it possible to simulate spin systems of considerably larger size. More precisely, the cost of each time-integration step drops from $\mathcal{O}(N^{3})$ to $\mathcal{O}(r^{2}N)$, where $N=2^{n}$ is the dimension for $n$ spins, and $r\ll N$ is the numerical rank of the solution matrix. Likewise, the memory requirement is reduced from $\mathcal{O}(N^{2})$ to $\mathcal{O}(rN)$, since we never form any full $N\times N$ (or nearly full) matrices during the computation. Among several technical improvements we applied, one key idea was to handle the computation of the action of the invariant subspace of the density matrix associated with its zero eigenvalues. This was accomplished by applying Householder reflectors constructed for the dominant eigenspace, thereby enabling the entire solution process to proceed without ever forming any large matrices.

In addition, similar to the original method, the LREI method preserves physical properties such as spectrum preservation, trace conservation, and positivity of the density matrix in time evolution.

In numerical experiments, we demonstrated that the LREI algorithm achieves the same level of accuracy as the original full-rank method, while scaling up efficiently to systems with up to 24 effective spins (35 periodic spins) on our personal laptop, which is well beyond the reach of the previous technique. For larger spin systems, as low-rank representations and massive (sparse) Hamiltonian components still require storing and manipulating sizable matrices, the bottleneck is memory, not computation time or lack of parallelization. The simulation can be extended beyond 24 spins by adapting and executing the code on computers with higher memory capacity or HPC platforms either through nodes with higher RAM (fat nodes) or by using distributed memory architectures across nodes.

Here we highlight two avenues for future research: (i) developing adaptive time-stepping schemes that adjust the step size using local error estimates to improve efficiency without sacrificing accuracy; and (ii) advancing beyond standard Krylov-type partial eigensolvers by employing state-of-the-art randomized methods to further accelerate computation of the required eigenmodes for large spin systems.

From a quantum physics perspective, our new algorithm, LREI, makes it possible to study dynamical quantum phenomena in many-body spin systems that were previously out of reach. It gives researchers a new tool to uncover subtle quantum effects, gain a deeper understanding of the quantum properties and behavior of complex magnetic systems, and push forward the development of quantum technologies that use spin-based phenomena.