Classical fractional time series from quantum XXZ spin chains

TL;DR

This paper investigates how quantum correlations in the XXZ spin chain can generate classical fractional stochastic processes with varying roughness, revealing the influence of quantum symmetries on the Hurst exponent.

Contribution

It demonstrates the connection between quantum symmetries in the XXZ chain and the fractional properties of classical processes derived from it, using numerical quantum methods.

Findings

Processes with quantum symmetries tend to have H=0 with logarithmic scaling.

Symmetry-breaking processes can produce H≥0.5, but not likely H<0.5.

Numerical methods like MERA and TEBD support these results.

Abstract

Entangled quantum mechanical states in one dimension can be used to represent and simulate classical stochastic processes with nontrivial statistical properties. Long-range quantum correlations translate into fractional processes with their asymptotic Hurst exponents characterizing roughness and persistence. We explore this analogy in the case of the spin-1/2 XXZ chain and investigate properties of four different classical two-state processes that this quantum system can generate. These processes show fractional characteristics with varying Hurst exponents. We argue that the continuous quantum symmetries such as U(1) or SU(2) of the XXZ chain give rise to $H=0$ with logarithmic scaling. Processes generated without these symmetries can produce $H \geq0.5$ but likely not $H < 0.5$ unless the dominant term responsible for $H=0.5$ gets canceled. This does not seem to happen for the XXZ model. We use standard quantum methods, including MERA and TEBD, to numerically substantiate our findings.

Tables

Table 1. Table 1: The details of the analogy between stochastic processes and quantum spin chains.

	Classical stochastic example	Quantum spin chain
Increment process	GBP or fGN	Spin-1/2 quantum chain
Integrated process	fBm of fBP	Spin domain magnetization
Correlations	Power law (in time)	Power law (in space)
Fractal characteristics	Hurst exponent	Correlation exponents
Probability measure	Conditional Bernoulli/normal	Quantum ground state implied
Sample	Process trajectory	Quantum measurement
Sampling algorithm	\pbox20cmCholesky, Circulant, Kernel
or by the definition	\pbox20cmTN (MPS, MERA) sampling
or physical measurement

Table 2. Table 2: The asymptotic forms of the autocovariance function, the exact scaling of the domain magnetizations variance in the Z direction, and the scaling according to long-range numerical calculations in the X direction with the exact formula for Δ = 0 \Delta=0 .

		$\Delta =0$
Z direction	C(n) $\sim n^{\zeta }$	$\zeta =-2$
	$M^2\left(l\right)$	$\log \left(l\right)$
	C(n) $\sim {\left(-1\right)}^n{n}^{\xi }$	$\xi =-2$
	$N^2\left(l\right)$	$l$
X direction	C(n) $\sim n^{\zeta }$	$\zeta =-2.5$
	$M^2\left(l\right)$	l
	C(n) $\sim {\left(-1\right)}^n{n}^{\xi }$	$\xi =-0.5$
	$N^2\left(l\right)$	$l^{1.5}$

Table 3. Table 3: The numerically fitted Hurst exponent for different magnetizations at selected Δ \Delta parameters.

$\Delta$	$M_z$	$N_z$	$M_x$	$N_x$
-0.9	0	0.5	0.5	0.89
-0.5	0	0.5	0.5	0.82
0	0	0.5	0.5	0.74
0.5	0	0.51	0.49	0.67
1	0	0.58	0.05	0.58

Table 4. Table 4: The scaling of the variance for the different magnetizations.

		$\Delta =1$
$M_z^2\propto$	$\log \left(l\right)$
$N_z^2\propto$	$l$
$M_x^2\propto$	$l$	$\log \left(l\right)$
$N_x^2\propto$	$l^{2-\frac{\arccos \left(-\Delta \right)}{\pi }}$

Equations

$$C(n)$$

$$C(n)\sim {\sigma }^{2}H(2H-1)n^{2H-2},\mathrm{for}n\gg 1.$$

$$B_t=\sum _{i=1}^t{X}_i.$$

$$\text{Var}(B_t)={\sigma }^2{t}^{2H}.$$

$$⟨X_t{X}_{t+n}⟩=c{n}^{2h-2},n>=1.$$

$$Var(B_t)\propto \{\begin{array}{ll}t& h\in (0,0.5)\\ t\ln t& h=0.5\\ t^{2h}& h\in (0.5,1)\end{array}$$

$$H=J\sum _{j=1}^L\left[\left(S_j^x{S}_{j+1}^x+S_j^y{S}_{j+1}^y\right)+\Delta S_j^z{S}_{j+1}^z\right],$$

$$C_x(n)$$

$$C_z(n)$$

$$\Delta =-\cos (\pi \eta )$$

$$M_z(l)$$

$$M_x(l)$$

$$\mathrm{Var}(B_t)=⟨B_t^2⟩$$

$$\sum _{j=1}^{k}C(j)$$

$$\mathrm{Var}(B_t)$$

$$\sum _{k=1}^n(-1{)}^k{k}^p$$

$$\sum _{n=1}^l\sum _{k=1}^n(-1{)}^k{k}^p$$

$$C(l)=A{l}^{-\zeta }+(-1{)}^l\tilde{A}{l}^{-\xi },$$

$$⟨M^2(l)⟩=B{l}^{-\zeta +2}+\tilde{B}{l}^{-\xi }+Cl$$

$$⟨N^2(l)⟩=D{l}^{-\zeta }+\tilde{D}{l}^{-\xi +2}+Fl,$$

$$C_z(n)=\frac{1}{2{\pi }^2{n}^2}\left[(-{)}^n-1\right],n\ne 0.$$

$$C_x(n)=\frac{(-1{)}^m}{4}\prod _{k=1}^{[\frac{m}{2}]}\frac{{\Gamma }^2(k)}{\Gamma (k-\frac{1}{2})\Gamma (k+\frac{1}{2})}\prod _{k=1}^{[\frac{m+1}{2}]}\frac{{\Gamma }^2(k)}{\Gamma (k-\frac{1}{2})\Gamma (k+\frac{1}{2})}.$$

$$⟨M^2(l)⟩$$

$$⟨N^2(l)⟩$$

$$M_z(L)|\mathrm{GS}⟩=\sum _{j=1}^L{S}_j^z|\mathrm{GS}⟩=0.$$

$$⟨\mathrm{GS}|M_z^2(L)|\mathrm{GS}⟩=\sum _{i=1}^L\sum _{j=1}^L⟨\mathrm{GS}|S_i^z{S}_j^z|\mathrm{GS}⟩=0.$$

$$0=\sum _{j=1}^L⟨(S_j^z{)}^2⟩+2\sum _{j=1}^L\sum _{n=1}^{L/2-1}⟨S_j^z{S}_{j+n}^z⟩,$$

$$0=⟨(S_j^z{)}^2⟩+2\sum _{n=1}^{\infty }{C}_z(n).$$

$$⟨M_z^2(l)⟩$$

$$⟨M_z^2(l)⟩=-2\left({\Xi }^{(1)}(l)+{\Xi }^{(2)}(l)\right),$$

$${\Xi }^{(1)}(l)$$

Full text

††preprint: APS/123-QED

I Introduction

Fractional processes are stochastic processes with long-range dependence, asymptotic self-similarity, and scale-invariant behavior. They have been applied to many fields, such as economics [1], physics and geophysics [2, 3], computer science [4, 5], financial modeling [6, 7], especially in stochastic volatility models [8, 9, 10]. Fractional processes are represented by scaling exponents, one or many, which determine the asymptotic scaling of their moments, such as their variance. Representing and sampling fractional processes is not straightforward and are an active field of research [11].

Power-law scaling of moments arises from slowly decaying correlations of the process increments. Certain quantum mechanical systems, such as quantum spin chains, show similar correlation functions at their critical point. This creates an interesting analogy between the two systems, which may be used to represent and generate fractional time series. The quantum analogy potentially opens a way to sample fractional processes directly from carefully engineered physical systems with quantum measurement, using quantum computers, or via quantum-inspired methods on classical hardware.

In this study, we explore this analogy through the lens of the Heisenberg XXZ quantum spin chain, which is a well-studied, standard model of one-dimensional strongly interacting discrete quantum degrees of freedom with a rich phase diagram containing an asymptotically self-similar (critical) region with power-law correlations. We demonstrate how the XXZ model can be used to generate fractional series of two states. Having parameter-dependent critical exponents makes the model a potential candidate to generate fractional processes with Hurst exponents both below and above 0.5. Note, however, that the Heisenberg chain is not self-similar in the strict sense, self-similarity only sets in asymptotically. This brings in a dependence on the specific short-range properties of the model. As we will argue, this interferes with and modifies the naively expected scaling behavior for the generated classical process. The analogy between time series and quantum spin chains was already suggested [12], while MPS has been used extensively in recent studies for machine learning tasks such as generative modeling [13, 14, 15, 16] and classification [17, 18, 19].

The paper is structured as follows: in Sections II and III, background information and known results about fractional time series and the XXZ model are introduced, then in Section IV, analytical considerations regarding the scaling of the different processes’ variance are discussed, while Section V, the results of the numerical calculations are presented. Section VI describes sampling techniques for fractional processes and tensor networks. Finally, Section VII summarizes our conclusions.

II Fractional time series

II.1 Fractional Gaussian noise and Brownian motion

The archetype of fractional processes is the fractional Brownian motion (fBM).[20] With our forthcoming line of reasoning in mind, we define it here in discrete time through its increments. We start with fractional Gaussian noise (fGN) $X_{t}$, which is a stationary Gaussian process with auto-covariance

[TABLE]

where $n\geq 0$ is an integer time lag, $\sigma^{2}=\langle X_{t}X_{t}\rangle$ is the constant variance of fGN, and $0<H<1$ is its Hurst exponent. For large time distances, $C(n)$ behaves like a power law:

[TABLE]

Fractional Brownian motion (fBM) is defined as the cumulative process of fGN increments

[TABLE]

As a result of the covariance structure in Eq. (1), fBM has a time-dependent variance that behaves exactly as a power law with exponent $2H$ (see Section IV.1 for a derivation of how the short-term properties of $C(n)$ play a role in this):

[TABLE]

Fractional Brownian motion is exactly self-similar (not just asymptotically) in terms of its variance.

II.2 Generalized Bernoulli process and fractional binomial process

Discrete state processes are better suited for the analogy with finite-level quantum systems. The generalized Bernoulli process (gBP) proposed by Lee [21] is a discretized version of the fractional Gaussian noise. The process built up of these increments is called the fractional binomial process (fBP).

The generalized Bernoulli process is a two-state stochastic process with $Prob(1)=p$, $Prob(0)=1-p$, defined through its auto-covariance function:

[TABLE]

with $c$ and $h$ constants. This mimics the asymptotic form of the fGN covariance in Eq. (1), but not its specific short-term features.

The variance of the fBP for large distances turns out to be [21]:

[TABLE]

Note that, contrary to the fBM, the fraction Binomial process does not produce rough scaling for $h<0.5$. The variance cannot grow slower than linear, its Hurst exponent is $H=1/2$. Sublinear scaling is suppressed by an $H=1/2$ scaling operator, which is not canceled as in the case of the fBM without an appropriately engineered short-term structure.

III Heisenberg XXZ 1/2-spin chain

The Heisenberg XXZ 1/2-spin chain with anisotropy parameter $\Delta$ is defined by the Hamiltonian:

[TABLE]

where $S^{a}_{j}$ is the spin-$\frac{1}{2}$ projection operator in direction $a$ acting on site $j$, $J$ is the coupling, chosen to be $1$ and $\Delta$ is the anisotropy parameter.

The system is in a gapless phase if $|\Delta|<1$ and experiences long-range correlations in its paramagnetic ground state. Asymptotic results are known for the correlation of spin-spin operators from field theoretical calculations. Exact results are hard to calculate but are known for small distances [22, 23]. For $\Delta=0$ the $\langle S^{z}_{i}S^{z}_{i+l}\rangle$ auto-covariance takes a simple form and the $\langle S^{x}_{i}S^{x}_{i+l}\rangle$ is also calculable in closed form [24].

The leading term in the critical region is an algebraically decaying correlation for all directions [25], similar to the auto-covariance of the GBP (Eq. (5)). The dominating terms are:

[TABLE]

where

[TABLE]

connects scaling exponents to spin chain anisotropy.

The asymptotic behavior is a sum of a monotonically decaying and an oscillating term. Hence, it is customary to define and study two spin aggregation functions, the ”uniform” and the ”staggered” domain magnetizations over length $l$. With longitudinal and transverse directions wrt the direction of anisotropy, we have four interesting quantities to study:

[TABLE]

Spin chains provide a natural analogy with fractional stochastic processes based on their stationary, zero-mean increments and asymptotic scaling of the autocovariance of the increments (see Table 1). In Ref [12] we showed how to leverage this analogy to generate classical time series from quantum states such as the ground state of the XXZ model. The fractal properties of the time series generated will reflect the statistical properties of the magnetic domains. Hence, in the next section, we will investigate these in detail.

IV Scaling of the variance of magnetic domains

In general, we can calculate the variance of any cumulative process using the increments’ covariance function. With $X_{t}$ being the increment process, we define the cumulative process as $B_{t}=\sum_{i=1}^{t}X_{i}$. When the increment process has zero mean $\langle B_{t}\rangle=0$, the variance of the cumulative process can be expressed as:

[TABLE]

where we have used that increments are stationary with $\left\langle X_{i}^{2}\right\rangle=\sigma^{2}$. The first, linear term is always present irrespectively of the autocovariance structure and in itself alone would give rise to $H=0.5$ scaling.

Asymptotic scaling of the variance can be derived from the scaling of the covariance function. To see this, look at the double sum as a double integral, which is a valid approximation in the large domain length limit. If there is a nonzero constant term coming from the first integral, it will be integrated again, leading to a term proportional to $l$. This $\mathcal{O}\left(l\right)$ term is irrelevant when $C(l)\sim l^{-\alpha},\quad\alpha<1$, i.e., decays slower than $1/l$, because then the dominant contribution from the double integral is $\mathcal{O}\left(l^{2-\alpha}\right)$ which increases faster than linear. This superlinear scaling leads to a Hurst exponent $H>0.5$ as expected.

On the other hand, when $C(l)\sim l^{-\alpha}$ with $\quad\alpha<1$, i.e., covariance decays faster than $1/l$, the double integral scales as $\mathcal{O}\left(l^{2-\alpha}\right)+\mathcal{O}\left(l\right)$ in which $\mathcal{O}\left(l\right)$ dominates unless the constant term of the first integral vanishes for any reason. Without this the Hurst exponent is $H=0.5$.

IV.1 Fractional Brownian motion

It is illuminating to evaluate the double sum term exactly for fractional Brownian motion and see how it interacts with the linear term. With the covariance function in Eq. (1), the inner sum of the second term is a telescopic series:

[TABLE]

Substituting back we get a telescopic series again in the outer sum. The whole expression becomes

[TABLE]

Notice that the linear term gets canceled exactly. This ”miraculous” cancellation is the consequence of the particular small-distance auto-covariance structure of the fractional Gaussian noise process. What remains for $\mathrm{Var}(B_{t})$ is a single power-law term whose exponent is $2H$ and thus behaves sub- or super-linearly below or above $H=0.5$.

IV.2 Generalized Bernoulli process

In this case the linear term does not cancel as was already found by J. Lee.[21]. The variance becomes dominated by this linear scaling term for $h<0.5$, and this suppresses sub-linear scaling. For $h>0.5$, super-linear scaling arises as expected, and the Hurst exponent is $H=h$.

IV.3 Heisenberg spin chain

For a Gaussian process, like the fGn or the fBm, the expectation value and the (auto)covariance function are sufficient to define the process. This is not true for non-Gaussian processes, like the fBP or gBP. Many processes with the same expectation value and covariance structure can be constructed.

A major property of fractional processes is the power-law scaling of the variance with the Hurst exponent. In this study, we will only concentrate on the variance (second central moment) of a given process when analyzing its fractional and self-similar nature, arguing that this can be easily obtained for real-world time series.

Without knowing the exact short-distance form of the auto-covariance, it is impossible to calculate the variance of the magnetization. Even the asymptotic behavior can depend on the unknown short-distance behavior.

It is important to see the contribution of the staggered term of the covariance to the uniform magnetization and vice-versa for any real $p$ (for derivation see Appendix A):

[TABLE]

In summary, if the autocovariance of the increments has the long-distance asymptotic of:

[TABLE]

then the long-distance asymptotics of the variance of the uniform cumulative process are:

[TABLE]

and

[TABLE]

if $\zeta,\xi\notin\{1,2\}$. Only if the constants $C$ or $F$ are [math], can we expect a behavior that is similar to the fBm with $H<0.5$. These calculations can be applied to the Heisenberg spin chain to predict the asymptotic behavior of the domain magnetization’s variance. The resulting exponents can be seen in Fig. 4 with a dashed line if $C$ or $F$ is [math] and with a solid line if not [math].

There are two cases where we can make stronger statements: a) when $\Delta=0$ and b) when there is a conserved quantum number behind the generating quantum system.

IV.4 Case of $\Delta=0$

At the XX point, $\Delta=0$, the exact form of the covariance function is known exactly [24]:

[TABLE]

and

[TABLE]

Using the above formulae, the variance of the uniform and staggered domain magnetization in the $z$ direction can be calculated analytically. The calculation details are presented in Appendix B, we only state the end results here. For the uniform magnetization we find

[TABLE]

and the staggered magnetization reads

[TABLE]

We refer to the Appendix to see how the linear term vanishes exactly in the uniform case.

For the $x$ direction, we cannot derive results analytically, but we can use the exact form of the covariance function and do the summation numerically for long distances. Results are shown in Fig. 1. Behavior of the uniform magnetization is consistent asymptotically with $H=1/2$, but we need long chains of hundreds of spins to see this with precision. This means that the linear term of the double summation does not vanish and dominates the behavior. On shorter distances, the empirically fitted exponent can appear misleadingly smaller.

The staggered magnetization can be fitted well with an asymptotic power law with $H\sim 0.75$. This is consistent with the exponent $\xi=1/2$ of the covariance function through Eqs. (9) and (18). Since this is superlinear, the linear term does not count.

IV.5 Case of a quantum symmetry

The $U(1)$ symmetry of the model, i.e., the invariance under the simultaneous rotation of spins around the $z$-axis, leads to the conservation of the $z$-component of the total spin, and this has important bearings on the scaling of the uniform magnetization $M_{z}$. Suppose, for a moment, a chain of length $L$ with periodic boundary conditions. Assume that $L$ is even, which makes the derivation simpler but has no relevance in the $L\rightarrow\infty$ limit.

For $|\Delta|\leq 1$, the gapless ground state of the model is an eigenstate of $M_{z}(L)$ with zero eigenvalue,

[TABLE]

Consequently, the variance of the magnetization is also zero,

[TABLE]

We now change variables in the sums,

[TABLE]

then use translational invariance and take the $L\rightarrow\infty$ limit to get the simple sum rule,

[TABLE]

Note that this involves the infinite sum of the covariance function and connects it to the single-site variance.

The variance of the magnetization of a finite segment reads as

[TABLE]

Using the sum rule Eq. (26) we get

[TABLE]

or

[TABLE]

with

[TABLE]

Thus, we have rewritten the double sum into simple sums. It is safe to insert the asymptotic forms. This leads to an unknown constant error in $\Xi^{(1)}(l)$ which depends on the short-distance behavior of $C_{z}(l)$, but this will become irrelevant. The second term $\Xi^{(2)}(l)$ looks dangerous due to the $l$ prefactor, but here the sum starts only at $l$, and we can safely use the asymptotic form if $l\gg 1$. By Eq. (LABEL:eq:spin_cov) and (52) the two expressions read:

[TABLE]

As $\eta<1$, the leading term is the logarithm, i.e., for $l\gg 1$ we expect:

[TABLE]

The $\Delta=1$ point is special in two senses: (I) not only $M_{z}(L)$ but also $M_{x}(L)$ is conserved and zero due to the full $SU(2)$ symmetry, leading to identical behavior in the two spin directions, (II) in this limit ($\eta\rightarrow 1$) the correlator $C_{z}(n)$ has a slightly different form [26],

[TABLE]

This, however, does not affect the leading term that is still $\propto\ln l$.

V Numerical results

We have seen above that the asymptotic scaling of (uniform and staggered) domain magnetization depends much on the short-term structure of the correlation functions and the symmetry of the quantum state. Consequently, any classical stochastic process defined from these quantum states inherits this dependency. The fractional properties that can be achieved via this quantum-classical mapping may be restricted.

Out of the 4 processes we have defined we expect $M_{z}$ to scale logarithmically slowly and thus effectively represent Hurst $H=0$ due to the $U(1)$ continuous symmetry. There is no such symmetry for the other three processes except at $\Delta=1$ for $M_{x}$. In the general cases $M_{x}$, $N_{z}$, $N_{x}$ can show sub- or super linear scaling with $H<1/2$ or $H>1/2$ depending on $\Delta$. We see no difficulty with super-linear scaling, but sub-linear scaling requires a constant to become exactly zero in the correlation (auto-covariance) structure. Whether this condition holds cannot be determined analytically in the XXZ model. Hence we turn to numerical methods.

We calculated the ground state of the Heisenberg XXZ chains with MERA and MPS ansatzes at different anisotropy parameters. We could then calculate two-point correlations and magnetic domains of different lengths.

The MPS was calculated using the infinite time-evolving block decimation (iTEBD) algorithm [27] with a bond dimension of 1,024 and exploiting one $U(1)$ symmetry, using the opeNAMPS package [28, 29]. The scale-invariant MERA was calculated using a gradient-based method [30] with a maximal bond dimension of 20 and two transitional layers. The number of parameters in these settings, accounting for the independent components of the tensors, is ($2\chi^{2}$) 2,097,152 in the case of the MPS and 334,416 for the MERA. This difference accounts for the higher accuracy of the MPS representation but is necessary because of the increased computational demands of working with the MERA ansatz.

V.1 Autocorrelation of increments

In Fig. 2 we show a comparison of the short-term behavior of spin chain correlators and asymptotically equivaent fraction Brownian motion and fractional Binomial processes. Deviations are small, except for very small distances, but these can be responsible for significant qualitative alterations in domain scaling behavior.

V.2 Scaling of the variance of the domain magnetization

We also calculated the variance of the cumulative process (domain magnetization) using Eq. (11). Looking for asymptotic scaling in the form of $l^{2H}$, we can determine the Hurst exponents $H$ numerically.

We calculated the variance of domain magnetization up to a distance of 1,000 sites. Going further would accumulate numerical errors as the tensor network calculations loose precision. We will describe what can be seen up to these distances based on Fig. 3 and the fitted parameters in Table 3. The behaviors are summarized in Table 4 and Fig. 4.

Process $M_{z}$

For $M_{z}$, the variance of the uniform magnetization scales logarithmically for all $\Delta$. This is visible by the straight lines on the lin-log plot. This is the only case where we see no power-law scaling. As discussed this is a consequence of the $U(1)$ symmetry.

Process $N_{z}$

If we consider the staggered magnetization in the $z$ direction, we see a slowly increasing scaling exponent from $0.5$ to $0.58$. According to our previous calculations, this is an artifact on short distances, and all curves will converge asymptotically to the exponent $0.5$. This was shown for $\Delta=0$ analytically earlier.

Process $M_{x}$

The uniform magnetization in the $x$ direction is an interesting case. Except for $\Delta=1$, we see close to linear scaling. However, for short distances, the scaling exponent tends to [math] as $\Delta$ is increased, which would correspond to rough behavior. Only with increasing distance does the linear term of equation (11) take over and lead to $H=1/2$. This result is in line with Ref. [25], who found that for $|\Delta|<1$, scaling is linear. We obtained a similar result for $\Delta=0$ in our exact calculation earlier. The $\Delta=1$ point is different: here we observe the crossover to logarithmic scaling due to the emerging SU(3) symmetry.

Process $N_{x}$

The staggered magnetization in the $x$ direction is the most straightforward to interpret, as we see fractional behavior with scaling exponents between $0.5$ and $1$, which is in agreement with the expectations. The reason why the exponent is not decreased below $0.58$ is likely the same artifact mentioned above: the interplay between the different subleading terms results in a behavior that looks like a higher-order scaling at short distances. The theoretical connection between the anisotropy parameter and the Hurst exponent is $H=1-\frac{\arccos(-\Delta)}{2\pi}$.

VI Series generation methods

Our quantum chain to classical time series mapping allows us to generate classical series from the ground state (or any other quantum state) of the quantum spin chain. If the quantum system has algebraically decaying long-range correlations and hence non-trivial full counting statistics with self-similarity for its magnetic domains we believe this structure will get transferred to the time series. It is important, however, how series generation can work in practice and how it scales with the length of the time series we want to sample. We discussed the general sampling methodology in an earlier publication [12], but have not focused on the specific complications of fractional series where sampling is known to be a difficult problem classically.

VI.1 Classical

For the simulation of fractional Brownian motion, several algorithms have been proposed [11, 31, 32]. They all exploit the Gaussian property of the stochastic process and are, therefore, inapplicable to the fractional binomial process. The Davies-Harte is the optimal known exact method [32] that scales as $\mathcal{O}(N\log N)$ with the predefined number of time points ($N$). If the sample size is not known in advance, one can use the Hosking method, which has $\mathcal{O}(N^{2})$ scaling.

VI.2 Quantum-inspired

The exponentially growing state space makes approximating methods necessary to get the ground state of a quantum spin chain longer than a couple of dozen spins.

One well-established approach uses tensor network states as an ansatz for the ground state of spin chains [33]. These try to approximate the high-dimensional quantum state by breaking it up into many interconnected smaller tensors with a specified bond dimension. The simplest type of tensor network states are matrix product states (MPS) [34], which can be calculated by powerful density-matrix renormalization group (DMRG) algorithms [35] or time-evolving block decimation (TEBD) [36], but can only approximate the long-range autocorrelation. Yet, another tensor network ansatz, the multiscale entanglement renormalization ansatz (MERA), was proposed to naturally represent algebraically decaying correlations [37, 38]. Both MPS and MERA can be optimized with automatic differentiation-based algorithms [39].

For practical stochastic time-series generation, we would need to sample these ground states. Fortunately, methods already exist for sampling from finite tensor network states [40, 41]. We extended these methods for infinite MPS and MERA states, see below.

Sampling from an MPS scales linearly with the number of sites. However, it scales cubically with the bond dimension and the required bond dimension to maintain precision at long distances (or the correlation length) scales with $N^{\kappa}$, where $\kappa=0.73$ for the Heisenberg model [42]. Sampling from a MERA has a $N^{3\log_{3}{\chi}}$ scaling, at worst, but is highly dependent on the contraction order of the tensors. The number of MERA tensors grows only linearly. In a previous work, we described in detail how an infinite MPS can be sampled [12]. Here, we extended these methods to be able to sample from infinite MERA states.

VI.3 Sampling an infinite MPS

The stochastic sampling of tensor networks was originally worked out by Ferris and Vidal [40]. They worked with finite chains and assumed a canonical form such as a unitary tensor network, e.g., uMPS. In case of infinite chains we can relax this criterion but maintain the same linear scaling of the computational cost. We need to calculate the dominant left and right eigenvectors of the MPS transfer matrix ($l$, $r$) and use them as the left and right MPS vectors to close the (very long) chain (see Fig. 5) making it effectively infinite. The sampling is done site-by-site. An example for the third site is shown in Fig. 6. The figure shows how the density matrix for the third site is calculated after the first two sites are already fixed in the $\uparrow\downarrow$ states. The density matrix for the third site can be used to sample its spin state with probabilities of up and down:

[TABLE]

For a detailed description, see our previous work [12].

VI.4 Sampling an infinite MERA

Sampling works analogously in other tensor networks. Suppose we have a scale-invariant (infinite) 1D ternary MERA with $T$ transitional layers. The isometries are ($w_{0},w_{1}\dots w_{T},w$), the disentanglers ($u_{0},u_{1}\dots u_{T},u$). We can calculate the average two-site density matrix ($\rho_{2}$) of the scale-invariant layers by finding the dominant eigenoperator of the average descending superoperator.

VI.4.1 The average three-site MERA density matrix

In a ternary MERA, every $N$-site operator will eventually be mapped to a two- or three-site operator after enough layers. Three-site operators are, however, mapped to two- or three-site operators in the next layer, depending on their position, making infinitely many layers necessary to map them to a two-site operator. This is done by the three three-site ascending superoperators $\mathcal{A}^{3\rightarrow 2}_{L}$, $\mathcal{A}^{3\rightarrow 3}_{C}$, $\mathcal{A}^{3\rightarrow 2}_{R}$, shown in Fig. 7. They all occur on $\frac{1}{3}$th of the positions in every layer.

We can derive an analytical expression for the average three-site density matrix to allow sampling and the calculation of expectation values for all $N$-site operators.

Suppose we calculate the expectation value of a three-site operator $O_{3}$ - according to the middle-ground contraction scheme - by applying the ascending superoperator on $O_{3}$ until the three-site operator is mapped to a two-site operator and applying the descending superoperator on the average two-site density matrix $\bar{\rho}_{2}$ to descend the density matrix to the same layer:

[TABLE]

where we exploited the fact that $\bar{\rho}_{2}$ is an eigenoperator of $\mathcal{D}$ with eigenvalue $1$. The condition for closed-form substitution of the infinite sum in the last line is that the dominant eigenvalue of $\mathcal{A}^{3\rightarrow 3}_{C}$ is smaller in absolute value than three.

If we define the average three-site ascending superoperator $\mathcal{A}^{3\rightarrow 2}$:

[TABLE]

and consider that:

[TABLE]

is true for arbitrary $O_{3}$. We see that the average three-site density matrix $\bar{\rho}_{3}$ is:

[TABLE]

Calculating the inverse in (36) can be an obstacle due to the $\chi^{12}$ size of $\mathcal{A}^{3\rightarrow 3}_{C}$. In this case, one can fall back to the infinite sum form of (35) and calculate it approximately.

VI.4.2 Sampling

If we want to sample $N$ sites from MERA, we build a MERA network with $N$ physical sites and with at least $T$ layers, but as many as we need to coarse-grain them to two or three sites and then connect them to the average two- or three-site density matrix.

The resulting network is an $N$-site density matrix with $2N$ free indices. Of course, this would be untreatable in the matrix form. For the sequential sampling from left to right, we first trace out the rightmost $2N-2$ indices and take a sample according to the probabilities of the resulting one-site density matrix. In the next step, we fix the first site, trace out $2N-4$ indices, and take a sample according to the resulting one-site density matrix. We repeat the steps until no free indices are left, and then we have a length $N$ sample from the MERA network.

For an exact result, one would need to average all possible layouts, in other words, over every different $N$-site density operator. As there are three possible layouts per layer, as shown in Fig. 8, this multiplies the computational cost by three with every renormalization layer. However, it was previously shown [43] that fluctuations due to the violation of translational invariance of the MERA network are of order $10^{-4}$ and diminish with growing bond dimension. For the sake of simplicity and demonstration, we used only a single layout, where the first free index is the middle leg. Such a MERA network for sampling $N=10$ sites is shown in Fig. 9.

Calculating expectation values of two-site operators like two-point correlation can be done in a very similar way, e.g., with the network in Fig. 9. We just put the two operators on the particular sites and trace the rest.

Finding the optimal contraction order is a hard task. However, the computational cost can be estimated with a systematic contraction order. Such a scheme is shown in Fig. 10. Suppose that the $n-1$ sites are already fixed and we want to sample the $n$th site. We start at the bottom layer in the middle of the network (Fig. 10 (a)), where no free indices are present. First, we contract the isometries with the disentanglers in the bottom layer. This is a $\chi^{4}$ operation. In Fig. 10 (b), we see the resulting network, with the additional generalization that the tensors on the bottom are connected with $k$ parallel bonds. This is because, as we will see, moving up layer by layer, only the effective bond dimension of this line will multiply. To see this, we contract all the tensors between the dashed lines to end up with the network in Fig. 10(c). This can be done with cost scaling as $\chi^{3k+4}$ in the order shown in the figure. After doing it, we end up with one less layer, and the new bottom layer will have a form similar to the previous one, only the tensors at the bottom will have one more parallel bond between them.

If we are contracting the part with the free index, that will only multiply the cost by the physical dimension. For the whole network with $\log_{3}(n)$ layers and $\frac{n}{3^{k}}$ tensors in the $k$th layer, the computational cost scales with:

[TABLE]

To obtain an $N$-site sample, one would need to do the above contraction for $n=1,2\dots N$. This leads to a computational cost that is upper bounded by $N^{3\log_{3}{\chi}+1}$.

VII Conclusions

Power-law scaling is a characteristic feature of both one-dimensional quantum systems at their critical point and classical fractional time series. This formal similarity can be made rigorous, and a mapping can be established that defines classical fractional series from quantum systems and creates a means to generate/sample these series in a sequential way. The concepts of this mapping were described in an earlier publication [12]. In this paper we focused on the mapping between the ground state of the spin-1/2 Heisenberg XXZ chain and four different fractional stochastic processes of two states built upon this quantum model.

The Heisenberg model has an extended critical region in which critical exponents depend on a parameter called anisotropy. This brings the hope that these two-state fractional processes can be created with different Hurst exponents. A priori, simply looking at the asymptotic decay of quantum correlations, we expected to see Hurst exponents above and below $1/2$. Our analysis substantiated this expectation for $H\geq 1/2$. It is easy to generate fractional processes with positive persistence from the spin chain. On the other hand, we identified a difficulty for $H<1/2$, i.e., anti-persistent (rough) behavior. We realized that our mapping directly introduces a scaling operator with linear scaling, which in itself leads to $H=1/2$. When this operator is present the anti-persistent behavior is suppressed and non-fractional scaling prevails. There is only one case when this operator vanishes: continuous quantum symmetries dictate a special sum rule of the quantum correlation function, and as a result, the linear operator is canceled. A $U(1)$ symmetry is present in one of the processes we defined for any value of the chain anisotropy, and an $SU(2)$ symmetry in a special point of the parameter space for two of the processes. In these cases we do see $H<1/2$. However, in the XXZ Heisenberg chain this can only lead to $H=0$ exactly with logarithmic scaling. It is an open question if there is any other one-dimensional quantum system where a quantum symmetry kills the dangerous linear operator, but it still possesses a scaling operator in its spectrum with $H>0$. This question is left for future work.

In the second half of the paper, we discussed how fractional series can be generated from a tensor network representation of the quantum state. The simplest generation method works with MPS (matrix product) states [44]. This has linear scaling in the chain length. Note, however, that the MPS is only an approximation of the quantum state, and its representation accuracy becomes worse close to the critical point. The MPS requires more and more parameters to keep accuracy as the critical point is approached, and this ruins linear scaling at the critical point. Other tensor networks have been proposed in the literature for better approximating quantum critical systems, such as MERA.[37, 38] Here we extended sampling from MERA states for infinite chains and discussed how MERA can be used to generate a classical fractional process in this context. Working with MERA, though, has its own complexity, and it is unclear if the complexity of the process is counterbalanced by the smaller number of parameters MERA requires for the same precision. Again, this is left for future studies.

VIII Acknowledgments

This work was supported by the Ministry of Innovation and Technology and the National Research, Development and Innovation Office of Hungary (NKFIH) within the National Quantum Technology Program (Grant No. 2022-2.1.1- NL-2022-00004) and Grant Nos. K134983, TKP2021-NVA-04. Ö.L. acknowledges support of the Hans Fischer Senior Fellowship programme funded by the Technical University of Munich - Institute for Advanced Study, and of the Center for Scalable and Predictive methods for Excitation and Correlated phenomena (SPEC), funded as part of the Computational Chemical Sciences Program FWP 70942 by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences at Pacific Northwest National Laboratory. M.A.W. has also been supported by the Janos Bolyai Research Scholarship of the Hungarian Academy of Sciences. Z.U. acknowledges the support of the Doctoral Student Scholarship Program of the Co-operative Doctoral Program of the Ministry of Innovation and Technology, financed by the National Research, Development and Innovation Fund of Hungary.

.1 Auxiliary functions

We list first some special functions that will be used in the following derivations. The Riemann zeta function $\zeta(a)$ and the Dirichlet eta function $\eta(a)$ are defined as

[TABLE]

We will need the partial sums

[TABLE]

where $H_{n}^{(a)}$ is called the generalized harmonic number. Easy to see that the two are connected; for $n$ even, we get

[TABLE]

which can be compared to the infinite sum result

[TABLE]

For positive integer powers $m\in\mathcal{N}$, these sums can also be expressed using the polygamma functions [45]:

[TABLE]

where $\gamma=0.5772...$ is the Euler-Mascheroni constant, and $\psi^{(m)}$ is the $m$-th derivative of the digamma function (aka polygamma). In particular for $m=2$ and $n$ even:

[TABLE]

Later, we will need for $n$ even

[TABLE]

where for the last equality see https://en.wikipedia.org/wiki/Trigamma_function.

For large $n$ the Euler-Maclaurin [46] sum formula can be used:

[TABLE]

and for the Dirichlet sum for $n$ even, we find

[TABLE]

The first we found empirically also valid for $a>0$ (using the analytic continuation of the zeta function for $a<1$). For $n$ odd we have a slightly different expansion.

Appendix A Alternatig power sum

A.1 Single sum

We can write the alternating partial sum as the difference between two non-alternating partial sums:

[TABLE]

We will use the series expansion of the generalized harmonic numbers $H_{n}^{p}$ (obtained with the help of the Euler-McLaurin formula) to obtain the asymptotic behavior of the sum. It is valid for any real $p\neq-1$.

[TABLE]

where $\zeta$ denotes the Riemann zeta function.

For $p=-1$ the behavior is similar:

[TABLE]

where $\gamma$ is the Euler-Mascheroni constant.

For even $l$:

[TABLE]

For odd $l$:

[TABLE]

Now if combine (50) and (A.1) we end up with the general formula:

[TABLE]

[TABLE]

A.2 Double sum

[TABLE]

[TABLE]

for every $p\neq 2$.

Appendix B Exact calculations of magnetization variance ($\Delta=0$)

At the XX point, $\Delta=0$, the exact form of the autocovariance function is known for the Heisenberg spin chain [24], see Eqs. (19) and (20). This makes it possible to calculate the variance of the magnetization without any uncertainty.

B.1 Uniform magnetization in $z$ direction ($M_{z}$)

Introduce

[TABLE]

and the difference

[TABLE]

Substituting in the autocovariance function:

[TABLE]

Using $\zeta(2)=\pi^{2}/6$, the constant term cancels, and we get

[TABLE]

Or alternatively, with the exact form above

[TABLE]

With this, the variance of the uniform domain magnetization in the $z$ direction reads

[TABLE]

The sums on the right-hand side have closed-form results as

[TABLE]

which gives

[TABLE]

As we see, it scales sublinearly. Note that

[TABLE]

B.2 Staggered magnetization in $z$ direction ($N_{z}$)

Introduce

[TABLE]

and the difference

[TABLE]

Substituting in the autocovariance function:

[TABLE]

Using $\zeta(2)=\pi^{2}/6$, now the constant term does not cancel, and we get

[TABLE]

Or alternatively, with the exact form above

[TABLE]

With this, the variance of the staggered domain magnetization in the $z$ direction reads

[TABLE]

Again, the sums on the right-hand side have the closed-form results

[TABLE]

which gives

[TABLE]

As we see, it scales linearly with a log correction as the subleading term. So the Hurst exponent is $H=1/2$.