* (54) P. P. Fiziev, _Teukolsky-Starobinsky Identities: A Novel Derivation and Generalizations_ , _Phys. Rev. D_ 80 (2009) 124001, [0906.5108].
  * (55) V. Cardoso, A. S. Miranda, E. Berti, H. Witek and V. T. Zanchin, _Geodesic stability, Lyapunov exponents and quasinormal modes_ , _Phys. Rev. D_ 79 (2009) 064016, [0812.1806].
  * (56) H. Yang, D. A. Nichols, F. Zhang, A. Zimmerman, Z. Zhang and Y. Chen, _Quasinormal-mode spectrum of Kerr black holes and its geometric interpretation_ , _Phys. Rev. D_ 86 (2012) 104006, [1207.4253].
  * (57) A. Gruzinov and M. Kleban, _Causality Constrains Higher Curvature Corrections to Gravity_ , _Class. Quant. Grav._ 24 (2007) 3521–3524, [hep-th/0612015].
  * (58) G. T. Horowitz, M. Kolanowski and J. E. Santos, _Almost all extremal black holes in AdS are singular_ , _JHEP_ 01 (2023) 162, [2210.02473].
  * (59) M. Casals, S. E. Gralla and P. Zimmerman, _Horizon Instability of Extremal Kerr Black Holes: Nonaxisymmetric Modes and Enhanced Growth Rate_ , _Phys. Rev. D_ 94 (2016) 064003, [1606.08505].
  * (60) G. Goon and R. Penco, _Universal Relation between Corrections to Entropy and Extremality_ , _Phys. Rev. Lett._ 124 (2020) 101103, [1909.05254].
  * (61) P. A. Cano, S. Chimento, R. Linares, T. Ortín and P. F. Ramírez, _$\alpha^{\prime}$ corrections of Reissner-Nordström black holes_, _JHEP_ 02 (2020) 031, [1910.14324].
  * (62) D. Li, A. Hussain, P. Wagle, Y. Chen, N. Yunes and A. Zimmerman, _Isospectrality breaking in the Teukolsky formalism_ , _Phys. Rev. D_ 109 (2024) 104026, [2310.06033].


# A restricted eigenvalue condition for unit-root non-stationary data

Etienne Wijler

(Vrije Universiteit Amsterdam
Department of Econometrics and Data Science
)

###### Abstract

In this paper, we develop a restricted eigenvalue condition for unit-root non-
stationary data and derive its validity under the assumption of independent
Gaussian innovations that may be contemporaneously correlated. The method of
proof relies on matrix concentration inequalities and offers sufficient
flexibility to enable extensions of our results to alternative time series
settings. As an application of this result, we show the consistency of the
lasso estimator on ultra high-dimensional cointegrated data in which the
number of integrated regressors may grow exponentially in relation to the
sample size.

Keywords: REC, random matrix theory, lasso, cointegration.

JEL-Codes: C32, C55

##  1 Introduction

Modern time series applications frequently involve the analysis of high-
dimensional datasets in which the number of time series $N$ is large in
relation to, and may exceed, the number of temporal observations $T$. A rather
successful approach to modelling such datasets has been to impose some form of
sparsity on the underlying data generating process (DGP) and to incorporate
$\ell_{1}$-based regularization such as the lasso to reduce model complexity
and/or recover the true sparsity pattern (e.g Basu and Michailidis,, 2015;
Medeiros and Mendes,, 2016; Masini et al.,, 2019). However, these approaches
heavily rely on stationarity of the considered time series, thereby either
excluding the presence of stochastic trends or requiring correct
identification and corrections of the order of integration. As noted by
Smeekes and Wijler, (2020), a key difficulty that hinders extensions of the
theory for lasso-type estimators to high-dimensional non-stationary settings
is the absence of a verifiable restricted eigenvalue condition (REC).
Accordingly, we aim to fill this gap in the literature by proposing an
appropriately scaled REC for unit root non-stationary time series and
verifying its validity directly on the gram matrix with the use of matrix
concentration inequalities. This result is then used to derive novel finite
sample error-bounds for the lasso and prove its asymptotic consistency on
(ultra) high-dimensional (co)integrated datasets.

The behaviour of $\ell_{1}$-based regularization on stationary time series
data is well-researched. Basu and Michailidis, (2015) and Kock and Callot,
(2015) consider lasso based estimation of high-dimensional VAR processes
driven by Gaussian innovations and derive tight error bounds for the
prediction and estimation error. These bounds show that the estimation error
deteriorates in the dimension at a logarithmic rate, thereby validating the
lasso in (ultra) high-dimensional time series settings. Medeiros and Mendes,
(2016) consider the adaptive lasso and derive both estimation and selection
consistency on weakly dependent and possibly heteroskedastic time series data
without imposing Gaussianity on the predictors or innovations. By replacing
the assumption of Gaussianity by a set of moment restrictions, they provide a
considerable generalization to previous results, albeit at the cost of error
bounds that deteriorate at faster rates in the dimension. Masini et al.,
(2019) generalize the aforementioned results on sparse VARs estimated with the
lasso by deriving error bounds under the assumption of weak sparsity, while
only assuming the innovation vector to be geometrically strong
($\alpha$-)mixing and a martingale difference process. Finally, Wong et al.,
(2020) provide further generalizations by avoiding the assumption of a
specific parametric form of the DGP and relying only on (strict) stationarity
and geometrically decaying $\beta$-mixing coefficients to establish error
bounds for the lasso for subweibull random vectors.

Owing to the lack of suitable restricted eigenvalue conditions, the literature
on the lasso in the non-stationary setting is substantially more scarce, with
current results being limited to fixed or moderate dimensional settings only.
For example, in the fixed dimensional setting, Liao and Phillips, (2015)
develop a penalized estimator for a vector error-correction model (VECM) and
derive its consistency and pointwise limit distribution. Lee et al., (2021)
consider a two-stage adaptively weighted lasso procedure applied to
(co)integrated data for which the oracle property is shown to hold. Koo et
al., (2020) consider the lasso applied to predictive regressions with
covariates that are integrated of order one and derive error bounds under the
assumption that the gram matrix satisfies an REC. The authors correctly note
that the resulting lower bound that is assumed on the restricted eigenvalue
may decrease to zero as the number of integrated time series grows, thereby
inflating the estimation error at an unknown rate. Accordingly, consistency is
again only derived for the case of a fixed number of integrated time series.
Moving on to the moderate dimensional setting, Liang and Schienle, (2019)
develop a sparse VECM estimator that is argued to be consistent in a setting
where the number of stochastic trends in the data increases at a rate slightly
slower than $T^{1/4}$. Smeekes and Wijler, (2020) consider lasso-type
estimation of a single-equation error-correction model and derive consistency
in an asymptotic framework in which the number of integrated time series is
allowed to grow at rate $T^{1/4}$. In addition, the authors show that the
minimum eigenvalue of the usual non-stationary gram matrix converges to zero
as the dimension increases, but can be bounded away from zero with probability
converging to one when scaled up by its dimension if the latter grows no
faster than rate $T^{1/2}$.

To the best of the author’s knowledge, this paper is the first to develop and
verify a restricted eigenvalue condition that enables theoretical analysis of
$\ell_{1}$-regularized estimators in true high-dimensional (co)integrated
datasets. Since the gram matrix based on integrated time series does not
converge to a deterministic matrix, the commonly used indirect approach of
imposing a restricted eigenvalue condition on the limit of the sample
covariance matrix is not well-suited to high-dimensional settings. Therefore,
we instead verify the condition directly on an appropriately scaled version of
the gram matrix. Outside the realm of time series, a broad line of research
exists that is concerned with direct verification of RE type conditions. In
particular, letting $\bm{X}$ be a $(T\times N)$-matrix with independent rows,
numerous authors have shown that $\frac{1}{T}\bm{X}^{\prime}\bm{X}$ satisfies
the RE condition under Gaussian or sub-Gaussian tail conditions (e.g. Zhou,,
2009; Raskutti et al.,, 2010; Rudelson and Zhou,, 2012; Tropp,, 2015;
Kasiviswanathan and Rudelson,, 2018). To the best of my knowledge, however, no
contributions that deal with dependent rows exist in the current literature,
let alone the case in which the rows are strongly persistent as is the case
when dealing with integrated time series.

The RE condition developed in this paper enables the theoretical analysis of
methods based on $\ell_{1}$-regularization, which we illustrate by deriving
tight error bounds for the lasso when applied to cointegrated regression. We
then use these error bounds to show that consistency can be maintained even
when the number of potentially relevant integrated regressors grows
exponentially in the sample size. Furthermore, we conjecture that the reliance
on matrix concentration inequalities in the proofs offers sufficient
flexibility for prospective future theoretical developments, which we motivate
by providing several suggestions for valuable extensions.

The paper proceeds as follows. Section 2 presents the assumptions on the data,
the restricted eigenvalue and the theoretical validity of the latter. Next, in
Section 3 we apply the results of the preceding section to derive error bounds
and asymptotic consistency for the lasso when applied to a high-dimension
(co)integrated dataset. We conclude in Section 4.

Finally, a word on notation. For an $N$-dimensional vector $\bm{x}$, we define
the $L_{p}$-norm as
$\left\lVert\bm{x}\right\rVert_{p}=\left(\sum_{j=1}^{N}x_{j}^{p}\right)^{1/p}$.
Then, for a $(T\times N)$ matrix $\bm{Y}$, the induced matrix norm is defined
as
$\left\lVert\bm{Y}\right\rVert_{p}=\sup_{\bm{x}\in\mathbb{R}^{N}\setminus\\{0\\}}\left\\{\frac{\left\lVert\bm{Y}\bm{x}\right\rVert_{p}}{\left\lVert\bm{x}\right\rVert_{p}}\right\\}$.
If $\bm{Y}$ is square, we use $\lambda_{i}\left(\bm{Y}\right)$ to denote its
$i$-th eigenvalue and adopt the ordering
$\lambda_{1}\left(\bm{Y}\right)\leq\ldots\leq\lambda_{N}\left(\bm{Y}\right)$.
Letting $A,B$ be two index sets, we define $\bm{x}_{A}$ as the sub-vector of
$\bm{x}$ consisting of the elements indexed by the set $A$ and $\bm{Y}_{A,B}$
corresponds to the sub-matrix of $\bm{Y}$ that contains the rows and columns
indexed by $A$ and $B$, respectively. To further simplify notation, we let
$\bm{Y}_{A}=\bm{Y}_{A,A}$.

##  2 Restricted eigenvalues for integrated time series

In this section, our main theoretical result is derived. In particular, we
derive a finite-sample bound on the probability that the restricted eigenvalue
of the gram matrix lies above a positive, universal constant.

Let $\bm{S}=(\bm{s}_{1},\ldots,\bm{s}_{T})^{\prime}$ denote a $(T\times
N)$-dimensional matrix with random walks of the form
$\bm{s}_{t}=\sum_{s=1}^{t}\bm{\epsilon}_{s}$. The restricted minimum
eigenvalue is defined as

$\kappa(\bm{S},f_{T},s,c_{0}):=\underset{\begin{subarray}{c}J_{0}\subseteq\\{1,\ldots,N\\},\\\
\left\lvert J_{0}\right\rvert\leq
s\end{subarray}}{\text{min}}\underset{\begin{subarray}{c}\bm{x}\neq 0,\\\
\left\lVert\bm{x}_{J_{0}^{c}}\right\rVert_{1}\leq
c_{0}\left\lVert\bm{x}_{J_{0}}\right\rVert_{1}\end{subarray}}{\text{min}}\frac{f_{T}\left\lVert\bm{S}\bm{x}\right\rVert_{2}}{\left\lVert\bm{x}_{J_{0}}\right\rVert_{2}}.$
(1)

We impose the following assumption on the innovations driving the random
walks.

######  Assumption 1.

Assume that
$\bm{\epsilon}_{s}\overset{i.i.d.}{\sim}\mathscr{N}\left(0,\bm{\varSigma}\right)$,
where
$0<c_{\sigma}\leq\lambda_{\min}\left(\bm{\varSigma}\right)\leq\lambda_{\max}\left(\bm{\varSigma}\right)\leq
C_{\sigma}<\infty$.

######  Remark 1.

Assumption 1 imposes the innovations to be independently and identically
distributed according to a multivariate Gaussian distribution. Somewhat
surprisingly, the choice for a Gaussian distribution is not motivated for its
exponential tail decay, although this property is exploited in the results to
follow. Instead, the Gaussian distribution’s invariance to orthonormal
transformations is the key ingredient that enables the derivation of
probabilistic bounds on restricted eigenvalues via matrix concentration
inequalities.

The main theoretical result brought forward in this paper concerns the
behaviour of the minimum restricted eigenvalue of $\bm{S}^{\prime}\bm{S}$.

######  Theorem 1.

Assume that
$\left\lceil\frac{c_{0}^{2}C_{\sigma}s}{c_{\sigma}}\right\rceil<N$. Then,
under Assumption 1, there exists $\kappa_{0}>0$ such that

$\operatorname{\mathbb{P}}\left(\kappa\left(\bm{S},\frac{s\log^{3/4}N}{T},s,c_{0}\right)\geq\kappa_{0}\right)\geq
1-(4T+C_{2}s)\text{e}^{-C_{1}s\log N},$ (2)

where $\kappa_{0}:=\kappa_{0}(c_{0},c_{\sigma},C_{\sigma},C_{1})$, whose
expression is given in (34), is a constant that is independent of $s,N,T$,
polynomially decreases in $c_{0}$, $C_{\sigma}$, $C_{1}$, and polynomially
increases in $c_{\sigma}$.

A direct implication of Theorem 1 is that
$\operatorname{\mathbb{P}}\left(\kappa\left(\bm{S},\frac{s\log^{3/4}N}{T},s,c_{0}\right)\geq\kappa_{0}\right)\to
1$ whenever $T\text{e}^{-s\log N}\to 0$ and $s\text{e}^{-s\log N}\to 0$.
Hence, perhaps somewhat counter-intuitively, the restricted eigenvalue
condition is more likely to hold when $s$ and $N$ are large. The reason
underlying this behaviour is the scaling factor of $s\log^{3/4}N$ included in
the restricted minimum eigenvalue. We stress, however, that this is no free
lunch, as the scaling factor will inversely affect the convergence rates of
$\ell_{1}$-regularized estimators, as will become apparent in the following
section.

##  3 Lasso estimation of cointegrated data

Suppose a researcher is interested in modelling a scalar-valued time series
$y_{t}$ based on an $N$-dimensional vector time series $\bm{x}_{t}$ by
considering the model

$\begin{split}\bm{y}_{t}=\bm{\beta}^{\prime}\bm{x}_{t}+\epsilon_{y,t},\\\
\bm{x}_{t}=\bm{x}_{t-1}+\epsilon_{x,t}.\end{split}$ (3)

The vector $\bm{\beta}$ is assumed to be sparse, with
$S_{\beta}=\\{j:\beta_{j}\neq 0\\}$ and $\left\lvert S_{\beta}\right\rvert=s$.
Hence, $\bm{y}_{t}$ is a unit root non-stationary time series that is
cointegrated with $s$ integrated time series. We stack the innovations in the
$(N+1)$ vector
$\bm{\epsilon}_{t}=(\epsilon_{y,t},\bm{\epsilon}_{x,t}^{\prime})^{\prime}$.

Given that $N$ may be (very) large, we consider estimating the linear model
(3) by the lasso. Define $\bm{y}=(y_{1},\ldots,y_{T})^{\prime}$ and
$\bm{X}=(\bm{x}_{1},\ldots,\bm{x}_{T})$. Then, the lasso estimator for
$\bm{\beta}$ is given by

$\hat{\bm{\beta}}=\ \underset{\bm{\beta}}{\text{arg
min}}\left\lVert\bm{y}-\bm{X}\bm{\beta}\right\rVert_{2}^{2}+\lambda_{T}\left\lVert\bm{\beta}\right\rVert_{1}.$
(4)

Based on Theorem 1, we are able to derive the following error bounds on the
prediction and estimation errors.

######  Theorem 2.

Under Assumption 1, it holds that

$\left\lVert\bm{X}(\hat{\bm{\beta}}-\bm{\beta})\right\rVert_{2}^{2}+\lambda_{T}\left\lVert\hat{\bm{\beta}}-\bm{\beta}\right\rVert_{1}\leq\frac{4\lambda_{T}^{2}s^{3}\log^{3/2}N}{T^{2}\phi_{0}^{2}},$

with probability at least

$\begin{split}&1-2N\left[\emph{exp}\left(-\frac{2\lambda_{T}^{2}}{T^{1+2m}}\right)+2T\emph{exp}\left(-\frac{T^{m-\frac{1}{2}}}{2C_{\sigma}}\right)+2\emph{exp}\left(-\frac{\lambda_{T}}{4C_{\sigma}}\right)\right]-(4T+C_{2}s)e^{-C_{1}s\log
N},\end{split}$ (5)

for all $m>0$, and universal constants $C_{1},C_{2}>0$.

The probabilistic finite-sample error bound in Theorem 2 elucidates the dual
role of the penalty parameter $\lambda_{T}$. Whilst a high degree of
penalization results in a looser error bound, it simultaneously increases the
probability with which this error bound holds. The attentive reader may note
that the probability in (5) contains the lower bound on the probability that
the REC holds from Theorem 1, minus an additional term. The latter term
resembles the probability that $\lambda_{T}$ fails to dominate the empirical
process, and is the decisive factor for determining the minimum amount of
penalization required to ensure consistency. Indeed, based on Theorem 2, the
consistency of the lasso estimator can be established under a suitable choice
of $\lambda_{T}$ and appropriate restrictions on the growth rates of $T,N$,
and $s$.

######  Corollary 1.

Assume that,as $T,N,s\to\infty$, $\frac{T^{1+\xi}\log^{1/2}N}{\lambda_{T}}\to
0$ and $\frac{\log N}{T^{\xi}}\to 0$ for some $\xi>0$, and $\frac{\log
T}{s\log N}\to 0$. Then, under Assumption 1, it holds that

$\left\lVert\hat{\bm{\beta}}-\bm{\beta}\right\rVert_{1}=O_{p}\left(\frac{s^{3}\log^{2}N}{T^{1-\xi}}\right),$

for any $\xi>0$.

Corollary 1 makes the rate of consistency, as well as the sufficient
restrictions on growth rates of $\lambda_{T},T,N$ and $s$ explicit. Most
noteworthy is that consistency may be attained when $N$ grows exponentially in
$T$, thereby providing a first justification for the use of the lasso in ultra
high-dimensional cointegrated applications. Furthermore, we note that when the
number of time series that cointegrate with the dependent variable ($s$) grows
at a slow polynomial rate in $T$, while the total number of candidate time
series $N$ grows at any polynomial rate in $T$, then $\hat{\bm{\beta}}$ almost
attains attains the same $T$-consistency as in the fixed-dimensional setting.

##  4 Conclusion

In this paper, we validate the restricted eigenvalue condition (REC) proposed
by Bickel et al., (2009) in a stationary setting to unit root nonstationary
data. The gram matrix is scaled by $\frac{s^{2}\log^{3/2}N}{T^{2}}$ instead of
the usual $\frac{1}{T^{2}}$ that is familiar from the fixed-dimensional case.
We apply our extended REC to derive finite-sample error bounds on the lasso
when applied to high-dimensional cointegrating regressions. The bounds reveal
that consistency can be attained, even under exponential growth of the number
of time series. In addition, when the number of potentially relevant
integrated regressors grow at any polynomial rate in $T$, while the number of
regressors that cointegrate with the dependent variable grows at a slow
polynomial rate, the estimator is almost $T$-consistent, as is the case in the
classic fixed-dimensional setting.

A natural question that may arise concerns the validity of these results in
the presence of multiple cointegrating relationships. We are optimistic that
the results in this paper can serve useful when deriving the properties of
system estimators such as penalized vector error correction models, and we
consider this an interesting avenue for future research.

##  Appendix A Proofs

###  A.1 Lemmas

######  Lemma 1 (Theorem 1.1 of Tropp, (2012)).

Consider a finite sequence $\\{\bm{X}_{t}\\}$ of independent, random, self-
adjoint matrices with dimension $N$. Assume that each random matrix satisfies

$\bm{X}_{t}\succeq 0,\ \text{ and }\ \lambda_{\max}\left(\bm{X}_{t}\right)\leq
R\ \text{ almost surely}.$

Define

$\mu_{\min}:=\lambda_{\min}\left(\sum_{t}\operatorname{\mathbb{E}}(\bm{X}_{t})\right)\quad\text{and}\quad\mu_{\max}:=\lambda_{\max}\left(\sum_{t}\operatorname{\mathbb{E}}(\bm{X}_{t})\right).$

Then,

$\displaystyle\operatorname{\mathbb{P}}\left(\lambda_{\min}\left(\sum_{t=1}^{T}\bm{X}_{t}\right)\leq(1-\delta)\mu_{\min}\right)\leq
N\left(\frac{e^{-\delta}}{(1-\delta)^{1-\delta}}\right)^{\mu_{\min}/R}\text{
for }\delta\in[0,1]\text{, and}$ (6)
$\displaystyle\operatorname{\mathbb{P}}\left(\lambda_{\max}\left(\sum_{t=1}^{T}\bm{X}_{t}\right)\geq(1+\delta)\mu_{\max}\right)\leq
N\left(\frac{e^{\delta}}{(1+\delta)^{1+\delta}}\right)^{\mu_{\min}/R}\text{
for }\delta\geq 0.$ (7)

######  Lemma 2.

Let $f_{T}=\frac{s\log^{3/4}N}{T}$ and define the set
$\mathscr{A}_{T}(a)=\left\\{f_{T}^{2}\left\lVert\bm{X}^{\prime}\bm{\epsilon}_{y}\right\rVert_{\infty}\geq
a\right\\}$. Then, under Assumption 1, it holds that

$\operatorname{\mathbb{P}}\left(\mathscr{A}_{T}(a)\right)\geq
1-2N\left[\emph{exp}\left(-\frac{2a^{2}T^{3-2m}}{s^{4}\log^{4}N}\right)+2T\emph{exp}\left(-\frac{T^{m-\frac{1}{2}}}{2C_{\sigma}}\right)+2\emph{exp}\left(-\frac{aT^{2}}{4C_{\sigma}s^{2}\log^{2}N}\right)\right],$
(8)

for any $m\in\mathbb{R}$.

###### Proof.

The main argument to prove Lemma 2 relies on a truncation argument in
combination with the Azuma-Hoeffding bound and is inspired by the proof of
Theorem 1 in Kock and Callot, (2015). First, by application of the union
bound,

$\begin{split}&\operatorname{\mathbb{P}}\left(\max_{1\leq i\leq
N}\left\lvert\sum_{t=1}^{T}x_{i,t}\epsilon_{y,t}\right\rvert\geq\frac{a}{f_{T}^{2}}\right)\\\
&\leq\operatorname{\mathbb{P}}\left(\max_{1\leq i\leq
N}\left\lvert\sum_{t=1}^{T}x_{i,t-1}\epsilon_{y,t}\right\rvert\geq\frac{a}{2f_{T}^{2}}\right)+\operatorname{\mathbb{P}}\left(\max_{1\leq
i\leq
N}\left\lvert\epsilon_{i,t}\epsilon_{y,t}\right\rvert\geq\frac{a}{2f_{T}^{2}}\right)\end{split}$
(9)

We bound each of the two RHS terms in (9) individually. For the first term, we
obtain

$\begin{split}&\operatorname{\mathbb{P}}\left(\max_{1\leq i\leq
N}\left\lvert\sum_{t=1}^{T}x_{i,t-1}\epsilon_{y,t}\right\rvert\geq\frac{a}{2f_{T}^{2}}\right)\\\
&\leq\operatorname{\mathbb{P}}\left(\max_{1\leq i\leq
N}\left\lvert\sum_{t=1}^{T}x_{i,t-1}\epsilon_{y,t}\right\rvert\geq\frac{a}{2f_{T}^{2}},\max_{1\leq
i\leq N}\max_{1\leq t\leq T}\left\lvert
x_{i,t-1}\epsilon_{y,t}\right\rvert\leq T^{m}\right)\\\
&\quad+\operatorname{\mathbb{P}}\left(\max_{1\leq i\leq N}\max_{1\leq t\leq
T}\left\lvert x_{i,t-1}\epsilon_{y,t}\right\rvert>T^{m}\right)\\\
&\leq\sum_{i=1}^{N}\left[\operatorname{\mathbb{P}}\left(\left\lvert\sum_{t=1}^{T}x_{i,t-1}\epsilon_{y,t}\right\rvert\geq\frac{a}{2f_{T}^{2}},\max_{1\leq
t\leq T}\left\lvert x_{i,t-1}\epsilon_{y,t}\right\rvert\leq
T^{m}\right)+\operatorname{\mathbb{P}}\left(\max_{1\leq t\leq T}\left\lvert
x_{i,t-1}\epsilon_{y,t}\right\rvert>T^{m}\right)\right]\end{split}$ (10)

To bound the first term on the RHS of (10), we aim to apply the Azuma-
Hoeffding inequality for bounded martingale difference sequences (e.g.
Wainwright,, 2019, Corollary 2.20) . To this end, let $\mathscr{F}_{t}$ denote
the natural filtration with respect to the stochastic process
$\\{\bm{\epsilon}_{j}\\}_{j=-\infty}^{t}$. By the proof of Lemma 4 in Kock and
Callot, (2015), it follows that

$\operatorname{\mathbb{E}}\left(x_{i,t-1}\epsilon_{y,t}\mathbbm{1}\\{\left\lvert
x_{i,t-1}\epsilon_{y,t}\right\rvert\leq T^{m}\\}|\mathscr{F}_{t-1}\right)=0.$

Then, by the Azuma-Hoeffding inequality,

$\begin{split}&\sum_{i=1}^{N}\operatorname{\mathbb{P}}\left(\left\lvert\sum_{t=1}^{T}x_{i,t-1}\epsilon_{y,t}\right\rvert\geq\frac{a}{2f_{T}^{2}},\max_{1\leq
t\leq T}\left\lvert x_{i,t-1}\epsilon_{y,t}\right\rvert\leq T^{m}\right)\\\
&\leq
2N\text{exp}\left(\frac{-2a^{2}}{T^{2m+1}f_{T}^{4}}\right)=2N\text{exp}\left(\frac{-2a^{2}T^{3-2m}}{s^{4}\log^{3}N}\right).\end{split}$
(11)

To bound the second term in (10), we first note that for any two random
variables $x,y$ and constant $c>0$, it holds that

$\operatorname{\mathbb{P}}(\left\lvert
xy\right\rvert>c)\leq\operatorname{\mathbb{P}}(\left\lvert x\right\rvert
c^{k}>c,\left\lvert y\right\rvert\leq
c^{k})+\operatorname{\mathbb{P}}(\left\lvert
y\right\rvert>c^{k})\leq\operatorname{\mathbb{P}}(\left\lvert
x\right\rvert>c^{1-k})+\operatorname{\mathbb{P}}(\left\lvert
y\right\rvert>c^{k}).$

In addition, note that $\bm{x}_{i,t-1}\sim\mathscr{N}(0,t\sigma_{ii}^{2})$ and
$\epsilon_{y,t}\sim\mathscr{N}(0,\sigma_{11}^{2})$ by Assumption 1. Then,
based on a concentration inequality for sub-Gaussian variables (e.g.
Wainwright,, 2019, equation (2.9)), we bound

$\begin{split}\operatorname{\mathbb{P}}\left(\left\lvert
x_{i,t-1}\epsilon_{y,t}\right\rvert>T^{m}\right)&\leq\operatorname{\mathbb{P}}\left(\left\lvert
x_{i,t-1}\right\rvert>T^{mk}\right)+\operatorname{\mathbb{P}}\left(\left\lvert\epsilon_{y,t}\right\rvert>T^{m(1-k)}\right)\\\
&\leq
2\text{exp}\left(-\frac{T^{2mk-1}}{2\sigma_{ii}^{2}}\right)+2\text{exp}\left(-\frac{T^{2m(1-k)}}{2\sigma_{11}^{2}}\right).\end{split}$
(12)

We optimize the rate at which (12) converges to zero by equating the
exponents:

$2mk-1=2m(1-k)\Rightarrow k=\frac{1}{2}+\frac{1}{4m}.$

Plugging in this value for $k$, together with the upper bound
$(\sigma_{ii}^{2},\vee\sigma_{11}^{2})\leq C_{\sigma}$ for all $1\leq i\leq
N+1$ and a union bound, it follows that

$\sum_{i=1}^{N}\operatorname{\mathbb{P}}\left(\max_{1\leq t\leq T}\left\lvert
x_{i,t-1}\epsilon_{y,t}\right\rvert>T^{m}\right)\leq
4NT\text{exp}\left(-\frac{T^{m-\frac{1}{2}}}{2C_{\sigma}}\right).$ (13)

Next, again using the sub-Gaussian concentration inequality, we bound the
second term of (9) as

$\begin{split}\operatorname{\mathbb{P}}\left(\max_{1\leq i\leq
N}\left\lvert\epsilon_{i,t}\epsilon_{y,t}\right\rvert\geq\frac{a}{2f_{T}^{2}}\right)&\leq\sum_{i=1}^{N}\left[\operatorname{\mathbb{P}}\left(\left\lvert\epsilon_{i,t}\right\rvert\geq\left(\frac{a}{2f_{T}^{2}}\right)^{1/2}\right)+\operatorname{\mathbb{P}}\left(\left\lvert\epsilon_{y,t}\right\rvert\geq\left(\frac{a}{2f_{T}^{2}}\right)^{1/2}\right)\right]\\\
&\leq
4N\text{exp}\left(-\frac{a}{4f_{T}^{2}C_{\sigma}}\right)=4N\text{exp}\left(-\frac{aT^{2}}{4C_{\sigma}s^{2}\log^{3/2}N}\right).\end{split}$
(14)

Finally, combining (11), (13) and (14), we bound the RHS of (9) as

$\begin{split}&\operatorname{\mathbb{P}}\left(\max_{1\leq i\leq
N}\left\lvert\sum_{t=1}^{T}x_{i,t}\epsilon_{y,t}\right\rvert\geq\frac{a}{f_{T}^{2}}\right)\\\
&\leq
2N\left[\text{exp}\left(-\frac{2a^{2}T^{3-2m}}{s^{4}\log^{3}N}\right)+2T\text{exp}\left(-\frac{T^{m-\frac{1}{2}}}{2C_{\sigma}}\right)+2\text{exp}\left(-\frac{aT^{2}}{4C_{\sigma}s^{2}\log^{3/2}N}\right)\right].\end{split}$

∎

###  A.2 Proofs of main results

###### Proof of Theorem 1.

We first show that the restricted minimum eigenvalue condition holds for a
transformed matrix. Rewrite $\bm{S}=\bm{U}\bm{E}$, where $\bm{U}$ is a lower
triangular matrix with ones on and below the diagonal and
$\bm{E}=(\bm{\epsilon}_{1},\ldots,\bm{\epsilon}_{T})^{\prime}$. Decompose
$\bm{U}^{\prime}\bm{U}=\bm{V}\bm{\varLambda}\bm{V}^{\prime}$, where
$\bm{\varLambda}=\operatorname{\text{diag}}(\lambda_{1},\ldots,\lambda_{T})$
are the eigenvalues of $\bm{U}^{\prime}\bm{U}$. By Lemma 1 in Akesson and
Lehoczky, (1998), it holds that

$\lambda_{t}^{-1}=4\sin^{2}\left(\frac{\omega_{t}}{2}\right)=2(1-\cos\omega_{t}),$
(15)

with $\omega_{t}=\frac{(2t-1)\pi}{2T+1}$.111The second equality in (15) is
based on the identity $\cos(2\alpha)=1-2\sin^{2}(\alpha)$. Both
$\lambda_{1}=O(T^{2})$ and $\sum_{t=1}^{T}\lambda_{t}=O(T^{2})$, which turns
out to be troublesome when deriving the behaviour of the spectrum.
Accordingly, we consider a transformed matrix

$\bm{S}^{(\phi)^{\prime}}\bm{S}^{(\phi)}=\sum_{t=1}^{T}\lambda_{\phi,t}\tilde{\bm{\epsilon}}_{t}\tilde{\bm{\epsilon}}t^{\prime},$

with

$\lambda^{-1}_{\phi,t}=2(1+\phi-\cos\omega_{t}),$

where $\phi$ is a positive sequence that decreases to zero as
$T,N,s\to\infty$. Note that

$\bm{S}^{\prime}\bm{S}-\bm{S}^{(\phi)^{\prime}}\bm{S}^{(\phi)}=\sum_{t=1}^{T}(\lambda_{t}-\lambda_{\phi,t})\tilde{\bm{\epsilon}}_{t}\tilde{\bm{\epsilon}}_{t}^{\prime}\succeq
0,$ (16)

thus implying that
$\kappa\left(\bm{S},s,c_{0}\right)>\kappa\left(\bm{S}^{(\phi)},s,c_{0}\right)$.
Define the sparse minimum and maximum eigenvalues as

$\begin{split}\phi_{\min}(\bm{S}^{(\phi_{u})},f_{T},u)&=\underset{\bm{x}\in\mathbb{R}^{N}:1\leq\mathscr{M}(\bm{x})\leq
u}{\text{min}}\frac{\left\lVert
f_{T}\bm{S}^{(\phi_{u})}\bm{x}\right\rVert_{2}^{2}}{\left\lVert\bm{x}\right\rVert_{2}^{2}},\\\
\phi_{\max}(\bm{S}^{(\phi_{u})},f_{T},u)&=\underset{\bm{x}\in\mathbb{R}^{N}:1\leq\mathscr{M}(\bm{x})\leq
u}{\text{max}}\frac{\left\lVert
f_{T}\bm{S}^{(\phi_{u})}\bm{x}\right\rVert_{2}^{2}}{\left\lVert\bm{x}\right\rVert_{2}^{2}}.\end{split}$
(17)

By Lemma 4.1(ii) in Bickel et al., (2009), it follows that

$\kappa\left(\bm{S}^{(\phi_{s})},f_{T},s,c_{0}\right)\geq\sqrt{\phi_{\min}(\bm{S}^{(\phi_{s})},f_{T},s+m)}-c_{0}\sqrt{\phi_{\max}(\bm{S}^{(\phi_{s})},f_{T},m)}\sqrt{\frac{s}{m}},$
(18)

where $m\geq s$ and $s+m\leq N$. Accordingly, we aim the provide a lower bound
for (18) by deriving probabilistic upper and lower bounds for the sparse
minimum and maximum eigenvalues, respectively.

Let the collection of all index sets of cardinality $s$ be denoted by

$\mathscr{A}(s):=\left\\{A\subseteq\\{1,\ldots,N\\}:\left\lvert
A\right\rvert=s\right\\},$

and note that its cardinality is upper bounded by

$\left\lvert\mathscr{A}(s)\right\rvert\leq\left(\frac{eN}{s}\right)^{s}=e^{s\log(eN/s)}\leq
e^{s(1-\log s+\log N)}\leq e^{s\log N},$ (19)

for $s\geq 3$. For any $A\in\mathscr{A}(s)$, let
$\bm{S}^{(\phi)}_{A}=\bm{U}^{(\phi)}\bm{E}_{A}$ represent the columns of
$\bm{S}^{(\phi)}$ indexed by $A$. Then, an equivalent definition of the sparse
eigenvalues in (17) is given by

$\begin{split}&\phi_{\min}(\bm{S}^{(\phi_{s})},f_{T},s)=\underset{A\in\mathscr{A}(s)}{\text{min}}\lambda_{\min}\left(f_{T}^{2}\bm{S}_{A}^{(\phi_{s})\prime}\bm{S}_{A}^{(\phi_{s})}\right),\\\
&\phi_{\max}(\bm{S}^{(\phi_{s})},f_{T},s)=\underset{A\in\mathscr{A}(s)}{\text{max}}\lambda_{\max}\left(f_{T}^{2}\bm{S}_{A}^{(\phi_{s})\prime}\bm{S}_{A}^{(\phi_{s})}\right).\end{split}$
(20)

We proceed to derive the stochastic order of the quantities in (20). For any
$A\in\mathscr{A}(s)$ it follows that

$\bm{S}_{A}^{(\phi_{s})\prime}\bm{S}_{A}^{(\phi_{s})}=\bm{E}_{A}^{\prime}\bm{V}\bm{\varLambda}^{(\phi_{s})}\bm{V}^{\prime}\bm{E}_{A}=\tilde{\bm{E}}_{A}^{\prime}\bm{\varLambda}^{(\phi_{s})}\tilde{\bm{E}}_{A}=\sum_{t=1}^{T}\lambda_{\phi_{s},t}\tilde{\bm{\epsilon}}_{A,t}\tilde{\bm{\epsilon}}_{A,t}^{\prime},$

where
$\tilde{\bm{\epsilon}}_{A,t}\overset{i.i.d.}{\sim}\mathscr{N}(\bm{0},\bm{\varSigma}_{A})$
by the rotational invariance of the multivariate normal distribution. Let

$\begin{split}\mu^{\phi_{s}}_{\min}(A)&:=\lambda_{\min}\left(\sum_{t=1}^{T}\lambda_{\phi_{s},t}\operatorname{\mathbb{E}}\left(\tilde{\bm{\epsilon}}_{A,t}\tilde{\bm{\epsilon}}_{A,t}^{\prime}\right)\right)=\lambda_{\min}\left(\bm{\varSigma}_{A}\right)\sum_{t=1}^{T}\lambda_{\phi_{s},t}\\\
\mu^{\phi_{s}}_{\max}(A)&:=\lambda_{\max}\left(\sum_{t=1}^{T}\lambda_{\phi_{s},t}\operatorname{\mathbb{E}}\left(\tilde{\bm{\epsilon}}_{A,t}\tilde{\bm{\epsilon}}_{A,t}^{\prime}\right)\right)=\lambda_{\max}\left(\bm{\varSigma}_{A}\right)\sum_{t=1}^{T}\lambda_{\phi_{s},t},\end{split}$

and the corresponding quantities

$\begin{split}\mu^{s}_{\min}&=\min_{A\in\mathscr{A}(s)}\mu^{\phi_{s}}_{\min}(A),\\\
\mu^{s}_{\max}&=\max_{A\in\mathscr{A}(s)}\mu^{\phi_{s}}_{\max}(A).\end{split}$

In addition, define

$R_{s}=C_{\sigma}^{2}\lambda_{\phi_{s},1}s\left((C_{1}+1)\log^{1/2}N+1\right),$
(21)

for an arbitrary constant $C_{1}>0$. First, we use the truncation argument

$\begin{split}&\operatorname{\mathbb{P}}\left(\min_{A\in\mathscr{A}(s)}\lambda_{\min}\left(\frac{1}{\mu_{\min}(A)}\bm{S}_{A}^{(\phi_{s})^{\prime}}\bm{S}_{A}^{(\phi_{s})}\right)\leq
1-\delta\right)\\\
&\leq\operatorname{\mathbb{P}}\left(\min_{A\in\mathscr{A}(s)}\lambda_{\min}\left(\frac{1}{\mu_{\min}(A)}\bm{S}_{A}^{(\phi_{s})^{\prime}}\bm{S}_{A}^{(\phi_{s})}\right)\leq
1-\delta,\sup_{A\in\mathscr{A}(s)}\sup_{1\leq t\leq
T}\lambda_{\max}\left(\lambda_{\phi_{s},t}\tilde{\bm{\epsilon}}_{A,t}\tilde{\bm{\epsilon}}_{A,t}^{\prime}\right)\leq
R_{s}\right)\\\
&\quad+\operatorname{\mathbb{P}}\left(\sup_{A\in\mathscr{A}(s)}\sup_{1\leq
t\leq
T}\lambda_{\max}\left(\lambda_{\phi_{s},t}\tilde{\bm{\epsilon}}_{A,t}\tilde{\bm{\epsilon}}_{A,t}^{\prime}\right)\geq
R_{s}\right).\\\ \end{split}$ (22)

The first term on the RHS can now be bounded with the use of the Chernoff
bound in Lemma 1. Define the constant

$K_{\delta}=\max\left(\frac{e^{-\delta}}{(1-\delta)^{1-\delta}},\frac{e^{\delta}}{(1+\delta)^{1+\delta}}\right)$

and note that $0<K_{\delta}<1$ for any $0<\delta<1$. Then, Lemma 1 states
that, for each $A\in\mathscr{A}(s)$,

$\begin{split}\operatorname{\mathbb{P}}\left(\lambda_{\min}\left(\frac{1}{\mu_{\min}(A)}\bm{S}_{A}^{(\phi_{s})^{\prime}}\bm{S}_{A}^{(\phi_{s})}\right)\leq
1-\delta,\sup_{1\leq t\leq
T}\lambda_{\max}\left(\lambda_{\phi_{s},t}\tilde{\bm{\epsilon}}_{A,t}\tilde{\bm{\epsilon}}_{A,t}^{\prime}\right)\leq
R_{s}\right)\leq sK_{\delta}^{\mu^{\phi_{s}}_{\min}(A)/R_{s}}.\end{split}$
(23)

For any $s\leq N$, we can bound $\mu^{s}_{\min}$ from below as

$\begin{split}2\mu^{s}_{\min}&\geq\sum_{t=1}^{T}\frac{c_{\sigma}}{1-\cos\omega_{t}+\phi_{s}}\geq\sum_{t=1}^{T}\frac{c_{\sigma}}{\left(\frac{(2t+1)\pi}{2T+1}\right)^{2}+\phi_{s}}\\\
&=\sum_{t=1}^{T}\frac{c_{\sigma}(2T+1)^{2}}{(2t+1)^{2}\pi^{2}+\phi_{s}(2T+1)^{2}}.\end{split}$
(24)

Since $\phi_{s}$ is defined as a positive sequence decreasing to zero in $T$,
we have that $\phi_{s}<1$ for large enough $T$. It follows that

$\begin{split}\mu^{s}_{\min}&\geq\sum_{t=1}^{T}\frac{c_{\sigma}(2T+1)^{2}}{2\left((2t+1)^{2}\pi^{2}+\phi_{s}(2T+1)^{2}\right)}\geq\sum_{t=1}^{T}\frac{c_{\sigma}(2T)^{2}}{2\left((3t)^{2}\pi^{2}+\phi_{s}(3T)^{2}\pi^{2}\right)}\\\
&=\sum_{t=1}^{T}\frac{4c_{\sigma}T^{2}}{18\pi^{2}\left(t^{2}+\phi_{s}T^{2}\right)}\geq\sum_{t=1}^{\left[T\phi_{s}^{1/2}\right]}\frac{4c_{\sigma}T^{2}}{18\pi^{2}\left(t^{2}+\phi_{s}T^{2}\right)}\\\
&\geq\frac{4c_{\sigma}\phi_{s}^{1/2}T^{3}}{36\pi^{2}\phi_{s}T^{2}}=\frac{c_{\sigma}T}{9\pi^{2}\phi_{s}^{1/2}}.\end{split}$
(25)

Next, bound $\mu^{s}_{\min}/R_{s}$ as

$\begin{split}\frac{\mu^{s}_{\min}}{R_{s}}&\geq\frac{c_{\sigma}T}{9\pi^{2}C_{\sigma}^{2}\phi_{s}^{1/2}\lambda_{\phi,1}s\left((C_{1}+1)\log^{1/2}N+1\right)}=\frac{c_{\sigma}T(1-\cos\omega_{1}+\phi_{s})}{9\pi^{2}C_{\sigma}^{2}\phi_{s}^{1/2}s\left((C_{1}+1)\log^{1/2}N+1\right)}\\\
&\overset{(*)}{\geq}\frac{c_{\sigma}T\phi_{s}^{1/2}}{9\pi^{2}C_{\sigma}^{2}(C_{1}+3)s\log^{1/2}N}\overset{(**)}{\geq}-\frac{(C_{1}+1)s\log
N}{\log K_{\delta}},\end{split}$ (26)

where (*) follows from the observation that $2\log^{1/2}N>1$ for any $N\geq 2$
and (**) follows after setting

$\phi_{s}^{1/2}=-\frac{9\pi^{2}(C_{1}+1)(C_{1}+3)C_{\sigma}^{2}s^{2}\log^{3/2}N}{c_{\sigma}T\log
K_{\delta}}.$ (27)

Note that $\phi_{s}^{1/2}$ is positive due to the negative term $\log
K_{\delta}$ in the denominator. Then, combining (25)-(27), we use the Chernoff
matrix-bound from Lemma 1 to bound the first RHS term of (22) as

$\begin{split}&\operatorname{\mathbb{P}}\left(\min_{A\in\mathscr{A}(s)}\lambda_{\min}\left(\frac{1}{\mu^{\phi_{s}}_{\min}(A)}\bm{S}_{A}^{(\phi)^{\prime}}\bm{S}_{A}^{(\phi)}\right)\leq
1-\delta,\sup_{A\in\mathscr{A}(s)}\sup_{1\leq t\leq
T}\lambda_{\max}\left(\lambda_{\phi_{s},t}\tilde{\bm{\epsilon}}_{t}\tilde{\bm{\epsilon}}_{t}^{\prime}\right)\leq
R_{s}\right)\\\
&\leq\sum_{A\in\mathscr{A}(s)}\operatorname{\mathbb{P}}\left(\min_{A\in\mathscr{A}(s)}\lambda_{\min}\left(\frac{1}{\mu^{\phi_{s}}_{\min}(A)}\bm{S}_{A}^{(\phi)^{\prime}}\bm{S}_{A}^{(\phi)}\right)\leq
1-\delta,\sup_{A\in\mathscr{A}(s)}\sup_{1\leq t\leq
T}\lambda_{\max}\left(\lambda_{\phi_{s},t}\tilde{\bm{\epsilon}}_{t}\tilde{\bm{\epsilon}}_{t}^{\prime}\right)\leq
R_{s}\right)\\\ &\leq s\cdot\text{exp}\left(s\log N+\frac{\mu^{s}_{\min}\log
K_{\delta}}{R_{s}}\right)\leq s\cdot\text{exp}\left(-C_{1}s\log
N\right).\end{split}$ (28)

Next, we proceed by bounding the second RHS term of (22). Let
$\bm{v}_{t}\sim\mathscr{N}\left(0,\bm{I}_{\left\lvert A\right\rvert}\right)$,
such that
$\tilde{\bm{\epsilon}}_{A,t}\overset{d}{=}\bm{\varSigma}_{A}^{1/2}\bm{v}_{t}$.
By application of the union bound, the fact that the
$\tilde{\bm{\epsilon}}_{t}$ are identically distributed, the definition of
$R_{s}$ and the triangle inequality, it holds that

$\begin{split}&\operatorname{\mathbb{P}}\left(\sup_{A\in\mathscr{A}(s)}\sup_{1\leq
t\leq
T}\lambda_{\max}\left(\lambda_{\phi_{s},t}\tilde{\bm{\epsilon}}_{A,t}\tilde{\bm{\epsilon}}_{A,t}^{\prime}\right)\geq
R_{s}\right)\leq\sum_{A\in\mathscr{A}(s)}\sum_{t=1}^{T}\operatorname{\mathbb{P}}\left(\lambda_{\max}\left(\tilde{\bm{\epsilon}}_{A,t}\tilde{\bm{\epsilon}}_{A,t}^{\prime}\right)\geq\frac{R_{s}}{\lambda_{\phi_{s},t}}\right)\\\
&\leq
T\sum_{A\in\mathscr{A}(s)}\operatorname{\mathbb{P}}\left(\lambda_{\max}\left(\tilde{\bm{\epsilon}}_{A,1}\tilde{\bm{\epsilon}}_{A,1}^{\prime}\right)\geq\frac{R_{s}}{\lambda_{\phi_{s},1}}\right)=T\sum_{A\in\mathscr{A}(s)}\operatorname{\mathbb{P}}\left(\left\lVert\tilde{\bm{\epsilon}}_{A,1}\right\rVert_{2}^{2}\geq\frac{R_{s}}{\lambda_{\phi_{s},1}}\right)\\\
&\leq
T\sum_{A\in\mathscr{A}(s)}\operatorname{\mathbb{P}}\left(\left\lVert\tilde{\bm{v}}_{1}\right\rVert_{2}^{2}\geq\frac{R_{s}}{\lambda_{\phi_{s},1}\left\lVert\bm{\varSigma}_{A}\right\rVert_{2}^{2}}\right)\leq
T\sum_{A\in\mathscr{A}(s)}\operatorname{\mathbb{P}}\left(\frac{1}{s}\sum_{j\in
A}v^{2}_{j,1}\geq(C_{1}+1)\log^{1/2}N+1\right)\\\ &\leq
T\sum_{A\in\mathscr{A}(s)}\operatorname{\mathbb{P}}\left(\left\lvert\frac{1}{s}\sum_{j\in
A}v^{2}_{j,1}-1\right\rvert\geq(C_{1}+1)\log^{1/2}N\right)\overset{(*)}{\leq}2T\exp\left(-C_{1}s\log
N\right),\end{split}$ (29)

where (*) follows by the sub-exponential concentration of property of
$\chi^{2}$ random variables (see Example 2.11 in Wainwright, (2019)). Finally,
plugging (28) and (29) into (22), we obtain the bound for the minimum sparse
eigenvalue as

$\begin{split}&\operatorname{\mathbb{P}}\left(\phi_{\min}\left(\bm{S}^{(\phi_{s})},\frac{1}{\sqrt{\mu^{s}_{\min}}},s\right)\leq
1-\delta\right)\leq\operatorname{\mathbb{P}}\left(\min_{A\in\mathscr{A}(s)}\lambda_{\min}\left(\frac{1}{\mu^{\phi_{s}}_{\min}(A)}\bm{S}_{A}^{(\phi_{s})^{\prime}}\bm{S}_{A}^{(\phi_{s})}\right)\leq
1-\delta\right)\\\ &\leq(2T+s)\text{exp}\left(-C_{1}s\log
N\right).\end{split}$ (30)

Finally, since the lower bounds in (25) and (26) remain valid for $\mu_{\max}$
and $\frac{\mu_{\max}}{R_{s}}$, respectively, it follows by analogous
reasoning that the maximum sparse eigenvalue is bounded as

$\begin{split}&\operatorname{\mathbb{P}}\left(\phi^{2}_{\max}\left(\bm{S}^{(\phi_{s})},\frac{1}{\sqrt{\mu^{s}_{\max}}},s\right)\geq
1+\delta\right)\leq\operatorname{\mathbb{P}}\left(\max_{A\in\mathscr{A}(s)}\lambda_{\max}\left(\frac{1}{\mu^{\phi_{s}}_{\max}(A)}\bm{S}_{A}^{(\phi_{s})^{\prime}}\bm{S}_{A}^{(\phi_{s})}\right)\geq
1+\delta\right)\\\ &\leq(2T+s)\text{exp}\left(-C_{1}s\log
N\right).\end{split}$ (31)

Having derived probabilistic bounds for the minimum and maximum sparse
eigenvalues, we proceed to prove the restricted eigenvalue condition in
Theorem 1. Define the set

$\mathscr{R}(s,m)=\left\\{\phi_{\min}\left(\bm{S}^{(\phi_{s+m})},\frac{1}{\sqrt{\mu^{s+m}_{\min}}},s+m\right)>1-\delta\right\\}\bigcap\left\\{\phi_{\max}\left(\bm{S}^{(\phi_{s})},\frac{1}{\sqrt{\mu^{s}_{\max}}},s\right)<1+\delta\right\\}.$

For a positive constant $\kappa>0$, define

$m_{\kappa}=\left\lceil\frac{c_{0}^{2}C_{\sigma}(1+\delta)s}{c_{\sigma}(\sqrt{1-\delta}-\kappa)^{2}}\right\rceil=:\left\lceil
C_{\kappa}s\right\rceil,$

with
$C_{\kappa}=\frac{c_{0}^{2}C_{\sigma}(1+\delta)}{c_{\sigma}(\sqrt{1-\delta}-\kappa)^{2}}$.
Then, continuing from (18), it holds on $\mathscr{R}\left(s,m_{\kappa}\right)$
that

$\begin{split}&\kappa\left(\bm{S}^{(\phi_{s})},\frac{1}{\sqrt{\mu^{s+m_{\kappa}}_{\min}}},s,c_{0}\right)\\\
&\geq\sqrt{\phi_{\min}\left(\bm{S}^{(\phi_{s})},\frac{1}{\sqrt{\mu^{s+m_{\kappa}}_{\min}}},s+m_{\kappa}\right)}-c_{0}\sqrt{\phi_{\max}\left(\bm{S}^{(\phi_{s})},\frac{1}{\sqrt{\mu^{s+m_{\kappa}}_{\min}}},s\right)}\sqrt{\frac{s}{m_{\kappa}}}\\\
&\overset{(*)}{\geq}\sqrt{\phi_{\min}\left(\bm{S}^{(\phi_{s+m_{\kappa}})},\frac{1}{\sqrt{\mu^{s+m_{\kappa}}_{\min}}},s+m_{\kappa}\right)}-c_{0}\sqrt{\phi_{\max}\left(\bm{S}^{(\phi_{s})},\frac{1}{\sqrt{\mu^{s+m_{\kappa}}_{\min}}},s\right)}\sqrt{\frac{s}{m_{\kappa}}}\\\
&\geq\sqrt{\phi_{\min}\left(\bm{S}^{(\phi_{s+m_{\kappa}})},\frac{1}{\sqrt{\mu^{s+m_{\kappa}}_{\min}}},s+m_{\kappa}\right)}-c_{0}\sqrt{\frac{\mu^{s}_{\max}}{\mu^{s+m_{\kappa}}_{\min}}}\sqrt{\phi_{\max}\left(\bm{S}^{(\phi_{s})},\frac{1}{\sqrt{\mu^{s}_{\max}}},s\right)}\sqrt{\frac{s}{m_{\kappa}}}\\\
&\overset{(**)}{\geq}\sqrt{1-\delta}-c_{0}\sqrt{\frac{C_{\sigma}(1+\delta)s}{c_{\sigma}m_{\kappa}}}\geq\kappa,\end{split}$
(32)

where ($*$) uses that $\lambda_{\phi_{s},t}>\lambda_{\phi_{s+m},t}$ and (**)
uses that

$\frac{\mu^{s}_{\max}}{\mu^{s+m_{\kappa}}_{\min}}\leq\frac{\mu^{s+m_{\kappa}}_{\max}}{\mu^{s+m_{\kappa}}_{\min}}\leq\frac{C_{\sigma}}{c_{\sigma}}.$

Note that by (25) and (27)

$\frac{1}{\mu^{s+m_{\kappa}}_{\min}}\leq\frac{9\pi^{2}\phi_{s+m_{\kappa}}^{1/2}}{c_{\sigma}T}=-\frac{81\pi^{4}(C_{1}+1)(C_{1}+3)C_{\sigma}^{2}(C_{\kappa}+1)^{2}s^{2}\log^{3/2}N}{c_{\sigma}^{2}T^{2}\log
K_{\delta}}=\frac{C_{\mu}s^{2}\log^{3/2}N}{T^{2}},$

with
$C_{\mu}=-\frac{81\pi^{4}(C_{1}+1)(C_{1}+3)C_{\sigma}^{2}(C_{\kappa}+1)^{2}}{c_{\sigma}^{2}\log
K_{\delta}}$. Then, we conclude that

$\kappa\left(\bm{S},\frac{s\log^{3/4}N}{T},s,c_{0}\right)\geq
C_{\mu}^{-1/2}\kappa\left(\bm{S}^{(\phi_{s})},\frac{1}{\sqrt{\mu^{s+m_{\kappa}}_{\min}}},s,c_{0}\right)\geq
C_{\mu}^{-1/2}\kappa=:\kappa_{0}$ (33)

on the set $\mathscr{R}\left(s,m_{\kappa}\right)$, where the fully expanded
expression for $\kappa_{0}$ reads

$\kappa_{0}:=\kappa_{0}(c_{0},c_{\sigma},C_{\sigma},C_{1})=-\frac{\kappa
c_{\sigma}\log^{1/2}K_{\delta}}{9\pi^{2}\sqrt{(C_{1}+1)(C_{1}+3)}C_{\sigma}\left(\frac{c_{0}^{2}C_{\sigma}(1+\delta)}{c_{\sigma}(\sqrt{1-\delta}-\kappa)^{2}}+1\right)},$
(34)

where $\kappa>0$ and $0<\delta<1$ are free constants chosen such that
$\left\lceil\frac{c_{0}^{2}C_{\sigma}(1+\delta)s}{c_{\sigma}(\sqrt{1-\delta}-\kappa)^{2}}\right\rceil<N$.

Finally, based on (30) and (31),

$\begin{split}\operatorname{\mathbb{P}}\left(\mathscr{R}(s,m_{\kappa})\right)&\geq
1-(2T+s)\text{e}^{-C_{1}s\log N}-(2T+m_{\kappa})\text{e}^{-C_{1}m_{\kappa}\log
N}\\\ &\geq 1-(4T+(C_{k}+1)s)\text{e}^{-C_{1}s\log
N}=1-(4T+C_{2}s)\text{e}^{-C_{1}s\log N},\end{split}$

with $C_{2}=C_{\kappa}+1$. This completes the proof of Theorem 1.

∎

###### Proof of Theorem 2.

The proof of the error bound is based on standard arguments, see for example
Theorem 6.1 in Bühlmann and Van De Geer, (2011), but the probability on which
it holds requires new arguments. First, let $f_{T}=\frac{s\log^{3/4}N}{T}$ and
define $\tilde{\lambda}_{T}=f_{T}^{2}\lambda_{T}$, such that the estimator can
equivalently be defined as the minimizer of

$\mathscr{L}(\bm{\beta})=f_{T}^{2}\left\lVert\bm{y}-\bm{X}\bm{\beta}\right\rVert_{2}^{2}+\tilde{\lambda}_{T}\left\lVert\bm{\beta}\right\rVert_{1}.$
(35)

Next, define the set

$\mathscr{A}_{T}\left(\tilde{\lambda}_{2}\right)=\left\\{f_{T}^{2}\left\lVert\bm{X}^{\prime}\bm{\epsilon}_{y}\right\rVert_{\infty}\leq\frac{\tilde{\lambda}_{T}}{4}\right\\}.$
(36)

By construction
$\mathscr{L}\left(\hat{\bm{\beta}}\right)\leq\mathscr{L}\left(\bm{\beta}\right)$.
Then, on the set $\mathscr{A}_{T}$, it holds that

$\begin{split}&f_{T}^{2}\left\lVert\bm{y}-\bm{X}\hat{\bm{\beta}}\right\rVert_{2}^{2}+\tilde{\lambda}_{T}\left\lVert\hat{\bm{\beta}}\right\rVert_{1}\leq
f_{T}^{2}\left\lVert\bm{y}-\bm{X}\bm{\beta}\right\rVert_{2}^{2}+\tilde{\lambda}_{T}\left\lVert\bm{\beta}\right\rVert_{1}\\\
&\begin{split}\Rightarrow
f_{T}^{2}\left\lVert\bm{X}(\hat{\bm{\beta}}-\bm{\beta})\right\rVert_{2}^{2}+\tilde{\lambda}_{T}\left\lVert\hat{\bm{\beta}}\right\rVert_{1}&\leq
2f_{T}^{2}\left\lVert\bm{X}^{\prime}\bm{\epsilon}_{y}\right\rVert_{\infty}\left\lVert\hat{\bm{\beta}}-\bm{\beta}\right\rVert_{1}+\tilde{\lambda}_{T}\left\lVert\bm{\beta}_{S}\right\rVert_{1}\\\
&\leq\frac{\tilde{\lambda}_{T}}{2}\left\lVert\hat{\bm{\beta}}-\bm{\beta}\right\rVert_{1}+\tilde{\lambda}_{T}\left\lVert\bm{\beta}_{S}\right\rVert_{1}\end{split}\end{split}.$
(37)

Using that

$\left\lVert\hat{\bm{\beta}}\right\rVert_{1}=\left\lVert\hat{\bm{\beta}}_{S}\right\rVert_{1}+\left\lVert\hat{\bm{\beta}}_{S}{{}^{c}}\right\rVert_{1}\geq\left\lVert\hat{\bm{\beta}}_{S}\right\rVert_{1}-\left\lVert\hat{\bm{\beta}}_{S}-\bm{\beta}_{S}\right\rVert_{1}+\left\lVert\hat{\bm{\beta}}_{S^{c}}\right\rVert_{1},$

it follows from (37) that

$\begin{split}2\tilde{\lambda}_{T}\left\lVert\hat{\bm{\beta}}\right\rVert_{1}\leq\tilde{\lambda}_{T}\left\lVert\hat{\bm{\beta}}-\bm{\beta}\right\rVert_{1}+2\tilde{\lambda}_{T}\left\lVert\bm{\beta}_{S}\right\rVert_{1}\Rightarrow\left\lVert\hat{\bm{\beta}}_{S^{c}}\right\rVert_{1}\leq
3\left\lVert\hat{\bm{\beta}}_{S}-\bm{\beta}_{S}\right\rVert_{1}.\end{split}$
(38)

Define the cone

$\mathscr{C}_{N}(S,c_{0}):=\left\\{\bm{x}\in\mathbb{R}^{N}:\left\lVert\bm{x}_{S^{c}}\right\rVert_{1}\leq
c_{0}\left\lVert\bm{x}_{S}\right\rVert_{1}\right\\}$ (39)

and note that (38) implies that
$\hat{\bm{\beta}}-\bm{\beta}\in\mathscr{C}_{N}(S,3)$. Recall the definition of
$\kappa\left(\hat{\bm{\varSigma}}_{T},s,c_{0}\right)$ in (1) and define the
set

$\mathscr{B}_{T}=\left(\kappa\left(\hat{\bm{\varSigma}}_{T},s,c_{0}\right)\geq\phi_{0}\right\\}$
(40)

assume that $\kappa\left(\hat{\bm{\varSigma}}_{T},s,3\right)\geq\phi_{0}^{2}$.
Then, combining (37) and (38),

$\begin{split}&2\left\lVert
f_{T}\bm{X}(\hat{\bm{\beta}}-\bm{\beta})\right\rVert_{2}^{2}+\tilde{\lambda}_{T}\left\lVert\hat{\bm{\beta}}-\bm{\beta}\right\rVert_{1}\overset{\eqref{eq:bound_step1}}{\leq}4\tilde{\lambda}_{T}\left\lVert\hat{\bm{\beta}}_{S}-\bm{\beta}_{S}\right\rVert_{1}\leq
4\sqrt{s}\tilde{\lambda}_{T}\left\lVert\hat{\bm{\beta}}_{S}-\bm{\beta}_{S}\right\rVert_{2}\\\
&\leq\frac{4\tilde{\lambda}_{T}\sqrt{s}\left\lVert
f_{T}\bm{X}\left(\hat{\bm{\beta}}-\bm{\beta}\right)\right\rVert_{2}}{\phi_{0}}\leq\frac{4\tilde{\lambda}_{T}^{2}s}{\phi_{0}^{2}}+\left\lVert
f_{T}\bm{X}\left(\hat{\bm{\beta}}-\bm{\beta}\right)\right\rVert_{2}^{2},\end{split}$
(41)

with the last inequality following from the elementary inequality $2uv\leq
u^{2}+v^{2}$. Hence, we are left to conclude that, on
$\mathscr{A}_{T}\left(\tilde{\lambda}_{T}\right)\cap\mathscr{B}_{T}$,

$\begin{split}&\left\lVert
f_{T}\bm{X}(\hat{\bm{\beta}}-\bm{\beta})\right\rVert_{2}^{2}+\tilde{\lambda}_{T}\left\lVert\hat{\bm{\beta}}-\bm{\beta}\right\rVert_{1}\leq\frac{4\tilde{\lambda}_{T}^{2}s}{\phi_{0}^{2}}\\\
&\Leftrightarrow\left\lVert\bm{X}(\hat{\bm{\beta}}-\bm{\beta})\right\rVert_{2}^{2}+\lambda_{T}\left\lVert\hat{\bm{\beta}}-\bm{\beta}\right\rVert_{1}\leq\frac{4f_{T}^{2}\lambda_{T}^{2}s}{\phi_{0}^{2}}=\frac{4\lambda_{T}^{2}s^{3}\log^{3/2}N}{T^{2}\phi_{0}^{2}}\end{split}$

Finally, we note that
$\operatorname{\mathbb{P}}\left(\mathscr{A}_{T}\left(\tilde{\lambda}_{T}\right)\cap\mathscr{B}_{T}\right)\leq
1-\operatorname{\mathbb{P}}\left(\mathscr{A}_{T}\left(\tilde{\lambda}_{T}\right)^{c}\right)-\operatorname{\mathbb{P}}\left(\mathscr{B}_{T}^{c}\right)$.
By Lemma 2, it holds that

$\begin{split}&\operatorname{\mathbb{P}}\left(\mathscr{A}_{T}^{c}\right)\leq
2N\left[\text{exp}\left(-\frac{2\tilde{\lambda}_{T}^{2}T^{3-2m}}{s^{4}\log^{3}N}\right)+2T\text{exp}\left(-\frac{T^{m-\frac{1}{2}}}{2C_{\sigma}}\right)+2\text{exp}\left(-\frac{\tilde{\lambda}_{T}T^{2}}{4C_{\sigma}s^{2}\log^{3/2}N}\right)\right]\\\
&=2N\left[\text{exp}\left(-\frac{2\lambda_{T}^{2}}{T^{1+2m}}\right)+2T\text{exp}\left(-\frac{T^{m-\frac{1}{2}}}{2C_{\sigma}}\right)+2\text{exp}\left(-\frac{\lambda_{T}}{4C_{\sigma}}\right)\right],\end{split}$

whereas it follows from Theorem 1 that

$\operatorname{\mathbb{P}}\left(\mathscr{B}_{T}^{c}\right)\leq(4T+C_{2}s)\text{e}^{-C_{1}s\log
N}.$

Combining the latter two inequalities completes the proof. ∎

###### Proof of Corollary 1.

We show that under the the rate restrictions in Corollary 1, and Assumption 1,
the error bound in Theorem 2 holds with probability converging to 1 as
$T,N,s\to\infty$. This requires us to show that

$\begin{split}&2N\left[\text{exp}\left(-\frac{2\lambda_{T}^{2}}{T^{1+2m}}\right)+2T\text{exp}\left(-\frac{T^{m-\frac{1}{2}}}{2C_{\sigma}}\right)+2\text{exp}\left(-\frac{\lambda_{T}}{4C_{\sigma}}\right)\right]+z(4T+C_{2}s)e^{-C_{1}s\log
N}\\\ &=A_{1}+A_{2}+A_{3}+A_{4}\to 0,\end{split}$

as $T,N,s\to\infty$. For $A_{1}$, we set $m=1/2+\xi$, and note that

$A_{1}=\text{exp}\left(-\frac{2\lambda_{T}^{2}}{T^{2+2\xi}}+\log N\right)\to
0$

if $\frac{T^{1+\xi}\log^{1/2}N}{\lambda_{T}}\to 0$ as $T,N,s\to\infty$. Next,
note that

$A_{2}=\text{exp}\left(-\frac{T^{\xi}}{2C_{\sigma}}+\log N+\log 2T\right)\to
0,$

if $\frac{\log N}{T^{\xi}}\to 0$. Moving to $A_{3}$, it holds that

$A_{3}=2\text{exp}\left(-\frac{\lambda_{T}}{4C_{\sigma}}+\log N\right)\to 0,$

if $\frac{\log N}{\lambda_{T}}\to 0$ as $T,N,s\to\infty$, which is also
implied by the combination of conditions that
$\frac{T^{1+\xi}\log^{1/2}N}{\lambda_{T}}\to 0$ and $\frac{\log N}{T^{\xi}}\to
0$ as $T,N,s\to\infty$. Finally, it holds that

$A_{4}=(4T+C_{2}s)\text{e}^{-C_{1}s\log N}=4\text{exp}\left(-C_{1}s\log N+\log
T\right)+C_{2}\text{exp}\left(-C_{1}s\log N+\log s\right)\to 0,$

if $\frac{s\log N}{\log T}\to\infty$ and $\frac{s\log N}{\log s}\to\infty$, as
$T,N,s\to\infty$. The latter condition holds trivially, while the former
condition is implied by the condition $\frac{\log T}{s\log N}\to 0$ as
$T,N,s\to\infty$. As a result, it follows that $\mathscr{P}(T,N,s)\to 1$ as
$T,N,s\to\infty$, with $\mathscr{P}(T,N,s)$ being defined in (5), thereby
completing the proof. ∎

## References

  * Akesson and Lehoczky, (1998) Akesson, F. and Lehoczky, J. (1998).  Discrete eigenfunction expansion of multidimensional Brownian motion and the Ornstein-Uhlenbeck process.  Preprint.
  * Basu and Michailidis, (2015) Basu, S. and Michailidis, G. (2015).  Regularized estimation in sparse high-dimensional time series models.  The Annals of Statistics, 43:1535–1567.
  * Bickel et al., (2009) Bickel, P. J., Ritov, Y., and Tsybakov, A. B. (2009).  Simultaneous analysis of lasso and dantzig selector.  The Annals of Statistics, 37:1705–1732.
  * Bühlmann and Van De Geer, (2011) Bühlmann, P. and Van De Geer, S. (2011).  Statistics for High-Dimensional Data: Methods, Theory and Applications.  Springer Science & Business Media.
  * Kasiviswanathan and Rudelson, (2018) Kasiviswanathan, S. P. and Rudelson, M. (2018).  Restricted eigenvalue from stable rank with applications to sparse linear regression.  In Conference On Learning Theory, pages 1011–1041. PMLR.
  * Kock and Callot, (2015) Kock, A. B. and Callot, L. (2015).  Oracle inequalities for high dimensional vector autoregressions.  Journal of Econometrics, 186:325–344.
  * Koo et al., (2020) Koo, B., Anderson, H. M., Seo, M. H., and Yao, W. (2020).  High-dimensional predictive regression in the presence of cointegration.  Journal of Econometrics, 219:456–477.
  * Lee et al., (2021) Lee, J. H., Shi, Z., and Gao, Z. (2021).  On lasso for predictive regression.  Journal of Econometrics.
  * Liang and Schienle, (2019) Liang, C. and Schienle, M. (2019).  Determination of vector error correction models in high dimensions.  Journal of Econometrics, 208(2):418–441.
  * Liao and Phillips, (2015) Liao, Z. and Phillips, P. C. B. (2015).  Automated estimation of vector error correction models.  Econometric Theory, 31:581–646.
  * Masini et al., (2019) Masini, R. P., Medeiros, M. C., and Mendes, E. F. (2019).  Regularized estimation of high-dimensional vector autoregressions with weakly dependent innovations.  arXiv e-prints 1912.09002, arXive.
  * Medeiros and Mendes, (2016) Medeiros, M. C. and Mendes, E. F. (2016).  $\ell_{1}$-regularization of high-dimensional time series models with non-gaussian and heteroskedastic errors.  Journal of Econometrics, 191:255–271.
  * Raskutti et al., (2010) Raskutti, G., Wainwright, M. J., and Yu, B. (2010).  Restricted eigenvalue properties for correlated gaussian designs.  The Journal of Machine Learning Research, 11:2241–2259.
  * Rudelson and Zhou, (2012) Rudelson, M. and Zhou, S. (2012).  Reconstruction from anisotropic random measurements.  In Conference on Learning Theory, pages 10–1. JMLR Workshop and Conference Proceedings.
  * Smeekes and Wijler, (2020) Smeekes, S. and Wijler, E. (2020).  An automated approach towards sparse single-equation cointegration modelling.  Journal of Econometrics.
  * Tropp, (2012) Tropp, J. A. (2012).  User-friendly tail bounds for sums of random matrices.  Foundations of Computational Mathematics, 12:389–434.
  * Tropp, (2015) Tropp, J. A. (2015).  Convex recovery of a structured signal from independent random linear measurements.  In Sampling Theory, a Renaissance, pages 67–101. Springer.
  * Wainwright, (2019) Wainwright, M. J. (2019).  High-dimensional statistics: A non-asymptotic viewpoint, volume 48.  Cambridge University Press.
  * Wong et al., (2020) Wong, K. C., Li, Z., and Tewari, A. (2020).  Lasso guarantees for $\beta$-mixing heavy-tailed time series.  Annals of Statistics, 48:1124–1142.
  * Zhou, (2009) Zhou, S. (2009).  Restricted eigenvalue conditions on subgaussian random matrices.  arxiv e-prints, arXive.


# Simulation of the Dissipative Dynamics of Strongly Interacting NV Centers
with Tensor Networks

Jirawat Saiphet<EMAIL_ADDRESS>Daniel Braun  Institut für
Theoretische Physik, Eberhard Karls Universität Tübingen, 72076 Tübingen,
Germany

###### Abstract

NV centers in diamond are a promising platform for highly sensitive quantum
sensors for magnetic fields and other physical quantities. The quest for high
sensitivity combined with high spatial resolution leads naturally to dense
ensembles of NV centers, and hence to strong, long-range interactions between
them. Hence, simulating strongly interacting NVs becomes essential. However,
obtaining the exact dynamics for a many-spin system is a challenging task due
to the exponential scaling of the Hilbert space dimension, a problem that is
exacerbated when the system is modelled as an open quantum system. In this
work, we employ the Matrix Product Density Operator (MPDO) method to represent
the many-body mixed state and to simulate the dynamics of an ensemble of NVs
in the presence of strong long-range couplings due to dipole-dipole forces. We
benchmark different time-evolution algorithms in terms of numerical accuracy
and stability. Subsequently, we simulate the dynamics in the strong
interaction regime, with and without dissipation.

##  I Introduction

Probing magnetic fields with high sensitivity and resolution is important in
frontier research applications. A single NV-center (NV for short) in diamond
has been proprosed as a nanoscale probe [1] and used for measuring a magnetic
field [2, 3]. Recently, controlled systems of double and triple NV centers
were successfully fabricated [4]. Using NV ensembles with many spins has the
potential to increase the sensitivity by having more spins in a probe. [5, 6,
7]. Requesting at the same time high spatial resolution leads to ensembles
with high density. The resulting strong dipole-dipole interactions between the
NVs lead, however, to a rapid population of sub-spaces of Hilbert space with
reduced total spin. This can be considered a form of intrinsic decoherence [8]
in addition to the remaining external decoherence mechanisms, resulting in
short coherence times [9] and hence less sensitivity. Decoupling the
interactions with specific control pulses, or engineering the alignments of
NVs to reduce interactions can increase the coherence time and sensitivity [8,
10, 11]. However, ideally one would like to profit from the inteactions for
generating entangled states that could highly enhance the sensitivity of the
probe.

To model and optimize entanglement generation in a dissipative system with
strong and long-range interactions, simulations of its dynamics in the
presence of microwave control pulses is necessary. However, given the many-
body nature of an ensemble, simulating its exact dynamics is intractable.
Exact simulations of closed spin-1/2 systems with long-range interaction have
been implemented up to 32 spins [12, 13, 8, 14], and up to 12 spins with
dissipation [15, 16]. In order to address this challenge, we use a tensor
network approach to capture the dynamics of the ensemble. Matrix Product
States (MPS) [17, 18, 19, 20] were proposed for the efficient simulation of
quantum metrology [21] and open quantum systems [22]. Simulation of non-
Markovian systems has been performed using Matrix Product Operators (MPO) with
a nearest-neighbers model [23, 24].

In this work, we consider an NV ensemble that consists of spin-1 particles.
All NVs interact with each other by long-range dipole-dipole interaction. We
simulate the dissipative dynamics by using the Matrix Product Density Operator
(MPDO) method. We investigate the efficiency of using MPDO in simulating
dynamics in the strong interaction limit and under dissipation. Operator
entanglement entropy (opEE) is computed and used to demonstrate the interplay
between strong interaction and dissipation for the capability of the MPDO to
approximate the exact states. We then address and quantify possible
sensitivity improvements from NV-NV interactions during the time evolution
using quantum Fisher information.

##  II Theory

###  II.1 Strongly interacting NVs

An NV-center is a spin-1 system that froms in a carbon lattice of diamonds if
two adjacent carbons are replaced by a nitrogen atom and a vacancy. Without an
external field applied, the diamond structure creates a Zero Field Splitting
(ZFS). Energy levels corresponding to states $|0\rangle$ and $|\pm 1\rangle$
are seperated by $D=2\pi\times 2870$ MHz.

Diamond’s crystalographic axes provide 4 possible orientations of the princial
axes of the NV, seperating them into different groups. Different groups have
different couplings to a magnetic field, which leads to different energy
levels that can be distinguished by optically detected magnetic resonance
(ODMR).

In an external magnetic field, the field component $B_{z}$ along the NV’s
princial axis creates an additional energy splitting between $|-1\rangle$ and
$|+1\rangle$ proportional to $B_{z}$. This allows selective transitions from
the ground state to one of the two exited states by applying a microwave field
at resonance frequency. Here we restrict ourselves to transitions to
$|-1\rangle$.

(a)

(b)

(c)

Figure 1: (1(a)) Energy levels of an NV center. (1(b)) 4 possible orientations
given by diamond’s crystalographic axes. (1(c)) ODMR spectra of 4 different
axes from measurements [4] and analytic calculations.

The Hamiltonian for an individual NV with microwave drive with Rabi frequency
$\Omega(t)$ is given by

$\hat{H}_{NV,i}=\hbar\left(D\hat{S}_{z}^{2}+\gamma
B^{(i)}_{z}\hat{S}_{z}+\Omega(t)\cos(\omega t)\hat{S}_{x}^{(i)}\right),$ (1)

where the spin operators for $S=1$ are

$\displaystyle\hat{S}_{x}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&1&0\\\ 1&0&1\\\
0&1&0\end{pmatrix},\hat{S}_{y}=\frac{i}{\sqrt{2}}\begin{pmatrix}0&-1&0\\\
1&0&-1\\\ 0&1&0\end{pmatrix},$
$\displaystyle\hat{S}_{z}=\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{pmatrix}.$ (2)

Applying the unitary $\hat{U}=e^{i\omega t\hat{S}_{z}^{2}}$ to Eq. 1 to
transform to a frame co-rotating with the microwave and using rotating-wave
approximation (RWA) yields

$\hat{H}_{\text{NV},i}=\hbar\left((D-\omega)\hat{S}_{z}^{2}+\gamma
B^{(i)}_{z}\hat{S}_{z}+\frac{\Omega(t)}{2}\hat{S}_{x}^{(i)}\right).$ (3)

At resonance frequency, $\omega=D+\gamma B^{(i)}_{z}$, the microwave drive
makes a transition between $|0\rangle\leftrightarrow|-1\rangle$ for the $i$-th
NV.

NVs interact with each other via dipole-dipole interaction. We consider a case
of strong interaction between NVs and hence ignore nuclear spin. The general
definition of dipole-dipole interaction between NV-$i$ and NV-$j$ is

$\displaystyle H_{\text{dip},ij}=-\frac{\mu_{0}\gamma^{2}\hbar^{2}}{4\pi
r_{ij}^{3}}(3(\vec{S}^{(i)}\dotproduct\hat{r})(\vec{S}^{(j)}\dotproduct\hat{r})-\vec{S}^{(i)}\vec{S}^{(j)}),$
(4)

where $\vec{S}^{(i)}=(\hat{S}_{x},\hat{S}_{y},\hat{S}_{z})^{(i)}$ are spin
operators of NV-$i$, and $\hat{r}_{ij}=(r_{x},r_{y},r_{z})_{ij}$ are the unit
vectors connecting the two NVs. Taking into account groups with different
princial axes, the dipole-dipole interaction in the rotating frame can be
tranformed into an effective Hamiltonian [16],

^H_eff,{ij} = {Cdip(^Sx(i)^Sx(j)\+ ^Sy(i)^Sy(j)\- ^Sz(i)^Sz(j)), same group-
Cdip^Sz(i)^Sz(j), different group

(5)

where $C_{\text{dip}}=\frac{J_{0}q_{ij}}{r_{ij}^{3}}$,
$J_{0}=\frac{\mu_{0}h^{2}\gamma^{2}}{4\pi}=2\pi\times 52$ MHz$\cdot$ nm3,
$q_{ij}=3(\hat{r}_{ij}\cdot\hat{z}_{i})(\hat{r}_{ij}\cdot\hat{z}_{j})-\hat{z}_{i}\cdot\hat{z}_{j}$.
It reduces to $q_{ij}=(3\cos^{2}\theta-1)$ and
$\cos\theta\equiv\hat{z}_{i}\cdot\hat{r}_{ij}$ in the case of the same group
($\hat{z}_{i}=\hat{z}_{j}$).

###  II.2 QFI and Cramér-Rao bound

According to the Cramér-Rao bound, sensitivity in the estimation of a
parameter $\theta$ using quantum probes is bounded by the inverse of the
Quantum Fisher Information (QFI), $(\delta\theta)^{2}\geq\frac{1}{MF_{Q}}$
where $M$ is the number of independent measurements [25]. Hence, sensitivities
can be improved by repeating the measurements or using more probes. For $N$
independent probes, we obtain the Standard Quantum Limit (SQL),
$\delta\theta\geq\mathcal{O}(1/\sqrt{N})$. However, the QFI can be increased
when the probes are highly engtangled. For example, under evolution with a
pure Zeeman term, and in the absence of decoherence and dissipation, probes
prepared in a GHZ state achieve optimal sensitivity for magnetic field
measurement that follows the ”Heisenberg limit” (HL),
$\delta\theta\geq\mathcal{O}(1/N)$ [26].

Thus, since interactions are necessary for the creation of entanglement, dense
NV ensembles harbor the potential for higher sensitivity compared to non-
interacting NVs. However, the lack of permutational symmetry leads to the
population of other irreducible representations of SU(2) starting from the one
with maximum spin, i.e. on average the total spin decays and sensitivity is
reduced. Optimal control is therefore necessary to harvest the entanglement
from strong interactions for achieving higher sensitivity.

###  II.3 Tensor network state

Due to the $r^{-3}$ scaling, the interactions within an ensemble of NVs
extends beyond the nearest-neighbers and becomes long-range. These all-to-all
interactions increase complexity and limit our capability to exactly simulate
the system to only a few spins.

(a) MPS

(b) MPO

(c) MPDO

Figure 2: Tensor diagrams representing (2(a)) a wavefunction for a pure state,
(2(b)) an operator, and (2(c)) a vectorized density operator constructed from
an MPS and its conjugate by contracting auxiliary indices, depicted as dash
lines, and combining corresponding physical (bonds) indices.

To simulate the many-body dynamics for an ensemble of NVs, we represent the
quantum state as a tensor network state. In Fig. 2(a) the state of a closed
system is decompsed into a 1-dimensional tensor network structure, called
Matrix Product State (MPS) [17, 18, 19, 20],

$|\Psi\rangle=\sum_{i_{1}...i_{N}=1}^{d}A^{i_{1}}A^{i_{2}}...A^{i_{N}}|i_{1}\cdot\cdot\cdot
i_{N}\rangle,$ (6)

where a matrix $A^{i_{n}}\equiv A^{i_{n}}_{\alpha_{(n-1)}\alpha_{n}}$. Each
$A^{i_{n}}$ in the MPS contains a physical index, $i_{n}$ representing the
local Hilbert space of the $n^{th}$ NV. The virtual index, or bond index
$\alpha_{n}$ having bond dimensions $d(\alpha_{n})=\chi_{n}$ labels links
between two tensors. Physically, the bond dimension contains information about
the entanglement entropy between the two parts of the tensor network that the
bond connects. Using MPS allows us to compress the bond dimensions and to
efficiently represent the ground state in compact, small Hibert spaces.

A similar tensor network structure can be adapted to simulate an open system
[27, 28]. To represent an operator we need a Matrix Product Operator (MPO) as
given in Fig. 2(b). This MPO, in an orthogonal basis, represents a density
operator of the system,

$\hat{\rho}=\sum_{i_{1}...i_{N}=1}^{d}\sum_{i^{\prime}_{1}...i^{\prime}_{N}=1}^{d}B^{i_{1}i^{\prime}_{1}}B^{i_{2}i^{\prime}_{2}}...B^{i_{N}i^{\prime}_{N}}|i_{1}\cdot\cdot\cdot
i_{N}\rangle\langle i^{\prime}_{1}\cdot\cdot\cdot i^{\prime}_{N}|.$ (7)

For pure states, this MPO can be contructed by contracting auxiliary indices
of two MPS and then combining corresponding physical (bond) indices.
Furthermore, we combine the physical indices to create a Matrix Product
Density Operator (MPDO) representing a vectorized density operator.

$|\rho\rangle\rangle=\sum_{j_{1}...j_{N}=1}^{d^{2}}B^{j_{1}}B^{j_{2}}...B^{j_{N}}|j_{1}\cdot\cdot\cdot
j_{N}\rangle\rangle.$ (8)

Here, $j_{n}=(i_{n},i^{\prime}_{n})$ and $|j_{n}\rangle\in\mathbb{C}^{d^{2}}$
is combined from two physical indices. The dimension of each resulting index
is doubled compared to the MPS. A graphical representation of this process is
shown in Fig. 2(c)

###  II.4 Simulation of dissipative dynamics

We simulate $\rho$ directly, including dissipation, by solving the vectorized
master equation

$\frac{\partial}{\partial
t}|\rho\rangle\rangle=\mathcal{L}|\rho\rangle\rangle,$ (9)

where $|\rho\rangle\rangle$ is the MPDO given by Eq. 8, and $\mathcal{L}$ is a
vectorized Lindblad operator defined as

$\displaystyle\mathcal{L}(t)=$
$\displaystyle-i(\hat{H}(t)\otimes\mathbb{I}-\mathbb{I}\otimes\hat{H}^{T}(t))$
$\displaystyle+\sum_{\nu}\gamma_{i}[L_{i}\otimes(L^{\dagger}_{i})^{T}-\frac{1}{2}(L^{\dagger}_{i}L_{i}\otimes\mathbb{I}+\mathbb{I}\otimes(L^{\dagger}_{i}L_{i})^{T})].$
(10)

Note that when $\gamma_{i}=0$, the dynamics are unitary. In this case when MPS
is sufficient to simulate the system, utilizing MPDO is unnecessary and
computationally more expensive due to the squared memory.

##  III Results

We model an ensemble of NVs as a 1-dimensional spin-1 chain with long-range
interaction. We assume uniform separation, $r_{i,i+1}=r$, for nearest neighbor
spins. Note that an ensemble with any real-space configuration can be mapped
to this model with different interaction strengths. For the very strong
interaction regime, this separation is set to be $r<2$ nm while an actual
sample in experiment has $r\sim 5$ nm [4]. In our simulations we use the
following conditions: i) The amplitude of external magnetic field
$|B_{\text{ext}}|=9$ mT. ii) $\hat{B}_{\text{ext}}$ has a direction calculated
from the ODMR data in Fig. 1(c). iii) Uniform Rabi frequency
$\Omega(t)=\Omega=2\pi\times 2.00$ MHz. iv) All NVs belong to the same
orientation group parallel to [111]. (v) Time step $dt=1$ ns.

###  III.1 Simulation algorithms in strong interaction regime

Firstly, we investigate the precision and numerical stability of tensor
network algorithms in the very strong interaction regime, for two different
algorithms that simulate time-evolution and support the long-range model [29]:
i) MPO $W^{II}$ [30] and ii) time-dependent variational principle (TDVP) [31,
32]. We simulate dynamics of ensembles up to 10 NVs and
$\chi_{\text{max}}=16$. Here $\chi_{\text{max}}$ is a maximum dimension for
all bonds. Fig. 3(a) shows that, when $r=2$ nm, both algorithms can simulate
the dynamics and have similar expectation values of observables, with only
small discrepency. However, the TDVP has significantly smaller errors in the
imaginary parts, as shown in Fig. 3(b). When the interaction strength becomes
large the $W^{II}$ also has numerical stability issues and diverging
expectation value of a local observable $\langle S^{(1)}_{z}\rangle$, while
the TDVP is still valid compared to exact diagonalization. Fig. 3(c) and Fig.
3(d) are plots of average magnetization of 3 NVs with $r=0.5$ nm. Hence, we
keep the TDVP as the main algorithm for our simulations.

(a)

(b)

(c)

(d)

Figure 3: Comparison of simulation results between $W^{II}$ and TDVP
algorithms for different $N$. (3(a)) Relative difference in real part of
$\langle\sigma_{z}\rangle$ from both algorithms. (3(b)) Numerical errors in
imaginary part of $W^{II}$ algorithm is 5 orders of magnitude larger than in
the TDVP algorithm. (3(c)) and (3(d)) Results from $N=3$. $W^{II}$ fails when
$r=0.5$ nm while TDVP can capture the dynamics compared to exact
diagonalization.

###  III.2 Dissipative dynamics with finite bond dimension

A major source of error when using MPS/MPDO is truncation error. Keeping
finite bond dimensions and truncating when the bond indices exceed
$\chi_{\text{max}}$ introduces an error. This error is given by the square
root of the sum of squares of the truncated singular values of all bonds,

$\epsilon=\sqrt{\sum_{k>\chi_{\text{max}}}(s_{k})^{2}}\,.$ (11)

Typically, these errors remain bounded if the state has entanglement that
follows an area-law, e.g. when the system has only nearest-neighbor
interactions. In such a case, the bond dimension grows with small singular
values during the time evolution, allowing the simulation of large systems
with small error. However, since we are considering a case of long-range
interaction, this argument should not hold. In Fig. 4 and Fig. 5 we plot
simulation results for different $\chi_{\text{max}}$ with $r=2.0$ nm and
$r=1.5$ nm, compared to a result from exact diagonalization for $N=4,7$. We
find that stronger interactions introduce bigger errors in bond truncation.
The results for smaller $\chi_{\text{max}}$ only follow the exact calculation
as long as entanglement entropy does still not reach $\chi_{\text{max}}$; i.e.
up to a certain number of timesteps before diverging and becoming inaccurate
due to truncation errors.

These results imply that it requires bigger $\chi_{\text{max}}$ to capture the
dynamics in the presence of long-range, strong interactions, thus less
efficiency of the MPS/MPDO method.

(a)

(b)

(c)

(d)

Figure 4: 4NVs with $\gamma=0$. The horizontal axis of each plots has a unit
of $dt=1$ ns.

(a)

(b)

(c)

(d)

Figure 5: 7NV with $\gamma=0$

Next, we add dissipation by setting $\gamma_{i}=\gamma\neq 0$. So each spin
has the same dissipation rate. Numerical values of $\gamma$ are given in
$(\mu\text{s})^{-1}$ units throughout. For the dissipation operator
$\hat{L}_{i}$ we choose the dephsasing operator, $\hat{L}_{i}=\hat{S}_{z}$,
relevant, e.g., for magnetic field noise. Fig. 6 shows the absolute errors
compared to exact calculations of dissipative dynamics for different
$\chi_{\text{max}}$. As a result, the discrepency of these dissipative
dynamics compared to the exact one is smaller than in the non-dissipative case
for the same $\chi_{\text{max}}$. When the growth of singular values of bi-
partition of the total state across the bond is limited by a larger $\gamma$,
bond truncation becomes more effective since it introduces less errors. This
demonstrates the interplay between the growth of entanglement entropy due to
interactions and dissipation that hinders it.

(a)

(b)

(c)

(d)

Figure 6: 4NV with $\gamma\neq 0$, compared with exact diagonalization.

To understand quantitatively how stronger interactions induce entanglement
entropy during time-evolution, we calculate the Operator Entanglement Entropy
(opEE) of the reduced density matrix determined by cutting the system into two
halves at the middle bond and tracing out the other half. We use the von
Neuman entropy defined as [33]

$S_{op}=-\sum_{i}^{\chi}(\lambda_{i})^{2}\log_{2}(\lambda_{i})^{2},$ (12)

where $\lambda_{i}$ are singular values of the cut bond with bond dimension
$\chi$ for the vectorized operator. Here the $\lambda_{i}$ are squared because
they are singular values of a vectorized operator $|\rho\rangle$. Note that in
general opEE does not give a direct measure of entanglement for mixed states.
The opEE for a pure state is twice the value of its standard entanglement
entropy, $S_{op}(|\rho\rangle)=2S(|\psi\rangle)$, when
$\rho=|\psi\rangle\langle\psi|$ [34]. Yet, the opEE provides insight into how
well the state can be approximated by an MPDO, indicating simulation errors
from a truncation [35].

After time-evolution, the middle bond index $\alpha_{(N/2+1)}$ which is the
bond that equally separates the MPDO into two halves has the largest
dimension. Fig. 7(d) shows the time evolution of opEE calculated from cutting
the middle bond. The opEE growth as function of time is accelerated as
interactions become stronger. This agrees with results of larger errors from
truncation shown earlier. Except at very strong interaction, i.e. for $r=1.5$
nm, opEE grows at the beginning and saturates rapidly. This is because the
interaction strengths $|C_{\text{dip}}|$ are comparable to the Rabi frequency
$\Omega$. However, we observe that this can be altered by increasing $\Omega$.

(a)

(b)

(c)

(d)

Figure 7: Operator entanglement entropy of middle bond. (7(a)) and (7(b)) show
trends of opEE to grow faster when interactions get stronger. (7(c)) and
(7(d))

dissipation suppresses generation of opEE.

###  III.3 Dynamics of QFI

To quantify the sensitivity of the ensemble of interacting NVs to a uniform
magnetic field $B_{z}$, we compute the Quantum Fisher Information of the mixed
state during time evolution. We use the tensor network approach for
calculating QFI as in [36]. Note that, unlike in the original work [36], we
only optimize the symmetric logarithmic derivative (SLD) but not the input
state $\rho$. In short, we iteratively search for the SLD that satisfies a
definition of QFI:

$F(\rho_{\phi},L)=\underset{L}{\text{sup}}[2\text{Tr}(\rho^{\prime}_{\phi}L)-\text{Tr}(\rho_{\phi}L^{2})].$
(13)

At the beginning, we randomly create the MPO approximation of an SLD, $L$,
given by

$L=\sum_{jk}\text{Tr}(S[1]_{k_{1}}^{j_{1}}...S[n]_{k_{n}}^{j_{n}})|\textbf{j}\rangle\langle\textbf{k}|.$
(14)

Here $S[l]_{k_{l}}^{j_{l}}$ is a Hermitian matrix. In searching for an optimal
$L$ we locally update $(S[1]_{k_{1}}^{j_{1}},...,S[n]_{k_{n}}^{j_{n}})$ by
sweeping from $S[1]$ to $S[n]$ and back to $S[1]$ again until $F$ has
converged. For example, when updating $S[l]_{k_{l}}^{j_{l}}$, all other
tensors $S_{k}^{j}$ are fixed and combined. Then after contractions, Eq. 13
becomes

$F(\rho,L)=2\sum_{\alpha}b_{\alpha}S[l]_{\alpha}-2\sum_{\alpha\beta}S[l]_{\alpha}A_{\alpha\beta}S[l]_{\beta}.$
(15)

Here, $b$ and $A$ are a vector and a matrix resulting from contracting the
fixed tensors and combining the un-contracted indices. The diagrammatic
explanation can be found on page 8 of [36]. Taking the derivative with respect
to $S[l]$, $\partial F(\rho,L)/\partial S[L]=0$, Eq. 15 yields
$\frac{1}{2}(A+A^{T})|S[l]\rangle=|b\rangle$. A solution to this equation
provides a local extremum for the QFI. Since this local update approach tends
to get stuck at local extrema, several repetitions with different initial $L$
are needed and can be computational expensive for large $N$.

Fig. 8 shows the dynamics of the QFI for a non-dissipative system with small
size and different $r$. At each time step, we find the optimal SLD from 10
independent realizations and select the one with the largest QFI. Furthermore,
we utilize the optimal SLD obtained from the previous time step as the
initialization for the current time step. Based on our experience, although
not guaranteed, implementing such a strategy can lead to a smoother result for
the QFI.

According to the plots, the QFI of $r=4.0$ nm nearly reaches values of $N$
non-interacting probes, denoted as $F_{|+1\rangle^{\otimes N}}=N$. The QFI
increases if the interactions are stonger, especially when $r\leq 2.5$ nm,
leading to a QFI that can be greater than $N$. This suggests that the
sensitivity for the probes is enhanced by the entanglement that is created by
interaction. Still, as for opEE, the QFI remains relatively small when $r=1.5$
nm, i.e. for very strong interaction. We find that this is due to an
interaction strength of the nearest NV pairs ($C_{dip,(i,i+1)}=2\pi\times
15.41$ MHz) that becomes substantially larger than the Rabi frequency
($\Omega=2\pi\times 2.00$ MHz). As shown Fig. 9, the QFI can be further
increased also for $r=1.5$ nm by increasing the Rabi frequency.

(a)

(b)

(c)

Figure 8: Dynamics of QFI with Rabi frequency = $\Omega$. All plots are the
best value computed from 10 different SLD initializations. The insets contain
moving averages over 10 time steps.

(a)

(b)

(c)

Figure 9: Dynamics of QFI when the Rabi frequency is doubled from Fig. 8 to be
$2\Omega=(2\pi)4.00$ MHz. The QFI when $r=1.5$ nm improvements can be
observed.

##  IV Conclusion

An ensemble of NV centers is a promising quantum metrology platform. While
adding more independent spins improves the sensitivity, the probe can in
principle also reach greater sensitivity due to entanglement created by
interactions. Also, having a dense ensemble means that the total size of the
probe can be smaller, resulting in better spatial resolution. However, without
external control of the entanglement and the states generated, strong dipole-
dipole interactions due to small separations between NV centers within a dense
ensemble are usually detrimental for the sensitivity as they lead to a rapid
decay of the total spin.

To investigate a strongly interacting NVs system, we simulated the master
equation with dissipation using MPDO. We benchmarked tensor network algorithms
for time evolution with a long-range interaction model. We found the $W^{II}$
algorithm to have bigger imaginary errors and to be less stable in a very
small separation regime than TDVP. Then, for small NV separations, we
investigated the effect of finite maximum bond dimensions on the simulation
accuracy. We determined the accuracy by comparing the obtained states to those
from exact diagonalization. Stronger interactions result in larger errors,
suggesting the need for larger bonds. For dissipative dynamics, we found
higher accuracy, which shows suppression of operator entanglement entropy
growth. These results indicate that the tensor network method can be suitable
for simulating an open system.

For applications in quantum sensing, we used the approach of locally
optimizing approximations of the SLD to obtain the quantum Fisher information.
The QFI of driven ensembles shows that strong interaction can create
entanglement-enhanced sensitivity compared to independent spins.

Finally, we observed that the boost in sensitivity can diminish when the
interaction becomes significantly larger than the Rabi frequency. In this
situation, it is necessary to utilize microwaves with higher frequencies for
driving in order to increase the sensitivity.

##  Appendix A Pure state entanglement entropy for a vector and an operator.

Any pure state has a Schmidt decomposition
$|\psi_{AB}\rangle=\sum_{i}^{r}\sqrt{\lambda}_{i}|i\rangle_{A}\otimes|i\rangle_{B}$
with Schmidt coefficients $\sqrt{\lambda}_{i}>0$ and
$\sum_{i}^{r}\lambda_{i}=1$. Then we can define the reduced density matrix for
a bipartite system by partial trace over the other parts:
$\rho_{A}=\Tr_{B}[\rho_{AB}]=\sum_{i}^{r}\lambda_{i}|i\rangle_{A}\langle
i|_{A}$ and
$\rho_{B}=\Tr_{A}[\rho_{AB}]=\sum_{j}^{r}\lambda_{j}|j\rangle_{B}\langle
j|_{B}.$ The density operator of the total system reads
$\rho_{AB}=|\psi_{AB}\rangle\langle\psi_{AB}|=\sum_{ij}^{r,r}\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}(|i\rangle_{A}\otimes|i\rangle_{B})(\langle
j|_{A}\otimes\langle j|_{B})$. The entanglement entropy for reduced density
operators is defined as

$\displaystyle
S(\rho_{A})=-\Tr[\rho_{A}\log(\rho_{A})]=-\sum_{i=1}^{r}\lambda_{i}\log{\lambda_{i}},$
(16)

where $S(\rho_{A})=S(\rho_{B})$. To find the operator entanglement entropy
(opEE), we vectorize $\rho_{AB}\rightarrow|\rho_{AB}\rangle\rangle$,

$\displaystyle|\rho_{AB}\rangle\rangle$
$\displaystyle=|\psi_{AB}\rangle|\psi_{AB}\rangle$ (17)
$\displaystyle=\sum_{\mu}^{r^{2}}\Lambda_{\mu}|i_{A}i_{B};j_{A}j_{B}\rangle.$
(18)

Here $\Lambda_{\mu}=\sqrt{\lambda_{i}}\sqrt{\lambda_{j}}$ and $\mu=(i,j)$. The
super density operator, denoted by $\rho^{\sharp}$, can be created using an
outer product:

$\displaystyle\rho^{\sharp}_{AB}$
$\displaystyle=|\rho_{AB}\rangle\langle\rho_{AB}|$ (19)
$\displaystyle=\sum_{\begin{subarray}{c}\mu=(i,j),\\\
\nu=(k,l)\end{subarray}}^{r^{2},r^{2}}\Lambda_{\mu}\Lambda_{\nu}|i_{A}i_{B};j_{A}j_{B}\rangle\langle
k_{A}k_{B};l_{A}l_{B}|.$ (20)

Similar to $\rho_{A}(\rho_{B})$, we can also have
$\rho_{A}^{\sharp}(\rho_{B}^{\sharp})$; e.g.
$\rho_{A}^{\sharp}=\sum_{\mu}^{r^{2}}\Lambda_{\mu}^{2}|i_{A};j_{A}\rangle\langle
i_{A};j_{A}|,$ Then we can find the opEE,

$\displaystyle S_{OP}(\rho_{A}^{\sharp})$
$\displaystyle=-\left[\sum_{\mu}\Lambda_{\mu}^{2}\log(\Lambda_{\mu}^{2})\right]$
(21)
$\displaystyle=-\left[\sum_{i,j}(\lambda_{i}\lambda_{j})\log(\lambda_{i}\lambda_{j})\right]$
(22)
$\displaystyle=-\left[\sum_{i,j}(\lambda_{i}\lambda_{j})\left(\log(\lambda_{i})+\log(\lambda_{j})\right)\right]$
(23)
$\displaystyle=-\left[\sum_{i}(\lambda_{i})\log(\lambda_{i})\sum_{j}(\lambda_{j})+\sum_{i}(\lambda_{i})\sum_{j}(\lambda_{j})\log(\lambda_{j})\right]$
(24) $\displaystyle=-\left[2\sum_{i}\lambda_{i}\log(\lambda_{i})\right]$ (25)
$\displaystyle=2S(\rho_{A}).$ (26)

##  Appendix B opEE for Rabi frequency $2\Omega$

The operator entanglement entropy when the Rabi frequency is doubled from
$\Omega=(2\pi)2.00$ MHz, hence $2\Omega=2\pi\times 4.00$ MHz.

(a)

(b)

Figure 10: Operator entanglement entropy of middle bond. Rabi frequency =
$2\Omega$.

###### Acknowledgements.

This work was supported by the Baden Württemberg Stiftung, project CDINQUA.
Tensor network codes are modified and implemented algorithms from TenPy [37].

## References

  * Taylor  _et al._ [2008] J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R. Hemmer, A. Yacoby, R. Walsworth, and M. D. Lukin, Nature Phys 4, 810 (2008), publisher: Nature Publishing Group.
  * Welter  _et al._ [2022] P. Welter, J. Rhensius, A. Morales, M. S. Wörnle, C. H. Lambert, G. Puebla-Hellmann, P. Gambardella, and C. L. Degen, Applied Physics Letters 120, 074003 (2022).
  * Dwyer  _et al._ [2022] B. L. Dwyer, L. V. Rodgers, E. K. Urbach, D. Bluvstein, S. Sangtawesin, H. Zhou, Y. Nassab, M. Fitzpatrick, Z. Yuan, K. De Greve, E. L. Peterson, H. Knowles, T. Sumarac, J.-P. Chou, A. Gali, V. Dobrovitski, M. D. Lukin, and N. P. de Leon, PRX Quantum 3, 040328 (2022), publisher: American Physical Society.
  * Haruyama  _et al._ [2019] M. Haruyama, S. Onoda, T. Higuchi, W. Kada, A. Chiba, Y. Hirano, T. Teraji, R. Igarashi, S. Kawai, H. Kawarada, Y. Ishii, R. Fukuda, T. Tanii, J. Isoya, T. Ohshima, and O. Hanaizumi, Nature Communications 2019 10:1 10, 1 (2019), publisher: Nature Publishing Group.
  * Balasubramanian  _et al._ [2019] P. Balasubramanian, C. Osterkamp, Y. Chen, X. Chen, T. Teraji, E. Wu, B. Naydenov, and F. Jelezko, Nano Lett. 19, 6681 (2019), publisher: American Chemical Society.
  * Zheng  _et al._ [2019] H. Zheng, J. Xu, G. Z. Iwata, T. Lenz, J. Michl, B. Yavkin, K. Nakamura, H. Sumiya, T. Ohshima, J. Isoya, J. Wrachtrup, A. Wickenbrock, and D. Budker, Phys. Rev. Appl. 11, 064068 (2019), publisher: American Physical Society.
  * Wang  _et al._ [2022] C. Wang, Q. Liu, Y. Hu, F. Xie, K. Krishna, N. Wang, L. Wang, Y. Wang, K. C. Toussaint Jr, J. Cheng, H. Chen, and Z. Wu, Realization of high-dynamic-range broadband magnetic-field sensing with ensemble nitrogen-vacancy centers in diamond (2022), arXiv:2209.11360 [quant-ph].
  * Zhou  _et al._ [2020] H. Zhou, J. Choi, S. Choi, R. Landig, A. M. Douglas, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, P. Cappellaro, H. S. Knowles, H. Park, and M. D. Lukin, Phys. Rev. X 10, 031003 (2020), publisher: American Physical Society.
  * Bauch  _et al._ [2020] E. Bauch, S. Singh, J. Lee, C. A. Hart, J. M. Schloss, M. J. Turner, J. F. Barry, L. M. Pham, N. Bar-Gill, S. F. Yelin, and R. L. Walsworth, PHYSICAL REVIEW B 102, 134210 (2020).
  * Farfurnik  _et al._ [2018] D. Farfurnik, Y. Horowicz, and N. Bar-Gill, Phys. Rev. A 98, 033409 (2018).
  * Osterkamp  _et al._ [2019] C. Osterkamp, M. Mangold, J. Lang, P. Balasubramanian, T. Teraji, B. Naydenov, and F. Jelezko, Sci Rep 9, 5786 (2019), publisher: Nature Publishing Group.
  * Sandvik [2010] A. W. Sandvik, Phys. Rev. Lett. 104, 137204 (2010).
  * Schiffer  _et al._ [2019] S. Schiffer, J. Wang, X.-J. Liu, and H. Hu, Phys. Rev. A 100, 063619 (2019), publisher: American Physical Society.
  * Cheng [2023] C. Cheng, Phys. Rev. B 108, 155113 (2023).
  * Essink [2018] S. Essink,  _Boundary-Driven XXZ Spin-1/2 Chain_ , Ph.D. thesis, Universität Bonn, Bonn (2018).
  * Kucsko  _et al._ [2018] G. Kucsko, S. Choi, J. Choi, P. Maurer, H. Zhou, R. Landig, H. Sumiya, S. Onoda, J. Isoya, F. Jelezko, E. Demler, N. Yao, and M. Lukin, Phys. Rev. Lett. 121, 023601 (2018).
  * Vidal [2003] G. Vidal, Phys. Rev. Lett. 91, 147902 (2003), publisher: American Physical Society.
  * Perez-Garcia  _et al._ [2007] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Quantum Info. Comput. 7, 401–430 (2007).
  * Schollwöck [2011] U. Schollwöck, Annals of Physics January 2011 Special Issue, 326, 96 (2011).
  * Orús [2014] R. Orús, Annals of Physics 349, 117 (2014).
  * Jarzyna and Demkowicz-Dobrzański [2013] M. Jarzyna and R. Demkowicz-Dobrzański, Phys. Rev. Lett. 110, 240405 (2013).
  * Finsterhölzl  _et al._ [2020] R. Finsterhölzl, M. Katzer, A. Knorr, and A. Carmele, Entropy 22, 984 (2020), number: 9 Publisher: Multidisciplinary Digital Publishing Institute.
  * Schlimgen  _et al._ [2022] A. W. Schlimgen, K. Head-Marsden, L. M. Sager, P. Narang, and D. A. Mazziotti, Phys. Rev. Research 4, 023216 (2022).
  * Fux  _et al._ [2023] G. E. Fux, D. Kilda, B. W. Lovett, and J. Keeling, Phys. Rev. Research 5, 033078 (2023).
  * Paris [2009] M. G. A. Paris, Int. J. Quantum Inform. 07, 125 (2009), publisher: World Scientific Publishing Co.
  * Giovannetti  _et al._ [2004] V. Giovannetti, S. Lloyd, and L. Maccone, Science 306, 1330 (2004), publisher: American Association for the Advancement of Science.
  * Weimer  _et al._ [2021] H. Weimer, A. Kshetrimayum, and R. Orús, Rev. Mod. Phys. 93, 015008 (2021).
  * Jaschke  _et al._ [2018] D. Jaschke, S. Montangero, and L. D. Carr, Quantum Science and Technology 4, 013001 (2018).
  * Paeckel  _et al._ [2019] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Annals of Physics 411, 167998 (2019).
  * Zaletel  _et al._ [2015] M. P. Zaletel, R. S. K. Mong, C. Karrasch, J. E. Moore, and F. Pollmann, Phys. Rev. B 91, 165112 (2015).
  * Haegeman  _et al._ [2016] J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016), publisher: American Physical Society.
  * Hubig  _et al._ [2018] C. Hubig, J. Haegeman, and U. Schollwöck, Phys. Rev. B 97, 045125 (2018), publisher: American Physical Society.
  * Murciano  _et al._ [2024] S. Murciano, J. Dubail, and P. Calabrese, J. Phys. A: Math. Theor. 57, 145002 (2024), publisher: IOP Publishing.
  * Preisser Beltrán [2023] G. J. Preisser Beltrán,  _Numerical studies of entanglement dynamics in open spin systems_ , These de doctorat, Strasbourg (2023).
  * Preisser  _et al._ [2023] G. Preisser, D. Wellnitz, T. Botzung, and J. Schachenmayer, Phys. Rev. A 108, 012616 (2023).
  * Chabuda  _et al._ [2020] K. Chabuda, J. Dziarmaga, T. J. Osborne, and R. Demkowicz-Dobrzański, Nat Commun 11, 250 (2020), number: 1 Publisher: Nature Publishing Group.
  * Hauschild and Pollmann [2018] J. Hauschild and F. Pollmann, SciPost Phys. Lect. Notes , 5 (2018).