Exact large deviation function of spin current for the one dimensional XX spin chain with domain wall initial condition

TL;DR

This paper derives an exact large deviation function for spin current fluctuations in a 1D XX spin chain with domain wall initial conditions, using determinant formulas and Coulomb gas methods, validated by DMRG.

Contribution

It provides the first exact analytical expression for the large deviation function of spin current in this model, advancing understanding of quantum spin transport.

Findings

Exact large deviation function derived

Determinant representation with Bessel kernel obtained

Results confirmed by DMRG simulations

Abstract

We investigate the fluctuations of the spin current for the one dimensional XX spin chain starting from the domain wall initial condition. The generating function of the current is shown to be written as a determinant with the Bessel kernel. An exact analytical expression for the large deviation function is obtained by applying the Coulomb gas method. Our results are also compared with DMRG calculations.

Figures

Equations

$$⟨N(t)⟩\simeq \frac{t}{\pi },{⟨N(t{)}^2⟩}_c=⟨{\left[N(t)-⟨N(t)⟩\right]}^2⟩\simeq \frac{1}{2{\pi }^2}(\log t+C),$$

$$H_{XX}=-J\sum _{m=-\infty }^{\infty }(s_m^x{s}_{m+1}^x+s_m^y{s}_{m+1}^y).$$

$$s_m^i={\sigma }_m^i/2,m\in \mathbb{Z},i=x,y,z.$$

$$|\mathrm{DW}⟩=|\cdots \uparrow \uparrow \uparrow \downarrow \downarrow \downarrow \cdots ⟩.$$

$$N(t):=\sum _{m=\alpha +1}^{\infty }\left[s_m^z(t)+\frac{1}{2}\right].$$

$$M(t):=\sum _{m=\alpha +1}^{\infty }{s}_m^z(t).$$

$$\frac{d}{dt}{s}_m^z(t)=i\left[H_{XX},s_m^z\right]=j_{m-1}-j_m,$$

$$j_m=J\left(s_{m+1}^x{s}_m^y-s_{m+1}^y{s}_m^x\right).$$

$$N(t)={\int }_0^t{j}_{\alpha }(\tau )d\tau .$$

$$\underset{m\to \pm \infty }{\lim }{j}_m(t)=0.$$

$$\chi (\lambda ,t):=⟨e^{\lambda N(t)}⟩.$$

$$⟨e^{\lambda N(t)}⟩=\sum _{n=0}^{\infty }{e}^{\lambda n}\mathrm{P}[N(t)=n].$$

$$\mathrm{P}[N(t)=n]=\sum _{\begin{array}{c}{s}_m^z=\pm 1/2,m\in \mathbb{Z}\\ \\ N:={\sum }_{m=\alpha +1}^{\infty }(s_m^z+1/2)=n\end{array}}⟨\mathrm{DW}⟩e^{it{H}_{XX}}\dots ,s_m^z,\dots ⟨\dots ,s_m^z,\dots ⟩e^{-it{H}_{XX}}\mathrm{DW},$$

$$\mathrm{P}[N(t)=aJt]\sim e^{-t^2{J}^2\Psi (a)}$$

$$F(r)=\frac{1}{4}(1-r)+a(r)\sqrt{r}\left[K\left(1-\frac{1}{r}\right)-E\left(1-\frac{1}{r}\right)\right]$$

$$a(r)=\{\begin{array}{ll}\frac{1}{\pi }\sqrt{1-r}\left[E\left(\frac{r}{r-1}\right)-K\left(\frac{r}{r-1}\right)\right],& \mathrm{for}0\le r<1,\\ \frac{1}{\pi }\sqrt{r}E\left(\frac{1}{r}\right),& \mathrm{for}1<r.\end{array}$$

$$s_m^{\pm }:=s_m^x\pm i{s}_m^y.$$

$$c_j^{†}$$

$$c_j$$

$$H_{XX}=-\frac{J}{2}\sum _{m=-\infty }^{\infty }\left(c_m^{†}{c}_{m+1}+c_{m+1}^{†}{c}_m\right)=\sum _{i,j=-\infty }^{\infty }{c}_i^{†}{H}_{ij}{c}_j,$$

$$H_{ij}=-\frac{J}{2}({\delta }_{i,j-1}+{\delta }_{i,j+1}).$$

$$c_m(t)=\sum _{n=-\infty }^{\infty }{U}_{mn}(t,0)c_n=\sum _{n=-\infty }^{\infty }{i}^{m-n}{J}_{m-n}(t)c_n,$$

$$J_{\nu }(z)=\sum _{k=0}^{\infty }\frac{(-1{)}^k}{k!\Gamma (k+\nu +1)}{\left(\frac{z}{2}\right)}^{2k+\nu }.$$

$$|\mathrm{DW}⟩=\prod _{m=-\infty }^0{c}_m^{†}|0⟩,$$

$$N(t)=\sum _{m=\alpha +1}^{\infty }{c}_m^{†}(t)c_m(t).$$

$$j_m=-\frac{1}{2i}\left(c_{m+1}^{†}{c}_m-c_m^{†}{c}_{m+1}\right).$$

$$\chi (\lambda ,t)=\det {\left[1+(e^{\lambda }-1)C(t)\right]}_{{\ell }^2({\mathbb{Z}}_{>\alpha })},$$

$$C_{mn}(t)=⟨c_m^{†}(t)c_n(t)⟩=\frac{1}{2\pi }{\int }_{-k_F}^{k_F}{\phi }_m^{\ast }(k;t){\phi }_n(k;t)dk.$$

$${\phi }_m(k;t)=\sum _{j=-\infty }^0{i}^{m-j}{J}_{m-j}(t)e^{ikj},$$

$$\chi (\lambda ,t)=\det {\left[1+(e^{\lambda }-1)K_{\mathrm{DBes}}(t)\right]}_{{\ell }^2({\mathbb{Z}}_{>\alpha })},$$

$$\begin{array}{rl}(K_{\mathrm{DBes}}{)}_{mn}(t)& =\sum _{k=0}^{\infty }{J}_{m+k}(t)J_{n+k}(t)\\ & =\frac{t}{2(m-n)}\left[J_{m-1}(t)J_n(t)-J_m(t)J_{n-1}(t)\right],m,n\in \mathbb{Z}.\end{array}$$

Full text

1 Introduction

Recently, quantum dynamics of many-body systems attract interests of many researchers, and there has been impressive progress in both experiments and theories [1, 2]. It is generally hard to study dynamical properties of quantum many-body systems analytically. But in the case of one dimension, some exact results have been obtained. They have not only provided us useful theoretical insights about quantum dynamics and thermalization but are also in many cases experimentally relevant.

For example, the famous experiment about a quantum version of Newton’s cradle with one-dimensional Bose gases [3] proposed an example of systems which do not thermalize due to their integrability. This experiment inspired the introduction of the generalized Gibbs ensemble (GGE) [4] as a state quantum integrable systems tend to relax to.

More recently, a hydrodynamic formulation called the Generalized Hydrodynamics (GHD) was proposed in [5, 6] and has been successfully applied to various quantum integrable systems such as XXZ Heisenberg chain, Lieb-Liniger model, for calculating their density, current profile and so on.

One of the most standard setups to study systems out of equilibrium is the so-called partitioning protocol in an isolated quantum system, in which initially two regions of a system are prepared in macroscopically different states, at $t=0$ they are connected abruptly, and then one studies the time evolution of the whole system.

In the context of spin chains, the simplest partitioning protocol is to take the domain wall initial state, in which spins in two regions are prepared in the opposite directions, and to study the spin dynamics after that. This inhomogeneous magnetization profile before the quantum quench takes place helped to reveal several non-equilibrium behaviors in different models (see for instance [6, 7, 8]). Domain wall initial state having a step-like spin configuration is pretty easy to make in real materials because we only have to apply the external strong magnetic field into each system.

Many results have been obtained for this setting. For example, for the case of the XX spin chain, the magnetization density profile has already been found in [9]. It is also reproduced by GHD in [6]. Time-dependent DMRG were also applied to spin chains’ dynamics initiated from this condition (for instance [10, 11, 12]) including different anisotropy parameter $\Delta$ of XXZ spin chain. But most results so far have been for average behaviors of the systems. Fluctuation properties have much less been addressed, though very recently there was important progress in the context of GHD [13].

In this paper, we focus on the XX chain and integrated spin current $N(t)$ which is given by the total number of up spins which moved from the left to the right subsystem for enough large time period $t$ under the domain wall initial condition at absolute zero temperature.

The large time behavior of the mean and the variance of this quantity have been already known in the previous works [9, 14] :

[TABLE]

where the constant $C$ is given as a sum of a few constants and definite integrals and numerically $C=2.963510026\dots$. The notation $\simeq$ means that both sides are asymptotically equal. Our aim in this paper is to study the full distribution of $N(t)$ for large $t$, in particular to give an explicit formula for the large deviation function of $N(t)$.

The XX model is well-known to be equivalent to a free fermionic system, for which the seminal Levitov-Lesovik formula [15] is available. The formula, in particular, the long-time approximated form of it (see eq. (9) in [15]) was established in 1993 by using the scattering matrix and has a wide application, but for our setting with zero temperature and perfect transmission starting from the domain wall state, a naive application of the Levitov-Lesovik formula would lead to the Dirac delta shape distribution and does not create the variance in eq. (1.1) which increases logarithmically in time as pointed out in [16, 17].

In this paper, we do not rely on the long-time approximated Levitov-Lesovik formula and study the statistics of $N(t)$, in particular the large deviation function (a.k.a. the rate function), by exact calculations. It contains far more information than the average and the variance in eq. (1.1). As is well known, the moment generating function and its large deviation function are mutually convertible through the Legendre transform [18]. The large deviation in non-equilibrium steady state (NESS) has been investigated for open systems with boundary driven condition [19, 20, 21], however, we directly treat the isolated system from the quench without assuming NESS.

Our arguments are based on the determinant formula for the generating function of $N(t)$ found in [22]. In [22], the authors mainly studied the quantum propagating front which is the edge of the melting up spins moving to the other subsystem. They found that the statistics at the quantum front is described by the Airy kernel [23] asymptotically and thus showed that the rightmost up spin’s existence probability is given by the Tracy-Widom distribution [24]. In this paper, we will rewrite the generating function in terms of the Bessel kernel, and apply the Coulomb gas technique to study the large deviation properties. We also found that time-dependent DMRG method can be extended to calculate the rate function through moment generating function by applying matrix product operators having the dependence on $\lambda$ successively.

The rest of the paper is organized as follows. In section 2, we introduce the model and some notations and also present our main result. In section 3, we explain the mapping from our dynamical problem into a static problem in terms of random matrix theory. Section 4 is to be devoted to a review of the role of the Bessel kernel [23] in random matrix theory which appeared in the previous section. In section 5, we will explain the techniques to study the large deviation properties of the Wishart matrix and apply it to our setting to find the exact formula for the large deviation function of $N(t)$. In section 6, by expanding the large deviation function around the average, we reproduce the variance which increases logarithmically in time. In section 7, we will apply the time-dependent DMRG calculation to the evaluation of the large deviation function. Section 8 is the conclusion.

2 The model and results

We consider the dynamics of the XX model on the one-dimensional lattice. A site at which a spin is located is designated by an integer. The Hamiltonian is given by

[TABLE]

Here, $J$ is the coupling constant and $s_{m}^{i}$ denotes the spin operator at site $m$ defined as half of the Pauli operator,

[TABLE]

As an initial condition at $t=0$, we employ the domain wall initial condition, in which for the left half of the lattice $m\leq 0$, all the spins are set to be up spins, while for the right half of the lattice $m\geq 1$, all the spins are set to be down spins (see Fig.1). We denote this state by

[TABLE]

For this model, the quantity we are interested in is the integrated spin current, which is defined as

[TABLE]

Here, the time dependence of an operator is from the Heisenberg picture. Without loss of generality, we can only consider the case $\alpha\geq 0$ because, as long as the domain wall initial condition is considered, we can use the symmetry of the system between up spins and down spins. Counting up spins, for example, in the region $m\geq-1$ or $\alpha=-2$ is equivalent to counting down spins in the region $m\leq-2$.

To see the meaning of this observable $N(t)$, let us consider the total magnetization operator in the area $m>\alpha$ as

[TABLE]

For our domain wall initial condition, this quantity is infinite, but the difference $M(t)-M(0)=N(t)$, which can be interpreted as the difference of the total magnetization in the area $m>\alpha$ between [math] and $t$ is well-defined. At time $t=0$, we can, of course, confirm that $N(0)=0$ from the definition of the initial condition.

When $\alpha=0$, we are supposed to count all the transported up spins in the right half of the chain $m\geq 1$ at $t$. Note that the dependence on $\alpha$ is abbreviated in the above notations of $N(t)$ and $M(t)$. From the time derivative of the local magnetization $s_{m}^{z}$,

[TABLE]

we can also consider the instantaneous spin current [25],

[TABLE]

Here and in the following, we set the reduced Planck constant $\hbar$ to be unity for simplicity.

In terms of the time-dependent spin current $j_{\alpha}(t)$, the integrated current $N(t)$ given in eq. (2.4) can also be written as

[TABLE]

This can be immediately derived by taking the sum of the local conservation law in eq. (2.6) from $m=\alpha+1$ to $\infty$ and integrating with regard to time from [math] to $t$. Here, we implicitly assumed the situation that the effect of the other boundary current away from the center of the chain is zero,

[TABLE]

To consider the statistics of the total number of transmitted up spins $N(t)$, we define the moment generating function as

[TABLE]

The moment generating function can be expanded as

[TABLE]

Here $\mathrm{P}[N(t)=n]$, defined as

[TABLE]

has the meaning of the probability that the total number of transported up spins $N(t)$ at $t$ is $n$. The letter $N(t)$ is slightly abused because here it is treated like a random variable, not an operator as defined in eq. (2.4), but there should be no confusion. Here we remark that the probability $\mathrm{P}[N(t)=n]$ is, in fact, independent of the sign of the coupling constant $J$, because changing the sign of $J$ is equivalent to reversing the time, the Hamiltonian $H_{XX}$ is real symmetric and the domain wall initial state has only real components. In other words, the distribution of the spin current $N(t)$ is actually the same for ferromagnetic and antiferromagnetic chains. Therefore in the following discussions, we set $J>0$ for convenience.

The probability $\mathrm{P}[N(t)=n]$ can be recovered from the moment generating function $\chi(\lambda,t)$. After the replacement of $\lambda$ by $i\lambda$, the Fourier coefficients of the characteristic function $\chi(i\lambda,t)$ give the probability $\mathrm{P}[N(t)=n]$ for each $n$.

In this paper, our aim is to calculate the large deviation function of the spin current exactly under the specific setup explained above. Our main result is the following. For the case of $\alpha=0$, the full distribution function of the integrated spin current behaves for long time as

[TABLE]

where the large deviation function can be described in a parametric way $\Psi(a)=F(r(a))$ as

[TABLE]

with

[TABLE]

and $r(a)$ is the inverse function of $a(r)$. Here $E$ and $K$ are the complete elliptic integrals of the first and second kind, see eq. (5.16) for their definitions, and $\sim$ means that both hand sides are equal after their logarithm taken, divided by $t^{2}$ and performing the limit $t\rightarrow\infty$.

To solve our problem, we use a well-known mapping from the XX spin chain to a free fermion. Let us introduce the raising and lowering operators acting on each Hilbert space $\mathbb{C}^{2}$ at site $m$,

[TABLE]

Performing the Jordan-Wigner transformation,

[TABLE]

one can map a spin chain into a fermionic chain. For the XX spin chain, one obtains the free fermionic Hamiltonian :

[TABLE]

where the single particle Hamiltonian is just given by

[TABLE]

In the following, the ferromagnetic coupling constant $J$ is set to unity. Using the eq. (2.19), the time evolution of fields are as follows [9],

[TABLE]

where $U$ is the time evolution operator whose matrix element is given in eq. (A.6). Here, $J_{n}(t)$ is the Bessel function of the first kind of order $n$ [26], generally defined by

[TABLE]

In the free fermion language, the state corresponding to the domain wall initial condition in eq. (2.3) is written as

[TABLE]

where $\ket{\mathrm{0}}$ means the vacuum state which does not contain any fermion. The total magnetization flowing into the region $m>\alpha$ in eq. (2.4) equals the number of the transmitted charges,

[TABLE]

The spin current in eq. (2.7) also has a fermionic representation,

[TABLE]

3 Full counting statistics and random matrix analogy

In the studies of quantum transport, especially in the field of mesoscopic physics, obtaining the full distribution of the transferred charges is of great interest [27]. This is referred to as full counting statistics (FCS) in literature. The moment generating function $\chi(\lambda,t)$ of the number of transferred fermions $N(t)$ in this context has been investigated in a series of studies [15, 28, 29, 16, 17] and many others. For the case of free fermion at zero temperature $T=0$, it is known to be written in a form of determinant,

[TABLE]

where $C(t)$ is an infinite dimensional matrix. When there is no fermion on the right chain, its element is of the form

[TABLE]

Here $\varphi_{m}(k;t)$ is the wave function at the position $m$ and at time $t$. Letter $k$ denotes the wave number in the single particle system and $k_{F}$ is its Fermi wave number. In our specific case, considering the unitary time evolution in eq. (2.21), the explicit form of $\varphi_{m}(k;t)$ is given by

[TABLE]

with $-\pi<k\leq\pi$ (i.e. $k_{F}=\pi$). After some calculations, as shown in [22], the moment generating function can be written as

[TABLE]

where the discrete Bessel kernel [30, 31] is given by

[TABLE]

The discrete Bessel kernel also appears in the studies of the polynuclear growth (PNG) model (e.g. [32]) and the longest increasing subsequence (LIS) in random permutation [33]. The determinant that appears in eq. (3.1) or eq. (3.4) is the Fredholm determinant, defined as

[TABLE]

By comparing the Taylor expansion at the point $z=1$ of this Fredholm determinant assuming $z=1-e^{\lambda}$, with the identity in eq. (2.11), we can see that

[TABLE]

In [22], it is also shown that the right hand side of eq. (3.7) has a meaning that exactly $n$ particles lie in the interval $m>\alpha$ and the analogous structures between FCS and Random Matrix Theory (RMT) was pointed out. For general instruction about RMT, see for instance [34, 35].

Now an important observation for our analysis in this paper is that the Fredholm determinant with the discrete Bessel kernel in eq. (3.4) can also be written as the Fredholm determinant with the celebrated Bessel kernel [23] for any parameter $z\in\mathbb{C}$,

[TABLE]

where

[TABLE]

Our proof is based on a direct calculation of the trace which appears in the expansion of the logarithm of the determinants, see Appendix A.

The relation in eq. (3.8) casts a fermion counting problem in the infinite and discrete spatial interval with fixed $t$ into the same problem but in the finite and continuous interval. By considering $[0,t^{2}]$ as if it were a spatial interval, the latter situation can be taken as the non-interacting fermions in a confining potential and there have been considerable works on the subject (for example, see [36]). Let $E(n;[0,t^{2}])$ denote the probability that $n$ fermions exist in the interval $[0,t^{2}]$ and write it with the Bessel kernel,

[TABLE]

By observing that the expansion of both sides of eq. (3.8) in terms of $z$ at $z=1$ are given by eq. (3.7) and eq. (3.10), we have

[TABLE]

In fact, we will study the large deviation of the left hand side of eq. (3.11) by analyzing that of the right hand side of eq. (3.11).

When $\alpha=0$, from eq. (3.7) we readily notice that the special case $n=0$ in eq. (3.7), is nothing but the return probability $R(t)$ (the quantum Loschmidt echo), which have been studied in [32, 37, 38, 39, 40, 41]. It is defined as the modulus square of the overlap of the initial state and the final state, which is taken as the same as the initial state :

[TABLE]

For the case of domain wall initial state, and when we observe the up spins in the whole right subsystem $m>\alpha=0$ at time $t$, we have

[TABLE]

All representations here are totally equivalent. The meaning of the first equality is clear. When the final state at time $t$ is completely the same as the initial state, there should be no particle in the right half of the chain $m>\alpha=0$.

The rightmost expression for the return probability as the Fredholm determinant is a special case ($z=1$) of the eq. (3.8)444The $z=1$ case can also be found by combining (8.96) in [35] and (4.3) in [42], as pointed out by P. Forrester. and seems to have close connection with the expression appearing in [40] by considering the calculation of the partition function using transfer matrix formalism in Euclidian time and some well-known determinant formula related to six-vertex model in two dimension with domain wall boundary conditions. But we emphasize that the relation in eq. (3.8), which holds for arbitrary $z$, and whose direct proof is given in Appendix A, will be crucial in our derivation of the large deviation function.

Actually, in the case of the XX spin chain with domain wall initial condition, a simple result for the return probability was obtained in [38, 43]

[TABLE]

for any time $t$. By taking $t$ large, we can calculate $\Psi(a=0)=1/4$. This also makes us expect that the tail of the distribution of the integrated spin current decays as Gaussian as time passes as we will see later.

We can also evaluate the gradient at the leftmost edge of the large deviation tail. Let us focus on the symmetric case $\alpha=0$ and write the following asymptotic formula for $r(n,t^{2})$ which has been known for fixed $n$ and large $t$, see eq. (1.28) in [44],

[TABLE]

[TABLE]

where $G(x)$ is the Barnes’ G-Function [26]. We now consider the large $n$ behavior of this. With the help of the asymptotic expansion of the Barnes’ G-Function $G(x)$ for large argument $z$ [26] (see chapter 5.17),

[TABLE]

where $A=1.282427\dots$ is the Glaisher-Kinkelin constant, and using Stirling’s approximation, we find

[TABLE]

for large $t$ and $n$, but with the condition $n/t\ll 1$, because we considered the large $n$ limit after taking large $t$ limit. When we scale the variable $n$ as $at$ due to the fact that mean current $\braket{N(t)}\simeq t/\pi$ is ballistic, we obtain

[TABLE]

for large $t$. From that asymptotics, we can read that the gradient near the origin would be given as $\Psi^{\prime}(0)=-2$.

We remark that, even though our main focus of this work is for a finite $\alpha$, in particular the $\alpha=0$ case, some formulas in this paper should also be useful for studying the fluctuation of the quantum front (by setting the position of measuring current as $\alpha=t+2^{-1/3}t^{1/3}s$) for large $t$ [45, 22, 38, 46]. Here $s$ is a local coordinate describing the wavefront. In [47], it was known that the regime described by the Bessel kernel with the interval $[0,Q_{\alpha}(s)]$, where $Q_{\alpha}(s)=(\alpha-2^{-1/3}\alpha^{1/3}s)^{2}$, tends to be that described by the Airy kernel [23] with the interval $[s,\infty)$, in the large $\alpha$ limit. By identifying $Q_{\alpha}(s)$ as $t^{2}$ and considering large $\alpha$ limit, we can focus on the statistics of up spins in the front regime.

When we evaluate the asymptotic behavior of the number variance of integrated current $N(t)$ in Appendix C, it is useful to work with the variable $\sqrt{x}\rightarrow x$ in eq. (3.9), which leads to a transformed Bessel kernel,

[TABLE]

In this case, the kernel $\tilde{K}_{\mathrm{Bes}}^{(\alpha)}$ in eq. (3.20) acts on $L^{2}[0,t]$. Thus, expanding the Fredholm determinant with regard to $\tilde{K}_{\mathrm{Bes}}^{(\alpha)}$ with the use of a continuous version of eq. (3.6), we have

[TABLE]

4 Bessel kernel in RMT

Here, we briefly review how the Bessel kernel emerges in RMT. Let $X$ be an $M\times N~{}(M\geq N)$ Gaussian matrix whose elements are complex and follow the Gaussian distribution having unit variance but zero mean. Then, the $N\times N$ matrix

[TABLE]

is called the Wishart matrix which is originally from [48]. The joint eigenvalue probability density function of the Wishart matrix is given by [49]

[TABLE]

where $\alpha=M-N$. Note that all the eigenvalues are positive $x_{i}\geq 0$.

If we see the eigenvalues as the positions of identical fermionic particles, the joint eigenvalue probability density function can be represented by $N$ point function of fermions as

[TABLE]

where the kernel $K_{N}(x,y)$ can be constructed in terms of wave functions $\{\phi_{k}(x)\}_{k=0}^{N-1}$,

[TABLE]

Since the weight function for this case is $x^{\alpha}e^{-x}$, the wave function can be written in terms of the generalized Laguerre polynomials $L_{k}^{(\alpha)}(x)$ as

[TABLE]

With the help of the Christoffel Darboux formula [50], we have

[TABLE]

It is worth checking that the one point function after a proper scaling [51] gives the Marčenko and Pastur distribution [52]

[TABLE]

When we study behaviors of eigenvalues near the hard edge at $x\simeq 0$, the Bessel kernel appears after the scaling [23],

[TABLE]

In order to get this kernel, the following large $N$ asymptotics formula of the generalized Laguerre polynomials $L_{k}^{(\alpha)}(x)$ [50] is helpful,

[TABLE]

Here, a shorthand notation $N^{\prime}=N+(\alpha+1)/2$ is used.

Thus, we can see that the Bessel kernel describes the statistics of eigenvalues at the hard edge where the smallest eigenvalue of $N\times N$ Wishart matrix is located. To evaluate the probability measure $E(n;[0,t^{2}])$ in an exponential form, we consider the corresponding probability measure $E_{N}(n;[0,t^{2}/4N])$ that $n$ eigenvalues lie in the interval $[0,t^{2}/4N]$ for a large but finite $N$ Wishart ensemble. The change of the interval as compared to eq. (3.10) follows from the scaling of the argument appearing in eq. (4.8).

For large $N$ with finite and fixed $\alpha$, the contribution from the one body potential $\alpha\log x$ can be ignored compared to another one body potential $x/2$. Therefore, we only consider the case $\alpha=0$ from now on. Going back to the joint eigenvalue probability density function in eq. (4.2) with $\alpha=0$ of the Wishart matrix, let us write this probability in an exponential form

[TABLE]

This can be interpreted as the Boltzmann factor and the exponent plays a role of the total energy of the $N$ particle system

[TABLE]

The second term implies that this energy is the energy of gases having two dimensional Coulomb interaction as an analogy.

5 Number statistics via coulomb gas method

Because of the analogy mentioned at the end of the previous section, the Coulomb gas method [53, 35] is widely used to determine certain properties of fermions confined in some potential at zero temperature (see, e.g. [54] and references therein). In [55, 56], large deviation properties of the Wishart matrix were studied with the method, with the normalization such that a scaling $x\rightarrow Nx$ is performed in eq. (4.2). The authors in [55, 56] showed that the probability $E_{N}(n;I)$, that there are $n$ eigenvalues in an arbitrary interval $I(\subseteq[0,\infty))$ for $N\times N$ Wishart matrix with this normalization, satisfies the following large deviation property in the large $N$ limit,

[TABLE]

We will utilize their result to our problem. Notice that, with the change of normalization of the Wishart matrix mentioned above, the discussions in the previous section tells us that the probability (3.10), associated with the Bessel kernel, is recovered in the $N\to\infty$ limit as

[TABLE]

Now we briefly review description of the large deviation function $\psi_{I}$ from [55, 56], for the case of the interval $I=[0,t^{2}/4N^{2}]$ of our interest. Given the interval, and the number of eigenvalues or the fermions $n$ in the interval $I$, the most probable distribution of the spectrum density $\rho^{\ast}(x)$ is determined as a result of minimization problem of the free energy under the constraint. Here the density $\rho^{\ast}(x)$ should be taken so that it satisfies the normalization condition :

[TABLE]

Since $\rho^{\ast}(x)$ which minimizes the total energy in the Boltzmann factor corresponding to the expression in (4.2) gives the main contribution to the probability distribution, finding such a density $\rho^{\ast}(x)$ is a crucial problem. Indeed, $\rho^{\ast}(x)$ exists and the following statement is known. The large deviation function $\psi_{I}(\kappa)$ for large $N$ is given as [55, 56]

[TABLE]

Here $\mu$ and $\eta$ in our case are given by

[TABLE]

where $\lambda_{\pm}$ denotes the two spectral edges of the four, as explained later. In fact, the supports of $\rho^{\ast}(x)$ will be found to be

[TABLE]

Additionally, the resolvent $G(x)$ is defined by the Stieltjes transformation of the density $\rho^{\ast}(x)$ as

[TABLE]

In fact, $G(x)$ satisfies the self-consistent and quadratic equation coming from the saddle point equation with regards to $\rho^{\ast}(x)$. The two solutions are given as

[TABLE]

The sign between the first and the second terms is to be determined by the density $\rho^{\ast}(x)$. Explicitly, it is given in eq. (B.5) and in eq. (B.6). In this resolvent, $\lambda_{+}>\lambda_{-}$ are the roots of the numerator of the fraction inside the square root. From the definition of the resolvent in eq. (5.8), it should decay like $G(z)\simeq 1/z$ for large $z$ as the density $\rho^{\ast}(x)$ is normalized as in eq. (5.3). Therefore, by eq. (5.9), we have

[TABLE]

We can extract the spectral density $\rho^{\ast}(x)$ from the imaginary part of this resolvent as

[TABLE]

Now we can identify the role of parameters $\lambda_{\pm}$, and $t^{2}/4N^{2}$. The spectral density $\rho^{\ast}(x)$ has two compact supports with four edges. The leftmost endpoint $x=0$ and the rightmost endpoint $x=\lambda_{+}$ are always unchanged for large $N$, while the order of the rest endpoints $x=\lambda_{-}$ and $x=t^{2}/4N^{2}$ could change depending on the value of $\lambda_{-}$.

The condition that this spectral density contains $n$ eigenvalues in the interval $I=[0,t^{2}/4N^{2}]$ is reflected in another constraint,

[TABLE]

The two constraints in eq. (5.10) and eq. (5.12) determine the two parameters $\lambda_{\pm}$ uniquely as functions of $n$, $t$ and $N$. Consequently these constraints implies that when $n<t/\pi$, we can see $\lambda_{-}<t^{2}/4N^{2}$ while when $n>t/\pi$, we can see $t^{2}/4N^{2}<\lambda_{-}$. In order to focus on to the hard edge regime $x\simeq 0$, we need to take the large $N$ limit as explained in the previous section.

Before we tackle the large deviation, let us see the mean value $\braket{N(t)}$ with eq. (5.12). For large $N$, the interval $I$ around $x\simeq 0$ having the order $\mathcal{O}(1/N^{2})$ shrinks to zero and the Marčenko and Pastur distribution [52] in $I$ diverges as $\mathcal{O}(N)$ as $N$ increases. Since the spectral density contains the total number of eigenvalues $N$ inside its support, multiplying $\rho^{\ast}(x)$ by $N$, we can extract a meaningful and finite quantity. It turns out to be

[TABLE]

The result is consistent with the previous result shown in eq. (1.1).

Since we are considering the deviation from the mean value $\braket{N(t)}\simeq t/\pi$, the number of transferred up spins $n$ should have the same order as the mean. Thus, it is reasonable to expect that the parameter $\lambda_{-}$ which adjusts $n$ and $t^{2}/4N^{2}$ have the same order. For this reason, we could introduce the ratio $r$ that satisfies

[TABLE]

When $r=1$, the spectral distribution $\rho^{\ast}$ becomes the Marčenko and Pastur distribution [52], since $\lambda_{+}=4$ and $\lambda_{-}=t^{2}/4N^{2}$.

Next, we will analyze the large deviation. As a preparation, we first expand the large deviation function $\psi_{I}(\kappa)$ as a series in $1/N$ with $t$ fixed. Details of the calculation will be given in Appendix B. As a result of expansion, the terms having the order $\mathcal{O}(1)$ and $\mathcal{O}(1/N)$ do not remain. We find

[TABLE]

where $K(k^{2})$ and $E(k^{2})$ are the complete elliptic integrals of the first and second kind [26] respectively. The definitions are as follows

[TABLE]

By definition, $K(0)=E(0)=\pi/2$ follows immediately. Furthermore, we can also check the asymptotic relation such as $K(-k^{2})\simeq(\log 4k)/k$ and $E(-k^{2})\simeq k$ for large $k>0$ from [26, 57]. The subleading terms $\mathcal{O}(1/N^{3})$ are not important because those should be discarded after the limit $N\rightarrow\infty$.

If we set $n=at$ with a new parameter $a$ in (5.15), the leading order terms are proportional to $t^{2}/N^{2}$. Multiplying $2N^{2}$ in front of the large deviation function $\psi_{I}(\kappa)$ as we have seen in (5.1) for large $N$, we get terms proportional to $t^{2}$ in the exponent. Though the expansion in $1/N$ in Appendix B is performed for a fixed $t$, one should consider large $t$ limit because the large deviation results in [55, 56] were obtained for large $n(=at)$. Namely, we should consider the limits of large $t$ and $N$ but with the condition $t\ll N$. In this limit, one can expect that, although there could be some subleading terms, the large deviation behaviors for large $t$ are captured by the terms proportional to $t^{2}$ which were found above.

Recalling eqs. (3.11) and (5.2), and substituting (5.15) into (5.1), we finally obtain

[TABLE]

where the desired large deviation function $\Psi(a)=F(r(a))$ is found as

[TABLE]

with the following condition from eq. (B.20) and eq. (B.21) :

[TABLE]

which is derived by considering large $N$ behaviors of the integral in eq. (5.12), see Appendix B. The shapes of $F(r)$ and $a(r)$ are shown in Fig. 3 and Fig. 3. The large deviation function $\Psi(a)$ is also shown as the solid curve in Fig. 4 below (together with the result by DMRG explained in section 7). The graph is slightly asymmetric with respect to the refection at $a=1/\pi$.

Restoring the ferromagnetic coupling constant in eq. (5.17), we obtain eq. (2.13), which is our main result in the paper.

As one sees in Fig. 3, as $a$ increases, $r$ also tends to increase. When $a$ takes a values from [math] to $1/\pi$, the parameter $r$ varies from [math] to $1$ and when $a$ takes a value more than $1/\pi$, $r$ also varies from $1$ to $\infty$. The function $a(r)$ is always invertible as $r(a)$.

We can check that the minimum of $\Psi(a)$ is attained at $a=1/\pi$ and this value is zero. First, the function $F(r)$ has indeed the minimum at $r=1$. Using the chain rule, the derivative of $\Psi(a)$ with regard to $a$ is

[TABLE]

Here, the derivative of $a$ with regards to $r$ is given by

[TABLE]

This derivative takes a positive value for all $r>0$. Therefore, the sign of the sum of the first and the third term in eq. (5.20) should be determined by the following terms,

[TABLE]

Substituting $a(r)$ in eq. (5.19) tells us that the terms above are monotonically increasing function in terms of $r$. Noticing that

[TABLE]

the sum of the first and third term in eq. (5.20) is negative for $0\leq r<1$ and positive for $r>1$. We can also confirm that the second term in eq. (5.20) takes negative value for $0\leq r<1$ and positive value for $r>1$ from the behavior of the elliptic integrals $E$ and $K$. Hence, the large deviation function $F(r)$ has an unique global minimum at $r=1$. Recalling that $r=1$ corresponds to $a=1/\pi$ from eq. (5.19), we can confirm that

[TABLE]

The Gaussian fluctuation around $a=1/\pi$, will be studied in the next section.

In section 3, we estimated that the large deviation function should behave as $\Psi(0)=1/4$ and $\Psi^{\prime}(0)=-2$ at the leftmost endpoint. Now, we can also check the correctness of the estimates by using the explicit representation of large deviation function derived above. In evaluating those, we can utilize the asymptotic behavior $K(-k^{2})\simeq(\log 4k)/k$ and $E(-k^{2})\simeq k$ for large $k>0$, mentioned above. As shown in Fig. 3, $r$ obviously goes to [math] as well as $a$. Therefore we only consider $r\rightarrow 0$ limit.

First, using $a(0)=0$ and $E(1-1/r)\simeq 1/\sqrt{r}$ for small $r$, we can easily confirm that

[TABLE]

Here we used the fact that $K(1-1/r)$ goes to zero when taking the limit $r\rightarrow 0$.

Second, we confirm that $\Psi^{\prime}(0)=-2$. The function $\Psi(a)$ has already differentiated in eq. (5.20). Here the following two terms in eq. (5.20) vanish in the limit $r\rightarrow 0$,

[TABLE]

Next, the derivative of $a$ with regards to $r$ for $0\leq r<1$ in eq. (5.21) gives $a^{\prime}(0)=1/4$ and thus the first term in eq. (5.20) is equal to $-1$. Considering the limit of the contribution from $-\sqrt{r}E(1-1/r)$ after using the asymptotics $E(-1/r)\simeq 1/\sqrt{r}$, we can see that $-\sqrt{r}E(1-1/r)$ also gives $-1$. Combining the above calculation, we obtain

[TABLE]

6 Gaussian fluctuation with variance of $\mathcal{O}\left(\log t\right)$

As already mentioned in the introduction, in [14], the variance was shown to behave as $\braket{N(t)^{2}}_{c}\simeq(\log{t}+C)/2\pi^{2}$ and the constant $C$ was written as a sum of a few mathematical constants and definite integrals. This behavior is also observed in the context of the hard edge statistics in RMT and it is proven that the normalized random variable $(N(t)-\braket{N(t)})/\sqrt{\braket{N(t)^{2}}_{c}}$ converges to the Gaussian distribution [58]. The proof of the convergence is based on the discussion in [59]. In this section, we show one can recover this behavior by expanding the large deviation function obtained in the previous section around the most probable point $a=1/\pi$. As a byproduct, we find a much simpler formula for the constant $C$, see eq. (6.12). To see the main contribution, we can follow the same procedure as [55, 56], and as in their case, we can see the property that the large deviation function at the minimum is not analytical inherit. Therefore, we need to see the increment which emerges when we slide the value $a$ towards the horizontal axis only by small difference $\delta a>0$. As the eq. (5.19) gives the relation between $\delta a>0$ and $\delta r>0$, we can access the increment $\delta F$.

In the case of $0<a<1/\pi$ or $0<r<1$, we have

[TABLE]

Decreasing the parameter by $\delta r>0$, expansion of the function $F(r)$ in terms of $\delta r$ is as follows

[TABLE]

The transformation $a\rightarrow a-\delta a$ can be read from the eq. (5.19) as

[TABLE]

In the case of $1/\pi<a$ or $1<r$, we have

[TABLE]

The same calculation as above gives

[TABLE]

The transformation $a\rightarrow a+\delta a$ can be read from the eq. (5.19)

[TABLE]

The relation between $\delta a$ and $\delta r$ in eq. (6.3) is, in fact, the same as that in eq. (6.6). That is why it is enough to consider one of them. Inverting this relation, we have

[TABLE]

where $W(x)$ is the Lambert W-function [26] (see chapter 4.13), which is defined to be the solution of the following equation

[TABLE]

For $x<0$, notice that $W(x)$ has a branch because $W(x)$ is a multi-valued function, and we should take the one which takes the range of the function in $(\infty,-1/e]$ in the calculations below. For $W(x)\leq-1/e$, this function has an expansion

[TABLE]

as $x\rightarrow-0$ where $\eta=\ln\left(-1/x\right)$.

By using this asymptotic expansion of the Lambert W-function $W(x)$ and substituting eqs. (6.3), (6.7) and eq. (6.6), (6.7) into eq. (6.2) and eqs. (6.5) respectively, we obtain

[TABLE]

Therefore, after we restore the scale of $\delta n=t\delta a$,

[TABLE]

Inserting this into eq. (5.17), we could derive the main contribution to the variance, which increases logarithmically in time with the coefficient $1/2\pi^{2}$, as shown in [14].

To determine the number variance up to the correction $\mathcal{O}(1)$ in eq. (6.11), we need to carry out the direct calculation by using the Bessel kernel. In Appendix C, we show

[TABLE]

see eq. (C.24). Therefore we have found that the constant $C$ in eq. (1.1) is, in fact, simply given by

[TABLE]

which numerically matches precisely with the value from the calculation in [14].

In principle, we can address the cumulants higher than the second one by using the formula shown in eq. (C.1) as we calculate the second cumulant in Appendix C though the full direct calculation of all the cumulants is cumbersome and remains as a future problem. We can see numerically that a few subsequent cumulants are suppressed compared to the second cumulant and are around zero. This observation also supports the fact that higher-order cumulants are smaller than the second cumulant which diverges logarithmically in time [58] and, as a consequence, that the fluctuation of the spin current is the Gaussian when focusing on the vicinity of the mean on a scale $\mathcal{O}(\sqrt{\log t})$.

7 Numerical Analysis

We numerically verified the above formula of the large deviation function $\Psi(a)$ by DMRG calculations, using the ITensor C++ library [60]. The result is shown in Fig. 4 exhibiting a striking agreement between analytical formula and the DMRG calculation.

We employed the matrix product states (MPS) time evolution method with matrix product operator (MPO) form of time-evolution operator $e^{-it{H}}$[61] not only to simulate the real dynamics of the spin chain but also to obtain the cumulant generating function. For general information about MPS and MPO in DMRG, see for instance the reviews [62, 63].

The large deviation function $\Psi(a)$ can be calculated from the cumulant generating function of the integrated spin current [64]. However, as we have seen, the tail of the distribution has an exponent denoted by time squared $t^{2}$. This decay rate prevents us from applying the usual formula in large deviation theory. For adapting our situation to the form available for the large deviation formalism, let us modify the moment generating function $\chi(\lambda,t)$. Noting that $\lambda$ is just a parameter, we change the parameter $\lambda$ to $\lambda t$ and write $\chi(\lambda t,t)$ as $\tilde{\chi}(\lambda,t)$ after evaluating $\braket{e^{\lambda N(t)}}$. In this way, we can obtain the large deviation function from the moment generating function :

[TABLE]

via the Legendre-Fenchel transform:

[TABLE]

Thus, our goal is to obtain the numerical data of the expectation value $\left\langle e^{\lambda N(t)}\right\rangle$ at each $t$ and $\lambda$. This quantity can be recast into the Schrödinger picture:

[TABLE]

The procedure to evaluate this value is the following:

Calculate the state $\ket{\phi(t)}$ as an MPS by cumulatively applying the time evolution MPO $e^{-i\Delta t{H}}$, from an initial state $\ket{\phi(t=0)}$ of interest. In this study, we chose the Domain Wall state in eq. (2.3) with finite length $2L$. 2. 2.

Prepare the MPO form of the operator $e^{\pm\Delta\lambda N(0)}$ as in the case of the time evolution operator $e^{-i\Delta t{H}}$, and cumulatively apply it to the state $\ket{\phi(t)}$ up to $\pm\lambda/2$, when the desired value is $\pm\lambda$. It is required to generate the states $e^{\pm\lambda{N}(0)/2}\ket{\phi(t)}$ separately from $\ket{\phi(t)}$ in order to obtain the value of $\chi(\lambda)$ within the interval $[-\lambda,\lambda]$. 3. 3.

Evaluate the norm of the resulting state $e^{\pm\lambda{N}(0)/2}\ket{\phi(t)}$.

We can access to only finite values of $t$ and $\lambda$ in this manner because it is based on successive application of time evolution operator with finite time step. The range of $\lambda$ determine the upper and lower bound of the gradient of the large deviation function. Additionally, accurate calculation of cumulant generating function for large $t$ requires much more bond dimension—the number of the singular vectors preserved in MPS, whose corresponding singular values are greater than or equal to the cutoff $\epsilon$—of the MPS $\ket{\phi(t)}$ and much smaller time step to reach it than other calculations, such as spin expectation values and correlation functions. We evaluated the reliability of the result by comparing the return probability $\left|\braket{\phi(0)}{\phi(t)}\right|^{2}$ with the analytical expression $e^{-t^{2}/4}$ in eq. (3.14), which corresponds to the value of the rate function at the origin: ${\Psi}(0)$.

We set our calculation parameters as follows; the length of the system: $2L=13+13$, the final time and the time step: $t=8,\Delta t=2\times 10^{-6}$, the upper limit of the parameter of the cumulant generating function and its step: $\lambda=36,\Delta\lambda=5\times 10^{-3}$, the cutoff of singular values for $t$ and $\lambda$: $\epsilon_{t}=1\times 10^{-28},\epsilon_{\lambda}=1\times 10^{-13}$.

For these parameter values, the range of $a$ which could be numerically evaluated was restricted to $0\leq a<0.75$. In the region, we find a good agreement between the theory and numerics in Fig. 4. In the same figure, we also observe some deviation as $a$ increases. We expect that the agreement becomes better as $t$ becomes large. For instance, in another numerical calculation for a larger time $t=10$, we observed a better agreement than $t=8$, but the range of available $a$ was narrower for the reason that we need more information of $\lambda$ as $t$ becomes large.

8 Conclusion

In this work, we studied the fundamental XX model and investigated the large deviation property of the spin current starting from the domain wall initial condition for absolute zero temperature and perfect transmission. For this setup, we found that the values of the rare event for spin current in the right and left tail have an exponential decay in the form $e^{-t^{2}\Psi(a)}$. The time squared decay seems inconsistent with the expectation from the approximated Levitov-Lesovik formula which predicts the decay rate $t$, but this is expected to be a robust behavior in our settings. In addition to the exponent, the analytical large deviation function is exactly obtained. The large deviation function $\Psi(a)$ has a non-symmetric shape, which would be a reflection of many-body effects in a non-equilibrium situation.

From the point of view of the experiment, one can construct the XX spin chain [65], using bosonic atoms in an optical lattice and adjusting the parameter of the system through changing of the laser intensity or the Feshbach resonance. By trapping the hard core bosons undergoing Bose-Einstein condensation only in a center of periodic lattice potential plus harmonic potential by optical laser beams [1], the initial condition corresponding to our domain wall initial state could be prepared. Indeed, characteristic density profile of the form $\braket{s_{m}(t)}+1/2=\pi^{-1}\arccos\xi$ with $\xi=m/t$ has been observed by L. Vidmar et al. as shown Fig.2 in [66]. As discussed below the eq. (2.12), our results apply to both ferromagnetic and antiferromagnetic coupling. The antiferromagnet XX spin chain is realized by compounds such as $\mathrm{PrCl_{3}}$ [67] or $\mathrm{CsCoCl_{3}}$ [68]. One can apply strong magnetic field to these materials for preparing the domain wall. It would be interesting to measure the fluctuation of $N(t)$ in experiments to compare with our theoretical predictions.

As a future prospect, analyzing the details of the dynamics of XX model with other initial conditions would be a most natural direction. The analytical behavior of the large deviation for other models, in particular the XXZ spin chain with the same or other initial condition, are most interesting and challenging problems. Furthermore, we would like to clarify the connection between this work and the description of GHD extended to higher fluctuation [13]. It may also be interesting to see if our results and methods in this paper reveal some properties of the hard edge statistics, on the contrary.

9 Acknowledgments

The authors are grateful to T. Yoshimura for helpful discussions and comments, and to P. Forrester, G. M. Schütz, L. Vidmar for pointing out a few related reference. T. S. is also grateful to M. Katori, T. Imamura, T. Nagao and T. Shirai for their comments. H. M. also acknowledges the financial support from Advanced Research Center for Quantum Physics and Nanoscience, Tokyo Institute of Technology. The work of T. S. is supported by JSPS KAKENHI Grant Numbers JP15K05203, JP16H06338, JP18H01141, JP18H03672, JP19K03665.

Note added.—At the stage of the proofreading, we have realized that the large deviation function of the particle number in the hard edge is studied in a different fashion in [69]. The expression of large deviation in [69] seems different from ours but numerically looks almost the same.

Appendix A Determinant Identity

In this appendix, we illustrate how we get the relation in eq. (3.8). The determinant expression of the moment generating function $\chi(\lambda,t)$ in eq. (3.1) contains a single particle density operator restricted in a specific area at time $t$ as a correlation kernel

[TABLE]

Where $P$ is a projection operator onto the area that measurement would be performed, $U(t,0)$ is an unitary operator from 0 to $t$ generated by the single particle Hamiltonian in eq. (2.20) and $\rho(t)$ is a density operator of the system at time $t$. Expanding the cumulant generating function $\log{\chi(\lambda,t)}$ as

[TABLE]

enables us to handle the correlation kernel as a single particle operator. When we focus on the trace part in eq. (A.2), that has the form

[TABLE]

In the case we are considering, the initial condition is just given by projection operator onto the left chain

[TABLE]

Another projection operator $P$ for measurement of up spins which takes projection onto the site $\{\alpha+1,\alpha+2\dots\}$ are to be written as $P_{\alpha}$.

We now proceed to calculate the trace

[TABLE]

The unitary operator $U(t,0)$ in the coordinate representation is in eq. (2.21). To bring the propagator to that close to the form appearing in the Coordinate Bethe Ansatz [70, 71] (c.f. Theorem 2.1 in [71] for $N=1$), let us move on to $z$ plane to carry out the sum over all the momentum $k$ of plane waves $e^{ikn}$. Then,

[TABLE]

In the second equality, the integral representation of the Bessel function of the first kind [26]

[TABLE]

is used. The last equality is derived from the equivalence of the hopping rate to the nearest neighbor. Instead of performing the Wick rotation, we replace $iz$ by $z$, which does not change the contour of $z$

[TABLE]

The general property of evolution operator,

[TABLE]

and changing the variable $z\rightarrow 1/z$ lead to

[TABLE]

Substituting eq. (A.8) and eq. (A.10) into the trace expression in eq. (A.5), we have

[TABLE]

Note that $z_{2m+1}=z_{1}$, $n_{2k+1}=n_{1}$ and the like are assumed in this notation because of the periodicity originated from the trace. By taking the summation over their indices $n_{k}~{}(k=0,\dots,2m)$, one can carry out the trace as

[TABLE]

We use the following identity that Derrida and Gerschenfeld employed in [72],

[TABLE]

Going back to the trace and switching the corresponding part so that this form appears in the integrand, we obtain

[TABLE]

Here we can see that

[TABLE]

Similarly, we get

[TABLE]

Since the integrands do not have any poles with respect to $z_{k}$, because of $\alpha>0$, these terms readily equal zero by the residue theorem. Therefore, the calculation becomes rather simple as the poles which behave singularly no longer exist

[TABLE]

Rearranging the order of the exponent,

[TABLE]

holds. Changing the variables $z_{k}$ as $\sqrt{\dfrac{t_{k-1}}{t_{k}}}z_{k}$ leads to

[TABLE]

Now let us recall that the Bessel function of the first kind of order $\alpha$ has a complex contour representation.

[TABLE]

Hence, we arrive at the following expression

[TABLE]

Moreover one can integrated out $m$ variables alternatively, by taking new variables and using a formula related to two product of the Bessel functions [26]

[TABLE]

From this formula, we can proceed further

[TABLE]

We also used recurrence relation of $J_{\alpha}(t)$ in the third equality. In this way, we obtain

[TABLE]

In other words,

[TABLE]

The proof has been completed.

Appendix B Derivation of the large deviation function

In this section we derive the formula for the rate function $\Psi(a)$ in eq. (5.18) and eq. (5.19). The starting point is the eq. (5.4). At first, we decompose the rate function $\Psi(a)$ into four parts and name the first three terms as $I_{1},I_{2},I_{3}$,

[TABLE]

Each term is given by

[TABLE]

The resolvent $G(x)$ in eq. (5.9) has a different sign depending on the value of the argument $x$. From the definition of the resolvent $G(x)$ in eq. (5.8), when $r<1$,

[TABLE]

Similarly, when $r>1$,

[TABLE]

follows. According to the formula in [55] or from the direct calculation following [57], we know that

[TABLE]

The expression above has nothing to do with the order of $r,1$ on $[0,\infty)$.

For the other two terms $I_{2},I_{3}$, we expand those into the series of $1/N$. In the case of $r<1$,

[TABLE]

where $K$ and $E$ are the elliptic integrals defined in eq. (5.16). In the first equality, we used, if $a<b$,

[TABLE]

and the definition of $r$ in eq. (5.14). In the case of $r>1$, following to the previous case, we obtain a similar expression,

[TABLE]

Next, we write the integral part in the third term $I_{3}$ in eq. (B.4).

[TABLE]

When $x<1$, the square root of the form $\sqrt{1-x}$ can be expanded in terms of $x$ as Taylor series,

[TABLE]

Substituting this formula in the above expression, the integrand in eq. (B.11) is written as

[TABLE]

Expanding this integrand and put together the contributions up to the order $\mathcal{O}(1/N^{2})$, this becomes

[TABLE]

Performing the integration from $\lambda_{+}$ to $\infty$, we obtain

[TABLE]

The two infinite series appearing in the first and the third term converge to the following values,

[TABLE]

[TABLE]

Now, we return back to the original form of $I_{3}$ and we can see that

[TABLE]

In the last equality, we used the definition of $\lambda_{-}$ in eq. (5.14). Combining all the terms in eq. (B.7), (B.10), and (B.18) evaluated separately into eq. (B.1), multiplying $\psi_{I}$ by $2N^{2}$ and taking the large $N$ limit, we finally arrive at the expression for the large deviation function $\Psi(a)$ in eq. (5.18). Here we can see that the term which has the order $\mathcal{O}(1)$ in $I_{1},~{}I_{3}$ and the remaining constant $3/4$ in eq. (B.1) are canceled.

In a similar way, evaluating the constraint in eq. (5.12) so that we do not leave the fictitious parameter $N$ in the large $N$ limit, we can obtain eq. (5.19) as a form which does not include the parameter $N$. Since the right hand side of eq. (5.12) has only $1/N$ dependence, we can ignore its higher correction smaller than $1/N$.

For $0\leq r<1$, equivalently for $0\leq a<1/\pi$, eq. (5.12) becomes

[TABLE]

Recalling that $\lambda_{+}$ has an expansion as $\lambda_{+}=4+\mathcal{O}(1/N^{2})$ from eq. (5.10), and taking the large $N$ limit, we have,

[TABLE]

For $1<r$, equivalently for $1/\pi<a$, we can apply the same procedure and get the corresponding result as

[TABLE]

Appendix C Constant in the Variance

In section 6, we investigated the main contribution to the variance of the integrated spin current $\braket{N(t)^{2}}_{c}$ for large $t$. In this appendix, we should determine the further correction up to $\mathcal{O}(1)$ order, by using the asymptotics of the Bessel kernel.

Making use of the expansion of the cumulant generating function $\log\chi(\lambda,t)$, we see

[TABLE]

where $W(m,n)=(-1)^{n}n!S(m+1,n+1)$ is the (signed) Worpitzky number and $S(m,n)$ is the Stirling number of the second kind [26] (see chapter 26.8). One of the exponential generating function of $S(m,n)$ [26] is given as follows,

[TABLE]

From the second coefficient, we have

[TABLE]

This kind of trace is explicitly written down as follows by using eq. (A.21),

[TABLE]

Moreover, using the representation of the Dirac delta function in terms of the Bessel function $J_{\nu}(t)$ [26],

[TABLE]

we can unite the two multiple integral terms as

[TABLE]

if one integrates the variables $t_{1}$ and $t_{3}$, this leads to the double integral,

[TABLE]

Applying asymptotic expansion of the Bessel function $J_{\nu}(t)$ and its derivative $J^{\prime}_{\nu}(t)$ for a large argument [26],

[TABLE]

we can perform the integration in eq. (C.7) with elementary functions. If one is interested in higher corrections, one could take subleading corrections of $J_{\nu}(t)$ and $J^{\prime}_{\nu}(t)$ at this point.

The large $t$ behavior of this integral in eq. (C.7) is reduced to the following integral,

[TABLE]

Expanding the square in the integrand, we decompose the integral into three parts. Let us name them as $J_{1},J_{2},J_{3}$ for simplicity,

[TABLE]

where

[TABLE]

respectively. We first integrate them with respect to $y$. In the following, we are to use sine integral and cosine integral defined by [26]

[TABLE]

[TABLE]

For $J_{1}$, we have

[TABLE]

Before the integration of $J_{2}$, let us rewrite the integrand to perform the integration easily as

[TABLE]

This leads to

[TABLE]

The calculation of $J_{3}$ is similar to that of $J_{1}$. Changing the variable as $y\rightarrow-y$ in the process of its calculation, we obtain

[TABLE]

Since the second term $J_{2}$ vanishes for large $t$, we integrate the remaining two terms $J_{1}$, $J_{3}$,

[TABLE]

Simplifying each term, we get

[TABLE]

where

[TABLE]

is the Euler-Mascheroni constant. As the time $t$ becomes large, the following terms go to zero as

[TABLE]

An asymptotic property such that $\lim_{x\to\infty}\mathrm{Ci}(x)=0$ implies that both $\mathrm{Ci}(4t)$ and $2\mathrm{Ci}(2t)$ also vanish.

In conclusion, we obtain the variance for large $t$ as

[TABLE]