$\displaystyle=\int_{\operatorname{supp}(Y)}\frac{m_{2}(y)^{2}}{m_{0}(y)}\mathrm{d}y-2\int_{\operatorname{supp}(Y)}m_{2}(y)\mathrm{d}y+\int_{\operatorname{supp}(Y)}m_{0}(y)\mathrm{d}y$
$\displaystyle=\int_{\operatorname{supp}(Y)}\frac{m_{2}(y)^{2}}{m_{0}(y)}\mathrm{d}y-1,$

since

$\displaystyle\begin{split}\int_{\operatorname{supp}(Y)}m_{2}(y)\mathrm{d}y&=\int_{\operatorname{supp}(Y)}\mathbb{E}\left[p(y|G)G^{2}\right]\mathrm{d}y=\mathbb{E}\left[G^{2}\right]=1,\\\
\int_{\operatorname{supp}(Y)}m_{0}(y)\mathrm{d}y&=\int_{\operatorname{supp}(Y)}\mathbb{E}\left[p(y|G)\right]\mathrm{d}y=1.\end{split}$
(D.2)

To show Equation D.1, it then suffices to show
$\int\frac{m_{2}^{2}}{m_{0}}\leq 3$. Using the Cauchy–Schwarz inequality, we
can bound $m_{2}^{2}$ as follows,

$\displaystyle m_{2}(y)^{2}$
$\displaystyle=\left(\int_{\mathbb{R}}f(g)p(y|g)g^{2}\mathrm{d}g\right)^{2}$
$\displaystyle=\left(\int_{\mathbb{R}}\sqrt{f(g)p(y|g)}g^{2}\cdot\sqrt{f(g)p(y|g)}\mathrm{d}g\right)^{2}$
$\displaystyle\leq\left(\int_{\mathbb{R}}f(g)p(y|g)g^{4}\mathrm{d}g\right)\cdot\left(\int_{\mathbb{R}}f(g)p(y|g)\mathrm{d}g\right)$
$\displaystyle=\left(\int_{\mathbb{R}}f(g)p(y|g)g^{4}\mathrm{d}g\right)\cdot
m_{0}(y),$

where $f(g)$ denotes the standard Gaussian density function. The desired bound
then follows:

$\displaystyle\int_{\operatorname{supp}(Y)}\frac{m_{2}(y)^{2}}{m_{0}(y)}\mathrm{d}y$
$\displaystyle\leq\int_{\operatorname{supp}(Y)}\int_{\mathbb{R}}f(g)p(y|g)g^{4}\mathrm{d}g\mathrm{d}y$
$\displaystyle=\int_{\mathbb{R}}\left[f(g)\left(\int_{\operatorname{supp}(Y)}p(y|g)\mathrm{d}y\right)g^{4}\right]\mathrm{d}g$
$\displaystyle=\mathbb{E}\left[G^{4}\right]$ $\displaystyle=3.$

##  Appendix E Evaluation of spectral threshold and overlap

Recall that the optimal threshold for weakly recovering $x_{1}^{*}$ and
$x_{2}^{*}$ using the spectral method are given by

$\displaystyle\delta^{*}_{1}$
$\displaystyle=\left(\alpha^{2}\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y\right)^{-1},$
(E.1a) $\displaystyle\delta^{*}_{2}$
$\displaystyle=\left((1-\alpha)^{2}\int_{\operatorname{supp}(Y)}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y\right)^{-1},$
(E.1b)

respectively. In what follows, we will evaluate the optimal spectral threshold
and the optimal overlap for two well-studied models: _mixed linear regression_
and _mixed phase retrieval_.

###  E.1 Mixed linear regression (proof of Corollary 3.8)

Consider the following _mixed linear regression_ problem, a special case of
mixed GLM,

$\displaystyle y_{i}$ $\displaystyle=\left\langle
a_{i},\eta_{i}x_{1}^{*}+(1-\eta_{i})x_{2}^{*}\right\rangle+\varepsilon_{i}$

for $1\leq i\leq n$. Here the signals
$x_{1}^{*},x_{2}^{*}\in{\mathbb{R}}^{d}$, the row of the design matrix
$a_{i}\in{\mathbb{R}}^{d}$, the mixing variable $\eta_{i}\in\\{0,1\\}$ have
been defined in Section 2.2 and are subject to Items (A1), (A2), (A4) and
(A5). We further assume that
$\varepsilon_{i}\overset{\mathrm{i.i.d.}}{\sim}{\mathcal{N}}(0,\sigma^{2})$
for some $\sigma>0$.

Let $g_{i}=\left\langle
a_{i},\eta_{i}x_{1}^{*}+(1-\eta_{i})x_{2}^{*}\right\rangle$ and let
$p(y_{i}|g_{i})$ denote the conditional density of $y_{i}$ given $g_{i}$. Due
to the Gaussianity of $\varepsilon_{i}$, we have

$\displaystyle p(y|g)$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(y-g)^{2}}{2\sigma^{2}}\right).$
(E.2)

To evaluate the integral in Equation E.1, let us first calculate
$\mathbb{E}\left[p(y|G)\right]$ and $\mathbb{E}\left[p(y|G)(G^{2}-1)\right]$.
We have

$\displaystyle\mathbb{E}\left[p(y|G)\right]$
$\displaystyle=\int_{\mathbb{R}}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{g^{2}}{2}\right)\cdot\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(y-g)^{2}}{2\sigma^{2}}\right)\mathrm{d}g$
$\displaystyle=(f_{G}*f_{\varepsilon})(y)$
$\displaystyle=f_{G+\varepsilon}(y)$
$\displaystyle=\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right),$
(E.3)

and

$\displaystyle\mathbb{E}\left[p(y|G)G^{2}\right]$
$\displaystyle=\int_{\mathbb{R}}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{g^{2}}{2}\right)\cdot\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(y-g)^{2}}{2\sigma^{2}}\right)g^{2}\mathrm{d}g$
$\displaystyle=\frac{1}{2\pi\sigma}\int_{\mathbb{R}}\exp\left(-\frac{1+\sigma^{2}}{2\sigma^{2}}g^{2}+\frac{y}{\sigma^{2}}g-\frac{y^{2}}{2\sigma^{2}}\right)g^{2}\mathrm{d}g$
$\displaystyle=\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\cdot\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}.$
(E.4)

Therefore,

$\displaystyle\mathbb{E}\left[p(y|G)(G^{2}-1)\right]$
$\displaystyle=\mathbb{E}\left[p(y|G)G^{2}\right]-\mathbb{E}\left[p(y|G)\right]$
$\displaystyle=\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\cdot\frac{y^{2}-(1+\sigma^{2})}{(1+\sigma^{2})^{2}},$
(E.5)

and

$\displaystyle\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}$
$\displaystyle=\left(\frac{y^{2}-(1+\sigma^{2})}{(1+\sigma^{2})^{2}}\right)^{2}\cdot\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right).$

We then see that the evaluation of integral in Equation E.1 can be reduced to
computing Gaussian moments:

$\displaystyle\int_{\mathbb{R}}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=\int_{\mathbb{R}}\left(\frac{y^{2}-(1+\sigma^{2})}{(1+\sigma^{2})^{2}}\right)^{2}\cdot\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\mathrm{d}y$
$\displaystyle=\mathop{\mathbb{E}}_{Y\sim{\mathcal{N}}(0,1+\sigma^{2})}\left[\left(\frac{Y^{2}-(1+\sigma^{2})}{(1+\sigma^{2})^{2}}\right)^{2}\right]$
$\displaystyle=\frac{1}{(1+\sigma^{2})^{4}}\left(\mathop{\mathbb{E}}_{Y\sim{\mathcal{N}}(0,1+\sigma^{2})}\left[Y^{4}\right]-2(1+\sigma^{2})\mathop{\mathbb{E}}_{Y\sim{\mathcal{N}}(0,1+\sigma^{2})}\left[Y^{2}\right]+(1+\sigma^{2})^{2}\right)$
$\displaystyle=\frac{1}{(1+\sigma^{2})^{4}}\left(3(1+\sigma^{2})^{2}-2(1+\sigma^{2})(1+\sigma^{2})+(1+2\sigma^{2})^{2}\right)$
$\displaystyle=\frac{2}{(1+\sigma^{2})^{2}}.$

This allows us to conclude that

$\displaystyle\delta^{*}_{1}$
$\displaystyle=\frac{(1+\sigma^{2})^{2}}{2\alpha^{2}},\quad\delta^{*}_{2}=\frac{(1+\sigma^{2})^{2}}{2(1-\alpha)^{2}}.$

Furthermore, according to Equation 3.25, the corresponding optimal
preprocessing functions specialize to

$\displaystyle{\mathcal{T}}_{1}^{*}(y)$
$\displaystyle=1-\frac{1}{\alpha\cdot\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+(1-\alpha)},\quad{\mathcal{T}}_{2}^{*}(y)=1-\frac{1}{(1-\alpha)\cdot\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+\alpha}.$

By Equations 3.22, 3.23 and 3.24, the above optimal preprocessing functions
yield the following optimal overlaps

$\displaystyle\frac{1}{\sqrt{\beta_{1}^{*}(\delta,\alpha,\sigma)+\alpha}},\quad\frac{1}{\sqrt{\beta_{2}^{*}(\delta,\alpha,\sigma)+(1-\alpha)}},$

where
$\beta_{1}^{*}(\delta,\alpha,\sigma),\beta_{2}^{*}(\delta,\alpha,\sigma)$ are
the unique positive solutions to the following equations

$\displaystyle(\beta_{1}^{*}(\delta,\alpha,\sigma)-(1-\alpha))\int_{\mathbb{R}}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\alpha\mathbb{E}\left[p(y|G)G^{2}\right]+\beta_{1}^{*}(\delta,\alpha,\sigma)\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=\frac{1}{\alpha^{2}\delta},$ (E.6)
$\displaystyle(\beta_{2}^{*}(\delta,\alpha,\sigma)-\alpha)\int_{\mathbb{R}}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{(1-\alpha)\mathbb{E}\left[p(y|G)G^{2}\right]+\beta_{2}^{*}(\delta,\alpha,\sigma)\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=\frac{1}{(1-\alpha)^{2}\delta},$ (E.7)

respectively. Using the formulas we just derived (Equations E.3, E.4 and E.5),
the fixed point equations (Equations E.6 and E.7) above can be simplified to

$\displaystyle(\beta_{1}^{*}(\delta,\alpha,\sigma)-(1-\alpha))\mathop{\mathbb{E}}_{G\sim{\mathcal{N}}(0,1)}\left[\frac{\left(\frac{G^{2}-1}{1+\sigma^{2}}\right)^{2}}{\alpha\cdot\frac{G^{2}+\sigma^{2}}{1+\sigma^{2}}+\beta_{1}^{*}(\delta,\alpha,\sigma)}\right]$
$\displaystyle=\frac{1}{\alpha^{2}\delta},$
$\displaystyle(\beta_{2}^{*}(\delta,\alpha,\sigma)-\alpha)\mathop{\mathbb{E}}_{G\sim{\mathcal{N}}(0,1)}\left[\frac{\left(\frac{G^{2}-1}{1+\sigma^{2}}\right)^{2}}{(1-\alpha)\cdot\frac{G^{2}+\sigma^{2}}{1+\sigma^{2}}+\beta_{2}^{*}(\delta,\alpha,\sigma)}\right]$
$\displaystyle=\frac{1}{(1-\alpha)^{2}\delta},$

since the integral in Equations E.6 and E.7 can be written as

$\displaystyle\int_{\mathbb{R}}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\alpha\mathbb{E}\left[p(y|G)G^{2}\right]+\beta\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=\int_{\mathbb{R}}\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\frac{\left(\frac{y^{2}-(1+\sigma^{2})}{(1+\sigma^{2})^{2}}\right)^{2}}{\alpha\cdot\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+\beta}\mathrm{d}y$
$\displaystyle=\mathop{\mathbb{E}}_{Y\sim{\mathcal{N}}(0,1+\sigma^{2})}\left[\frac{\left(\frac{Y^{2}-(1+\sigma^{2})}{(1+\sigma^{2})^{2}}\right)^{2}}{\alpha\cdot\frac{Y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+\beta}\right]$
$\displaystyle=\mathop{\mathbb{E}}_{G\sim{\mathcal{N}}(0,1)}\left[\frac{\left(\frac{G^{2}-1}{1+\sigma^{2}}\right)^{2}}{\alpha\cdot\frac{G^{2}+\sigma^{2}}{1+\sigma^{2}}+\beta}\right]$
$\displaystyle=-\frac{\alpha+\beta}{\alpha^{2}}+\left(\frac{\alpha+\beta}{\alpha}\right)^{2}\sqrt{\frac{\pi(1+\sigma^{2})^{2}}{2\alpha(\sigma^{2}\alpha+(1+\sigma^{2})\beta)}}\exp\left(\frac{\sigma^{2}\alpha+(1+\sigma^{2})\beta}{2\alpha}\right)\mathrm{erfc}\left(\sqrt{\frac{\sigma^{2}\alpha+(1+\sigma^{2})\beta}{2\alpha}}\right).$

As a final remark, we observe that as $\sigma\to 0$, $\delta^{*}_{1}$ and
$\delta^{*}_{2}$ scale as

$\displaystyle\delta^{*}_{1}$
$\displaystyle=\frac{1+2\sigma^{2}}{2\alpha^{2}}+{\mathcal{O}}(\sigma^{4}),\quad\delta^{*}_{2}=\frac{1+2\sigma^{2}}{2(1-\alpha)^{2}}+{\mathcal{O}}(\sigma^{4}),$
(E.8)

respectively.

###  E.2 Mixed phase retrieval (proof of Corollary 3.9)

Consider another special case of mixed GLM: the _mixed phase retrieval_
problem defined below

$\displaystyle y_{i}$ $\displaystyle=\left|\left\langle
a_{i},\eta_{i}x_{1}^{*}+(1-\eta_{i})x_{2}^{*}\right\rangle\right|+\varepsilon_{i},$

for each $1\leq i\leq n$. Again, we assume
$\varepsilon_{i}\overset{\mathrm{i.i.d.}}{\sim}{\mathcal{N}}(0,\sigma^{2})$
for some $\sigma>0$.

The conditional density $p(y|g)$ for this model is given by

$\displaystyle p(y|g)$
$\displaystyle\coloneqq\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(y-|g|)^{2}}{2\sigma^{2}}\right).$

Analogously, let us start by computing $\mathbb{E}\left[p(y|G)\right]$ and
$\mathbb{E}\left[p(|G)G^{2}\right]$:

$\displaystyle\mathbb{E}\left[p(y|G)\right]$
$\displaystyle=\int_{\mathbb{R}}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{g^{2}}{2}\right)\cdot\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(y-|g|)^{2}}{2\sigma^{2}}\right)\mathrm{d}g$
$\displaystyle=2\int_{0}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{g^{2}}{2}\right)\cdot\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(y-g)^{2}}{2\sigma^{2}}\right)\mathrm{d}g$
$\displaystyle=\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)\right],$

and

$\displaystyle\mathbb{E}\left[p(y|G)G^{2}\right]$
$\displaystyle=\int_{\mathbb{R}}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{g^{2}}{2}\right)\cdot\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{(y-|g|)^{2}}{2\sigma^{2}}\right)g^{2}\mathrm{d}g$
$\displaystyle=\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)\right]\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}$
$\displaystyle\phantom{=}+\frac{\sigma
y}{\pi(1+\sigma^{2})^{2}}\exp\left(-\frac{y^{2}}{2\sigma^{2}}\right).$

Thus,

$\displaystyle\mathbb{E}\left[p(y|G)(G^{2}-1)\right]$
$\displaystyle=\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)\right]\frac{y^{2}-(1+\sigma^{2})}{(1+\sigma^{2})^{2}}$
$\displaystyle\phantom{=}+\frac{\sigma
y}{\pi(1+\sigma^{2})^{2}}\exp\left(-\frac{y^{2}}{2\sigma^{2}}\right),$

and

$\displaystyle\phantom{=}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}$
$\displaystyle=\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)\right]\cdot\left(\frac{y^{2}-(1+\sigma^{2})}{(1+\sigma^{2})^{2}}\right)^{2}$
(E.9)
$\displaystyle\phantom{=}+2\cdot\frac{y^{2}-(1+\sigma^{2})}{(1+\sigma^{2})^{2}}\cdot\frac{\sigma
y}{\pi(1+\sigma^{2})^{2}}\exp\left(-\frac{y^{2}}{2\sigma^{2}}\right)$ (E.10)
$\displaystyle\phantom{=}+\frac{\left[\frac{\sigma
y}{\pi(1+\sigma^{2})^{2}}\exp\left(-\frac{y^{2}}{2\sigma^{2}}\right)\right]^{2}}{\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)\right]}.$
(E.11)

Let

$\displaystyle\texttt{I}(y)$
$\displaystyle\coloneqq\text{\lx@cref{creftypecap~refnum}{eqn:pr-eval-
term1}},\quad\texttt{II}(y)\coloneqq\text{\lx@cref{creftypecap~refnum}{eqn:pr-
eval-
term2}},\quad\texttt{III}(y)\coloneqq\text{\lx@cref{creftypecap~refnum}{eqn:pr-
eval-term3}}.$

We will evaluate the integrals of
$\texttt{I}(y),\texttt{II}(y),\texttt{III}(y)$ term by term. The first two
terms admit closed-form integrals:

$\displaystyle\int_{\mathbb{R}}\texttt{I}(y)\mathrm{d}y$
$\displaystyle=\frac{2}{(1+\sigma^{2})^{2}},\quad\int_{\mathbb{R}}\texttt{II}(y)\mathrm{d}y=0.$

As for the third term, we have

$\displaystyle\int_{\mathbb{R}}\texttt{III}(y)\mathrm{d}y$
$\displaystyle=\frac{\sqrt{2}\sigma^{2}}{\pi^{3/2}(1+\sigma^{2})^{7/2}}\int_{\mathbb{R}}\frac{\exp\left(-\frac{2+\sigma^{2}}{2\sigma^{2}(1+\sigma^{2})}y^{2}\right)y^{2}}{1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)}\mathrm{d}y$
$\displaystyle=\frac{4\sigma^{5}}{\pi^{3/2}(1+\sigma^{2})^{2}}\int_{\mathbb{R}}\frac{\exp\left(-(2+\sigma^{2})z^{2}\right)z^{2}}{1+\mathrm{erf}(z)}\mathrm{d}z.$

In the last equality, we apply the change of variable
$z=\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}$.

Finally, defining

$\displaystyle h(s)$
$\displaystyle\coloneqq\int_{\mathbb{R}}\frac{\exp\left(-(2+s)z^{2}\right)z^{2}}{1+\mathrm{erf}(z)}\mathrm{d}z,$

we obtain

$\displaystyle\int_{\mathbb{R}}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=\int_{\mathbb{R}}\texttt{I}(y)\mathrm{d}y+\int_{\mathbb{R}}\texttt{II}(y)\mathrm{d}y+\int_{\mathbb{R}}\texttt{III}(y)\mathrm{d}y$
$\displaystyle=\frac{2}{(1+\sigma^{2})^{2}}+\frac{4\sigma^{5}h(\sigma^{2})}{\pi^{3/2}(1+\sigma^{2})^{2}}.$

We conclude that

$\displaystyle\delta^{*}_{1}$
$\displaystyle=\frac{1}{\alpha^{2}}\left(\frac{2}{(1+\sigma^{2})^{2}}+\frac{4\sigma^{5}h(\sigma^{2})}{\pi^{3/2}(1+\sigma^{2})^{2}}\right)^{-1},\quad\delta^{*}_{2}=\frac{1}{(1-\alpha)^{2}}\left(\frac{2}{(1+\sigma^{2})^{2}}+\frac{4\sigma^{5}h(\sigma^{2})}{\pi^{3/2}(1+\sigma^{2})^{2}}\right)^{-1}.$

Furthermore, according to Equation 3.25, the corresponding optimal
preprocessing functions specialize to

$\displaystyle{\mathcal{T}}_{1}^{*}(y)$
$\displaystyle=1-\frac{1}{\alpha\Delta(y)+(1-\alpha)},\quad{\mathcal{T}}_{2}^{*}(y)=1-\frac{1}{(1-\alpha)\Delta(y)+\alpha},$

where

$\displaystyle\Delta(y)$
$\displaystyle\coloneqq\frac{\mathbb{E}\left[p(y|G)G^{2}\right]}{\mathbb{E}\left[p(y|G)\right]}$
$\displaystyle=\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+\frac{\frac{\sigma
y}{\pi(1+\sigma^{2})^{2}}\exp\left(-\frac{y^{2}}{2\sigma^{2}}\right)}{\frac{1}{\sqrt{2\pi(1+\sigma^{2})}}\exp\left(-\frac{y^{2}}{2(1+\sigma^{2})}\right)\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)\right]}$
$\displaystyle=\frac{y^{2}+\sigma^{2}+\sigma^{4}}{(1+\sigma^{2})^{2}}+\sqrt{\frac{2}{\pi}}\cdot\frac{\sigma
y\exp\left(-\frac{y^{2}}{2\sigma^{2}(1+\sigma^{2})}\right)}{(1+\sigma^{2})^{3/2}\left[1+\mathrm{erf}\left(\frac{y}{\sqrt{2\sigma^{2}(1+\sigma^{2})}}\right)\right]}.$

We observe that as $\sigma\downarrow 0$, the asymptotic behaviours of
$h(\sigma)$ and
$\int_{\mathbb{R}}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
are easy to understand. Indeed,

$\displaystyle h(\sigma^{2})$
$\displaystyle=\int_{\mathbb{R}}\frac{e^{-2z^{2}}z^{2}}{1+\mathrm{erf}(z)}\cdot
e^{-\sigma^{2}z^{2}}\mathrm{d}z\xrightarrow{\sigma\to
0}\int_{\mathbb{R}}\frac{e^{-2z^{2}}z^{2}}{1+\mathrm{erf}(z)}\mathrm{d}z\approx
1.22564\in(0,\infty).$

Hence

$\displaystyle\int_{\mathbb{R}}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=2-4\sigma^{2}+{\mathcal{O}}(\sigma^{4}),$

and

$\displaystyle\delta^{*}_{1}$
$\displaystyle=\frac{1+2\sigma^{2}}{2\alpha^{2}}+{\mathcal{O}}(\sigma^{4}),\quad\delta^{*}_{2}=\frac{1+2\sigma^{2}}{2(1-\alpha)^{2}}+{\mathcal{O}}(\sigma^{4}).$
(E.12)

If $\sigma=0$, $\Delta(y)=y^{2}$, so the optimal preprocessing functions take
the form

$\displaystyle{\mathcal{T}}_{1}^{*}(y)$ $\displaystyle=1-\frac{1}{\alpha
y^{2}+(1-\alpha)},\quad{\mathcal{T}}_{2}^{*}(y)=1-\frac{1}{(1-\alpha)y^{2}+\alpha}.$

Consequently, the optimal overlaps are given by

$\displaystyle\frac{1}{\sqrt{\beta_{1}^{*}(\delta,\alpha)+\alpha}},\quad\frac{1}{\sqrt{\beta_{2}^{*}(\delta,\alpha)+(1-\alpha)}},$

where $\beta_{1}^{*}(\delta,\alpha),\beta_{2}^{*}(\delta,\alpha)$ are the
unique positive solutions to the fixed point equations (3.23). The integral in
Equation 3.23 can be evaluated explicitly:

$\displaystyle\int_{0}^{\infty}\frac{\mathbb{E}\left[p(y|G)(G^{2}-1)\right]^{2}}{\alpha\mathbb{E}\left[p(y|G)G^{2}\right]+\beta\mathbb{E}\left[p(y|G)\right]}\mathrm{d}y$
$\displaystyle=\int_{\mathbb{R}}\sqrt{\frac{2}{\pi}}\exp\left(-\frac{y^{2}}{2}\right)\frac{(y^{2}-1)^{2}}{\alpha
y^{2}+\beta}\mathrm{d}y$
$\displaystyle=-\frac{\alpha+\beta}{\alpha^{2}}+\sqrt{\frac{\pi}{2\alpha\beta}}\left(\frac{\alpha+\beta}{\alpha}\right)^{2}\exp\left(\frac{\beta}{2\alpha}\right)\mathrm{erfc}\left(\sqrt{\frac{\beta}{2\alpha}}\right).$

Therefore in the case of mixed noiseless phase retrieval, the fixed point
equations (3.23) can be written explicitly as:

$(\beta_{1}^{*}(\delta,\alpha)-(1-\alpha))\left[-\frac{\alpha+\beta_{1}^{*}(\delta,\alpha)}{\alpha^{2}}+\sqrt{\frac{\pi}{2\alpha\beta_{1}^{*}(\delta,\alpha)}}\left(\frac{\alpha+\beta_{1}^{*}(\delta,\alpha)}{\alpha}\right)^{2}\right.\\\
\left.\exp\left(\frac{\beta_{1}^{*}(\delta,\alpha)}{2\alpha}\right)\mathrm{erfc}\left(\sqrt{\frac{\beta_{1}^{*}(\delta,\alpha)}{2\alpha}}\right)\right]=\frac{1}{\alpha^{2}\delta},$
(E.13)

and

$(\beta_{2}^{*}(\delta,\alpha)-\alpha)\left[-\frac{(1-\alpha)+\beta_{2}^{*}(\delta,\alpha)}{(1-\alpha)^{2}}+\sqrt{\frac{\pi}{2(1-\alpha)\beta_{2}^{*}(\delta,\alpha)}}\left(\frac{(1-\alpha)+\beta_{2}^{*}(\delta,\alpha)}{1-\alpha}\right)^{2}\right.\\\
\left.\exp\left(\frac{\beta_{2}^{*}(\delta,\alpha)}{2(1-\alpha)}\right)\mathrm{erfc}\left(\sqrt{\frac{\beta_{2}^{*}(\delta,\alpha)}{2(1-\alpha)}}\right)\right]=\frac{1}{(1-\alpha)^{2}\delta}.$
(E.14)

##  Appendix F Random matrix and free probability theory background

For a finite positive Borel measure (including probability measures) $\mu$ on
${\mathbb{R}}$, the _Stieltjes transform_
$s_{\mu}\colon{\mathbb{C}}\setminus\operatorname{supp}(\mu)\to{\mathbb{C}}$ of
$\mu$ is defined as

$\displaystyle s_{\mu}(z)$
$\displaystyle\coloneqq\int_{\mathbb{R}}\frac{1}{x-z}\mathrm{d}\mu(x).$ (F.1)

The _$R$ -transform_ [Voi86] of $\mu$ is the convergent power series defined
by

$\displaystyle R_{\mu}(z)$ $\displaystyle\coloneqq
s_{\mu}^{-1}(-z)-\frac{1}{z},$ (F.2)

where $s_{\mu}^{-1}$ is the functional inverse of $s_{\mu}$. The _free
additive convolution_ , denoted by $\mu\boxplus\nu$, of two compactly
supported probability measures $\mu$ and $\nu$ is a compactly supported
probability measure characterized by the following relation:

$\displaystyle R_{\mu\boxplus\nu}$ $\displaystyle=R_{\mu}+R_{\nu}.$ (F.3)

For a symmetric matrix $M\in{\mathbb{R}}^{p\times p}$, the _empirical spectral
distribution_ of $M$ is defined as

$\displaystyle\mu_{M}$
$\displaystyle\coloneqq\frac{1}{p}\sum_{i=1}^{p}\delta_{\lambda_{i}(M)}.$

Suppose $A,B\in{\mathbb{R}}^{p\times p}$ are two (sequences of) deterministic
matrices whose empirical spectral distributions converge weakly to two
compactly supported distributions $\mu$ and $\nu$, respectively. Let $O$ be a
random orthogonal matrix distributed according to the Haar distribution on the
orthogonal group ${\mathbb{O}}(p)$. Then $A$ and $OBO^{\top}$ are almost
surely _asymptotically free_. The same conclusion continues to hold if $A$ and
$B$ are random matrices independent of $O$. As a result [Voi91, Spe93], the
empirical spectral distribution of $D\coloneqq A+OBO^{\top}$ converges weakly
to $\mu\boxplus\nu$ almost surely. The free additive convolution
$\mu\boxplus\nu$ is again a compactly supported distribution on
${\mathbb{R}}$. We caution that the weak convergence of the empirical spectral
distribution of $D$ does _not_ imply that all eigenvalues of $D$ are close to
$\operatorname{supp}(\mu\boxplus\nu)$. The spectrum of $D$ may contain
outliers, though they must not affect the limiting spectral distribution
$\mu\boxplus\nu$.

Given a pair of Borel probability measures $\mu$ and $\nu$ on ${\mathbb{R}}$,
there exists a unique pair of analytic functions
$w_{1},w_{2}\colon{\mathbb{H}}\to{\mathbb{H}}$ satisfying the following three
properties, where ${\mathbb{H}}\coloneqq\\{z\in{\mathbb{C}}:\Im(z)>0\\}$
denotes upper complex half-plane.

  1. 1.

For $j\in\\{1,2\\}$,

$\displaystyle\lim_{z\to\infty}\frac{w_{j}(\iota z)}{\iota z}$
$\displaystyle=1,$

where $\iota$ denotes the imaginary unit.

  2. 2.

For every $z\in{\mathbb{H}}$,

$\displaystyle s_{\mu\boxplus\nu}(z)$
$\displaystyle=s_{\mu}(w_{1}(z))=s_{\nu}(w_{2}(z))=\frac{1}{z-(w_{1}(z)+w_{2}(z))}.$
(F.4)

  3. 3.

For any $z\in{\mathbb{H}}\cup{\mathbb{R}}$ such that $w_{1}$ is analytic at
$z$, $w_{1}(z)$ is the attracting fixed point of the self-map of
${\mathbb{H}}$ defined by

$\displaystyle w$ $\displaystyle\mapsto
F_{\nu}(F_{\mu}(w)-w+z)-(F_{\mu}(w)-w)$

where

$\displaystyle F_{\mu}(z)$ $\displaystyle\coloneqq-\frac{1}{s_{\mu}(z)},\quad
F_{\nu}(z)\coloneqq-\frac{1}{s_{\nu}(z)}.$

A similar statement, with $\mu$ and $\nu$ interchanged, holds for $w_{2}$.

The functions $w_{1},w_{2}$ defined through the above functional equations are
called _subordination functions_ [Voi93, Bia98, BB07]. We emphasize that each
of the subordination functions $w_{1},w_{2}$ depends on _both_ distributions
$\mu,\nu$. Though $w_{1},w_{2}$ are defined on the upper half-plane, they can
be continuously extended to the real line.

######  Lemma F.1 ([BBCF17, Lemma 3.1]).

Let $\mu$ and $\nu$ be two compactly supported Borel probability measures on
${\mathbb{R}}$ neither of which is a point mass. Then the subordination
functions $w_{1},w_{2}\colon{\mathbb{H}}\to{\mathbb{H}}$ associated with
$\mu,\nu$ can be extended to ${\mathbb{H}}\cup{\mathbb{R}}\cup\\{\infty\\}$
with the following properties being preserved.

  1. 1.

$w_{1},w_{2}\colon{\mathbb{H}}\cup{\mathbb{R}}\cup\\{\infty\\}\to{\mathbb{H}}\cup{\mathbb{R}}\cup\\{\infty\\}$
are continuous.

  2. 2.

If $x\in{\mathbb{R}}\setminus\operatorname{supp}(\mu\boxplus\nu)$, then
$w_{1},w_{2}$ continue meromorphically to a neighbourhood of $x$:

$\displaystyle w_{1}(x)$
$\displaystyle=F_{\nu}(w_{2}(x))-w_{2}(x)+x\in({\mathbb{R}}\cup\\{\infty\\})\setminus\operatorname{supp}(\mu),$
$\displaystyle w_{2}(x)$
$\displaystyle=F_{\mu}(w_{1}(x))-w_{1}(x)+x\in({\mathbb{R}}\cup\\{\infty\\})\setminus\operatorname{supp}(\nu).$

If $w_{1}(x)=\infty$, then

$\displaystyle w_{2}(x)$
$\displaystyle=x-\int_{\mathbb{R}}z\mathrm{d}\mu(z)\in{\mathbb{R}},$

and if $w_{2}(x)=\infty$, then

$\displaystyle w_{2}(x)$
$\displaystyle=x-\int_{\mathbb{R}}z\mathrm{d}\nu(z)\in{\mathbb{R}}.$

  3. 3.

Conversely, suppose that $w_{1}$ continues meromorphically to a neighbourhood
of $x\in{\mathbb{R}}$ and that $w_{1}(z)\in{\mathbb{R}}$ when
$z\in(x-\delta,x+\delta)\setminus\\{x\\}$ for some $\delta>0$. If
$x\in\operatorname{supp}(\mu\boxplus\nu)$, then $x$ is an isolated atom for
$\mu\boxplus\nu$.

##  Appendix G Theorem 2.1 of [BBCF17]

For the reader’s convenience, we transcribe below Theorem 2.1 in [BBCF17]. To
this end, we need a few steps of preparation.

###  G.1 Assumptions

Let $\mu_{1},\mu_{2}$ be two compactly supported Borel probability measures on
${\mathbb{R}}$. Let $\ell_{1},\ell_{2}$ be two positive integers. Let
$O\sim\operatorname{Haar}({\mathbb{O}}(d))$ be a (sequence of) random
orthogonal matrix distributed according to the Haar distribution. Let
$\theta_{1}\geq\theta_{2}\geq\cdots\theta_{\ell_{1}}$ and
$\gamma_{1}\geq\gamma_{2}\geq\cdots\geq\gamma_{\ell_{2}}$ be two fixed
sequences of real numbers which do not depend on
$\operatorname{supp}(\mu_{1})$ and $\operatorname{supp}(\mu_{2})$, and denote
$\Theta\coloneqq(\theta_{i})_{i=1}^{\ell_{1}},\Gamma=(\gamma_{i})_{i=1}^{\ell_{2}}$
as the _multiset_ of these two sequences of numbers, respectively. Let
$M=M_{1}+OM_{2}O^{\top}$ where $M_{1}$, $M_{2}$ and $O$ are matrices
satisfying the following assumptions.

  1. 1.

The matrix (sequence) $M_{1}\in{\mathbb{R}}^{d\times d}$ is deterministic and
symmetric for each $d\in{\mathbb{Z}}_{\geq 1}$ and satisfies:

    1. (a)

$\mu_{M_{1}}\to\mu_{1}$ weakly as $d\to\infty$;

    2. (b)

for every $d\geq\ell_{1}$ and $\theta\in\Theta$,

$\displaystyle\left|\left\\{i\in[d]:\lambda_{i}(M_{1})=\theta\right\\}\right|$
$\displaystyle=\left|\left\\{i\in[\ell_{1}]:\theta_{i}=\theta\right\\}\right|,$

in words, the number of eigenvalues of $M_{1}$ that are equal to $\theta$
equals the multiplicity of $\theta$ in $\Theta$;

    3. (c)

the set of eigenvalues of $M_{1}$ that do not equal any $\theta\in\Theta$
converges uniformly to $\operatorname{supp}(\mu_{1})$ as $d\to\infty$, that
is,

$\displaystyle\lim_{d\to\infty}\max_{\begin{subarray}{c}i\in[d]\\\
\lambda_{i}(M_{1})\notin\Theta\end{subarray}}\operatorname{dist}(\lambda_{i}(M_{1}),\operatorname{supp}(\mu_{1}))$
$\displaystyle=0.$

  2. 2.

The matrix (sequence) $M_{2}\in{\mathbb{R}}^{d\times d}$ is deterministic and
symmetric for each $d\in{\mathbb{Z}}_{\geq 1}$ and satisfies:

    1. (a)

$\mu_{M_{2}}\to\mu_{2}$ weakly as $d\to\infty$;

    2. (b)

for every $d\geq\ell_{2}$ and $\gamma\in\Gamma$,

$\displaystyle\left|\left\\{i\in[d]:\lambda_{i}(M_{2})=\gamma\right\\}\right|$
$\displaystyle=\left|\left\\{i\in[\ell_{1}]:\gamma_{i}=\gamma\right\\}\right|,$

in words, the number of eigenvalues of $M_{2}$ that are equal to $\gamma$
equals the multiplicity of $\gamma$ in $\Gamma$;

    3. (c)

the set of eigenvalues of $M_{2}$ that do not equal any $\gamma\in\Gamma$
converges uniformly to $\operatorname{supp}(\mu_{2})$ as $d\to\infty$, that
is,

$\displaystyle\lim_{d\to\infty}\max_{\begin{subarray}{c}i\in[d]\\\
\lambda_{i}(M_{2})\notin\Gamma\end{subarray}}\operatorname{dist}(\lambda_{i}(M_{2}),\operatorname{supp}(\mu_{2}))$
$\displaystyle=0.$

######  Remark G.1.

It was mentioned in [BBCF17, Remark 2.7] that the same results below (Theorem
G.1) hold if $M_{1}$ and $M_{2}$ are random matrices such that $M_{1},M_{2},U$
are mutually independent and the outlier eigenvalues of $M_{1}$ (resp.
$M_{2}$) converge to (instead of being exactly equal to) $\theta_{i}$ (resp.
$\gamma_{i}$) almost surely as $d\to\infty$.

######  Remark G.2.

It was mentioned in [BBCF17, Remark 2.8] that analogous results below (Theorem
G.1) can be obtained if $M$ is a sum of $k\geq 3$ matrices. Indeed, let
$k\in{\mathbb{Z}}_{\geq 2}$,
$\\{A_{i}\\}_{i\in[k]}\subset{\mathbb{R}}^{d\times d}$ be deterministic
symmetric matrices444Again, this condition can be relaxed so that $A_{i}$ are
independent random matrices with desired limiting spectral distribution and
outlier eigenvalues. See Remark G.1.  and
$(U_{i})_{i\in[k-1]}\sim\operatorname{Haar}({\mathbb{O}}(d))^{\otimes(k-1)}$
be independent Haar matrices. Suppose that for $i\in[k]$, the empirical
spectral distribution of $A_{i}$ converges weakly to $\mu_{i}$ as $d\to\infty$
and $A_{i}$ has outlier eigenvalues subject to the assumptions in Section G.1.
Let

$\displaystyle M$
$\displaystyle\coloneqq\sum_{i=1}^{k-1}U_{i}A_{i}U_{i}^{\top}+A_{k}.$

Then the limiting spectral distribution of $M$ equals
$\bigboxplus_{i=1}^{k}\mu_{i}$, and the outlier eigenvalues of $M$ can be
described in an analogous manner as in Theorem G.1 below.

###  G.2 Theorem statement

We are now ready to state Theorem 2.1 of [BBCF17].

######  Theorem G.1 ([BBCF17, Theorem 2.1]).

Let $w_{1}$ and $w_{2}$ be the subordination functions associated with the
pair of probability measures $\mu_{1}$ and $\mu_{2}$. Define

$\displaystyle\overline{\mathcal{K}}\coloneqq\operatorname{supp}(\mu_{1}\boxplus\mu_{2}),\quad\widehat{\mathcal{K}}\coloneqq\bigcup_{i=1}^{\ell_{1}}\\{w_{1}^{-1}(\theta_{i})\\}\cup\bigcup_{j=1}^{\ell_{2}}\\{w_{2}^{-1}(\gamma_{j})\\},\quad{\mathcal{K}}\coloneqq\overline{\mathcal{K}}\cup\widehat{\mathcal{K}}.$
(G.1)

Then, the following statements hold.

  1. 1.

For any $\varepsilon>0$,

$\displaystyle\Pr\left[\exists d_{0},\,\forall d>d_{0},\,\forall
i\in[d],\,\lambda_{i}(M)\in{\mathcal{K}}_{\varepsilon}\right]$
$\displaystyle=1,$

where for a set of real numbers ${\mathcal{K}}\subset{\mathbb{R}}$, we use
${\mathcal{K}}_{\varepsilon}$ to denote its _$\varepsilon$ -enlargement_:

$\displaystyle{\mathcal{K}}_{\varepsilon}\coloneqq\left\\{\rho\in{\mathbb{R}}:\inf_{\rho^{\prime}\in{\mathcal{K}}}|\rho-\rho^{\prime}|\leq\varepsilon\right\\}.$

In words, for all sufficiently large dimensions, all eigenvalues of $M$
converge into ${\mathcal{K}}$.

  2. 2.

Fix some $\rho\in\widehat{\mathcal{K}}$, let $\varepsilon>0$ be sufficiently
small so that
$(\rho-2\varepsilon,\rho+2\varepsilon)\cap\widehat{\mathcal{K}}=\\{\rho\\}$.
Set
$k\coloneqq\left|\left\\{i\in[\ell_{1}]:w_{1}(\rho)=\theta_{i}\right\\}\right|$
and
$\ell=\left|\left\\{i\in[\ell_{2}]:w_{2}(\rho)=\gamma_{i}\right\\}\right|$.
That is, $k$ and $\ell$ are the numbers of the outliers of $M_{1}$ and
$M_{2}$, respectively, that go to the outlier $\rho$ of $M$. Then

$\displaystyle\Pr\left[\exists d_{0},\,\forall
d>d_{0},\,\left|\left\\{i\in[d]:\lambda_{i}(M)\in(\rho-\varepsilon,\rho+\varepsilon)\right\\}\right|=k+\ell\right]$
$\displaystyle=1.$

In words, the multiplicity of the outlier eigenvalue $\rho$ of $M$ is $k+\ell$
if there are $k$ outliers of $M_{1}$ and $\ell$ outliers of $M_{2}$ that
correspond to $\rho$.

##  Appendix H State evolution of GAMP for mixed GLMs (proof of Proposition
6.1)

Recall the GAMP iteration in Equation 6.2:

$\displaystyle\begin{split}u^{t}&=\frac{1}{\sqrt{\delta}}\overline{A}f_{t}(v^{t};\,\overline{x}_{1}^{*},\overline{x}_{2}^{*})-{\sf
b}_{t}\,g_{t-1}(u^{t-1};y),\qquad
v^{t+1}=\frac{1}{\sqrt{\delta}}\overline{A}^{\top}g_{t}(u^{t};y)-{\sf
c}_{t}\,f_{t}(v^{t};\,\overline{x}_{1}^{*},\overline{x}_{2}^{*}),\end{split}$
(H.1)

where the memory coefficients are ${\sf
b}_{t}=\frac{1}{n}\sum_{i=1}^{d}f_{t}^{\prime}(v_{i}^{t};\,\overline{x}_{1,i}^{*},\overline{x}_{2,i}^{*})$
and ${\sf c}_{t}=\frac{1}{n}\sum_{i=1}^{n}g_{t}^{\prime}(u_{i}^{t};y_{i})$.
Using the initialization ${\tilde{v}}^{0}$, the iteration starts with
$u^{0}=\frac{1}{\sqrt{\delta}}\overline{A}{\tilde{v}}^{0}$.

To prove the proposition, we rewrite Equation H.1 as an AMP iteration with
matrix-valued iterates. This matrix-valued AMP is designed to be a special
case of an abstract AMP iteration for which a state evolution result has been
established in [JM13, FVRS22]. This state evolution result is then translated
to obtain the results in Equations 6.12 and 6.13. (See also [GB21] for an
analysis of a general graph-based AMP iteration that includes AMP with matrix-
valued iterates.) Given the iteration in Equation H.1, for $t\geq 1$, let

$\displaystyle
e^{t}\coloneqq\begin{bmatrix}g_{1}&g_{2}&u^{t}\end{bmatrix},\qquad
h^{t+1}\coloneqq\,$ $\displaystyle
v^{t+1}-\chi_{1,t+1}\overline{x}_{1}^{*}-\chi_{2,t+1}\overline{x}_{2}^{*}\,,$
(H.2)

where we recall that $g_{1}=\overline{A}\overline{x}_{1}^{*}$,
$g_{2}=\overline{A}\overline{x}_{2}^{*}$, and $\chi_{1,t},\chi_{2,t}$ are the
state evolution parameters computed via the recursion in Equations 6.4, 6.5
and 6.6. We also introduce the functions
${\breve{f}}_{t}:{\mathbb{R}}^{3}\to{\mathbb{R}}^{3}$ and
${\breve{g}}_{t}:{\mathbb{R}}^{5}\to{\mathbb{R}}$, defined as:

$\displaystyle{\breve{f}}_{t}(h^{t};\,\overline{x}_{1}^{*},\overline{x}_{2}^{*})=\begin{bmatrix}\sqrt{\delta}\overline{x}_{1}^{*},&\sqrt{\delta}\overline{x}_{2}^{*},&f_{t}(h^{t}+\chi_{1,t}\overline{x}_{1}^{*}+\chi_{2,t}\overline{x}_{2}^{*};\overline{x}_{1}^{*},\overline{x}_{2}^{*})\end{bmatrix},$
(H.3)
$\displaystyle{\breve{g}}_{t}(e^{t};\,{\underline{\eta}},{\underline{\varepsilon}})=g_{t}(e^{t}_{3};\,q({\underline{\eta}}\odot
e^{t}_{1}+(1-{\underline{\eta}})\odot e^{t}_{2},{\underline{\varepsilon}})).$
(H.4)

Here, $\odot$ denotes element-wise multiplication, ${\breve{f}}_{t}$ and
${\breve{g}}_{t}$ act row-wise on their matrix-valued inputs and $e^{t}_{j}$
denotes the $j$-th column of $e^{t}\in{\mathbb{R}}^{n\times 3}$. We also
recall the notation ${\underline{\eta}}=(\eta_{1},\ldots,\eta_{n})$,
${\underline{\varepsilon}}=(\varepsilon_{1},\ldots,\varepsilon_{n})$ and that
$y=q({\underline{\eta}}\odot g_{1}+(1-{\underline{\eta}})\odot
g_{2},\,{\underline{\varepsilon}})=q({\underline{\eta}}\odot
e^{t}_{1}+(1-{\underline{\eta}})\odot e^{t}_{2},\,{\underline{\varepsilon}})$.
With these definitions, we claim that the AMP iteration in Equation H.1 is
equivalent to the following one:

$\begin{split}&e^{t}=\frac{1}{\sqrt{\delta}}\overline{A}{\breve{f}}_{t}(h^{t};\,\overline{x}_{1}^{*},\overline{x}_{2}^{*})\,-\,{\breve{g}}_{t-1}(e^{t-1};\,{\underline{\eta}},{\underline{\varepsilon}}){\sf
B}_{t}^{\top},\\\
&h^{t+1}=\frac{1}{\sqrt{\delta}}\overline{A}^{\top}{\breve{g}}_{t}(e^{t};\,{\underline{\eta}},{\underline{\varepsilon}})\,-\,{\breve{f}}_{t}(h^{t};\,\overline{x}_{1}^{*},\overline{x}_{2}^{*}){\sf
C}_{t}^{\top},\end{split}$ (H.5)

where ${\sf B}_{t}\in{\mathbb{R}}^{3\times 1}$ and ${\sf
C}_{t}\in{\mathbb{R}}^{1\times 3}$ are given by:

$\displaystyle{\sf
B}_{t}=\begin{bmatrix}0&0&\frac{1}{n}\sum_{i=1}^{d}f_{t}^{\prime}(h_{i}^{t}+\chi_{1,t}\overline{x}^{*}_{1,i}+\chi_{2,t}\overline{x}^{*}_{2,i};\overline{x}_{1,i}^{*},\overline{x}_{2,i}^{*})\end{bmatrix}^{\top},$
$\displaystyle{\sf
C}_{t}=\begin{bmatrix}{\mathbb{E}}[\partial_{1}g_{t}(U_{t};\,q(\eta
G_{1}+(1-\eta)G_{2},\varepsilon))]\\\
{\mathbb{E}}[\partial_{2}g_{t}(U_{t};\,q(\eta
G_{1}+(1-\eta)G_{2},\varepsilon))]\\\
\frac{1}{n}\sum_{i=1}^{n}g_{t}^{\prime}(u_{i}^{t};\,q(\eta_{i}g^{1}_{i}+(1-\eta_{i})g^{2}_{i},\varepsilon_{i}))]\end{bmatrix}^{\top}.$
(H.6)

Here $\partial_{1}g_{t}$ and $\partial_{2}g_{t}$ refer to the partial
derivatives of $g_{t}(u;\,q(\eta g_{1}+\eta g_{2},\varepsilon))$ with respect
to $g_{1}$ and $g_{2}$, respectively,
$(G_{1},G_{2},\eta,\varepsilon)\sim{\mathcal{N}}(0,1)\otimes{\mathcal{N}}(0,1)\otimes\operatorname{Bern}(\alpha)\otimes
P_{\varepsilon}$, and $U_{t}$ is defined as in Equation 6.4. The iteration in
Equation H.5 is initialized with $e^{0}=[g_{1}\quad
g_{2}\quad\frac{1}{\sqrt{\delta}}\overline{A}{\tilde{v}}^{0}]$ .

The equivalence between Equations H.1 and H.5 can be seen by substituting in
Equation H.5 the definitions of $e^{t}$ and $h^{t+1}$ from Equation H.2, and
the fact that by Stein’s lemma, $\chi_{1,t+1}$ and $\chi_{2,t+1}$ defined in
Equations 6.4, 6.5 and 6.6 can be expressed as [FVRS22, Lemma 4.1]:

$\chi_{1,t+1}=\sqrt{\delta}\,{\mathbb{E}}[\partial_{1}g_{t}(U_{t};q(\eta
G_{1}+(1-\eta)G_{2},\varepsilon))],\;\chi_{2,t+1}=\sqrt{\delta}\,{\mathbb{E}}[\partial_{2}g_{t}(U_{t};q(\eta
G_{1}+(1-\eta)G_{2},\varepsilon))].$ (H.7)

The recursion in Equation H.5 is a special case of the abstract AMP recursion
with matrix-valued iterates for which a state evolution result has been
established [JM13], [FVRS22, Sec. 6.7]. The standard form of the abstract AMP
recursion uses empirical estimates (instead of expected values) for the first
two entries of ${\sf C}_{t}$ in Equation H.6. However, the state evolution
result remains valid for the recursion in Equation H.5 (see Remark 4.3 of
[FVRS22]). This result states that the empirical distributions of the rows of
$e^{t}$ and $h^{t+1}$ converge to the Gaussian distributions
${\mathcal{N}}(0,\Sigma_{t})$ and ${\mathcal{N}}(0,\Omega_{t+1})$,
respectively. The covariances $\Sigma_{t}\in{\mathbb{R}}^{3\times 3}$ and
$\Omega_{t+1}\in{\mathbb{R}}$ are defined by the following state evolution
recursion:

$\displaystyle\Sigma_{t}$
$\displaystyle=\frac{1}{\delta}{\mathbb{E}}[{\breve{f}}_{t}(G^{\omega}_{t};\,X_{1},X_{2}){\breve{f}}_{t}(G^{\omega}_{t};\,X_{1},X_{2})^{\top}]$
(H.8)
$\displaystyle=\begin{bmatrix}1&0&\frac{{\mathbb{E}}[X_{1}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})]}{\sqrt{\delta}}\\\
0&1&\frac{{\mathbb{E}}[X_{2}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})]}{\sqrt{\delta}}\\\
\frac{{\mathbb{E}}[X_{1}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})]}{\sqrt{\delta}}&\frac{{\mathbb{E}}[X_{2}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})]}{\sqrt{\delta}}&\frac{{\mathbb{E}}[(f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2}))^{2}]}{\delta}\end{bmatrix},$
$\displaystyle\Omega_{t+1}$
$\displaystyle={\mathbb{E}}[({\breve{g}}_{t}(G^{\sigma}_{t};\,\eta,\varepsilon))^{2}]={\mathbb{E}}[(g_{t}(G^{\sigma}_{t,3};\,q(G^{\sigma}_{t,1},G^{\sigma}_{t,2},\eta,\varepsilon)))^{2}],$
(H.9)

where
$G^{\sigma}_{t}\equiv(G^{\sigma}_{t,1},G^{\sigma}_{t,2},G^{\sigma}_{t,3})\sim{\mathcal{N}}(0,\Sigma_{t})$
is independent of $(\eta,\varepsilon)\sim\operatorname{Bern}(\alpha)\otimes
P_{\varepsilon}$, and
$(G^{\omega}_{t},X_{1},X_{2})\sim{\mathcal{N}}(0,\Omega_{t})\otimes{\mathcal{N}}(0,1)\otimes{\mathcal{N}}(0,1)$.
The recursion is initialized with

$\Sigma_{0}=\begin{bmatrix}1&0&\frac{1}{\sqrt{\delta}}{\mathbb{E}}[F_{0}(X_{1},X_{2})X_{1}]\\\
0&1&\frac{1}{\sqrt{\delta}}{\mathbb{E}}[F_{0}(X_{1},X_{2})X_{2}]\\\
\frac{1}{\sqrt{\delta}}{\mathbb{E}}[F_{0}(X_{1},X_{2})X_{1}]&\frac{1}{\sqrt{\delta}}{\mathbb{E}}[F_{0}(X_{1},X_{2})X_{2}]&\frac{1}{\delta}{\mathbb{E}}[(F_{0}(X_{1},X_{2}))^{2}]\end{bmatrix}.$
(H.10)

The sequences $(G^{\sigma}_{t})_{t\geq 0}$ and $(G^{\omega}_{t+1})_{t\geq 0}$
are each jointly Gaussian with the following covariance structure:

$\begin{split}&G^{\sigma}_{t,1}=G_{1},\qquad
G^{\sigma}_{t,2}=G_{2},\quad\forall\,t\geq 0\quad\text{ where
}(G_{1},G_{2})\sim{\mathcal{N}}(0,1)\otimes{\mathcal{N}}(0,1),\\\
&\mathbb{E}\left[G^{\sigma}_{0,3}G^{\sigma}_{t,3}\right]=\frac{1}{\delta}\mathbb{E}\left[F_{0}(X_{1},X_{2})f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})\right],\quad
t\geq 1,\end{split}$ (H.11)

and for $r,t\geq 1$:

$\mathbb{E}\left[G^{\omega}_{r}G^{\omega}_{t}\right]={\mathbb{E}}[g_{r-1}(G^{\sigma}_{r-1,3};\,q(\eta
G^{\sigma}_{r-1,1}+(1-\eta)G^{\sigma}_{r-1,2},\varepsilon)\\\ \times
g_{t-1}(G^{\sigma}_{t-1,3};\,q(\eta
G^{\sigma}_{t-1,1}+(1-\eta)G^{\sigma}_{t-1,2},\varepsilon)],$ (H.12a)
$\mathbb{E}\left[G^{\sigma}_{r,3}G^{\sigma}_{t,3}\right]=\frac{1}{\delta}\mathbb{E}\left[f_{r}(G^{\omega}_{r}+\chi_{1,r}X_{1}+\chi_{2,r}X_{2};X_{1},X_{2})\,f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})\right].$
(H.12b)

The following proposition follows from the state evolution result in [JM13],
[FVRS22, Sec. 6.7] for AMP with matrix-valued iterates.

######  Proposition H.1 (State Evolution).

With the setup and assumptions of Proposition 6.1, consider the AMP recursion
in Equation H.5. The following holds almost surely for any ${\rm PL}(2)$
functions $\Psi:{\mathbb{R}}^{t+1}\to{\mathbb{R}}$,
$\Phi:{\mathbb{R}}^{t+3}\to{\mathbb{R}}$, for $t\geq 0$:

$\displaystyle\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\Psi(e^{t}_{i},e^{t-1}_{i},\ldots,e^{0}_{i})={\mathbb{E}}[\Psi(G^{\sigma}_{t},\,G^{\sigma}_{t-1},\,\ldots,G^{\sigma}_{0})],$
(H.13)
$\displaystyle\lim_{d\to\infty}\frac{1}{d}\sum_{i=1}^{d}\Phi(\overline{x}^{*}_{1,i},\overline{x}^{*}_{2,i},h^{t+1}_{i},h^{t}_{i},\ldots,h^{1}_{i})={\mathbb{E}}[\Phi(X_{1},X_{2},G^{\omega}_{t+1},\,G^{\omega}_{t},\,\ldots,G^{\omega}_{1})],$
(H.14)

where the joint distributions of
$(G^{\sigma}_{t},\,G^{\sigma}_{t-1},\,\ldots,G^{\sigma}_{0})$ and
$(G^{\omega}_{t+1},\,G^{\omega}_{t},\,\ldots,G^{\omega}_{1}))$ are as given in
Equations H.8, H.9, H.10, H.11 and H.12b.

Recalling the definitions of $e^{t},h^{t+1}$ from Equation H.2, the results in
Equations H.13 and H.14 imply that for any ${\rm PL}(2)$ function
$\Phi:{\mathbb{R}}^{t+3}\to{\mathbb{R}}$, we have for $t\geq 0$:

$\displaystyle\lim_{d\to\infty}\frac{1}{d}\sum_{i=1}^{d}\Phi(\overline{x}^{*}_{1,i},\overline{x}^{*}_{2,i},v^{t+1}_{i},\ldots,v^{1}_{i})$
$\displaystyle\qquad={\mathbb{E}}[\Phi(X_{1},X_{2},\chi_{1,t+1}X_{1}+\chi_{2,t+1}X_{2}+G^{\omega}_{t+1},\,\ldots,\,\chi_{1,1}X_{1}+\chi_{2,1}X_{2}+G^{\omega}_{1})]$
(H.15)
$\displaystyle\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\Phi(g_{1,i},g_{2,i},u^{t}_{i},\ldots,u^{0}_{i})=\mathbb{E}\left[\Phi(G_{1},G_{2},G^{\sigma}_{t,3},\ldots,G^{\sigma}_{0,3})\right].$
(H.16)

Recalling $G^{\sigma}_{t,1}=G_{1}$ and $G^{\sigma}_{t,2}=G_{2}$ for $t\geq 0$,
using Equation H.8 we can write

$\displaystyle G^{\sigma}_{t,3}$
$\displaystyle=\mathbb{E}\left[G^{\sigma}_{t,3}\,|\,G^{\sigma}_{t,1},G^{\sigma}_{t,2}\right]\,+\,W^{\sigma}_{t}$
$\displaystyle=\frac{1}{\sqrt{\delta}}\mathbb{E}\left[X_{1}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})\right]G_{1}$
$\displaystyle\phantom{=}\,+\,\frac{1}{\sqrt{\delta}}\mathbb{E}\left[X_{2}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})\right]G_{2}\,+\,W^{\sigma}_{t},$

where $W^{\sigma}_{t}$ is a zero-mean Gaussian independent of $(G_{1},G_{2})$
with variance

$\mathbb{E}\left[(W^{\sigma}_{t})^{2}\right]=\frac{1}{\delta}{\mathbb{E}}[(f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2}))^{2}]\\\
-\frac{1}{\delta}(\mathbb{E}\left[X_{1}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})\right])^{2}-\frac{1}{\delta}(\mathbb{E}\left[X_{2}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})\right])^{2}.$
(H.17)

To complete the proof, for $t\geq 0$, we define:

$\begin{split}&W_{U,t}\coloneqq W^{\sigma}_{t},\quad W_{V,t+1}\coloneqq
G^{\omega}_{t+1},\\\
&\mu_{1,t}\coloneqq\frac{1}{\sqrt{\delta}}{\mathbb{E}}[X_{1}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})],\quad\mu_{2,t}\coloneqq\frac{1}{\sqrt{\delta}}{\mathbb{E}}[X_{2}f_{t}(G^{\omega}_{t}+\chi_{1,t}X_{1}+\chi_{2,t}X_{2};X_{1},X_{2})],\\\
&U_{t}\coloneqq G^{\sigma}_{3,t}=\mu_{1,t}G_{1}+\mu_{2,t}G_{2}+W_{U,t},\quad
V_{t+1}=\chi_{1,t}X_{1}+\chi_{2,t}X_{2}+W_{V,t+1}.\end{split}$

Using these in the convergence statements in Equations H.15 and H.16 gives the
result of Proposition 6.1.

##  Appendix I Auxiliary lemmas

######  Lemma I.1 (Conditional distribution of Gaussians).

Let $d\geq 2$ and $1\leq p\leq d-1$ be integers. Let

$\displaystyle\begin{bmatrix}G_{1}\\\ G_{2}\end{bmatrix}$
$\displaystyle\sim{\mathcal{N}}\left(\begin{bmatrix}\mu_{1}\\\
\mu_{2}\end{bmatrix},\begin{bmatrix}\Sigma_{1,1}&\Sigma_{1,2}\\\
\Sigma_{1,2}^{\top}&\Sigma_{2,2}\end{bmatrix}\right)$

be a $d$-dimensional Gaussian random vector, where
$G_{1}\in{\mathbb{R}}^{p},G_{2}\in{\mathbb{R}}^{d-p},\mu_{1}\in{\mathbb{R}}^{p},\mu_{2}\in{\mathbb{R}}^{d-p},\Sigma_{1,1}\in{\mathbb{R}}^{p\times
p},\Sigma_{1,2}\in{\mathbb{R}}^{p\times(d-p)},\Sigma_{2,2}\in{\mathbb{R}}^{(d-p)\times(d-p)}$.
Then for any $g_{2}\in{\mathbb{R}}^{d-p}$, the distribution of $G_{1}$
conditioned on $G_{2}=g_{2}$ is given by
$G_{1}\,|\,\\{G_{2}=g_{2}\\}\sim{\mathcal{N}}(\mu_{1}^{\prime},\Sigma_{1}^{\prime})$
where

$\displaystyle\mu_{1}^{\prime}$
$\displaystyle=\mu_{1}+\Sigma_{1,2}\Sigma_{2,2}^{-1}(g_{2}-\mu_{2})\in{\mathbb{R}}^{p},\quad\Sigma_{1}^{\prime}=\Sigma_{1,1}-\Sigma_{1,2}\Sigma_{2,2}^{-1}\Sigma_{1,2}^{\top}\in{\mathbb{R}}^{p\times
p}$

and $\Sigma_{2,2}^{-1}$ denotes the generalized inverse of $\Sigma_{2,2}$.

######  Lemma I.2.

Consider the setting of Section 2.4 and let Item (A7) hold. Then the following
properties of $\varphi(\lambda),\psi(\lambda;\Delta)$ hold.

  1. 1.

$\varphi(\cdot)$ is strictly decreasing;

  2. 2.

For any $\Delta>0$, $\psi(\cdot;\Delta)$ is strictly convex in the first
argument;

  3. 3.

For any $\lambda>\sup\operatorname{supp}(Z)$, $\psi(\lambda;\cdot)$ is
strictly decreasing in the second argument.

###### Proof.

The proof follows by checking the derivatives. For $\varphi$, we have

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\lambda}\varphi(\lambda)$
$\displaystyle=\mathbb{E}\left[\frac{ZG^{2}}{\lambda-Z}\right]-\lambda\mathbb{E}\left[\frac{ZG^{2}}{(\lambda-Z)^{2}}\right]=-\mathbb{E}\left[\left(\frac{ZG}{\lambda-Z}\right)^{2}\right]<0,$
(I.1)

The last strict inequality holds since $Z$ is not almost surely zero, i.e.,
$\Pr\left[Z=0\right]<1$ in Item (A7). Item 1 of Lemma I.2 then follows. For
$\psi$, we have

$\displaystyle\frac{\partial}{\partial\lambda}\psi(\lambda;\Delta)$
$\displaystyle=\frac{1}{\Delta}+\mathbb{E}\left[\frac{Z}{\lambda-Z}\right]-\lambda\mathbb{E}\left[\frac{Z}{(\lambda-Z)^{2}}\right]=\frac{1}{\Delta}-\mathbb{E}\left[\frac{Z^{2}}{(\lambda-Z)^{2}}\right],$
(I.2)

and

$\displaystyle\frac{\partial^{2}}{\partial\lambda^{2}}\psi(\lambda;\Delta)$
$\displaystyle=2\mathbb{E}\left[\frac{Z^{2}}{(\lambda-Z)^{3}}\right].$

Since $\lambda>\sup\operatorname{supp}(Z)$ and $Z$ is not almost surely zero,
$\psi^{\prime\prime}(\cdot;\Delta)>0$ and therefore Item 2 of Lemma I.2 holds.
Finally, Item 3 of Lemma I.2 is obvious, provided
$\lambda>\sup\operatorname{supp}(Z)>0$ where the second inequality is
guaranteed by Item (A7). ∎

######  Lemma I.3 (Explicit formulas).

Fix any $\Delta>0$. The parameter $\overline{\lambda}(\Delta)$ satisfies

$\displaystyle\mathbb{E}\left[\left(\frac{Z}{\overline{\lambda}(\Delta)-Z}\right)^{2}\right]$
$\displaystyle=\frac{1}{\Delta}.$ (I.3)

If $\lambda^{*}(\Delta)>\overline{\lambda}(\Delta)$, the parameter
$\lambda^{*}(\Delta)$ satisfies

$\displaystyle\mathbb{E}\left[\frac{Z(G^{2}-1)}{\lambda^{*}(\Delta)-Z}\right]$
$\displaystyle=\frac{1}{\Delta}.$ (I.4)

###### Proof.

Since $\overline{\lambda}(\Delta)$ (cf. Equation 3.3) is the minimum point of
$\psi(\cdot;\Delta)$, it satisfies

$\displaystyle\left.\frac{\partial}{\partial\lambda}\psi(\lambda;\Delta)\right|_{\lambda=\overline{\lambda}(\Delta)}$
$\displaystyle=0.$

This gives Equation I.3 according to Equation I.2.

Under the condition $\lambda^{*}(\Delta)>\overline{\lambda}(\Delta)$, we have
$\zeta(\lambda^{*}(\Delta);\Delta)=\psi(\lambda^{*}(\Delta);\Delta)$ and
$\lambda^{*}(\Delta)$ satisfies the fixed point equation
$\psi(\lambda^{*}(\Delta);\Delta)=\varphi(\lambda^{*}(\Delta))$ which in turn
can be written more explicitly as Equation I.4. ∎

######  Remark I.1.

Equation I.4 is often used with $\Delta=\delta_{i}$ ($i\in\\{1,2\\}$) under
the condition $\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$. This is
legitimate since $\overline{\lambda}(\delta)>\overline{\lambda}(\delta_{i})$
(see, e.g., Equation 5.33) and the latter condition is stronger than the one
in Lemma I.3.

Lemmas I.4 and I.5 below follow by chasing the definitions of various
functions and parameters defined at the beginning of Section 3. We omit the
proofs.

######  Lemma I.4 (Equivalent characterization of spectral threshold).

For $i\in\\{1,2\\}$, the following are equivalent.

  1. 1.

$\lambda^{*}(\delta_{i})>\overline{\lambda}(\delta)$;

  2. 2.

$\psi(\overline{\lambda}(\delta);\delta_{i})<\varphi(\overline{\lambda}(\delta))$;

  3. 3.

$\zeta(\lambda^{*}(\delta_{i});\delta_{i})>\zeta(\overline{\lambda}(\delta);\delta_{i})$;

  4. 4.

$\psi^{\prime}(\lambda^{*}(\delta_{i});\delta)>0$.

######  Lemma I.5 (Alternative formulas of overlaps, spectral).

If $\lambda^{*}(\delta_{1})>\overline{\lambda}(\delta)$, then

$\displaystyle\sqrt{\frac{1}{\alpha}\cdot\frac{\frac{1}{\delta}-\mathbb{E}\left[\left(\frac{Z}{\lambda^{*}(\delta_{1})-Z}\right)^{2}\right]}{\frac{1}{\delta_{1}}+\mathbb{E}\left[\left(\frac{Z}{\lambda^{*}(\delta_{1})-Z}\right)^{2}(G^{2}-1)\right]}}$
$\displaystyle=\sqrt{\frac{\psi^{\prime}(\lambda^{*}(\delta_{1});\delta)}{\alpha(\psi^{\prime}(\lambda^{*}(\delta_{1});\delta_{1})-\varphi^{\prime}(\lambda^{*}(\delta_{1})))}}.$

If $\lambda^{*}(\delta_{2})>\overline{\lambda}(\delta)$, then

$\displaystyle\sqrt{\frac{1}{1-\alpha}\cdot\frac{\frac{1}{\delta}-\mathbb{E}\left[\left(\frac{Z}{\lambda^{*}(\delta_{2})-Z}\right)^{2}\right]}{\frac{1}{\delta_{2}}+\mathbb{E}\left[\left(\frac{Z}{\lambda^{*}(\delta_{2})-Z}\right)^{2}(G^{2}-1)\right]}}$
$\displaystyle=\sqrt{\frac{\psi^{\prime}(\lambda^{*}(\delta_{2});\delta)}{(1-\alpha)(\psi^{\prime}(\lambda^{*}(\delta_{2});\delta_{2})-\varphi^{\prime}(\lambda^{*}(\delta_{2})))}}.$

######  Theorem I.6 ([MM19, Lemma 2]).

Let $x^{*}\sim\operatorname{Unif}({\mathbb{S}}^{d-1})$ be a signal vector. Let
$A\in{\mathbb{R}}^{n\times d}$ be a random design matrix with i.i.d. entries
distributed according to ${\mathcal{N}}(0,1)$. Let
${\mathcal{T}}\colon{\mathbb{R}}\to{\mathbb{R}}$ be a preprocessing function
and $T=\mathrm{diag}({\mathcal{T}}(q(Ax^{*},{\underline{\varepsilon}})))$,
where $q\colon{\mathbb{R}}\times{\mathbb{R}}\to{\mathbb{R}}$ is a link
function and ${\underline{\varepsilon}}\sim P_{\varepsilon}^{\otimes n}$ is a
noise sequence. Suppose $n/d\to\delta$ for a constant $\delta>0$. Define
random variables $(G,\varepsilon)\sim{\mathcal{N}}(0,1)\otimes
P_{\varepsilon}$, $Y=q(G,\varepsilon)$ and $Z={\mathcal{T}}(Y)$. Assume
$\Pr\left[Z=0\right]<1$ and $Z$ has bounded support.

Define functions
$\varphi\colon(\sup\operatorname{supp}(Z),\infty)\to{\mathbb{R}}$ and
$\psi\colon(\sup\operatorname{supp}(Z),\infty)\times(0,\infty)\to{\mathbb{R}}$
as

$\displaystyle\varphi(\lambda)$
$\displaystyle=\lambda\,\mathbb{E}\left[\frac{ZG^{2}}{\lambda-Z}\right],\quad\psi(\lambda;\delta)=\lambda\left(\frac{1}{\delta}+\mathbb{E}\left[\frac{Z}{\lambda-Z}\right]\right),$

respectively. Define $\overline{\lambda}(\delta)$ and
$\zeta\colon(\sup\operatorname{supp}(Z),\infty)\times(0,\infty)\to{\mathbb{R}}$
as

$\displaystyle\overline{\lambda}(\delta)$
$\displaystyle=\operatorname*{argmin}_{\lambda>\sup\operatorname{supp}(Z)}\psi(\lambda;\delta),\quad\zeta(\lambda;\delta)=\psi(\max\\{\lambda,\overline{\lambda}(\delta)\\};\delta).$

Consider the matrix $D=\frac{1}{n}A^{\top}TA\in{\mathbb{R}}^{d\times d}$. Let
$v_{1}(D),\cdots,v_{d}(D)$ and $\lambda_{1}(D)\geq\cdots\geq\lambda_{d}(D)$
denote the eigenvectors and the corresponding eigenvalues of $D$,
respectively. Then, the following results hold:

  1. 1.

The equation $\zeta(\lambda;\delta)=\varphi(\lambda)$ has a unique solution in
$\lambda\in(\sup\operatorname{supp}(Z),\infty)$ which is denoted by
$\lambda^{*}(\delta)$.

  2. 2.

We have that, almost surely,

$\displaystyle\lim_{d\to\infty}\frac{\left|\left\langle
v_{1}(D),x^{*}\right\rangle\right|}{\left\|v_{1}(D)\right\|_{2}\left\|x^{*}\right\|_{2}}$
$\displaystyle=\begin{cases}\sqrt{\frac{\psi^{\prime}(\lambda^{*}(\delta);\delta)}{\psi^{\prime}(\lambda^{*}(\delta);\delta)-\varphi^{\prime}(\lambda^{*}(\delta))}},&\psi^{\prime}(\lambda^{*}(\delta);\delta)>0,\\\
0,&\psi^{\prime}(\lambda^{*}(\delta);\delta)\leq 0.\end{cases}$

  3. 3.

We have, almost surely,

$\displaystyle\lim_{d\to\infty}\lambda_{1}(D)$
$\displaystyle=\zeta(\lambda^{*}(\delta);\delta),\quad\lim_{d\to\infty}\lambda_{2}(D)=\zeta(\overline{\lambda}(\delta);\delta).$

######  Remark I.2.

The above theorem is a generalization of [LL20, Theorem 1] which requires
additionally that the preprocessing function can only take non-negative
values.

######  Theorem I.7 ([MM19, Lemma 3]).

Suppose that the following conditions hold:

  1. 1.

$A\in{\mathbb{C}}^{n\times d}$ is a random matrix whose entries are i.i.d.
such that $\mathbb{E}\left[A_{i,j}\right]=0$,
$\mathbb{E}\left[|A_{i,j}|^{2}\right]=1$ and
$\mathbb{E}\left[|A_{i,j}|^{4}\right]<\infty$;

  2. 2.

The empirical spectral distribution of $T\in{\mathbb{C}}^{n\times n}$
converges weakly to a probability measure $\mu$ as $n\to\infty$;

  3. 3.

$n/d\to\delta\in(0,\infty)$ as $n,d\to\infty$;

  4. 4.

The spectral norm of $T$ is bounded for every $n\in{\mathbb{Z}}_{\geq 0}$.

Then, the empirical spectral distribution of $\frac{1}{d}A^{\top}TA$ converges
weakly to a measure whose inverse Stieltjes transform is given by

$\displaystyle s^{-1}(z)$
$\displaystyle=-\frac{1}{z}+\delta\mathbb{E}\left[\frac{Z}{1+zZ}\right],$

for $z\in{\mathbb{H}}$.

######  Remark I.3.

The above theorem was first proved in [MP67] (see also [BS10, Chapter 4])
under an additional assumption that $T\in{\mathbb{C}}^{n\times n}$ is PSD for
all $n\in{\mathbb{Z}}_{\geq 0}$. It was recently showed in [MM19] that the
same conclusion holds without the PSD condition on $T$. For a general diagonal
$T\in{\mathbb{R}}^{n\times n}$ (not necessarily PSD), the idea in [MM19] is to
decompose $T$ into $T_{+}$, consisting of only non-negative entries, and
$-T_{-}$, consisting of only negative entries, and then apply free probability
tools to study how they jointly contribute to the matrix product
$\frac{1}{d}A^{\top}TA$.

######  Lemma I.8 (Lemma 4.5 in [DSS11]).

Let $(Q_{t})_{t\geq 1}$ be a sequence of distributions converging weakly to
some distribution $Q$, and let $h$ be a non-negative continuous function such
that

$\lim_{t\to\infty}\int h\,{\rm d}Q_{t}=\int h\,{\rm d}Q.$ (I.5)

Then, for any continuous function $\Psi$ such that $\left|\Psi\right|/(1+h)$
is bounded,

$\lim_{t\to\infty}\int\Psi\,{\rm d}Q_{t}=\int\Psi\,{\rm d}Q.$ (I.6)


# Real-time detection and control of correlated charge tunneling in a quantum
dot

Johannes C. Bayer  Institut für Festkörperphysik, Leibniz Universität
Hannover, Hannover, Germany  Fredrik Brange  Department of Applied Physics,
Aalto University, 00076 Aalto, Finland  Adrian Schmidt  Institut für
Festkörperphysik, Leibniz Universität Hannover, Hannover, Germany  Timo Wagner
Institut für Festkörperphysik, Leibniz Universität Hannover, Hannover, Germany
Eddy P. Rugeramigabo  Institut für Festkörperphysik, Leibniz Universität
Hannover, Hannover, Germany  Christian Flindt<EMAIL_ADDRESS>Department of Applied Physics, Aalto University, 00076 Aalto, Finland  RIKEN
Center for Quantum Computing, Wakoshi, Saitama 351-0198, Japan  Rolf J. Haug
<EMAIL_ADDRESS>Institut für Festkörperphysik, Leibniz Universität
Hannover, Hannover, Germany

###### Abstract

The tunneling of interacting charges in nanostructures is a random and noisy
process, which gives rise to a variety of intriguing physical phenomena,
including Coulomb blockade, the Kondo effect, and Fermi-edge singularities.
Such correlation effects are visible in low-frequency measurements of the
electric current or the shot noise, however, they are also expected to leave
clear fingerprints in the tunneling times of individual charges. Here, we
experimentally demonstrate the real-time detection and control of correlated
charge tunneling in a dynamically driven quantum dot. To this end, we measure
the joint distribution of waiting times between tunneling charges, and we show
that subsequent waiting times are strongly correlated due to the Coulomb
interactions between the tunneling charges and the external periodic drive.
Our measurements are in excellent agreement with a theoretical model that
allows us to develop a detailed understanding of the correlated tunneling
events. We also demonstrate that the degree of correlations can be controlled
by the external drive. Our experiment paves the way for systematic real-time
investigations of correlated electron transport in low-dimensional
nanostructures with potential applications in metrology and sensing.

Introduction. Correlations in random processes are ubiquitous in nature and
are important across a wide range of sciences and technologies [1]. A
prominent example is that of entanglement in quantum physics, where the
measurement on one particle can affect a subsequent measurement on another
particle, even if they are far apart [2]. Correlations can also occur due to
Coulomb interactions in nanostructures as the tunneling of one electron blocks
the tunneling of the next one, which is visible in noise measurements [3]. The
significance of correlations can best be appreciated with a simple example as
in Fig. 1: If a fair coin is tossed, it has an equal chance of landing on
either side up, and there are no correlations between the tosses with each
coin toss being independent of the others. As a comparison, one can imagine an
unfair coin that in each toss alternates between heads and tails. Such an
unfair coin would yield heads half of the time and tails the other half, just
as a fair coin. However, if one investigates the correlations between
subsequent coin tosses, one would find that heads is always followed by tails,
and tails is always followed by heads.

Figure 1: Correlations in random processes. a, Tossing a fair coin leads to a
random sequence of heads and tails. One half of the tosses yield heads and the
other half tails. Also, there are no correlations between the tosses, such
that heads and tails are equally likely to be followed by heads or tails. b,
One can instead consider an unfair coin that produces an alternating sequence
of heads and tails. Also in that case, one half of the tosses yield heads and
the other half tails. However, in this case, there are correlations between
the coin tosses, since heads is always followed by tails and vice versa.
Figure 2: Measurements of waiting times. a, Colored SEM image of the device. A
quantum dot (yellow circle) is formed by negative voltages on the gates
$V_{\mathrm{T1,P,T2,C1,C2}}$. A quantum point contact (QPC) serves as a charge
detector, which can be tuned by the gate $V_{\mathrm{q}}$. Ohmic contacts to
the quantum dot are indicated by the crossed squares, and single electrons
tunnel on and off the quantum dot from two electrodes. A small bias voltage
$V_{\mathrm{qd}}=$-1\text{\,}\mathrm{m}\mathrm{V}$$ drives a current through
the quantum dot. The green gates $V_{\mathrm{T1,P,T2}}$ are used to modulate
the tunneling rates periodically in time. b, Charge stability diagram of the
device. The quantum dot is driven periodically around the operating point on
the first charging line, where the charge states with $N_{e}=0,1$ electrons on
the quantum dot are energetically degenerate (red circle). Here, we use
$\Delta V_{\mathrm{T}1}=\Delta V_{\mathrm{T}2}=\Delta V_{\mathrm{T}}$. c, The
time-dependent current through the QPC switches between two values
corresponding to the number of electrons on the quantum dot, $N_{e}=0,1$. The
waiting times between tunneling events are indicated. d, Time-dependent
tunneling rates in and out of the quantum dot. The rate for tunneling out is
roughly constant, while the rate for tunneling in depends exponentially on the
gate voltages as in Eq. (3). e, Waiting time distribution for electrons
tunneling out of the quantum dot. The experimental data are well-captured by
Eq. (4) corresponding to a device with constant tunneling rates. f, Waiting
time distribution for electrons tunneling into the quantum dot or,
equivalently, holes tunneling out of the quantum dot. In this case, the
experimental data are well-captured by Eq. (5), which contains an oscillatory
factor due to the time-dependent rate at which holes tunnel.

In the context of quantum transport in nanostructures, several statistical
tools have been used to describe the tunneling of electrons [4]. In one
approach, one considers the number of electrons that have been transmitted
through a mesoscopic conductor, which is known as the full counting statistics
[5, 6, 7, 8, 9, 10]. Alternatively, one may investigate the first-passage time
distribution, which characterizes the time it takes for a given number of
electrons to have tunneled [11, 12, 13]. In a complementary approach, one
measures the distribution of waiting times between tunneling events [14, 15,
16]. In the past, photon waiting times have been investigated theoretically in
quantum optics, however, they have rarely been measured since highly efficient
detectors are needed [17, 18]. By contrast, measurements of electron waiting
times have become possible thanks to the development of accurate charge
detectors. The time between tunneling charges can be described by a
probability density, $\mathcal{W}(\tau)$, which has been measured in several
experiments [19, 20, 21, 22, 23]. However, to identify correlations between
waiting times, it is necessary to measure their joint distribution. The joint
distribution is the probability density of first observing a waiting time of
duration $\tau_{1}$ followed by a second waiting time of duration $\tau_{2}$.
It factorizes as
$\mathcal{W}(\tau_{1},\tau_{2})=\mathcal{W}(\tau_{1})\mathcal{W}(\tau_{2})$ in
the absence of correlations. Thus, to quantify correlations, one may consider
the difference

$\Delta\mathcal{W}(\tau_{1},\tau_{2})=\mathcal{W}(\tau_{1},\tau_{2})-\mathcal{W}(\tau_{1})\mathcal{W}(\tau_{2}),$
(1)

and the correlation ratio, which we define as

$\mathcal{R}(\tau_{1},\tau_{2})=\frac{\Delta\mathcal{W}(\tau_{1},\tau_{2})}{\mathcal{W}(\tau_{1})\mathcal{W}(\tau_{2})}.$
(2)

The correlation ratio is positive (or negative) for positively (or negatively)
correlated waiting times and vanishes in the absence of correlations. Joint
waiting time distributions have not been measured before, however, it has been
predicted that the waiting times for non-interacting electrons should be
correlated due to the Pauli exclusion principle [24], and a periodic drive may
also induce correlations between waiting times [25].

Figure 3: Correlated waiting times. a, Measurements of the joint waiting time
distribution for the tunneling of electrons together with model calculations.
b, The corresponding correlation ratio shows that there essentially are no
correlations between subsequent waiting times. c, Measurements of the joint
waiting time distribution for the tunneling of holes together with model
calculations. d, For the tunneling of holes, there are strong correlations
between the waiting times.

Quantum dot device. Figure 2a shows our device, which is based on a
GaAs/AlGaAs heterostructure with patterned metallic top gates and operated at
a temperature of $1.5\,\text{K}$. Negative voltages applied to the gate
electrodes lead to the formation of a quantum dot in the two-dimensional
electron gas about $110\,\text{nm}$ below the surface. In addition to the
constant offset voltages, we modulate the gate voltages in time to modify the
rate at which electrons tunnel back and forth between the quantum dot and two
external electrodes. To detect the electrons on the quantum dot, we use an
adjacent quantum point contact (QPC), whose conductance changes as electrons
tunnel in and out of the quantum dot [26, 27]. As shown in Fig. 2b, we can map
out the charge stability diagram of the quantum dot by monitoring the
conductance of the QPC as a function of the gate voltages. In the following,
we operate the quantum dot close to the point, where the charge states with
zero or one electrons on the quantum dot are energetically degenerate. As
shown in Fig. 2c, the time-resolved current in the QPC then switches back and
forth between two distinct values as electrons tunnel in and out of the
quantum dot. From such time-traces of the current, we can extract the rates at
which electrons tunnel, and how they are modulated by the time-dependent gate
voltages [28]. In Fig. 2d, we show the rates for tunneling in and out of the
dot. The rate for tunneling out of the dot is roughly constant with magnitude
$\Gamma_{o}\simeq 3.4$ kHz. On the other hand, the rate for tunneling into the
dot depends exponentially on the gate voltages as

$\Gamma_{i}(t)=\Gamma_{i}e^{\alpha\sin(2\pi ft)}/I_{0}[\alpha],$ (3)

where $f=8$ kHz is the driving frequency, and $\alpha\simeq 1.33$
characterizes the influence of the gate voltages on the tunneling rate [29,
22]. In the expression above, we have included the Bessel function,
$I_{0}[\alpha]=\int_{0}^{2\pi}dxe^{\alpha\sin x}/2\pi$, so that
$\Gamma_{i}\simeq 3.3$ kHz is the average tunneling rate. The parameter values
are found by fitting Eq. (3) to the experimental results in Fig. 2d, and these
tunneling rates form the basis for the modelling of our experiment.

Distribution of waiting times. From the time-traces of the current in the QPC,
we can extract the distribution of waiting times between subsequent tunneling
events. In Fig. 2e, we show the distribution of waiting times between
electrons tunneling out of the quantum dot and into one of the two electrodes.
This waiting time distribution depends only weakly on the periodic drive, and
it is well-approximated by the simple expression

$\mathcal{W}_{e}(\tau)\simeq\mathcal{W}_{s}(\tau)=\frac{\Gamma_{i}\Gamma_{o}}{\Gamma_{i}-\Gamma_{o}}(e^{-\Gamma_{o}\tau}-e^{-\Gamma_{i}\tau}).$
(4)

which is the distribution for a static device with constant tunneling rates
[14]. We obtain this expression by assuming that the tunneling rates are
roughly constant between each tunneling event (see Methods), and the
experimental data are well-captured by this approximation.

In Fig. 2f, we show the distribution of waiting times between electrons
tunneling into the quantum dot. We can equivalently think of these results as
the waiting times between holes tunneling into the electrodes. For the hole
waiting times, we find the expression

$\mathcal{W}_{h}(\tau)\simeq\mathcal{W}_{s}(\tau)I_{0}[2\alpha\cos(\pi
f\tau)]/I_{0}^{2}[\alpha],$ (5)

which again is given by the static distribution from Eq. (4), however, this
time, it is multiplied by an additional oscillatory function, which accounts
for the time-dependent rate at which holes tunnel (see Methods). As seen in
Fig. 2f, this prediction agrees well with the measurements of the waiting time
distribution.

Correlated waiting times. Although the results in Fig. 2e,f fully characterize
the distributions of waiting times, they do not capture possible correlations
between subsequent waiting times. Similarly, in the example of coin tosses,
both the fair and the unfair coin yield heads (or tails) half of the time.
However, for the unfair coin, that does not capture the correlations between
the coin tosses. To identify correlations between the waiting times, we need
to consider their joint distribution and the correlation ratio in Eq. (2).
Thus, in Fig. 3, we show measurements of the joint distributions of subsequent
waiting times. Specifically, in Fig. 3a,b, we show the joint waiting time
distribution for electrons and the correlation ratio. The experimental data
are well-captured by our model calculations, and the results resemble what one
would expect for a static device with no correlations. Indeed, turning to the
correlation ratio in Fig. 3b, we see that there essentially are no
correlations, as the correlation ratio is vanishingly small and mostly
dominated by experimental noise. Theoretically, we also expect that the
correlations are small as seen in the figure. Moreover, if we assume that the
tunneling rates are constant in between tunneling events, we predict that

$\mathcal{R}_{e}(\tau_{1},\tau_{2})\simeq 0,$ (6)

which is in good agreement with the measurements.

Figure 4: Joint distributions and correlation ratio. a, Joint distribution
along the first dashed line in Fig. 2c with $\tau_{1}=0.375$ ms. We also show
the distribution without correlations, given by the single waiting time
distribution. b, Correlation ratio along the first dashed line in Fig. 2d with
the theory given by Eq. (7). c, Joint distribution along the second dashed
line in Fig. 2c with $\tau_{1}=0.4375$ ms. We also show the distribution
without correlations, given by the single waiting time distribution. d,
Correlation ratio along the second dashed line in Fig. 2d with the theory
given by Eq. (7).  Figure 5: Control of correlations. a, Charge stability
diagram showing the direction of the periodic drive with respect to the
single-electron charging line. b, Corresponding correlation ratios for the
electron waiting times, which increase as we change the direction of the
periodic drive. The left panel corresponds to the results in Fig. 3b, however,
they are shown on a different scale and for a different range of times.

In Fig. 3c, a different picture emerges for the joint waiting time
distribution of holes. The joint distribution exhibits an oscillatory behavior
due to the time-dependent rate at which holes tunnel, and it is well-captured
by our model calculations. Moreover, the correlation ratio in Fig. 3d reveals
that there are strong correlations between subsequent waiting times, which can
be both positive (red) and negative (blue), such that the measurement of a
particular waiting time will likely be followed by a waiting time that is
longer or shorter than one would have expected. We can again calculate the
correlation ratio with the simplification that the tunneling rates are
constant in between tunneling events. We then find

$\mathcal{R}_{h}(\tau_{1},\tau_{2})\simeq
I_{0}[\alpha]\frac{I_{0}\big{[}\alpha\sqrt{3+2\cos(2\pi f\tau_{1})+4\cos(\pi
f\tau_{1})\cos(\pi f(\tau_{1}+2\tau_{2}))}\big{]}}{I_{0}[2\alpha\cos(\pi
f\tau_{1})]I_{0}[2\alpha\cos(\pi f\tau_{2})]}-1,$ (7)

which, surprisingly perhaps, is independent of the tunneling rates. It
captures the measurements in Fig. 3d, which can also be seen in Fig. 4, where
we show the joint waiting time distribution and the correlation ratio along
the two lines indicated in Fig. 3c,d, corresponding to specific measurement
outcomes for the first waiting time. Along the first line, the joint
distribution in Fig. 4a is oscillating, and it resembles the distribution for
uncorrelated waiting times. For that reason, the correlation ratio in Fig. 4b
is rather small with periods of positive and negative correlations. By
contrast, in Fig. 4c, the joint waiting time distribution along the second
line shows that the measurement outcome for the first waiting time strongly
affects the following waiting time. Indeed, Figure 4d reveals much stronger
correlations between the waiting times.

Control of correlations. Finally, we show that we can control the degree of
correlations between the waiting times. To this end, we change the direction
of the periodic drive as indicated in Fig. 5a. Specifically, we apply periodic
gate voltages of the form

$\begin{bmatrix}V_{\text{P}}(t)\\\
V_{\text{T}}(t)\end{bmatrix}=\begin{bmatrix}V_{\text{P}}^{0}\\\
V_{\text{T}}^{0}\end{bmatrix}+\begin{bmatrix}\Delta V_{\text{P}}\\\ \Delta
V_{\text{T}}\end{bmatrix}\sin{\left(2\pi ft\right)},$ (8)

where the first vector contains a constant offset, while the second vector
defines a direction in the charge stability diagram as indicated in Fig. 5a.
The corresponding correlation ratios for electrons are shown in Fig. 5b, and
we see how the correlations increase as we change the drive across the charge-
degeneracy line and eventually apply the drive in a perpendicular direction to
the charge-degeneracy line. These results demonstrate how the degree of
correlations can be controlled by the drive.

Summary. We have measured the joint distribution of waiting times between
tunneling charges and identified strong correlations, which arise because of
Coulomb interactions and an external periodic drive. We have also shown how
the degree of correlations can be controlled by changing the periodic drive.
We have performed real-time measurements of the tunnel processes in a low-
dimensional structure, however, future experiments may seek to detect
correlations between charge pulses that propagate in elongated electronic wave
guides [30]. Our experiment establishes a systematic approach for real-time
investigations of correlated quantum transport, which is not limited to
electrons and holes, but may also be extended to other discrete quanta such as
single photons [31] or phonons [32]. Our findings may also be important for
metrology and sensing applications, where accurate correlation measurements
are critical.

## Methods

Experiment. Our QD device is based on a two-dimensional electron gas (2DEG),
which is formed in a GaAs/AlGaAs heterostructure approximately
$110\text{\,}\mathrm{n}\mathrm{m}$ below the surface. The QD is formed in the
2DEG by applying negative voltages to metallic gates, which are patterned on
top of the heterostructure electron beam lithography. The measurements were
carried out in a 4He cryostat at $1.5\text{\,}\mathrm{K}$ in a dc transport
set-up. A $100\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ bandwidth low-noise
FEMTO transimpedance amplifier
($100\text{\,}\mathrm{M}\mathrm{V}\mathrm{A}^{-1}$ gain) was used to amplify
the QPC current at room temperature. Signals were applied and recorded with a
sampling rate of $400\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ using an ADwin
Pro2 real-time controller ($1\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ ADSP
T12).

The rates for tunneling in and out of the dot were extracted from the time-
traces of the QPC current [28]. For small time bins $\Delta t$, the rates can
be found as

$\Gamma_{\text{o}}(t)=\frac{1}{\Delta
t}\frac{N_{\text{out}}(t)}{N_{1}(t)},\quad\Gamma_{\text{i}}(t)=\frac{1}{\Delta
t}\frac{N_{\text{in}}(t)}{N_{0}(t)},$ (9)

where $N_{\text{in/out}}(t)$ is the number of tunneling events in and out of
the QD during the time $\Delta t$, and $N_{0/1}(t)$ measures how long the QD
was occupied by zero or one electron.

Model. The system is described by the rate equation

$\frac{d}{dt}|p(t)\rangle=\mathbf{L}(t)|p(t)\rangle,$ (10)

where the vector $|p(t)\rangle=[p_{0}(t),p_{1}(t)]^{T}$ contains the
probabilities for the quantum dot to be empty or occupied. The dynamics is
governed by the rate matrix

$\mathbf{L}(t)=\begin{bmatrix}-\Gamma_{i}(t)&\Gamma_{o}(t)\\\
\Gamma_{i}(t)&-\Gamma_{o}(t)\end{bmatrix},$ (11)

where $\Gamma_{i}(t)$ and $\Gamma_{o}(t)$ are the time-dependent tunneling
rates. For the results in Figs. 2, 3, and 4, the rate for tunneling in is
given by Eq. (3), while the rate for tunneling out to a good approximation is
constant. We also introduce the jump operators $\mathbf{J}_{i}(t)$ and
$\mathbf{J}_{o}(t)$, which describe the tunneling of electrons in and out of
the quantum dot,

$\mathbf{J}_{i}(t)=\begin{bmatrix}0&0\\\
\Gamma_{i}(t)&0\end{bmatrix},\quad\mathbf{J}_{o}(t)=\begin{bmatrix}0&\Gamma_{o}(t)\\\
0&0\end{bmatrix}.$ (12)

The distribution of waiting times between electrons ($x=o$) or holes ($x=i$)
can then be expressed as

$\mathcal{W}_{x}(\tau)=\int_{0}^{T}dt\frac{\langle
1|\mathbf{J}_{x}(t+\tau)\mathbf{U}_{x}(t+\tau,t)\mathbf{J}_{x}(t)|t\rangle}{\langle
1|\mathbf{J}_{x}(t)|t\rangle}P_{x}(t),$ (13)

where we have defined the time-evolution operator

$\mathbf{U}_{x}(t_{1},t_{0})=\mathcal{T}\left\\{e^{\int_{t_{0}}^{t_{1}}dt\left[\mathbf{L}(t)-\mathbf{J}_{x}(t)\right]}\right\\}$
(14)

without tunneling events of type $x=o,i$, and $\mathcal{T}$ is the time-
ordering operator. We have also defined the cyclic state at time $t$, given by
the eigenproblem $\mathbf{U}_{x}(t+T,t)|t\rangle=|t\rangle$, together with
$\langle 1|=(1,1)$. Moreover, the probability density for the time of a
tunneling event reads

$P_{x}(t)=\langle 1|\mathbf{J}_{x}(t)|t\rangle/\int_{0}^{T}dt^{\prime}\langle
1|\mathbf{J}_{x}(t^{\prime})|t^{\prime}\rangle.$ (15)

Equation (13) can be understood by considering two tunneling events that are
separated by the waiting time $\tau$ with no other events in between. We then
integrate over all possible times within the period of the drive that the
first tunneling event could have occurred.

Following a similar line of arguments, the joint distribution of subsequent
waiting times reads

$\mathcal{W}(\tau_{1},\tau_{2})=\int_{0}^{T}dt\frac{\langle
1|\mathbf{J}_{x}(t+\tau_{1}+\tau_{2})\mathbf{U}_{x}(t+\tau_{1}+\tau_{2},t+\tau_{1})\mathbf{J}_{x}(t+\tau_{1})\mathbf{U}_{x}(t+\tau_{1},t)\mathbf{J}_{x}(t)|t\rangle}{\langle
1|\mathbf{J}_{x}(t)|t\rangle}P_{x}(t).$ (16)

This expression describes three tunneling events that are separated by the
waiting times $\tau_{1}$ and $\tau_{2}$ with no other tunneling events
occuring in between. Again, we integrate over all times that the first
tunneling event could have occurred. In the main text, we evaluate Eqs.
(13,16) numerically to obtain the theoretical results in Figs. 2 and 3.
However, below we evaluate them analytically using a simple approximation that
is valid at high driving frequencies.

Fast driving. For high driving frequencies, the time-evolution operator
$\mathbf{U}_{x}(t+\tau,t)$ can be approximated by

$\mathbf{U}_{x}(t+\tau,t)\simeq
e^{\left(\mathbf{\bar{L}}-\mathbf{\bar{J}}_{x}\right)\tau},$ (17)

where $\mathbf{\bar{L}}$ and $\mathbf{\bar{J}}_{x}$ are the rate matrix and
the jump operator with the time-dependent rates $\Gamma_{i}(t)$ and
$\Gamma_{o}(t)$ replaced by their time-averages,
$\Gamma_{i}=\int_{0}^{T}\Gamma_{i}(t)dt/T$ and
$\Gamma_{o}=\int_{0}^{T}\Gamma_{o}(t)dt/T$. From Eq. (13), we then find

$\mathcal{W}_{x}(\tau)\simeq\Gamma_{\bar{x}}\frac{e^{-\Gamma_{\bar{x}}\tau}-e^{-\Gamma_{x}\tau}}{\Gamma_{x}-\Gamma_{\bar{x}}}\frac{\int_{0}^{T}\Gamma_{x}(t+\tau)\Gamma_{x}(t)dt}{\int_{0}^{T}\Gamma_{x}(t)dt}.$
(18)

Plugging in the expressions for the time-dependent rate in Eq. (3), we obtain
Eqs. (4) and (5) in the main text.

Using the same approximation as above, we find the joint waiting time
distributions as

$\mathcal{W}_{x}(\tau_{1},\tau_{2})\simeq\Gamma_{\bar{x}}^{2}\frac{e^{-\Gamma_{\bar{x}}\tau_{1}}-e^{-\Gamma_{x}\tau_{1}}}{\Gamma_{x}-\Gamma_{\bar{x}}}\frac{e^{-\Gamma_{\bar{x}}\tau_{2}}-e^{-\Gamma_{x}\tau_{2}}}{\Gamma_{x}-\Gamma_{\bar{x}}}\frac{\int_{0}^{T}\Gamma_{x}(t+\tau_{1}+\tau_{2})\Gamma_{x}(t+\tau_{1})\Gamma_{x}(t)dt}{\int_{0}^{T}\Gamma_{x}(t)dt}.$
(19)

With a constant rate for tunneling out, $\Gamma_{o}(t)=\Gamma_{o}$, and the
time-dependent rate in Eq. (3), we find

$\mathcal{W}_{e}(\tau_{1},\tau_{2})\simeq\mathcal{W}_{s}(\tau_{1})\mathcal{W}_{s}(\tau_{2}),$
(20)

and

$\mathcal{W}_{h}(\tau_{1},\tau_{2})\simeq\mathcal{W}_{s}(\tau_{1})\mathcal{W}_{s}(\tau_{2})I_{0}\left(\alpha\sqrt{3+2\cos(2\pi
f\tau_{1})+4\cos(\pi f\tau_{1})\cos(\pi
f[\tau_{1}+2\tau_{2}])}\right)/I^{3}_{0}(\alpha),$ (21)

which then lead us to Eqs. (6,7) for the corresponding correlation ratios.

## Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC 2123
QuantumFrontiers – 390837967, the State of Lower Saxony of Germany via the
Hannover School for Nanotechnology, the Jane and Aatos Erkko Foundation,
Research Council of Finland through the Finnish Centre of Excellence in
Quantum Technology (project number 352925), and the Japan Society for the
Promotion of Science through an Invitational Fellowship for Research in Japan.

## Author contributions

JCB and TW fabricated the device on a heterostructure grown by EPR. TW carried
out the measurements and JCB, AS, and TW analyzed the experimental data. FB
and CF performed the theoretical analysis of the experiment. JCB, FB, CF, and
RJH wrote the manuscript.

## References

  * Kampen [2007] N. G. V. Kampen,  _Stochastic Processes in Physics and Chemistry_ , 3rd ed. (Elsevier, 2007).
  * Horodecki  _et al._ [2009] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).
  * Blanter and Büttiker [2000] Y. M. Blanter and M. Büttiker, Shot noise in mesoscopic conductors, Phys. Rep. 336, 1 (2000).
  * Nazarov and Blanter [2009] Y. V. Nazarov and Y. M. Blanter,  _Quantum Transport: Introduction to Nanoscience_ (Cambridge University Press, 2009).
  * Levitov and Lesovik [1993] L. S. Levitov and G. B. Lesovik, Charge distribution in quantum shot noise, JETP Lett. 58, 230 (1993).
  * Levitov  _et al._ [1996] L. S. Levitov, H. Lee, and G. B. Lesovik, Electron counting statistics and coherent states of electric current, J. Math. Phys. 37, 4845 (1996).
  * Bagrets and Nazarov [2003] D. A. Bagrets and Y. V. Nazarov, Full counting statistics of charge transfer in Coulomb blockade systems, Phys. Rev. B 67, 085316 (2003).
  * Gustavsson  _et al._ [2006] S. Gustavsson, R. Leturcq, B. Simovič, R. Schleser, T. Ihn, P. Studerus, K. Ensslin, D. C. Driscoll, and A. C. Gossard, Counting Statistics of Single Electron Transport in a Quantum Dot, Phys. Rev. Lett. 96, 076605 (2006).
  * Flindt  _et al._ [2009] C. Flindt, C. Fricke, F. Hohls, T. Novotný, K. Netočný, T. Brandes, and R. J. Haug, Universal oscillations in counting statistics, Proc. Natl. Acad. Sci. 106, 10116 (2009).
  * Ubbelohde  _et al._ [2012] N. Ubbelohde, C. Fricke, C. Flindt, F. Hohls, and R. J. Haug, Measurement of finite-frequency current statistics in a single-electron transistor, Nat. Commun. 3, 612 (2012).
  * Saito and Dhar [2016] K. Saito and A. Dhar, Waiting for rare entropic fluctuations, EPL 114, 50004 (2016).
  * Singh  _et al._ [2019] S. Singh, P. Menczel, D. S. Golubev, I. M. Khaymovich, J. T. Peltonen, C. Flindt, K. Saito, E. Roldán, and J. P. Pekola, Universal First-Passage-Time Distribution of Non-Gaussian Currents, Phys. Rev. Lett. 122, 230602 (2019).
  * Landi  _et al._ [2024] G. T. Landi, M. J. Kewming, M. T. Mitchison, and P. P. Potts, Current Fluctuations in Open Quantum Systems: Bridging the Gap Between Quantum Continuous Measurements and Full Counting Statistics, PRX Quantum 5, 020201 (2024).
  * Brandes [2008] T. Brandes, Waiting times and noise in single particle transport, Ann. Phys. 520, 477 (2008).
  * Albert  _et al._ [2011] M. Albert, C. Flindt, and M. Büttiker, Distributions of Waiting Times of Dynamic Single-Electron Emitters, Phys. Rev. Lett. 107, 086805 (2011).
  * Albert  _et al._ [2012] M. Albert, G. Haack, C. Flindt, and M. Büttiker, Electron Waiting Times in Mesoscopic Conductors, Phys. Rev. Lett. 108, 186806 (2012).
  * Vyas and Singh [1988] R. Vyas and S. Singh, Waiting-time distributions in the photodetection of squeezed light, Phys. Rev. A 38, 2423 (1988).
  * Carmichael  _et al._ [1989] H. J. Carmichael, S. Singh, R. Vyas, and P. R. Rice, Photoelectron waiting times and atomic state reduction in resonance fluorescence, Phys. Rev. A 39, 1200 (1989).
  * Gorman  _et al._ [2017] S. K. Gorman, Y. He, M. G. House, J. G. Keizer, D. Keith, L. Fricke, S. J. Hile, M. A. Broome, and M. Y. Simmons, Tunneling Statistics for Analysis of Spin-Readout Fidelity, Phys. Rev. Appl. 8, 034019 (2017).
  * Jenei  _et al._ [2019] M. Jenei, E. Potanina, R. Zhao, K. Y. Tan, A. Rossi, T. Tanttu, K. W. Chan, V. Sevriuk, M. Möttönen, and A. Dzurak, Waiting time distributions in a two-level fluctuator coupled to a superconducting charge detector, Phys. Rev. Res. 1, 033163 (2019).
  * Ranni  _et al._ [2021] A. Ranni, F. Brange, E. T. Mannila, C. Flindt, and V. F. Maisi, Real-time observation of Cooper pair splitting showing strong non-local correlations, Nat. Commun. 12, 6358 (2021).
  * Brange  _et al._ [2021] F. Brange, A. Schmidt, J. C. Bayer, T. Wagner, C. Flindt, and R. J. Haug, Controlled emission time statistics of a dynamic single-electron transistor, Sci. Adv. 7, eabe0793 (2021).
  * Hänze  _et al._ [2021] M. Hänze, G. McMurtrie, S. Baumann, L. Malavolti, S. N. Coppersmith, and S. Loth, Quantum stochastic resonance of individual Fe atoms, Sci. Adv. 7, abg2616 (2021).
  * Dasenbrook  _et al._ [2015] D. Dasenbrook, P. P. Hofer, and C. Flindt, Electron waiting times in coherent conductors are correlated, Phys. Rev. B 91, 195420 (2015).
  * Hofer  _et al._ [2016] P. P. Hofer, D. Dasenbrook, and C. Flindt, Electron waiting times for the mesoscopic capacitor, Physica E 82, 3 (2016).
  * Wagner  _et al._ [2016] T. Wagner, P. Strasberg, J. C. Bayer, E. P. Rugeramigabo, T. Brandes, and R. J. Haug, Strong suppression of shot noise in a feedback-controlled single-electron transistor, Nat. Nanotech. 12, 218 (2016).
  * Bayer  _et al._ [2017] J. C. Bayer, T. Wagner, E. P. Rugeramigabo, and R. J. Haug, Charge reconfiguration in arrays of quantum dots, Phys. Rev. B 96, 235305 (2017).
  * Wagner  _et al._ [2019] T. Wagner, P. Talkner, J. C. Bayer, E. P. Rugeramigabo, P. Hänggi, and R. J. Haug, Quantum stochastic resonance in an ac-driven single-electron quantum dot, Nat. Phys 15, 330 (2019).
  * MacLean  _et al._ [2007] K. MacLean, S. Amasha, I. P. Radu, D. M. Zumbühl, M. A. Kastner, M. P. Hanson, and A. C. Gossard, Energy-Dependent Tunneling in a Quantum Dot, Phys. Rev. Lett. 98, 036802 (2007).
  * Bocquillon  _et al._ [2013] E. Bocquillon, V. Freulon, F. D. Parmentier, J. Berroir, B. Plaçais, C. Wahl, J. Rech, T. Jonckheere, T. Martin, C. Grenier, D. Ferraro, P. Degiovanni, and G. Fève, Electron quantum optics in ballistic chiral conductors, Ann. Phys. 526, 1 (2013).
  * Aharonovich  _et al._ [2016] I. Aharonovich, D. Englund, and M. Toth, Solid-state single-photon emitters, Nat. Photon. 10, 631 (2016).
  * Cohen  _et al._ [2015] J. D. Cohen, S. M. Meenehan, G. S. MacCabe, S. Gröblacher, A. H. Safavi-Naeini, F. Marsili, M. D. Shaw, and O. Painter, Phonon counting and intensity interferometry of a nanomechanical resonator, Nature 520, 522 (2015).


# Star Formation History of Two Fields in the Halo of NGC 5128

Sima T. Aghdam1 Atefeh Javadi1 Seyedazim Hashemi2,3 Jacco Th. van Loon4 Habib
Khosroshahi1 Roya H. Golshan5 Elham Saremi1,6,7 Maryam Saberi8 1School of
Astronomy, Institute for Research in Fundamental Sciences (IPM), Tehran,
19568-36613, Iran
email<EMAIL_ADDRESS>
2Department of Physics and Astronomy, University of California Riverside, CA
92521, The USA
3Department of Physics, Sharif University of Technology, Tehran, 11155-9161,
Iran
4Lennard-Jones Laboratories, Keele University, ST5 5BG, UK
5Physikalisches Institut der Universität zu Köln, Zülpicher Str. 77, 50937,
Köln, Germany
6Instituto de Astrofìsica de Canarias, C/ Vìa Làctea s/n, 38205 La Laguna,
Tenerife, Spain
7Departamento de Astrofìsica, Universidad de La Laguna, 38205 La Laguna,
Tenerife, Spain 8Rosseland Centre for Solar Physics, University of Oslo, P.O.
Box 1029, Blindern, NO-0315, Oslo, Norway

(2022)

###### Abstract

NGC 5128 galaxy is a giant elliptical galaxy located in the Centaurus group of
galaxies at $3.8$ Mpc. We aim to study the star formation history (SFH) of two
different fields of the galaxy. The northeastern field (Field 1) is located at
a distance of $18.8$ kpc, while the southern field (Field 2) is at $9.9$ kpc.
We use a photometric method that is based on identifying long period variable
(LPV) stars and asymptotic giant branch (AGB) stars, as they are strong
tracers of star formation and galaxy evolution due to their luminosity and
variability; $395$ LPVs in Field 1 and $671$ LPVs in Field 2 have been
identified. These two fields present similar SFHs, although the SF rate of
Field 2 is more enhanced. We find that the galaxy has three major star
formation episodes t $\sim$ $800$ Myr ago, t $\sim$ $3.2$ Gyr ago, and t
$\sim$ $10$ Gyr ago, where t is look-back time. The rate of star formation at
$\sim$ 800 Myr ago agrees with previous studies suggesting that the galaxy
experienced a merger around that time. Furthermore, NGC 5128 has experienced a
lower star formation rate in its recent history which could have been driven
by jet-induction star formation and multiple outbursts of AGN activity in this
galaxy, as well as a minor merger around $400$ Myr ago.

###### keywords:

stars: AGB and LPV – stars: formation – galaxies: jets – galaxies: nuclei –
galaxies: halos – galaxies: evolution – galaxies: star formation – galaxies:
individual: NGC 5128

††volume: 373††journal: Resolving the Rise and Fall of Star Formation in
Galaxies††editors: Kim W.-T. & Wong T., eds.

##  1 Introduction

There is no giant elliptical galaxy (GE) in the local group (LG), so to study
the various features of GE galaxies, we must look at the nearest elliptical
galaxy in other groups. NGC 5128 (or Centaurus A), with a distance of 3.8 Mpc
($\mu=27.87\pm 0.16$ mag; Rejkuba et al., 2004a, and $E(B-V)=0.15\pm 0.05$
mag; Rejkuba et al., 2001) provides for us a distinctive opportunity to study
the nearest GE galaxy (Harris et al., 1999; Charmandaris et al., 2000;
Rejkuba, 2004b; Rejkuba et al., 2005), located in the Centaurus group of
galaxies (Karachentsev, 2005).

We aim to find the SFH of two small fields in the halo of NGC 5128 using long
period variable stars (LPVs) to better understand the relation between the SFH
of the halo and likely recent merger. LPV stars are evolved asymptotic giant
branch stars (AGBs) (e.g. Fraser et al., 2005, 2008; Soszyński et al., 2009)
which are identified according to high luminosity ($1000$ to $60,000$
$L_{\odot}$), as well as their variability ($\sim$ 100 to 1300 days). In
addition, red supergiant (RSG) stars are massive stars with a mass up to
$\sim$ $30$ M⊙, known as a tracer of recent star formation (about $10-30$ Myr)
in the galaxy with a strong radial pulsation (Vassiliadis & Wood, 1993; Javadi
et al., 2011b; van Loon et al., 1999, 2005).

##  2 Data

We used the data that was obtained by the ISAAC near-IR imaging spectrometer
at the ESO Paranal UT1 Antu 8.2m Telescope (VLT) in two different fields in
the halo of the galaxy that were published by Rejkuba et al. (2003), which
have been separated into two fields in the halo of the galaxy and identified
as Field 1 and Field 2 by Rejkuba et al. (2001) and Rejkuba et al. (2003).
Field 1 is on the distinguished north-eastern part of the halo, at a distance
of $\sim$ 18.8 kpc from the center of the galaxy whereas Field 2 is located at
a distance of $\sim$ 9.9 kpc from the center (Rejkuba et al., 2001). The
detected evolved stars to investigate SFH are the LPV stars with periods
longer than $70$ days (from Rejkuba et al., 2003). Therefore, $395$ and $671$
LPVs are selected in Field 1 and Field 2, respectively.

##  3 Method

We used a method to study the SFH based on LPVs and it was developed by Javadi
et al. (2011b) which was explained in the previous papers in details (M33:
Javadi et al., 2011b, 2011c, 2016, 2017; LMC & SMC: Rezaei et al., 2014; NGC
147 & NGC 185: Hamedani Golshan et al., 2017; IC 1613: Hashemi et al., 2019;
Andromeda VII: Navabi et al., 2021; Andromeda I: Saremi et al., 2021). In the
first step, we should obtain mass, age, and pulsation duration of LPVs (the
duration that stars are in the LPV phase) by applying Padova evolutionary
models (Marigo et al., 2017) for various constant metallicities. Furthermore,
for studying the star formation history of the galaxy, we consider $Z=0.039$
(Woodley et al., 2010; Rejkuba et al., 2011). By assuming all mentioned
parameters, the SFRs for different bins with specified intervals in age and
mass as follows:

$\xi(t)=\frac{dn^{\prime}(t)}{\delta t}\ \frac{\int_{\rm min}^{\rm max}f_{\rm
IMF}(m)m\ dm}{\int_{m(t)}^{m(t+dt)}f_{\rm IMF}(m)\ dm}$ (1)

where m is mass, $f_{\rm IMF}(m)$ is Kroupa initial mass function (IMF)
(Kroupa, 2001), $dn^{\prime}$ is the observed LPVs in each bin, and $\delta t$
is the pulsation duration.

The statistical error bars for each bin come from the Poisson distribution:

$\sigma_{\xi(t)}=\frac{\sqrt{N}}{N}\xi(t)$ (2)

where N is the number of stars in each age bin.

Stellar evolution in the AGB phase causes that LPV stars inject dust into the
interstellar medium (ISM). Therefore, the light from the LPVs will be detected
with a fainter magnitude (Javadi et al., 2011b; van Loon et al., 1999, 2005).
To correct the circumstellar extinction, we should do a de-reddened process by
plotting the color-magnitude diagram and overlaying the theoretical isochrones
by Marigo et al. (2017).

Furthermore, since NGC 5128 is a distant galaxy, it is unlikely to detect all
of the LPVs in the galaxy. According to Rejkuba et al. (2003) study, they
simulated variable stars by using 3 crucial parameters which are the mean
magnitude of each star, period, and amplitude. By considering their
simulation, we applied the detection probability of each star to our results.
It should be noted that it is a key point in the presented method due to its
dependency on the number of detected LPVs.

##  4 Results and Discussion

By considering metallicity Z = $0.039$ and applying the probability functions,
the SFHs of Field 1 and Field 2 are presented in Fig. 1. Despite being located
in very different parts of the galaxy, these two fields have very similar
SFHs, as well as the higher rates of star formation for Field 2. There is a
similar supply of gas to stars during their evolution. In each field, the
consistent patterns of SFRs is seen that are comparable for Field 1 and Field
2 in increasing in log t(yr) $\sim$ $8.9$ (t $\sim$ $800$ Myr), log t(yr)
$\sim$ $9.5$ (t $\sim$ $3.2$ Gyr), and log t(yr) $\sim$ $10$ (t $\sim$ $10$
Gyr) where $t$ is look-back time.

Figure 1: The SFH of Field 1 and Field 2 for Z = $0.039$. The blue represents
the SFH of Field 1, while the red one indicates the SFH of Field 2.

Additionally, an increase in star formation in Field 1 and Field 2 around
$800$ Myr ago is consistent with the idea of recent merger of NGC 5128 with a
small gas-rich galaxy (Israel, 1998; Fassett & Graham, 2000; Mould et al.,
2000; Rejkuba et al., 2001). After log t(yr) $\sim$ $8.5$, the rate of recent
star formation stays constant according to the recent merger which happened
around $400$ Myr ago (Peng et al., 2002). Furthermore, NGC 5128 has
experienced AGN activity based on the central supermassive black hole (SMBH).
There is an acceptable scenario leading to starbursts and AGN activity which
is the interaction of gas-rich galaxies during mergers (Mo, Van den Bosch &
White, 2010).

## References

  * Charmandaris et al. (2000) Charmandaris V., Combes F., van der Hulst J. M., 2000, A&A, 356, L1
  * Fassett & Graham (2000) Fassett C. I., & Graham J. A., 2000, ApJ, 538, 594
  * Fraser et al. (2005) Fraser O. J., Hawley S. L., Cook K. H., Keller S. C., 2005, AJ, 129, 768
  * Fraser et al. (2008) Fraser O. J., Hawley S. L., Cook K. H., 2008, Astron. J, 136, 1242
  * Hamedani Golshan et al. (2017) Hamedani Golshan R., Javadi A., van Loon J. Th., et al. 2017, MNRAS, 466, 1764
  * Hashemi et al. (2019) Hashemi S. A., Javadi A., van Loon J. Th., 2019, MNRAS, 483, 4751
  * Harris et al. (1999) Harris G. L. H., Harris W. E., Poole G. B., 1999, Astron. J, 117, 855
  * Israel (1998) Israel F., 1998, A&A Rev. 8, 237
  * Javadi et al. (2011a) Javadi A., van Loon J. Th., Mirtorabi M. T., 2011a, MNRAS, 411, 263
  * Javadi et al. (2011b) Javadi A., van Loon J. Th., Mirtorabi M. T., 2011b, MNRAS, 414, 3394
  * Javadi et al. (2011c) Javadi A., van Loon J. Th., Mirtorabi M. T., 2011c, in: Why Galaxies Care About AGB Stars II, eds. F. Kerschbaum, T. Lebzelter & R. F. Wing, ASPC, 445, 497
  * Javadi et al. (2013) Javadi A., van Loon J. Th., Khosroshahi H., Mirtorabi M. T., 2013, MNRAS, 432, 2824
  * Javadi et al. (2016) Javadi A., van Loon J. Th., Khosroshahi H., 2016, Memorie della Società Astronomica Italiana, 87, 278
  * Javadi et al. (2017) Javadi A., van Loon J. Th., Khosroshahi H., et al., 2017, MNRAS, 464, 2103
  * Karachentsev (2005) Karachentsev I. D., 2005, Astron. J, 129, 178
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Marigo et al. (2017) Marigo P., Girardi L., Bressan A., et al. 2017, ApJ, 835, 77
  * Mo, Van den Bosch & White (2010) Mo H., Van den Bosch F., White S. D. M., 2010, Galaxy formation and evolution. Cambridge University Press
  * Mould et al. (2000) Mould J. R., Ridgewell A., Gallagher J. S., et al. 2000, ApJ, 536, 266
  * Navabi et al. (2021) Navabi M., Saremi E., Javadi A., et al. 2021, ApJ, 910, 127
  * Peng et al. (2002) Peng E. W., Ford H. C., Freeman K. C., White R. L., 2002, Astron. J, 124, 3144
  * Rejkuba et al. (2001) Rejkuba M., Minniti D., Silva D., Bedding T., 2001, A&A, 379, 781
  * Rejkuba et al. (2003) Rejkuba M., Minniti D., Silva D., 2003, A&A, 406, 75
  * Rejkuba et al. (2004a) Rejkuba M., Greggio L., Zoccali M., 2004, A&A, 415, 915
  * Rejkuba (2004b) Rejkuba M., 2004, A&A, 413, 903
  * Rejkuba et al. (2005) Rejkuba M., Greggio L., Harris W. E., et al. 2005, ApJ, 631, 262
  * Rejkuba et al. (2011) Rejkuba M., Harris W. E., Greggio L., Harris G. L. H., 2011, A&A, 526, A123
  * Rezaei et al. (2014) Rezaei Kh S., Javadi A., Khosroshahi H., van Loon J. Th., 2014, MNRAS, 445, 2214
  * Saremi et al. (2021) Saremi E., Javadi A., Navabi M., et al. 2021, ApJ, 923, 164
  * Soszyński et al. (2009) Soszyński I., Udalski A., Szymański M. K., et al. 2009, ACTA, 59, 239
  * van Loon et al. (1999) van Loon J. Th., Groenewegen M. A. T., de Koter A., 1999, A&A, 351, 559
  * van Loon et al. (2005) van Loon J. Th., Cioni M.-R. L., Zijlstra A. A., Loup C., 2005, A&A, 438, 273
  * Vassiliadis & Wood (1993) Vassiliadis E., & Wood P. R., 1993, ApJ, 413, 641
  * Woodley et al. (2010) Woodley K. A., Harris W. E., Puzia T. H., et al. 2010, ApJ, 708, 1335