# Tensor-Network Approaches to Counting Statistics for the Current in a
Boundary-Driven Diffusive System

Jiayin Gu<EMAIL_ADDRESS>Fan Zhang<EMAIL_ADDRESS>
###### Abstract

We apply tensor networks to counting statistics for the stochastic particle
transport in an out-of-equilibrium diffusive system. This system is composed
of a one-dimensional channel in contact with two particle reservoirs at the
ends. Two tensor-network algorithms, namely, Density Matrix Renormalization
Group (DMRG) and Time Evolving Block Decimation (TEBD), are respectively
implemented. The cumulant generating function for the current is numerically
calculated and then compared with the analytical solution. Excellent agreement
is found, manifesting the validity of these approaches in such an application.
Moreover, the fluctuation theorem for the current is shown to hold.

## 1 Introduction

Away from equilibrium, open systems composed of atoms and molecules manifest
macroscopic currents dissipating energy and producing thermodynamic entropy
[Prigogine_1967, Nicolis_RepProgPhys_1979, deGroot_1984, Callen_1985]. These
currents are intrinsically fluctuating with the probabilistic nature fully
characterized by their large deviation properties [Ellis_2006, Vulpiani_2014,
Touchette_PhysRep_2009]. The theoretical progress in the last three decades
reveals that, for nonequilibrium systems in steady state, there exists a
strikingly simple and general relation for the large deviation properties of
the fluctuating currents. This relation is nowadays known as the fluctuation
theorem [Evans_PhysRevLett_1993, Gallavotti_PhysRevLett_1996,
Kurchan_JPhysA_1998, Lebowitz_JStatPhys_1999, Gaspard_JChemPhys_2004,
Andrieux_JChemPhys_2004, Andrieux_JStatMech_2006, Derrida_JStatMech_2007]. At
the mesoscopic level of description, the evolution of the coarse-grained state
of the concerned system can be modeled as a Markov jump process. Accordingly,
the probability distribution for the system states obeys a master equation. To
evaluate the large deviation function for the current fluctuations, a direct
method is to numerically simulate the Markov jump process with the Gillespie
algorithm [Gillespie_JComputPhys_1976, Gillespie_JPhysChem_1977], which
generates a sample of random trajectories. Then, we perform counting
statistics for the integrated current, and the probability distribution of the
current can be finally constructed from the statistics [Touchette_arXiv_2012].
However, such direct sampling method encounters problem as it can not
accurately estimate probability of rare events or large deviations given
limited number of data. An alternative method is to first calculate the
cumulant generating function for the currents, and then obtain the large
deviation function through the Legendre-Fenchel transform. The cumulant
generating function turns out to be the leading eigenvalue of the tilted
generator of master equation including counting parameters
[Touchette_arXiv_2012]. Thus, it has transformed into an eigenvalue problem.

For a spatially extended system, a coarse-grained description at the
mesoscopic level requires partitioning the system in space. This leads to the
characterization of the states with many random variables. Accordingly, the
dimension of the tilted generator is very large, growing exponentially with
the number of these variables. The same issue arises in quantum many-body
systems for their underlying Hilbert space. Over the past few decades, tensor
networks [Orus_NatRevPhys_2019, Orus_AnnPhys_2014, Verstraete_AdvPhys_2008,
Schollwock_RevModPhys_2005, Schollwock_AnnPhys_2011, Cirac_RevModPhys_2021,
Montangero_2018, Ran_2020] have emerged as one of the most powerful tools for
handling such numerical complexity. Although originally developed in the
context of condensed matter physics, tensor networks are presently finding
applications in stochastic dynamics [Hieida_JPhysSocJpn_1998,
Carlon_EurPhysJB_1999, Degenhard_MultiscaleModelSimul_2004,
Helms_PhysRevE_2019, Helms_PhysRevLett_2020, Nagy_JStatPhys_2002,
Temme_PhysRevLett_2010, Johnson_PhysRevE_2010, Gorissen_PhysRevE_2009,
Gorissen_PhysRevE_2012, Gorissen_JPhysA_2011, Johnson_PhysRevLett_2015,
Banuls_PhysRevLett_2019, Causer_arXiv_2022, Strand_arXiv_2022a,
Strand_arXiv_2022b]. Indeed, there is a close analogy between the classical
stochastic systems in nonequilibrium steady state and the quantum many-body
systems in groud states. In this regard, two tensor-network approaches stand
out, namely, Density Matrix Renormalization Group (DMRG)
[White_PhysRevLett_1992, White_PhysRevB_1993, Schollwock_RevModPhys_2005,
Schollwock_AnnPhys_2011] and Time Evolving Block Decimation (TEBD)
[Vidal_PhysRevLett_2003, Vidal_PhysRevLett_2004]. For quantum many-body
systems, the former is a variational approach that enables the systematic
search for ground states, while the latter an approach that evolves the
states. If implemented in imaginary time, the latter can also be used to find
the ground states. Correspondingly, for classical stochastic systems, these
two approaches are naturally well suited for obtaining the nonequilibrium
steady states.

Here, our purpose is to provide a pedagogical demonstration of how DMRG and
TEBD are used to calculate the cumulant generating function for the current in
a boundary-driven diffusive system. The rest of this paper is structured as
follows. The stochastic description for the diffusive process is introduced in
Section 2 where we discretize the system in space and then establish the
master equation. The counting parameter for the particle transitions is
included in the master equation, and consequently the cumulant generating
function for the current is analytically obtained as the leading eigenvalue.
Section 3 is devoted to the DMRG approach to counting statistics. The
probability distribution is represented as a matrix product state (MPS). The
generator for the master equation is formulated in terms of the local
creation, annihilation, and particle number operators, and is then constructed
as a matrix product operator (MPO). In Section 4, the TEBD approach to
counting statistics is concerned. The small-time evolution operator for the
tilted master equation is approximately decomposed into a ordered product of
local ones. The numerical details and results are presented in Section 5,
followed by the concluding remarks given in Section 6. The cumulant generating
function is formally derived in Appendix A.

## 2 Boundary-Driven Diffusive System

Figure 1: Schematic representation of a transport channel in contact with two
reservoirs.

The system we consider here is modeled as a one-dimensional channel in contact
with two reservoirs at the ends, as illustrated in Figure 1. Distributed
inside the channel are mobile particles with their densities expressed as a
function of position, $n(x)$. At the boundaries, because the reservoirs are
supposed to be arbitrarily large, the densities of particles are constant in
time with $n_{\rm L}$ for the left reservoir and $n_{\rm R}$ for the right
reservoir. In the case that $n_{\rm L}\neq n_{\rm R}$, the system is driven
out of equilibrium and a directional current flows from one reservoir to the
other. The mobile particles undergo the Brownian motion [Einstein_1956,
Uhlenbeck_PhysRev_1930, Wang_RevModPhys_1945, Mazo_2002], so that the time
evolution of the densities is given by the diffusion equation

$\displaystyle\frac{\partial n(x,t)}{\partial
t}=D\frac{\partial^{2}n(x,t)}{\partial x^{2}}\text{,}$ (1)

where $D$ is the diffusion coefficient determined phenomenally by direct
observation of the random walk of a tagged Brownian particle according to the
formula

$\displaystyle
D=\lim_{t\to\infty}\frac{1}{2t}\left\langle(x_{t}-x_{0})^{2}\right\rangle\text{.}$
(2)

Here, $x_{0}$ (respectively, $x_{t}$) is the position at the initial time
(respectively, time $t$), and $\langle\cdot\rangle$ denotes the average taken
over the trajectories.

### 2.1 Stochastic Description

At the mesoscopic level, the time evolution of the distribution of mobile
particles in the channel can be described as a Markov jump process, which is
formulated in terms of a master equation (see Ref. [Gaspard_NewJPhys_2005] and
references therein). The fluctuations down to the mesoscopic scale can be
fully characterized by such a stochastic description. For this purpose, the
channel is spatially discretized into $L$ cells, each of width $\Delta x$ and
containing some number of mobile particles. These cells are labeled with
indices $i=1,2,\cdots,L$. For notational convenience, the indices $i=0$ and
$i=L+1$ are respectively used to refer to the cells for left and right
reservoirs. The number of particles in both cells for reservoirs are
maintained constant in time, and are given by $N_{0}\equiv\bar{N}_{\rm L}$ and
$N_{L+1}\equiv\bar{N}_{\rm R}$. The system state is specified by the numbers
${\bf N}=\\{N_{i}\\}_{i=1}^{L}$ of particles in the cells. These numbers
change in time due to the random transitions according to the kinetic network

$\leavevmode\resizebox{469.75499pt}{}{$\begin{array}[]{ccccccccccccc}\bar{N}_{\rm
L}&\autorightleftharpoons{$\scriptstyle W_{0}^{(+)}$}{$\scriptstyle
W_{0}^{(-)}$}&N_{1}&\autorightleftharpoons{$\scriptstyle
W_{1}^{(+)}$}{$\scriptstyle
W_{1}^{(-)}$}&N_{2}&\autorightleftharpoons{$\scriptstyle
W_{2}^{(+)}$}{$\scriptstyle
W_{2}^{(-)}$}&\cdots&\autorightleftharpoons{$\scriptstyle
W_{L-2}^{(+)}$}{$\scriptstyle
W_{L-2}^{(-)}$}&N_{L-1}&\autorightleftharpoons{$\scriptstyle
W_{L-1}^{(+)}$}{$\scriptstyle
W_{L-1}^{(-)}$}&N_{L}&\autorightleftharpoons{$\scriptstyle
W_{L}^{(+)}$}{$\scriptstyle W_{L}^{(-)}$}&\bar{N}_{\rm
R}\text{,}\end{array}$}$

with the transition rates linearly depending on the particle number in the
departure cell, $W_{i}^{(+)}=kN_{i}$, $W_{i}^{(-)}=kN_{i+1}$. Here, $k$ is the
rate constant related to $D$ by $k=D/\Delta x^{2}$ which can be obtained
through the continuum limit111The mean current is expressed as  $\displaystyle
J=-k\langle N_{i+1}\rangle+k\langle N_{i}\rangle=-k\Omega\Delta
x\frac{\partial n}{\partial x}=j\Sigma=-D\frac{\Omega}{\Delta x}\frac{\partial
n}{\partial x}\text{,}$ where $\Omega$ and $\Sigma$ denote the volume and
section area of a cell, respectively. From this expression, $k=D/\Delta x^{2}$
is immediately arrived.. The probability ${\cal P}({\bf N},t)$ that the cells
contain the particle numbers ${\bf N}$ for time $t$ obeys the master equation

$\displaystyle\frac{{\rm d}{\cal P}}{{\rm d}t}=\hat{L}{\cal P}=$
$\displaystyle k\bar{N}_{\rm L}{\cal P}(N_{1}-1,\cdots,t)-k\bar{N}_{\rm
L}{\cal P}(N_{1},\cdots,t)$ $\displaystyle+k(N_{1}+1){\cal
P}(N_{1}+1,\cdots,t)-kN_{1}{\cal P}(N_{1},\cdots,t)$ $\displaystyle\hskip
28.45274pt\vdots$ $\displaystyle+k(N_{i}+1){\cal
P}(\cdots,N_{i}+1,N_{i+1}-1,\cdots,t)-kN_{i}{\cal
P}(\cdots,N_{i},N_{i+1},\cdots,t)$ $\displaystyle+k(N_{i+1}+1){\cal
P}(\cdots,N_{i}-1,N_{i+1}+1,\cdots,t)-kN_{i+1}{\cal
P}(\cdots,N_{i},N_{i+1},\cdots,t)$ $\displaystyle\hskip 28.45274pt\vdots$
$\displaystyle+k(N_{L}+1){\cal P}(\cdots,N_{L}+1,t)-kN_{L}{\cal
P}(\cdots,N_{L},t)$ $\displaystyle+k\bar{N}_{\rm R}{\cal
P}(\cdots,N_{L}-1,t)-k\bar{N}_{\rm R}{\cal P}(\cdots,N_{L},t)\text{,}$ (3)

where $\hat{L}$ denotes the generator for the time evolution. The stationary
solution of the above master equation is given by a product of Poisson
distributions,

$\displaystyle{\cal P}_{\rm
st}(N_{1},\cdots,N_{L})=\prod_{i=1}^{L}\frac{\langle N_{i}\rangle_{\rm
st}^{N_{i}}}{N_{i}!}{\rm e}^{{-\langle N_{i}\rangle}_{\rm st}}\text{,}$ (4)

with the mean number of particles corresponding to a linear profile of
concentration,

$\displaystyle\langle N_{i}\rangle_{\rm st}=\bar{N}_{\rm L}+\frac{\bar{N}_{\rm
R}-\bar{N}_{\rm L}}{L+1}i\text{.}$ (5)

### 2.2 Counting Statistics

To investigate the current fluctuations for the particle transport, we now
perform counting statistics. For this purpose, the master equation (3) is
modified to include a counting parameter $\lambda$ for particle transitions
between the left-reservoir cell and the first one for the discretized channel,
yielding

$\displaystyle\frac{{\rm d}F}{{\rm d}t}=\hat{L}_{\lambda}F=$ $\displaystyle
k\bar{N}_{\rm L}{\rm e}^{-\lambda}F(N_{1}-1,\cdots,t)-k\bar{N}_{\rm
L}F(N_{1},\cdots,t)$ $\displaystyle+k(N_{1}+1){\rm
e}^{+\lambda}F(N_{1}+1,\cdots,t)-kN_{1}F(N_{1},\cdots,t)$ $\displaystyle\hskip
28.45274pt\vdots$
$\displaystyle+k(N_{i}+1)F(\cdots,N_{i}+1,N_{i+1}-1,\cdots,t)-kN_{i}F(\cdots,N_{i},N_{i+1},\cdots,t)$
$\displaystyle+k(N_{i+1}+1)F(\cdots,N_{i}-1,N_{i+1}+1,\cdots,t)-kN_{i+1}F(\cdots,N_{i},N_{i+1},\cdots,t)$
$\displaystyle\hskip 28.45274pt\vdots$
$\displaystyle+k(N_{L}+1)F(\cdots,N_{L}+1,t)-kN_{L}F(\cdots,N_{L},t)$
$\displaystyle+k\bar{N}_{\rm R}F(\cdots,N_{L}-1,t)-k\bar{N}_{\rm
R}F(\cdots,N_{L},t)\text{.}$ (6)

The eigensolution of this tilted master equation has the form

$\displaystyle F(N_{1},\cdots,N_{L},t)={\rm
e}^{-Q(\lambda)t}\prod_{i=1}^{L}\frac{C_{i}^{N_{i}}}{N_{i}!}\text{,}$ (7)

where $\\{C_{i}\\}$ are constants depending on $\lambda$, and $Q(\lambda)$ is
the cumulant generating function. Substituting this solution into the tilted
master equation and after some algebraic manipulations, we obtain

$\displaystyle Q(\lambda)=\frac{k\bar{N}_{\rm L}}{L+1}\left(1-{\rm
e}^{-\lambda}\right)+\frac{k\bar{N}_{\rm R}}{L+1}\left(1-{\rm
e}^{\lambda}\right)\text{.}$ (8)

If we define $A=\ln\left(\bar{N}_{\rm L}/\bar{N}_{\rm R}\right)$ as the
affinity [Donder_1936] for driving the system out of equilibrium, Then, the
cumulant generating function (8) satisfies the Gallavotti-Cohen symmetry
[Gallavotti_PhysRevLett_1995, Gallavotti_JStatPhys_1995,
Lebowitz_JStatPhys_1999]

$\displaystyle Q(\lambda)=Q(A-\lambda)\text{,}$ (9)

which is also known under the name fluctuation theorem. The statistical
cumulants can be obtained by taking successive derivatives with respect to the
counting parameter $\lambda$. For example, the mean current and its
diffusivity are given by

$\displaystyle J=\left.\frac{\partial
Q}{\partial\lambda}\right|_{\lambda=0}=\frac{k\left(\bar{N}_{\rm
L}-\bar{N}_{\rm R}\right)}{L+1}\text{,}$ (10) $\displaystyle
D=-\frac{1}{2}\left.\frac{\partial^{2}Q}{\partial\lambda^{2}}\right|_{\lambda=0}=\frac{k\left(\bar{N}_{\rm
L}+\bar{N}_{\rm R}\right)}{2(L+1)}\text{.}$ (11)

For nonequilibrium systems described by master equation, the affinity can be
systematically obtained through Schnakenberg graph analysis
[Schnakenberg_RevModPhys_1976]. More details about the counting statistics for
this diffusive process can be found in Ref. [Andrieux_JStatMech_2006]. A
formal derivation of the cumulant generating function (8) is provided in
Appendix A.

From the numerical perspective, the determination of cumulant generating
function requires to calculate the leading eigenvalue of the tilted generator
$\hat{L}_{\lambda}$ defined in Eq. (6). In principle, this can be done with
various numerical routines for diagonalizing a matrix. However, this is often
very difficult in practice. Since the dimension of the tilted generator grows
exponentially with the number of the discretized cells. This motivates the
need of tensor networks to tackle such a problem.

## 3 DMRG Approach to Counting Statistics

In this section, we present in detail the DMRG procedure to perform counting
statistics for the stochastic process in the preceding section. The key
differences from the standard usage for quantum lattice systems are
highlighted.

### 3.1 Notations

First, it is helpful to introduce the quantum-mechanical notation for
subsequent description convenience. We use $\ket{N_{i}}$ for the state that
the $i$-th cell contains $N_{i}$ particle and $\ket{\bf N}$ for the system
state specified by the numbers ${\bf N}=\\{N_{i}\\}_{i=1}^{L}$. Any
distribution $F({\bf N})$ is written as

$\displaystyle\ket{F({\bf N})}=\sum_{\bf N}F(\bf N)\ket{\bf N}\text{.}$ (12)

Accordingly, $\bra{N_{i}}$, $\bra{\bf N}$, and $\bra{F^{\dagger}(\bf N)}$ are
used to denote their adjoints. Following the conventions in quantum mechanics,
2-norm and orthogonality are assumed,

$\displaystyle\braket{N^{\prime}}{N}=\delta_{N^{\prime},N}\text{,}$ (13)
$\displaystyle\braket{{\bf N}^{\prime}}{{\bf N}}=\delta_{{\bf N}^{\prime},{\bf
N}}\text{,}$ (14) $\displaystyle\braket{F^{\dagger}({\bf N})}{F({\bf
N})}=\sum_{\bf N}F^{\dagger}({\bf N})F({\bf N})=1\text{.}$ (15)

Now, an issue arises that the probability distribution is of 1-norm,
$\sum_{\bf N}{\cal P}({\bf N})=1$. We can introduce the 2-norm probability
distribution,

$\displaystyle\tilde{\cal P}({\bf N})\equiv\frac{{\cal P}({\bf
N})}{\sqrt{\sum_{\bf N}{\cal P}({\bf N}){\cal P}({\bf N})}}\text{,}$ (16)

so that $\braket{\tilde{\cal P}^{\dagger}({\bf N})}{\tilde{\cal P}({\bf
N})}=1$. The 1-norm probability distribution can be recovered,

$\displaystyle{\cal P}({\bf N})=\frac{\tilde{\cal P}({\bf N})}{\sum_{\bf
N}\tilde{\cal P}({\bf N})}\text{.}$ (17)

There is a one-to-one correspondence between these two kinds of probability
distributions. In practical implementation in code, we always use the 2-norm
version, since many tensor-network techniques already developed for quantum-
mechanical systems can be utilized.

### 3.2 MPS and MPO

Figure 2: Graphical representation of (a) an MPS and (b) an MPO with the
counting parameter $\lambda$ included in the first tensor.

In quantum many-body systems, the presence of the celebrated area laws
[Eisert_RevModPhys_2010] for the entanglement entropy enables an efficient
numerical representation of the ground state as an MPS
[PerezGarcia_arXiv_2007] (see Figure 2a for illustration). Here, MPS can also
be used to represent the stationary solution for classical stochastic systems
in nonequilibrium steady states. Moreover, we notice that the feature of joint
product in solutions (4) and (7) imply that the bond dimensions of their MPS
representations are both one. This is very remarkable that the computer memory
used to store the stationary solution can be very low, scaling linearly with
the number of cells. The physical indices denote the numbers $\\{N_{i}\\}$.

Next, we show how to construct an MPO (see Figure 2b) for the tilted generator
$\hat{L}_{\lambda}$. This constitutes the most crucial step in applying the
DMRG approach to counting statistics for the stochastic process. We adopt the
formalism similar to the second quantization. The particle transitions between
two neighboring cells can be associated with an annihilation operator acting
on one cell and an creation operator acting the other. Following this
reasoning, we introduce the local annihilation operator $a_{i}^{-}$ and
creation operator $a_{i}^{+}$, which are respectively defined by

$\displaystyle\braket{N_{i}}{a_{i}^{-}}{N_{i}^{\prime}}=kN_{i}^{\prime}\delta_{N_{i},,N_{i}^{\prime}-1}\text{,}$
(18)
$\displaystyle\braket{N_{i}}{a_{i}^{+}}{N_{i}^{\prime}}=\delta_{N_{i},N_{i}^{\prime}+1}\text{.}$
(19)

It should be noted that the operator definitions by Eqs. (18)-(19) only apply
to the cells for the channel, $1\leq i\leq L$. For the cells corresponding to
the reservoirs, $i=0,L+1$, these operators should be defined differently,

$\displaystyle a_{0}^{-}=k\bar{N}_{\rm L}\text{,}$ (20) $\displaystyle
a_{0}^{+}=a_{L+1}^{+}=1\text{,}$ (21) $\displaystyle a_{L+1}^{-}=k\bar{N}_{\rm
R}\text{.}$ (22)

Because the particle numbers in the reservoir cells are kept constant. To
account for the probability loss in the tilted master equation, we should
still introduce another operator $a_{i}$ defined by

$\displaystyle\braket{N_{i}}{a_{i}}{N_{i}^{\prime}}=kN_{i}^{\prime}\delta_{N_{i},N_{i}^{\prime}}\text{,}$
(23)

which can be intuitively termed as particle number operator. This definition
applies to any cell. In this way, the tilted generator is expressed as

$\displaystyle\hat{L}_{\lambda}=a_{0}^{-}\otimes a_{1}^{+}{\rm
e}^{-\lambda}+a_{0}^{+}\otimes a_{1}^{-}{\rm
e}^{+\lambda}+\sum_{i=1}^{L}\left(a_{i}^{-}\otimes
a_{i+1}^{+}+a_{i}^{+}\otimes
a_{i+1}^{-}\right)-\sum_{i=0}^{L}\left(a_{i}+a_{i+1}\right)\text{,}$ (24)

where $\otimes$ denotes the tensor product and the counting parameter
$\lambda$ is included for particle transitions between the left-reservoir cell
and the first one in the discretized channel. This tilted generator is now
split into two parts

$\displaystyle\hat{L}_{\lambda}=\hat{L}_{\lambda}^{(1)}+\hat{L}^{(2)}\text{,}$
(25)

with each constructed as an MPO,

$\displaystyle\hat{L}_{\lambda}^{(1)}$ $\displaystyle=a_{0}^{-}\otimes
a_{1}^{+}{\rm e}^{-\lambda}+a_{0}^{+}\otimes a_{1}^{-}{\rm
e}^{+\lambda}+\sum_{i=1}^{L}\left(a_{i}^{-}\otimes
a_{i+1}^{+}+a_{i}^{+}\otimes a_{i+1}^{-}\right)$
$\displaystyle=\begin{pmatrix}a_{0}^{-}{\rm
e}^{-\lambda}&a_{0}^{+}&0&1\end{pmatrix}\otimes\begin{pmatrix}0&0&a_{1}^{+}&0\\\
0&0&a_{1}^{-}{\rm e}^{+\lambda}&0\\\ 0&0&1&0\\\
a_{1}^{-}&a_{1}^{+}&0&1\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}0&0&a_{i}^{+}&0\\\
0&0&a_{i}^{-}&0\\\ 0&0&1&0\\\
a_{i}^{-}&a_{i}^{+}&0&1\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}a_{L+1}^{+}\\\
a_{L+1}^{-}\\\ 1\\\ 0\end{pmatrix}$
$\displaystyle=\begin{pmatrix}a_{1}^{-}&a_{1}^{+}&a_{1}^{+}k\bar{N}_{\rm
L}{\rm e}^{-\lambda}+a_{1}^{-}{\rm
e}^{+\lambda}&1\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}0&0&a_{i}^{+}&0\\\
0&0&a_{i}^{-}&0\\\ 0&0&1&0\\\
a_{i}^{-}&a_{i}^{+}&0&1\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}a_{L}^{+}\\\
a_{L}^{-}\\\ 1\\\ a_{L}^{-}+a_{L}^{+}k\bar{N}_{\rm R}\end{pmatrix}$ (26)

and

$\displaystyle\hat{L}^{(2)}$
$\displaystyle=\sum_{i=0}^{L}\left(-a_{i}-a_{i+1}\right)$
$\displaystyle=\begin{pmatrix}1&-a_{0}&1\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}0&-a_{i}&0\\\
0&1&0\\\ 1&-a_{i}&1\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}-a_{L+1}\\\
1\\\ 0\end{pmatrix}$ $\displaystyle=\begin{pmatrix}1&-k\bar{N}_{\rm
L}-2a_{1}&1\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}0&-a_{i}&0\\\
0&1&0\\\ 1&-a_{i}&1\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}-a_{L}\\\
1\\\ -a_{L}-k\bar{N}_{\rm R}\end{pmatrix}\text{,}$ (27)

where the boundary tensors are explicitly expressed according to Eqs.
(20)-(22) and contracted with their neighboring ones. So, there are $L$
tensors in each MPO as expected. If we denotes

$\displaystyle\hat{L}_{\lambda}^{(1)}=l_{1,\lambda}^{(1)}\otimes\cdots\otimes
l_{i}^{(1)}\otimes\cdots\otimes l_{L}^{(1)}$ (28)

and

$\displaystyle\hat{L}^{(2)}=l_{1}^{(2)}\otimes\cdots\otimes
l_{i}^{(2)}\otimes\cdots\otimes l_{L}^{(2)}\text{,}$ (29)

then the tilted generator is constructed as an MPO,

$\displaystyle\hat{L}_{\lambda}=l_{1,\lambda}\otimes\cdots\otimes
l_{i}\otimes\cdots\otimes
l_{L}=\begin{pmatrix}l_{1,\lambda}^{(1)}&l_{1}^{(2)}\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}l_{i}^{(1)}&0\\\
0&l_{i}^{(2)}\end{pmatrix}\otimes\cdots\otimes\begin{pmatrix}l_{L}^{(1)}\\\
l_{L}^{(2)}\end{pmatrix}\text{,}$ (30)

whose bond dimension is 7. It should be noticed that the tilted generator
$\hat{L}_{\lambda}$ is intrinsically non-symmetric or non-Hermitian, implying
that (i) the right eigenvectors are not identical to the corresponding left
eigenvectors; and (ii) the right (left) eigenvectors are not mutually
orthogonal. This is the basic difference from the case of quantum Hamiltonian.

### 3.3 Optimization

Figure 3: Graphical representation of the contraction between an normalized
MPS and its adjoint. Both are in the canonical form with the centers lying at
$F_{2}$ and $F_{2}^{\dagger}$, respectively.

The aim of DMRG approach is to variationally find the eigensolution from which
the leading eigenvalue can be computed. This can be formulated as an
optimization problem

$\displaystyle Q(\lambda)=-{\rm max}\left[\braket{F^{\dagger}({\bf
N})}{\hat{L}_{\lambda}}{F({\bf N})}\right]\text{,}$ (31)

under the condition $\braket{F^{\dagger}({\bf N})}{F({\bf N})}=1$. The
graphical representation of this normalization is shown in Figure 3 where the
canonical form of MPS is presented due to the gauge freedom. The advantage of
this canonical form lies in the efficient contraction of an MPS with its
adjoint; we only need to contract the canonical centers. The shift of the
canonical center can be realized with the singular value decomposition. The
problem (31) can be further reformulated into a constrained optimization
problem using a Lagrangian multiplier $-Q$,

$\displaystyle\braket{F^{\dagger}({\bf N})}{\hat{L}_{\lambda}}{F({\bf
N})}+Q\braket{F^{\dagger}({\bf N})}{F({\bf N})}\text{.}$ (32)

We now adopt the single-tensor optimization procedure, where only one tensor
$F_{i}^{\dagger}$ is varied at once while keeping all others fixed, so we have

$\displaystyle\frac{\partial}{\partial
F_{i}^{\dagger}}\left[\braket{F^{\dagger}({\bf N})}{\hat{L}_{\lambda}}{F({\bf
N})}+Q\braket{F^{\dagger}({\bf N})}{F({\bf N})}\right]=0\text{.}$ (33)

Figure 4 shows the graphical representation of the consequence after this
variational procedure. Because two MPSs are both in canonical form with the
centers at $F_{i}$ and $F_{i}^{\dagger}$, Eq. (33) simplifies to

$\displaystyle\tilde{\hat{L}}_{\lambda}\cdot F_{i}+QF_{i}=0\text{,}$ (34)

where $\tilde{\hat{L}}_{\lambda}$ the operator indicated as the red part in
Figure 5. This now represents a standard eigenvalue problem. The dimension of
the operator $\tilde{\hat{L}}_{\lambda}$ is greatly reduced, making it
possible to be handled with normal numerical routines.

Figure 4: The graphical representation of the variational procedure (33) with,
e.g., $F_{2}$ being the canonical center.

We now need to find the right eigenvector $F_{i}$ of
$\tilde{\hat{L}}_{\lambda}$ corresponding to the largest eigenvalue $-Q$. This
can be achieved through the iterative method. Since ${\rm
e}^{\tilde{\hat{L}}_{\lambda}}>0$, the Perron-Frobenius theorem applies, and
the leading eigenvalue $-Q$ of $\tilde{\hat{L}}_{\lambda}$ corresponds to the
real maximum eigenvalue ${\rm e}^{-Qt}$ of ${\rm
e}^{\tilde{\hat{L}}_{\lambda}t}$ in magnitude (for some $t>0$). Consequently,
the right eigenvector can be asymptotically evaluated as

$\displaystyle F_{i}\sim_{t\to\infty}{\rm
e}^{\tilde{\hat{L}}_{\lambda}t}F_{i}(0)\text{,}$ (35)

which is then normalized

$\displaystyle F_{i}\leftarrow\frac{F_{i}}{\sqrt{F^{\dagger}_{i}\cdot
F_{i}}}\text{.}$ (36)

Here, $F_{i}(0)$ is an arbitrary vector assumed to have the component of the
desired eigenvector. The matrix exponential can be computed with many
numerical routines [Moler_SIAMRev_1978]. Then, we have

$\displaystyle Q=-F_{i}^{\dagger}\cdot\tilde{\hat{L}}_{\lambda}\cdot
F_{i}\text{.}$ (37)

After completing the optimization on $F_{i}$, we move the canonical center to
a neighboring one and continue such procedure. This operation is performed
back and forth between the first and last tensors until the leading eigenvalue
converges. With different values of $\lambda$, we can numerically construct
the desired cumulant generating function $Q(\lambda)$.

Figure 5: Graphical representation of the eigenvalue problem (34).

## 4 TEBD Approach to Counting Statistics

In this section, we present another tensor-network approach to counting
statistics for the stochastic current. Now that we want to obtain the
eigensolution of the tilted generator $\hat{L}_{\lambda}$ corresponding to the
leading eigenvalue, can we achieve this with direct iteration? The answer to
this question is yes. This can be done with TEBD which is an algorithm relying
on Trotter-Suzuki decomposition [Suzuki_CommunMathPhys_1976]. The desired
eigensolution is asymptotically evaluated as

$\displaystyle\ket{F({\bf N})}\sim_{t\to\infty}{\rm
e}^{\hat{L}_{\lambda}t}\ket{F({\bf N},0)}\text{,}$ (38)

where $\ket{F({\bf N},0)}$ is an arbitrary initial distribution. At this
point, the dimension of the evolution operator is exponentially large. For
small-time step $\delta$, however, ${\rm e}^{\hat{L}_{\lambda}\delta}$ can be
decomposed into an ordered product of local ones of a much smaller dimension.
Because the particles hop between neighboring cells, the tilted generator
$\hat{L}_{\lambda}$ can be split into two parts,

$\displaystyle\hat{L}_{\lambda}=\hat{L}_{\rm odd}+\hat{L}_{\rm even}\text{,}$
(39)

with

$\displaystyle\hat{L}_{\rm odd}=l_{1,2,\lambda}+l_{3,4}+\cdots\text{,}$ (40)
$\displaystyle\hat{L}_{\rm even}=l_{2,3}+l_{4,5}+\cdots\text{,}$ (41)

where $l_{i,i+1}$ represents the local operator acting simultaneously on the
$i$-th and $(i+1)$-th cells. The local operators in each part can be
diagonalized efficiently and mutually commuting. It is noted here that the
counting parameter is included in $l_{1,2,\lambda}$. The exact evolution
operator can be decomposed to any order. Here, we give the second-order one,

$\displaystyle{\rm e}^{\delta\hat{L}_{\lambda}}\approx{\rm
e}^{\frac{\delta}{2}\hat{L}_{\rm odd}}{\rm e}^{\delta\hat{L}_{\rm even}}{\rm
e}^{\frac{\delta}{2}\hat{L}_{\rm odd}}+O(\delta^{3})\text{,}$ (42)

as it is commonly used. If the transport channel is discretized into $5$
cells, then the local operators are explicitly given by

$\displaystyle l_{1,2,\lambda}=k\bar{N}_{\rm L}{\rm
e}^{-\lambda}a_{1}^{+}\otimes{\sf I}_{2}+{\rm
e}^{+\lambda}a_{1}^{-}\otimes{\sf I}_{2}+a_{1}^{-}\otimes
a_{2}^{+}+a_{1}^{+}\otimes a_{2}^{-}-k\bar{N}_{\rm L}{\sf I}_{1}\otimes{\sf
I}_{2}-2a_{1}\otimes{\sf I}_{2}-{\sf I}_{1}\otimes a_{2}\text{,}$ (43)
$\displaystyle l_{3,4}=a_{3}^{-}\otimes a_{4}^{+}+a_{3}^{+}\otimes
a_{4}^{-}-a_{3}\otimes{\sf I}_{4}-{\sf I}_{3}\otimes a_{4}\text{,}$ (44)
$\displaystyle l_{2,3}=a_{2}^{-}\otimes a_{3}^{+}+a_{2}^{+}\otimes
a_{3}^{-}-a_{2}\otimes{\sf I}_{3}-{\sf I}_{2}\otimes a_{3}\text{,}$ (45)
$\displaystyle l_{4,5}=a_{4}^{-}\otimes a_{5}^{+}+a_{4}^{+}\otimes
a_{5}^{-}+{\sf I}_{4}\otimes a_{5}^{-}+k\bar{N}_{\rm R}{\sf I}_{4}\otimes
a_{5}^{+}-a_{4}\otimes{\sf I}_{5}-2{\sf I}_{4}\otimes a_{5}-k\bar{N}_{\rm
R}{\sf I}_{4}\otimes{\sf I}_{5}\text{,}$ (46)

where ${\sf I}_{i}$ stands for the identity operator acting on the $i$-th
cell. The local evolution operators are alternatively applied to the MPS (see
Figure 6 for illustration). Once we obtain the desired eigensolution
$\ket{F({\bf N})}$, the cumulant generating function can be computed,

$\displaystyle Q(\lambda)=-\frac{1}{\delta}\ln\braket{F^{\dagger}({\bf
N})}{{\rm e}^{\hat{L}_{\lambda}\delta}}{F({\bf N})}\text{,}$ (47)

with $\delta$ being a very small real number.

Figure 6: Graphical representation of the TEBD algorithm implemented with the
second-order Trotter-Suzuki decomposition. The local evolution operators are
alternatively applied to the MPS. The counting parameter $\lambda$ is included
in $l_{1,2,\lambda}$.

## 5 Numerical Results

Table 1: The values of quantities and parameters specified in numerical computations. diffusion coefficient, $D$ |  $0.01$
---|---
width of each cell, $\Delta x$ |  $0.1$
number of cells in discretized channel, $L$ |  $5$
number of particle in left-reservoir cell, $\bar{N}_{\rm L}$ |  $9$
number of particle in right-reservoir cell, $\bar{N}_{\rm R}$ |  $3$
dimension of $a_{i}^{+}$, $a_{i}^{-}$, and $a_{i}$ after truncation, $N_{\rm max}$ |  $20$

Figure 7: Mean particle number $\langle N_{i}\rangle$. Dash line plots the
exact linear profile solution (5). Asterisks are the numerical results from
DMRG approach, while the circles from TEBD approach. The values for numerical
computations are listed in Table 1.

With the tensor-network approaches having already exhibited in the preceding
sections, we now dive into the numerical details and present the results.
First of all, we define the system with specific numbers which are listed in
Table 1. We realize that an arbitrary number of particles can concentrate in a
discretized cell of the transport channel. However, the probability of the
event that the particle numbers $\\{N_{i}\\}$ far exceeds their mean value is
exponentially small. This means that we can truncate the state space. The
particle numbers in a cell only take the values $N_{i}=0,1,2,\cdots,N_{\rm
max}-1$. Accordingly, the operators $a_{i}^{-}$, $a_{i}^{-}$, $a$ defined by
Eqs. (18)-(23) have finite dimension $N_{\rm max}\times N_{\rm max}$, as also
specified in Table 1. Choosing a proper value of $N_{\rm max}$ is very crucial
in numerical computation. Smaller value leads to inaccuracy of the results
while larger value costs more computational resources.

Figure 8: Probability distribution of particle number in the third cell ${\cal
P}(N_{3})$. Dash line plots the exact solution which is a Poisson distribution
with the mean $\langle N_{3}\rangle=6$. The state space for a single cell is
20, stipulated by dimension $N_{\rm max}$ for truncation. Asterisks are the
numerical results from DMRG approach, while the circles from TEBD approach.
The values for numerical computations are listed in Table 1.

We are now interested in the stationary distribution of the system in steady
state. Numerically, this can be obtained with both DMRG and TEBD approaches
for the generator $\hat{L}$. The 2-norm probability that the $i$-th cell have
$N_{i}$ particles is calculated as $\tilde{\cal P}(N_{i})=\sum_{{\bf
N}\backslash N_{i}}\tilde{\cal P}({\bf N})$. This probability can be obtained
by contracting $\ket{\tilde{\cal P}({\bf N})}$ represented as an MPS with an
auxiliary MPS given by

$\displaystyle\bra{{\rm auxilary\;MPS}}={\bf 1}_{1}\otimes\cdots\otimes{\bf
1}_{i-1}\otimes{\bf V}_{i}\otimes{\bf 1}_{i+1}\otimes\cdots\otimes{\bf
1}_{L}\text{,}$ (48)

where ${\bf 1}_{j}$ is a vector with all elements one, and ${\bf V}_{i}$ a
vector with $N_{i}$-th element one and others zero. The tensor product
$\otimes$ creates a bond of dimension 1, linking the vectors to form an MPS.
Then, the 1-norm probability distribution ${\cal P}(N_{i})$ can be recovered
through Eq. (17). The mean value can also be calculated, $\langle
N_{i}\rangle=\sum_{N_{i}}N_{i}{\cal P}(N_{i})$. Figures 7 and 8 respectively
plot the profile of mean particle numbers and the probability distribution of
particle number in a given cell. The profile is linear, in agreement with the
theory (Eq. 5). The probability distribution is Poissonian, also in agreement
with the theory. We notice in Figure 8 that the probability ${\cal P}(N_{\rm
max}-1)$ is indeed vanishingly small, manifesting that there is no over-
truncation of the state space. Next, we numerically calculate the cumulant
generating function $Q(\lambda)$. The counting parameter $\lambda$ takes
several discrete values from $0$ to the affinity $A=\ln\left(\bar{N}_{\rm
L}/\bar{N}_{\rm R}\right)$. The result is presented in Figure 9. The agreement
between the numerical results and the theory (Eq. 8) is also striking. These
agreements strongly manifest validity of both DMRG and TEBD approaches in such
an application.

Although DMRG and TEBD both give nice results for the system specified with
the numbers in Table 1. DMRG is in general much more efficient than TEBD.
Because, in the adopted version of DMRG we deals with one tensor at a time,
whereas in TEBD we are simultaneously deals two tensors at a time. On the
other hand, DMRG also gives more accurate result than TEBD does. Because, in
TEBD there is an error associated with the Trotter-Suzuki decomposition.

The computer program for numerical computation is coded in C++
[Stroustrup_2013] with the ITensor library [Fishman_arXiv_2020]. Readers are
referred to Ref. [Psarras_arXiv_2021] where a comprehensive and up-to-date
snapshot of software for tensor computations is assembled.

Figure 9: Cumulant generating function $Q(\lambda)$. Dash line plots the exact
solution (8). Asterisks are the numerical results from DMRG approach, while
the circles from TEBD approach. The values for numerical computations are
listed in Table 1.

## 6 Conclusion and Perspectives

In this paper, we studied the diffusive process of Brownian particles in a
channel, with a particular interest in the current fluctuation. A
nonequilibrium constraint is imposed by placing the channel in between two
particle reservoirs with different concentrations. To describe the randomness
of the particle transport, the channel is spatially discretized into cells so
that a master equation can be written down for the time evolution of the
system state specified by the particle numbers in the cells. The counting
statistics is performed by modifying the master equation to include a counting
parameter for particle transitions between the cells. This gives a tilted
generator for the time evolution of the corresponding master equation. The
cumulant generating function for the current can be obtained from the leading
eigenvalue of the tilted generator. Because the transition rates are linear
with the local particle concentration, an analytical solution for the cumulant
generating function is obtained, exhibiting the Gallavotti-Cohen symmetry.

The focus of this paper is on the application of tensor networks to counting
statistics. We presented two approaches, DMRG and TEBD, in details. The tilted
generator is written in terms of the local annihilation, creation, particle
number operators associated with random particle transitions. The annihilation
and creation operators for the reservoir cells are defined differently to make
the nonequilibrium constraint sustainable. We also emphasized the differences
from their conventional applications for quantum many-body systems. We
introduced the so-called 2-norm probability distribution in order to make use
of the tensor-network techniques already developed in the quantum-mechanical
setting. The conventional probability distribution can be recovered once the
corresponding 2-norm probability distribution is obtained. In the DMRG
approach, the key point is to construct an MPO for the tilted generator. This
is done by constructing two MPOs for different parts of the tilted generator,
and then performing tensor summation of these two MPOs. The resulting MPO is
of bond dimension 7. We formulated the DMRG approach in such a way that a
single tensor is optimized at a time. In the TEBD approach, the tilted
evolution operator is approximately decomposed into an ordered product of
local ones. Then these local evolution operators are alternatively applied to
the an arbitrary distribution represented as an MPS. In both approaches, what
we finally obtain is the eigensolution which is then used to calculate the
cumulant generating function. When the counting parameter is set to zero, the
leading eigensolution of the generator for the master equation represents the
stationary distribution for the system in steady state. We compared the
profile of mean particle numbers and the probability distribution of particles
in a single cell with those given by the analytic solution. We remark here
that the joint probability distribution of particles in several cells can also
be obtained. We also compared the numerically obtained cumulant generating
function with its analytical counterpart. Striking agreements are found for
all these comparisons between the numerical results and theory.

The tensor-network approaches provides a good means to investigate the
fluctuations for a vast majority of realistic systems, e.g., diffusion-
reaction systems with many species of particles, systems with long-range
electrostatic interactions between particles [Andrieux_JStatMech_2009,
Gu_PhysRevE_2018, Gu_PhysRevE_2019], or systems of complex geometries with
more than two particle reservoirs. Moreover, we notice that the tensor-network
approaches, especially DMRG, enable a very accurate determination of the
cumulant generating function. This makes it possible to test Onsager
reciprocal relations [Onsager_PhysRev_1931a, Onsager_PhysRev_1931b,
Casimir_RevModPhys_1945] and their generalizations to nonlinear regimes
[Andrieux_JChemPhys_2004, Andrieux_JStatMech_2007, Gaspard_NewJPhys_2013,
Barbier_JPhysA_2018, Barbier_JPhysA_2019, Barbier_JPhysA_2020] for spatially
extended systems. It also makes it possible to study the mathematical
properties for statistical cumulants. For example, whether the limited number
of cumulants can be used to infer the affinity [Gu_JStatMech_2020]. To
conclude, an accurate determination of the cumulant generating function
possibly reveals many underlying physics.

## Acknowledgments

Financial support from National Science Foundation of China under the Grants
No. 11775001, 11825501, 12147162 is acknowledged.

##  Appendix A Formal Derivation of the Cumulant Generating Function

In this appendix, we provide a formal derivation of the cumulant generating
function (8) with the method introduced in Ref. [Gaspard_JStatMech_2018a]. We
now consider the probability ${\cal P}(Z,N_{1},\cdots,N_{L},t)$ that the cells
contain given particle numbers and that the signed cumulated number $Z$ of
particles is transferred from the $I$-th to the $(I+1)$-th cells during time
interval $[0,\,t]$. The time evolution of this probability is ruled by the
following master equation,

$\displaystyle\frac{{\rm d}{\cal P}}{{\rm d}t}=$ $\displaystyle k\bar{N}_{\rm
L}\left({\rm e}^{-\partial_{N_{1}}}-1\right){\cal P}+k\left({\rm
e}^{+\partial_{N_{1}}}-1\right)N_{1}{\cal P}+k\left({\rm
e}^{+\partial_{N_{1}}}{\rm e}^{-\partial_{N_{2}}}-1\right)N_{1}{\cal P}$
$\displaystyle+\sum_{i=2}^{I-1}\left[k\left({\rm e}^{+\partial_{N_{i}}}{\rm
e}^{-\partial_{N_{i-1}}}-1\right)N_{i}{\cal P}+k\left({\rm
e}^{+\partial_{N_{i}}}{\rm e}^{-\partial_{N_{i+1}}}-1\right)N_{i}{\cal
P}\right]$ $\displaystyle+k\left({\rm e}^{+\partial_{N_{I}}}{\rm
e}^{-\partial_{N_{I-1}}}-1\right)N_{I}{\cal P}+k\left({\rm
e}^{+\partial_{N_{I}}}{\rm e}^{-\partial_{N_{I+1}}}{\rm
e}^{-\partial_{Z}}-1\right)N_{I}{\cal P}$ $\displaystyle+k\left({\rm
e}^{+\partial_{N_{I+1}}}{\rm e}^{-\partial_{N_{I}}}{\rm
e}^{+\partial_{Z}}-1\right)N_{I+1}{\cal P}+k\left({\rm
e}^{+\partial_{N_{I+1}}}{\rm e}^{-\partial_{N_{I+2}}}-1\right)N_{I+1}{\cal P}$
$\displaystyle+\sum_{i=I+2}^{L-1}\left[k\left({\rm e}^{+\partial_{N_{i}}}{\rm
e}^{-\partial_{N_{i-1}}}-1\right)N_{i}{\cal P}+k\left({\rm
e}^{+\partial_{N_{i}}}{\rm e}^{-\partial_{N_{i+1}}}-1\right)N_{i}{\cal
P}\right]$ $\displaystyle+k\bar{N}_{\rm R}\left({\rm
e}^{-\partial_{N_{L}}}-1\right){\cal P}+k\left({\rm e}^{+\partial_{N_{L}}}{\rm
e}^{-\partial_{N_{L-1}}}-1\right)N_{L}{\cal P}+k\left({\rm
e}^{+\partial_{N_{L}}}-1\right)N_{L}{\cal P}\text{.}$ (49)

We then define the moment generating function

$\displaystyle
G(\eta,s_{1},\cdots,s_{L},t)=\sum_{Z,N_{1},\cdots,N_{L}}\eta^{Z}\left(\prod_{i}s_{i}^{N_{i}}\right){\cal
P}(Z,N_{1},\cdots,N_{L},t)\text{.}$ (50)

where $\eta\equiv{\rm e}^{-\lambda}$ with $\lambda$ being the counting
parameter for the particle transfers $Z$. From the master equation (49), we
can easily obtain the following first-order partial differential equation,

$\displaystyle\partial_{t}G$
$\displaystyle+\left[k(s_{1}-1)+k(s_{1}-s_{2})\right]\partial_{s_{1}}G$
$\displaystyle+\sum_{i=2}^{I-1}\left[k(s_{i}-s_{i-1})+k(s_{i}-s_{i+1})\right]\partial_{s_{i}}G$
$\displaystyle+\left[k(s_{I}-s_{I-1})+k(s_{I}-\eta
s_{I+1})\right]\partial_{s_{I}}G$
$\displaystyle+\left[k(s_{I+1}-\eta^{-1}s_{I})+k(s_{I+1}-s_{I+2})\right]\partial_{s_{I+1}}G$
$\displaystyle+\sum_{i=I+1}^{L-1}\left[k(s_{i}-s_{i-1})+k(s_{i}-s_{i+1})\right]\partial_{s_{i}}G$
$\displaystyle+\left[k(s_{L}-s_{L-1})+k(s_{L}-1)\right]\partial_{s_{L}}G$
$\displaystyle=\left[k\bar{N}_{\rm L}(s_{1}-1)+k\bar{N}_{\rm
R}(s_{L}-1)\right]G\text{,}$ (51)

which, in vectorial notations, can be written in the following form,

$\displaystyle\partial_{t}G+({\boldsymbol{\sf L}}\cdot{\bf s}+{\bf
f})\cdot\partial_{\bf s}G=({\bf g}\cdot{\bf s}+h)G$ (52)

where

$\displaystyle{\boldsymbol{\sf L}}\equiv k\begin{pmatrix}2&-1&&&&&&\\\
-1&2&-1&&&&&\\\ &\ddots&\ddots&\ddots&&&&\\\ &&-1&2&-\eta&&&\\\
&&&-\eta^{-1}&2&-1&&\\\ &&&&\ddots&\ddots&\ddots&\\\ &&&&&-1&2&-1\\\
&&&&&&-1&2\end{pmatrix}\text{,}$ (53) $\displaystyle{\bf
s}\equiv\begin{pmatrix}s_{1}\\\ s_{2}\\\ \vdots\\\ s_{L-1}\\\
s_{L}\end{pmatrix}\text{,}\hskip 28.45274pt{\bf f}\equiv-\begin{pmatrix}k\\\
0\\\ \vdots\\\ 0\\\ k\end{pmatrix}\text{,}\hskip 56.9055pt{\bf
g}\equiv\begin{pmatrix}k\bar{N}_{\rm L}\\\ 0\\\ \vdots\\\ 0\\\ k\bar{N}_{\rm
R}\end{pmatrix}\text{,}$ (54)

and

$\displaystyle h\equiv-k\bar{N}_{\rm L}-k\bar{N}_{\rm R}\text{.}$ (55)

By setting $\eta=1$ in the matrix ${\boldsymbol{\sf L}}$, we can define
${\boldsymbol{\sf L}}_{0}$. So, we have the relations

$\displaystyle{\bf f}=-{\boldsymbol{\sf L}}_{0}\cdot{\bf 1}\text{,}$ (56)
$\displaystyle h=-{\bf g}\cdot{\bf 1}\text{,}$ (57)

where ${\bf 1}$ denotes the vector with all elements being equal to one.
Besides, the stationary values of mean particle numbers are given by

$\displaystyle{\bf\Gamma}_{0}={\boldsymbol{\sf L}}_{0}^{-1{\rm T}}\cdot{\bf
g}\text{.}$ (58)

Here, ${\rm T}$ denotes the transpose. The first-order partial differential
equation (52) can be solved by the method of characteristics [Gardiner_2004].
The solution is given by

$\displaystyle G=G_{0}\exp\left[{\bf g}\cdot{\boldsymbol{\sf
L}}^{-1}\cdot\left({\boldsymbol{\sf I}}-{\rm e}^{-{\boldsymbol{\sf
L}}t}\right)\cdot\left({\bf s}+{\boldsymbol{\sf L}}^{-1}\cdot{\bf
f}\right)+\left(h-{\bf g}\cdot{\boldsymbol{\sf L}}^{-1}\cdot{\bf
f}\right)t\right]\text{,}$ (59)

where $G_{0}$ denotes the initial condition, and ${\boldsymbol{\sf I}}$ the
identity matrix. The cumulant generating function for the current can be
obtained,

$\displaystyle
Q(\lambda)\equiv\lim_{t\to\infty}-\frac{1}{t}\ln\left[G(\eta={\rm
e}^{-\lambda},{\bf 1},t)\right]={\bf g}\cdot\left({\bf 1}+{\boldsymbol{\sf
L}}^{-1}\cdot{\bf f}\right)\text{.}$ (60)

We observe that

$\displaystyle{\boldsymbol{\sf L}}=\boldsymbol{\sf M}\cdot{\boldsymbol{\sf
L}}_{0}\cdot\boldsymbol{\sf M}^{-1}\text{,}$ (61)

where

$\displaystyle\boldsymbol{\sf M}=\eta\boldsymbol{\sf P}_{\rm
L}+\boldsymbol{\sf P}_{\rm R}$ (62)

with the projection matrices

$\displaystyle\boldsymbol{\sf P}_{\rm
L}=\begin{pmatrix}1&\cdots&0&0&\cdots&0\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\ 0&\cdots&1&0&\cdots&0\\\
0&\cdots&0&0&\cdots&0\\\ \vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&\cdots&0&0&\cdots&0\\\ \end{pmatrix}\hskip 28.45274pt\text{and}\hskip
28.45274pt\boldsymbol{\sf P}_{\rm R}=\begin{pmatrix}0&\cdots&0&0&\cdots&0\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\ 0&\cdots&0&0&\cdots&0\\\
0&\cdots&0&1&\cdots&0\\\ \vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
0&\cdots&0&0&\cdots&1\\\ \end{pmatrix}\text{.}$ (63)

The identity matrix in $\boldsymbol{\sf P}_{\rm L}$ is of dimension $I\times
I$, while the identity matrix in $\boldsymbol{\sf P}_{\rm R}$ is
$(L-I)\times(L-I)$. Since the projection matrices satisfy the condition
$\boldsymbol{\sf P}_{\rm L}+\boldsymbol{\sf P}_{\rm R}={\boldsymbol{\sf I}}$,
we thus have

$\displaystyle\boldsymbol{\sf M}={\boldsymbol{\sf I}}+(\eta-1)\boldsymbol{\sf
P}_{\rm L}\text{,}$ (64) $\displaystyle\boldsymbol{\sf
M}^{-1}={\boldsymbol{\sf I}}+(\eta^{-1}-1)\boldsymbol{\sf P}_{\rm L}\text{.}$
(65)

From the above related expressions, the cumulant generating function (60) can
be written in the following form,

$\displaystyle Q(\lambda)$ $\displaystyle={\bf g}\cdot\left({\boldsymbol{\sf
I}}-\boldsymbol{\sf M}\cdot{\boldsymbol{\sf L}}_{0}^{-1}\cdot\boldsymbol{\sf
M}^{-1}\cdot{\boldsymbol{\sf L}}_{0}\right)=W_{+}\left(1-{\rm
e}^{-\lambda}\right)+W-\left(1-{\rm e}^{+\lambda}\right)\text{,}$ (66)

with the global transition rates given by

$\displaystyle W_{+}={\bf\Gamma}_{0}\cdot{\boldsymbol{\sf
L}}_{0}\cdot\boldsymbol{\sf P}_{\rm L}\cdot{\boldsymbol{\sf
L}}_{0}^{-1}\cdot\boldsymbol{\sf P}_{\rm R}\cdot{\boldsymbol{\sf
L}}_{0}\cdot{\bf 1}\text{,}$ (67) $\displaystyle
W_{-}={\bf\Gamma}_{0}\cdot{\boldsymbol{\sf L}}_{0}\cdot\boldsymbol{\sf P}_{\rm
R}\cdot{\boldsymbol{\sf L}}_{0}^{-1}\cdot\boldsymbol{\sf P}_{\rm
L}\cdot{\boldsymbol{\sf L}}_{0}\cdot{\bf 1}\text{.}$ (68)

These two global transition rates can be further developed as

$\displaystyle W_{+}=k^{2}\bar{N}_{\rm L}\left({\boldsymbol{\sf
L}}_{0}^{-1}\right)_{1L}\hskip 28.45274pt\text{and}\hskip
28.45274ptW_{-}=k^{2}\bar{N}_{\rm R}\left({\boldsymbol{\sf
L}}_{0}^{-1}\right)_{L1}\text{.}$ (69)

Inverting the matrix ${\boldsymbol{\sf L}}_{0}$, we get

$\displaystyle\left({\boldsymbol{\sf
L}}_{0}^{-1}\right)_{ij}=\begin{cases}\frac{i}{k(L+1)}(L+1-j)&\text{if}\quad
i\leq j\text{,}\\\ &\vspace{-0.3cm}\textstyle\\\
\frac{j}{k(L+1)}(L+1-i)&\text{if}\quad i>j\text{.}\end{cases}$ (70)

So, the two global transition rates are calculated as

$\displaystyle W_{+}=\frac{k\bar{N}_{\rm L}}{L+1}\hskip
28.45274pt\text{and}\hskip 28.45274ptW_{-}=\frac{k\bar{N}_{\rm
R}}{L+1}\text{,}$ (71)

finally giving

$\displaystyle Q(\lambda)=\frac{k\bar{N}_{\rm L}}{L+1}\left(1-{\rm
e}^{-\lambda}\right)+\frac{k\bar{N}_{\rm R}}{L+1}\left(1-{\rm
e}^{\lambda}\right)\text{.}$ (72)


Comparing this with (192), we obtain

$\displaystyle T_{SD}=T_{p,1}\cdot[1+o_{p}(1)]+o_{p}(1).$ (193)

So to finish the proof, by using the Slutsky lemma, it suffices to prove that
$T_{p,1}\to(1-\sum_{i=1}^{\infty}\rho_{i}^{2})^{1/2}\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$
in distribution.

Set $m=n-1$. Then

$\displaystyle\sqrt{n}\bar{{\boldsymbol{X}}}\sim
N_{p}(\mathbb{0},{\boldsymbol{\Sigma}}),\ n\hat{{\bf S}}\sim
W_{p}(m,{\boldsymbol{\Sigma}}),\ \mbox{and}\ \bar{{\boldsymbol{X}}}\
\mbox{and}\ \hat{{\bf S}}\ \mbox{are independent},$ (194)

where $W_{p}(m,{\boldsymbol{\Sigma}})$ is the Wishart distribution defined
after (165); see, for example, Theorem 3.1.2 from Muirhead (1982). This
implies

$\displaystyle\hat{{\bf D}}\,\overset{d}{=}\,\mbox{the diagonal matrix of}\
\frac{1}{n}W_{p}(m,{\boldsymbol{\Sigma}}).$ (195)

In particular, $\bar{{\boldsymbol{X}}}$ is independent of $\hat{{\bf D}}$. By
Lemma 4.1 and assumption $\lim_{p\to\infty}\frac{p}{m\|\mathbb{R}\|_{F}}=0$,

$\displaystyle\frac{n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}}{\sqrt{2\,\mbox{tr}(\mathbb{R}^{2})}}=\frac{n\bar{{\boldsymbol{X}}}^{T}{\bf
D}^{-1}\bar{{\boldsymbol{X}}}}{\sqrt{2\,\mbox{tr}(\mathbb{R}^{2})}}+o_{p}(1)$

as $p\to\infty$, where ${\bf D}$ is the diagonal matrix of the population
covariance matrix ${\boldsymbol{\Sigma}}$. Now we will analyze the behavior of
$n\bar{{\boldsymbol{X}}}^{T}{\bf D}^{-1}\bar{{\boldsymbol{X}}}$ in two steps.
Once they are established, the limiting distribution of $T_{p,1}$ will be
identified quickly.

Step 1. The exact distribution of $n\bar{{\boldsymbol{X}}}^{T}{\bf
D}^{-1}\bar{{\boldsymbol{X}}}$. By assumption,
$\lambda_{1},\cdots,\lambda_{p}$ are the eigenvalues of $\mathbb{R}={\bf
D}^{-1/2}{\boldsymbol{\Sigma}}{\bf D}^{-1/2}$. We claim that

$\displaystyle n\bar{{\boldsymbol{X}}}^{T}{\bf
D}^{-1}\bar{{\boldsymbol{X}}}\overset{d}{=}\sum_{i=1}^{p}\lambda_{i}\xi_{i}^{2},$
(196)

where $\xi_{1},\cdots,\xi_{p}$ are i.i.d. $N(0,1)$ and where
“$\overset{d}{=}$” means both sides of “$=$” have the same distribution. We
show (196) next.

Since $\sqrt{n}\bar{{\boldsymbol{X}}}\sim
N_{p}(\mathbb{0},{\boldsymbol{\Sigma}})$, we are able to write
$\sqrt{n}\bar{{\boldsymbol{X}}}={\boldsymbol{\Sigma}}^{1/2}\boldsymbol{\xi}$,
where ${\boldsymbol{\Sigma}}^{1/2}$ is a non-negative definite matrix
satisfying
${\boldsymbol{\Sigma}}^{1/2}\cdot{\boldsymbol{\Sigma}}^{1/2}={\boldsymbol{\Sigma}}$
and $\boldsymbol{\xi}\sim N_{p}(\mathbb{0},\mathbb{I}_{p}).$ It follows that

$\displaystyle n\bar{{\boldsymbol{X}}}^{T}{\bf
D}^{-1}\bar{{\boldsymbol{X}}}=\boldsymbol{\xi}^{T}\big{(}{\boldsymbol{\Sigma}}^{1/2}{\bf
D}^{-1}{\boldsymbol{\Sigma}}^{1/2}\big{)}\boldsymbol{\xi}.$ (197)

Since $\mathbb{A}\mathbb{B}$ and $\mathbb{B}\mathbb{A}$ have the same
eigenvalues for any square matrix $\mathbb{A}$ and $\mathbb{B}$, it is easy to
see that $\lambda_{1},\cdots,\lambda_{p}$ are also the eigenvalues of
${\boldsymbol{\Sigma}}^{1/2}{\bf D}^{-1}{\boldsymbol{\Sigma}}^{1/2}$. Write
${\boldsymbol{\Sigma}}^{1/2}{\bf
D}^{-1}{\boldsymbol{\Sigma}}^{1/2}=\mathbb{O}^{T}\mbox{diag}(\lambda_{1},\cdots,\lambda_{p})\mathbb{O}$
for some orthogonal matrix $\mathbb{O}$. Then by the orthogonal invariance
property, we know $\mathbb{O}\boldsymbol{\xi}\overset{d}{=}\boldsymbol{\xi}$.
Consequently,

$\displaystyle\boldsymbol{\xi}^{T}\big{(}{\boldsymbol{\Sigma}}^{1/2}{\bf
D}^{-1}{\boldsymbol{\Sigma}}^{1/2}\big{)}\boldsymbol{\xi}=(\mathbb{O}\boldsymbol{\xi})^{T}\mbox{diag}(\lambda_{1},\cdots,\lambda_{p})(\mathbb{O}\boldsymbol{\xi})\overset{d}{=}\sum_{i=1}^{p}\lambda_{i}\xi_{i}^{2}$
(198)

where $\xi_{1},\cdots,\xi_{p}$ are i.i.d. $N(0,1)$. We then get (196).

Step 2: the limiting distribution of $n\bar{{\boldsymbol{X}}}^{T}{\bf
D}^{-1}\bar{{\boldsymbol{X}}}$. By (196) and the fact
$\mbox{tr}(\mathbb{R})=\lambda_{1}+\cdots+\lambda_{p}=p$, we know

$\displaystyle\frac{n\bar{{\boldsymbol{X}}}^{T}{\bf
D}^{-1}\bar{{\boldsymbol{X}}}-p}{\sqrt{2\,\mbox{tr}(\mathbb{R}^{2})}}\overset{d}{=}\frac{1}{\sqrt{2}}\sum_{i=1}^{p}\frac{\lambda_{i}}{\|\mathbb{R}\|_{F}}(\xi_{i}^{2}-1).$

Set $a_{p,i}=\frac{\lambda_{i}}{\|\mathbb{R}\|_{F}}$ for $1\leq i\leq p$. Then
$a_{p,1}\geq\cdots\geq a_{p,p}\geq 0$ and $a_{p,1}^{2}+\cdots+a_{p,p}^{2}=1$.
By assumption (a), $\lim_{p\to\infty}a_{p,i}=\rho_{i}\geq 0$ for each $i\geq
1$. From Lemma 4.2,

$\displaystyle\frac{1}{\sqrt{2}}\sum_{i=1}^{p}\frac{\lambda_{i}}{\|\mathbb{R}\|_{F}}(\xi_{i}^{2}-1)\to\Big{(}1-\sum_{i=1}^{\infty}\rho_{i}^{2}\Big{)}^{1/2}\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$

in distribution as $p\to\infty$, where $\xi_{0}\sim N(0,1)$ and $\xi_{0}$ is
independent of $\\{\xi_{i};\,i\geq 1\\}$. Therefore,

$\displaystyle\frac{n\bar{{\boldsymbol{X}}}^{T}{\bf
D}^{-1}\bar{{\boldsymbol{X}}}-p}{\sqrt{2\,\mbox{tr}(\mathbb{R}^{2})}}\to\Big{(}1-\sum_{i=1}^{\infty}\rho_{i}^{2}\Big{)}^{1/2}\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$
(199)

in distribution.

With Step 1 and Step 2 completed, let us now proceed to finish the proof. In
fact, by assumption (b),

$\displaystyle\frac{1}{\sqrt{2}\|\mathbb{R}\|_{F}}\big{[}p-pn(n-3)^{-1}\big{]}=-\frac{3p}{\sqrt{2}(n-3)\|\mathbb{R}\|_{F}}\to
0.$

Summing this and (199), we see from the Slutsky lemma that

$\displaystyle\frac{n\bar{{\boldsymbol{X}}}^{T}{\bf
D}^{-1}\bar{{\boldsymbol{X}}}-pn(n-3)^{-1}}{\sqrt{2}\|\mathbb{R}\|_{F}}\to\Big{(}1-\sum_{i=1}^{\infty}\rho_{i}^{2}\Big{)}^{1/2}\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$

in distribution. Combine this with (191) to see

$\displaystyle
T_{p,1},\to\Big{(}1-\sum_{i=1}^{\infty}\rho_{i}^{2}\Big{)}^{1/2}\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$

in distribution as $p\to\infty.$ The proof is completed. $\square$

Now we begin to prove Proposition 1. As far as the proofs go, parts (i) and
(ii) from the lemma have different natures, to make the presentation clear, we
will handle the two parts separately.

Proof of Proposition 1(i). Assume $\boldsymbol{\mu}=0$. Then
$\mathbb{X}_{1},\cdots,\mathbb{X}_{n}$ are i.i.d. with distribution
$N_{p}(\mathbb{0},{\boldsymbol{\Sigma}})$, where all of the $p^{2}$ entries of
${\boldsymbol{\Sigma}}$ are identical to $\sigma^{2}$. For this reason, we
write ${\boldsymbol{X}}_{i}=\xi_{i}(1,\cdots,1)^{T}\in\mathbb{R}^{p}$ for
$i=1,\cdots,n$, where $\xi_{1},\cdots,\xi_{n}$ are i.i.d. $N(0,\sigma^{2})$.
By (1) and (7),
$\bar{{\boldsymbol{X}}}=\frac{1}{n}\sum_{i=1}^{n}{\boldsymbol{X}}_{i}$,

$\displaystyle\hat{{\bf
S}}=\frac{1}{n}\sum_{i=1}^{n}({\boldsymbol{X}}_{i}-\bar{{\boldsymbol{X}}})({\boldsymbol{X}}_{i}-\bar{{\boldsymbol{X}}})^{T}~{}~{}~{}~{}~{}\mbox{and}~{}~{}~{}~{}~{}\hat{{\bf
R}}=\hat{{\bf D}}^{-1/2}\hat{{\bf S}}\hat{{\bf D}}^{-1/2},$

where $\hat{{\bf D}}$ is the diagonal matrix of $\hat{{\bf S}}$. Set
$\bar{\xi}=\frac{1}{n}\sum_{i=1}^{n}\xi_{i}$ and
$W=\frac{1}{n}\sum_{i=1}^{n}(\xi_{i}-\bar{\xi})^{2}$. Then
$\sqrt{n}\bar{\xi}/\sigma\sim N(0,1)$, $nW/\sigma^{2}\sim\chi^{2}(n-1)$, and
$\bar{\xi}$ and $W$ are independent. Since
$\bar{{\boldsymbol{X}}}=\bar{\xi}(1,\cdots,1)^{T}\in\mathbb{R}^{p}$, we see
${\boldsymbol{X}}_{i}-\bar{{\boldsymbol{X}}}=(\xi_{i}-\bar{\xi})^{T}(1,\cdots,1)^{T}\in\mathbb{R}^{p}$.
Hence, $\hat{{\bf S}}=W\cdot\mathbb{J}$, where $\mathbb{J}$ is a $p\times p$
matrix whose entries are all equal to $1$. It follows $\hat{{\bf
D}}=W\cdot\mathbb{I}_{p}$ and $\hat{{\bf R}}=\mathbb{J}$. In particular,
$\mathrm{tr}(\hat{{\bf R}}^{2})=p^{2}$. Consequently,

$\displaystyle n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}=\frac{np}{n-1}\cdot\frac{(\sqrt{n}\bar{\xi}/\sigma)^{2}}{nW/[(n-1)\sigma^{2}]}\overset{d}{=}\frac{np}{n-1}\cdot
F_{1,n-1}.$ (200)

By (2),

$\displaystyle T_{SD}=\frac{n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}-p(n-1)(n-3)^{-1}}{\sqrt{2[\mathrm{tr}(\hat{{\bf
R}}^{2})-p^{2}(n-1)^{-1}]}}$ $\displaystyle\overset{d}{=}$
$\displaystyle\frac{\frac{n}{n-1}\cdot
F_{1,n-1}-\frac{n-1}{n-3}}{\sqrt{2[1-(n-1)^{-1}]}}.$ (201)

Notice $F_{1,n-1}\to\chi^{2}(1)$ as $n\to\infty$. We see from the Slutsky
lemma that $T_{SD}\to\frac{1}{\sqrt{2}}\cdot[\chi^{2}(1)-1]$ as $p\to\infty$.
Finally, by (9) and (200),

$\displaystyle T_{p,1}=\frac{n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}-pn(n-3)^{-1}}{\sqrt{2\big{|}\mathrm{tr}(\hat{{\bf
R}}^{2})-p(p-1)(n-1)^{-1}\big{|}}}\overset{d}{=}\frac{\frac{np}{n-1}\cdot
F_{1,n-1}-pn(n-3)^{-1}}{\sqrt{2\big{|}p^{2}-p(p-1)(n-1)^{-1}\big{|}}}.$

As a consequence,

$\displaystyle T_{p,1}$ $\displaystyle\overset{d}{=}$
$\displaystyle\frac{\frac{n}{n-1}\cdot
F_{1,n-1}-n(n-3)^{-1}}{\sqrt{2\\{1-(p-1)/[(n-1)p]\\}}}\to\frac{1}{\sqrt{2}}\cdot\big{[}\chi^{2}(1)-1\big{]}.$

The proof is completed. $\square$

Proof of Proposition 1(ii). Recall

$\displaystyle T_{SD}=\frac{n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}-p(n-1)(n-3)^{-1}}{\sqrt{2[\mathrm{tr}(\hat{{\bf
R}}^{2})-p^{2}(n-1)^{-1}]}}.$

Since $T_{SD}$ is scale-invariant, without loss of generality, we assume
${\boldsymbol{\Sigma}}=\mathbb{I}_{p}$. By (194),

$\displaystyle\sqrt{n}\bar{{\boldsymbol{X}}}\sim
N_{p}(\mathbb{0},\mathbb{I}_{p}),\ n\hat{{\bf S}}\sim
W_{p}(m,\mathbb{I}_{p}),\ \mbox{and}\ \bar{{\boldsymbol{X}}}\ \mbox{and}\
\hat{{\bf S}}\ \mbox{are independent},$

where $m:=n-1$ and $W_{p}(m,{\boldsymbol{\Sigma}})$ is the Wishart
distribution defined after (165). Also, by (195), $\hat{{\bf D}}$ and the
diagonal matrix of $\frac{1}{n}W_{p}(m,\mathbb{I}_{p})$ have the same
distribution. By definition of $W_{p}(m,\mathbb{I}_{p})$, its $p$ diagonal
entries are i.i.d. $\chi^{2}(m)$. Therefore,
$n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf D}}^{-1}\bar{{\boldsymbol{X}}}$ is a sum
of $p$ i.i.d. random variables with distribution
$n\cdot\frac{N(0,1)^{2}}{\chi^{2}(m)}$, where the numerator and denominator
are independent. Thus, we are able to write

$\displaystyle n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}=\frac{n}{m}\sum_{i=1}^{p}X_{p,i},$ (202)

where $X_{p,1},\cdots,X_{p,p}$ are i.i.d. with distribution $F(1,m)$. By Lemma
6.36,

$\displaystyle
U_{p}:=\frac{1}{\sqrt{2p}}\Big{(}\sum_{i=1}^{p}X_{p,i}-\frac{mp}{m-2}\Big{)}\to
N(0,1)$ (203)

as $p\to\infty$, regardless of the speeds of $n$ and $p$ going to infinity.
Solve $\sum_{i=1}^{p}X_{p,i}$ in terms of $U_{p}$ to see
$\sum_{i=1}^{p}X_{p,i}=\sqrt{2p}U_{p}+\frac{mp}{m-2}$. By plugging this into
(202), we obtain

$\displaystyle n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}=\frac{n}{m}\sqrt{2p}U_{p}+\frac{np}{m-2}.$
(204)

Consequently,

$\displaystyle n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}-p(n-1)(n-3)^{-1}$ $\displaystyle=$
$\displaystyle\frac{n}{m}\sqrt{2p}U_{p}+\frac{np}{m-2}-\frac{mp}{m-2}$
$\displaystyle=$
$\displaystyle\sqrt{2p}U_{p}+\frac{\sqrt{2p}}{m}U_{p}+\frac{p}{m-2}.$

It follows that

$\displaystyle\frac{1}{\sqrt{2p}}\big{[}n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}-p(n-1)(n-3)^{-1}\big{]}=U_{p}+\frac{1}{m}U_{p}+\frac{\sqrt{p}}{\sqrt{2}(m-2)}.$
(205)

On the other hand, $\hat{{\bf R}}=\hat{{\bf D}}^{-1/2}\hat{{\bf S}}\hat{{\bf
D}}^{-1/2}$ by (7). Write $n\hat{{\bf
S}}=(\boldsymbol{\xi}_{1},\cdots,\boldsymbol{\xi}_{p})^{T}(\boldsymbol{\xi}_{1},\cdots,\boldsymbol{\xi}_{p})$,
where $\boldsymbol{\xi}_{1},\cdots,\boldsymbol{\xi}_{p}$ are i.i.d.
$N_{m}(0,\mathbb{I}_{m})$. Set
$\mathbb{e}_{i}=\frac{\boldsymbol{\xi}_{i}}{\|\boldsymbol{\xi}_{1}\|}$ for
$i=1,\cdots,p$. Then $\\{\mathbb{e}_{1},\cdots,\mathbb{e}_{p}\\}$ are i.i.d.
uniformly distributed on the $m$-dimensional sphere; see, for example, p. 38
from Muirhead (1982). Also, $\hat{{\bf
R}}=(\mathbb{e}_{i}^{T}\mathbb{e}_{j})$. Set $V_{p}:=\frac{1}{p}\sum_{1\leq
i<j\leq p}\big{[}m(\mathbb{e}_{i}^{\prime}\mathbb{e}_{j})^{2}-1\big{]}.$ We
claim

$\displaystyle\sup_{p\geq 2}\mbox{Var}(V_{p})\leq\frac{3}{2}.$ (206)

In fact, by Lemma 11 from Cai et al. (2013) or (22) in Lemma 4.1 from Cai and
Jiang (2012), $\\{\mathbb{e}_{i}^{T}\mathbb{e}_{j};\,1\leq i<j\leq p\\}$ are
pairwise i.i.d. Thus,
$\mbox{Var}(V_{p})=\frac{1}{p^{2}}\cdot\frac{1}{2}p(p-1)\cdot
m^{2}\cdot\mbox{Var}((\mathbb{e}_{1}^{\prime}\mathbb{e}_{2})^{2}).$ Since
$\mathbb{e}_{1}$ and $\mathbb{e}_{2}$ are independent, by Theorems 1.5.6 and
1.5.7 from Muirhead (1982), we know
$\mathbb{e}_{1}^{T}\mathbb{e}_{2}\overset{d}{=}\xi_{1}(\xi_{1}^{2}+\cdots+\xi_{m}^{2})^{-1/2}$,
where $\xi_{1},\cdots,\xi_{m}$ are i.i.d. $N(0,1)$. In particular, if $n=2$
then $m=1$, $\mathbb{e}_{1}^{T}\mathbb{e}_{2}$ is a symmetric Bernoulli random
variable, and hence $(\mathbb{e}_{1}^{T}\mathbb{e}_{2})^{2}=1$ and
$\mbox{Var}((\mathbb{e}_{1}^{T}\mathbb{e}_{2}))=0$. Now, for $n\geq 3$,

$\displaystyle\mbox{Var}((\mathbb{e}_{1}^{\prime}\mathbb{e}_{2})^{2})=\mbox{Var}\Big{(}\frac{\xi_{1}^{2}}{\xi_{1}^{2}+\cdots+\xi_{m}^{2}}\Big{)}\leq
E\frac{\xi_{1}^{4}}{(\xi_{1}^{2}+\cdots+\xi_{m}^{2})^{2}}=\frac{3}{m(m+2)}$

by taking $a_{1}=2$ and other $a_{i}=0$ in Lemma 2.4 from Jiang (2012). The
inequality is true for all $m\geq 1$. Combining this and the earlier
expression of $\mbox{Var}(V_{p})$, we obtain (206). In particular, by the
Chebyshev inequality, (206) indicates

$\displaystyle\frac{V_{p}}{m}\to 0$ (207)

as $p\to\infty$ regardless of the speeds of $n$ and $p$ going to infinity.
Write

$\displaystyle\sum_{1\leq i<j\leq
p}(\mathbb{e}_{i}^{\prime}\mathbb{e}_{j})^{2}=\frac{p}{m}V_{p}+\frac{p(p-1)}{2m}.$

Since $\hat{{\bf R}}$ is a symmetric matrix whose diagonal entries are all
equal to $1$. The above implies

$\displaystyle\mathrm{tr}(\hat{{\bf
R}}^{2})=p+\frac{2p}{m}V_{p}+\frac{p(p-1)}{m}.$ (208)

Consequently, we have from (207) that

$\displaystyle\frac{2}{p}\cdot\Big{[}\mathrm{tr}(\hat{{\bf
R}}^{2})-\frac{p^{2}}{m}\Big{]}=2+\frac{4}{m}V_{p}-\frac{2}{m}\to 2$ (209)

in probability. It follows from this and (205) that

$\displaystyle
T_{SD}=\frac{\sqrt{2p}}{\sqrt{p(2+4m^{-1}V_{p}-2m^{-1})}}\cdot\Big{[}U_{p}+\frac{1}{m}U_{p}+\frac{\sqrt{p}}{\sqrt{2}(m-2)}\Big{]}.$

By (207) and the Slutsky lemma,
$\frac{\sqrt{2p}}{\sqrt{p(2+4m^{-1}V_{p}-2m^{-1})}}\to 1$ in probability.
Therefore, by (203),

$\displaystyle
U_{p}+\frac{1}{m}U_{p}+\frac{\sqrt{p}}{\sqrt{2}(m-2)}\to\begin{cases}\xi_{0},&\text{if
$p/n^{2}\to 0$};\\\ \xi_{0}+\sqrt{h/2},&\text{if $p/n^{2}\to h$};\\\
\infty,&\text{if $p/n^{2}\to\infty$}\end{cases}$

in distribution, where $\xi_{0}\sim N(0,1)$. So we get part (ii) of
Proposition 1.

Finally, as for $T_{p,1},$ in lieu of (9),

$\displaystyle T_{p,1}=\frac{n\bar{{\boldsymbol{X}}}^{T}\hat{{\bf
D}}^{-1}\bar{{\boldsymbol{X}}}-pn(n-3)^{-1}}{\sqrt{2\big{|}\mathrm{tr}(\hat{{\bf
R}}^{2})-p(p-1)(n-1)^{-1}\big{|}}}=\frac{\frac{n}{m}\sqrt{2p}U_{p}}{\sqrt{2p+\frac{4p}{m}V_{p}}}=\frac{\frac{n}{m}U_{p}}{\sqrt{1+\frac{2}{m}V_{p}}}$

by (204) and (208). This implies $T_{p,1}\to N(0,1)$ as $p\to\infty$ by (203)
and (207), regardless of the speeds of $n$ and $p$ going to infinity.
$\square$

Proof of Theorem 2. Recall

$\displaystyle\bar{\mathbb{X}}_{i}=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\mathbb{X}_{ij}\
\mbox{for}\ i=1,2;$ $\displaystyle\hat{{\bf
S}}=\frac{1}{n_{1}+n_{2}}\Big{[}\sum_{j=1}^{n_{1}}(\mathbb{X}_{1j}-\bar{\mathbb{X}}_{1})(\mathbb{X}_{1j}-\bar{\mathbb{X}}_{1})^{T}+\sum_{j=1}^{n_{2}}(\mathbb{X}_{2j}-\bar{\mathbb{X}}_{2})(\mathbb{X}_{2j}-\bar{\mathbb{X}}_{2})^{T}\Big{]}.$

The $p\times p$ matrix $\hat{{\bf D}}$ is the diagonal matrix of $\hat{{\bf
S}}$ and $\hat{{\bf R}}=\hat{{\bf D}}^{-1/2}\hat{{\bf S}}\hat{{\bf D}}^{-1/2}$
is the $p\times p$ pooled sample correlation matrix. In particular, all of the
diagonal entries of $\hat{{\bf R}}$ are $1$. Thus,
$\|\mathbb{R}\|_{F}\geq\sqrt{p}$. Now, $\bar{\mathbb{X}}_{i}$ is independent
of
$\sum_{j=1}^{n_{i}}(\mathbb{X}_{ij}-\bar{\mathbb{X}}_{i})(\mathbb{X}_{ij}-\bar{\mathbb{X}}_{i})^{T}$
for each $i=1,2$; see, for example, Theorem 3.1.2 from Muirhead (1982).
Second, by assumption, the two samples are independent. Therefore, the three
random vectors $\bar{\mathbb{X}}_{1}$, $\bar{\mathbb{X}}_{2}$ and $\hat{{\bf
S}}$ are independent. In particular, $\bar{\mathbb{X}}_{1}$ and
$\bar{\mathbb{X}}_{2}$ are independent of $\hat{{\bf D}}$. By assumption
${\boldsymbol{\mu}}_{1}={\boldsymbol{\mu}}_{2}$, we obtain

$\displaystyle\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\sim
N_{p}\Big{(}\mathbb{0},\Big{(}\frac{1}{n_{1}}+\frac{1}{n_{2}}\Big{)}{\boldsymbol{\Sigma}}\Big{)}.$

Hence,

$\displaystyle\sqrt{\frac{n_{1}n_{2}}{n_{1}+n_{2}}}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}\sim
N_{p}(\mathbb{0},{\boldsymbol{\Sigma}}).$ (210)

Furthermore, for each $i$,
$\sum_{j=1}^{n_{i}}(\mathbb{X}_{ij}-\bar{\mathbb{X}}_{i})(\mathbb{X}_{ij}-\bar{\mathbb{X}}_{i})^{T}\sim
W_{p}(n_{i}-1,{\boldsymbol{\Sigma}})$; see, for example, Theorem 3.1.2 from
Muirhead (1982). Therefore, $(n_{1}+n_{2})\hat{{\bf S}}\sim
W_{p}(n_{1}+n_{2}-2,{\boldsymbol{\Sigma}})$ by independence. Consequently,

$\displaystyle\hat{{\bf D}}_{1}:=\frac{n_{1}+n_{2}}{n_{1}+n_{2}-1}\hat{{\bf
D}}\ \overset{d}{=}\,\mbox{the diagonal matrix of}\
\frac{1}{m}W_{p}\big{(}m-1,{\boldsymbol{\Sigma}}\big{)},$ (211)

where $m:=n_{1}+n_{2}-1.$ As explained earlier, the left hand side of (210) is
independent of $\hat{{\bf D}}_{1}$. In particular, by replacing
“$\boldsymbol{\eta}$” and “$n$” from (174) with
“$\sqrt{\frac{n_{1}n_{2}}{n_{1}+n_{2}}}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}$”
and “$m$”, respectively, we obtain

$\displaystyle
E\Big{[}\frac{n_{1}n_{2}}{n_{1}+n_{2}-1}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}^{T}\hat{{\bf
D}}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})\Big{]}$ (212)
$\displaystyle=$ $\displaystyle
E\Big{[}\frac{n_{1}n_{2}}{n_{1}+n_{2}}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}^{T}\hat{{\bf
D}}_{1}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})\Big{]}$
$\displaystyle=$ $\displaystyle\frac{(n_{1}+n_{2}-1)p}{n_{1}+n_{2}-4}.$

This is the reason how come up with the numerator of $T_{p,2}$ defined in
(16). The assumption $\lim_{p\to\infty}\frac{p}{(n_{1}+n_{2})^{a}}=0$ implies
$n_{1}+n_{2}\to\infty$ as $p\to\infty.$ Thus we know from assumption (b) that
$\lim_{p\to\infty}\frac{p}{m\|\mathbb{R}\|_{F}}=0$. Replacing “$n$” with “$m$”
in Lemma 4.1, we have from (210) and (211) that

$\displaystyle\frac{\frac{n_{1}n_{2}}{n_{1}+n_{2}}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}^{T}\hat{{\bf
D}}_{1}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})-p}{\sqrt{2\mathrm{tr}({\bf
R}^{2})}}$ (213) $\displaystyle=$
$\displaystyle\frac{\frac{n_{1}n_{2}}{n_{1}+n_{2}}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}^{T}{\bf
D}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})-p}{\sqrt{2\,\mbox{tr}(\mathbb{R}^{2})}}+o_{p}(1),$

where ${\bf D}$ is the diagonal matrix of ${\boldsymbol{\Sigma}}$, ${\bf
R}={\bf D}^{-1/2}{\boldsymbol{\Sigma}}{\bf D}^{-1/2}$ is the population
correlation matrix. By (196) and (210),

$\displaystyle\frac{n_{1}n_{2}}{n_{1}+n_{2}}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}^{T}{\bf
D}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})\overset{d}{=}\sum_{i=1}^{p}\lambda_{i}\xi_{i}^{2},$
(214)

where $\lambda_{1}\geq\cdots\geq\lambda_{p}\geq 0$ are the eigenvalues of
$\mathbb{R}$ and $\xi_{1},\cdots,\xi_{p}$ are i.i.d. $N(0,1).$ Similar to the
argument as in Step 2 in the proof of Theorem 1, we see

$\displaystyle\frac{1}{\|\mathbb{R}\|_{F}}\Big{(}-p+\sum_{i=1}^{p}\lambda_{i}\xi_{i}^{2}\Big{)}\to\sqrt{2}\Big{(}1-\sum_{i=1}^{\infty}\rho_{i}^{2}\Big{)}^{1/2}\xi_{0}+\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$

in distribution, where $\xi_{0}\sim N(0,1)$ and $\xi_{0}$ is independent of
$\xi_{1},\xi_{2},\cdots$. Combining this with (213) and (214), we have

$\displaystyle\frac{\frac{n_{1}n_{2}}{n_{1}+n_{2}-1}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}^{T}\hat{{\bf
D}}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})-p}{\sqrt{2\mathrm{tr}({\bf
R}^{2})}}$ (215) $\displaystyle=$
$\displaystyle\frac{\frac{n_{1}n_{2}}{n_{1}+n_{2}}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}^{T}\hat{{\bf
D}}_{1}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})-p}{\sqrt{2\mathrm{tr}({\bf
R}^{2})}}$ $\displaystyle\to$
$\displaystyle\Big{(}1-\sum_{i=1}^{\infty}\rho_{i}^{2}\Big{)}^{1/2}\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$

in distribution. Evidently,

$\displaystyle\frac{1}{\sqrt{2\mathrm{tr}({\bf
R}^{2})}}\Big{[}p-\frac{(n_{1}+n_{2}-1)p}{n_{1}+n_{2}-4}\Big{]}=\frac{1}{\sqrt{2}}\frac{-3p}{(n_{1}+n_{2}-4)\|\mathbb{R}\|_{F}}\to
0$

by the assumption
$\lim_{p\to\infty}\frac{p}{(n_{1}+n_{2})\|\mathbb{R}\|_{F}}=0$. Add the left
hand sides of the above two assertions to obtain

$\displaystyle\frac{\frac{n_{1}n_{2}}{n_{1}+n_{2}-1}\big{(}\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2}\big{)}^{T}\hat{{\bf
D}}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})-\frac{(n_{1}+n_{2}-1)p}{n_{1}+n_{2}-4}}{\sqrt{2\mathrm{tr}({\bf
R}^{2})}}$ (216) $\displaystyle\to$
$\displaystyle\Big{(}1-\sum_{i=1}^{\infty}\rho_{i}^{2}\Big{)}^{1/2}\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$

in distribution. Next we will replace $\mathrm{tr}({\bf R}^{2})$ by its ratio-
unbiased-estimator.

Recall $(n_{1}+n_{2})\hat{{\bf S}}\sim
W_{p}(n_{1}+n_{2}-2,{\boldsymbol{\Sigma}})$, and $\hat{{\bf D}}$ is the
diagonal matrix of $\hat{{\bf S}}$ and $\hat{{\bf R}}=\hat{{\bf
D}}^{-1/2}\hat{{\bf S}}\hat{{\bf D}}^{-1/2}$. Set
$\kappa=\frac{n_{1}+n_{2}}{n_{1}+n_{2}-1}$. Then

$\displaystyle\kappa\hat{{\bf
S}}\sim\frac{1}{n_{1}+n_{2}-1}W_{p}(n_{1}+n_{2}-2,{\boldsymbol{\Sigma}});$
$\displaystyle\hat{{\bf R}}=(\kappa\hat{{\bf D}})^{-1/2}(\kappa\hat{{\bf
S}})(\kappa\hat{{\bf D}})^{-1/2}\ \ \mbox{and}\ \ \kappa\hat{{\bf D}}\
\mbox{is the diagonal matrix of}\ \kappa\hat{{\bf S}}.$

The essential assumption from Theorem 3 is that, in its own notation,
$\hat{\mathbb{R}}$ is the sample correlation matrix obtained from $\hat{{\bf
S}}\sim\frac{1}{n}W_{p}(n-1,{\boldsymbol{\Sigma}})$. Replace “$\hat{{\bf S}}$”
and “$n$” from Theorem 3 with “$\kappa\hat{{\bf S}}$” and “$n_{1}+n_{2}-1$”,
respectively. Assumption (b) indicates
$\lim_{p\to\infty}\frac{p}{(n_{1}+n_{2}-1)\|\mathbb{R}\|_{F}}=0$ and
$\lim_{p\to\infty}\frac{p}{(n_{1}+n_{2}-1)^{a}}=0$. Then by Theorem 3,

$\displaystyle\frac{1}{\mbox{tr}(\mathbb{R}^{2})}\Big{[}\mbox{tr}(\hat{\mathbb{R}}^{2})-\frac{p(p-1)}{n_{1}+n_{2}-2}\Big{]}\to
1$ (217)

in probability. This together with (216) and the Slutsky lemma yields

$\displaystyle T_{p,2}$ $\displaystyle=$
$\displaystyle\frac{\frac{n_{1}n_{2}}{n_{1}+n_{2}-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})^{T}\hat{{\bf
D}}^{-1}(\bar{{\boldsymbol{X}}}_{1}-\bar{{\boldsymbol{X}}}_{2})-\frac{(n_{1}+n_{2}-1)p}{n_{1}+n_{2}-4}}{\sqrt{2\big{|}\mathrm{tr}(\hat{{\bf
R}}^{2})-\frac{p(p-1)}{n_{1}+n_{2}-2}\big{|}}}$ $\displaystyle\to$
$\displaystyle\Big{(}1-\sum_{i=1}^{\infty}\rho_{i}^{2}\Big{)}^{1/2}\xi_{0}+\frac{1}{\sqrt{2}}\sum_{i=1}^{\infty}\rho_{i}(\xi_{i}^{2}-1)$

in distribution as $p\to\infty$.

Finally, by the same argument as deriving (193), we see
$T_{SD}^{\prime}=T_{p,2}\cdot[1+o_{p}(1)]+o_{p}(1).$ So the conclusion for
$T_{SD}^{\prime}$ follows from the above and the Slutsky lemma. The proof is
completed. $\square$

###  6.4 A Lemma and Verification of (11), (12) and (13)

The following lemma is used in the discussion after Theorem 2.

######  LEMMA 6.37

Let $\\{\mathbb{R}_{p};\,p\geq 1\\}$ be non-negative definite matrices whose
diagonal entries are all equal to $1$. Let
$\lambda_{p,1}\geq\cdots\geq\lambda_{p,p}\geq 0$ be the eigenvalues of
$\mathbb{R}_{p}$. Assume condition “$C4$” stated in Theorem 2 from Zhang et
al. (2020b) hold, that is,
$\lim_{p\to\infty}\frac{\lambda_{p,i}}{\|\mathbb{R}_{p}\|_{F}}=\rho_{i}$ for
all $i\geq 1$ with $\rho_{1}>0$ and
$\lim_{p\to\infty}\sum_{i=1}^{p}\frac{\lambda_{p,i}}{\|\mathbb{R}_{p}\|_{F}}=\sum_{i=1}^{\infty}\rho_{i}<\infty$.
Then $\sum_{i=1}^{\infty}\rho_{i}^{2}=1$.

Proof of Lemma 6.37. Notice
$\sum_{i=1}^{p}\lambda_{p,i}=\mbox{tr}(\mathbb{R}_{p})=p$ because all of the
diagonal entries of $\mathbb{R}_{p}$ are identical to $1$. The assumption
$\lim_{p\to\infty}\sum_{i=1}^{p}\frac{\lambda_{p,i}}{\|\mathbb{R}_{p}\|_{F}}=\sum_{i=1}^{\infty}\rho_{i}<\infty$
implies that

$\displaystyle\frac{p}{\|\mathbb{R}_{p}\|_{F}}\to\sum_{i=1}^{\infty}\rho_{i}$
(218)

as $p\to\infty$. Set $a_{p,i}=\frac{\lambda_{p,i}}{\|\mathbb{R}_{p}\|_{F}}.$
Then $\lim_{p\to\infty}a_{p,i}=\rho_{i}$ and
$a_{p,1}^{2}+\cdots+a_{p,p}^{2}=1$ by the definition of
$\|\mathbb{R}_{p}\|_{F}$. For any $K\geq 1$, write

$\displaystyle a_{p,1}^{2}+\cdots+a_{p,K}^{2}=1-\sum_{i=K+1}^{p}a_{p,i}^{2}.$
(219)

We claim that
$\lim_{K\to\infty}\limsup_{p\to\infty}\sum_{i=K+1}^{p}a_{p,i}^{2}=0$. If this
is true, by letting $p\to\infty$ first and then sending $K\to\infty$, then
$\sum_{i=1}^{\infty}\rho_{i}^{2}=1$. We now prove the claim. In fact, write

$\displaystyle\sum_{i=K+1}^{p}a_{p,i}=\frac{p}{\|\mathbb{R}_{p}\|_{F}}-\sum_{i=1}^{K}a_{p,i}.$

For fixed $K\geq 1$, let $p\to\infty$ and use (218) to have

$\displaystyle\lim_{p\to\infty}\sum_{i=K+1}^{p}a_{p,i}=\Big{(}\sum_{i=1}^{\infty}\rho_{i}\Big{)}-\sum_{i=1}^{K}\rho_{i}=\sum_{i=K+1}^{\infty}\rho_{i}.$

Use the assumption $\sum_{i=1}^{\infty}\rho_{i}<\infty$ and let $K\to\infty$
to see $\lim_{K\to\infty}\limsup_{p\to\infty}\sum_{i=K+1}^{p}a_{p,i}=0$. The
claim is then verified since $a_{p,i}\leq 1$ for all $1\leq i\leq p$.
$\square$.

In the following the notation “$A_{p}\sim B_{p}$” means that $A_{p}/B_{p}\to
1$ as $p\to\infty$.

The Verification of (11). Review Example 2.1. We have

$\displaystyle\mathbb{R}=\begin{pmatrix}1&r&r&\cdots&r\\\ r&1&r&\cdots&r\\\
\vdots&&&&\\\ r&r&r&\cdots&1\end{pmatrix}_{p\times p}.$

Easily, $\lambda_{1}=1+(p-1)r$ and $\lambda_{2}=\cdots=\lambda_{p}=1-r$. To
maintain $\mathbb{R}$ to be non-negative definite, all eigenvalues have to be
non-negative, that is, $-\frac{1}{p-1}\leq r\leq 1.$ Recall $r=r_{p}$
satisfies that $\lim_{p\to\infty}\sqrt{p}\cdot r=c$. Now we consider three
cases: $c=0$, $c\in(0,\infty)$ and $c=\infty$.

Case 1: $c=0$. In this case
$\|\mathbb{R}\|_{F}^{2}=[1+(p-1)r]^{2}+(p-1)(1-r)^{2}\sim p$ as $p\to\infty$,
thus $\|\mathbb{R}\|_{F}\sim\sqrt{p}$. Easily,
$\frac{p}{n\|\mathbb{R}\|_{F}}=O(\frac{\sqrt{p}}{n})\to 0$ provided
$p=o(n^{2})$. Also,
$\rho_{1}=\lim_{p\to\infty}\frac{\lambda_{1}}{\|\mathbb{R}\|_{F}}=0$, hence
$\lambda_{i}=0$ for every $i\geq 1$. By Theorem 1, both $T_{SD}$ and $T_{p,1}$
converge to $N(0,1)$ in distribution.

Case 2: $c\in(0,\infty)$. In this case
$\|\mathbb{R}\|_{F}^{2}=[1+(p-1)r]^{2}+(p-1)(1-r)^{2}\sim(c^{2}+1)p$ as
$p\to\infty$, hence $\|\mathbb{R}\|_{F}\sim\sqrt{c^{2}+1}\cdot\sqrt{p}$.
Readily,

$\displaystyle\frac{p}{n\|\mathbb{R}\|_{F}}=O\Big{(}\frac{\sqrt{p}}{n}\Big{)}\to
0~{}~{}~{}~{}\mbox{and}~{}~{}~{}~{}\rho_{1}=\lim_{p\to\infty}\frac{\lambda_{1}}{\|\mathbb{R}\|_{F}}\to\frac{c}{\sqrt{c^{2}+1}};$
$\displaystyle\rho_{i}=\lim_{p\to\infty}\frac{\lambda_{i}}{\|\mathbb{R}\|_{F}}=0~{}~{}\mbox{for
every}~{}~{}i\geq 2$

provided $p=o(n^{2})$. By Theorem 1, under condition $p=o(n^{2})$, we know
both $T_{SD}$ and $T_{p,1}$ converge to

$\displaystyle\frac{1}{\sqrt{c^{2}+1}}\xi_{0}+\frac{c}{\sqrt{2(c^{2}+1)}}(\xi_{1}^{2}-1)$

in distribution, where $\xi_{0}$ and $\xi_{1}$ are i.i.d. $N(0,1).$

Case 3: $c=\infty$. In this case
$\|\mathbb{R}\|_{F}^{2}=[1+(p-1)r]^{2}+(p-1)(1-r)^{2}\sim(pr)^{2}$ as
$p\to\infty$. Therefore, $\|\mathbb{R}\|_{F}\sim pr$. Notice
$\lim_{p\to\infty}p/(n\|\mathbb{R}\|_{F})=\lim_{p\to\infty}1/(nr_{p})=0$,
$\rho_{1}=\lim_{p\to\infty}\lambda_{1}/\|\mathbb{R}\|_{F}=1$ and $\rho_{i}=0$
for $i\geq 2$ provided $nr_{p}\to\infty$. So, under conditions
$nr_{p}\to\infty$ and $p=o(n^{a})$ for some constant $a>0$, we have
$T_{p,1}\to\frac{1}{\sqrt{2}}[\chi^{2}(1)-1]$. $\square$

The Verification of (12). Recall $r\in(0,1)$, $m=[p^{r}]$ and $\mathbb{A}_{m}$
from Example 2.1. By definition,

$\displaystyle\mathbb{R}=\begin{pmatrix}\mathbb{A}_{m}&\mathbb{0}\\\
\mathbb{0}&\mathbb{I}_{p-m}\end{pmatrix}.$

Then the largest eigenvalues of $\mathbb{R}$ is $\lambda_{1}=1+(m-1)r$ and the
rest of them are $\lambda_{2}=\cdots=\lambda_{m}=1-r$ and
$\lambda_{m+1}=\cdots=\lambda_{p}=1$. Then

$\displaystyle\|\mathbb{R}\|_{F}^{2}=[1+(m-1)r]^{2}+(m-1)(1-r)^{2}+(p-m)=m^{2}r^{2}+p+O(m).$

Thus,

$\displaystyle\|\mathbb{R}\|_{F}=\begin{cases}\sqrt{p}\cdot[1+o(1)],&\text{if
$0<r<\frac{1}{2}$};\\\ \sqrt{5p/4}\cdot[1+o(1)],&\text{if $r=\frac{1}{2}$};\\\
rp^{r}\cdot[1+o(1)],&\text{if $\frac{1}{2}<r<1$.}\end{cases}$

Also,

$\displaystyle\rho_{1}=\lim_{p\to\infty}\frac{\lambda_{1}}{\|\mathbb{R}\|_{F}}=\begin{cases}0,&\text{if
$0<r<\frac{1}{2}$};\\\ 1/\sqrt{5},&\text{if $r=\frac{1}{2}$};\\\ 1,&\text{if
$\frac{1}{2}<r<1$.}\end{cases}$

Obviously, $\rho_{i}=0$ for $i\geq 2.$ At last,

$\displaystyle\frac{p}{n\|\mathbb{R}\|_{F}}=\begin{cases}O(\sqrt{p}/n),&\text{if
$0<r\leq\frac{1}{2}$};\\\ O(p^{1-r}/n),&\text{if
$\frac{1}{2}<r<1$.}\end{cases}$

Assume ${\boldsymbol{\mu}}=\boldsymbol{0}$. Then, we have from Theorem 1 that
both $T_{SD}$ and $T_{p,1}$ go to

$\displaystyle\begin{cases}N(0,1),&\text{if $0<r<\frac{1}{2}$};\\\
\frac{2}{\sqrt{5}}N(0,1)+\frac{1}{\sqrt{10}}\cdot[\chi^{2}(1)-1],&\text{if
$r=\frac{1}{2}$};\\\ \frac{1}{\sqrt{2}}\cdot[\chi^{2}(1)-1],&\text{if
$\frac{1}{2}<r<1$}\end{cases}$

in distribution under condition $p=o(n^{2})$ for $0<r\leq\frac{1}{2}$ and
$p=o(n^{1/(1-r)})$ for $\frac{1}{2}<r<1$, where the random variables $N(0,1)$
and $\chi^{2}(1)$ are independent as $r=\frac{1}{2}$. $\square$

The Verification of (13). Recall $m=[\log(p+2)].$ Given $\tau\geq 0$, we have
$p^{\prime}=p-[\tau\sqrt{p}]$, $\lambda_{i}=1+\tau
2^{-i}\big{(}1-2^{-m}\big{)}^{-1}\sqrt{p}$ for $1\leq i\leq m$ and
$\lambda_{i}=1$ for $i=m+1,\cdots,p^{\prime}-1$,
$\lambda_{p^{\prime}}=1+[\tau\sqrt{p}]-\tau\sqrt{p}\in[0,1]$ and the rest of
$\lambda_{i}$ are zero. Obviously, $\lambda_{1}+\cdots+\lambda_{k}\geq k$ for
each $i=1,\cdots,p^{\prime}-1$ and $\lambda_{1}+\cdots+\lambda_{p^{\prime}}$
is identical to

$\displaystyle\sum_{i=1}^{m}\big{[}1+\tau
2^{-i}\big{(}1-2^{-m}\big{)}^{-1}\sqrt{p}\big{]}+(p^{\prime}-m-1)\cdot
1+\big{(}1+[\tau\sqrt{p}]-\tau\sqrt{p}\big{)}$ $\displaystyle=$
$\displaystyle(m+\tau\sqrt{p})+\big{(}p-[\tau\sqrt{p}]-m-1\big{)}+1+[\tau\sqrt{p}]-\tau\sqrt{p}$
$\displaystyle=$ $\displaystyle p,$

where we use the fact $\sum_{i=1}^{m}2^{-i}\big{(}1-2^{-m}\big{)}^{-1}=1$.
Lemma 4.3 tells us that there exists a correlation matrix $\mathbb{R}$ such
that $\mathbb{R}$ has eigenvalues $\lambda_{i},\,1\leq i\leq p.$ Now,

$\displaystyle\|\mathbb{R}\|_{F}^{2}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\big{[}1+\tau
2^{-i}\big{(}1-2^{-m}\big{)}^{-1}\sqrt{p}\big{]}^{2}+(p^{\prime}-m-1)\cdot
1^{2}+\big{(}1+[\tau\sqrt{p}]-\tau\sqrt{p}\big{)}^{2}$ $\displaystyle=$
$\displaystyle\tau^{2}p\big{(}1-2^{-m}\big{)}^{-2}\sum_{i=1}^{m}\frac{1}{4^{i}}+p+O\big{(}\sqrt{p}\big{)}$
$\displaystyle=$ $\displaystyle\Big{(}\frac{\tau^{2}}{3}+1\Big{)}p+o(p)$

since $\sum_{i=1}^{m}\frac{1}{4^{i}}\to\frac{1}{3}$ as $p\to\infty$. In
particular, $\lim_{p\to\infty}\frac{p}{n\|\mathbb{R}\|_{F}}=0$ provided
$p=o(n^{2})$. Easily, for each $i\geq 1$,

$\displaystyle\rho_{i}:=\lim_{p\to\infty}\frac{\lambda_{i}}{\|\mathbb{R}\|_{F}}=\lim_{p\to\infty}\frac{1+\tau
2^{-i}\big{(}1-2^{-m}\big{)}^{-1}\sqrt{p}}{\sqrt{(\frac{\tau^{2}}{3}+1)p+o(p)}}\to\sqrt{\frac{3\tau^{2}}{\tau^{2}+3}}\cdot\frac{1}{2^{i}}$

as $p\to\infty$. Now,

$\displaystyle\sum_{i=1}^{\infty}\rho_{i}^{2}=\frac{3\tau^{2}}{\tau^{2}+3}\sum_{i=1}^{\infty}\frac{1}{4^{i}}=\frac{\tau^{2}}{\tau^{2}+3}.$

Assuming ${\boldsymbol{\mu}}=\boldsymbol{0}$, by Theorem 1, we see that
$T_{p,1}\to\sqrt{\frac{3}{\tau^{2}+3}}\xi_{0}+\sqrt{\frac{3\tau^{2}}{2(\tau^{2}+3)}}\sum_{i=1}^{\infty}\frac{1}{2^{i}}(\xi_{i}^{2}-1)$
in distribution, where $\xi_{0},\xi_{1},\cdots,\xi_{d}$ are i.i.d. $N(0,1).$

## References

  * Anderson (2003) Theodore W Anderson.  _An Introduction to Multivariate Statistical Analysis_.  John Wiley & Sons, 2003.
  * Bai and Saranadasa (1996) Zhidong Bai and Hewa Saranadasa.  Effect of high dimension: by an example of a two sample problem.  _Statistica Sinica_ , pages 311–329, 1996.
  * Bai and Silverstein (2010) Zhidong Bai and Jack W Silverstein.  _Spectral analysis of large dimensional random matrices_ , volume 20.  Springer, 2010.
  * Bao et al. (2012) Zhigang Bao, Guangming Pan, and Wang Zhou.  Tracy-widom law for the extreme eigenvalues of sample correlation matrices.  _Electronic Journal of Probability_ , 17:1–32, 2012.
  * Biswas and Ghosh (2014) Munmun Biswas and Anil K Ghosh.  A nonparametric two-sample test applicable to high dimensional data.  _Journal of Multivariate Analysis_ , 123:160–171, 2014\.
  * Boucheron et al. (2013) Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  _Concentration inequalities: A nonasymptotic theory of independence_.  Oxford university press, 2013.
  * Brockwell et al. (2016) Peter J Brockwell, Peter J Brockwell, Richard A Davis, and Richard A Davis.  _Introduction to time series and forecasting_.  Springer, 2016.
  * Cai and Jiang (2012) T Tony Cai and Tiefeng Jiang.  Phase transition in limiting distributions of coherence of high-dimensional random matrices.  _Journal of Multivariate Analysis_ , 107:24–39, 2012.
  * Cai et al. (2013) T Tony Cai, Jianqing Fan, and Tiefeng Jiang.  Distributions of angles in random packing on spheres.  _Journal of Machine Learning Research_ , 14:1837, 2013.
  * Cai et al. (2014) T Tony Cai, Weidong Liu, and Yin Xia.  Two-sample test of high dimensional means under dependence.  _Journal of the Royal Statistical Society: Series B: Statistical Methodology_ , pages 349–372, 2014.
  * Chakraborty and Chaudhuri (2017) Anirvan Chakraborty and Probal Chaudhuri.  Tests for high-dimensional data based on means, spatial signs and spatial ranks.  _The Annals of Statistics_ , 45(2):771–799, 2017\.
  * Chang et al. (2017) Jinyuan Chang, Chao Zheng, Wen-Xin Zhou, and Wen Zhou.  Simulation-based hypothesis testing of high dimensional means under covariance heterogeneity.  _Biometrics_ , 73(4):1300–1310, 2017.
  * Chen et al. (2011) Lin S Chen, Debashis Paul, Ross L Prentice, and Pei Wang.  A regularized hotelling’s t 2 test for pathway analysis in proteomic studies.  _Journal of the American Statistical Association_ , 106(496):1345–1360, 2011.
  * Chen and Qin (2010) Song Xi Chen and Ying-Li Qin.  A two-sample test for high-dimensional data with applications to gene-set testing.  _The Annals of Statistics_ , 38(2):808–835, 2010\.
  * Chen et al. (2014) Song Xi Chen, Jun Li, and Ping-Shou Zhong.  Two-sample tests for high dimensional means with thresholding and data transformation.  _arXiv preprint arXiv:1410.2848_ , 2014.
  * Chen et al. (2019) Song Xi Chen, Jun Li, and Ping-Shou Zhong.  Two-sample and anova tests for high dimensional means.  _The Annals of Statistics_ , 47(3):1443–1474, 2019.
  * Chow and Teicher (1997) Yuan Shih Chow and Henry Teicher.  Probability theory: Independence, interchangeability, martingales.  1997\.
  * Dempster (1958) Arthur P Dempster.  A high dimensional two sample significance test.  _The Annals of Mathematical Statistics_ , pages 995–1010, 1958.
  * Dempster (1960) Arthur P Dempster.  A significance test for the separation of two highly multivariate small samples.  _Biometrics_ , 16(1):41–50, 1960.
  * Dong et al. (2016) Kai Dong, Herbert Pang, Tiejun Tong, and Marc G Genton.  Shrinkage-based diagonal hotelling’s tests for high-dimensional small sample size data.  _Journal of Multivariate Analysis_ , 143:127–142, 2016\.
  * Durrett (2019) Rick Durrett.  _Probability: theory and examples_ , volume 49.  Cambridge university press, 2019.
  * Fan and Jiang (2019) Jianqing Fan and Tiefeng Jiang.  Largest entries of sample correlation matrices from equi-correlated normal populations.  _The Annals of Probability_ , 47(5):3321–3374, 2019.
  * Fan et al. (2020) Jianqing Fan, Jianhua Guo, and Shurong Zheng.  Estimating number of factors by adjusted eigenvalues thresholding.  _Journal of the American Statistical Association_ , pages 1–10, 2020\.
  * Feng et al. (2016) Long Feng, Changliang Zou, and Zhaojun Wang.  Multivariate-sign-based high-dimensional tests for the two-sample location problem.  _Journal of the American Statistical Association_ , 111(514):721–735, 2016.
  * Gao et al. (2017) Jiti Gao, Xiao Han, Guangming Pan, and Yanrong Yang.  High dimensional correlation matrices: The central limit theorem and its applications.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 79(3):677–693, 2017.
  * Ghosh (1966) B. K. Ghosh.  Asymptotic expansions for the moments of the distribution of correlation coefficient.  _Biometrika_ , 53(1/2):258–262, 1966.
  * Gregory et al. (2015) Karl Bruce Gregory, Raymond J Carroll, Veerabhadran Baladandayuthapani, and Soumendra N Lahiri.  A two-sample test for equality of means in high dimension.  _Journal of the American Statistical Association_ , 110(510):837–849, 2015.
  * Gretton et al. (2012) Arthur Gretton, Dino Sejdinovic, Heiko Strathmann, Sivaraman Balakrishnan, Massimiliano Pontil, Kenji Fukumizu, and Bharath K Sriperumbudur.  Optimal kernel choice for large-scale two-sample tests.  In _Advances in neural information processing systems_ , pages 1205–1213. Citeseer, 2012.
  * Guiard (1986) V Guiard.  A general formula for the central mixed moments of the multivariate normal distribution.  _Statistics: A Journal of Theoretical and Applied Statistics_ , 17(2):279–289, 1986.
  * Guo and Chen (2016) Bin Guo and Song Xi Chen.  Tests for high dimensional generalized linear models.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ , pages 1079–1102, 2016.
  * Horn and Johnson (2012) Roger A Horn and Charles R Johnson.  _Matrix analysis_.  Cambridge university press, 2012.
  * Hotelling (1931) Harold Hotelling.  The generalization of student’s ratio.  _The Annals of Mathematical Statistics_ , 2(3):360–378, 1931.
  * Hu and Bai (2016) Jiang Hu and Zhidong Bai.  A review of 20 years of naive tests of significance for high-dimensional mean vectors and covariance matrices.  _Science China Mathematics_ , 59(12):2281–2300, 2016.
  * Jiang (2004) Tiefeng Jiang.  The limiting distributions of eigenvalues of sample correlation matrices.  _Sankhyā: The Indian Journal of Statistics_ , pages 35–48, 2004\.
  * Jiang (2009) Tiefeng Jiang.  A variance formula related to a quantum conductance problem.  _Physics Letters A_ , 373(25):2117–2121, 2009\.
  * Jiang and Yang (2013) Tiefeng Jiang and Fan Yang.  Central limit theorems for classical likelihood ratio tests for high-dimensional normal distributions.  _The Annals of Statistics_ , 41(4):2029–2074, 2013.
  * Johansen (1980) Søren Johansen.  The Welch-James approximation to the distribution of the residual sum of squares in a weighted linear regression.  _Biometrika_ , 67(1):85–92, 1980.
  * Krishnamoorthy and Yu (2004) K Krishnamoorthy and Jianqi Yu.  Modified nel and van der merwe test for the multivariate behrens–fisher problem.  _Statistics & probability letters_, 66(2):161–169, 2004.
  * Ledoux (2001) Michel Ledoux.  _The concentration of measure phenomenon_.  Number 89. American Mathematical Soc., 2001.
  * Li et al. (2020) Haoran Li, Alexander Aue, Debashis Paul, Jie Peng, and Pei Wang.  An adaptable generalization of hotelling’s $t^{2}$ test in high dimension.  _The Annals of Statistics_ , 48(3):1815–1847, 2020.
  * Lopes et al. (2011) Miles Lopes, Laurent Jacob, and Martin J Wainwright.  A more powerful two-sample test in high dimensions using random projection.  _Advances in Neural Information Processing Systems_ , 24:1206–1214, 2011.
  * Morales-Jimenez et al. (2021) David Morales-Jimenez, Iain M Johnstone, Matthew R McKay, and Jeha Yang.  Asymptotics of eigenstructure of sample correlation matrices for high-dimensional spiked models.  _Statistica Sinica_ , 31(2):571, 2021.
  * Muirhead (1982) Robb J Muirhead.  _Aspects of multivariate statistical theory_.  John Wiley & Sons, 1982.
  * Nel and Van der Merwe (1986) DG Nel and CA Van der Merwe.  A solution to the multivariate behrens-fisher problem.  _Communications in Statistics-Theory and Methods_ , 15(12):3719–3735, 1986.
  * Park and Ayyala (2013) Junyong Park and Deepak Nag Ayyala.  A test for the mean vector in large dimension and small samples.  _Journal of Statistical Planning and Inference_ , 143(5):929–943, 2013.
  * Srivastava (2009) Muni S Srivastava.  A test for the mean vector with fewer observations than the dimension under non-normality.  _Journal of Multivariate Analysis_ , 100(3):518–532, 2009.
  * Srivastava and Du (2008) Muni S Srivastava and Meng Du.  A test for the mean vector with fewer observations than the dimension.  _Journal of Multivariate Analysis_ , 99(3):386–402, 2008.
  * Srivastava and Fujikoshi (2006) Muni S Srivastava and Yasunori Fujikoshi.  Multivariate analysis of variance with fewer observations than the dimension.  _Journal of Multivariate Analysis_ , 97(9):1927–1940, 2006.
  * Srivastava et al. (2013) Muni S Srivastava, Shota Katayama, and Yutaka Kano.  A two sample test in high dimensional data.  _Journal of Multivariate Analysis_ , 114:349–358, 2013\.
  * Srivastava et al. (2016) Radhendushka Srivastava, Ping Li, and David Ruppert.  RAPTT: An exact two-sample test in high dimensions using random projections.  _Journal of Computational and Graphical Statistics_ , 25(3):954–970, 2016.
  * Wang et al. (2015) Lan Wang, Bo Peng, and Runze Li.  A high-dimensional nonparametric multivariate test for mean vector.  _Journal of the American Statistical Association_ , 110(512):1658–1669, 2015.
  * Welch (1947) Bernard L Welch.  The generalization of ‘student’s’problem when several different population varlances are involved.  _Biometrika_ , 34(1-2):28–35, 1947.
  * Wu et al. (2006) Yujun Wu, Marc G Genton, and Leonard A Stefanski.  A multivariate two-sample mean test for small sample size and missing data.  _Biometrics_ , 62(3):877–885, 2006.
  * Xiao and Zhou (2010) Han Xiao and Wang Zhou.  Almost sure limit of the smallest eigenvalue of some sample correlation matrices.  _Journal of Theoretical Probability_ , 23(1):1–20, 2010.
  * Xu et al. (2016) Gongjun Xu, Lifeng Lin, Peng Wei, and Wei Pan.  An adaptive two-sample test for high-dimensional means.  _Biometrika_ , 103(3):609–624, 2016.
  * Xue and Yao (2020) Kaijie Xue and Fang Yao.  Distribution and correlation-free two-sample test of high-dimensional means.  _The Annals of Statistics_ , 48(3):1304–1328, 2020.
  * Yao (1965) Ying Yao.  An approximate degrees of freedom solution to the multivariate behrens fisher problem.  _Biometrika_ , 52(1/2):139–147, 1965.
  * Zhang et al. (2017) Jin-Ting Zhang, Jia Guo, and Bu Zhou.  Linear hypothesis testing in high-dimensional one-way MANOVA.  _Journal of Multivariate Analysis_ , 155:200–216, 2017\.
  * Zhang et al. (2020a) Jin-Ting Zhang, Jia Guo, Bu Zhou, and Ming-Yen Cheng.  A simple two-sample test in high dimensions based on l 2-norm.  _Journal of the American Statistical Association_ , 115(530):1011–1027, 2020a.
  * Zhang et al. (2020b) Liang Zhang, Tianming Zhu, and Jin-Ting Zhang.  A simple scale-invariant two-sample test for high-dimensional data.  _Econometrics and Statistics_ , 14:131–144, 2020b.
  * Zhong et al. (2013) Ping-Shou Zhong, Song Xi Chen, and Minya Xu.  Tests alternative to higher criticism for high-dimensional means under sparsity and column-wise dependence.  _The Annals of Statistics_ , 41(6):2820–2851, 2013.