#  Restarted Nonnegativity Preserving Tensor Splitting Methods via Relaxed
Anderson Acceleration for Solving Multi-linear Systems

Dongdong Liu School of Mathematics and Statistics, Guangdong University of
Technology, Guangzhou, China. (ddliu@gdut.edu.cn). Xifu Liu Corresponding
author, School of Mathematical Sciences, Chongqing Normal University,
Chongqing<EMAIL_ADDRESS>
###### Abstract

Multilinear systems play an important role in scientific calculations of
practical problems. In this paper, we consider a tensor splitting method with
a relaxed Anderson acceleration for solving multilinear systems. The new
method preserves nonnegativity for every iterative step and improves the
existing ones. Furthermore, the convergence analysis of the proposed method is
given. The new algorithm performs effectively for numerical experiments.

Key words. Tensor splitting, multi-linear systems, Anderson acceleration,
strong $\mathcal{M}$-tensor MSC 15A48 15A69 65F10 65H10

##  1 Introduction

Let $\mathcal{A}=(a_{i{i_{2}}\cdots{i_{l}}})$ be an order $l$ dimension $n$
tensor and ${\bf{b}}\in\mathbb{R}^{n}$. The multi-linear system was first
given in [5]:

$\mathcal{A}{\bf{x}}^{l-1}={\bf{b}},$ (1.1)

where ${\bf{x}}$ is the unknown vector to be solved. The $i$-th entry of the
vector $\mathcal{A}{\bf{x}}^{l-1}$ (see [22, 23]) is defined as follows:

${(\mathcal{A}{{\bf{x}}^{l-1}})_{i}}\equiv\sum\limits_{{i_{2}},\cdots,{i_{l}}=1}^{n}{{a_{i{i_{2}}\cdots{i_{l}}}}{x_{{i_{2}}}}}\cdots{x_{{i_{l}}}},{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}i=1,2,...,n,$ (1.2)

where $x_{i}$ denote the $i$-th component of a vector ${\bf{x}}$. Since there
are many applications such as tensor complementarity problems(e.g., see [5,
32]), numerical partial differential equations(e.g., see [5, 14]) and Bellman
equation(e.g., see [1, 26]), the nonlinear equation (1.1) has attracted many
researchers’ attention, e.g., theoretical analysis of solutions(e.g.,see [5,
12, 15, 27, 28]), numerical algorithms such as splitting type methods(e.g.,see
[9, 13, 29]), the Gradient-based iteration(e.g., [2, 5, 6, 7, 14, 18, 19,
32]), dynamical approaches(e.g.,see [26, 30]) and so on. In [9], several
tensor splitting methods were given via setting
$\mathcal{A}=\mathcal{E}-\mathcal{F}$:

${{\bf{x}}_{k}}={(M{(\mathcal{E})^{-1}}\mathcal{F}{\bf{x}}_{k-1}^{l-1}+M{(\mathcal{E})^{-1}}{\bf{b}})^{[\frac{1}{{l-1}}]}},$
(1.3)

the iterative solution of which is a fixed point of the following iterative
function

$g_{\mathcal{E}}({{\bf{x}}})={(M{(\mathcal{E})^{-1}}\mathcal{F}{\bf{x}}^{l-1}+M{(\mathcal{E})^{-1}}{\bf{b}})^{[\frac{1}{{l-1}}]}},$
(1.4)

where $M{(\mathcal{E})^{-1}}\mathcal{F}$ was the iterative tensor of the
splitting method (1.3). Some preconditioned tensor splitting methods were
further studied(e.g., see [3, 10, 11, 34]). It is noted that tensor splitting
type methods have a fine numerical performance in terms of the running time
because of avoiding calculation of the Jacobian matrix(e.g., see [10]).
Therefore, an improved tensor splitting method combing with the efficient
acceleration technique is still an attractive topic.

Anderson acceleration is a powerful extrapolation technique to improve the
fixed-point iteration for solving the linear or nonlinear systems(e.g., see
[25, 31]). However, there is no work yet on solving the multilinear systems by
this technique. By this motivation, we propose a tensor splitting method with
the relaxed Anderson acceleration and the nonnegativity is preserved of the
new iteration in every iterative step for solving the multilinear systems with
the coefficient tensor being a strong $\mathcal{M}$-tensor. Furthermore, the
convergence analysis is given for preserving the nonnegativity.

The rest of this paper is given as follows. We introduce some definitions and
lemmas in Section 2. Besides, the tensor splitting method with the relaxed
Anderson acceleration is proposed in Section 3. Furthermore, the convergence
analysis is considered in Section 4. We report numerical experiments to
present the efficient performance of the proposed method in Section 5. The
final section is the conclusion.

##  2 Preliminaries

An order $l$ and dimension $n$ tensor $\mathcal{A}$ with $n^{l}$ elements:

$\mathcal{A}=(a_{i_{1}\cdots i_{l}}),a_{i_{1}\cdots
i_{l}}\in\mathbb{R},i_{j}\in\langle n\rangle,j=1,\cdots,l,$

where $\mathbb{R}$ denotes the real field, $\langle n\rangle=\\{1,\cdots,n\\}$
and $n$ is a positive integer. We use $\mathbb{R}^{[l,n]}$ to denote the set
of all order $l$ dimension $n$ real tensors. If $m=1$, $\mathbb{R}^{[1,n]}$ is
reduces to $\mathbb{R}^{n}$, which is commonly used to denote the set of all
dimension $n$ real vectors. More specifically,
$\mathbb{R}_{+}^{n}$($\mathbb{R}_{++}^{n}$) is indicted as the set of all
nonnegative(positive) vectors. Let ${{0}}$ denote a null real number, a null
vector, a null matrix or a null tensor. Let
$\mathcal{A},~{}\mathcal{B}\in\mathbb{R}^{[l,n]}$. The notation
$\mathcal{A}\geq\mathcal{B}(>\mathcal{B})$ refers that every element of
$\mathcal{A}$ is no less than (strictly greater than) the corresponding one of
$\mathcal{B}$. Furthermore, $||\cdot||$ denotes 2-norm of a vector or a
matrix.

Qi [22] and Lim [17] independently gave the definition of tensor eigenvalues
and eigenvectors, i.e., a complex number $\lambda$ and a nonzero complex
vector ${\bf{x}}$ is called as an eigenvalue and eigenvector of
$\mathcal{A}\in\mathbb{R}^{[l,n]}$, respectively, if the following equation
hold:

$\mathcal{A}{\bf{x}}^{l-1}=\lambda{\bf{x}}^{[l-1]},$

where ${\bf{x}}^{[l-1]}=(x_{1}^{l-1},...,x_{n}^{l-1})^{T}$. Besides,
$(\lambda,{\bf{x}})$ is called as an $H$-eigenpair if both $\lambda$ and
${\bf{x}}$ are real. Furthermore, the spectral radius of $\mathcal{A}$ is
signified by
$\rho(\mathcal{A})=\max\\{\left|\lambda\right||\lambda\in\sigma(\mathcal{A})\\}$
with $\sigma(\mathcal{A})$ being the set of all eigenvalues of $\mathcal{A}$.
We call $\mathcal{I}_{l}\in\mathbb{R}^{[l,n]}$ an unit tensor if its entries
are given by:

${\delta_{{i_{1}}\cdots{i_{l}}}}=\left\\{\begin{array}[]{l}1,{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{i_{1}}=\cdots={i_{l}},\\\
0,{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern 1.0pt}{\kern
1.0pt}{\kern 1.0pt}{\kern 1.0pt}\mbox{else.}\\\ \end{array}\right.$

If $\mathcal{A}\in\mathbb{R}^{[l,n]}$ is a $\mathcal{Z}$-tensor if its off-
diagonal entries are non-positive. Besides, $\mathcal{A}$ is an
$\mathcal{M}$-tensor(see the works of [4, 33]) if there exists a nonnegative
tensor $\mathcal{B}$ and a positive real number $\eta\geq\rho(\mathcal{B})$
such that $\mathcal{A}=\eta\mathcal{I}_{l}-\mathcal{B}.$ Furthermore, if
$\eta>\rho(\mathcal{B})$, $\mathcal{A}$ is a strong $\mathcal{M}$-tensor. Let
$A\in{\mathbb{R}^{[2,n]}}$ and $\mathcal{B}\in{\mathbb{R}^{[k,n]}}$. A product
[8] $\mathcal{C}=A\mathcal{B}\in\mathbb{R}^{[k,n]}$ is introduced by

${c_{j{i_{2}}\cdots{i_{k}}}}=\sum\limits_{{j_{2}}=1}^{{n}}{{a_{j{j_{2}}}}{b_{{j_{2}}{i_{2}}\cdots{i_{k}}}}}.$

The majorization matrix $M(\mathcal{A})$ of $\mathcal{A}$ is defined by
$M(\mathcal{A})_{ij}=a_{ij...j}$ with $i,j\in\langle n\rangle$. In [9], it is
pointed out that if $\mathcal{A}$ is a left-nonsingular tensor if and only if
$\mathcal{A}=M(\mathcal{A})\mathcal{I}_{m}$ and $M(\mathcal{A})$ is
nonsingular. If $\mathcal{A}$ is a strong $\mathcal{M}$-tensor, then
$M(\mathcal{A})$ is a nonsingular $\mathcal{M}$-matrix. Besides, the some
concepts on tensor splittings were listed as follows.

######  Definition 2.1.

[9] Let $\mathcal{A},{\kern 2.0pt}\mathcal{E},{\kern
2.0pt}\mathcal{F}\in\mathbb{R}^{[l,n]}$.

(1)

$\mathcal{A}=\mathcal{E}-\mathcal{F}$ is a splitting of $\mathcal{A}$ if and
only if $\mathcal{E}$ is a left-nonsingular.

(2)

$\mathcal{A}=\mathcal{E}-\mathcal{F}$ is a regular splitting of $\mathcal{A}$
if and only if $\mathcal{E}$ is left-nonsingular with $M(\mathcal{E})^{-1}\geq
0$, and $\mathcal{F}\geq 0$.

(3)

$\mathcal{A}=\mathcal{E}-\mathcal{F}$ is a weak regular splitting of
$\mathcal{A}$ if and only if $\mathcal{E}$ is left-nonsingular with
$M(\mathcal{E})^{-1}\geq 0$ and $M(\mathcal{E})^{-1}\mathcal{F}\geq 0$.

(4)

$\mathcal{A}=\mathcal{E}-\mathcal{F}$ is a convergent splitting if and only if
$\rho(M(\mathcal{E})^{{-1}}\mathcal{F})<1$.

Let $\mathcal{A}$ be a $\mathcal{Z}$-tensor. By the works of [9], it is noted
that $\mathcal{A}$ is a strong $\mathcal{M}$-tensor if and only if
$\mathcal{A}$ has a convergent (weak) regular splitting(or all (weak) regular
splittings of $\mathcal{A}$ are convergent). The following lemmas on
multilinear systems in the existing works will be used in Section 4.

######  Lemma 2.2.

[5] Let $\mathcal{A}$ be a strong $\mathcal{M}$-tensor and
${\bf{b}}\in\mathbb{R}^{n}_{++}$. Thus, the multi-linear system (1.1) has a
unique positive solution.

######  Lemma 2.3.

[9] Let $\mathcal{A}$ be a strong $\mathcal{M}$-tensor and
${\bf{b}}\in\mathbb{R}^{n}_{++}$. Thus, there always exists a
${\bf{x}}_{0}\in\mathbb{R}^{n}_{++}$ with
${\bf{0}}<\mathcal{A}{{\bf{x}}_{0}^{m-1}}\leq{\bf{b}}$.

##  3 Tensor splitting methods with Anderson acceleration

Let $f_{\mathcal{E}}({\bf{x}})=g_{\mathcal{E}}({\bf{x}})-{\bf{x}}.$ The tensor
splitting methods with the relaxed Anderson acceleration is proposed for
solving the multilinear systems (1.1) as follows.

Algorithm 1 Restarted tensor splitting methods with the relaxed Anderson
acceleration(RTSMRAA)

1: Require: Given a positive vector ${\bf{b}}$, a strong $\mathcal{M}$-tensor
$\mathcal{A}$ with its (weak) regular splitting
$\mathcal{A}=\mathcal{E}-\mathcal{F}$, maximum $k_{{\bf{max}}}$, a positive
initial vector ${\bf{z}}_{0}$ and $m\geq 1$.

2: Update ${\bf{z}}_{1}=g_{\mathcal{E}}({\bf{z}}_{0})$,
${\boldsymbol{\mu}}_{1}={\bf{z}}_{1}$, ${\bf{y}}_{1}={\bf{z}}_{1}$ and
$f_{0}={\bf{z}}_{1}-{\bf{z}}_{0}$

3: for k=1,2,…, $k_{{\bf{max}}}$ do

4: Set $m_{k}=\min\\{m,k\\}$

5: Update
$\begin{split}{\boldsymbol{\mu}}_{k+1}&=g_{\mathcal{E}}({\bf{z}}_{k}),\\\
f_{k}&={\boldsymbol{\mu}}_{k+1}-{\bf{z}}_{k}.\end{split}$ (3.1)

6: Set $Q_{k}=(f_{k-m_{k}},...,f_{k})$.

7: Determine
$\boldsymbol{\alpha}^{(k)}=(\alpha_{0}^{(k)},...,\alpha_{m_{k}}^{(k)})^{T}$
that solves
$\min_{\alpha=(\alpha_{0},...,\alpha_{m_{k}})^{T}}||Q_{k}\alpha||_{2}~{}\textit{s.t.}~{}\sum\limits_{i=0}^{m_{k}}\alpha_{i}=1.$
(3.2)

8: Update
${\bf{y}}_{k+1}=\sum_{i=0}^{m_{k}}\alpha_{i}^{(k)}g_{\mathcal{E}}(z_{k-m_{k}+i}).$

9: if ${\bf{y}}_{k+1}\geq 0$ and
$\sum\limits_{i=0}^{m_{k}}|\alpha_{i}^{(k)}|\leq\kappa_{\boldsymbol{\alpha}}$
then

10: Choose $\theta_{k}\in[0,1]$ to update
${\bf{z}}_{k+1}=\theta_{k}{\bf{y}}_{k+1}+(1-\theta_{k}){{\boldsymbol{\mu}}}_{k+1}$.

11: else

12: Update ${\bf{z}}_{k+1}=\boldsymbol{\mu}_{k+1}$.

13: end if

14: Update $k\leftarrow k+1$.

15: end for

16: return ${\bf{z}}_{k}$

######  Remark 3.1.

The optimization problem (3.2) can be solved by the following equivalent the
square problem[31]: Give
$\boldsymbol{\zeta}^{(k)}=(\zeta_{0}^{(k)},\zeta_{1}^{(k)},...,\zeta_{m_{k}-1}^{(k)})$
by solving

$\min_{\boldsymbol{\zeta}=(\zeta_{0}^{(k)},\zeta_{1}^{(k)},...,\zeta_{m_{k}-1}^{(k)})^{T}}||f_{k}-F_{k}\boldsymbol{\zeta}^{(k)}||_{2},$
(3.3)

where $\Delta F_{k}=(\Delta f_{k-m_{k}},...,\Delta f_{k-1})$ and $\Delta
f_{i}=f_{i+1}-f_{i}$. Furthermore, Step 8 in RTSMRAA is written as follows.

${\bf{y}}_{k+1}={\boldsymbol{\mu}}_{k+1}-G_{k}\boldsymbol{\zeta}^{(k)},$ (3.4)

where $G_{k}=(\Delta g_{k-m_{k}},...,\Delta g_{k-1})$ and $\Delta
g_{i}=g_{\mathcal{E}}({\bf{z}}_{i+1})-g_{\mathcal{E}}({\bf{z}}_{i}).$ The
solution $\boldsymbol{\zeta}$ of the least square problem (3.3) and the
solution $\alpha$ of the optimization problem 3.2 are related by
$\alpha_{0}=\zeta_{0},\alpha_{i}=\zeta_{i}-\zeta_{i-1}$ for $i=1,...,m_{k}-1$
and $\alpha_{m_{k}}=1-\zeta_{m_{k}-1}$.

##  4 Convergence analysis

In this section, we will consider the convergence analysis of Algorithm 3. Let
$\mathbb{U}_{\delta}\equiv\\{{\bf{x}}\big{|}||{\bf{x}}-{\bf{x}}^{*}||\leq\delta,\delta>0\\}$.
We write $a_{i_{1}i_{2}\cdots i_{l}}$ as $a_{i_{1}\omega}$, where
$\omega=i_{2}\cdots i_{l}$. Notice that for any
$\mathcal{A}=(a_{i_{1}i_{2}\cdots i_{l}})\in{\mathbb{R}}^{[l,n]}$, it is easy
to check that

${\mathcal{A}}{\bf{x}}^{l-1}=\overline{\mathcal{A}}{\bf{x}}^{l-1}$

where $\overline{\mathcal{A}}=(\hat{a}_{i_{1}i_{2}\cdots i_{m}})$ is given by:

${\overline{a}_{ii_{2}\cdots
i_{l}}}=\frac{1}{(l-1)!}\sum\limits_{\omega\in\mathbb{G}}{{a_{i\omega}}},$

and $\mathbb{G}$ denotes the set of all permutations in
$({i_{2}}\cdots{i_{l}})$. Then $\overline{\mathcal{A}}$ is symmetric on the
last $m-1$ indices ($\overline{\mathcal{A}}$ is called semi-symmetric[20]).
For the sake of clarity, we always denote the semi-symmetric of a tensor by
the notation $\overline{{\cdot}}$(e.g., $\overline{\mathcal{A}}$). Firstly, we
give some lemmas as follows.

######  Lemma 4.1.

Let $\mathcal{A}$ be a strong $\mathcal{M}$-tensor,
${\bf{b}}\in\mathbb{R}^{n}_{++}$ and $\mathcal{A}=\mathcal{E}-\mathcal{F}$ be
a (weak) regular splitting. Then the function $g_{\mathcal{E}}({\bf{x}})$ has
the following properties:

(1)

$g_{\mathcal{E}}({\bf{x}})\in\mathbb{R}^{n}_{++}$ for all
${\bf{x}}\in\mathbb{R}^{n}_{+}$.

(2)

$g_{\mathcal{E}}({\bf{x}})$ is monotone increasing in $\mathcal{R}^{n}_{+}$,
i.e., for any $0\leq{\bf{x}}\leq{\bf{y}}$, $g_{\mathcal{E}}({\bf{x}})\leq
g_{\mathcal{E}}({\bf{y}})$.

Proof: Since $\mathcal{A}=\mathcal{E}-\mathcal{F}$ is a (weak) regular
splitting, the second assertion is obvious. Because
${\bf{b}}\in\mathbb{R}^{n}_{++}$, it is obtained that
$M(\mathcal{E})^{-1}{\bf{b}}\in\mathbb{R}^{n}_{++}$. Thus, the first
affirmation is proved.

######  Proposition 4.2.

Under the assumption of Lemma 4.1, for any initial vector
${\bf{z}}_{0}\in\mathbb{R}^{n}_{++}$ with
${{0}}<\mathcal{A}{{\bf{x}}_{0}^{l-1}}\leq{\bf{b}}$, we get
${\bf{z}}_{k}\in\mathbb{R}^{n}_{++}$ generated by RTSMRAA if
$\theta_{k}\in[0,1)$ in Step 10 of RTSMAA for all $k\geq 1$.

Proof: Since $\mathcal{A}{\bf{z}}_{0}^{l-1}\leq{\bf{b}}$, it is true that
${\bf{z}}_{0}\leq{\bf{z}}_{1}$ under the assumption of Lemma 4.1. By Lemma
4.1, we know $\boldsymbol{\mu}_{k}\in\mathbb{R}^{n}_{++}$ for any $k\geq 1$
via the induction. Then the proof is completed.

######  Remark 4.3.

If $\theta_{k}=1$ in Step 10 of RTSMRAA for some $k$, it is note that TSM is
just improved by the restarted Anderson acceleration for the $k$-th iteration,
i.e., ${\bf{z}}_{k}={\bf{y}}_{k}$. In this case, we still set
${\bf{y}}_{k+1}\geq 0$ rather than ${\bf{y}}_{k+1}>0$ in Step 9 of RTSMRAA.
That’s because, by Lemma 4.1, we can still get a vector
${\bf{z}}_{k+1}\in\mathbb{R}^{n}_{++}$.

######  Lemma 4.4.

Let ${\bf{x}}^{*}$ be the solution of the multilinear systems (1.1) with
$\mathcal{A}$ being a strong $\mathcal{M}$-tensor and
${\bf{b}}\in\mathbb{R}^{n}_{++}$. For every (weak) regular splitting
$\mathcal{A}=\mathcal{E}-\mathcal{F}$, there exists a $\mathbb{U}_{\delta}$
such that the mapping
$g_{\mathcal{E}}({\bf{x}}):\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ has the
following properties for a sufficiently small $\delta$:

(i)

For any ${\bf{x}}\in\mathbb{U}_{\delta}$,
$g_{\mathcal{E}}({\bf{x}})\in\mathbb{U}_{\delta}$.

(ii)

For any ${\bf{x}}$, ${\bf{y}}\in\mathbb{U}_{\delta}$, there exists a real
number $c\in[0,1)$ such that

$||g_{\mathcal{E}}({\bf{x}})-g_{\mathcal{E}}({\bf{y}})||\leq
c||{\bf{x}}-{\bf{y}}||.$

(iii)

For any ${\bf{x}}\in\mathbb{U}_{\delta}$,

$\nabla
g_{\mathcal{E}}({\bf{x}})=\big{(}M(\mathcal{E})g_{\mathcal{E}}({\bf{x}})^{[m-2]}\big{)}^{-1}\overline{F}{\bf{x}}^{m-2},~{}\nabla
f_{\mathcal{E}}({\bf{x}})=\nabla g_{\mathcal{E}}({\bf{x}})-I.$

(iv)

$\nabla g_{\mathcal{E}}({\bf{x}}^{*})$ is Lipschitz continuously, i.e., for
any ${\bf{x}}$, ${\bf{y}}\in\mathbb{U}_{\delta}$, there exists a real number
$\gamma>0$ such that

$||\nabla f_{\mathcal{E}}({\bf{x}})-\nabla
f_{\mathcal{E}}({\bf{y}})||\leq\gamma||{\bf{x}}-{\bf{y}}||.$

Proof: Since $\mathcal{A}$ is a strong $\mathcal{M}$-tensor and ${\bf{b}}$ is
positive, by Lemma 2.2, ${\bf{x}^{*}}$ is the unique positive solution. By the
proof of Theorem 5.7 in [9], we have that
$M(\mathcal{E})g_{\mathcal{E}}({\bf{x}^{*}})^{[m-2]}=M(\mathcal{E})({\bf{x}^{*}})^{[m-2]}$
is nonsingular and $\nabla
g_{\mathcal{E}}({\bf{x}})=\big{(}M(\mathcal{E})g_{\mathcal{E}}({\bf{x}})^{[m-2]}\big{)}^{-1}\overline{F}{\bf{x}}^{m-2}$
satisfying $|\rho(\nabla g_{\mathcal{E}}({\bf{x}}^{*}))|\leq c<1$. Because
$g_{\mathcal{E}}({\bf{x}})$ and $\nabla g_{\mathcal{E}}({\bf{x}})$ are
continuous, there exists $\mathbb{U}_{\delta}$ such that $|\rho\big{(}\nabla
g_{\mathcal{E}}({\bf{x}})\big{)}|\leq c<1$ and
$M(\mathcal{E})g_{\mathcal{E}}({\bf{x}})^{[m-2]}$ is nonsingular for any
${\bf{x}}\in\mathbb{U}_{\delta}$. Thus, it is easy to check that the
affirmation (ii) holds by differential mean value theorem and the term (iii)
is given. Besides, for any ${\bf{x}}\in\mathbb{U}_{\delta}$, we have

$\displaystyle||g_{\mathcal{E}}({\bf{x}})-{\bf{x}}^{*}||=||g_{\mathcal{E}}({\bf{x}})-g_{\mathcal{E}}({\bf{x}}^{*})||\leq
c||{\bf{x}}-{\bf{x}}^{*}||<||{\bf{x}}-{\bf{x}}^{*}||\leq\delta.$

Then the desired result (i) is obtained. Furthermore, by the result 3.3.5 of
[21] and the second G-derivative $\nabla^{2}f_{\mathcal{E}}({\bf{x}})$ of
$\nabla f_{\mathcal{E}}({\bf{x}})$ at ${\bf{x}}$ is continuous, we have

$||\nabla f_{\mathcal{E}}({\bf{x}})-\nabla
f_{\mathcal{E}}({\bf{y}})||\leq\max\limits_{{\bf{x},{\bf{y}}}\in\mathbb{R}^{n}_{+}}\max\limits_{t\in[0,1]}\big{(}||\nabla^{2}f_{\mathcal{E}}({\bf{x}}+t({\bf{y}}-{\bf{x}}))||\big{)}||{\bf{x}}-{\bf{y}}||.$

Since $\mathbb{R}^{n}_{+}$ is a closed set, the term (iv) is shown by
$\gamma=\max\limits_{{\bf{x},{\bf{y}}}\in\mathbb{R}^{n}_{+}}\max\limits_{t\in[0,1]}\big{(}||\nabla^{2}f_{\mathcal{E}}({\bf{x}}+t{\bf{y}}-{\bf{x}})||\big{)}.$

######  Proposition 4.5.

Under the assumption of Lemma 4.4, we have
$\big{(}\mathbb{I}\setminus\\{{\bf{x}}^{*}\\}\big{)}\cap\big{(}\mathbb{U}_{\delta}\setminus\\{{\bf{x}}^{*}\\})\neq\emptyset$,
where
$\mathbb{I}=\big{\\{}{\bf{x}}\big{|}{\bf{x}}>0,0<\mathcal{A}{\bf{x}}^{l-1}\leq
b\big{\\}}$ and $\mathbb{U}_{\delta}$ is defined in Lemma 4.4.

Proof: By Lemma 2.3, $\mathbb{I}\setminus\\{{\bf{x}}^{*}\\}\neq\emptyset$. Let
$p({\bf{x}})=\mathcal{A}{\bf{x}}^{l-1}$ which is a continuous function from
$\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$. Thus, there exists a real number
$\delta^{{}^{\prime}}>0$ such that,
$0<\mathcal{A}{\widehat{{\bf{x}}}}^{l-1}\leq b$ for any
$\widehat{{\bf{x}}}\in\mathbb{U}_{\delta^{{}^{\prime}}}\setminus\\{{\bf{x}}^{*}\\}$.
Let
$\delta^{{}^{\prime\prime}}=\min\big{\\{}\delta^{{}^{\prime}},\delta\big{\\}}$.
It is easy to check that
$\widehat{{\bf{x}}}\in\mathbb{U}_{\delta^{{}^{\prime\prime}}}\subseteq\mathbb{U}_{\delta}$.
Therefore, the proof is completed.

######  Lemma 4.6.

Under the assumption of Lemma 4.4, there exists
$\mathbb{U}_{\hat{\delta}}\subseteq\mathbb{U}_{\delta}$ with sufficiently
small $\hat{\delta}$ such that for any ${\bf{x}}\in\mathbb{R}^{n}_{++}$,

$||f_{\mathcal{E}}({\bf{x}})-\nabla
f_{\mathcal{E}}({\bf{x}}^{*})({\bf{x}}-{\bf{x}}^{*})||\leq\frac{\gamma}{2}||{\bf{x}}-{\bf{x}}^{*}||^{2}$
(4.1)

and

$(1-c)||{\bf{x}}-{\bf{x}}^{*}||\leq||f_{\mathcal{E}}({\bf{x}})||\leq(1+c)||{\bf{x}}-{\bf{x}}^{*}||.$
(4.2)

Proof: The formula (4.1) is a direct result by 3.2.12 of [21] and (iii) of
Lemma 4.4. The second formula can be show by (ii) of Lemma 4.4 and

$\displaystyle||f_{\mathcal{E}}({\bf{x}})||=||f_{\mathcal{E}}({\bf{x}})-f_{\mathcal{E}}({\bf{x}}^{*})||=||g_{\mathcal{E}}({\bf{x}})-g_{\mathcal{E}}({\bf{x}}^{*})+({\bf{x}}-{\bf{x}}^{*})||.$

Let $\boldsymbol{\epsilon}_{k}={\bf{z}}_{k}-{\bf{z}}^{*}$,
$\hat{\boldsymbol{\epsilon}}_{k}={\bf{y}}_{k}-{\bf{z}}^{*}$and
${\bf{z}}^{*}={\bf{x}}^{*}.$ Then we give the convergence analysis of the
proposed algorithm as follows.

######  Theorem 4.7.

Under the assumption of Lemma 4.4, for any initial vector
${\bf{z}}_{0}\in\mathbb{I}\cap\mathbb{U}_{\delta}$, there exists
$\mathbb{U}_{\hat{\delta}}\subseteq\mathbb{I}\cap\mathbb{U}_{\delta}$ with a
small enough $\hat{\delta}$ such that
${\bf{y}}_{k},~{}{\bf{z}}_{k}\in\mathbb{U}_{\hat{\delta}}$($k\geq 1$), and

$||f_{\mathcal{E}}({\bf{z}}_{k})||\leq\tau\overline{c}^{k}||f_{\mathcal{E}}({\bf{z}}_{0})||,$
(4.3)

where ${\bf{y}}_{k},~{}{\bf{z}}_{k}$ is generated by RTSMRAA,
$\forall\overline{c}\in(c,1)$ and $\tau=\frac{1+c}{1-c}$.

Proof: We can always choose the ${\bf{z}}_{0}\in\mathbb{U}_{\hat{\delta}}$
with a sufficiently small $\hat{\delta}$ such that
$\hat{\delta}<2(1-c)/\gamma$,

$\frac{\kappa_{\boldsymbol{\alpha}}\tau(2c+\gamma\hat{\delta})}{2(1-c)}\overline{c}^{-m}||f_{\mathcal{E}}({\bf{z}}_{0})||\leq\frac{\kappa_{\boldsymbol{\alpha}}\tau(2c+\gamma\hat{\delta})(1+c)}{2(1-c)}\overline{c}^{-m}||{\boldsymbol{\epsilon}}_{0}||\leq\hat{\delta}$
(4.4)

and

$\frac{\frac{c}{\overline{c}}+\frac{\kappa_{\boldsymbol{\alpha}}\gamma\hat{\delta}}{2(1-c)}\overline{c}^{-m-1}}{1-\frac{\hat{\delta}\gamma}{2(1-c)}}\leq
1.$ (4.5)

For any initial vector
${\bf{z}}_{0}\in\mathbb{U}_{\hat{\delta}}\subseteq\mathbb{I}\cap\mathbb{U}_{\delta}$
with a sufficiently small $\hat{\delta}$ such that (4.1) and (4.2) hold, by
Lemma 4.4 and the proof of Proposition 4.5, we know that there exists
$\mathbb{U}_{\hat{\delta}}\subseteq\mathbb{I}\cap\mathbb{U}_{\delta}$ such
that
${\bf{z}}_{1},~{}{\bf{y}}_{1},~{}{\boldsymbol{\mu}}_{1},~{}{\boldsymbol{\mu}_{2}}\in\mathbb{U}_{\hat{\delta}}$
and the inequality (4.3) holds. Assume that
${\bf{y}}_{k},~{}{\bf{z}}_{k}\in\mathbb{U}_{\hat{\delta}}$ and
$||f_{\mathcal{E}}({\bf{z}}_{k})||\leq\tau\overline{c}^{k}||f_{\mathcal{E}}({\bf{z}}_{0})||$
for any $0\leq k\leq N$. Now we will prove
${\bf{z}}_{N+1}\in\mathbb{U}_{\hat{\delta}}$ and the inequality
$||f_{\mathcal{E}}({\bf{z}}_{N+1})||\leq\tau\overline{c}^{N+1}||f_{\mathcal{E}}({\bf{z}}_{0})||$
by the induction. By Lemma 4.6, we have

$f_{\mathcal{E}}({\bf{z}}_{k})=\nabla
f_{\mathcal{E}}({\bf{z}}^{*}){{\boldsymbol{\epsilon}}_{k}}+{\varOmega}_{k},$
(4.6)

and
$||{\varOmega_{k}}||\leq\frac{\gamma}{2}||{\boldsymbol{\epsilon}_{k}}||^{2}$.
Therefore, it is easy to check that

$g_{\mathcal{E}}({\bf{z}}_{k})={\bf{z}}^{*}+\nabla
g_{\mathcal{E}}({\bf{z}}^{*})\boldsymbol{\epsilon}_{k}+\varOmega_{k}.$ (4.7)

Firstly, we will prove (4.3) and ${\bf{y}}_{N+1}\in\mathbb{U}_{\hat{\delta}}$
from the following two cases.

Case I: If ${\bf{y}}_{N+1}\leq 0$ or
$\sum\limits_{i=0}^{m_{k}}|\alpha_{i}^{(k)}|>\kappa_{\boldsymbol{\alpha}}$, by
Lemma 4.4, Lemma 4.6 and the hypothesis, we imply
${\bf{z}}_{N+1}=\boldsymbol{\mu}_{N+1}\in\mathbb{U}_{\hat{\delta}}$ and

$\displaystyle||f_{\mathcal{E}}({\bf{z}}_{N+1})||$
$\displaystyle\leq(1+c)||{\boldsymbol{\epsilon}}_{N+1}||$
$\displaystyle\leq(1+c)||g_{\mathcal{E}}({\bf{z}}_{N})-g_{\mathcal{E}}({\bf{z}}^{*})||$
$\displaystyle\leq c(1+c)||\boldsymbol{\epsilon}_{N}||$
$\displaystyle\leq\frac{c(1+c)}{1-c}||f_{\mathcal{E}}({\bf{z}}_{N})||$
$\displaystyle\leq\frac{1+c}{1-c}\overline{c}^{N+1}||f_{\mathcal{E}}({\bf{z}}_{0})||=\tau\overline{c}^{N+1}||f_{\mathcal{E}}({\bf{z}}_{0})||.$

This gives the inequality (4.3).

Case II: If ${\bf{y}}_{N+1}\geq 0$ and
$\sum\limits_{i=0}^{m_{k}}|\alpha_{i}^{(k)}|\leq\kappa_{\boldsymbol{\alpha}}$,
note that the fact $\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}=1$ combing with
(4.7) gives

${\bf{y}}_{N+1}={\bf{z}}^{*}+\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\nabla
g_{\mathcal{E}}({\bf{z}}^{*})\boldsymbol{\epsilon}_{N-m_{N}+i}+\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\varOmega_{N-m_{N}+i}.$
(4.8)

By Lemma 4.6, we can obtain

$||\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\varOmega_{N-m_{N}+i}||\leq\sum\limits_{i=0}^{m_{N}}|\alpha_{i}^{(N)}|||\varOmega_{N-m_{N}+i}||\leq\gamma\sum\limits_{i=0}^{m_{N}}|\alpha_{i}^{(N)}|||\boldsymbol{\epsilon}_{N-m_{N}+i}||^{2}/2.$
(4.9)

Since $m_{N}=\min\\{m,N\\}$, we get $N-m_{N}+i\geq N-m$. Thus, by Lemma 4.6
and the assumption, it is not difficult to get that

$\displaystyle||\boldsymbol{\epsilon}_{N-m_{N}+i}||\leq\frac{1}{1-c}||f_{\mathcal{E}}({\bf{z}}_{N-m_{N}+i})||\leq\frac{\tau}{1-c}\overline{c}^{N-m_{N}+i}||f_{\mathcal{E}}({\bf{z}}_{0})||\leq\frac{\tau}{1-c}\overline{c}^{-m}||f_{\mathcal{E}}({\bf{z}}_{0})||$
(4.10)

and

$\displaystyle\begin{split}||\boldsymbol{\epsilon}_{N-m_{N}+i}||^{2}&\leq\frac{1}{1-c}||\boldsymbol{\epsilon}_{N-m_{N}+i}||||f_{\mathcal{E}}({\bf{z}}_{N-m_{N}+i})||\\\
&\leq\frac{\hat{\delta}\tau}{1-c}\overline{c}^{N-m_{N}+i}||f_{\mathcal{E}}({\bf{z}}_{0})||\\\
&\leq\frac{\hat{\delta}\tau}{1-c}\overline{c}^{N-m}||f_{\mathcal{E}}({\bf{z}}_{0})||,\end{split}$
(4.11)

together with
$\sum\limits_{i=0}^{m_{N}}|\alpha_{i}^{(N)}|\leq\kappa_{\boldsymbol{\alpha}}$,
the inequalities (4.9) and (4.10), we imply that

$||\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\varOmega_{N-m_{N}+i}||\leq\frac{\kappa_{\boldsymbol{\alpha}}\gamma\hat{\delta}\tau}{2(1-c)}\overline{c}^{N-m}||f_{\mathcal{E}}({\bf{z}}_{0})||\leq\frac{\kappa_{\boldsymbol{\alpha}}\gamma\hat{\delta}\tau}{2(1-c)}\overline{c}^{-m}||f_{\mathcal{E}}({\bf{z}}_{0})||$
(4.12)

and

$||\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\nabla
g_{\mathcal{E}}({\bf{z}}^{*})\boldsymbol{\epsilon}_{N-m_{N}+i}||\leq\frac{\kappa_{\boldsymbol{\alpha}}\tau
c}{1-c}\overline{c}^{-m}||f_{\mathcal{E}}({\bf{z}}_{0})||.$ (4.13)

From (4.8), we obtain

${\hat{\boldsymbol{\epsilon}}_{N+1}}=\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\nabla
g_{\mathcal{E}}({\bf{z}}^{*})\boldsymbol{\epsilon}_{N-m_{N}+i}+\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\varOmega_{N-m_{N}+i},$

combing with (4.12), (4.13) and the condition (4.4) yields

$\displaystyle||\hat{\boldsymbol{\epsilon}}_{N+1}||\leq\frac{\kappa_{\boldsymbol{\alpha}}\tau(2c+\gamma\hat{\delta})}{2(1-c)}\overline{c}^{-m}||f_{\mathcal{E}}({\bf{z}}_{0})||\leq\hat{\delta},$
(4.14)

and hence ${\bf{y}}_{N+1}\in\mathbb{U}_{\hat{\delta}}$. Because
${\boldsymbol{\mu}_{N+1}}\in\mathbb{U}_{\hat{\delta}}$ by Lemma 4.4, it is
found that ${\bf{z}}_{N+1}\in\mathbb{U}_{\hat{\delta}}.$ Therefore, by Lemma
4.6, we can give the Taylor’s formula of $f_{\mathcal{E}}({\bf{z}})$ at the
true solution ${\bf{z}}^{*}$ as follows.

$f_{\mathcal{E}}({\bf{z}}_{N+1})=(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I){\boldsymbol{\epsilon}_{N+1}}+\varOmega_{N+1},$
(4.15)

with
$||\varOmega_{N+1}||\leq\frac{\gamma}{2}||\boldsymbol{\epsilon}_{N+1}||^{2}$.
This combing with (4.8) and the fact $\nabla
g_{\mathcal{E}}({\bf{z}}^{*})(\nabla g_{\mathcal{E}}({\bf{z}}^{*})-I)=(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I)\nabla g_{\mathcal{E}}({\bf{z}}^{*})$ yields

$\displaystyle\begin{split}f_{\mathcal{E}}({\bf{z}}_{N+1})=&(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I)\big{(}\theta_{k}\hat{\boldsymbol{\epsilon}}_{N+1}+(1-\theta_{k})(g_{\mathcal{E}}({\bf{z}}_{N})-{\bf{z}}^{*})\big{)}+\varOmega_{N+1}\\\
=&\theta_{k}\nabla
g_{\mathcal{E}}({\bf{z}}^{*})\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I){\boldsymbol{\epsilon}}_{N-m_{N}+i}+\theta_{k}(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I)\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\varOmega_{N-m_{N}+i}\\\
&+(1-\theta_{k})(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I)(g_{\mathcal{E}}({\bf{z}}_{N})-{\bf{z}}^{*})+\varOmega_{N+1}.\end{split}$
(4.16)

By (4.9)-(4.12), (4.15) the inductive assumption and the fact
$||\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}f_{\mathcal{E}}({\bf{z}}_{N-m_{N}+i})||\leq||f_{\mathcal{E}}({\bf{z}}_{N})||$,
it is obtained that

$\displaystyle\begin{split}&||\nabla
g_{\mathcal{E}}({\bf{z}}^{*})\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I){\boldsymbol{\epsilon}}_{N-m_{N}+i}+(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I)\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\varOmega_{N-m_{N}+i}||\\\
=&||\nabla
g_{\mathcal{E}}({\bf{z}}^{*})\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\big{(}f_{\mathcal{E}}({\bf{z}}_{N-m_{N}+i})-\varOmega_{N-m_{N}+i}\big{)}+(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I)\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\varOmega_{N-m_{N}+i}||\\\
=&||\nabla
g_{\mathcal{E}}({\bf{z}}^{*})\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}f_{\mathcal{E}}({\bf{z}}_{N-m_{N}+i})-\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\varOmega_{N-m_{N}+i}||\\\
\leq&c||f_{\mathcal{E}}({\bf{z}}_{N})||+||\sum\limits_{i=0}^{m_{N}}\alpha_{i}^{(N)}\varOmega_{N-m_{N}+i}||\\\
\leq&(\frac{c}{\overline{c}}+\frac{\kappa_{\boldsymbol{\alpha}}\gamma\hat{\delta}}{2(1-c)}\overline{c}^{-m-1})\tau\overline{c}^{N+1}||f_{\mathcal{E}}({\bf{z}}_{0})||\end{split}$
(4.17)

and

$\displaystyle\begin{split}||(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I)(g_{\mathcal{E}}({\bf{z}}_{N})-{\bf{z}}^{*})||=&||(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I)f_{\mathcal{E}}({\bf{z}}_{N})+(\nabla
g_{\mathcal{E}}({\bf{z}}^{*})-I)\boldsymbol{\epsilon}_{N}||\\\ \leq&||\nabla
g_{\mathcal{E}}({\bf{z}}^{*})f_{\mathcal{E}}({\bf{z}}_{N})||+||\varOmega_{N}||\\\
\leq&(\frac{c}{\overline{c}}+\frac{\gamma\hat{\delta}}{2(1-c)}\overline{c}^{-m-1})\tau\overline{c}^{N+1}||f_{\mathcal{E}}({\bf{z}}_{0})||.\end{split}$
(4.18)

Note that
$||\varOmega_{N+1}||\leq\frac{\gamma\hat{\delta}}{2(1-c)}||f_{\mathcal{E}}({\bf{z}}_{N+1})||$.
By (4.16)-(4.18) and the hypothesis (4.5) and the fact
$\kappa_{\boldsymbol{\alpha}}\geq 1$, it is not difficulty to get that

$\displaystyle\begin{split}||f_{\mathcal{E}}({\bf{z}}_{N+1})||&\leq\theta_{k}\frac{\frac{c}{\overline{c}}+\frac{\kappa_{\boldsymbol{\alpha}}\gamma\hat{\delta}}{2(1-c)}\overline{c}^{-m-1}}{1-\frac{\gamma\hat{\delta}}{2(1-c)}}\tau\overline{c}^{N+1}||f_{\mathcal{E}}({\bf{z}}_{0})||+(1-\theta_{k})\frac{\frac{c}{\overline{c}}+\frac{\gamma\hat{\delta}}{2(1-c)}\overline{c}^{-m-1}}{1-\frac{\gamma\hat{\delta}}{2(1-c)}}\tau\overline{c}^{N+1}||f_{\mathcal{E}}({\bf{z}}_{0})||\\\
&\leq\tau\overline{c}^{N+1}||f_{\mathcal{E}}({\bf{z}}_{0})||.\end{split}$
(4.19)

The proof is completed.

From Lemma 4.6 and Theorem 4.7, we can get the following corollary
immediately.

######  Corollary 4.8.

Under the assumption of Lemma 4.4, for any initial vector
${\bf{z}}_{0}\in\mathbb{I}\cap\mathbb{U}_{\delta}$, there exists
$\mathbb{U}_{\hat{\delta}}\subseteq\mathbb{I}\cap\mathbb{U}_{\delta}$ with
small enough $\hat{\delta}$ such that

$||\boldsymbol{\epsilon}_{k}||\leq\tau^{2}\overline{c}^{k}||\boldsymbol{\epsilon}_{0}||,$
(4.20)

where $\tau$ and $\overline{c}$ is given in Theorem 4.7.

##  5 Numerical examples

All the numerical experiments will be done in MATLAB R2018a with the
configuration: Intel(R) Core(TM) i9-11900K CPU 3.50GHz and 64.00G RAM. For all
numerical examples, we use ’CPU(s)’, ’IT’ and ’RES’ denote the running time,
the iterative step and the norm of a residual vector(i.e.,
$\textrm{RES}=||f_{\mathcal{E}}({\bf{z}}_{k})||$) of the tested algorithms,
respectively. We set the maximum iterative number as $1000$,
$\kappa_{\boldsymbol{\alpha}}=1000$ and the stopping criteria as
$||f_{\mathcal{E}}({\bf{z}}_{k})||<\varepsilon$ with $\varepsilon=10^{-11}$.
Let $\mathcal{A}=\mathcal{E}-\mathcal{F}$, $\mathcal{D}=D\mathcal{I}_{m}$ and
$\mathcal{L}=L\mathcal{I}_{m}$, where $D$ and $-L$ are the diagonal part and
the strictly lower triangle part of $M(\mathcal{A})$, respectively. By the
proposed acceleration technique, we give several illustrations on the tested
splitting methods listed as in Table 1.

Table 1: The proposed splitting methods with the related Anderson acceleration

The proposed methods | $\mathcal{E}$ | Abbreviations
---|---|---
Jacobi-Anderson type method | $\mathcal{D}$ | Jacobi-RAA(m)
G-S-Anderson type method | $\mathcal{D}-\mathcal{L}$ | G-S-RAA(m)
SOR-Anderson type method |  $\frac{1}{\omega}(\mathcal{D}-\omega\mathcal{L}$) | SOR-RAA(m)

######  Remark 5.1.

If $m=0$, it is noted that the numerical results are given by TSM.

We select the following two common examples in the existing literature(see [5,
10]). In numerical experiments, let $\theta=\theta_{k}$ for all $k$. We set
the right-hand side ${\bf{b}}={\bf{e}}$ and the initial value
${\bf{z}}_{0}={\bf{b}}$, where ${\bf{e}}=({1,1,...,1})^{T}$.

######  Example 5.2.

Taking order-$3$ dimension $n$ nonsingular $\mathcal{M}$-tensors
$\mathcal{A}=s\mathcal{I}-\mathcal{B}$, where $\mathcal{B}$ is generated
randomly by MATLAB, and
$s=(1+\varepsilon)\mathop{\max}\limits_{i=1,2,...,n}{(B{{\bf{e}}^{2}})_{i}},\varepsilon>0$,
and $\varepsilon=1$.

######  Example 5.3.

Let $\mathcal{A}\in\mathbb{R}^{[3,n]}$. Then entries $a_{111}=a_{nnn}=4.$ For
$i\in\langle n\rangle$, $a_{iii}=4$, $a_{i11}=a_{ii-1i}=a_{iii+1}=-1/3$.

It is easy to check that $\mathcal{A}$ in Examples 5.2 and 5.3 is a strong
$\mathcal{M}$-tensor and the corresponding multilinear system (1.1) has a
unique solution. For all numerical experiments, we set $n=100$. Firstly,
taking $\theta=0.7$, we draw Figures 1(a)-2(b) to show the running time with
the different $m$ in Step 4 of RTSMRAA for Examples 5.2-5.3. Besides, we also
report numerical results in Table 2 and the optimal parameters are listed in
Table 3, where ’-’means that it is unavailable. It is seen that if $m=1$ and
$m=3$, the SOR-type methods perform best for Examples 5.2-5.3, respectively.
The proposed algorithm is valid for a just small $m$ and increasing the value
of $m$ does not improve the algorithm further any better. This is a good
numerical behavior because a small value of $m$ means that there is not much
computation involved for solving the least square problem (3.3). Furthermore,
we give the variation of the running time with the different $\theta$ for the
tested methods when $m=3$. We search the optimal parameter $\theta$ from $0.1$
to $1$ in the interval of $0.1$ and the parameter $\omega$ from $0.1$ to $2$
in the interval of $0.1$ for the SOR type methods. The numerical results are
shown in Figures 2(a)-2(b), which show that the SOR-RAA($m$) method performs
best if the parameters $\theta$ and $\omega$ are chosen appropriately.

(a) Example 5.2

(b) Example 5.3

Figure 1: The running time of the tested algorithms with the variation of $m$.

(a) Example 5.2

(b) Example 5.3

Figure 2: The running time of the tested algorithms with the variation of
$\theta$. Table 2: The running time of the proposed algorithms with the
variation of $m$.

|  |  | Jacobain-RAA($m$) |  |  | GS-RAA($m$) |  |  | SOR-RAA($m$) |
---|---|---|---|---|---|---|---|---|---|---
| $m$ | CPU(s) | IT | RES | CPU(s) | IT | RES | CPU(s) | IT | RES
Example 5.2 | 0 | 0.082037 | 47 | 2.42E-12 | 0.080918 | 47 | 2.29E-12 | 0.053604 | 30 | 4.68E-12
| 1 | 0.018290 | 9 | 5.11E-14 | 0.019070 | 9 | 4.96E-14 | 0.019348 | 9 | 4.96E-14
| 2 | 0.019090 | 9 | 4.26E-15 | 0.019435 | 9 | 2.70E-15 | 0.019040 | 9 | 2.70E-15
| 3 | 0.023908 | 12 | 3.43E-12 | 0.022639 | 11 | 5.44E-15 | 0.023014 | 11 | 5.44E-15
| 4 | 0.024765 | 12 | 2.52E-12 | 0.024244 | 12 | 4.14E-12 | 0.024886 | 12 | 4.14E-12
Example 5.3 | 0 | 0.029323 | 16 | 3.77E-12 | 0.029483 | 16 | 2.43E-12 | 0.022033 | 12 | 2.15E-12
| 1 | 0.022518 | 11 | 3.74E-12 | 0.021555 | 11 | 4.32E-12 | 0.021878 | 11 | 5.85E-13
| 2 | 0.021770 | 11 | 7.41E-13 | 0.020660 | 10 | 4.58E-12 | 0.020009 | 10 | 2.23E-12
| 3 | 0.021846 | 11 | 4.40E-13 | 0.020020 | 10 | 1.83E-12 | 0.01971 | 10 | 1.31E-12
| 4 | 0.020400 | 10 | 3.10E-12 | 0.020933 | 10 | 1.31E-12 | 0.019767 | 10 | 4.73E-12

Table 3: The optimal parameters

| $m$ | Jacobain-RAA($m$) | GS-RAA($m$) | SOR-RAA($m$)
---|---|---|---|---
|  | $\theta$ | $\theta$ | ($\theta,~{}\omega$)
Example 5.2 | 0 | - | - | (-, 1.4)
| 1 | 1 | 1 | (1, 1)
| 2 | 1 | 1 | (1, 1)
| 3 | 0.9 | 1 | (1, 1)
| 4 | 1 | 0.9 | (0.9, 1)
Example 5.3 | 0 | - | - | (-, 1.1)
| 1 | 0.6 | 0.7 | (0.9, 1.1)
| 2 | 0.8 | 0.8 | (0.9, 1.1)
| 3 | 0.9 | 0.9 | (0.6, 1.1)
| 4 | 0.9 | 0.9 | (0.4, 1.1)

##  6 Conclusion

The main contribution in this paper is aimed at proposing the restarted
nonnegativity preserving tensor splitting methods by the relaxed Anderson
acceleration for solving the multilinear systems with a strong
$\mathcal{M}$-tensor and the local convergence analysis is also given. From
the numerical experiments, it is seen that the proposed methods with a
suitable parameter perform better than the tensor splitting methods without
acceleration. Besides, the increasing value of $m$ does not improve the given
algorithm better, and the proposed algorithm is valid just for a small value
of $m.$ However, how to find an optimal parameter $\theta_{k}$ is still an
open problem, which is one of our future works.
Declarations
Funding The first author was supported in part by National Natural Science
Foundation of China (No. 12101136), the Project of Science and Technology of
Guangzhou (No.202102020273), the Guangdong Basic and Applied Basic Research
Foundations (Nos. 2020A1515110967, 2021SFKC030), the Opening Project of
Guangdong Province Key Laboratory of Computational Science at the Sun Yat sen
University(No. 2021004), the Open Project of Key Laboratory, School of
Mathematical Sciences, Chongqing Normal University (No. CSSXKFKTQ202002). The
second author was funded by the Science and Technology Research Program of
Chongqing Municipal Education Commission (Grant No. KJQN202100505), the
Natural Science Foundation Project of Chongqing (Grant No. cstc2021jcyj-
msxmX0195), and the Program of Chongqing Innovation Research Group Project in
University (Grant No. CXQT19018).
Conflict of interest The authors declare that they have no conflict of
interest.

## References

  * [1] P. Azimzadeh, E. Bayraktar, High order Bellman equations and weakly chained diagonally dominant tensors, SIAM Journal on Matrix Analysis and Applications, 40(1)(2019): 276-298.
  * [2] X. Bai, H. He, C. Ling, G. Zhou, A nonnegativity preserving algorithm for multilinear systems with nonsingular ${\mathcal{M}}$-tensors, Numerical Algorithms, 87(3)(2021): 1301-1320.
  * [3] L. Cui, M. Li, Y. Song, Preconditioned tensor splitting iterations method for solving multi-linear systems, Applied Mathematics Letters, 96(2019): 89-94.
  * [4] W. Ding, L. Qi, Y. Wei, M-tensors and nonsingular M-tensors, Linear Algebra and Its Applications, 439(2013): 3264-3278.
  * [5] W. Ding, Y. Wei, Solving multi-linear system with M-tensors, Journal of Scientific Computing, 68 (2) (2016): 689-715.
  * [6] L. Han, A homotopy method for solving multilinear systems with M-tensors, Applied Mathematics Letters, 69(2017): 49-54.
  * [7] H. He, C. Ling, L. Qi, G. Zhou, A globally and quadratically convergent algorithm for solving multilinear systems with M-tensors, Journal of Scientific Computing, 76(3)(2018): 1718-1741.
  * [8] T. G. Kolda, B.W. Bader, Tensor decompositions and applications, SIAM Review, 51(2009): 455-500.
  * [9] D. Liu, W. Li, S. W. Vong, The tensor splitting with application to solve multi-linear systems, Journal of Computational and Applied Mathematics, 330 (2018), 75-94.
  * [10] D. Liu, W. Li, S. W. Vong, A new preconditioned SOR method for solving multi-linear systems with an M-tensor, Calcolo, 57(2)(2020): 1-18.
  * [11] W. Li, D. Liu, S.W. Vong, Comparison results for splitting iterations for solving multi-linear systems, Applied Numerical Mathematics, 134(2018): 105-121.
  * [12] D. Li, H.B. Guan, X.Z. Wang, Finding a nonnegative solution to an M-tensor equation, Pacific Journal of Optimization, 16(3)(2020): 419-440.
  * [13] D. Li, S. Xie, H. Xu, Splitting methods for tensor equations, Numerical Linear Algebra with Applications, 24(5)(2017): e2102.
  * [14] D. Li, J. Xu, H. Guan, Newton’s Method for M-Tensor Equations, Journal of Optimization Theory and Applications, 190(2)(2021): 628-649.
  * [15] L. Liu, X. Li, S. Liu, Further study on existence and uniqueness of positive solution for tensor equations, Applied Mathematics Letters, 124(2022): 107686.
  * [16] M. Liang, B. Zheng, Y. Zheng, R. Zhao, A two-step accelerated Levenberg–Marquardt method for solving multilinear systems in tensor-train format, Journal of Computational and Applied Mathematics, 382(2021): 113069.
  * [17] L.H. Lim, Singular values and eigenvalues of tensors: a variational approach, in: Proceedings of the IEEE International Workshop on ComputationalAdvances in Multi-Sensor Adaptive Processing, CAMSAP 05, vol.1, IEEE Computer Society Press, Piscataway, NJ. 2005: 129-C132.
  * [18] M. Liang, B. Zheng, R. Zhao, Alternating iterative methods for solving tensor equations with applications, Numerical Algorithms, 80(4)(2019):1437-1465.
  * [19] C. Q. Lv, C. F. Ma, A Levenberg-Marquardt method for solving semi-symmetric tensor equations, Journal of Computational and Applied Mathematics, 332(2018): 13-25.
  * [20] Q. Ni, L. Qi, A quadratically convergent algorithm for finding the largest eigenvalue of a nonnegative homogeneous polynomial map, Journal of Global Optimization, 61(4)(2015): 627-641.
  * [21] J. M. Ortega, W. C. Rheinboldt, Iterative solution of nonlinear equations in several variables. Classics Appl. Math. 30, Society for Industrial and Applied Mathematics, Philadelphia, 2000.
  * [22] L. Qi, Eigenvalues of a real supersymmetric tensor, Journal of Symbolic Computation, 40(2005): 1302-1324.
  * [23] L. Qi, Z. Luo, Tensor analysis: Spectral theory and special tensors, Society for Industrial and Applied Mathematics, 2017.
  * [24] J. Shao, L. You, On some properties of three different types of triangular blocked tensors, Linear Algebra and its Applications, 511(2016): 110-140.
  * [25] A. Toth, C. T. Kelley, Convergence analysis for Anderson acceleration, SIAM Journal on Numerical Analysis, 53(2)(2015): 805-819.
  * [26] X. Wang, M. Che, Y. Wei, Neural networks based approach solving multi-linear systems with M-tensors, Neurocomputing, 351(2019): 33-42.
  * [27] X. Wang, M. Che, Y. Wei, Global uniqueness and solvability of tensor complementarity problems for $\mathcal{H}_{+}$-tensors, Numerical Algorithms, 84(2)(2020): 567-590.
  * [28] X. Wang, M. Che, Y. Wei, Existence and uniqueness of positive solution for H+-tensor equations, Applied Mathematics Letters, 98(2019): 191-198.
  * [29] X. Wang, M. Che, Y. Wei, Preconditioned tensor splitting AOR iterative methods for $\mathcal{H}$-tensor equations, Numerical Linear Algebra with Applications, 27(6)(2020): e2329.
  * [30] X. Wang, C. Mo, M. Che , Y. Wei, Accelerated dynamical approaches for finding the unique positive solution of $\mathcal{K}\mathcal{S}$-tensor equations. Numerical Algorithms, 88(4)(2021): 1787-1810.
  * [31] H. F. Walker, P. Ni, Anderson acceleration for fixed-point iterations, SIAM Journal on Numerical Analysis, 49(4)(2011): 1715-1735.
  * [32] Z. Xie, X. Q. Jin, Y. Wei, Tensor methods for solving symmetric ${\mathcal{M}}$-tensor systems, Journal of Scientific Computing, 74(1)(2018): 412-425.
  * [33] L. Zhang, L. Qi, G. Zhou, M-tensors and some applications, SIAM Journal on Matrix Analysis and Applications, 35(2014): 437-452.
  * [34] Y. Zhang, Q. Liu, Z. Chen, Preconditioned Jacobi type method for solving multi-linear systems with M-tensors, Applied Mathematics Letters, 104(2020): 106287.


# Characterization of resistive Micromegas detectors for the upgrade of the
T2K Near Detector Time Projection Chambers

D. Attié  M. Batkiewicz-Kwasniak  P. Billoir  A. Blanchet  A. Blondel  S.
Bolognesi  D. Calvet  M.G. Catanesi  M. Cicerchia  G. Cogo  P. Colas  G.
Collazuol  A. Delbart  J. Dumarchez  S. Emery-Schrenk  M. Feltre  C. Giganti
F. Gramegna  M. Grassi  M. Guigue  P. Hamacher-Baumann  S.
<EMAIL_ADDRESS>F. Iacob  C. Jesús-Valls  R. Kurjata  M.
Lamoureux  M. Lehuraux  A. Longhin  T. Lux  L. Magaletti  T. Marchi  A. Maurel
L. Mellet  M. Mezzetto  L. Munteanu  Q. V. Nguyen  Y. Orain  M. Pari  J.-M.
Parraud  C. Pastore  A. Pepato  E. Pierre  B. Popov  H. Przybiliski  T.
Radermacher  E. Radicioni  M. Riallot  S. Roth  A. Rychter  L. Scomparin  J.
Steinmann  S<EMAIL_ADDRESS>J. Swierblewski  D.
Terront  N. Thamm  F. Toussenel  V. Valentino  G. Vasseur  U. Yevarouskaya  M.
Ziembicki  M. Zito  IRFU, CEA, Université Paris-Saclay, Gif-sur-Yvette, France
H. Niewodniczanski Institute of Nuclear Physics PAN, Cracow, Poland  LPNHE,
Sorbonne Université, Université de Paris, CNRS/IN2P3, Paris; France  INFN
sezione di Bari, Università di Bari e Politecnico di Bari, Italy  INFN:
Laboratori Nazionali di Legnaro (LNL), Padova , Italy  INFN Sezione di Padova
and Università di Padova, Dipartimento di Fisica e Astronomia, Padova, Italy
RWTH Aachen University, III. Physikalisches Institut, Aachen, Germany
Institut de Fisica d’Altes Energies (IFAE), The Barcelona Institute of Science
and Technology, Bellate Spain  Warsaw University of technology, Warsaw, Poland
Institute for Nuclear Research of the Russian Academy of Sciences, Moscow,
Russia

###### Abstract

The second phase of the T2K experiment is expected to start data taking in
autumn 2022. An upgrade of the Near Detector (ND280) is under development and
includes the construction of two new Time Projection Chambers called High-
Angle TPC (HA-TPC). The two endplates of these TPCs will be paved with eight
Micromegas type charge readout modules. The Micromegas detector charge
amplification structure uses a resistive anode to spread the charges over
several pads to improve the space point resolution. This innovative technique
is combined with the bulk-Micromegas technology to compose the ”Encapsulated
Resistive Anode Micromegas” detector. A prototype has been designed, built and
exposed to an electron beam at the DESY II test beam facility.
The data have been used to characterize the charge spreading and to produce a
RC map. Spatial resolution better than 600 µm and energy resolution better
than 9% are obtained for all incident angles. These performances fulfil the
requirements for the upgrade of the ND280 TPC.

###### keywords:

Resistive Micromegas, T2K Near Detector Time Projection Chambers

††journal: Nuclear Instruments and Methods

###### Contents

  1. 1 Introduction
  2. 2 Resistive Micromegas
  3. 3 Experimental setup
  4. 4 Collected data
  5. 5  Charge spreading characterisation
  6. 6 Track reconstruction
  7. 7 Spatial resolution
    1. 7.1 The Pad Response Function method
    2. 7.2 Spatial resolution estimation
    3. 7.3 Spatial resolution dependence on the drift distance, momentum, high voltage
    4. 7.4 Spatial resolution dependence on the track inclination
    5. 7.5 Bias measurements
  8. 8 Deposited energy resolution
    1. 8.1 Definition of the cluster charge
    2. 8.2 Deposited energy resolution for inclined tracks
    3. 8.3 Dependence of the $dE/dx$ resolution on the number of clusters
  9. 9 RC map calculation
  10. 10 Conclusions

##  1 Introduction

T2K [1] is a long-baseline neutrino oscillation experiment exploiting a muon
neutrino beam produced by the JPARC accelerator complex in Japan. The T2K
experiment includes the secondary neutrino beamline, a set of near detectors
(INGRID and ND280) and the far detector SuperKamiokande.

T2K provided the first evidence of non-zero mixing angle $\theta_{13}$ [2] and
discovered the appearance of electron neutrinos in a muon neutrino beam [2, 3,
4]. Combining T2K data with precise $\theta_{13}$ measurement from reactor
experiments, T2K has recently reported hints of large Charge-Parity (CP)
violation in the leptonic sector [5], excluding CP-conservation at about the
2$\sigma$ Confidence Level.

The T2K collaboration is now preparing for the second phase of the experiment
(T2K-II), starting in Fall 2022, which will exploit the upgrade of the beam
from 500 kW to 750 kW. T2K-II will collect in total more than $10^{22}$
Protons-On-Target (POT), including the $3.6\times 10^{21}$ POT already
collected, thus enabling 3$\sigma$ sensitivity on CP-Violation, in case of
maximally violated CP, as currently indicated by the T2K results [5]. In order
to cope with such increased statistics, improved control of the relevant
systematic uncertainties is needed. To this aim, an upgrade to ND280 is being
constructed.

ND280 is a magnetized multi-purpose detector that measures, before the
oscillation, the neutrino differential rate versus the kinematics variables of
the particles, mostly muons, that result from the neutrino interaction with
nuclei in dense parts of the detector. This allows one to constrain the
neutrino flux and the neutrino-nucleus interaction cross-section. The present
ND280 consists of two main parts: an upstream $\pi^{0}$ detector (P0D) and a
downstream tracker which includes two Fine Grained Scintillators (FGD)
interleaved with three vertical Time Projection Chambers (TPCs). The P0D and
the tracker are surrounded by an electromagnetic calorimeter (ECAL) and by the
UA1 magnet providing a 0.2 T magnetic field. The magnet yoke is further
instrumented with a Side Muon Range Detector (SMRD).

The role of ND280 in the T2K oscillation analysis is crucial, allowing one to
constrain the uncertainty on the expected number of events at the far detector
to 4-5%. ND280 measurements are performed on different targets (Carbon and
Oxygen) and rely on the precise measurement of the muon momentum measured by
the TPCs, with a momentum resolution of 10% at 1 GeV [6]. Better resolution is
not needed since the determination of the incoming neutrino energy from the
outgoing muon momentum in $\nu_{\mu}$ charged-current interactions is limited
by the uncertainty induced by the Fermi motion of the nucleon in the nucleus.

Another important result of ND280 is the measurement of the $\nu_{e}$
contamination in the beam [7, 8], that constitutes one of the two main
backgrounds to the $\nu_{e}$ appearance in the electron-like sample selected
at Super-Kamiokande. This analysis, based on the Particle Identification (PID)
capabilities of the TPCs and of ECAL, is possible thanks to the deposited
energy resolution of 8% in the TPCs that allows sufficient e-$\mu$ separation
between a few hundred MeV and $\sim$2 GeV.

The upgrade of ND280 consists in substituting the P0D with a new tracker,
similar to the existent one but with a horizontal orientation parallel to the
neutrino beam, and able to detect particles transverse to the beam, for which
the acceptance of the current tracker is small. The new tracker includes a
3-dimensional scintillator target for neutrino interactions(Super-FGD) made of
about 2 million 1 $\text{cm}{{}^{3}}$ cubes, read out by wavelength shifting
fibers in the 3 directions [9]. On the top and the bottom of the Super-FGD,
two High-Angle TPCs (HA-TPCs) will be installed. The new tracker system will
be surrounded by six Time-of-Flight modules. This new detector configuration
will allow one to improve the angular acceptance of ND280, being close to the
full 4$\pi$ phase-space accessible at SuperKamiokande. In addition, the better
tracking performances of the super-FGD will allow to improve the
reconstruction of the hadronic part of the neutrino interactions, that will be
exploited in combination with the muon kinematics. Also, more target mass for
neutrino interactions yields more data statistics.

Each endplate of the new TPCs will be instrumented with eight Micromegas type
charge readout modules. The Micromegas detector charge amplification structure
uses a resistive anode to ”spread” the charges over several pads to improve
the space point resolution: a signal is induced in pads adjacent to the pad
receiving directly electrons from the avalanche, and the information from
those additional signals allows one to better measure the track position. This
effect is different from the charge sharing between pads resulting from the
diffusion of the primary electrons in the TPC gas. This innovative technique
is combined with the bulk-Micromegas technology [6] to compose the
”Encapsulated Resistive Anode Micromegas” (ERAM) detector. Performances of a
prototype of an ERAM detector exposed to a test beam at CERN were shown in
[10]. This paper describes the performance of one prototype ERAM module, built
with the same design that will be used for the HA-TPCs, exposed to an electron
beam at DESY. With respect to our previous studies reported in [10], in this
paper we characterize both spatial and dE/dx resolution as a function of the
angle of the track projection reconstructed in the ERAM plane, with the side
of the ERAM plane, corresponding to the sides of the ERAM pads.

The HA-TPC should match the performance of existing TPCs in terms of momentum
and dE/dx resolution. In addition, the super-FGD will enable the
reconstruction of low momentum protons and neutrons [11], thus allowing one to
reconstruct the incoming neutrino energy more precisely and effectively
correct for the Fermi motion: an improved track momentum resolution in the HA-
TPC, even beyond the previous specifications, is therefore useful to improve
the neutrino energy resolution and, ultimately, the precision on neutrino
oscillation measurements.

##  2 Resistive Micromegas

The ERAM detector uses the bulk-Micromegas technology invented in 2004 by a
collaboration between the CERN/EP-DT-EF PCB workshop and CEA-IRFU [12]. It was
developed and used for the construction of the 72 bulk-Micromegas modules (9
$\textrm{m}^{2}$ total area) which equip the three T2K/ND280 vertical TPCs.
The woven micromesh of the Micromegas detector is sandwiched between two
layers of the same insulating material (pyralux) and exposed to UV radiation
at the location where the pyralux is kept on top of the pad-segmented anode
PCB after chemical development.

The ERAM detector is a 128 µm amplification gap bulk-Micromegas using the
standard SD45/18 304L woven micromesh built on top of a resistive anode PCB.
When a track crosses the gas volume of a TPC, it generates a cloud of ionized
electrons. These electrons are drifted to the anode readout plane of the TPC
under a uniform electric field. On the readout plane, an avalanche is
generated by a high electric field in the Micromegas amplification region,
called DLC voltage in the following of this paper. The 3D position of the
track is then reconstructed from the arrival time and the position of the
avalanche on the readout plane. In the case of the bulk-Micromegas, the
electron avalanche in the amplification gap is quite narrow with respect to
the pad size and therefore the position resolution is often limited by the pad
size. The resistive anode Micromegas, first introduced and extensively studied
by the ILC-TPC collaboration [13], provides a way to induce a signal on a
larger number of pads allowing for a better reconstruction of the track
position. A sketch of the bulk and the resistive Micromegas concepts are
presented in Figure 1.

Figure 1:  A sketch of the Standard Bulk Micromegas layout (left) and the
Resistive anode Micromegas layout (right).

The anode, segmented in pads, is covered by a foil of insulating material,
which has a thin resistive layer on its top side. The ERAM detector uses a
Diamond-Like Carbon (DLC) thin layer sputtered on a 50 ${\rm\text{\textmu}m}$
thick APICAL (Kapton) insulator sheet. The avalanche is then naturally
quenched because the potential difference locally drops in presence of a high
charge density. The resistive layer acts like a 2-D RC network and the charge
deposited by the avalanche induces a signal also on the adjacent pads, looking
as if the charge spreads naturally with time following a Gaussian behaviour.
For a point charge deposited at $r=0$ and $t=0$, the charge density as a
function of radius $r$ and time $t$ is given by:

$\rho(r,t)=\frac{RC}{4\pi t}e^{\frac{-r^{2}RC}{4t}}$ (1)

where $R$ is the resistivity per unit area and $C$ is the capacitance per unit
area. For this structure, the capacitance $C$ is defined by the distance
between the resistive layer and the grounded pads (glue thickness plus APICAL
foil). The insulating layer thickness determines the width of the induced
charge spread $\sigma_{t}$ at a given time. When measured, due to integration
of the charge collected by a front-end electronics of shaping time $t$, this
spread is of the order of $\sigma_{t}=\sqrt{\frac{2t}{RC}}$.

The resistive anode provides mainly two advantages: by spreading the charge
between neighbouring pads, it improves greatly the resolution with respect to
the pitch$/\sqrt{12}$ provided by a mere hodoscope, and it suppresses the
formation of sparks and limits their intensity. A further and novel
improvement of this technique is a new High Voltage powering scheme, where the
mesh is set to ground and the anode to a positive amplification voltage. The
insulation of the resistive anode from the pads, hence from the electronics,
ensures a safe operation by a capacitive coupling readout and thus allows us
to remove the cumbersome anti-spark protection circuitry of electronic read-
out boards necessary in the case of the standard bulk readout.

At the end of 2017, the first series of prototypes were produced and tested to
assess the feasibility of large charge spreading by a low resistivity anode.
The vertical TPCs PCB, with an active area of 36$\times$34 cm2 covered by
0.97$\times$0.69 cm2 pads, was adapted to build an ERAM structure with ”on-
shelf” 2.5 MOhm/$\Box$ DLC foils. This first prototype was tested with a
particle beam at CERN in 2018 and its performances are summarized in [10].

In fall 2018, the global design of the ND280 upgrade detector was fixed. The
sub-detector envelops were defined and the size of the ERAM modules fixed to
be $420\times 340$ mm2 with $32\times 36$ rectangular pads of size
$10.09\times 11.18$ mm2. The ERAM module studied in this paper has a
resistivity close to the required one of 200 kOhm/$\Box$ using DLC foils stack
on a 75 ${\rm\text{\textmu}m}$ glue layer. The resistivity maps measured at
critical steps of the manufacturing of the DLC foil shows a non uniformity of
20%. The detector is readout with two analogue 576 channels AFTER based Front-
End Cards (eight 72 channels AFTER ASIC per card). The AFTER chip, already
used for the existing TPCs of T2K and for the test beam described in [10]
allows tuning several parameters such as the gain, the shaping time and the
sampling time. This prototype was tested in a test beam at DESY in June 2019.

##  3 Experimental setup

The prototype has been exposed to an electron beam at the DESY II test beam
facility [14]. DESY II provides electron beams of $1-6~{}\rm{GeV/c}$ at a rate
of up to several kHz, depending on the chosen momentum. In the test beam area
TB24/1, a large-bore superconducting solenoid, called PCMAG, provides a
magnetic field of up to $1.25~{}\rm{T}$. The magnet is mounted on a movable
platform, which allows the setup to be moved horizontally and vertically,
perpendicular to the beam line, as well as rotated by
$\pm$45\text{\,}\mathrm{\SIUnitSymbolDegree}$$ in the horizontal plane. The
platform can position the device under test with a precision of about
$0.2\text{\,}\mathrm{mm}$ horizontally, $0.1\text{\,}\mathrm{mm}$ vertically,
and within $0.1\text{\,}\mathrm{\SIUnitSymbolDegree}$ in angle.

Inside the bore of the magnet, a rail system is installed on which test
devices can be mounted at different positions within the magnet. The TPC
prototype is supported on a sled, which can move in and out of the magnet and
can be used to rotate the chamber around the magnetic field axis as
illustrated in Figure 2.

Usually, the magnet is positioned perpendicular to the beam, so the magnetic
field is also perpendicular to the beam, in a horizontal plane. The ERAM
module is a vertical plane parallel to the beam, but perpendicular to the
magnetic field. The walls of the magnet present about
$20\text{\,}\mathrm{\char 37\relax}$ of a radiation length, so that an
electron beam of $6~{}\rm{GeV/c}$ easily penetrates the magnet and the device
under test. A set of four consecutive scintillation counters, of which each
has an area of approximately $2.5\text{\times}2.5\text{\,}\mathrm{cm}$, is
mounted about $1.5\text{\,}\mathrm{m}$ in front of the magnet. The coincidence
between them is used as a beam trigger. In addition, a second set of
scintillation counters above and below the magnet provides a cosmic trigger
for tests without beam.

The prototype has been placed inside the magnetic field of 0.2 T provided by
the PCMAG magnet and has been operated for these measurements with a gas
mixture of 95% argon, 3% tetrafluoromethane (CF4), and 2% isobutane (iC4H10)
that is the same gas mixture as the one used in the existing TPCs of ND280.
The gas quality was constantly monitored during the measurement. For the
results reported in this paper, the oxygen contamination was around
$30~{}\rm{ppm}$ at a gas flow rate of $30~{}\ell/h$. The chamber was operated
at atmospheric pressure. Ambient temperature and pressure were constantly
monitored.

Figure 2: View of the setup.

##  4 Collected data

The collected data is used to analyze in detail the role of all the relevant
parameters (electronics peaking time, DLC voltage, drift distance) to fully
characterize the charge spreading, the resistive foil uniformity, and to
ensure a performance satisfying the ND280 upgrade requirements. The tests at
DESY were done in a short chamber with a 15 cm drift distance. The prototype
was operated at a voltage of 360 V. The settings chosen for the AFTER chip
were a sampling time of 40 ns, a peaking time of 412 ns or 200 ns, and a gain
of 120 fC.

The results presented in this paper were obtained with electrons with momenta
varying from 0.5 to 5 GeV/c. We have carried out drift distance scans with
seven points, spaced by 2 cm, at $B=0~{}\textrm{T}$ and $B=0.2~{}\textrm{T}$,
peaking time of 200 ns or 412 ns, and high voltage of 370 V and 360 V. In
addition a scan of the DLC voltage, varied from 330 V to 400 V, was performed.

We have calibrated the $T_{0}$, the moving table position, and the drift
velocities at different drift fields with a 4 GeV/c electron beam and a short
peaking time of 116 ns. At the standard T2K field of 275 V/cm, we obtained a
drift velocity of $7.68\pm 0.03~{}{\rm cm}/\text{\textmu}\textrm{s}$ by drift
distance scans with the accelerator beam. A Gas Monitoring Chamber (GMC),
identical to the ones deployed at T2K’s ND280 detector [6], monitored the
exhaust gas for the duration of the beam-time. The GMC measured a drift
velocity of $7.81\pm 0.02~{}{\rm cm}/text{\textmu}\textrm{s}$ at the T2K
field, in agreement with the beam scans, see Figure 3. Under the electric
field of 140 V/cm, which is associated with the minimum transverse diffusion,
we found a drift velocity of $5.84~{}{\rm cm}/\mu\textrm{s}$. Despite the
impact of the gas bottle changes on the gas properties, the drift velocity
$V_{\rm drift}$ under the electric field of $E=275{\rm V/cm}$ varies less than
$6\text{\textperthousand}$.

a

b

Figure 3:  The drift velocity under a constant field of 275 V/cm as a function
of days (a) and as a function of electric field (b) with density corrections
applied. The vertical lines on (a) correspond to the moments when the bottles
of premixed gas were replaced. Transition regions at ramp-up and between
bottle changes can clearly be seen. The measurements at different fields in
(b) correspond to a full day during the first bottle (14.06).

##  5  Charge spreading characterisation

As explained before, the resistive Micromegas technology produces a spreading
of the collected charge into neighbouring pads. The charge spreading
phenomenon, which drives the waveform shape on the pads adjacent the ”leading”
ones receiving the inital charge from the avalanche, is described in [13]. The
signal induced by the resistive layer is smaller in amplitude, delayed and
longer in time compared to the signal from direct charge deposition from the
track in the leading pad. The time delay of the waveform increases with the
distance of the adjacent pad to the track. Hence, the charge spreading is
significant in the transverse direction w.r.t. the track projection
reconstructed in the ERAM plane, while in the longitudinal direction it is
masked by the direct charge. To study the phenomenon of interest, we define a
“cluster” as a group of pads in the perpendicular direction to the track. A
schematic view of such a cluster for horizontal tracks and waveforms in the
leading and in adjacent pads are shown in Figure 4.

For each pad, the maximum of the waveform is used as an estimator of the
charge, which corresponds to different times for the leading and adjacent
pads. In Figure 5, the pad multiplicity per cluster and the ratio of the
charge measured in the pad with the largest signal with respect to the sum of
the charges of the pads in the cluster
($q_{\textup{max}}/q_{\textup{cluster}}$) are shown.

Most of the clusters are formed by more than two pads and the pad with the
largest signal contains typically 80% of the total charge of the cluster. The
effect of the high voltage on the multiplicity is also clearly seen in Figure
5. The cluster multiplicity increases with the high voltage, as the
probability for the smaller signals in some pads to pass the threshold
increases with the gain.

Figure 4: A schematic view of a column cluster (left) and the waveforms of
each pad composing this cluster (right) for 412 ns shaping time.

Figure 5:  The number of pads in a cluster versus the DLC high voltage (left)
and the fraction of the cluster charge collected in the pad with largest
signal at 360 V (right).

For tracks projections in the pad plane which are parallel to the sides of the
pads, transverse spreading is defined precisely within the given column. While
for inclined tracks and large square pads, the separation between longitudinal
(along the track) and transverse spreading is more complicated. In order to
distinguish these two spreading topologies, we define more sophisticated
cluster patterns shown in Figure 6. These clusters are repeated to pave the
whole ERAM. In the case of square pads, these patterns are optimized for
angles with tangents 0 (column), 1 (diagonal), 0.5 (2 by 1), 0.3 (3 by 1)
respectively. In the case of rectangular pads, the optimal angles are slightly
different.

a

b

c

d

Figure 6: The different cluster patterns that can be used depending on the
track (red line) angle with the pads sides : (a) column, (b) diagonal, (c) 2
by 1, (d) 3 by 1. For each case, the colored pads correspond to one cluster,
associated to the colored leading pad crossed by the track.

##  6 Track reconstruction

In the test beam analysis, we focus on the studies of the through-going tracks
as more complicated typologies (e.g. showering, curved low-energy tracks) are
difficult to interpret. Hence simple reconstruction algorithms like DBSCAN
[15] are sufficient in our case.

We select a track if it is crossing the whole detector without breaks or
splits. A split is defined as the case where there is more than one cluster in
a given column. An event containing a split is thus a multiple track candidate
and is rejected in our analysis. However, with the resistive spreading, two
close parallel tracks may not be separated by a gap and thus misreconstructed
as one single track. To reject such a topology, a cut on the pad multiplicity
in each cluster was implemented. The cut value was optimised for each
Micromegas voltage and electronics shaping time. Examples of the accepted and
rejected events are shown in Figure 7 (a) and (b) respectively.

a

b

Figure 7: Event displays of the (a) single track and (b) multi-track in the
prototype.

##  7 Spatial resolution

As described in section 1, the momentum resolution has to be better than 10%
at 1 GeV/c but, to fully exploit the ND280 upgrade capabilities, an even
better momentum resolution would be desirable. The momentum resolution is
directly connected to the spatial resolution obtained for one cluster through
the Gluckstern formula [16]. For tracks with 70 point measurements (using
clusters), a maximum drift distance of 1 m, and a magnetic field of 0.2 T, a
spatial resolution of $\sim 800\text{\textmu}\textrm{m}$ would be sufficient
to reach a momentum resolution of 10% at 1 GeV/c.

As we will show in this section, the resistive Micromegas technology allows
one to significantly improve the spatial resolution with respect to the bulk
Micromegas, even in presence of slightly larger pads thus allowing one to
reduce by $\sim 30\%$ the total number of electronic channels for the same
active surface. The test beam data have been used to characterize the ERAM
module performances for electrons with different angles with the pads sides.

###  7.1 The Pad Response Function method

The charge spread described in the section 5 results in charge detection in a
few pads around the avalanche arrival point. The charge measurements are
discrete with the finite pad size, while the spreading in the RC layer is
continuous. Moreover, the signal induced in the adjacent pads does not result
from a real sharing of an initial charge, and it is measured at different
times, even if it is related to the charge measured in the leading pad. Thus
the barycentric method (Center of Charge, weighted mean) that assigns all the
collected charge to the pad centre doesn’t provide a precise position
reconstruction. Instead, we used the so-called Pad Response Function (PRF)
which characterizes the relation between observed charge ratios and track
position w.r.t. the pad (Equation 2). This method improves the spatial
resolution compared to the barycentric method for TPCs with the resistive
anode [17]. The PRF is defined as:

$PRF(x_{track}-x_{pad})=Q_{pad}/Q_{cluster}$ (2)

where $x_{track}$ is the reconstructed position of the track, $x_{pad}$ is the
centre of the pad, $Q_{pad}$ is the charge collected on a given pad and
$Q_{cluster}$ is the charge collected on the whole cluster. The definition of
the cluster is the same as described in section 5: it’s a group of pads where
one receives a charge from the initial avalanche and the others detect the
charge spread in the resistive foil.

To parametrize the PRF we used the empirical ratio of two symmetric $4^{th}$
order polynomials proposed in [18]:

$PRF(x,\Gamma,\Delta,a,b)=A\times\dfrac{1+a_{2}x^{2}+a_{4}x^{4}}{1+b_{2}x^{2}+b_{4}x^{4}}$
(3)

where the parameters $a_{i}$ and $b_{i}$ can be related to the more physical
parameters: the full width at half maximum $\Gamma$, the base width $\Delta$,
and two scaling parameters $a$ and $b$.

###  7.2 Spatial resolution estimation

This parametrized function of the PRF allows one to determine the parameters
using an iterative method. To get a first guess of the PRF parameters, we use
the track position reconstruction obtained from the barycentric method. After
having defined the position of all the clusters, the global track is fit with
a parabola. The fit is based on many measurements using the different clusters
along the track ($\geqslant 34$) thus it is considered as a true track
position and the PRF scattered plot is filled (Figure 8 (a)). The scattered
plot is profiled along Y axis to form a graph that is further fit with the
analytical function from Equation 3 (Figure 8 (b)).

a

b

Figure 8: The Pad Response Function (PRF) obtained with (a) scattered plot and
(b) results of the profile and fit with analytical function. The pad borders
are represented by vertical lines.

After having determined the PRF, the track position in each cluster is
obtained with the following $\chi^{2}$ minimization:

$\chi^{2}=\sum_{pads}\frac{Q_{pad}/Q_{cluster}-PRF\left(x_{track}-x_{pad}\right)}{\sigma_{Q_{pad}/Q_{cluster}}}$
(4)

where $\sigma_{Q}$ is the uncertainty on the charge measurements. In our
analysis, we assume that charge measurement probability follows a Poisson
distribution, hence:
$\sigma_{Q_{pad}/Q_{cluster}}=\sqrt{Q_{pad}/Q_{cluster}}$. We proceed through
the iterative process of the PRF estimation until the track fit quality is not
improving anymore. Typically this procedure converges after few iterations and
for the results shown in this paper 10 iterations were used.

The spatial resolution is defined as the difference between the reconstructed
position in a given cluster and the track global fit (residual). The
particular cluster where the resolution is studied is excluded from the fit to
prevent biases. The residuals distribution is fit with a Gaussian function
whose standard deviation defines the spatial resolution.

###  7.3 Spatial resolution dependence on the drift distance, momentum, high
voltage

The spatial resolution was studied over different samples. The beam position
was varied within the drift distance of the field cage keeping the tracks
parallel to the Micromegas plane. For resistive Micromegas, the spatial
resolution is expected to degrade slightly for a larger drift distance
affected by the transverse and longitudinal diffusion. The observed
dependence, for a DLC voltage of 360 V, is shown in Figure 9 (a). A resolution
between 200 and 250 µm is observed for the whole drift length.

In addition, the high voltage applied to the Micromegas mesh was varied to
study the detector performance in different regimes. Higher voltage is
expected to enhance the initial avalanche, thus increase the smaller charge
spreading signal as well. Signals in the neighbour pads are then more likely
to pass above the threshold and are less affected by statistical fluctuations.
As track position reconstruction relies on the charge spreading measurements
we expect better performance with higher voltage. The results are shown in
Figure 9 (b).

a

b

Figure 9: Spatial resolution w.r.t. (a) beam injection position and (b)
Micromegas voltage for horizontal tracks parallel to the MM plane. Points
represent the mean value over detector columns and errors represent the
fluctuations (RMS).

The DESY beamline allows changing the momentum of the electrons delivered to
the test beam area. In this way, we studied spatial resolution as a function
of the track kinematics. We find no significant changes for the position
accuracy reconstruction for the tracks in the range between 1 and 5 GeV/c
(Figure 10).

Figure 10: Spatial resolution as a function of the momentum of the electron
beam.

###  7.4 Spatial resolution dependence on the track inclination

Inclined tracks are expected to be reconstructed less precisely compared to
horizontal ones. As an example, in the current ND280 TPCs, the resolution
degrades as a function of the track angle from
$600~{}\text{\textmu}\textrm{m}$ to $\sim 1~{}\textrm{mm}$ [6]. It is then
particularly interesting to investigate the behaviour of the spatial
resolution in the ERAM detector as a function of the angle of the
reconstructed track on the ERAM plane with respect to the ERAM side (track
inclination).

In order to do this, the different cluster patterns described in section 5
were used for different track inclination. The $\chi^{2}$ fit (Equation 4) is
applied to each cluster to extract the track position. As for horizontal
tracks, the positions in given clusters are fit together to form a global fit
and the iterative analysis is applied: the barycentric estimation is used as a
prior following with PRF calibration. The results of the inclined track
spatial resolution estimations are shown in Figure 11.

Figure 11: Spatial resolution w.r.t. track angle within the MM side obtained
with different cluster definitions.

As expected, the column clustering method leads to a severe spatial resolution
degradation with the track slope, reaching a maximum of $1.4~{}mm$ for 45∘
tracks. For these angles, the diagonal pattern provides a significant
performance improvement. In the intermediate regions (20∘, 30∘, 60∘, 70∘), the
best result is achieved with the more complex patterns: “2 by 1” and “3 by 1”.
The asymmetry w.r.t. 45∘ is caused by the rectangular pad shape
11.3$\times$10.2 mm. Thus, the diagonal pattern is considered a better choice
for 48∘ tracks than for 45∘ tracks. Hence, tracks inclined with 50∘ are
reconstructed more accurately compared to 40∘ tracks. Similar behavior is
observed for all the other patterns.

By taking the best clustering algorithm we observe a spatial resolution better
than $600~{}\text{\textmu}\textrm{m}$ for all the angles. We understand the
difference between horizontal tracks and inclined tracks as due to the larger
effective pad size for diagonal clustering and to the rectangular shape of the
pads, while diagonal clustering would work better for square pads. For a
spatial resolution of $600~{}\text{\textmu}\textrm{m}$, for 70 point
measurements and a magnetic field of 0.2 T we expect a momentum resolution of
6% at 1 GeV/c that scales linearly with the spatial resolution.

###  7.5 Bias measurements

As described above, we define the spatial resolution as a standard deviation
of the difference between the reconstructed position in a given cluster and a
global track fit. Meanwhile, the mean value of the residuals is also an
important characteristic that shows the bias of our measurements. In
particular, it is interesting to study the biases with respect to the track
position in the pad. For that, we use the natural beam spread. The electron
beam profile is nearly Gaussian with a standard deviation $\approx$ 1 cm. We
sample the residuals with the reconstructed track position in the pad. Thus,
we can study the resolution and biases in the different pad regions. In Figure
12, we show both the spatial resolution and bias per column. Individual PRFs
were used for each column to analyze the behaviour in the different regions of
the detector independently. The resolution undergoes some oscillations because
if a track is close to one of the leading pad borders, the neighbouring pads
see a larger signal and we thus have a more reliable input for the position
reconstruction.

Figure 12: The resolution and the bias of the track reconstruction over the
position in the pad for various columns. The pad borders are represented with
vertical lines. In this coordinate system, the beam is centered around
$Y=-0.5cm$, but tracks are also measured on the neighbour pads due to the beam
spread.

We generalize the bias study for the whole detector. Figure 13 represents the
fluctuations of the resolution and biases in the given column. We conclude
that for most of the detector the biases are under control and smaller
compared to the spatial resolution. In the downstream part of the detector, we
found larger biases that could be related to the non-uniformities in the
resistivity of this ERAM detector that will be described in section 9.

Figure 13: The spatial resolution and bias fluctuations observed for the
different position of the track in the pad. Dots represent the mean value in a
given column, filled areas correspond to the RMS, and error bars represent
minimum and maximum values.

##  8 Deposited energy resolution

One of the main goals of a TPC is to perform particle identification (PID)
based on the measurement of the ionization produced by charged particles
crossing the gas volume. The PID capabilities depend on the resolution in the
ionization energy loss measurements.

In the case of T2K TPCs, the PID is mainly used to distinguish electrons
(produced by $\nu_{e}$) from muons (produced by $\nu_{\mu}$). In the momentum
range studied by T2K, the amount of ionization between electrons and muons
differs by $\sim 40\%$. Therefore, a resolution of less than 10% is needed to
efficiently distinguish these two particles. In general, the resolution
depends on the number of independent ionization measurements (i.e. the number
of clusters) and on the amount of ionization in each cluster. For the existing
TPCs, a resolution of $8\%$ was obtained by combining measurements in 2
Micromegas detectors (72 independent measurements of ionization).

In this section, we will describe the performances observed with one single
ERAM detector (36$\times$32 pads).

The method used to estimate the energy loss of a given track is called the
truncated mean method: the charges contained in each cluster of the track are
sorted by increasing order and only a fraction of the lowest charges is kept
to compute the mean deposited energy per track. Such a method allows to reject
clusters with a large amount of charge, due to fluctuations in the ionization
processes, that would degrade the relative resolution on the mean value, and
thus the power to separate different types of particles.

The dependence of the $dE/dx$ resolution on the truncation fraction is shown
in Figure 14. The best resolution is obtained for values of truncation
fraction between 50% and 70%. Therefore, for all the results presented in this
section, a truncation factor of 70% is used.

Figure 14: The dependence of the $dE/dx$ resolution on the truncation factor
$\alpha$ for the beam parallel to the pad sides with a magnetic field of 0.2 T
and electronics peaking times of either 200 ns and 412 ns

###  8.1 Definition of the cluster charge

As explained before, the basic ingredient of the $dE/dx$ resolution is the
amount of charge seen in each cluster. This quantity can be defined in
different ways.

In the existing ND280 TPCs or in the results published in [10], the cluster
charge is defined by summing the waveform maximum seen in each pad composing
the cluster. This definition was shown to be a good estimator of the charge in
the case of bulk Micromegas, or when the drift distance is large enough to
allow the predominance of transverse diffusion over the charge spreading
induced in the ERAM module.

However, in the present configuration, the limited size of the TPC implies
that the transverse diffusion is small and that charge spreading dominates,
except when the track is close to the pad border. Therefore, by summing the
maximum of the waveform seen in each pad of the cluster, we are double
counting the leading pad charge: first, it is seen in the leading pad before
the spread, and then seen by the neighboring pads after the spread.

To remove this double counting effect, the charge per cluster can be defined
in a different way. For each cluster, we build a summed waveform defined as
the sum of the waveform amplitude seen in each pad of the cluster at all
times. We then take as a cluster charge estimator the maximum of this summed
waveform. We refer to this method as WFsum.

A comparison of the results of the two methods as a function of the drift
distance is shown in Figure 15 for the two different peaking times and for
data taken with 0.2 T magnetic field and with the beam parallel to the ERAM
side. It can be seen that with both methods, the resolution is well below 10%
for all the drift distances and, as expected, the WFsum method gives a better
resolution comprised between 8.5% and 9.0%.

a

b

Figure 15:  dE/dx resolution with respect to the drift distance for the beam
parallel to the pad side with a magnetic field of 0.2 T and peaking times of
200 ns (a) and 412 ns (b). Qsum method consists in summing the maximum of the
waveform in each pad of the cluster while WFsum corresponds to maximum of the
summed waveforms in a cluster.

###  8.2 Deposited energy resolution for inclined tracks

As explained in section 5, in order to reconstruct inclined tracks, different
clustering algorithms are used. In the case of the deposited energy
resolution, the usage of such algorithms has two impacts: a larger number of
clusters per track are reconstructed but the track will have different paths
in different clusters.

In the column or row clusters (defined in Figure 6), tracks with the same
angle with respect to the pad sides have the same $dx$ in each cluster
(neglecting the curvature induced by the magnetic field). This is not true for
diagonal clusters in which the $dx$ can vary between 0 and the diagonal of the
pad ($\sim 1.5~{}\textrm{cm}$). The distribution of the charge as a function
of $dx$ for tracks inclined by 45 degrees with respect to the pad plane is
shown in Figure 16. It is clear from this figure that there is a dependence,
although this dependence is not linear as one would expect from the simple
consideration that the deposited energy should be proportional to the path.

The non-linearity is due to the fact that each cluster sees not only the
direct charge due to the primary ionization, but it also sees some charge due
to the spread on the resistive plane and to the transverse diffusion in the
gas.

Figure 16: The distribution of the charge with respect to the $dx$ of the
track inclined by 45 degrees with respect to the pad plane

In order to correct the $dx$ in each cluster, we fit the charge in each slice
of $dx$ with a Landau function and we take the Most Probable Value (MPV). The
distribution of MPV as a function of $dx$ is then parametrized with a third
degree polynomial. In each cluster, the charge is corrected to take into
account the real path length $dx$ and then the truncated mean is computed.

The $dE/dx$ distributions for diagonal clustering with and without the $dx$
correction are shown in Figure 17.

a

b

Figure 17: The $dE/dx$ distribution for 40 deg inclined tracks at an
electronics peaking time of 200 ns (a) and 412 ns (b) with and without the
correction for $dx$.

The deposited energy resolution as a function of the beam inclination w.r.t.
the pad plane is shown in Figure 18. As expected, the diagonal clustering,
after proper $dx$ correction, provides the best resolution thanks to the
larger amount of clusters in which the track is sampled.

a

b

Figure 18: $dE/dx$ resolution versus the angle w.r.t the pad plane using
column/row clustering or diagonal clustering. Column clustering is used from 0
to 40 deg and after 40 deg row clustering is used. Runs at a peaking time of
200 ns (a) and 412 ns (b), with a magnetic field of 0.2 T applied to the TPC
prototype. Diagonal clustering is corrected for the $dx$ as described in the
text.

###  8.3 Dependence of the $dE/dx$ resolution on the number of clusters

In this test beam, only one ERAM module was used. In the HA-TPC that will be
installed at ND280, most of the tracks will cross two ERAM modules before
exiting the TPC resulting in a larger number of clusters (72 for tracks
parallel to the pad plane). The observed dependence of the deposited energy on
the number of clusters can then be used to extrapolate the expected resolution
in the HA-TPCs.

This can be done by computing the truncated mean using only a fraction of the
available clusters. The dependence of the truncated mean on the number of
clusters used for different samples is shown in Figure 19.

a

b

Figure 19: $dE/dx$ resolution versus the number of clusters for runs at a
peaking time of 412 ns with a magnetic field of 0.2T applied to the TPC
prototype. Diagonal clustering (including $dx$ correction) is used for 45 deg
inclined tracks (a). In (b), column clustering is employed for 0 deg
horizontal tracks and row clustering for 90 deg vertical tracks (b).

The resulting distribution of the deposited energy resolution as a function of
the number of clusters $N$ is then fit with the function: $f(N)=AN^{-B}$. $B$
equal to 0.5 would correspond to a simple $\sqrt{N}$ dependence.

In the case of horizontal or vertical tracks, we observe similar behavior: the
horizontal tracks being slightly better because of the larger pad size in the
horizontal direction. Extrapolated to two ERAM modules, we obtain a deposited
energy resolution of the order of 6% for all the angles.

##  9 RC map calculation

The quantity controlling the charge spreading over time is the product $RC$
(section 2. To have a better understanding of our detector, we reconstruct the
map of $RC$ using horizontal tracks in data where a scan in y-direction was
performed at peaking time of 412 ns. This map is crucial to characterize our
detector and its uniformity and is also needed for a detailed detector
simulation.
To extract these $RC$ values, we use an analytical model of the charge
dispersion, adjusted to the waveforms measured in the pads. The induced charge
on a rectangular pad below the resistive layer can be calculated by
integrating the charge density function over the pad area [19] :

$\begin{split}\mathcal{Q}(t)=\frac{q_{e}}{4}\left[erf(\frac{x_{\textup{high}}-x_{0}}{\sqrt{2}\sigma(t)})-erf(\frac{x_{\textup{low}}-x_{0}}{\sqrt{2}\sigma(t)})\right]\times\\\
\left[erf(\frac{y_{\textup{high}}-y_{0}}{\sqrt{2}\sigma(t)})-erf(\frac{y_{\textup{low}}-y_{0}}{\sqrt{2}\sigma(t)})\right]\end{split}$
(5)

with $q_{e}$ is the initial charge, ($x_{0}$,$y_{0}$) the track position,
$x_{\textup{high}}$, $x_{\textup{low}}$, $y_{\textup{high}}$,
$y_{\textup{low}}$ the pad boundaries. In the denominator
$\sigma(t)=\sqrt{(2t/\tau)+\omega^{2}}$, the term $\tau=RC$ where $R$ is the
surface resistivity of the layer and $C$ the capacitance determined by the
spacing between the anode and readout planes. Finally, $\omega$ is associated
to the transverse diffusion term.

To compare to data, the characteristics of the front-end charge preamplifiers
need also to be included. Longitudinal diffusion increases the size of
electron charge clusters in the drift direction. The longitudinal diffusion is
neglected here since we have only 15 cm drift distance. The parameterization
of the electronics shaping time effects $\mathcal{I}(t)$ is obtained from the
simulation. The convolution of $\mathcal{I}(t)$ and $\mathcal{Q}(t)$ results
in the full theoretical model, which is compared to the data. This convolution
is handled numerically. The fit is based on clusters of pads perpendicular to
the track. Each cluster consists of a so-called leading pad collecting
essentially the initial charge deposit, and in so-called neighbour pads
sensitive mostly to the induced charge due to the resistive effect, as shown
in Figure 4.

The fit procedure is as follows: we first fit the leading pad waveform with
the electronics response function $\mathcal{I}(t)$, then we fit simultaneously
the two neighbouring pads waveform with a convolution of $\mathcal{I}(t)$ with
$\mathcal{Q}(t)$ in order to extract $RC$. Then a simultaneous fit of two
neighbouring pads waveform use separate $\mathcal{Q}(t)$ functions, as the
distance to the track can be different for the two pads, but we consider $RC$
and the electronics response parameters as common in the fit. Since we are
using horizontal tracks, we can only fit the $y_{0}$-position of the track and
we do not have any constraints on its $x$-position. The track position $y_{0}$
is obtained with the PRF $\chi^{2}$ minimization method (see equation (4)).
Figure 20 shows an example of waveform fit results for the leading pad and its
neighbours.

a

b

c

Figure 20: Example of waveform fit results for the leading pad (b) and its
neighbors (a) and (c) in a given cluster.

Another method can be used to cross-check the $RC$ values obtained with the
analytical model. This alternative method studies two parameters related to
the signal propagation: the time at which the signal in neighbour pads is
maximal and the ratio of amplitudes of the neighbour pads and the leading pad.
The time difference between the leading pad and one of its neighbour pad is
found to be proportional to $RC$ as can be seen in the following formula:
$\Delta t_{1}-\Delta t_{2}=RC\times L\times y_{0}$, with $L$ the pad length
and $y_{0}$ the track position and $\Delta t_{1,2}=t_{\textup{Leading
Pad}}-t_{\textup{Neighbor}}$ the time difference between the leading pad and
the neighbor pad.

The $RC$ maps obtained using the analytical model and the time difference
method described above are shown in Figure 21. The $RC$ values are given in
unit of $\textrm{ns}/\textrm{mm}^{2}$. Both methods give results of the same
order of magnitude though the second method is less precise than the fit
procedure. A non-uniformity of $RC$ up to 30% is observed using both methods.

This non-uniformity is confirmed with the observation of the charge collected
in the second to leading pad. The charge observed in the leading pad is
uniform across the detector while the measurements in the neighbours
$Q_{\textup{second}}/Q_{\textup{leading}}$ demonstrate fluctuations,
especially in the downstream detector region (Figure 22). Lower charge
fraction is consistent with the higher $RC$ value.

Figure 21: The RC ($\textrm{ns/mm}^{2}$) map obtained using the fit from the
analytical model (left) and the time difference method (right) described in
the text. Figure 22: The distribution of the charge in the leading pad (left)
and the charge fraction in the second to leading pad (right) versus the column
for horizontal beam tracks.

##  10 Conclusions

We measured the performance of the ERAM prototype with beam particles at DESY.
We studied both spatial and $dE/dx$ resolution as a function of the angle of
the track with respect to the ERAM plane. We also characterized charge
spreading and produced a RC map of the prototype.
Spatial resolution better than 600 µm is obtained for all the angles using a
dedicated clustering algorithm which is adapted to the track angle. Energy
resolution better than 9% is obtained for all the angles. We expect ionization
energy loss resolution to be better than 7% for tracks crossing two ERAMs.
Such performances fully satisfies the requirements for the upgrade of the
ND280 TPC.

## Acknowledgements

The measurements leading to these results have been performed at the Test Beam
Facility at DESY Hamburg (Germany), a member of the Helmholtz Association. The
authors would like to thank the technical team at the DESY II accelerator and
test beam facility for the smooth operation of the test beam and the support
during the test beam campaign.
We acknowledge the support of CEA and CNRS/IN2P3, France; DFG, Germany; INFN,
Italy; National Science Centre (NCN) and Ministry of Science and Higher
Education (Grant No. DIR/WK/2017/05), Poland; MINECO and ERDF funds, Spain.
In addition, participation of individual researchers and institutions has been
further supported by H2020 Grant No. RISE-GA822070-JENNIFER2 2020, MSCA-
COFUND-2016 No.754496, ANR-19-CE31-0001, RFBR grants #19-32-90100, the
Ministry of Science and Higher Education of Russia (contract #075-15-2020-778)
the Spanish Ministerio de Economía y Competitividad (SEIDI - MINECO) under
Grants No. PID2019-107564GB-I00 and SEV-2016-0588. IFAE is partially funded by
the CERCA program of the Generalitat de Catalunya.

## References

  * [1] K. Abe, et al., The T2K Experiment, Nucl. Instrum. Meth. A659 (2011) 106–135.  arXiv:1106.1238, doi:10.1016/j.nima.2011.06.067.
  * [2] K. Abe, et al., Indication of Electron Neutrino Appearance from an Accelerator-produced Off-axis Muon Neutrino Beam, Phys. Rev. Lett. 107 (2011) 041801.  arXiv:1106.2822, doi:10.1103/PhysRevLett.107.041801.
  * [3] K. Abe, et al., Evidence of Electron Neutrino Appearance in a Muon Neutrino Beam, Phys. Rev. D 88 (3) (2013) 032002.  arXiv:1304.0841, doi:10.1103/PhysRevD.88.032002.
  * [4] K. Abe, et al., Observation of Electron Neutrino Appearance in a Muon Neutrino Beam, Phys. Rev. Lett. 112 (2014) 061802.  arXiv:1311.4750, doi:10.1103/PhysRevLett.112.061802.
  * [5] K. Abe, et al., Constraint on the matter–antimatter symmetry-violating phase in neutrino oscillations, Nature 580 (7803) (2020) 339–344, [Erratum: Nature 583, E16 (2020)].  arXiv:1910.03887, doi:10.1038/s41586-020-2177-0.
  * [6] N. Abgrall, et al., Time Projection Chambers for the T2K Near Detectors, Nucl. Instrum. Meth. A637 (2011) 25–46.  arXiv:1012.0865, doi:10.1016/j.nima.2011.02.036.
  * [7] K. Abe, et al., Measurement of the intrinsic electron neutrino component in the T2K neutrino beam with the ND280 detector, Phys. Rev. D 89 (2014) 092003\.  arXiv:1403.2552, doi:10.1103/PhysRevD.89.092003.
  * [8] K. Abe, et al., Measurement of the charged-current electron (anti-)neutrino inclusive cross-sections at the T2K off-axis near detector ND280, JHEP 10 (2020) 114.  arXiv:2002.11986, doi:10.1007/JHEP10(2020)114.
  * [9] A. Blondel, M. Bogomilov, S. Bordoni, F. Cadoux, D. Douqa, K. Dugas, T. Ekelof, Y. Favre, S. Fedotov, K. Fransson, R. Fujita, E. Gramstad, A. Ichikawa, S. Ilieva, K. Iwamoto, C. Jesús-Valls, C. Jung, S. Kasetti, M. Khabibullin, A. Khotjantsev, A. Korzenev, A. Kostin, Y. Kudenko, T. Kutter, T. Lux, L. Maret, T. Matsubara, A. Mefodiev, A. Minamino, O. Mineev, G. Mitev, M. Nessi, L. Nicola, E. Noah, S. Parsa, G. Petkov, F. Sanchez, D. Sgalaberna, W. Shorrock, K. Skwarczynski, S. Suvorov, A. Teklu, R. Tsenov, Y. Uchida, G. Vankova-Kirilova, N. Yershov, M. Yokoyama, J. Zalipska, Y. Zou, W. Zurek, The SuperFGD prototype charged particle beam tests, Journal of Instrumentation 15 (12) (2020) P12003–P12003.  doi:10.1088/1748-0221/15/12/p12003.
URL https://doi.org/10.1088/1748-0221/15/12/p12003

  * [10] D. Attié, et al., Performances of a resistive Micromegas module for the Time Projection Chambers of the T2K Near Detector upgrade, Nucl. Instrum. Meth. A 957 (2020) 163286.  arXiv:1907.07060, doi:10.1016/j.nima.2019.163286.
  * [11] K. Abe, et al., T2K ND280 Upgrade - Technical Design Report (2019).  arXiv:1901.03750.
  * [12] I. Giomataris, R. De Oliveira, S. Andriamonje, S. Aune, G. Charpak, P. Colas, A. Giganon, P. Rebourgeard, P. Salin, Micromegas in a bulk, Nucl. Instrum. Meth. A560 (2006) 405–408.  arXiv:physics/0501003, doi:10.1016/j.nima.2005.12.222.
  * [13] M. S. Dixit, J. Dubeau, J. P. Martin, K. Sachs, Position sensing from charge dispersion in micropattern gas detectors with a resistive anode, Nucl. Instrum. Meth. A518 (2004) 721–727.  arXiv:physics/0307152, doi:10.1016/j.nima.2003.09.051.
  * [14] R. Diener, et al., The DESY II Test Beam Facility, Nucl. Instrum. Meth. A 922 (2019) 265–286.  arXiv:1807.09328, doi:10.1016/j.nima.2018.11.133.
  * [15] M. Ester, H.-P. Kriegel, J. Sander, X. Xu, A density-based algorithm for discovering clusters in large spatial databases with noise, in: Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, AAAI Press, 1996, pp. 226–231.
  * [16] R. Gluckstern, Uncertainties in track momentum and direction, due to multiple scattering and measurement errors, Nuclear Instruments and Methods 24 (1963) 381–389.  doi:https://doi.org/10.1016/0029-554X(63)90347-1.
URL https://www.sciencedirect.com/science/article/pii/0029554X63903471

  * [17] D. Attie, Beam tests of Micromegas LC-TPC large prototype, JINST 6 (2011) C01007.  doi:10.1088/1748-0221/6/01/C01007.
  * [18] K. Boudjemline, M. S. Dixit, J. P. Martin, K. Sachs, Spatial resolution of a GEM readout TPC using the charge dispersion signal, Nucl. Instrum. Meth. A574 (2007) 22–27.  arXiv:physics/0610232, doi:10.1016/j.nima.2007.01.017.
  * [19] M. S. Dixit, J. Dubeau, J. P. Martin, K. Sachs, Position sensing from charge dispersion in micropattern gas detectors with a resistive anode, Nucl. Instrum. Meth. A518 (2004) 721–727.  arXiv:physics/0307152, doi:10.1016/j.nima.2003.09.051.