\section{Introduction}

In many problems in quantum physics one is led to study Schr\"odinger operators of the type
$$H_V=-\Delta+V,$$ where $V$ describes a complex system. For example, in the theory of disordered quantum systems it is common to apply a one body approximation to reduce the many body Schr\"odinger operator to a Schr\"odinger operator of the form
\beq\label{complex_Schr}
H_\Omega=-\Delta+\sum_{\omega\in\Omega}V(\cdot-\omega),
\eeq
where $\Omega$ denotes the set of positions of the obstacles, and, to simplify the analysis, one usually assumes the potential $V$ to be continuous and compactly supported (although more general potentials may be considered).

In 1958, Anderson demonstrated, in the case of a lattice model, that in a regime where the disorder is sufficiently strong the quantum eigenfunctions become exponentially localized -- a phenomenon known as Anderson localization \cite{A}. Anderson localization has been proved rigorously in various models, discrete as well as continuous \cite{GMP,FS,GHK,BK}. An important question concerns the existence of a transition from a localized to a delocalized regime (Anderson transition). The scaling theory of Anderson, Abrahams, Licciardello and Ramakrishnan \cite{AALR} predicts that such a delocalization transition should occur in dimension $d\geq 3$, whereas  the presence of disorder should suffice to localize all eigenfunctions in dimensions $d=1,2$. The case $d=2$ is critical in the sense that the localization length tends to be very large so that in practice (physical experiments and numerical simulations) a delocalization transition is often observed, because the localization length exceeds the physical dimensions of the system.

Whereas Anderson localization is by now well-understood at the mathematical level \cite{BK,GMP,FS,GHK}, a proof of Anderson delocalization remains elusive. Progress has, however, been made by Aizenman-Warzel \cite{AW} who constructed extended states on trees, and by Erd\"os-Salmhofer-Yau \cite{ESY} who prove quantum diffusion of the time-dependent random Schr\"odinger equation in a scaling limit. Moreover, Erd\"os-Yau and Erd\"os-Eng \cite{EY,EE} derived that the quantum dynamics is governed by a linear Boltzmann equation in the weak coupling and low density limits.

In the present article, we prove a quantitative equidistribution theorem for the eigenfunctions of a Schr\"odinger operator with an $L^2$ potential on a rectangular $2$-torus. Although we state our main results for the case of real-valued potentials on rectangular tori, which is the most interesting case for our purposes, our proof should work for any rectangular domain with general boundary conditions. In particular, it may be applied to Schr\"odinger operators with quasi-periodic boundary conditions which arise via Bloch-Floquet theory and play an important role in the study of quantum transport in graphene, as discussed by Fefferman-Weinstein in \cite{FW}. We discuss such applications in section \ref{sec-BC}.

We apply our equidistribution theorem to the case of Schr\"odinger operators with potentials describing a disordered system, such as \eqref{complex_Schr}. We study the thermodynamic limit, as the number of potentials $N$ on the fixed torus $\T^2$ tends to infinity, where the support of the potentials has diameter $O(N^{-1/2})$. Under a weak disorder hypothesis on the set of positions of the potentials, we derive an equidistribution theorem which is scale-invariant in the thermodynamic limit, as $N\to\infty$. This application of our theorem is discussed in section \ref{sec-complex}. The proof of the main theorem is given in section \ref{sec-equi}.

We derive a corollary which deals with the case of random potentials, where we prove equidistribution of eigenfunctions which holds almost surely with respect to the joint distribution of the random positions and is scale-invariant in the thermodynamic limit. Another corollary deals with the general (possibly strong disorder) regime, where we derive a scale-invariant equidistribution theorem in the thermodynamic limit, provided that energy and coupling strength of the potential satisfy certain constraint. This corollary gives rise to a lower bound for the localization length of a general random Schr\"odinger operator as a function of energy and coupling parameter.

\subsection*{Acknowledgements}
The author's work was supported by the grant ANR-17-CE40-0011-01 (project SpInQS) of the French Nantional Research Agency (ANR) and a d\'el\'egation CNRS. The author is grateful to Jens Marklof and Gabriel Rivi\`ere for many helpful comments.

\subsection{Notation}

Let $\T^2=\R^2/2\pi\cL_0$, where $\cL_0=\Z(a,0)\oplus\Z(0,1/a)$, $a>0$, is a rectangular unimodular lattice. We denote by $\Delta=\partial_x^2+\partial_y^2$ the Laplacian on $\T^2$. So the Laplace eigenvalues are values of the quadratic form $Q(\xi)=|\xi|^2=a^2x^2+a^{-2}y^2$, $(x,y)\in\Z^2$.

We denote the set of ordered Laplace eigenvalues, counted without multiplicity, as follows
\beq
\cQ=\{0=n_0<n_1<n_2<\cdots<\cdots<n_k<\cdots\to+\infty\}.
\eeq

If $V\in L^2(\T^2,\R)$, then the Schr\"odinger operator $-\Delta+V$ on the domain $H^2(\T^2)$ is essentially self-adjoint and has discrete spectrum (cf. section 2.2 in \cite{BBZ} for a detailed discussion). Throughout this paper we will assume that $V$ is real-valued which is the most interesting case for our purposes. However, our theorem can easily be adapted for complex-valued potentials and/or general boundary conditions. We briefly discuss this in subsection \ref{sec-BC}.

We will consider the eigenvalue problem
\beq
(-\Delta+V)\psi_\lambda=\lambda\psi_\lambda, \quad \|\psi_\lambda\|_{L^2(\T^2)}=1.
\eeq

We have the following $L^2$ representation for the eigenfunctions
\beq
\psi_\lambda=\sum_{\xi\in\cL}\hat{\psi}(\xi)e_\xi, \quad e_\xi(x):=\frac{1}{2\pi}e^{i\left\langle\xi,x\right\rangle}
\eeq
where $\cL$ denotes the dual lattice of $\cL_0$.

Moreover, for $\lambda\in (n_k,n_{k+1})$ let us define the following functions which are truncated outside an annulus
\beq
\psi_{\lambda,L}(x)=\sum_{\xi\in A(n_k,L)}\hat\psi_\lambda(\xi)e_\xi(x)
\eeq
where $L>0$ and $$A(n_k,L):=\{\xi \in \cL \mid ||\xi|^2-n_k|\leq L\}.$$

\subsection{Main theorem}
For a subset of the ordered set of Laplace eigenvalues, $\cQ'\subset \cQ$, we denote the following subset of the real line $$\Sigma:=\bigcup_{n_k\in \cQ'}(n_k,n_{k+1}).$$
We prove the following quantitative equidistribution theorem for the eigenfunctions of the Schr\"odinger operator $H_V=-\Delta+V$ on $\T^2$, where we recall $V\in L^2(\T^2)$.
\begin{thm}\label{equi_thm}
There exists a subset of full density $\cQ'\subset \cQ$ such that for any orthonormal basis of eigenfunctions of $H_V$, denoted by $\{\psi_k\}_{k\geq0}$, $H_V\psi_k=\lambda_k\psi_k$, and for any test function $a\in C^\infty(\T^2)$, we have for any $k$ such that $\lambda_k\in\Sigma$ the bound
\footnote{For two functions $f,g\geq 0$ we denote by $f\lesssim g$ that there exists a positive constant $C$ s. t. $f\leq Cg$.}
\beq\label{equi_bound}
\Big|\int_{\T^2} a |\psi_k|^2 d\mu-\frac{1}{4\pi^2}\int_{\T^2}a d\mu\Big| \lesssim_{a,\epsilon}
 \|V\|_{L^2(\T^2)}\lambda_k^{-(1-3\theta)/2+\epsilon}
\eeq
where the implied constant depends on a certain number of derivatives of $a$, $d\mu$ denotes the Lebesgue measure on $\T^2$ and $\theta$ denotes the exponent in the error term of the circle law. \footnote{The Gauss circle law asserts $\#\{\xi\in\cL \mid |\xi|^2\leq T\} = \pi T + O_\epsilon(T^{\theta+\epsilon})$.}
\end{thm}
\begin{remark}
Under the assumption of the Gauss circle problem, which asserts that $\theta=1/4$, we get an optimal equidistribution rate of $O_\epsilon(\lambda_k^{-1/8+\epsilon})$. The best known exponent to date is due to Bourgain-Watt \cite{BW}, $\theta=517/1648$, so our theorem gives a rate of $O_\epsilon(\lambda_k^{-97/3296+\epsilon})$.
\end{remark}

\subsection{General boundary conditions, non-self-adjoint operators}\label{sec-BC}
We state theorem \ref{equi_thm} in the, for our purposes, most natural case of a self-adjoint operator on the torus $\T^2$. However, our proof may easily be adapted to deal with any rectangular domain and boundary condition as well as complex-valued potentials. In this case the same equidistribution theorem holds for eigenfunctions which are solutions to the eigenvalue problem
\beq
(-\Delta+V)\psi_\lambda=\lambda\psi_\lambda
\eeq
on a rectangular domain $\scrR=[0,2\pi a]\times[0,2\pi/a]$, provided the eigenvalue satisfies the condition $\re\lambda\in\Sigma$. The boundary conditions are in fact of no importance in the proof.

In order to adapt the proof one simply replaces the bound \eqref{Fourier_bound} with
$$|\hat\psi_\lambda(\xi)|^2\leq \|V\|_{L^2(\scrR)}^2(|\xi|^2-\re\lambda)^{-2}$$
and uses the assumption $\re\lambda\in\Sigma$ throughout the proof.

We thus have for any test function $a\in C^\infty(\T^2)$ and for any $k$ such that $\re\lambda_k\in\Sigma$ the bound
\beq
\Big|\int_{\T^2} a |\psi_k|^2 d\mu-\frac{1}{4\pi^2}\int_{\T^2}a d\mu\Big| \lesssim_{a,\epsilon} \|V\|_{L^2(\T^2)}(\re\lambda_k)^{-(1-3\theta)/2+\epsilon}.
\eeq

\subsection{Eigenfunctions in disordered systems}
Let us now consider the case of a potential which models a disordered system of $N$ scatterers of diameter $O(N^{-1/2})$ on $\T^2$
\beq\label{thermo_pot}
V_N(x)=\sum_{\omega_j\in\Omega_N} V(N^{1/2}(x-\omega_j))
\eeq
where $V\in L^2(\T^2)$, $\supp V\subset B(0,r)\subset\T^2$ and
$$\Omega_N:=\{\omega_1,\cdots,\omega_N\}\subset\T^2.$$

For fixed $N$ one may apply Theorem \ref{equi_thm} to the Schr\"odinger operator $H_N=-\Delta+V_N$ and one will obtain an error term which depends on $N$ (cf. section \ref{sec-complex}). So the equidistribution theorem is in general not valid in the thermodynamic limit, as $N\to+\infty$.

The key question is whether, under a {\em weak disorder hypothesis} on the configuration of the set $\Omega_N$, one may strengthen the theorem in such a way that the error term obtained is {\em uniform with respect to the parameter $N$}. This means that, if the disorder in the system is sufficiently weak and the energy is above a certain threshold, then {\em delocalized eigenfunctions do exist} and their delocalization persists in the limit as $N\to\infty$, which means the delocalization is scale-invariant (we can think of the limit $N\to\infty$ as ``zooming'' out to larger and larger scales).

This question is answered in the affirmative by the following equidistribution theorem for the eigenfunctions of the Schr\"odinger operator $H_N=-\Delta+V_N$ which is valid in the thermodynamic limit if a certain weak disorder hypothesis is satisfied. The proof is given in section
\ref{subsec-thermo}.
\begin{thm}\label{thermo}
Let us suppose that the set $\Omega_N$ satisfies
\beq\label{rigid}
\#\Omega_N\cap B(x_0,\frac{R}{N^{1/2}}) \lesssim \max(1,R^2).
\eeq
for any $x_0\in\T^2$ and $0<R\lesssim N^{1/2}$.
Let $H_N=-\Delta+V_N$, where $V_N$ is given by \eqref{thermo_pot}. Then for any orthonormal basis of eigenfunctions of $H_V$ denoted by $\{\psi_{N,k}\}_{k\geq0}$, $H_N\psi_k=\lambda_{N,k}\psi_{N,k}$, and for any test function $a\in C^\infty(\T^2)$, we have for any $k$ such that $\lambda_{N,k}\in\Sigma$ the bound
\beq\label{equi_bound}
\Big|\int_{\T^2} a |\psi_{N,k}|^2 d\mu-\frac{1}{4\pi^2}\int_{\T^2}a d\mu\Big| \lesssim_{a,\epsilon} r^2\|V\|_{L^2(\T^2)}(\lambda_{N,k})^{-(1-3\theta)/2+\epsilon}
\eeq
where we recall that $\supp V\subset B(0,r)$.
\end{thm}

\begin{remark}
Note that condition \eqref{rigid} is satisifed by distorted lattices. If we consider, for simplicity, $\Z^2\subset \R^2$ and we introduce finite displacements at each lattice point, we obtain the distorted lattice
$$\Omega=\{\xi+\omega_\xi \mid \xi\in\Z^2, \; |\omega_\xi|\leq r_0\}.$$
So the set $\Omega_N:=N^{-1/2}\Omega$ satisfies our condition for any ball $B(x_0,R/N^{1/2})$, because
$$|N^{-1/2}(\xi+\omega_\xi)-x_0|\leq R/N^{1/2}$$ implies
$$|\xi-N^{1/2}x_0|\leq |(\xi+\omega_\xi)-N^{1/2}x_0|+|\omega_\xi|\leq R+r_0$$
and the number of lattice points in this ball is $\sim\pi (R+r_0)^2$.
\end{remark}

\subsection{Random potentials}
We may apply Theorem \ref{thermo} to a popular model of disordered quantum systems, so-called random displacement models (RDM). Whereas Anderson localization was proved for these models at the bottom of the spectrum \cite{KLNS}, we note that the RDM for $N$ obstacles on $\T^2$ which we study here, when rescaled to a torus of size $N^{1/2}$, corresponds to a weak coupling strength $1/N$ (as opposed to the model considered in \cite{KLNS}), as well as a low energy $E/N$. We discuss the case of general random Schr\"odinger operators below.

We will show that delocalization occurs for suffiently high energy in our RDM, in the sense that we prove an equidistribution theorem for the eigenfunctions which is scale-invariant in the thermodynamic limit. We point out that random Schr\"odinger operators such as those arising in RDM will satisfy condition \eqref{rigid} almost surely if the random variables satisfy a weak disorder assumption.
Let $\PP$ be a compactly supported probability density on $\R^2$ where $\supp\PP\subset B(0,r_1)$ for some $r_1>0$. We consider a sequence $\{\xi_j\}_{j=1}^{+\infty}$ of i.i.d. random variables with probability densities $\PP_j=\PP$. So the sequence $\{\xi_j\}_{j=1}^{+\infty}$ has a joint probability density $\PP_\infty=\bigotimes_{j=1}^\infty\PP_j$. The proof of the following corollary of Theorem \ref{thermo} will be given in section \ref{subsec-random}.
\begin{coro}
Let $\bar{\Omega}_N=\{\bar\omega_1,\cdots,\bar\omega_N\}\subset\T^2$ be a fixed point set which satisfies \eqref{rigid}. $\Omega_N=\{\omega_1,\cdots,\omega_N\}$ be a random set, where $\omega_j=\bar\omega_j+\xi_j$. Suppose $V\in L^2(\T^2)$ with $\supp V\subset B(0,r_2)$.

We consider the random Schr\"odinger operator
\beq
H_N=-\Delta+\sum_{\omega_j\in\Omega_N}V(N^{1/2}(x-\omega_j))
\eeq
Under the assumptions above we have $\PP_\infty$-a.s. that, for any $E>0$, any orthornormal basis of eigenfunctions of the random Schr\"odinger operator $H_N$ and for any test function $a\in C^\infty(\T^2)$, for all eigenvalues $\lambda_k^N$ of $H_N$ such that $\lambda_k^N\in\Sigma\cap[E,2E]$, the associated eigenfunctions satisfy
\beq\label{equi_bound}
\Big|\int_{\T^2} a |\psi_k^N|^2 d\mu-\frac{1}{4\pi^2}\int_{\T^2}a d\mu\Big| \lesssim_{a,\epsilon} (r_1+r_2)^2\|V\|_{L^2(\T^2)}E^{-(1-3\theta)/2+\epsilon}.
\eeq
\end{coro}

Let us now consider the case of a Schr\"odinger operator on $\T^2_L=\R^2/(2\pi L\Z)^2$, where $L>0$ is a large parameter,
\beq
H_L=-\Delta+\alpha\sum_{j=1}^N V(x-x_j), \quad H_L\Psi=E\Psi,
\eeq
$E$ denotes the energy, $\alpha=\alpha(L)$ a coupling parameter, and $\Omega_N=\{x_1,\cdots,x_N\}\subset\T^2_L$.
If we scale back to the torus $\T^2=\R^2/(2\pi\Z)^2$, then we obtain
\beq
h_L=-\Delta+\alpha L^2 \sum_{j=1}^N V(Lx-x_j), \quad h_L\psi=\lambda\psi, \quad \lambda=EL^2.
\eeq

Let $\Omega\subset\R^2$ be a point set which is of positive density in the sense that
$$L^{-2}\#\Omega\cap[-L,L]^2 \asymp \rho$$ for some $\rho>0$. If we now apply periodic boundary conditions and take $\Omega_N=\Omega\cap 2\pi[-L,L]^2$, then this gives rise to a Schr\"odinger operator $H_L$ on $\T^2_L$, and in turn $h_L$ on $\T^2$. We then have the following corollary of Theorem \ref{equi_thm}, which we prove in section \ref{sec-strong}. In particular, me may take $\Omega$ to be a random point set which is generated by a suitable stochastic process (e. g. a Poisson process or a RDM).
\begin{coro}\label{strong}
For any $E>0$, any orthornormal basis of eigenfunctions of the Schr\"odinger operator $h_L$ and for any test function $a\in C^\infty(\T^2)$, for all eigenvalues $\lambda_k^L$ of $h_L$ such that $\lambda_k^L\in\Sigma\cap [EL^2,2EL^2]$, the associated eigenfunctions satisfy
\beq\label{equi_bound}
\Big|\int_{\T^2} a |\psi_k^L|^2 d\mu-\frac{1}{4\pi^2}\int_{\T^2}a d\mu\Big|
\lesssim_{a,\epsilon} \rho|\alpha|\|V\|_{L^2(\T^2)}\frac{L^{2+3\theta+\epsilon}}{E^{(1-3\theta)/2-\epsilon}}.
\eeq
\end{coro}
\begin{remark}
We note that in a regime, where\footnote{
We denote by $\ll$ ``much smaller than'', and by $\gg$ ``much large than''.
}
\beq\label{equi_cond}
\rho|\alpha|\|V\|_{L^2(\T^2)}L^{2+3\theta+\epsilon}\ll E^{(1-3\theta)/2-\epsilon},
\eeq
we have equidistribution in the thermodynamic limit, as $L\to+\infty$.

In particular, this gives a lower bound for the length scale of Anderson localization as a function of the energy $E$, the density $\rho$ and the coupling parameter $\alpha$:
\beq
L_{loc}(\alpha,E,\rho,V)\gg \left(\frac{E^{(1-3\theta)/2-\epsilon}}{\rho|\alpha|\|V\|_{L^2(\T^2)}}\right)^{\frac{1}{2+3\theta+\epsilon}}
\eeq

We stress that one may, of course, improve this lower bound significantly by exploiting the particular properties of the stochastic process which generates the point set $\Omega$. For example, in the case of RDM, we used the weak disorder assumption \eqref{rigid} to improve the term $L^{2+3\theta+\epsilon}$ to $L^{1+3\theta+\epsilon}$. Moreover, we were considering a weak coupling/low energy regime, where $\alpha=L^{-2}$, $E=L^{-2}E_1$, where $E_1$ is the fixed energy on $\T^2$. So, condition \eqref{equi_cond} becomes $\rho\|V\|_{L^2(\T^2)}\ll E_1^{(1-3\theta)/2-\epsilon}$ which shows that scale-invariant equidistribution holds, when $E_1$ is large enough with respect to the density of the scatterers and the $L^2$ norm of the potential.
\end{remark}

\section{Application to disordered systems}\label{sec-complex}

\subsection{General case}\label{sec-strong}
Let us consider Schr\"odinger operators of the type
\beq
H_N=-\Delta+V_N
\eeq
where $V\in L^2(\T^2)$ and
\beq
V_N(x)=\sum_{\omega_j\in\Omega_N} V(x-\omega_j).
\eeq
where $$\Omega_N:=\{\omega_1,\cdots,\omega_N\}\subset\T^2.$$

Given $\Omega_N$ let $\{\psi_j^N\}_{j\geq0}$, $H_\omega\psi_j^N=\lambda_j\psi_j^N$, be an orthonormal basis of eigenfunctions of the operator $H_N$. Our equidistribution theorem yields, for any $a\in C^\infty(\T^2)$ with $\int_{\T^2} a d\mu=0$, and for any $j$ such that $\lambda_j^N\in \Sigma$,
\beq
|\bra a\psi_j^\omega,\psi_j^\omega\ket| \lesssim_a \|V_N\|_{L^2(\T^2)}(\lambda_j^N)^{-\beta}.
\eeq

Moreover, we have, for any $\Omega_N\subset\T^2$,
\beq\label{uni_bound}
\|V_N\|_{L^2(\T^2)}\leq N\|V\|_{L^2(\T^2)}.
\eeq

If we recall the setting of corollary \ref{strong}, we have $N\asymp \rho L^2$ and our potential $V_N$ is of the form
$$V_N=\alpha L^2\sum_{j=1}^N V(Lx-x_j)$$
so we get $$\|V_N\|_{L^2(\T^2)}\leq N|\alpha|L^2\|V(L\cdot)\|_{L^2(\T^2)}=N|\alpha|L\|V\|_{L^2(\T^2)}\ll \rho|\alpha|L^3\|V\|_{L^2(\T^2)}$$
which yields the error term in the corollary in view of $\lambda=EL^2$.

\subsection{Weak disorder}\label{subsec-thermo}
An important question concerns the validity of our equidistribution theorem in the limit as $N\to+\infty$. Let us consider potentials whose support is of diameter $O(N^{-1/2})$. We assume $\supp V\subset B(0,r)$, for some $r>0$, and consider the potential
\beq
V_N(x)=\sum_{\omega_j\in\Omega_N}V(N^{1/2}(x-\omega_j))
\eeq

Let us assume that the set $\Omega_N$ satisfies the following {\bf weak disorder hypothesis}. We suppose that for any $x_0\in \T^2$ and $0<R\leq N^{1/2}$ we have
\beq
\Omega_N\cap B(x_0,\frac{R}{N^{1/2}}) \lesssim \max(1,R^2).
\eeq

Under this assumption, we have the following bound on the $L^2$ norm of $V_N$.
First of all, we may rewrite
$$
\int_{\T^2}|V_N|^2 d\mu =\sum_{\omega_j\in\Omega_N} \sum_{|\omega_k-\omega_j|\leq 2 r/N^{1/2}} \int_{\T^2} V(N^{1/2}x)\overline{V(N^{1/2}(x-\omega_j+\omega_k))}dx
$$
where we note that $|\omega_j-\omega_k|>2r/N^{1/2}$ implies
$$|N^{1/2}(x-\omega_j+\omega_k)|\geq N^{1/2}|\omega_j-\omega_k|-N^{1/2}|x|>r$$
where we used $N^{1/2}x\in\supp V\subset B(0,r)$.

Now
\begin{align*}
&\Big|\int_{\T^2} V(N^{1/2} x)\overline{V(N^{1/2} (x-\omega_j+\omega_k))}dx\Big|\\
\leq& (\int_{\T^2} |V(N^{1/2} x)|^2dx)^{1/2}
(\int_{\T^2} |V(N^{1/2} (x-\omega_j+\omega_k))|^2dx)^{1/2}\\
=&\int_{\T^2} |V(N^{1/2} x)|^2dx=N^{-1}\|V\|_{L^2(\T^2)}^2
\end{align*}
because of translation invariance.

The weak disorder condition \eqref{rigid} now yields
\beq
\int_{\T^2}|V_N|^2 d\mu\lesssim \#\Omega_N N^{-1}\|V\|_{L^2(\T^2)}^2 r^2 = \|V\|_{L^2(\T^2)}^2 r^2
\eeq
so that the conditions of Theorem \ref{equi_thm} are satisfied, and, moreover, this yiels an equidistribution bound \eqref{equi_bound} which is uniform with respect to $N$.

\subsection{Random displacement models}\label{subsec-random}
If we take the positions $\omega_j$ to be independent random variables with probability densities whose compact supports are disjoint and of diameter $O(N^{-1/2})$, then we may ensure that the weak disorder hypothesis \eqref{rigid} is satisfied almost surely, and we may, thus, apply Theorem \ref{thermo} to obtain an equidistribution bound which is uniform with respect to $N$, and, therefore, valid in the thermodynamic limit, as $N\to+\infty$.

We first construct an infinite sequence of random variables on $\R^2$. Let $\PP$ be a compactly supported probability density on $\R^2$ where $\supp\PP\subset B(0,r_1)$ for some $r_1>0$. We consider a sequence $\{\xi_j\}_{j=1}^{+\infty}$ of i.i.d. random variables with probability densities $\PP_j=\PP$. So the sequence $\{\xi_j\}_{j=1}^{+\infty}$ has a joint probability density $\PP_\infty=\bigotimes_{j=1}^{+\infty}\PP_j$.

Let us suppose that $\bar\Omega_N\subset\T^2$ is a set of $N$ points on $\T^2$ which satisfy the weak disorder condition \eqref{rigid}. Let us introduce the random set $\Omega_N=\{\omega_j\}_{j=1}^N$, $\omega_j=\bar\omega_j+N^{-1/2}\xi_j$.

It follows that $\Omega_N$ satisifes the following bound almost surely for any $x_0\in\T^2$, $R>0$:
$$\#\Omega_N\cap B(x_0,R/N^{1/2})\lesssim \max(1,(r_1+R)^2).$$
This follows from the observation that if $|\omega_j-x_0|\leq R/N^{1/2}$, then, almost surely,
$$|\bar\omega_j-x_0|\leq |\bar\omega_j-\omega_j|+|\omega_j-y_0|\leq (r_1+R)N^{-1/2}$$
and, therefore, almost surely,
$$\#\{j\mid |\omega_j-x_0|\leq R/N^{1/2}\}\leq \#\{j\mid |\bar\omega_j-x_0|\leq r_1+R\}\lesssim \max(1,(r_1+R)^2).$$

So $\Omega_N$ satisfies our weak disorder hypthosesis almost surely.
Then, we have for an orthonormal basis of eigenfunctions $\{\psi_j^\omega\}_{j\geq0}$, for any $a\in C^\infty(\T^2)$ with $\int_{\T^2} a d\mu=0$, and for any $j$ such that $\lambda_j\in \Sigma$,
\beq
|\bra a\psi_j^\omega,\psi_j^\omega\ket| \lesssim_{a,\epsilon} (r_1+R)^2\|V\|_{L^2(\T^2)}\lambda_j^{-(1-3\theta)/2+\epsilon}.
\eeq

\section{Proof of Theorem \ref{equi_thm}}\label{sec-equi}
Let $a\in C^\infty(\T^2)$ with $\int_{\T^2}a d\mu=0$ and $L:=n_k^\delta$.

We have
\beq
\bra a\psi_\lambda,\psi_\lambda\ket=\bra a\psi_{\lambda,L},\psi_{\lambda,L}\ket+R(\psi_\lambda,L)
\eeq
where the remainder is of the form
$$R(\psi_\lambda,L):=\bra a\psi_{\lambda,L},\psi_{\lambda,L}^R\ket+\bra a\psi_{\lambda,L}^R,\psi_{\lambda,L}\ket+\bra a\psi_{\lambda,L}^R,\psi_{\lambda,L}^R\ket$$ where we denote $\psi_{\lambda,L}^R:=\psi_\lambda-\psi_{\lambda,L}$.

\subsection{Construction of $\cQ'$}

We will construct a full density subset $\cQ_1\subset \cQ$ such that for certain small $\delta>0$ and any $n\in \cQ_1$ the thin annuli $A(n,n^\delta)$ have the property that their lattice points are well-spaced. More precisely, we have the following proposition which we will prove in section \ref{sec-good}.
\begin{prop}\label{good_ann}
Let $\delta\in(\theta/2,1/2-\theta)$ and $\epsilon<\tfrac{1}{2}-\theta-\delta$. There exists a full density subset $\cQ_1\subset \cQ$ such that for any $n\in \cQ_1$ we have
that $\xi\in A(n,n^\delta)$ implies that for any $\zeta\in\cL\setminus\{0\}$, $|\zeta|\leq n^\epsilon$, we have $||\xi+\zeta|^2-n|\gg n^\delta$.
\end{prop}

We recall (cf. \cite{RU}, Lemma 2.1) that the subset
$$\cQ_2:=\{n_k \in S \mid |n_{k+1}-n_k|\leq n_k^{o(1)}\}$$
is of full density in $\cQ$, where $n_k^{o(1)}\to+\infty$, as $n_k\to+\infty$.

We define $\cQ':=\cQ_1\cap \cQ_2$.

\subsection{Bound on $R(\psi_\lambda,L)$}
Clearly,
$$R(\psi_\lambda,L)\leq 2 \|\psi_{\lambda,L}^R\|_{L^2(\T^2)}+\|\psi_{\lambda,L}^R\|_{L^2(\T^2)}^2,$$
so we have to estimate $\|\psi_{\lambda,L}^R\|_{L^2(\T^2)}$.

We have
$$\|\psi_{\lambda,L}^R\|_{L^2(\T^2)}^2=\sum_{\xi\notin A(n_k,L)}|\hat\psi_\lambda(\xi)|^2.$$
Recall that $\lambda\in(n_k,n_{k+1})$. From the identity $\psi_\lambda=(\Delta+\lambda)^{-1}(V\psi_\lambda)$ we deduce
\beq
\hat\psi_\lambda(\xi)=(-|\xi|^2+\lambda)^{-1}\int_{\T^2}e_\xi(-y)V(y)\psi_\lambda(y)dy
\eeq
and by Cauchy-Schwarz
\beq\label{Fourier_bound}
|\hat\psi_\lambda(\xi)|^2 \leq \|V\|_{L^2(\T^2)}^2(|\xi|^2-\lambda)^{-2}.
\eeq

In view of this estimate, we have
\beq
\|\psi_{\lambda,L}^R\|_{L^2(\T^2)}^2=\sum_{\xi\notin A(n_k,L)}|\hat\psi_\lambda(\xi)|^2
\leq \|V\|_{L^2(\T^2)}^2\sum_{\xi\notin A(n_k,L)}(|\xi|^2-\lambda)^{-2}.
\eeq
Because $\xi\notin A(n_k,L)$, we have, by definition,
$||\xi|^2-\lambda|\geq ||\xi|^2-n_k|-|n_k-n_{k+1}|>L-O(n_k^{o(1)})\gg n_k^\delta$
and, thus,
\beq
\|\psi_{\lambda,L}^R\|_{L^2(\T^2)}^2\leq \|V\|_{L^2(\T^2)}^2\sum_{||\xi|^2-\lambda|\gg n_k^\delta}(|\xi|^2-\lambda)^{-2}.
\eeq
This lattice sum can be bounded (cf. \cite{RU}, p. 11., eq. (5.9)) by a summation by parts which yields
\beq
\sum_{||\xi|^2-\lambda|\gg n_k^\delta}(|\xi|^2-\lambda)^{-2}\lesssim \frac{1}{L}+\frac{n_k^\theta}{L^2}
\eeq
where we recall $L=n_k^\delta$, so that $\alpha=\min(\delta,2\delta-\theta)=2\delta-\theta$, because $\delta<1/2-\theta<\theta$ (recall $\theta>1/4$), which gives $$\|\psi_{\lambda,L}^R\|_{L^2(\T^2)}^2\lesssim n_k^{-2\delta+\theta}.$$
Finally, $$R(\psi_\lambda,L)\lesssim n_k^{-\delta+\theta/2}.$$

\subsection{Bound on $\bra a\psi_{\lambda,L},\psi_{\lambda,L}\ket$}
Let us, first of all, prove the result for a monomial $a=e_\zeta$, $\zeta\in\cL\setminus\{0\}$.

We have
\beq
\bra e_\zeta \psi_{\lambda,L},\psi_{\lambda,L}\ket
=\sum_{\xi\in A(n_k,L)}\hat\psi_\lambda(\xi)\overline{\hat\psi_\lambda(\xi-\zeta)}
\eeq
and, thus, we obtain the bound
\beq
|\bra e_\zeta \psi_{\lambda,L},\psi_{\lambda,L}\ket|
\leq \|V\|_{L^2(\T^2)}\left(\sum_{\xi\in A(n_k,L)}(|\xi-\zeta|^2-\lambda)^{-2}\right)^{1/2},
\eeq
where we used
$$\sum_{\xi\in A(n_k,L)}|\hat\psi_\lambda(\xi)|^2\leq \|\psi_\lambda\|_{L^2(\T^2)}^2=1$$
as well as the bound \eqref{Fourier_bound}.

If $n_k\in \cQ'$, meaning that $A(n_k,L)$ is a good annulus, then, by definition, $\xi-\zeta\notin A(n_k,L)$, which yields the lower bound
$$||\xi-\zeta|^2-\lambda|\geq ||\xi-\zeta|^2-n_k|-O(n_k^{o(1)})\gg n_k^\delta.$$

And, hence, we obtain the bound
\beq
\begin{split}
|\bra e_\zeta \psi_{\lambda,L},\psi_{\lambda,L}\ket|
&\leq \|V\|_{L^2(\T^2)} n_k^{-\delta} \# A(n_k,L)^{1/2}\\
&\lesssim \|V\|_{L^2(\T^2)} n_k^{-\delta+\theta/2}.
\end{split}
\eeq

One may now readily generalize this bound for trigonometric polynomials of the form
$$a_\epsilon=\sum_{|\zeta|\leq n_k^\epsilon}\hat a(\zeta) e_\zeta,$$
namely,
\beq
\begin{split}
|\bra a_\epsilon \psi_{\lambda,L},\psi_{\lambda,L} \ket|
\leq& \sum_{|\zeta|\leq n_k^\epsilon}|\hat a(\zeta)||\bra e_\zeta \psi_{\lambda,L},\psi_{\lambda,L}\ket| \\
\lesssim& \|\hat{a}\|_{\ell^1} \|V\|_{L^2(\T^2)} n_k^{-\delta+\theta/2}.
\end{split}
\eeq

\subsection{General test functions}\label{sec-general}
It remains to generalize our argument for $C^\infty$ test functions.
We have
\beq
|\bra a\psi_{\lambda,L},\psi_{\lambda,L}\ket|
\leq |\bra a_\epsilon\psi_{\lambda,L},\psi_{\lambda,L}\ket|+|\bra(a-a_\epsilon)\psi_{\lambda,L},\psi_{\lambda,L}\ket|
\eeq
The first term is of order $O(n_k^{-\delta+\theta/2})$, as we saw above.

For the second term, we have
$$|\bra(a-a_\epsilon)\psi_{\lambda,L},\psi_{\lambda,L}\ket|
\leq \sum_{|\zeta|>n_k^\epsilon}|\hat a(\zeta)|
\lesssim_{a,K} \sum_{|\zeta|>n_k^\epsilon} |\zeta|^{-K}
\lesssim_K n_k^{\epsilon(2-K)}$$
for any $K>2$.

In summary, we have (choose $\delta=1/2-\theta-\epsilon$, for some $\epsilon>0$)
\beq
|\bra a\psi_\lambda,\psi_\lambda\ket|\lesssim_{a,\epsilon} n_k^{-\delta+\theta/2}=n_k^{(-1+3\theta)/2+\epsilon}.
\eeq
where the implied constant depends on a certain number of derivatives of $a$.

\section{Proof of Proposition \ref{good_ann}}\label{sec-good}
In this section we will prove the following proposition.
\begin{prop}\label{good}
Let $\delta\in(\theta/2,1/2-\theta)$ and $\epsilon<\tfrac{1}{2}-\theta-\delta$. There exists a full density subset $\cQ_1\subset \cQ$ such that for any $n\in \cQ_1$ we have
that $\xi\in A(n,n^\delta)$ implies that for any $\zeta\in\cL\setminus\{0\}$, $|\zeta|\leq n^\epsilon$, we have $||\xi+\zeta|^2-n|\gg n^\delta$.
\end{prop}
\begin{proof}
Let us introduce the set of ``bad'' lattice points
\beq
S_\zeta:=\{\xi\in\cL \mid |\bra\xi,\zeta\ket|\leq |\xi|^{2\delta}\}.
\eeq

We define
\beq
B_\zeta:=\{n\in S \mid  A(n,n^\delta)\cap S_\zeta\neq\emptyset\}
\eeq
and
\beq\label{S1}
\cQ_1=\bigcap_{\substack{\zeta\in\cL\setminus\{0\} \\ |\zeta|\leq n^\epsilon}}B_\zeta^c
\eeq

We will adapt an argument from \cite{RU} to show that $B_\zeta$ is of density zero in $\cQ$ in the sense that\footnote{We note that in \cite{RU} the bad set was defined as a subset of a given interlacing sequence of the Laplace eigenvalues. Here the bad set is a subset of the Laplace eigenvalues. But the proof works in an analogous fashion.}
\beq\label{bad_ev_bound}
\# \{n\in B_\zeta \mid n\leq X\} \lesssim \frac{X^{1/2+\theta+\delta}}{|\zeta|}
\eeq
where $\delta\in(\theta/2,1/2-\theta)$.

To see this, we first of all recall the following bound (cf. Lemma 6.2 in \cite{RU}) on the number of bad lattice vectors in a ball of radius $X^{1/2}$
\beq\label{badset_bound}
\#\{\eta\in S_\zeta \mid |\eta|^2\leq X\} \lesssim \frac{X^{1/2+\delta}}{|\zeta|}.
\eeq
We will now use this bound to obtain a bound on the number of elements $n\in B_\zeta$, $n\leq X$.

Let us define the set of distinct values of the quadratic form $Q(\xi)=|\xi|^2$, $\xi\in\cL$, which arise from bad lattice vectors. We denote
\beq
\cQ_\zeta:=\{q\in\R_+ \mid \exists\eta\in S_\zeta: q=|\eta|^2\}
\eeq

To this end, we define the projection map $\pi: B_\zeta\to \cQ_\zeta$ which assigns to $n\in B_\zeta$ the element of $\cQ_\zeta$ which is closest to $n$. If there are two such elements, then we define $\pi(n)$ as the larger of these two elements.

For a given $m\in \cQ_\zeta$ we have that the preimage satisfies
$$\pi^{-1}(m)\subset \{n\in B_\zeta \mid \exists\eta\in S_\zeta\cap A(n,n^\delta): |\eta|^2=m\}.$$
To see the inclusion, suppose $n\in\pi^{-1}(m)$. Recall that by definition of $\cQ_\zeta$ we have some $\eta_0\in S_\zeta$ s. t. $|\eta_0|^2=m$. Suppose for a contradiction that $\eta_0\notin S_\zeta\cap A(n,n^\delta)$, so $\eta_0\notin A(n,n^\delta)$, which means $||\eta_0|^2-n|>n^\delta$. On the other hand, $m=\pi(n)$ and so any other lattice vector $\eta\in S_\zeta$ must satisfy $||\eta|^2-n|\geq|m-n|>n^\delta$, by definition of the map $\pi$. This means $\eta\notin A(n,n^\delta)$, hence $S_\zeta\cap A(n,n^\delta)$ is empty, a contradiction to $n\in B_\zeta$. It follows that $\eta_0\in S_\zeta\cap A(n,n^\delta)$ which proves the inclusion.

Moreover, we have (for $m$ large enough)
\beq\label{incl}
\{n\in B_\zeta \mid \exists\eta\in S_\zeta\cap A(n,n^\delta): |\eta|^2=m\}\subset \cQ\cap[m-2m^\delta,m+2m^\delta]
\eeq
which easily follows, because for any $n\in B_\zeta$ in the set above we have $\eta\in A(n,n^\delta)$ with $|\eta|^2=m$. This means $|n-m|=|n-|\eta|^2|\leq n^\delta$ which implies $|n-m|\leq 2m^\delta$ for $m$ large enough.

We have
$$\#\{n\in B_\zeta \mid n\leq X\}= \sum_{m\in \cQ_\zeta}\#\pi^{-1}(m)\cap[0,X].$$

Now if $n\in\pi^{-1}(m)\cap[0,X]$, then $m=\pi(n)$ and $|n-m|\leq n^\delta$, because $S_\zeta\cap A(n,n^\delta)$ is non-empty, because $n\in B_\zeta$, and by definition $\pi(n)$ is the element of $\cQ_\zeta$ closest to $n$. So $m\leq n+n^\delta\leq X+X^\delta$, because $n\leq X$.

Hence,
$$\sum_{m\in\cQ_\zeta}\#\pi^{-1}(m)\cap[0,X] \leq \sum_{\substack{ m\in\cQ_\zeta \\ m\leq X+X^\delta}}\#\pi^{-1}(m)
\leq \sum_{\substack{ m\in\cQ_\zeta \\ m\leq X+X^\delta}}\#\cQ\cap [m-2m^\delta,m+2m^\delta]$$
where we used the inclusion \eqref{incl}.

The asymptotic $$\sum_{n\leq X} r_\cL(n)=\pi X+O(X^\theta)$$ yields (recall $\delta<\theta$)
$$\#\cQ\cap [m-2m^\delta,m+2m^\delta]\leq \sum_{|n-m|\leq 2m^\delta}r_\cL(n)\lesssim m^\delta+m^\theta\lesssim m^\theta.$$

Therefore,
$$\sum_{\substack{ m\in\cQ_\zeta \\ m\leq X+X^\delta}}\#\cQ\cap [m-2m^\delta,m+2m^\delta]
\lesssim X^\theta \sum_{\substack{ m\in\cQ_\zeta \\ m\leq X+X^\delta}}1\lesssim \frac{X^{1/2+\delta+\theta}}{|\zeta|}$$
where we used \eqref{badset_bound}.

This yields the bound \eqref{bad_ev_bound}.

Let us now use \eqref{bad_ev_bound} to prove that the complement of $\cQ_1$ is of density zero in $\cQ$:
\beq\label{compl-bound}
\begin{split}
&\#\{m\in \cQ_1^c \mid m\leq X\}\\
=&\#\{m\in S, m\leq X
\mid \exists \zeta\in\cL\setminus\{0\}, |\zeta|\leq m^\epsilon: A(m,m^\delta)\cap S_\zeta\neq\emptyset\}\\
\leq&\sum_{\substack{\zeta\in\cL\setminus\{0\} \\ |\zeta|\leq X^\epsilon}}
\#\{m\in B_\zeta \mid |m|\leq X\}\\
\lesssim& X^{1/2+\delta+\theta} \sum_{\substack{\zeta\in\cL\setminus\{0\}\\|\zeta|\leq X^\epsilon}}\frac{1}{|\zeta|}
\lesssim X^{1/2+\delta+\theta+\epsilon}=o(X).
\end{split}
\eeq

Let us now assume $n\in \cQ_1$ and $\xi\in A(n,n^\delta)$. We have to show that for any $\zeta\in\cL\setminus\{0\}$, $|\zeta|\leq n^\epsilon$, we have $||\xi+\zeta|^2-n|\gg n^\delta$.

Indeed,
\beq
\begin{split}
||\xi+\zeta|^2-n|&\geq 2|\bra\xi,\zeta\ket|-||\xi|^2-n|-|\zeta|^2\\
&\geq 2|\xi|^{2\delta}-n^\delta-O(n^{2\epsilon})\\
&\gg n^{\delta}
\end{split}
\eeq
because $\xi\in A(n,n^\delta)$.

\end{proof}


\section{\label{sec:Intro}Introduction}
The torsion pendulum, or torsion balance, is a sensitive device used to detect weak forces. \cite{Gil93} It was used in the experiments of Cavendish and Coulomb, conducted in the 18th century to measure gravitational and electrostatic constants. The torsion pendulum is still one of the most sensitive measurement methods in modern physics and engineering, with applications ranging from low-frequency gravitational wave detection \cite{Shi20} to studies of the mechanical properties of materials. \cite{Kes21}

The torsion pendulum has been used to detect small mechanical torques exerted on macroscopic objects from quantum mechanical spins. The first experiment was performed in 1915 by Einstein and de Haas, \cite{Ein15} who demonstrated that the mechanical rotation of a ferromagnetic material is caused by flipping of its magnetization or internal atomic spins. The second experiment, by Beth, was carried out in 1936 \cite{Beth36} and demonstrated that flipping the photon spin from $+\hbar$ to $-\hbar$ ($\hbar$: Planck's constant divided by $2\pi$) induced the rotation of a transmitting half-wave plate. The results of both experiments are understood in terms of angular momentum conservation of the system of atomic or photonic spins, and the rotation of the macroscopic object.
Spin-induced torque has received renewed attention in recent spintronics and optomechanics studies. Angular momentum transmission due to a spin wave has been mechanically detected with an yttrium ion garnet cantilever. \cite{Har19} Gigahertz rotation of optically levitated nanoparticles has been demonstrated via spin transfer arising from circularly polarized trapping light. \cite{Rei18, Jon18}

In our previous report, \cite{Hat19} we discussed the development of a torsion pendulum to study spin transfer from an optically spin-polarized atomic gas to a solid gas container through atom-surface interactions. The experiment was similar to that of Beth, in which an object suspended with a thin wire was irradiated from the bottom with circularly polarized light. Beth's experiment has been described in many textbooks. \cite{Hec, Joa} However, to our knowledge, the careful and patient measurements conducted therein have not been replicated using the original forced oscillation method, partly due to experimental difficulties including radiation pressure torque and the so-called radiometer effects, as well as the multiple steps involved in deriving the torque from the measured amplitude of the forced oscillation, as pointed out in Ref.~\cite{Del05,Emi18}. In Ref.~\cite{Del05}, the authors directly measured photon-spin-induced torque on the order of $10^{-12}$~N m from the acceleration of a suspended half-wave plate with a high-power infrared laser.

In this study, we analyzed our developed pendulum system, including measurement of the photon-spin-induced torque. Our experimental strategy was similar to the strategy applied in Beth's experiment: forced torsional oscillation induced by periodic application of external torque at the resonance frequency of the pendulum. The main differences between our system and that of Beth's study are the use of a double pendulum, in which the first pendulum works as a vibration isolator, and a higher operating frequency of 0.1~Hz relative to 0.002~Hz in Beth's experiment.
Through this simple passive isolation from external vibration noise, the torque sensitivity was $2 \times 10^{-17}$~N m (in 10$^4$ -s measurement time), which is close to the thermal noise limit and one order smaller than the minimum torque measured in Beth's experiment. The observed spin-induced torque exerted on the light-absorbing optics is consistent with the angular momentum transfer of $\hbar$ in the direction of light propagation per photon for left-circularly polarized light. \cite{Hec,Cor}
Therefore, we concluded that our torsion pendulum properly detected the photon-spin-induced torque.

This paper is organized as follows. We first explain the experimental apparatus, with special attention given to the design of the double pendulum in Sec.~\ref{sec:App}. In Sec.~\ref{sec:Exp}, we describe the experimental procedure, including a forced torsional oscillation approach and a lock-in data analysis method. We then present our experimental results in Sec.~\ref{sec:Res}, in which measurement of photon-spin-induced torque and estimation of the torque sensitivity  are discussed. The paper is concluded in Sec.~\ref{sec:Con}.

\section{\label{sec:App}Experimental apparatus}

\begin{figure}
\includegraphics[width=85mm]{Fig1}
\caption{\label{app}Schematic diagram of the experimental apparatus. PBS: polarization beam splitter; PSD: position-sensitive detector; $\lambda/4$: quarter-wave plate; $\lambda/2$: half-wave plate; LCVR: liquid crystal variable retarder.}
\end{figure}

The experimental apparatus, shown in Fig.~\ref{app}, is a slightly modified version of our previously reported one. \cite{Hat19}
The double torsion pendulum, discussed in detail below, was suspended in a vacuum chamber fixed on an optical table. The pressure inside was maintained at $1\times10^{-5}$~Pa with an ion pump, which was lower than the radiation pressure exerted on a suspended object by a typical torque laser beam used in the experiment ($1 \times10^{-4}$~Pa for a laser intensity of $3\times10^4~\rm{W/m^2}$). Thus, the radiometer effect from the background gas was negligible.

The angle of the suspended object (``optics holder'') was monitored using an optical lever method with a red probe laser (CPS650F; Thorlabs Inc.) and a position-sensitive detector (PSD) (S3932; Hamamatsu Photonics K.K.) located $L=0.38$~m from the torsional axis. The angle displacement $d\phi$ was derived from the position displacement $dx$ measured with the PSD as $d\phi=dx/(2L)$.
Light from the torque laser, a volume-holographic-grating-stabilized laser with a wavelength of 852~nm (LD852-SEV600; Thorlabs Inc.), was incident on the object from the bottom. The light polarization was controlled with a set of waveplates, including a liquid crystal variable retarder (LCVR) (LCC1223-B; Thorlabs Inc.), to periodically flip the polarization of the torque laser.
Signal processing and data recording were systematically performed with PC-based digital and analog input/output devices controlled by LabVIEW programs (NI).

Compared to Beth's experiment, which used a quartz fiber 0.25 m in length, the optics holder to be rotated was suspended with a short tungsten (W) wire approximately 5~mm in length and 10~$\mu$m in diameter. We refer to this portion of our double pendulum as the second pendulum.
The short wire resulted in a higher torsional resonance frequency of 0.1 Hz, compared to 0.002 Hz in Beth's experiment. To compensate for this, the second pendulum was mounted at the end of another pendulum, which had a resonance frequency of 0.0008 Hz to provide better isolation from external vibration noise, as discussed below.
The higher resonance frequency also provided shorter measurement times, which minimized spurious effects of drifting. The optics holder accommodated half-inch optics, on which torque was exerted from the torque laser beam; a mirror reflected the torque laser beam when it transmitted the optics. A cubic mirror was attached to the top of the optics holder for the optical lever. The moment of inertia of the optics holder was estimated to be $(9.4 \pm 0.9)\times10^{-8}$~kg m$^2$ from the drawing. The uncertainty originated from inaccuracy in the manufacturing process; it was estimated by comparing the designed and measured weights of the optics holder.

The optics holder suspended with the short wire was attached to another torsion pendulum (the first pendulum), which consisted of an object (the ``wheel'') with a large moment of inertia and a long wire. The moment of inertia of the wheel was estimated to be $7\times 10^{-6}$~kg m$^2$, and the W wire was approximately 1~m in length and 10~$\mu$m in diameter. The first pendulum had a low resonance frequency of torsional oscillation at 0.0008~Hz and served as a vibration isolator for the second pendulum, as discussed below.

The long wire was attached to a vacuum chamber flange (VRS-70M; Shinku-Kogaku Corp.), which was rotatable with an electric stepper motor (0.00004$^{\circ}$/pulse). This rotatable flange was used to dampen the torsional oscillation of the first pendulum. Figure \ref{damp} shows the angle of the rotatable flange controlled electronically, as well as the angle of the optics holder measured with the optical lever. The curve representing the angle of the optics holder was low-pass-filtered to clarify the oscillation of the first pendulum at its resonance period (1240~s). At 0~s in Fig. \ref{damp}, proportional integral control to the flange angle was activated toward the target angle of 0~mrad for the optics holder to dampen the resonance oscillation of the first pendulum. This damping technique was quite useful for readying the system for forced oscillation measurements for the second pendulum, after the apparatus experienced large external perturbations such as earthquakes, which often caused large and slowly damping torsional oscillation of the first pendulum.

\begin{figure}[h]
\includegraphics[width=85mm]{Fig2}
\caption{\label{damp}Damping of the torsional oscillation of the first pendulum with the rotatable flange. The angle of the rotatable flange was controlled electronically, while the angle of the optics holder was measured with the optical lever. The curve representing the angle of the optics holder was low-pass-filtered to clarify the oscillation of the first pendulum at its resonance period (1240~s). At 0 s, the automatic proportional integral control to the flange angle was activated toward the target angle of 0 mrad for the optics holder.}
\end{figure}

Figure \ref{sim} shows the Fourier transform spectrum of a free oscillation, recorded for $2\times 10^4$~s after damping of the first pendulum oscillation had been completed. The peak at 0.0008~Hz corresponds to the resonant torsional oscillation of the first pendulum. The amplitude decreased as the frequency increased beyond the resonance frequency. The peak at 0.1~Hz corresponds to the torsional oscillation of the second pendulum and is the focus of this study.

\begin{figure}[h]
\includegraphics[width=85mm]{Fig3}
\caption{\label{sim}Fourier transform spectrum of free torsional oscillation of the pendulum for $2\times10^4$-s measurement. Numerical simulations for the double and single pendulums are also shown.}
\end{figure}

To quantitatively examine the effect of vibration isolation with the first pendulum, we numerically solved the coupled equations that describe the torsional motions of the two pendulums with application of white noise between 0 and 2~Hz to the fixed end of the first pendulum wire; the freedom of motion in the simulation was the rotation of the first and second pendulums around the fixed wire axis. With appropriate experimental parameters, as shown in Fig.~\ref{sim}, the Fourier transform spectrum of the simulation reproduced the experimental observations fairly well, particularly the decrease in baseline oscillation amplitude at higher frequencies than the resonance of the first pendulum. The effect of vibration isolation can be clearly understood by comparing the measured spectrum with the single pendulum simulation performed on the assumption that the second pendulum was attached directly to the flange without the first pendulum.

\section{\label{sec:Exp}Experiments}
We adopted a forced oscillation method at the resonance frequency of the second pendulum.
A numerical simulation using the coupled equations for the double pendulum confirmed that the second pendulum could be regarded as a single pendulum around its resonance frequency, while the first pendulum (to which the second pendulum was attached) oscillated at a much lower resonance frequency. Therefore, we can use basic knowledge regarding the forced oscillation of the single pendulum, as described below.

In a steady state, the frequency dependence of the quadrature component of the forced oscillation with respect to the periodic external torque $N'\cos\omega t$ ($\omega$: angular frequency; $t$: time) corresponded to the following Lorentzian function around the resonance frequency:
\begin{equation}
\frac{\mu N'}{4I\omega_R}\frac{1}{(\omega-\omega_R)^2+(\mu/2)^2}, \label{quadrature}
\end{equation}
where $I$ is the moment of inertia, $\omega_R$ is the resonance angular frequency, and $\mu$ is the full width at half maximum of the resonance, determined by friction proportional to the angular velocity. Therefore, the peak amplitude of the resonance $A$ is $N'/(\mu I\omega_R)$.
The torque modulation $N'$ was then derived from the quadrature component of the torsional oscillation at the resonance, as follows:
\begin{equation}
N'=\mu I \omega_R A=\frac{I\omega_R^2A}{Q},\label{N’}
\end{equation}
where $Q=\omega_R/\mu$ is the Q-factor of the resonance.

In this study, the optical torque dependent on the photon spin was modulated in a square wave by changing the torque laser polarization with the LCVR, rather than in a sinusoidal wave. According to the Fourier expansion, $\pm N$ square modulation in the torque corresponds to $N'=\frac{4}{\pi} N$ sinusoidal modulation. Therefore, the factor $4/\pi$ was considered when deriving the torque modulation $N$ from the measured quadrature component as follows:
\begin{equation}
N=\frac{\pi}{4}\mu I \omega_R A=\frac{\pi}{4}\frac{I\omega_R^2A}{Q}.\label{N}
\end{equation}
We evaluated $\omega_R$ and $Q$ in every set of measurements to derive $N$ from the measured $A$.

The quadrature and in-phase components of the oscillation (i.e., the oscillation amplitude and phase delay with respect to the reference signal, which was synchronized with the modulation of the torque laser polarization) were derived by lock-in analysis of the forced oscillation using Igor Pro data analysis software (WaveMetrics, Inc.).

The spin-induced torque modulated by flipping the polarization of the torque laser between the left- and right- circular polarizations is ${P}_{\rm{laser}}\lambda/(2\pi c)= {P}_{\rm{laser}}[W]\times 4.52 \times 10^{-16}$~N m, where ${P}_{\rm{laser}}$ and $\lambda$ are the power and the wavelength of the torque laser, respectively, and $c$ is the speed of light, because one photon carries $+\hbar$ or $-\hbar$ spin component along the laser propagation.
Special precautions were required against torque caused by radiation pressure.
The torque exerted on the optics by the radiation pressure of the torque laser beam $\bm{N}_{\rm{rad}}$ can be calculated by the following equation:
\begin{equation}
\bm{N}_{\rm{rad}}= \bm{r} \times \bm{F}_{\rm{rad}},
\end{equation}
where $\bm{r}$ is the position at the laser beam enters the optics with respect to the wire-fixed position, and $\bm{F}_{\rm{rad}}$ is the radiation force.
If the torque laser beam is incident either perfectly parallel to or onto the wire-fixed point, the torque around the rotation axis is zero. Assuming imperfect incidence of the torque laser with a deviation of $a$ and an incline of $\alpha$ from the rotation axis, the torque around the rotation axis is $\alpha |\bm{F}_{\rm{rad}}| a=\alpha P_{\rm{laser}} a/c$.
When we measure optical torque in polarization modulation mode, the ratio of modulation of the radiation pressure torque to that of the spin-induced torque is $2\pi(a/\lambda)\alpha\beta$, where $\beta$ is a factor representing the modulation of the radiation pressure torque. This modulation may have been caused by variation in the laser power. Careful adjustment of the optical system to minimize the imperfect incidence and power variation associated with the polarization modulation is required to make the factor $2\pi(a/\lambda)\alpha\beta$ much less than $1$, such that the spin-induced torque modulation dominated the radiation pressure torque modulation. It was not difficult to fulfill this requirement, assuming $a=1$~mm, $\alpha=10^{-2}$~rad, and $\beta=10^{-3}$ for $\lambda=852$~nm. The effect of radiation pressure was examined by intentionally modulating the power of the torque laser in the experiment.

\section{\label{sec:Res}Results}
Because the optics holder was able to accommodate various types of half-inch optics, we examined several optics.
We mainly focused on a neutral density (ND) filter (NE560B-B; Thorlabs Inc.) as a light absorber, because our future target, a solid container of an atomic gas, can be regarded as a light absorber. \cite{Hat19} The absorption of the torque light by the ND filter was greater than 99.9\%. The spatial mode of the laser beam was filtered to be TEM$_{00}$ with a single-mode fiber. The laser power incident on the ND filter was 132~mW. The polarization of the torque laser was modulated between left-circular and right-circular polarization using the LCVR.

Figure \ref{ND_spc} shows the quadrature component of the induced torsional oscillation as a function of the modulation frequency around the resonance frequency of the second pendulum. One data point was the average of data collected for $10^4$~s in a steady state. The estimation of the error bar is discussed below. Lorentzian fitting yielded a peak amplitude of $0.78 \pm 0.08$~$\mu$rad at the resonance frequency of $0.08982 \pm 0.00002$~Hz, and a full width at half maximum value of $(3.1 \pm 0.5) \times 10^{-4}$~Hz. The Q-factor was derived as $(2.9 \pm 0.5)\times 10^2$. Therefore, the measured torque modulation was $(6.3 \pm 1.3)\times 10^{-17}$~N m, derived using Eq.~(\ref{N}). The systematic uncertainty, mainly originating from conversion of the position signal measured with the PSD to the angle of the optics holder, was less than 5\%, smaller than this statistical uncertainty of $1.3\times 10^{-17}$~N m. The measured torque was consistent with the theoretically estimated value of $(6.0 \pm 0.2) \times 10^{-17}$~N~m derived from the assumption that one photon transfers $\hbar$ angular momentum to the ND filter; the main origin of the uncertainty in the theoretical estimation originated from the laser power.

\begin{figure}[h]
\includegraphics[width=85mm]{Fig4}
\caption{\label{ND_spc}Resonance spectrum of the second pendulum with the ND filter irradiated with counterrotating circularly polarized light. Each error bar was estimated to be $\pm0.1~\mu$rad (see the text).}
\end{figure}

The error bars in Fig. \ref{ND_spc} were estimated from different measurements. We obtained 10 sets of $10^4$-s free torsional oscillation data for the pendulum, and applied lock-in analysis to them with imaginary reference signals at the resonance frequency. Figure~\ref{Noise} shows the obtained quadrature components of the oscillation with the standard deviations. Each set of data had a standard deviation similar to the case of the measurements in Fig. \ref{ND_spc}. The 10 measurements resulted in an average quadrature component of $0.11~\rm{\mu rad}$ and a standard deviation of $0.11~\rm{\mu rad}$. We therefore estimated that the uncertainty in detection of the quadrature component was $0.1~\rm{\mu rad}$ in a single $10^4$-s measurement; we adopted this value as the error bar for each data point shown in Fig. \ref{ND_spc}.

\begin{figure}[h]
\includegraphics[width=85mm]{Fig5}
\caption{\label{Noise} Fluctuation of quadrature components of the free torsional oscillation. Each error bar represents the standard deviation for each $10^4$-s measurement.}
\end{figure}

For the second pendulum used to estimate the uncertainty mentioned above, we found that $\omega_R=0.110317\pm0.00007$~Hz and $Q=173\pm5$. Using Eq.~(\ref{N}) with these parameters, the uncertainty in the quadrature component detection of $0.1~\rm{\mu rad}$ was converted to the uncertainty in the torque detection of $2 \times10^{-17}$~N m. We regarded this as the torque sensitivity of our system over a measurement time of 10$^4$ s.
The thermal noise limit of the torque $N_{\rm{ther}}$ can be estimated using the following formula, \cite{Hai07}
\begin{equation}
N_{\rm{ther}}=\frac{\pi}{4}\sqrt{\frac{4k_BTI\omega_R}{Q\Delta t}} \label{min_torque},
\end{equation}
where $k_B$ is the Boltzmann constant, $T$ is temperature, and $\Delta t$ is the measurement time. The estimated thermal noise limit was $2 \times10^{-17}$~N m, which is close to the experimentally estimated sensitivity.

 We also confirmed that the observed oscillation phase was consistent with the expected rotation direction; specifically, left-circularly polarized light carries $+\hbar$ angular momentum along the direction of propagation \cite{Hec,Cor} and exerts a positive torque on the ND filter.

To confirm that the observed torque originated from photon spin, rather than from other effects such as radiation pressure torque, we measured the light polarization dependence. At a fixed modulation frequency, we changed the polarization modulation by rotating the LCVR axis from a pair of counterrotating circular polarizations to a pair of counterrotating elliptical polarizations. Figure \ref{ND_po} shows the measured torque as a function of the spin component along the light propagation direction. The spin component of $0~\hbar$ corresponds to nominally unmodulated linear polarization, whereas the spin component of $1~\hbar$ corresponds to circular polarizations, with elliptical polarizations between the two extremes. The observed torque is reasonably proportional to the photon spin component of the propagation direction.

\begin{figure}[h]
\includegraphics[width=85mm]{Fig6}
\caption{\label{ND_po} Polarization dependence of the photon-spin-induced torque. Each error bar was estimated to be $\pm2 \times 10^{-17}$~N~m (see the text).}
\end{figure}

We also modulated the power of the torque laser without changing the polarization, to estimate the torque modulation originating from the radiation pressure; specifically, a $\pm$7\% change in power induced a torque modulation of $1 \times10^{-16}$~N~m. The power modulation associated with the polarization change with the LCVR was 0.1\%. Therefore, we estimated that the torque modulation originating from the radiation pressure was $1\times10^{-18}$~N~m, which was negligible in the observed torque modulation.
Taken together, we concluded that the observed torque was induced by photon spin transfer.

\section{\label{sec:Con}Conclusions}
In this paper, we describe the characterization of our developed double torsion pendulum system using optical torque with a forced oscillation method similar to Beth's experiment. In the double pendulum approach, the first pendulum operates as a vibration isolator. With this simple passive isolation from external vibration noise, the torque sensitivity reached $2 \times 10^{-17}$~N m in 10$^4$ -s measurement time, which is close to the thermal noise limit and one order smaller than the minimum torque measured in Beth's experiment.
We evaluated the photon-spin-induced torque using a ND filter as a light absorber. The observed torque was consistent with the torque expected from the angular momentum transfer of $\hbar$ (or $-\hbar$) per photon in left- (or right-) circularly polarized light. We concluded that our torsion pendulum properly detected the photon-spin-induced torque. Our demonstration can be regarded as a recapitulation of Beth's experiment.

Finally, the developed double pendulum could be used in future detailed studies of spin transfer from gases to solids.

\begin{acknowledgments}
We thank Jean Jeener for his useful comments.
This study was supported by the Matsuo Foundation.
\end{acknowledgments}

\section*{Data Availability Statement}
The data that support the findings of this study are available from the corresponding author upon reasonable request.

\providecommand{\noopsort}[1]{}\providecommand{\singleletter}[1]{#1}%


\section{Introduction}

The fifth-generation (5G) wireless networks are under deployment over the world. As one of the key technologies in 5G, massive multiple-input-multiple-output (mMIMO)  makes use of a large number of antennas (e.g., $\geq 64$) to increase the spectral efficiency (SE) by, at least, an order of magnitude compared to 4G \cite{marzetta:10a}. Despite this, both academia and industry continue their efforts to allow meeting ever more stringent requirements, such as ultra high data rate and energy efficiency, global coverage and connectivity, as well as extremely high reliability and low latency \cite{Saad2020Vision}. To solve these challenges, different avenues are considered, such as to further densify the network, using even ultra-massive MIMO arrays or utilizing the large and available bandwidth of higher frequency bands such as (sub-) terahertz frequencies \cite{Boccardi2014Five} \cite{Rajatheva2020White}. However, all these methods lead to higher energy consumption. Moreover, one of the major weaknesses for higher frequency bands is the limited network coverage, which is due to the high penetration loss and the lack of scattering and diffraction.

In parallel, a novel concept named reconfigurable intelligent surface has received more and more attention from academia and industry \cite{Qingqing2020Towards}. It is also known as intelligent reflecting surfaces, large intelligent surfaces or software-controlled metasurfaces. Briefly speaking, RIS is a two-dimensional (i.e., with near-zero thickness) surface consisting of a large array of passive scattering elements, each of which can be controlled in a software-defined manner to impose an amplitude change and/or phase shift to the incident signals. By tuning the reflecting phases of the incident signals independently, the reflected signals can be added coherently to improve the received signal power or combined destructively to mitigate interference when arriving at the user equipment (UE) \cite{Shiming2020Toward}. In general, RIS is a prospective paradigm to realize a smart and reconfigurable wireless environment and improve the transmission performance with low energy consumption and hardware cost.

So far, numerous investigations have been performed on different aspects of RIS, such as prototyping \cite{Tang2019Wireless}\cite{Dai2020} \cite{Pei2021Ris-aided}, channel estimation  \cite{Zheng2020Intelligent}, promising application scenarios \cite{Liang2019Large}, etc. Among them, the modeling of received power for RIS-aided wireless communication systems is a fundamental yet challenging work, which is of great significance to the understanding of the RIS technology itself. Although several received power models have been proposed, most of them are lacking in experimental validations. The authors of \cite{Emil2020Power} make use of electromagnetic arguments to derive an exact channel gain expression that captures three essential near-field behaviors.
\cite{Ozdogan2020Intelligent} takes advantage of the physical optics techniques to derive the far-field pathloss model that is inversely proportional to the square of the product of distances. In \cite{Ellingson2019Path}, the received signal power is formulated as the sum of each scattering element's received power which is related to the distances between the RIS elements and the transceivers, the element gains, and the complex-valued coefficients of the elements. The authors in \cite{Tang2021Wireless} develop a free-space pathloss model for RIS and experimental validations were also performed. Their proposed pathloss model is related to the distances from transceivers to the RIS, the size of the RIS, and the radiation patterns of antennas and elements.

In this paper, we propose a novel radar cross section (RCS)-based received power model.
To the best of our knowledge, our model is the first that takes angles of incidence/reflection into account when describing the phase shifts of the reflection elements. More specifically, the reflection phase in our model is not only a function of the control voltage imposed on the RIS element, but also a function of the angles. This is an important fact to consider when calculating the desired reflection coefficients.
We validate the proposed model on our hardware platform and show that the simulation results based on our proposed model is very close to measurements.

\section{System Model}
We consider an RIS-aided wireless communication system as shown in Fig.~\ref{Geometric}. For simplicity, the transmitter and receiver are both equipped with a single antenna. However, the results of this paper can be easily generalized to the case of MIMO.

Since the model of line of sight (LOS) propagation  between the transmitter and the receiver is well known, we focus on modeling the signals reflected by the RIS in this paper only. In case the LOS path exists, we just need to add the signal of the LOS path and the reflected paths.
Referring to Fig.~\ref{Geometric}, the RIS with $M \times N$ elements is placed on the $x$-$y$ plane of a Cartesian coordinate system and each row is parallel to the $x$-axis. The origin of the coordinate system is located at the geometric center of the RIS. Let $\text{E}_{m,n}$ denote the element in the $m$-th row and $n$-th column. The length and width of an element are denoted by $dx$ and $dy$, respectively. The location of $\text{E}_{m,n}$ is expressed as
\begin{equation}
    \label{eq:pmn}
    \boldsymbol{p}_{m,n}=\left( \left( n-\frac{N+1}{2} \right)dx, \left( \frac{M+1}{2}-m \right)dy, 0 \right)
\end{equation}
in the Cartesian coordinate system.

Furthermore, we use the spherical coordinates to represent the positions of the transmitter/receiver. Specifically, $\left(d_{1},\theta_{1},\varphi_{1}\right)$ and $\left(d_{2},\theta_{2},\varphi_{2}\right)$ denote respectively the positions of the transmitter and receiver seen from the geometric center of the RIS plane, where $d,\theta,\varphi$ represent the distance, zenith angle and azimuth angle, respectively. If seen from $\text{E}_{m,n}$, the positions of the transmitter and receiver can be expressed as $\left(d_{m,n}^{t},\theta_{m,n}^{t},\varphi_{m,n}^{t}\right)$ and  $\left(d_{m,n}^{r},\theta_{m,n}^{r},\varphi_{m,n}^{r}\right)$. The spherical coordinates can be transformed into Cartesian coordinates by
\begin{equation}
	\label{eq:shptocar}
	\left\{
		\begin{array}{l}
			C_x = d \sin \theta \cos \varphi,\\
			C_y = d \sin \theta \sin \varphi,\\
			C_z = d \cos \theta.
		\end{array}
	\right.
\end{equation}

Denote the transmit signal as $x$. The received signal is written as
\begin{equation}
	\label{eq:received_signal}
	y = \left(\sum_{m=1}^{M} \sum_{n=1}^{N} h_{m,n} g_{m,n}  \Gamma_{m,n}\right) x + z,
\end{equation}
where $h_{m,n}$ is the channel coefficient between the transmitter and the reflecting element $\text{E}_{m,n}$, and $g_{m,n}$ is the channel coefficient between $\text{E}_{m,n}$ and the receiver. $ z $ is the noise. $\Gamma_{m,n}$ is the reflection coefficient of the element $\text{E}_{m,n}$.
\begin{figure}[!htbp]
	\centering
	\includegraphics[width=1\linewidth]{RIS2.eps}
	\caption{The geometric illustration of an RIS-aided communication system.}
	\label{Geometric}
\end{figure}

\section{The Received Power Model}
In this section, we propose a general model for the received power in the RIS-aided communication system.
Denote the transmit power by $P_{t}$ and transmit antenna gain by $G_{t}\left(\textbf{r}_{m,n}^{t}\right)$, where $\textbf{r}_{m,n}^{t}$ represents the direction from the transmitter to the element $\text{E}_{m,n}$ and is related to $\theta_{m,n}^{t}$ and $\varphi_{m,n}^{t}$. Then the spatial power density generated by the transmitter on the element $\text{E}_{m,n}$ is
\begin{equation}
	\label{eq:Smnt}
	S_{m,n}^{t}=P_{t} G_{t}\left(\textbf{r}_{m,n}^{t}\right) \frac{1}{4 \pi \left( {d_{m,n}^{t}}\right) ^{2}},
\end{equation}
where $d_{m,n}^{t}$ represents the distance from the transmitter to $\text{E}_{m,n}$. With the transmitter located at $\left(C_{x},C_{y},C_{z}\right)$, $d_{m,n}^{t}$ can be calculated by
\begin{equation}
	\label{eq:rmnt}
	\resizebox{.86\hsize}{!}{$d_{m,n}^{t}=\sqrt{\left[ C_{x}\!-\! \left( n\!-\!\frac{N\!+\!1}{2} \right)\!dx \right] ^{2} + \left[ C_{y}\!-\! \left( \frac{\!M\!+1}{2}\!-\!m \right)\!dy \right] ^{2} +C_{z}^{2}}.$}
\end{equation}

The power captured by $\text{E}_{m,n}$ is
\begin{equation}
	\label{eq:Pmnt}
	P_{m,n}^{t}=S_{m,n}^{t}A_{e}\left(-\textbf{r}_{m,n}^{t}\right),
\end{equation}
where $A_{e}\left(-\textbf{r}_{m,n}^{t}\right)$ is the effective area of the $m,n$-th element in the direction $-\textbf{r}_{m,n}^{t}$. Similar to \eqref{eq:Smnt}, the spatial power density generated by the $m,n$-th element is
\begin{equation}
	\label{eq:Smnr}
	S_{m,n}^{r}=P_{m,n}^{t} G_{e}\left(\textbf{r}_{m,n}^{r}\right) \frac{1}{4 \pi \left( {d_{m,n}^{r}}\right) ^{2}},
\end{equation}
where $d_{m,n}^{r}$ represents the distance from $\text{E}_{m,n}$ to the receiver, which can be calculated similarly to \eqref{eq:rmnt}. $G_{e}\left(\textbf{r}_{m,n}^{r}\right)$ is the gain of the element $\text{E}_{m,n}$ in the direction $\textbf{r}_{m,n}^{r}$. The relation between the effective area and element gain is generally expressed as \cite{Gustafsson2019Maximum}
\begin{equation}
	\label{eq:G_A}
	G_{e}\left( \textbf{r}_{m,n}^{r} \right) = \frac{4 \pi }{\lambda^2} A_{e}\left( \textbf{r}_{m,n}^{r}\right).
\end{equation}
Note that we do not consider the power loss of the reflection in \eqref{eq:G_A}. However, the efficiency of the reflection will later be taken into account in the reflection coefficient.
The power captured by the receiver that is reflected by the $m,n$-th element is
\begin{equation}
	\label{eq:Pmnr}
	P_{m,n}^{r}=S_{m,n}^{r}A_{r}\left(-\textbf{r}_{m,n}^{r}\right),
\end{equation}
where $A_{r}\left(-\textbf{r}_{m,n}^{r}\right)$ is the effective area of the receiving antenna in the direction $-\textbf{r}_{m,n}^{r}$, which has the similar relation with antenna gain $G_{r}\left(-\textbf{r}_{m,n}^{r}\right)$, namely
\begin{equation}
	\label{eq:A_G}
	A_{r}\left( -\textbf{r}_{m,n}^{r} \right) = \frac{\lambda^2}{4 \pi \eta_{r}} G_{r}\left(-\textbf{r}_{m,n}^{r} \right),
\end{equation}
where $\eta_{r} \in \left[ 0,1\right] $ is the efficiency of the receiving antenna. According to \eqref{eq:Smnt} $\thicksim$ \eqref{eq:A_G}, the received power reflected from the $m,n$-th element is
\begin{equation}
	P_{m,n}^{r}=P_{t} G_{t}\left(\textbf{r}_{m,n}^{t}\right)G_{r}\left(-\textbf{r}_{m,n}^{r} \right) \frac{A_{e}\left(-\textbf{r}_{m,n}^{t}\right) A_{e}\left( \textbf{r}_{m,n}^{r}\right)}{\left( 4 \pi d_{m,n}^{t} d_{m,n}^{r} \right)^{2}\eta_{r}}.
\end{equation}

Some recent experiments in \cite{Pei2021Ris-aided} have indicated that both the amplitude and phase shift of the reflection coefficient are related to the angle of incidence, the angle of reflection, and the control signal, e.g., the bias voltage of a varactor. An example is that, for a given control signal, the reflection coefficient also varies with the angles of arrival and departure.  As a result, we model the reflection coefficient as a function of the three factors mentioned above. For the element $\text{E}_{m,n}$, the reflection coefficient is modeled as
\begin{align}
	\label{eq:Gamma_mn}
	& \Gamma_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right)  \\
	& = \mu_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right) \mathrm{e}^{\text{j} \phi_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right)}, \nonumber
\end{align}
where $\mu_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right)$, $\phi_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right)$ represent the amplitude  attenuation and phase shift, respectively. $u_{m,n}$ is the control signal imposed on the element $\text{E}_{m,n}$.

The total received power reflected by the RIS is then written as
\begin{equation}
	\label{eq:Pr}
	P_{r} = \left| \sum_{m=1}^{M} \sum_{n=1}^{N} \Gamma_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right)  \sqrt{P_{m,n}^{r}} \mathrm{e}^{\text{j} \Phi_{m,n}} \right| ^{2},
\end{equation}
where $\Phi_{m,n}$ is the phase difference of the $m,n$-th element caused by the propagation delay:
\begin{equation}
    \Phi_{m,n} =\frac{2\pi}{\lambda}\left(d_{m,n}^{t}+d_{m,n}^{r} \right).
\end{equation}

In order to make the expression more intuitive, the concept of radar cross section (RCS) $\sigma_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right)$ is introduced to describe the joint impact of the amplitude attenuation and the effective area of the element $\text{E}_{m,n}$, namely
\begin{align}
		\label{eq:sigma_mn}
		&\sigma_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right)  \\
		& = \mu_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right) \sqrt{A_{e}\left(-\textbf{r}_{m,n}^{t}\right) A_{e}\left( \textbf{r}_{m,n}^{r}\right)}. \nonumber
\end{align}
The proposed received power model can now be written as
\begin{equation}
    \label{eq:Pr_appr}
	P_{r} =  \frac{P_{t}}{16 \pi^{2} \eta_{r}} \left| \sum_{m=1}^{M} \sum_{n=1}^{N} \frac{\sqrt{G_{t} G_{r}} \sigma_{m,n} \mathrm{e}^{\text{j} \left(\phi_{m,n} + \Phi_{m,n}\right)}}{d_{m,n}^{t}d_{m,n}^{r}} \right| ^{2}.
\end{equation}
Note that we temporarily drop the arguments of $G_{t}$, $G_{r}$, $\sigma_{m,n}$ and $\Phi_{m,n}$ in order to make the above expression more concise.

We now study the impact of angles on the RCS, which depends on multiple parameters, e.g., the incident/reflection angles and the efficiency of reflection. In practice, such a multi-dimensional data can be obtained by electromagnetic simulate software.
For simplicity however, we let the angles of $\theta^{t}$ and $\varphi^{t}$ both equal to zero degree(also the case in our experiments), and characterize the impact of zenith angle of reflection $\theta^{r}$ on the RCS $\hat{\sigma}$.
we adopt the model of \cite{Meyers2003Encyclopedia} here:
\begin{equation}
    \label{eq:RCS_fit}
	\hat{\sigma}\left(\theta^{r}\right) = \frac{4\pi A^{2}}{\lambda^{2}} \left[ \frac{\sin{\left( k\sqrt{A}\sin{(\theta^{r})}\right)}}{k\sqrt{A}\sin{(\theta^{r})}} \right]^{2} + c,
\end{equation}
where $A = dxdy$ is the geometry area of each element and $k=\frac{2\pi}{\lambda}$ is the wave number. $c$ is a constant which is related to other parameters, e.g., the reflection efficiency and/or the angle $\varphi^{r} $.

We now characterize the impact of  $\theta^{r}$ on the phase shift. We propose to model such a relationship with a cosine function:
\begin{equation}
    \label{eq:phi_fit}
	\hat{\phi}\left(\theta^{r}\right) = a \cos{(\theta^{r})} + b,
\end{equation}
where $a,b$ are the design parameters. An example of their values is given below.

Fig.~\ref{fit_curve} illustrates the reflection coefficient as a function of $\theta^{r}$ with azimuth angle $\varphi^{r}$ fixed to $45^{\circ}$. We can observe that our models fit the CST simulations well.
\begin{figure}[!htbp]
    \centering
    \subfigure[]{
	\label{RCS_fit}
d	\includegraphics[width=0.45\linewidth]	{RCS_fit.eps}}
	\subfigure[]{
		\label{phase_fit}
		\includegraphics[width=0.45\linewidth] {phase_fit.eps}}
    \caption{The reflection coefficient versus the zenith angle $\theta^{r}$ with $\phi = 45^{\circ}$. (a) $\hat{\sigma}$ versus $\theta^{r}$ with $c\!=\!1.42\!\times\!10^{-5}~\text{m}^{2}$. (b) $\hat{\phi}$ versus $\theta^{r}$ with $a\!=\!90^{\circ},b\!=\!180^{\circ}$. }
    \label{fit_curve}
\end{figure}

\section{Measurements and Validations}
In this section, we validate the proposed model with measurements and compare it with an existing model.
Some recent works have modeled the role of RIS as the specular reflection approximately when in the near field or the RIS is infinitely large \cite{Basar2019Wireless,Tang2021Wireless}. It means that the RIS acts as a plane mirror that reflects the incident signal with the same amplitude attenuation coefficient and perfect phase shift. Similar to the Friis transmission equation, the received power of the specular reflection model is expressed as
\begin{equation}
	\label{eq:nearfield_Pr}
	P_{r} = P_{t} G_{t} G_{r} \left[ \frac{\lambda \bar{\mu}}{4 \pi \left( d_{1} + d_{2} \right)} \right]^{2},
\end{equation}
where $\bar{\mu}$ is the average amplitude of the reflection coefficient
\begin{equation}
	\label{eq:barmu}
	\bar{\mu} = \frac{1}{MN}\sum_{m=1}^{M} \sum_{n=1}^{N} \mu_{m,n}\left( -\textbf{r}_{m,n}^{t}, \textbf{r}_{m,n}^{r}, u_{m,n}\right),
\end{equation}
and $G_{t}, G_{r}$ are the antenna gains.

From \eqref{eq:nearfield_Pr}, one may observe that the received signal power is inversely proportional to $(d_{1} + d_{2})^2$. In the following we will validate this model and our proposed model \eqref{eq:Pr_appr}.

The basic parameters of our experiments are listed in Table~\ref{basic}. Other details of our experimental platform can be found in \cite{Pei2021Ris-aided}.
\begin{table}[!htbp]
	\renewcommand\arraystretch{1.2}
	\caption{Basic Parameters Of the experiments}
	\begin{center}
		\begin{tabular}{c|c}
			\hline
			\textbf{Parameters}& \textbf{Value} \\
			\hline
			Operating frequency & $f = 5.8$ GHz  \\
			Antenna peak gain & $G_{t} = G_{r} = 17.1$ dB  \\
			Efficiency of antenna & $\eta_{r}  \approx 54.29 \%$  \\
			Rows and columns & $M = 20, N=55$  \\
			Element size & $dx = 14.3$ mm, $dy = 10.27$ mm  \\
			\hline
		\end{tabular}
		\label{basic}
	\end{center}
\end{table}

Fig.~\ref{Environments} illustrates the indoor and outdoor measurement environments.
The reflection coefficients are 1-bit quantized, i.e., the two control states have a phase difference of around 180 degrees, which is enabled by adjusting the bias voltage $u_{m,n}$ imposed to the varactor diodes. The reflection coefficients are adjusted according to Algorithm 1 in \cite{Pei2021Ris-aided}.

\begin{figure}[!htbp]
	\centering
	\subfigure[]{
	\label{nearfield}
	\includegraphics[width=0.95\linewidth]	{nearfield.eps}}
	\subfigure[]{
		\label{farfield}
		\includegraphics[width=0.95\linewidth] {farfield.eps}}
	\caption{Measurement environments. (a) Indoor near field testing scene. \\
	(b) Outdoor far field testing scene.}
	\label{Environments}
\end{figure}

We carry on near-field experiments indoor and far-field experiments outdoor. The boundary of near field and far field is defined as \cite{Tang2021Wireless}
\begin{equation}
	\label{eq:df}
	df = \frac{2MNdxdy}{\lambda^{2}}.
\end{equation}

For the RIS in our tests, the boundary $df \approx 6~\text{m}$. Thus, we select the distances from $1~\text{m}$ to $5~\text{m}$ in indoor near-field experiments and from $5~\text{m}$ to $50~\text{m}$ in outdoor far-field experiments. The antennas are placed in a way that the azimuth angles of the transmit/receive antennas seen from the origin of the RIS are always equal to zero degree(i.e., $\varphi_{1}=\varphi_{2}=0$).
\begin{figure}[!htbp]
	\centering
	\includegraphics[width=1\linewidth]{attenuation_distance_nearfield_d2_2.eps}
	\caption{Power attenuation versus the TX-RIS distance $d_{1}$ with $d_{2}=2, \theta_{1}=0, \theta_{2}=30^{\circ}, \varphi_{1}=0, \varphi_{2}=0$. }
	\label{nearfield_distance}
\end{figure}

Fig.~\ref{nearfield_distance} illustrates the indoor near-field measurement results and the model results of the power attenuation $P_{r}/P_{t}$ as a function of the distance between the transmitter and the geometric center of the RIS (i.e., $d_{1}$).
We can observe that the measurements have the same trend with our proposed RX power model, however with a gap of about $2$ dB. This difference might be due to the reflection of indoor metal objects and walls that makes the measurement results higher.
When it comes to the curve of the specular reflection model, it is also close to the measurement results in such a near-field case. These measurements validate the applicability of the received power model we proposed and the specular reflection model in this case.

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=1\linewidth]{attenuation_distance_farfield_d2_2.eps}
	\caption{Power attenuation versus the TX-RIS distance $d_{1}$ with $d_{2}=2, \theta_{1}=0, \theta_{2}=30^{\circ}, \varphi_{1}=0, \varphi_{2}=0$. }
	\label{farfield_distance}
\end{figure}

Fig.~\ref{farfield_distance} illustrates the outdoor far-field results of the power attenuation. The horizontal axis represents the distance between the transmitter and the geometric center of the RIS (i.e., $d_{1}$).
We can observe that there is an increasing gap between the specular reflection model and the measurement results. However, simulation results based on our model are consistent with the measurements, which validates the accuracy of the proposed received power model in far-field case.

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=1\linewidth]{attenuation_angle_nearfield_d2_2.eps}
	\caption{Power attenuation versus the zenith angle $\theta_{2}$ of RX with $d_{1}=3, d_{2}=2,\theta_{1}=0, \varphi_{1}=0, \varphi_{2}=0$. }
	\label{nearfield_angle}
\end{figure}

Fig.~\ref{nearfield_angle} illustrates the indoor power attenuation results as a function of the zenith angle of the receiver (i.e., $\theta_{2}$). We observe that the RX received power model we proposed matches the measurements with a gap of about $2$ dB. Since the specular reflection model does not consider the influence of angles, the received power model is a constant value, which is not aligned with the measured results.
\section{Conclusion}
In this work, we have proposed a novel RCS-based received power model for an RIS-aided wireless communication system. This model introduces an angular-dependent phase shift behavior of the RIS elements, which has not yet been considered in previous models. We then validated it through indoor and outdoor measurements.
The results show that the specular reflection model is able to describe the received power in near-field case when the incident/reflection angles are relatively small. However it fails when in the far-field case or in case the angles are large. By contrast, our model can hold for both near-field and far-field cases and it can characterize the impact of the angles well.

\section*{Acknowledgment}
This work was supported by the National Natural Science
Foundation of China under Grant 62071191. The work of D. Gesbert was partially funded via the HUAWEI France supported Chair on Future Wireless Networks at EURECOM.

\bibliographystyle{IEEEtran}


\section{Introduction}
 Two large common Monte Carlo (MC) simulation samples were produced for the study of the ILD detector concept, which was established in 2009~\cite{ILD}.
One was used for the ILC physics case study as ILC DBD~\cite{DBD} with center of mass (CM) energies of 250 GeV, 350 GeV, 500 GeV, and 1 TeV. The other was produced for the ILD detector optimization study, ILD IDR~\cite{IDR},  using the DBD generator samples at the CM energy of 500 GeV in 2018.\\
~~Since generator, simulation and reconstruction tools are well advanced since the DBD era,
a new full set of high statistics MC samples with the latest software, named MC-2020,  are requested from the physics working group for precise studies at the CM energy of 250 GeV.
Table~\ref{tab:request_statistics} shows the summary of requested statistics of MC samples from ILD physics working group.
In this request, $1~\mathrm{ab^{-1}}$ or $5~\mathrm{ ab^{-1}}$ are required for most of channels and $100~\mathrm{ k~events}$ for each of individual higgs decay channels. It is more than an order larger than the previously produced samples.
In order to meet the request of this huge production, we employ the DIRAC~\cite{DIRAC} and ILCDirac~\cite{ILCDIRAC} distributed computing system for mass production which was fully utilized in the ILD group since IDR production period.
DIRAC is originally developed for the LHCb experiment~\cite{LHCb} providing high level interface for users to use heterogeneous GRID resources, composed of job management, file catalog, job monitoring, control via web interface. ILCDirac is an extension of the DIRAC for linear collider experiments which is developed and maintained by the CLICdp group~\cite{CLIC}.

\begin{table}[htbp]
\begin{center}
\caption{Requested integrated luminosity or number of events for each process. Here
($e_{L}.p_{R}$, $e_{R}.p_{L}$, $e_{L}.p_{L}$, $e_{R}.p_{R}$) denote 100\% polarized (Left/Right, Right/Left, Left/Left, Right/Right)-handed electron/positron beams.}
\label{tab:request_statistics}
\begin{tabular}{|c|c|}
\hline
Process & Statistics \\
\hline
\hline
\multirow{2}{*}{$2f\to \ell\ell,~ qq$} & $e_L.p_R/e_R.p_L$ : $\mathrm{ 5~ab^{-1}}$\\
& $e_L.p_L/e_R.p_R$ : $\mathrm{1~ab^{-1}}$\\
\hline
\multirow{2}{*}{All $4f$} & $e_L.p_R/e_R.p_L$ : $\mathrm{ 5~ab^{-1}}$\\
& $e_L.p_L/e_R.p_R$ : $\mathrm{ 1~ab^{-1}}$\\
\hline
All $6f$ & $\mathrm{ 10~k~events}$\\
\hline
$e\gamma$, $\gamma e$, $\gamma\gamma$ process ($3f$, $5f$, $\gamma\gamma\to 2f$, $\gamma\gamma \to 4f$) & $\mathrm{1~ab^{-1}}$\\
\hline
$2f \to ee\gamma$ & $\mathrm{1~ab^{-1}}$\\
\hline
$h\to$ inclusive & $\mathrm{1~ab^{-1}}$\\
\hline
$h\to$ each decay ($5\times 9$ channels) & $\mathrm{100~k~events}$\\
\hline
$Z \to qq$, $Zh\to \nu\nu qq$ for LCFI & $\mathrm{50~k~events}$\\
\hline
$Z\to qq$ (91 GeV) for LCFI & $\mathrm{50~k~events}$\\
\hline
\end{tabular}
\end{center}
\end{table}

For MC-2020, ILD generator group produces the event samples using the latest Whizard 2.8.5~\cite{Whizard}.
In this production, high statistics generator samples are newly produced even with small cross section channels.
With the latest Whizard, generator outputs are directly stored in the LCIO~\cite{LCIO} format instead of previous stdhep format in the DBD sample.
The generated samples are simulated with the ILD large detector model (ILD\_l5\_v02) using the dd4hep~\cite{dd4hep}.
ILD\_l5\_v02 is a hybrid calorimeter model composed of the silicon and scintillator tungsten electromagnetic calorimeter and analog scintillator and semi-digital RPC hadron calorimeter.
The simulated data (SIM) are reconstructed by the Marlin package~\cite{Marlin} with the silicon electromagnetic and analog scintillator hadron calorimeter option (ILD\_l5\_o1\_v02), creating REC and DST files.
REC files contain full reconstruction information while DST files contain only partial information suitable for physics studies.
The software version is ILCSoft-v02-02~\cite{ILCSoft}.
In the reconstruction process, low $p_T$ hadrons by $\gamma\gamma$ interactions and $e^{+}e^{-}$ pairs by bremsstrahlung backgrounds are overlaid to the simulated events.
According to an estimation, required storage space for the requested statistics is $4.6~\mathrm{PB}$ by SIM and $5.5~\mathrm{PB}$ by REC, which exceeds the available storage capacity of ILD.
Thus we decided not to keep all the SIM and REC files.
In this production, DST files are only all kept but only 10\% or at least 500 jobs of REC files are partially kept.
Similar rule applied to the SIM files.

\section{MC mass production workflow in ILD}
ILD MC production workflow is conducted by the ILDProd tool~\cite{ILDProd} utilizing ILCDirac as shown in Fig.~{\ref{Fig:workflow}}.
At first, a list of production request is prepared in an excel format.
Generator file splitter and production transformation submission scripts are created from the input excel file,
and a json format file is created for controlling subsequent production workflow by python,
then a series of production is launched by a start script by operator.
After launching production,
at the first step, generator samples are split into small chunk by GenSplit. It is processed at the KEKCC local system.
Split generator files are then uploaded to the KEK-DISK storage element (SE), registered into a DIRAC file catalog,
and meta information is attached.
ILCDirac production transformations are submitted for simulation and reconstruction soon after launching the GenSplit process.
In the MC-2020 production, following ILCDirac transformations are used in order to cope with the difference of keeping files:
SIM, REC partially-kept, DST only-kept.  REC and DST files are fully kept in the case of calibration samples.
Job submission and rescheduling in the case of job failure are managed by ILCDirac. Produced files are uploaded to KEK or DESY disk only SEs.
In order to reduce the heavy load by retrieving the background files for overlay,
$\gamma\gamma$ low $p_T$ hadron and $e^{+}e^{-}$ pair background samples are replicated to KEK, DESY, and CERN SEs.\\
~~The completion of transformation is decided by the operator when all input files without error are processed.
After terminating the production transformations, DST merge and log retrieve steps are initiated by operator.
Since DST file size is relatively small, merged-DST files are created by merging 50 to 100 DST files.
The DST merge is done by utilizing ILCDirac user job and rescheduling in the case of job failure is managed by the ILDProd tool.
All job output logs are also retrieved and archived on the SEs complying the ILD production convention.
Production status is monitored via a web-based monitor and jobs are added periodically with a cron task.
In MC-2020 production, an automation  of rescheduling of DST-merge jobs and retrying of GenSplit registration step
is enhanced by ILDProd, and a need of manual intervention of production process is reduced compared to the previous production.
All the production step statuses are recorded automatically on ELOG~\cite{ELOG} system and
a production summary page~\cite{ProdSummary} is updated automatically when the DST merge and log retrieve
steps are completed.
This workflow is repeated until all the input groups are processed.

\begin{figure}[htbp]
\begin{center}
 \includegraphics[width=0.5\textwidth]{figures/workflow.eps}
 \caption{ILD production workflow with ILCDirac managed by ILDProd tool.}
 \label{Fig:workflow}
 \end{center}
 \end{figure}

 \section{MC production progress and status}

In October 2020, MC-2020 production has been launched
and single particle, higgs, and small or middle cross section channels were processed by the end of November 2020
as ``1st-stage'' period shown in Fig.~\ref{Fig:ProdProgress}.
After this period, some failure jobs, and LCFIPlus~\cite{LCFIPlus} vertex finder training samples were processed as of  ``left-over'' period.
In these period, all the SIM files were kept. In the ``2nd-stage'' period,
SIM files were not kept any more and they were removed after the completion of merged-DST production.
In this period, more than 3.5 M jobs have been processed and process is still on-going.

 \begin{figure}[htbp]
 \begin{center}
 \includegraphics[width=0.5\textwidth]{figures/stuckJobType.eps}
 \caption{Accumulated number of production jobs for simulated (SIM) and reconstructed (REC) production.}
\label{Fig:ProdProgress}
\end{center}
 \end{figure}

 Figure~\ref{Fig:JobType} and \ref{Fig:JobStatus} show the number of processed jobs clustered by job types and final status.
 In the period of this figure, 3 to 4 k jobs run concurrently in average, about 1 k of them were REC jobs.
 At the early stage of the production, we faced high job failure rate by ``Stalling'' of REC jobs,
 which could be caused by the slow or timeout of downloading background overlay files.
 Even though we replicated BG overlay files to several SEs, we found that jobs tried to download BG overlay file only from KEK-DISK SE, resulting high network overload.
 Figure~\ref{Fig:JobSites} shows the job executed sites of this period, showing that jobs were mostly processed at European sites.
Therefore, we were forced to limit the number of concurrent running REC jobs up to 1.5 k in this period.

 \begin{figure}[htbp]
 \begin{minipage}{0.48\textwidth}
 \begin{center}
 \includegraphics[width=0.9\textwidth]{figures/jobType.eps}
 \caption{Executed production jobs by job types.}
 \label{Fig:JobType}
  \end{center}
 \end{minipage}
  \begin{minipage}{0.48\textwidth}
  \begin{center}
 \includegraphics[width=0.9\textwidth]{figures/jobStatus.eps}
 \caption{Executed production jobs by final statuses.}
 \label{Fig:JobStatus}
 \end{center}
 \end{minipage}
 \end{figure}

 The problem was notified to the ILCDirac administrator through JIRA issue tracking system.
 Thanks to the ILCDirac administrators effort, an implementation of SE selection mechanism was improved in the ILCDirac.
 Also we restricted the SE for output files to DESY-SRM in order to avoid the heavy load to KEK-DISK.
 After applying these treatment, high failure rate of the REC overlay process has been reduced and productions are much stabilized.

\begin{figure}[htbp]
\begin{center}
 \includegraphics[width=0.5\textwidth]{figures/jobSites.eps}
 \caption{Summary of the job executed sites.}
 \label{Fig:JobSites}
 \end{center}
 \end{figure}

 Figure~\ref{Fig:StorageSpace} shows the occupancy of the storage used in the ILD production.
 A capacity of the disk storage is 300 TB for KEK-DISK SE and 600 TB for DESY-SRM.
 SIM files, which were stored in disk storage in early period were moved to the KEK-SRM (TAPE)  SE
 and removed from disk storage in a timely manner.
 Unfortunately, due to a miss operation by a production team at the end of January 2021,
 SIM and REC files were removed, but merged-DST files were not affected.
The lost files included background files for overlay.  The production were restarted after
 reproducing background files. The 2nd stage $4f$ and $2f$ channels were produced
 in this period and it continued rather smoothly.

 \begin{figure}[htbp]
 \begin{center}
 \includegraphics[width=0.5\textwidth]{figures/storage.eps}
 \caption{Occupancy of storage space used in the ILD production; KEK-DISK, KEK-SRM, and DESY-SRM.}
 \label{Fig:StorageSpace}
 \end{center}
 \end{figure}

 At the middle of March 2021,  ILCDirac and its base DIRAC version have been updated on the ILCDirac server.
 For this update, production jobs have been drained before the downtime and production was paused.
 After resuming the ILCDirac servers, we applied small test production with updating ILCDirac client version
 and small issues have observed but quickly fixed.
 Meanwhile, we reproduced some forgotten samples; single particles, LCFI calibration, and higgs samples with semi-digital hadron calorimeter option (ILD\_l5\_o2\_v02).
 After confirming the production workflow with updated ILCDirac,
 production has been resumed and $2f$ large cross section channel production is launched.

 \section{Conclusion}
 We have conducted mass production of new 250 GeV MC samples fully utilizing the ILCDirac system since autumn 2020.
 According to the initial plan, more than 2 Giga events of about 240 processes are processed creating about 130 TB of DST files for high statistics physics studies.
 By the end of 2021 April, the production of all initial processes has been completed except $2f$ hadronic channels.
It is expected to take a couple of more months to be completed remaining $2f$ channels.
In order to speed up, the ILDProd parameters are being tuned for larger number of concurrent jobs with lower job failure rate.
Additional channels, such as $6f$ and $ee\gamma$, will be produced when generator samples are ready to use.

 \section*{Acknowledgments}

We would like to thank ILCDirac administrators, especially Andre Sailer,
 for the support and maintenance of ILCDirac system
 and ILD software group conveners for the continuous support and discussion regarding the production.
 We would also acknowledge all the computing sites who provide the resources and support for ILC Virtual Organization.

  \begin{footnotesize}