$\displaystyle=\sum_{l_{1}=1}^{m_{1}}\sum_{l_{2}=1}^{m_{2}+1}...\sum_{l_{n}=1}^{m_{n}+1}e^{i\xi_{1,l_{1}}}p\pi_{l_{1},...,l_{n}}^{(0)}\prod_{k=2}^{n}z_{k,l_{k}}$
$\displaystyle\quad+\sum_{l_{2}=1}^{m_{2}+1}...\sum_{l_{n}=1}^{m_{n}+1}e^{i\xi_{1,m_{1}+1}}(1-p)\left(\sum_{l_{1}=1}^{m_{1}}\pi_{l_{1},...,l_{n}}^{(0)}\right)\prod_{k=2}^{n}z_{k,l_{k}}$
$\displaystyle\quad+\sum_{l_{2}=1}^{m_{2}+1}...\sum_{l_{n}=1}^{m_{n}+1}e^{i\xi_{1,m_{1}+1}}\pi_{m_{1}+1,l_{2},...,l_{n}}^{(0)}\prod_{k=2}^{n}z_{k,l_{k}}$
$\displaystyle=\sum_{l_{1}=1}^{m_{1}+1}...\sum_{l_{n=1}}^{m_{n}+1}e^{i\xi_{1,l_{1}}}\pi_{l_{1},...,l_{n}}^{(1)}\prod_{k=2}^{n}z_{k,l_{k}},$

where the coefficients $\pi_{l_{1},...,l_{n}}^{(1)}$ are given by

$\pi_{l_{1},...,l_{n}}^{(1)}=(1-\delta_{l_{1},m_{1}+1})p\pi_{l_{1},...,l_{n}}^{(0)}+\delta_{l_{1},m_{1}+1}\left(\pi_{m_{1}+1,l_{2},...,l_{n}}^{(0)}+(1-p)\sum_{l_{1}^{\prime}=1}^{m_{1}}\pi_{l_{1}^{\prime},l_{2},...,l_{n}}^{(0)}\right).$
(92)

Using the same splitting-sum arguments as shown in the sequence of equations
(91), using $z_{2,l_{2}}=pe^{i\xi_{2,l_{2}}}+(1-p)e^{i\xi_{2,m_{2}+1}}$, one
can find coefficients $\pi_{l_{1},...,l_{n}}^{(2)}$ such that

$\sum_{l_{1}=1}^{m_{1}+1}...\sum_{l_{n=1}}^{m_{n}+1}\pi_{l_{1},...,l_{n}}^{(0)}\prod_{k=1}^{n}z_{k,l_{k}}=\sum_{l_{1}=1}^{m_{1}+1}...\sum_{l_{n=1}}^{m_{n}+1}e^{i(\xi_{1,l_{1}}+\xi_{2,l_{2}})}\pi_{l_{1},...,l_{n}}^{(2)}\prod_{k=3}^{n}z_{k,l_{k}}.$
(93)

And repeating such arguments $n$ times one obtains

$\sum_{l_{1}=1}^{m_{1}+1}...\sum_{l_{n=1}}^{m_{n}+1}\pi_{l_{1},...,l_{n}}^{(0)}\prod_{k=1}^{n}z_{k,l_{k}}=\sum_{l_{1}=1}^{m_{1}+1}...\sum_{l_{n=1}}^{m_{n}+1}e^{i(\xi_{1,l_{1}}+...+\xi_{n,l_{n}})}\pi_{l_{1},...,l_{n}}^{(n)}.$
(94)

For sake of completeness, we provide the recursive formula for the
coefficients $\pi_{l_{1},...,l_{n}}^{(n)}$:

$\pi_{l_{1},...,l_{n}}^{(k)}:=(1-\delta_{l_{k},m_{k}+1})p\pi_{l_{1},...,l_{n}}^{(k-1)}+\delta_{l_{k},m_{k}+1}\left(\pi_{l_{1},...,l_{k-1},m_{k}+1,l_{k+1},...,l_{n}}^{(k-1)}+(1-p)\sum_{l_{k}^{\prime}=1}^{m_{k}}\pi_{l_{1},...,l_{k-1},l_{k}^{\prime},l_{k+1},...,l_{n}}^{(k-1)}\right),$
(95)

for $k\geq 1$. We conclude hence

$\varphi_{\bm{Q}}(\bm{\xi})=\left[\sum_{l_{1}=1}^{m_{1}+1}...\sum_{l_{n=1}}^{m_{n}+1}e^{i(\xi_{1,l_{1}}+...+\xi_{n,l_{n}})}\pi_{l_{1},...,l_{n}}^{(n)}\right]^{N},$
(96)

which is indeed the characteristic function of an ECM random arrangement, with
full path probabilities given by the coefficients
$\pi_{l_{1},...,l_{n}}^{(n)}$. $\blacksquare$

###  B.8 Proof of Proposition 4.2

For the expectation, we use the properties of the conditional expectation and
we obtain

$\displaystyle\mathbb{E}\left(Q_{t}(A)\right)$
$\displaystyle=\mathbb{E}\left(\mathbb{E}\left(Q_{t}(A)\ \big{|}\
\Phi\right)\right)$ (97)
$\displaystyle=\mathbb{E}\left(p\Phi_{t}(A)\right)=p\mathbb{E}\left(\Phi_{t}(A)\right),$

where we have used that $Q_{t}(A)$, conditioned on $\Phi$, is a binomial
random variable with size $\Phi_{t}(A)$ and success probability $p$ (rule ii).
The same argument can be used for the variance:

$\displaystyle\mathbb{V}ar(Q_{t}(A))$
$\displaystyle=\mathbb{E}\left(Q_{t}^{2}(A)\right)-\mathbb{E}\left(Q_{t}(A)\right)^{2}$
(98) $\displaystyle=\mathbb{E}\left(\mathbb{E}\left(Q_{t}^{2}(A)\ \big{|}\
\Phi\right)\right)-p^{2}\mathbb{E}\left(\Phi_{t}(A)\right)^{2}$
$\displaystyle=\mathbb{E}\left(p(1-p)\Phi_{t}(A)+p^{2}\Phi_{t}^{2}(A)\right)-p^{2}\mathbb{E}\left(\Phi_{t}(A)\right)^{2}$
$\displaystyle=p(1-p)\mathbb{E}\left(\Phi_{t}(A)\right)+p^{2}\mathbb{V}ar(\Phi_{t}(A)).$

Now, if we consider $A,B$ and $t,s$ such that $A\cap B=\emptyset$ or such that
$t\neq s$, we have

$\displaystyle\mathbb{C}ov(Q_{t}(A),Q_{s}(B))$
$\displaystyle=\mathbb{E}\left(Q_{t}(A)Q_{s}(B)\right)-\mathbb{E}\left(Q_{t}(A)\right)\mathbb{E}\left(Q_{s}(B)\right)$
(99) $\displaystyle=\mathbb{E}\left(\mathbb{E}\left(Q_{t}(A)Q_{s}(B)\ \big{|}\
\Phi\right)\right)-p^{2}\mathbb{E}\left(\Phi_{t}(A)\right)\mathbb{E}\left(\Phi_{s}(B)\right)$
$\displaystyle=\mathbb{E}\left(\mathbb{E}\left(Q_{t}(A)\ \big{|}\
\Phi\right)\mathbb{E}\left(Q_{s}(B)\ \big{|}\
\Phi\right)\right)-p^{2}\mathbb{E}\left(\Phi_{t}(A)\right)\mathbb{E}\left(\Phi_{s}(B)\right)$
$\displaystyle=\mathbb{E}\left(p\Phi_{t}(A)p\Phi_{s}(B)\right)-p^{2}\mathbb{E}\left(\Phi_{t}(A)\right)\mathbb{E}\left(\Phi_{s}(B)\right)$
$\displaystyle=p^{2}\mathbb{C}ov(\Phi_{t}(A),\Phi_{s}(B)),$

where we have used that $Q_{t}(A)$ and $Q_{s}(B)$ are independent binomial
random variables conditioned on $\Phi$, according to rules i and ii. This
proves the desired result. $\blacksquare$

###  B.9 Proof of Theorem 5.1

Suppose $\tilde{X}$ satisfies conditions a and b and that there exists an
extended-real random variable $T_{c}$ satisfying the required axioms. Let
$\mu_{T_{c}}$ be its probability law over $(t_{0},\infty]$. Let $t>t_{0}$ and
$\Delta t>0$. Then,

$\displaystyle\mu_{T_{c}}\left([t,t+\Delta t]\right)$
$\displaystyle=\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\right)$ (100)
$\displaystyle=\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\ ,\ \tilde{X}(t)\in
D_{c}\right)+\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\ ,\ \tilde{X}(t)\notin
D_{c}\right).$

Let us study each term of this sum. For the first, we have

$\displaystyle\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\ ,\ \tilde{X}(t)\in
D_{c}\right)$ $\displaystyle=\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\ ,\
T_{c}\geq t\ ,\ \tilde{X}(t)\in D_{c}\right)$ (101)
$\displaystyle=\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\ \big{|}\ T_{c}\geq t\
,\ \tilde{X}(t)\in D_{c}\right)\mathbb{P}\left(T_{c}\geq t\ ,\ \tilde{X}(t)\in
D_{c}\right).$

Axiom ii implies in particular that $\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\
\big{|}\ T_{c}\geq t\ ,\ \tilde{X}(t)\in D_{c}\right)\to 0$ as $\Delta t\to
0$. Thus, expression (101) also goes to $0$ as $\Delta t\to 0$. On the other
hand we have

$\displaystyle\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\ ,\ \tilde{X}(t)\notin
D_{c}\right)$ $\displaystyle=\int_{[t,t+\Delta
t]}\mathbb{P}\left(\tilde{X}(t)\in D_{c}^{c}\ |\
T_{c}=u\right)d\mu_{T_{c}}(u)$ (102) $\displaystyle=\int_{[t,t+\Delta
t]}\mu_{\tilde{X}(t)|\tilde{X}(u)\in D_{c}}(D_{c}^{c})d\mu_{T_{c}}(u)$ (axiom
iii) $\displaystyle=\int_{(t,t+\Delta t]}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c}^{c})d\mu_{T_{c}}(u)$
$\displaystyle(\mu_{\tilde{X}(t)|\tilde{X}(t)\in D_{c}}(D_{c}^{c})=0)$
$\displaystyle\leq\mu_{T_{c}}\left((t,t+\Delta t]\right)\to 0\ \hbox{as}\
\Delta t\to 0.$

Since (101) and (102) go to $0$ as $\Delta t\to 0$, it is concluded that
$\mu_{T_{c}}\left([t,t+\Delta t]\right)\to\mu_{T_{c}}\left(\\{t\\}\right)=0$
for every $t>t_{0}$. Thus, $\mu_{T_{c}}$ has null discrete part and $T_{c}$ is
indeed a continuous random variable.

Let us now prove the existence of the limit $\lim_{\Delta t\to
0}\frac{\mu_{T_{c}}\left([t,t+\Delta t]\right)}{\Delta t}$. For this we divide
by $\Delta t$ in equations (100), (101), and (102) to obtain

$\displaystyle\frac{\mu_{T_{c}}\left([t,t+\Delta t]\right)}{\Delta t}$
$\displaystyle=\frac{\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\ \big{|}\
T_{c}\geq t\ ,\ \tilde{X}(t)\in D_{c}\right)}{\Delta
t}\mathbb{P}\left(T_{c}\geq t\ ,\ \tilde{X}(t)\in D_{c}\right)$ (103)
$\displaystyle\quad+\frac{1}{\Delta t}\int_{[t,t+\Delta
t]}\mu_{\tilde{X}(t)|\tilde{X}(u)\in D_{c}}(D_{c}^{c})d\mu_{T_{c}}(u)$

The first term in the sum in the right-side of (103) converges to
$\alpha\mathbb{P}\left(T_{c}\geq t\ ,\ \tilde{X}(t)\in D_{c}\right)$ due to
axiom ii. We claim that the second term converges to $0$. To verify this, we
first consider that the positivity and continuity conditions a, b imply that
the functions

$t\mapsto\mu_{\tilde{X}(t)}(D_{c})\quad,\quad\quad(t,u)\mapsto\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})\quad,\quad(t,u)\mapsto\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c}^{c}),$ (104)

are all continuous over their respective open domains. Indeed, the function
$t\mapsto\mu_{\tilde{X}(t)}(D_{c})$ is the restriction of the continuous
function $(t,u)\mapsto\mu_{\tilde{X}(t),\tilde{X}(u)}(D_{c}\times D_{c})$ to
the diagonal $\\{t=u\\}$ so it is continuous; the function
$(t,u)\mapsto\mu_{\tilde{X}(t)|\tilde{X}(u)\in D_{c}}(D_{c})$ is the division
of two continuous functions with never-null denominator; and
$(t,u)\mapsto\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c}^{c})=1-\mu_{\tilde{X}(t)|\tilde{X}(u)\in D_{c}}(D_{c})$, so it is
also continuous. Now, since $\mu_{\tilde{X}(t)|\tilde{X}(t)\in
D_{c}}(D_{c}^{c})=0$, by continuity one has

$\lim_{\Delta t\to 0}\sup_{u\in[t,t+\Delta
t]}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c}^{c})=\mu_{\tilde{X}(t)|\tilde{X}(t)\in D_{c}}(D_{c}^{c})=0.$
(105)

Let $\epsilon>0$ (smaller than $1$). Suppose $\Delta t$ is small enough so

$\frac{\mathbb{P}\left(T_{c}\in[t,t+\Delta t]\ \big{|}\ T_{c}\geq t\ ,\
\tilde{X}(t)\in D_{c}\right)}{\Delta
t}\leq\alpha+\epsilon\quad\hbox{and}\quad\sup_{u\in[t,t+\Delta
t]}\mu_{\tilde{X}(t)|\tilde{X}(u)\in D_{c}}(D_{c}^{c})\leq\epsilon.$ (106)

Then, we can bound (103) through

$\displaystyle\frac{\mu_{T_{c}}\left([t,t+\Delta t]\right)}{\Delta t}$
$\displaystyle\leq(\alpha+\epsilon)\mathbb{P}\left(T_{c}\geq t\ ,\
\tilde{X}(t)\in D_{c}\right)+\frac{\mu_{T_{c}}\left([t,t+\Delta
t]\right)}{\Delta t}\sup_{u\in[t,t+\Delta t]}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c}^{c})$ (107)
$\displaystyle\leq\alpha+\epsilon+\epsilon\frac{\mu_{T_{c}}\left([t,t+\Delta
t]\right)}{\Delta t},$

and so

$\frac{\mu_{T_{c}}\left([t,t+\Delta t]\right)}{\Delta
t}\leq\frac{\alpha+\epsilon}{1-\epsilon}.$ (108)

$\frac{\mu_{T_{c}}\left([t,t+\Delta t]\right)}{\Delta t}$ is thus bounded as
$\Delta t\to 0$. We obtain thus,

$\frac{1}{\Delta t}\int_{[t,t+\Delta t]}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c}^{c})d\mu_{T_{c}}(u)\leq\frac{\mu_{T_{c}}\left([t,t+\Delta
t]\right)}{\Delta t}\sup_{u\in[t,t+\Delta t]}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c}^{c})\to 0.$ (109)

Doing $\Delta t\to 0$ in (103) we conclude

$\lim_{\Delta t\to 0}\frac{\mu_{T_{c}}\left([t,t+\Delta t]\right)}{\Delta
t}=\alpha\mathbb{P}\left(T_{c}\geq t\ ,\ \tilde{X}(t)\in D_{c}\right).$ (110)

This implies that the cumulative probability function of $T_{c}$ has a right-
derivative at every $t>t_{0}$, which we call $f_{T_{c}}$:

$f_{T_{c}}(t):=\lim_{\Delta t\to 0}\frac{\mu_{T_{c}}\left([t,t+\Delta
t]\right)}{\Delta t}=\alpha\mathbb{P}\left(T_{c}\geq t\ ,\ \tilde{X}(t)\in
D_{c}\right).$ (111)

We can develop expression (111) using the axiom iii and basic probability
computations

$\displaystyle f_{T_{c}}(t)$
$\displaystyle=\alpha\left(\mathbb{P}\left(\tilde{X}(t)\in
D_{c}\right)-\mathbb{P}\left(\tilde{X}(t)\in D_{c},T_{c}<t\right)\right)$
(112)
$\displaystyle=\alpha\left(\mu_{\tilde{X}(t)}(D_{c})-\int_{(t_{0},t)}\mathbb{P}\left(\tilde{X}(t)\in
D_{c}\ \big{|}\ T_{c}=u\right)d\mu_{T_{c}}(u)\right)$
$\displaystyle=\alpha\left(\mu_{\tilde{X}(t)}(D_{c})-\int_{(t_{0},t)}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})d\mu_{T_{c}}(u)\right).$

Now we claim that $f_{T_{c}}$ is continuous at every $t>t_{0}$. Since the
function $t\mapsto\mu_{\tilde{X}(t)}(D_{c})$ is continuous, we only need to
prove the continuity of
$t\mapsto\int_{(t_{0},t)}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})d\mu_{T_{c}}(u)$. This can be done using a typical application of
Lebesgue Dominated Convergence Theorem, which applies since the integrated is
bounded and the measure is finite. Using the continuity of
$(t,u)\mapsto\mu_{\tilde{X}(t)|\tilde{X}(u)\in D_{c}}(D_{c})$, we have

$\lim_{t_{n}\to
t^{-}}\int_{(t_{0},t_{n})}\mu_{\tilde{X}(t_{n})|\tilde{X}(u)\in
D_{c}}(D_{c})d\mu_{T_{c}}(u)=\int_{(t_{0},t)}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})d\mu_{T_{c}}(u),$ (113)

which proves the left-continuity. On the other hand

$\displaystyle\lim_{t_{n}\to
t^{+}}\int_{(t_{0},t_{n})}\mu_{\tilde{X}(t_{n})|\tilde{X}(u)\in
D_{c}}(D_{c})d\mu_{T_{c}}(u)$
$\displaystyle=\int_{(t_{0},t]}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})d\mu_{T_{c}}(u)$ (114)
$\displaystyle=\int_{(t_{0},t)}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})d\mu_{T_{c}}(u)+\underbrace{\mu_{\tilde{X}(t)|\tilde{X}(t)\in
D_{c}}(D_{c})}_{=1}\underbrace{\mu_{T_{c}}(\\{t\\})}_{=0}$
$\displaystyle=\int_{(t_{0},t)}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})d\mu_{T_{c}}(u),$

where we have used that $\mu_{T_{c}}$ has no atoms. This proves the right-
continuity of $f_{T_{c}}$ and thus $f_{T_{c}}$ is continuous over
$(t_{0},\infty)$.

Finally, since the cumulative distribution function of $T_{c}$,
$F_{T_{c}}(t):=\mu_{T_{c}}((t_{0},t])$ is continuous ($\mu_{T_{c}}$ has no
atoms) and has continuous right-derivative, it is differentiable555See for
example the Math-Exchange discussion
https://math.stackexchange.com/questions/418737/continuous-right-derivative-
implies-differentiability on the subject, consulted for the last time on July
the $16$th $2024$.. The measure $\mu_{T_{c}}$ has thus a density over
$(t_{0},\infty)$ which coincides with $f_{T_{c}}=F_{T_{c}}^{\prime}$. From
(112), $f_{T_{c}}$ satisfies the integral equation

$f_{T_{c}}(t)=\alpha\left(\mu_{\tilde{X}(t)}(D_{c})-\int_{t_{0}}^{t}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})f_{T_{c}}(u)du\right).$ (115)

The continuity of the functions $u\mapsto\mu_{\tilde{X}(u)}(D_{c})$ and
$(t,u)\mapsto\mu_{\tilde{X}(t)|\tilde{X}(u)\in Dc}(D_{c})$ guarantees the
existence and the uniqueness of a continuous function solution to (115) (well
known result on Volterra integral equations, see (Volterra, 1913; Brunner,
2017)). It follows that such solution must coincide with $f_{T_{c}}$ which is
then the only possible density for the probability distribution of $T_{c}$
over $(t_{0},\infty)$. $\blacksquare$

###  B.10 Proof of Proposition 5.1

We express the problem as a fixed-point problem involving a contraction.
Consider the following Lemma.

######  Lemma B.2.

Let $E$ be a real Banach space with norm $\|\cdot\|_{E}$. Let
$\mathcal{L}:E\to E$ be a contraction mapping with constant $r<1$. Let
$\tilde{E}\subset E$ be a non-empty closed subset of $E$. Suppose in addition
that $\mathcal{L}(\tilde{E})\subset\tilde{E}$. Then, the unique solution of
the fixed-point problem $\mathcal{L}(x)=x$ is in $\tilde{E}$.

Proof of Lemma B.2: This is trivial following the same principles as in Banach
fixed-point Theorem (Burkill Burkill, 2002, Theorem 3.46): take
$x_{0}\in\tilde{E}$, and define the sequence $x_{n}=\mathcal{L}^{(n)}(x_{0})$.
This sequence is contained in $\tilde{E}$ and it converges to the unique
solution of the fixed-point problem $x=\mathcal{L}(x)$. Since $\tilde{E}$ is
closed, $x\in\tilde{E}$. $\blacksquare$

Proof of Proposition 5.1: Consider the Banach space of continuous functions
over $[t_{0},t_{H}]$, $E=C([t_{0},t_{H}])$ endowed with the supremum norm
$\|\cdot\|_{\infty}$. Consider the operator $\mathcal{L}:E\to E$ defined by

$\mathcal{L}(f)(t)=\alpha\mu_{\tilde{X}(t)}(D_{c})-\alpha\int_{t_{0}}^{t}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})f(u)du.$ (116)

Continuity and positivity conditions a, b allow to ensure that
$\mathcal{L}(f)$ is indeed a continuous function over $[t_{0},t_{H}]$. The
Volterra equation (36) can be thus re-written as the fixed-point problem

$f=\mathcal{L}(f).$ (117)

Since we know that Volterra equation (36) has a unique continuous solution,
the only task is to verify if such solution is indeed positive with total
integral over $[t_{0},t_{H}]$ smaller than $1$. For that, we consider the
subset of $E$:

$\tilde{E}=\\{\ f\in C([t_{0},t_{H}])\ \big{|}\ 0\leq
f\leq\alpha\mu_{\tilde{X}(\cdot)}(D_{c})\\}.$ (118)

Since the function $t\mapsto\mu_{\tilde{X}(t)}(D_{c})$ is strictly positive
over $(t_{0},t_{H}]$, then $\tilde{E}$ is a non-empty set (the same function
$\alpha\mu_{\tilde{X}(\cdot)}(D_{c})$ belongs to it). It is trivial to verify
that $\tilde{E}$ is closed under the supremum norm.

Now, let us verify $\mathcal{L}$ is a contraction. Set
$r:=\alpha(t_{H}-t_{0})<1$. Let $f_{1},f_{2}\in E$. Then,

$\displaystyle|\mathcal{L}(f_{1}-f_{2})(t)|$
$\displaystyle=\left|-\alpha\int_{t_{0}}^{t}\underbrace{\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})}_{\leq 1}(f_{1}-f_{2})(u)du\right|$ (119)
$\displaystyle\leq\alpha(t-t_{0})\|f_{1}-f_{2}\|_{\infty}$ $\displaystyle\leq
r\|f\|_{\infty},$ $\displaystyle\forall t\in[t_{0},t_{H}].$

This proves $\mathcal{L}$ is indeed a contraction with constant $r<1$. Now,
let $f\in\tilde{E}$. Since $f\geq 0$, it is immediate that
$\mathcal{L}(f)(t)\leq\alpha\mu_{\tilde{X}(t)}(D_{c})$ for every
$t\in[t_{0},t_{H}]$. On the other hand,

$\displaystyle\mathcal{L}(f)(t)$
$\displaystyle=\alpha\mu_{\tilde{X}(t)}(D_{c})-\alpha\int_{t_{0}}^{t}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})f(u)du$ (120)
$\displaystyle\geq\alpha\mu_{\tilde{X}(t)}(D_{c})-\alpha\int_{t_{0}}^{t}\mu_{\tilde{X}(t)|\tilde{X}(u)\in
D_{c}}(D_{c})\alpha\mu_{\tilde{X}(u)}(D_{c})du$
$\displaystyle=\alpha\mu_{\tilde{X}(t)}(D_{c})-\alpha^{2}\int_{t_{0}}^{t}\mu_{\tilde{X}(t),\tilde{X}(u)}(D_{c}\times
D_{c})du$
$\displaystyle\geq\alpha\mu_{\tilde{X}(t)}(D_{c})-\alpha^{2}\mu_{\tilde{X}(t)}(D_{c})(t-t_{0})$
$\displaystyle\geq\alpha\mu_{\tilde{X}(t)}(D_{c})-\alpha\mu_{\tilde{X}(t)}(D_{c})r$
$\displaystyle=\alpha\mu_{\tilde{X}(t)}(D_{c})(1-r)\geq 0.$

Thus $\mathcal{L}(f)\geq 0$ for every $f\in\tilde{E}$. We conclude that
$\mathcal{L}(\tilde{E})\subset\tilde{E}$, and from Lemma B.2, the solution
$f_{T_{c}}$ of the Volterra fixed-point problem (117) must be in $\tilde{E}$.
Finally, since $f_{T_{c}}\in\tilde{E}$, we have $f_{T_{c}}\geq 0$ and also

$\int_{t_{0}}^{t_{h}}f_{T_{c}}(t)dt\leq\int_{t_{0}}^{t_{h}}\alpha\mu_{\tilde{X}(t)}(D_{c})dt\leq\alpha(t_{H}-t_{0})=r<1.\quad\blacksquare$
(121)

###  B.11 Proof of Proposition 5.2

We verify that the distribution of the triplet
$\tilde{X},T_{c}^{(1)},\tilde{T}_{c}^{(2)}$ is correctly defined
hierarchically. The distribution of $T_{c}^{(1)}$ is determined by rule i.
Rule ii and the condition $\tilde{T}_{c}^{(2)}=\infty$ if $T_{c}^{(1)}>t_{L}$
fully determine the distribution of $\tilde{T}_{c}^{(2)}$ given $T_{c}^{(1)}$:

$\mu_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}=\mu_{T_{c}|T_{c}>u,\tilde{X}(u)\in
D_{c}}\mathbbm{1}_{(t_{0},t_{L}]}(u)+\delta_{\infty}\mathbbm{1}_{(t_{L},\infty]}(u).$
(122)

Let $u\in(t_{0},\infty)$ and let $B\subset(t_{0},t_{L}]$. Then,

$\displaystyle\mu_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(B)$
$\displaystyle=\mu_{T_{c}|T_{c}>u,\tilde{X}(u)\in D_{c}}(B)$ (123)
$\displaystyle=\frac{\mathbb{P}\left(T_{c}\in B\cap(u,\infty),\tilde{X}(u)\in
D_{c}\right)}{\mathbb{P}\left(T_{c}>u,\tilde{X}(u)\in D_{c}\right)}$
$\displaystyle=\frac{\int_{B}\mathbbm{1}_{(u,\infty)}(v)\mathbb{P}\left(\tilde{X}(u)\in
D_{c}\ |\ T_{c}=v\right)f_{T_{c}}(v)dv}{\mathbb{P}\left(\tilde{X}(u)\in
D_{c}\right)-\int_{t_{0}}^{u}\mathbb{P}\left(\tilde{X}(u)\in D_{c}\ |\
T_{c}=v\right)dv}$ (total probabilities)
$\displaystyle=\frac{\int_{B}\mathbbm{1}_{(u,\infty)}(v)\mu_{\tilde{X}(u)|\tilde{X}(v)\in
D_{c}}(D_{c})f_{T_{c}}(v)dv}{\underbrace{\mu_{\tilde{X}(u)}(D_{c})-\int_{t_{0}}^{u}\mu_{\tilde{X}(u)|\tilde{X}(v)\in
D_{c}}(D_{c})dv}_{=\frac{1}{\alpha}f_{T_{c}}(u)}}$
$\displaystyle=\int_{B}\mathbbm{1}_{(u,\infty)}(v)\alpha\mu_{\tilde{X}(u)|\tilde{X}(v)\in
D_{c}}(D_{c})\frac{f_{T_{c}}(v)}{f_{T_{c}}(u)}dv.$

It follows that the conditional distribution
$\mu_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}$ for $u\in(t_{0},t_{L}]$ has a
density over $(t_{0},t_{L})$ given by

$f_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(u)=\alpha\mu_{\tilde{X}(u)|\tilde{X}(v)\in
D_{c}}(D_{c})\frac{f_{T_{c}}(v)}{f_{T_{c}}(u)}\mathbbm{1}_{(u,\infty)}(v)=\alpha\mu_{\tilde{X}(v)|\tilde{X}(u)\in
D_{c}}(D_{c})\frac{\phi_{T_{c}}(v)}{\phi_{T_{c}}(u)}\mathbbm{1}_{(u,\infty)}(v),$
(124)

which proves expression (50).

Concerning $\tilde{X}$ given $(T_{c}^{(1)},\tilde{T}_{c}^{(2)})=(u,v)$, the
rule iii equation (48) covers the case $u\leq t_{L},u<v<\infty$. The cases
$u>t_{L}$ and $u\leq t_{L},v=\infty$ are still left to be made precise.

Let $t_{1},...,t_{n}\in(t_{0},\infty)$ be arbitrary time-points. Consider the
case $u>t_{L}$. Since $T_{c}^{(1)}>t_{L}$ implies
$\tilde{T}_{c}^{(2)}=\infty$, one must have

$\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}=\infty}=\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=u},\quad\forall
u>t_{L}.$ (125)

Now, rule iii Eq. (47) provides the case $u<\infty$. For the case $u=\infty$
we assume $\mathbb{P}(T_{c}^{(1)}=\infty)=\mu_{T_{c}}(\\{\infty\\})>0$, since
if not the case $T_{c}^{(1)}=\infty$ does not worth to be considered in the
conditioning. Then,

$\mu_{T_{c}}(\\{\infty\\})\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=\infty}=\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})}-\int_{t_{0}}^{\infty}\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=u}f_{T_{c}}(u)du,$
(126)

so we conclude

$\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}=\infty}=\begin{cases}\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|\tilde{X}(u)\in
D_{c}}&\hbox{ if }u<\infty,\\\
\frac{\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})}-\int_{t_{0}}^{\infty}\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|\tilde{X}(u)\in
D_{c}}f_{T_{c}}(u)du}{\mu_{T_{c}}(\\{\infty\\})}&\hbox{ if
}u=\infty.\end{cases}$ (127)

Now, for the case $u\leq t_{L}$, $v=\infty$, consider $B\subset(t_{0},t_{L}]$.
Then,

$\displaystyle\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n}),T_{c}^{(1)},\tilde{T}_{c}^{(2)}}(\
\cdot\ \times B\times\\{\infty\\})$
$\displaystyle=\int_{B}\int_{\\{\infty\\}}\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}=v}d\mu_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(v)f_{T_{c}}(u)du$
(128)
$\displaystyle=\int_{B}\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}=\infty}\mu_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(\\{\infty\\})f_{T_{c}}(u)du.$

On the other hand, using complementary sets, we have

$\displaystyle\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n}),T_{c}^{(1)},\tilde{T}_{c}^{(2)}}(\
\cdot\ \times B\times\\{\infty\\})$ (129)
$\displaystyle=\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n}),T_{c}^{(1)}}(\
\cdot\ \times
B)-\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n}),T_{c}^{(1)},\tilde{T}_{c}^{(2)}}(\
\cdot\ \times B\times(t_{0},\infty))$
$\displaystyle=\int_{B}\left[\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=u}-\int_{u}^{\infty}\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}=v}f_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(v)dv\right]f_{T_{c}(u)}du$
$\displaystyle=\int_{B}\left[\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|\tilde{X}(u)\in
D_{c}}-\int_{u}^{\infty}\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|\tilde{X}(u)\in
D_{c},\tilde{X}(v)\in
D_{c}}f_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(v)dv\right]f_{T_{c}(u)}du.$

Since (128) and (129) must hold for every $B\subset(t_{0},t_{L}]$, we obtain

$\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}=\infty}=\frac{\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|\tilde{X}(u)\in
D_{c}}-\int_{u}^{\infty}\mu_{\tilde{X}(t_{1}),...,\tilde{X}(t_{n})|\tilde{X}(u)\in
D_{c},\tilde{X}(v)\in
D_{c}}f_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(v)dv}{\mu_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(\\{\infty\\})}$
(130)

for every $u\in(t_{0},t_{L}]$ such that
$\mu_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(\\{\infty\\})>0$. This proves that
the conditional distribution of $\tilde{X}$ given
$T_{c}^{(1)},\tilde{T}_{c}^{(2)}$ can be fully specified given the rules i,ii,
iii. $\blacksquare$

###  B.12 Proof of Proposition 5.4

We consider just the case $N=1$, the other cases follow immediately by
multiplying correspondingly by $N$ considering the independence assumption on
the trajectories (see Eq. (15)). With $N=1$, the variable $Q_{t}(A)$ is a
binary random variable which equals $1$ when the individual is captured in $A$
before the time $t$:

$\\{Q_{t}(A)=1\\}=\\{\tilde{X}(T_{c}^{(1)})\in A,T_{c}^{(1)}\leq t\\}.$ (131)

Thus, its expectation for $t\leq t_{L}$ is given by

$\displaystyle\mathbb{E}(Q_{t}(A))$ $\displaystyle=\mathbb{P}(Q_{t}(A)=1)$
(132) $\displaystyle=\mathbb{P}(T_{c}^{(1)}\leq
t_{L},\tilde{X}(T_{c}^{(1)})\in A)$
$\displaystyle=\int_{t_{0}}^{t}\mathbb{P}(\tilde{X}(T_{c}^{(1)})\in
A|T_{c}^{(1)}=u)f_{T_{c}}(u)du$
$\displaystyle=\int_{t_{0}}^{t}\mu_{\tilde{X}(u)|\tilde{X}(u)\in
D_{c}}(A)f_{T_{c}}(u)du$
$\displaystyle=\int_{t_{0}}^{t}\mu_{\tilde{X}(u)}(A)\phi_{T_{c}}(u)du,\quad
t\leq t_{L}.$

If $t>t_{L}$, $Q_{t}(A)=1$ when the individual is either captured in $A$ at a
first time $T_{c}^{(1)}\in(t_{L},t]$ (no capture before $t_{L}$), or if it has
been captured for a second time at a time $T_{c}^{(2)}\in(t_{L},t]$ at a point
in $A$. Equivalently,

$\\{Q_{t}(A)=1\\}=\\{\tilde{X}(\tilde{T}_{c}^{(2)})\in A,T_{c}^{(1)}\leq
t_{L},\tilde{T}_{c}^{(2)}\in(T_{c}^{(1)},T_{c}^{(1)}+t-t_{L}]\\}\cup\\{\tilde{X}(T_{c}^{(1)})\in
A,T_{c}^{(1)}\in(t_{L},t]\\}.$ (133)

Hence,

$\displaystyle\mathbb{E}(Q_{t}(A))$
$\displaystyle=\mathbb{P}\left(\tilde{X}(\tilde{T}_{c}^{(2)})\in
A,T_{c}^{(1)}\leq
t_{L},\tilde{T}_{c}^{(2)}\in(T_{c}^{(1)},T_{c}^{(1)}+t-t_{L}]\right)+\mathbb{P}\left(\tilde{X}(T_{c}^{(1)})\in
A,T_{c}^{(1)}>t_{L}\right)$ (134)
$\displaystyle=\int_{t_{0}}^{t_{L}}\int_{u}^{u+t-t_{L}}\mu_{\tilde{X}(v)|T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}=v}(A)f_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(v)dvf_{T_{c}}(u)du+\int_{t_{L}}^{t}\mu_{\tilde{X}(u)|T_{c}^{(1)}=u}(A)f_{T_{c}}(u)du$
$\displaystyle=\int_{t_{0}}^{t_{L}}\int_{u}^{u+t-t_{L}}\mu_{\tilde{X}(v)|\tilde{X}(u)\in
D_{c},\tilde{X}(v)\in D_{c}}(A)\alpha\mu_{\tilde{X}(v)|\tilde{X}(u)\in
D_{c}}(D_{c})\frac{\phi_{T_{c}}(v)}{\phi_{T_{c}}(u)}dv\mu_{\tilde{X}(u)}(D_{c})\phi_{T_{c}}(u)du$
$\displaystyle\quad+\int_{t_{L}}^{t}\mu_{\tilde{X}(u)|\tilde{X}(u)\in
D_{c}}(A)f_{T_{c}}(u)du$
$\displaystyle=\alpha\int_{t_{0}}^{t_{L}}\int_{u}^{u+t-t_{L}}\mu_{\tilde{X}(u),\tilde{X}(v)}(D_{c}\times
A)\phi_{T_{c}}(v)dvdu+\int_{t_{L}}^{t}\mu_{\tilde{X}(u)}(A)\phi_{T_{c}}(u)du,\quad
t>t_{L}.$

For two times $t,s$ and two subsets $A,B\subset D_{c}$, the product
$Q_{t}(A)Q_{s}(B)$ is also a binary random variable which equals $1$ only if
both $Q_{t}(A)$ and $Q_{s}(B)$ equal one. We apply thus the same idea
following the equivalences (131) and (133). If $t,s\leq t_{L}$ then

$\displaystyle\mathbb{E}\left(Q_{t}(A)Q_{s}(B)\right)$
$\displaystyle=\mathbb{P}\left(Q_{t}(A)Q_{s}(B)=1\right)$ (135)
$\displaystyle=\mathbb{P}\left(T_{c}^{(1)}\leq t\wedge
s,\tilde{X}(T_{c}^{(1)})\in A\cap B\right)$
$\displaystyle=\int_{t_{0}}^{t\wedge s}\mu_{\tilde{X}(u)}(A\cap
B)\phi_{T_{c}}(u)du.$

If $t,s>t_{L}$, we have

$\displaystyle\mathbb{E}\left(Q_{t}(A)Q_{s}(B)\right)$
$\displaystyle=\mathbb{P}\left(\tilde{X}(\tilde{T}_{c}^{(2)})\in A\cap
B,T_{c}^{(1)}\leq t_{L},\tilde{T}_{c}^{(2)}\in(T_{c}^{(1)},T_{c}^{(1)}+t\wedge
s-t_{L}]\right)$ (136)
$\displaystyle\quad+\mathbb{P}\left(\tilde{X}(T_{c}^{(1)})\in A\cap
B,T_{c}^{(1)}\in(t_{L},t\wedge s]\right)$
$\displaystyle=\int_{t_{0}}^{t_{L}}\int_{u}^{u+t\wedge
s-t_{L}}\mu_{\tilde{X}(u),\tilde{X}(v)}(D_{c}\times(A\cap
B))\phi_{T_{c}}(v)dvdu+\int_{t_{L}}^{t\wedge s}\mu_{\tilde{X}(u)}(A\cap
B)\phi_{T_{c}}(u)du.$

Finally, for $t\leq t_{L}<s$ we have

$\displaystyle\mathbb{E}\left(Q_{t}(A)Q_{s}(B)\right)$
$\displaystyle=\mathbb{P}\left(\tilde{X}(T_{c}^{(1)})\in A,T_{c}^{(1)}\leq
t,\tilde{X}(\tilde{T}_{c}^{(2)})\in
B,\tilde{T}_{c}^{(2)}\in(T_{c}^{(1)},T_{c}^{(1)}+s-t_{L}]\right)$ (137)
$\displaystyle=\int_{t_{0}}^{t}\int_{u}^{u+t\wedge
s-t_{L}}\mu_{\tilde{X}(u),\tilde{X}(v)}(A\times
B)\phi_{T_{c}}(v)dvdu.\quad\blacksquare.$

###  B.13 Proof of Proposition 6.1

For number i, we remark that the case $T_{c}^{(1)}=u>t_{L}$ implies
$T_{c}^{(2)}=\tilde{T}_{c}^{(2)}=\infty>t_{H}$. Thus (see Eq. (125) in the
Proof of Proposition 5.2),

$\mu_{\tilde{X}(T_{c}^{(1)})|T_{c}^{(1)}=u}(A)=\mu_{\tilde{X}(u)|T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}=\infty}(A)=\mu_{\tilde{X}(u)|\tilde{X}(u)\in
D_{c}}(A)=\frac{\mu_{\tilde{X}(u)}(A)}{\mu_{\tilde{X}(u)}(D_{c})},\quad
u>t_{L},A\subset D_{c}.$ (138)

Now, for the case $T_{c}^{(1)}=u\leq t_{L}$ but $T_{c}^{(2)}>t_{H}$, we need
to consider the fact that the individual has not been captured twice within
the time horizon. Let $A\subset D_{c}$. Then,

$\displaystyle\mu_{\tilde{X}(T_{c}^{(1)})|T_{c}^{(1)}=u,T_{c}^{(2)}>t_{H}}(A)=\mathbb{P}\left(\tilde{X}(T_{c}^{(1)})\in
A\ |\ T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}>u+t_{H}-t_{L}\right)$ (139)
$\displaystyle=\frac{\mathbb{P}\left(\tilde{X}(T_{c}^{(1)})\in
A,\tilde{T}_{c}^{(2)}>u+t_{H}-t_{L}\ |\
T_{c}^{(1)}=u\right)}{\mathbb{P}\left(\tilde{T}_{c}^{(2)}>u+t_{H}-t_{L}\ |\
T_{c}^{(1)}=u\right))}$
$\displaystyle=\frac{\mathbb{P}\left(\tilde{X}(T_{c}^{(1)})\in A\ |\
T_{c}^{(1)}=u\right)-\mathbb{P}\left(\tilde{X}(T_{c}^{(1)})\in
A,\tilde{T}_{c}^{(2)}\leq u+t_{H}-t_{L}\ |\
T_{c}^{(1)}=u\right)}{1-\mathbb{P}\left(\tilde{T}_{c}^{(2)}\leq u+t_{H}-t_{L}\
|\ T_{c}^{(1)}=u\right))}$
$\displaystyle=\frac{\mu_{\tilde{X}(u)|\tilde{X}(u)\in
D_{c}}(A)-\int_{u}^{u+t_{H}-t_{L}}\mu_{\tilde{X}(u)|T_{c}^{(1)}=u,\tilde{T}_{c}^{(2)}=v}(A)f_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(v)dv}{1-\int_{u}^{u+t_{H}-t_{L}}f_{\tilde{T}_{c}^{(2)}|T_{c}^{(1)}=u}(v)dv}$
$\displaystyle=\frac{\frac{\mu_{\tilde{X}(u)}(A)}{\mu_{\tilde{X}(u)}(D_{c})}-\int_{u}^{u+t_{H}-t_{L}}\mu_{\tilde{X}(u)|\tilde{X}(u)\in
D_{c},\tilde{X}(v)\in D_{c}}(A)\alpha\mu_{\tilde{X}(v)|\tilde{X}(u)\in
D_{c}}(D_{c})\frac{\phi_{T_{c}}(v)}{\phi_{T_{c}}(u)}dv}{1-\int_{u}^{u+t_{H}-t_{L}}\alpha\mu_{\tilde{X}(v)|\tilde{X}(u)\in
D_{c}}(D_{c})\frac{\phi_{T_{c}}(v)}{\phi_{T_{c}}(u)}dv}$
$\displaystyle=\frac{\mu_{\tilde{X}(u)}(A)-\frac{\alpha}{\phi_{T_{c}}(u)}\int_{u}^{u+t_{H}-t_{L}}\mu_{\tilde{X}(u),\tilde{X}(v)}(A\times
D_{c})\phi_{T_{c}}(v)dv}{\mu_{\tilde{X}(u)}(D_{c})-\frac{\alpha}{\phi_{T_{c}}(u)}\int_{u}^{u+t_{H}-t_{L}}\mu_{\tilde{X}(u),\tilde{X}(v)}(D_{c}\times
D_{c})\phi_{T_{c}}(v)dv},$

which coincides with expression (59) for $u\leq t_{L}$, for every $A\subset
D_{c}$.

For number ii we consider simply $T_{c}^{(1)}=u$, $\tilde{T}_{c}^{(2)}=v$ with
$u\leq t_{L},u<v\leq u+t_{H}-t_{L}$ and we obtain following rule iii:

$\displaystyle\mu_{\tilde{X}(T_{c}^{(1)}),\tilde{X}(T_{c}^{(2)})|\tilde{T}_{c}^{(2)}=v,T_{c}^{(1)}=u}(A\times
B)$
$\displaystyle=\mu_{\tilde{X}(u),\tilde{X}(v)|\tilde{T}_{c}^{(2)}=v,T_{c}^{(1)}=u}(A\times
B)$ (140) $\displaystyle=\mu_{\tilde{X}(u),\tilde{X}(v)|\tilde{X}(u)\in
D_{c},\tilde{X}(v)\in D_{c}}(A\times B)$
$\displaystyle=\frac{\mu_{\tilde{X}(u),\tilde{X}(v)}(A\times
B)}{\mu_{\tilde{X}(u),\tilde{X}(v)}(D_{c}\times D_{c})},$

for every $A,B\subset D_{c}$. $\blacksquare$

##  Appendix C Pseudo-codes for simulation algorithms

Algorithm 1 Algorithm to simulate the Snapshot model

1:$N,\theta_{X},p$, the times of snapshot takes $t_{1},...,t_{n}$, and for
every $k=1,...,n$, the disjoint regions $A_{k,1},...,A_{k,m_{k}}$ over which
the snapshots are taken at time $t_{k}$.

2:for $j=1,...,N$ do

3: Simulate the individual trajectory
$(x_{j}(t_{1}),...,x_{j}(t_{n}))\sim\mu_{X(t_{1}),...,X(t_{n})|\theta_{X}}$.

4:end for

5:Define $\bm{Q}=(\vec{Q}_{k})_{k\in\\{1,...,n\\}}$ a list of $n$ null
vectors, each vector $\vec{Q}_{k}=(Q_{k,1},...,Q_{k,m_{k}})$ being of
dimension $m_{k}$ ($\bm{Q}$ can be a null matrix if $m_{1}=...=m_{k}$).

6:for $k=1,...,n$ do

7: for $l=1,...,m_{k}$ do

8: $Q_{k,l}\leftarrow\sum_{j=1}^{N}\mathbbm{1}_{A_{k,l}}(x_{j}(t_{k}))$

9: $Q_{k,l}\leftarrow Binomial(Q_{k,l},p)$

10: end for

11:end for

12:return $\bm{Q}$

Algorithm 2 Ad-hoc simulation method of a discretized version of the Capture
model in a time horizon $(t_{0},t_{H}]$ with one liberation time $t_{L}$,
using a capture-time based simulation.

1:$N,\theta_{\tilde{X}},\alpha,t_{0},t_{L},t_{H}$, $\Delta t$, and a
constructive-target-partition of $D_{c}$, $A_{1},...,A_{m}$.

2:$n_{L}=\frac{t_{L}-t_{0}}{\Delta t}$ and $n_{AL}=\frac{t_{H}-t_{L}}{\Delta
t}$ are integers.

3:Define $t_{1},...,t_{n_{L}+n_{AL}}$ a regular discretization time grid with
$t_{k}=t_{k-1}+\Delta t$ for $k=1,...,n_{L}+n_{AL}.$

4:Define $PP1S$ a matrix of dimensions $(n_{L}+n_{AL})\times m$ such that
$PP1S_{k,l}\leftarrow\mu_{\tilde{X}(t_{k})|\theta_{\tilde{X}}}(A_{l})$.

5:Define $PP1$ a vector of dimension $n_{L}+n_{AL}$ such that
$PP1_{k}\leftarrow\mu_{\tilde{X}(t_{k})|\theta_{X}}(D_{c})=\sum_{l=1}^{m}\mu_{\tilde{X}(t_{k})|\theta_{\tilde{X}}}(A_{l}).$

6:Verify $0<PP1_{k}\leq 1$ for all $k=1,...,n_{L}+n_{AL}$.

7:Define $PP2S$ a $4$-dimensional array with dimensions
$(n_{L}+n_{AL})\times(n_{L}+n_{AL})\times m\times m$, such that
$PP2S_{k,k^{\prime},l,l^{\prime}}\leftarrow\mu_{\tilde{X}(t_{k}),\tilde{X}(t_{k^{\prime}})|\theta_{\tilde{X}}}(A_{l}\times
A_{l^{\prime}})$. It is enough to define it for $k^{\prime}\leq k$.

8:Define $PP2$ a matrix of dimensions $(n_{L}+n_{AL})\times(n_{L}+n_{AL})$
such that
$PP2_{k,k^{\prime}}\leftarrow\mu_{\tilde{X}(t_{k}),\tilde{X}(t_{k^{\prime}})|\theta_{\tilde{X}}}(D_{c}\times
D_{c})=\sum_{l=1}^{m}\sum_{l^{\prime}=1}^{m}\mu_{\tilde{X}(t_{k}),\tilde{X}(t_{k^{\prime}})|\theta_{\tilde{X}}}(A_{l}\times
A_{l^{\prime}}).$ It is enough to define it for $k^{\prime}\leq k$.

9:Verify $0\leq PP2_{k,k^{\prime}}\leq 1$ for every $k,k^{\prime}$, and
$\displaystyle\max_{k^{\prime}=1,...,n_{L}+n_{AL}}PP2_{k,k^{\prime}}\leq
PP1_{k}$ for every $k$.

10:Define $CPP$ a lower-triangular matrix of dimensions
$(n_{L}+n_{AL})\times(n_{L}+n_{AL})$ such that
$CPP_{k,k^{\prime}}\leftarrow\frac{PP2_{k,k^{\prime}}}{PP1_{k^{\prime}}}$ for
$k^{\prime}\leq k$.

11:Define the vector $\vec{f}=(f_{1},...,f_{n_{L}+n_{AL}})$ as the solution of
the triangular linear system obtained from the discretization of the Volterra
equation (36):  $\left(I+\alpha\Delta tCPP\right)\vec{f}=\alpha PP1,$ (141)
where $I$ is the identity matrix of dimension $n_{L}+n_{AL}$.

12:Define the vector of probabilities $PTc1=(\Delta tf_{1},...,\Delta
tf_{n_{L}+n_{AL}},1-\Delta t\sum_{k=1}^{n_{L}+n_{AL}}f_{k})$.

13:Define $\vec{\phi}$ a vector of dimension $n_{L}+n_{AL}$ given by
$\phi_{k}=\frac{f_{k}}{PP1_{k}}$.

14:Define $FT2$ a matrix of dimensions $n_{AL}\times n_{L}$ given by
$FT2_{k,k^{\prime}}\leftarrow\alpha
CPP_{k^{\prime}+k,k^{\prime}}\frac{\phi_{k^{\prime}+k}}{\phi_{k^{\prime}}}$.

15:Define $\bm{Q}=(Q_{k,l})_{k,l}$ as a zero matrix of dimensions
$(n_{L}+n_{AL})\times m.$

16:for $j=1,...,N$ do

17: Simulate the first-capture-time index of individual $j$, $k_{c1}^{(1)}$,
according to  $k_{c1}^{(j)}=\begin{cases}k&\hbox{with probability
}PTc1_{k},\hbox{ for }k=1,...,n_{L}+n_{AL}\\\ \infty&\hbox{with probability
}PTc1_{n_{L}+n_{AL}+1}.\end{cases}$ (142)

18: Define both the first-time-position index $l_{c1}^{(j)}$ and the second-
time-position index $l_{c2}^{(j)}$ of individual $j$ as NULL.

19: if $k_{c1}^{(j)}\leq n_{L}$ then

20: Define the vector $PTc2^{(j)}=(\Delta tFT2_{1,k_{c1}^{(j)}},...,\Delta
tFT2_{n_{AL},k_{c1}^{(j)}},1-\Delta
t\sum_{k=1}^{n_{AL}}FT2_{k,k_{c1}^{(j)}}).$

21: Simulate the speculative second-capture-time-lag index of individual $j$,
$\tilde{k}_{c2}^{(j)}$, according to
$\tilde{k}_{c2}^{(j)}=\begin{cases}\tilde{k}&\hbox{with probability
}PTc2^{(j)}_{\tilde{k}}\hbox{ for }\tilde{k}=1,...,n_{AL},\\\ \infty&\hbox{
with probability }PTc2^{(j)}_{n_{AL}+1}.\end{cases}$ (143)

22: if $\tilde{k}_{c2}^{(j)}\leq n_{AL}$ then

23: Simulate $l_{c1}^{(j)},l_{c2}^{(j)}$ according to
$(l_{c1}^{(j)},l_{c2}^{(j)})=(l^{\prime},l)\hbox{ with probability
}\frac{PP2S_{k_{c1}^{(j)}+\tilde{k}_{c2}^{(j)},k_{c1}^{(j)},l,l^{\prime}}}{PP2_{k_{c1}^{(j)}+\tilde{k}_{c2}^{(j)},k_{c1}^{(j)}}},\hbox{
for }l,l^{\prime}=1,...,m.$ (144)

24: $Q_{k_{c1}^{(j)},l_{c1}^{(j)}}\leftarrow Q_{k_{c1}^{(j)},l_{c1}^{(j)}}+1.$

25: $Q_{n_{L}+\tilde{k}_{c2}^{(j)},l_{c2}^{(j)}}\leftarrow
Q_{n_{L}+\tilde{k}_{c2}^{(j)},l_{c2}^{(j)}}+1.$

26: else

27: Simulate $l_{c1}^{(j)}$ according to  $l_{c1}^{(j)}=l\hbox{ with
probability }\frac{PP1S_{k_{c1}^{(j)},l}-\frac{\alpha\Delta
t}{\phi_{k_{c1}^{(j)}}}\sum_{\tilde{k}=1}^{n_{AL}}\sum_{l^{\prime}=1}^{m}PP2S_{k_{c1}^{(j)}+\tilde{k},k_{c1}^{(j)},l^{\prime},l}\phi_{k_{c1}^{(j)}+\tilde{k}}}{PP1_{k_{c1}^{(j)}}PTc2^{(j)}_{n_{AL}+1}},\hbox{
for }l=1,...,m.$ (145)

28: $Q_{k_{c1}^{(j)},l_{c1}^{(j)}}\leftarrow Q_{k_{c1}^{(j)},l_{c1}^{(j)}}+1.$

29: end if

30: else if $n_{L}<k_{c1}^{(j)}\leq n_{L}+n_{AL}$ then

31: Simulate $l_{c1}^{(j)}$ according to  $l_{c1}^{(j)}=l\hbox{ with
probability }\frac{PP1S_{k_{c1}^{(j)},l}}{PP1_{k_{c1}^{(j)}}},\hbox{ for
}l=1,...,m.$ (146)

32: $Q_{k_{c1}^{(j)},l_{c1}^{(j)}}\leftarrow Q_{k_{c1}^{(j)},l_{c1}^{(j)}}+1.$

33: end if

34:end for

35:for $l=1,...,m$ do

36: for $k=1,...,n_{L}+n_{AL}$ do

37: if $k\neq 1$ and $k\neq n_{L}+1$ then

38: $Q_{k,l}\leftarrow Q_{k,l}+Q_{k-1,l}$

39: end if

40: end for

41:end for

42:return $\bm{Q}$

## References

  * Andrewartha Birch (1954) andrewartha1954distributionAndrewartha, HG. Birch, LC.  1954\.  The distribution and abundance of animals The distribution and abundance of animals ( Edn 1).
  * Bachl . (2019) bachl2019inlabruBachl, FE., Lindgren, F., Borchers, DL. Illian, JB.  2019\.  inlabru: an R package for Bayesian spatial modelling from ecological survey data inlabru: an r package for bayesian spatial modelling from ecological survey data.  Methods in Ecology and Evolution106760–766.
  * Brunner (2017) brunner2017volterraBrunner, H.  2017\.  Volterra integral equations: an introduction to theory and applications Volterra integral equations: an introduction to theory and applications ( 30).  Cambridge University Press.
  * Burkill Burkill (2002) burkill2002secondBurkill, JC. Burkill, H.  2002\.  A second course in mathematical analysis A second course in mathematical analysis.  Cambridge University Press.
  * Cao . (2022) cao2022tlrmvnmvtCao, J., Genton, MG., Keyes, DE. Turkiyyah, GM.  2022\.  tlrmvnmvt: Computing high-dimensional multivariate normal and student-t probabilities with low-rank methods in r tlrmvnmvt: Computing high-dimensional multivariate normal and student-t probabilities with low-rank methods in r.  Journal of Statistical Software1011–25.
  * Chilès Delfiner (1999) chiles1999geostatisticsChilès, J. Delfiner, P.  1999\.  Geostatistics: Modeling Spatial Uncertainty Geostatistics: Modeling Spatial Uncertainty.  John Wiley & Sons.
  * Cressie Wikle (2015) cressie2015statisticsCressie, N. Wikle, CK.  2015\.  Statistics for spatio-temporal data Statistics for spatio-temporal data.  John Wiley & Sons.
  * Daskalakis . (2015) daskalakis2015structureDaskalakis, C., Kamath, G. Tzamos, C.  2015\.  On the structure, covering, and learning of poisson multinomial distributions On the structure, covering, and learning of poisson multinomial distributions.  2015 IEEE 56th annual symposium on foundations of computer science 2015 ieee 56th annual symposium on foundations of computer science ( 1203–1217).
  * Diggle . (2010) diggle2010partialDiggle, PJ., Kaimi, I. Abellana, R.  2010\.  Partial-likelihood analysis of spatio-temporal point-process data Partial-likelihood analysis of spatio-temporal point-process data.  Biometrics662347–354.
  * Doi (1976) doi1976stochasticDoi, M.  1976\.  Stochastic theory of diffusion-controlled reaction Stochastic theory of diffusion-controlled reaction.  Journal of Physics A: Mathematical and General991479.
  * Edelsparre . (2021) edelsparre2021scalingEdelsparre, AH., Hefley, TJ., Rodríguez, MA., Fitzpatrick, MJ. Sokolowski, MB.  2021\.  Scaling up: understanding movement from individual differences to population-level dispersal Scaling up: understanding movement from individual differences to population-level dispersal.  bioRxiv2021–01.
  * Eisaguirre . (2021) eisaguirre2021multistateEisaguirre, JM., Booms, TL., Barger, CP., Goddard, SD. Breed, GA.  2021\.  Multistate Ornstein–Uhlenbeck approach for practical estimation of movement and resource selection around central places Multistate Ornstein–Uhlenbeck approach for practical estimation of movement and resource selection around central places.  Methods in Ecology and Evolution123507–519.
  * Fleming . (20141) fleming2014fineFleming, CH., Calabrese, JM., Mueller, T., Olson, KA., Leimgruber, P. Fagan, WF.  20141\.  From fine-scale foraging to home ranges: a semivariance approach to identifying movement modes across spatiotemporal scales From fine-scale foraging to home ranges: a semivariance approach to identifying movement modes across spatiotemporal scales.  The American Naturalist1835E154–E167.
  * Fleming . (20142) fleming2014nonFleming, CH., Calabrese, JM., Mueller, T., Olson, KA., Leimgruber, P. Fagan, WF.  20142\.  Non-Markovian maximum likelihood estimation of autocorrelated movement processes Non-Markovian maximum likelihood estimation of autocorrelated movement processes.  Methods in Ecology and Evolution55462–472.
  * Gardiner (2009) gardiner2009stochasticGardiner, C.  2009\.  Stochastic methods Stochastic methods ( 4).  Springer Berlin Heidelberg.
  * Gardner . (2022) gardner2022integratedGardner, B., McClintock, BT., Converse, SJ. Hostetter, NJ.  2022\.  Integrated animal movement and spatial capture–recapture models: Simulation, implementation, and inference Integrated animal movement and spatial capture–recapture models: Simulation, implementation, and inference.  Ecology10310e3771.
  * Gelfand Schliep (2018) gelfand2018bayesianGelfand, AE. Schliep, EM.  2018\.  Bayesian inference and computing for spatial point patterns Bayesian inference and computing for spatial point patterns.  NSF-CBMS regional conference series in probability and statistics Nsf-cbms regional conference series in probability and statistics ( 10, i–125).
  * Glennie . (2021) glennie2021incorporatingGlennie, R., Buckland, ST., Langrock, R., Gerrodette, T., Ballance, L., Chivers, S. Scott, M.  2021\.  Incorporating animal movement into distance sampling Incorporating animal movement into distance sampling.  Journal of the American Statistical Association116533107–115.
  * Grimm Railsback (2005) grimm2005individualGrimm, V. Railsback, S.  2005\.  Individual-based modeling and ecology. Princeton, NJ: Princeton University Press. Individual-based modeling and ecology. princeton, nj: Princeton university press.  428 p.
  * Gurarie . (2017) gurarie2017correlatedGurarie, E., Fleming, CH., Fagan, WF., Laidre, KL., Hernández-Pliego, J. Ovaskainen, O.  2017\.  Correlated velocity models as a fundamental unit of animal movement: synthesis and applications Correlated velocity models as a fundamental unit of animal movement: synthesis and applications.  Movement ecology51–18.
  * Hefley Hooten (2016) hefley2016hierarchicalHefley, TJ. Hooten, MB.  2016\.  Hierarchical species distribution models Hierarchical species distribution models.  Current Landscape Ecology Reports187–97.
  * Hefley . (2017) hefley2017mechanismHefley, TJ., Hooten, MB., Russell, RE., Walsh, DP. Powell, JA.  2017\.  When mechanism matters: Bayesian forecasting using models of ecological diffusion When mechanism matters: Bayesian forecasting using models of ecological diffusion.  Ecology Letters205640–650.
  * Hooten . (2013) hooten2013computationallyHooten, MB., Garlick, MJ. Powell, JA.  2013\.  Computationally efficient statistical differential equation modeling using homogenization Computationally efficient statistical differential equation modeling using homogenization.  Journal of agricultural, biological, and environmental statistics18405–428.
  * Hooten . (2017) hooten2017animalHooten, MB., Johnson, DS., McClintock, BT. Morales, JM.  2017\.  Animal movement: statistical models for telemetry data Animal movement: statistical models for telemetry data.  CRC press.
  * Horne . (2007) horne2007analyzingHorne, JS., Garton, EO., Krone, SM. Lewis, JS.  2007\.  Analyzing animal movements using Brownian bridges Analyzing animal movements using Brownian bridges.  Ecology8892354–2363.
  * Hotelling (1927) hotelling1927differentialHotelling, H.  1927\.  Differential equations subject to error, and population estimates Differential equations subject to error, and population estimates.  Journal of the American Statistical Association22159283–314.
  * Johnson . (2013) johnson2013estimatingJohnson, DS., Hooten, MB. Kuhn, CE.  2013\.  Estimating animal resource selection from telemetry data using point process models Estimating animal resource selection from telemetry data using point process models.  Journal of Animal Ecology8261155–1164.
  * Kéry . (2016) kery2016modelingKéry, M., Royle, J. .  2016\.  Modeling static occurrence and species distributions using site-occupancy models Modeling static occurrence and species distributions using site-occupancy models.  Applied hierarchical modelling in ecology1551–630.
  * Kranstauber (2019) kranstauber2019modellingKranstauber, B.  2019\.  Modelling animal movement as Brownian bridges with covariates Modelling animal movement as Brownian bridges with covariates.  Movement ecology71–10.
  * Krebs (2009) krebs2009ecologyKrebs, C.  2009\.  Ecology: the experimental analysis of distribution and abundance. University of British Columbia, Vancouver Ecology: the experimental analysis of distribution and abundance. University of British Columbia, Vancouver.  British Columbia655.
  * Lewis . (2013) lewis2013dispersalLewis, MA., Maini, PK. Petrovskii, SV.  2013\.  Dispersal, individual movement and spatial ecology Dispersal, individual movement and spatial ecology.  Lecture notes in mathematics (mathematics bioscience series)2071.
  * Lin . (2022) lin2022poissonLin, Z., Wang, Y. Hong, Y.  2022\.  The poisson multinomial distribution and its applications in voting theory, ecological inference, and machine learning The poisson multinomial distribution and its applications in voting theory, ecological inference, and machine learning.  arXiv preprint arXiv:2201.04237.
  * Louvrier . (2020) louvrier2020mechanisticLouvrier, J., Papaix, J., Duchamp, C. Gimenez, O.  2020\.  A mechanistic–statistical species distribution model to explain and forecast wolf (Canis lupus) colonization in South-Eastern France A mechanistic–statistical species distribution model to explain and forecast wolf (canis lupus) colonization in south-eastern france.  Spatial Statistics36100428.
  * Malchow . (2024) malchow2024fittingMalchow, AK., Fandos, G., Kormann, UG., Grüebler, MU., Kéry, M., Hartig, F. Zurell, D.  2024\.  Fitting individual-based models of spatial population dynamics to long-term monitoring data Fitting individual-based models of spatial population dynamics to long-term monitoring data.  Ecological Applications344e2966.
  * Michelot . (2019) michelot2019langevinMichelot, T., Gloaguen, P., Blackwell, PG. Étienne, MP.  2019\.  The Langevin diffusion as a continuous-time model of animal movement and habitat selection The Langevin diffusion as a continuous-time model of animal movement and habitat selection.  Methods in ecology and evolution10111894–1907.
  * Newman . (2014) newman2014modellingNewman, K., Buckland, S., Morgan, BJ., King, R., Borchers, D., Cole, DJ.Thomas, L.  2014\.  Modelling population dynamics Modelling population dynamics.  Modelling Population Dynamics: Model Formulation, Fitting and Assessment using State-Space Methods. Springer New York, New York, USA169–195.
  * Okubo . (2001) okubo2001diffusionOkubo, A., Levin, SA. .  2001\.  Diffusion and ecological problems: modern perspectives Diffusion and ecological problems: modern perspectives ( 14).  Springer.
  * Roques (2013) roques2013modelesRoques, L.  2013\.  Modèles de réaction-diffusion pour l’écologie spatiale: Avec exercices dirigés Modèles de réaction-diffusion pour l’écologie spatiale: Avec exercices dirigés.
  * Roques . (2022) roques2022spatialRoques, L., Allard, D. Soubeyrand, S.  2022\.  Spatial statistics and stochastic partial differential equations: A mechanistic viewpoint Spatial statistics and stochastic partial differential equations: A mechanistic viewpoint.  Spatial Statistics50100591.
  * Royle Dorazio (2008) royle2008hierarchicalRoyle, JA. Dorazio, RM.  2008\.  Hierarchical modeling and inference in ecology: the analysis of data from populations, metapopulations and communities Hierarchical modeling and inference in ecology: the analysis of data from populations, metapopulations and communities.  Elsevier.
  * Schnoerr . (2016) schnoerr2016coxSchnoerr, D., Grima, R. Sanguinetti, G.  2016\.  Cox process representation and inference for stochastic reaction–diffusion processes Cox process representation and inference for stochastic reaction–diffusion processes.  Nature communications7111729.
  * Schoenberg . (2002) schoenberg2002pointSchoenberg, FP., Brillinger, DR. Guttorp, P.  2002\.  Point processes, spatial-temporal Point processes, spatial-temporal.  Encyclopedia of environmetrics31573–1577.
  * GA. Seber (1986) seber1986reviewSeber, GA.  1986\.  A review of estimating animal abundance A review of estimating animal abundance.  Biometrics267–292.
  * GA. Seber Schofield (2023) seber2023estimatingSeber, GA. Schofield, MR.  2023\.  Estimating Presence and Abundance of Closed Populations Estimating Presence and Abundance of Closed Populations.  Springer.
  * GAF. Seber (1982) seber1982estimationSeber, GAF.  1982\.  The estimation of animal abundance and related parameters The estimation of animal abundance and related parameters.
  * GAF. Seber Schofield (2019) seber2019captureSeber, GAF. Schofield, MR.  2019\.  Capture-recapture: Parameter estimation for open animal populations Capture-recapture: Parameter estimation for open animal populations.  Springer.
  * Sisson . (2018) sisson2018handbookSisson, SA., Fan, Y. Beaumont, M.  2018\.  Handbook of approximate Bayesian computation Handbook of approximate bayesian computation.  CRC press.
  * Skellam (1951) skellam1951randomSkellam, JG.  1951\.  Random dispersal in theoretical populations Random dispersal in theoretical populations.  Biometrika381/2196–218.
  * Soubeyrand Roques (2014) soubeyrand2014parameterSoubeyrand, S. Roques, L.  2014\.  Parameter estimation for reaction-diffusion models of biological invasions Parameter estimation for reaction-diffusion models of biological invasions.  Population ecology56427–434.
  * Sullivan . (2009) sullivan2009ebirdSullivan, BL., Wood, CL., Iliff, MJ., Bonney, RE., Fink, D. Kelling, S.  2009\.  eBird: A citizen-based bird observation network in the biological sciences eBird: A citizen-based bird observation network in the biological sciences.  Biological conservation142102282–2292.
  * Turchin (1998) turchin1998quantitativeTurchin, P.  1998\.  Quantitative analysis of movement: measuring and modeling population redistribution in animals and plants Quantitative analysis of movement: measuring and modeling population redistribution in animals and plants.
  * Vilk . (2022) vilk2022phaseVilk, O., Campos, D., Méndez, V., Lourie, E., Nathan, R. Assaf, M.  2022\.  Phase transition in a non-Markovian animal exploration model with preferential returns Phase transition in a non-Markovian animal exploration model with preferential returns.  Physical review letters12814148301.
  * Volterra (1913) volterra1913lecconsVolterra, V.  1913\.  Leçons sur les équations intégrales et les équations intégro-différentielles: Leçons professées à la Faculté des sciences de Rome en 1910 Leçons sur les équations intégrales et les équations intégro-différentielles: Leçons professées à la faculté des sciences de rome en 1910.  Gauthier-Villars.
  * Wackernagel (2003) wackernagel2003multivariateWackernagel, H.  2003\.  Multivariate Geostatistics: an introduction with applications Multivariate Geostatistics: an introduction with applications (3 ).  Springer.
  * Wikle (2003) wikle2003hierarchicalWikle, CK.  2003\.  Hierarchical Bayesian models for predicting the spread of ecological processes Hierarchical bayesian models for predicting the spread of ecological processes.  Ecology8461382–1394.
  * Williams . (2002) williams2002analysisWilliams, BK., Nichols, JD. Conroy, MJ.  2002\.  Analysis and management of animal populations Analysis and management of animal populations.  Academic press.
  * Zamberletti . (2022) zamberletti2022understandingZamberletti, P., Papaïx, J., Gabriel, E. Opitz, T.  2022\.  Understanding complex spatial dynamics from mechanistic models through spatio-temporal point processes Understanding complex spatial dynamics from mechanistic models through spatio-temporal point processes.  Ecography20225e05956.


111‡ These authors contributed to this work equally.

# Superresolving second-order correlation imaging using synthesized colored
noise speckles

Zheng Li‡ Email<EMAIL_ADDRESS>State key Laboratory for Mesoscopic
Physics, School of Physics, Peking University, Beijing 100871, China  Xiaoyu
Nie‡ Department of Physics and Astronomy, Texas A&M University, College
Station, Texas, 77843, USA  Physics Department, Xi’an Jiaotong University,
Xi’an, Shaanxi 710049, China  Fan Yang‡ Department of Physics and Astronomy,
Texas A&M University, College Station, Texas, 77843, USA  Physics Department,
Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China  Xiangpei Liu
Department of Physics and Astronomy, Texas A&M University, College Station,
Texas, 77843, USA  Hefei National Laboratory for Physics Science at the the
Microscale and Department of Modern Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China  Dongyu Liu  State key
Laboratory for Mesoscopic Physics, School of Physics, Peking University,
Beijing 100871, China  Xiaolong Dong  State key Laboratory for Mesoscopic
Physics, School of Physics, Peking University, Beijing 100871, China  Xingchen
Zhao  Department of Physics and Astronomy, Texas A&M University, College
Station, Texas, 77843, USA  Tao Peng  Correspondence to: Tao Peng, Email:
<EMAIL_ADDRESS>Department of Physics and Astronomy, Texas A&M University,
College Station, Texas, 77843, USA  M. Suhail Zubairy  Department of Physics
and Astronomy, Texas A&M University, College Station, Texas, 77843, USA
Marlan O. Scully  Department of Physics and Astronomy, Texas A&M University,
College Station, Texas, 77843, USA  Baylor Research and Innovation
Collaborative, Baylor University, Waco, 76706, USA  Princeton University,
Princeton, New Jersey 08544, USA

We present a novel method to synthesize non-trivial speckles that can enable
superresolving second-order correlation imaging. The speckles acquire a unique
anti-correlation in the spatial intensity fluctuation by introducing the blue
noise spectrum to the input light fields through amplitude modulation.
Illuminating objects with the blue noise speckle patterns can lead to a sub-
diffraction limit imaging system with a resolution more than three times
higher than first-order imaging, which is comparable to the resolving power of
ninth order correlation imaging with thermal light. Our method opens a new
route towards non-trivial speckle generation by tailoring amplitudes of the
input light fields and provides a versatile scheme for constructing
superresolving imaging and microscopy systems without invoking complicated
higher-order correlations.

##  I Introduction

The Rayleigh diffraction limit specifies the minimum separation between two
incoherent point sources that can be resolved into distinct objects
Abbe:1873aa ; rayleigh1879xxxi . Over the decades, there has been vibrant
research to develop superresolving imaging techniques that circumvent the
Rayleigh limit by using quantum optical $N$-photon correlation d2001two ;
li2008optical ; boto2000quantum , structured illumination
gustafsson2000surpassing ; zeng2014nanometer ; classen2017superresolution ,
Rabi oscillation kiffner2008resonant ; liao2010quantum , fluorescence
saturation hell1994breaking , photoswitching hess2006ultra ; rust2006sub ;
betzig2006imaging , and higher-order detection of classical light
hemmer2006quantum ; guerrieri2010sub . For the case of classical light
illumination, the resolution can be improved by a factor of $\sqrt{2}$ with
second-order speckle correlation of thermal light oh2013sub , and can be
further enhanced by spatial filtering sprigg2016super ; chen2017sub . Higher
$N$th-order correlation provides further enhancement of resolution by
$\sqrt{N}$ times zhou2012resolution ; oppel2014directional ; zhang2016high .

As a fundamental component of the modern optics toolbox, laser speckles appear
when coherent light impinges upon a scattering sample and are typically
generated by modulating a laser beam with rotating ground glass
martienssen1964coherence , or spatial light modulator (SLM)
shapiro2008computational . The amplitudes of the speckles are distributed
following the Rayleigh statistics, resulting in a negative exponential
intensity probability density function (PDF) goodman2005introduction , known
as pseudo-thermal light chen2013100 ; li2020photon . Recently, interest has
been aroused in tailoring and generating non-trivial speckle patterns, e.g.,
by phase modulation of the input light fields at the Fourier plane of SLM,
speckles exhibiting non-Rayleigh statistics can be produced and used to
optimize speckle illumination imaging waller2012phase ; takasaki2014phase ;
bromberg2014generating ; bender2018customizing ; kondakci2016sub ;
li2016generation ; zhou2017superbunching ; bender2018customizing .

This work presents a novel method to synthesize speckles that possess non-
trivial spatial correlation via amplitude modulation of the input light fields
and achieve superresolving second-order correlation imaging with the obtained
speckle illumination. The speckle patterns are generated by tailoring the
light fields’ amplitudes using a digital micromirror device (DMD). They have a
blue noise power spectrum with spectral intensity increasing for higher
spatial frequencies ulichney1988dithering ; ulichney1988dithering . In
addition, the PDF of the blue noise speckles is made to obey the Rayleigh
statistics, which is the same as the white noise speckles. The unique property
of our blue noise speckles is the negative correlation in intensity
fluctuation between neighboring spatial pixels, which endows the second-order
correlation imaging with superresolving power when the object is illuminated
by such speckles. We experimentally demonstrate that the resulting image shows
$>3$ times higher resolution than the first-order imaging system, which
essentially surpasses the $\sqrt{2}$ times enhancement by conventional second-
order imaging using thermal light.

##  II Materials and Methods

###  II.1 Synthesized colored noise speckles and the spatial correlation

Figure 1:  (a) A typical 2D blue noise speckle pattern, (b) its probability
distribution, and (c) its spectrum distribution. (d) A typical 2D white noise
speckle pattern, (e) its probability distribution, and (f) its spectrum
distribution. Both the blue and white noise speckle patterns obey the same
Rayleigh distribution. The red dashed lines in (b), (c), (e), and (f) are the
corresponding plots based on averaging over all the 3000 blue and white noise
patterns used in the experiment, respectively. The blue noise power spectrum
scales asymptotically as $\omega^{4}$ for spatial frequency $\omega$.

The noise speckle patterns are generated and imprinted into the amplitude of
the input light field by DMD. We take the spectral power distribution of the
noise speckles to be $I(\omega)\simeq C_{0}\delta(\omega)+C\omega^{n}$ for
spatial frequency $\omega$. White noise ($n=0$), and blue noise with $n=4$ are
used in our experiment. After applying random conjugated symmetric phases to
the power spectrum, an inverse Fourier transform is performed to obtain the
noise speckle patterns in grayscale. The generated patterns obey Gaussian
statistics with low spatial correlation, which makes the imaging system suffer
from measurement noise. Thus we redistribute the speckle PDF to satisfy
Rayleigh statistics with the local intensity transformation Bender:19

$\displaystyle\int^{I_{0}}_{0}P(I)dI=\int^{I^{\prime}_{0}}_{0}P(I^{\prime})dI^{\prime},$
(1)

where $P(I^{\prime})$ is the PDF of the original created Gaussian patterns.
The patterns are then one-to-one mapped to the negative Rayleigh distribution

$\displaystyle P(I)=(1/\bar{I})\exp{(-I/\bar{I})}$ (2)

The resulting Rayleigh blue noise speckle pattern, its PDF, and its power
spectrum are shown in Figs. 1(a), (b), and (c), respectively. The white noise
Rayleigh speckle pattern, its PDF, and power spectrum are also shown in Figs.
1(d), (e), and (f), respectively. Unlike previous works that obtain non-
trivial second-order correlations via tailoring the speckles’ statistics, we
show that the blue noise speckle has Rayleigh statistics resulting in the same
autocorrelation for each pixel, same as the white noise speckle. However, a
non-trivial correlation between one pixel and its neighboring pixel appears
due to the modified spectrum distribution.

Figure 2: 2D spatial intensity fluctuation correlation of speckle patterns for
an arbitrary pixel and its neighboring pixels. (a) white noise, (b) blue noise
($n=1$), and (c) blue noise ($n=4$). It is shown that the white noise has a
zero correlation between different pixels, as expected. For blue noise, there
exists a significantly negative spatial correlation between neighboring
pixels; the fluctuation correlation between adjacent pixels is $-0.21$ for
$n=1$ blue noise speckles and reaches $-0.44$ for $n=4$ blue noise speckles.

We then examine the spatial correlation of the white noise and blue noise
Rayleigh speckles. To simplify the calculation without loss of generality, we
start with one dimension. The spatial intensity fluctuation correlation, which
is central to the image formation, is defined as

$\displaystyle\gamma^{(2)}(\tilde{x})=\langle\Delta I(x_{1})\Delta
I(x_{2})\rangle=\frac{\langle I(x_{1})I(x_{2})\rangle}{\langle
I(x_{1})\rangle\langle I(x_{2})\rangle}-1,$ (3)

where $\tilde{x}=x_{2}-x_{1}$. Given the power spectrum density of the
patterns, the intensity fluctuation correlation is

$\displaystyle\gamma^{(2)}(\tilde{x})=\mathcal{F}^{-1}\\{|C\omega^{n}|^{2}\\}(\tilde{x}).$
(4)

For white noise speckles, the intensity of each pixel fluctuates randomly and
independently:

$\displaystyle\gamma^{(2)}_{\mathrm{w}}(\tilde{x})=\mathcal{F}^{-1}\\{|C_{\mathrm{w}}|^{2}\\}(\tilde{x})\sim\delta(\tilde{x}).$
(5)

There is no correlation between adjacent pixels.

For blue noise speckles, we have the intensity fluctuation correlation as

$\displaystyle\gamma^{(2)}_{\mathrm{b}}(\tilde{x})=\mathcal{F}^{-1}\\{|C_{\mathrm{b}}\omega^{4}|^{2}\\}(\tilde{x}).$
(6)

This results in a rather complicated form as compared with the white noise
case. The Rayleigh blue noise speckles possess significantly negative spatial
correlation of intensity fluctuation between neighboring pixels, with the
slope and negativity of $\gamma^{(2)}_{\mathrm{b}}(\tilde{x})$ essentially
enhanced for larger blue noise order $n$. The correlation oscillates as
distance between two pixels increases, and reaches zero when they are well
separated. We plot the 2D correlation functions of the intensity fluctuation
of white noise speckles (Fig. 2 (a)), blue noise speckles for $n=1$ (Fig.
2(b)), and $n=4$ (Fig. 2(c)), using 3000 patterns for each noise type. The
fluctuation correlation between adjacent pixels reaches $-0.44$ for $n=4$ blue
noise speckles. The striking anti-correlation of intensity fluctuation
correlation readily provides the capability to resolve objects beyond the
Rayleigh diffraction limit.

###  II.2 Measurement method of second-order correlation imaging

To demonstrate the corresponding resolving power of second-order correlation
imaging with blue noise speckle illumination, we employ the experimental setup
shown in Fig. 3.

Figure 3:  Schematic setup of the second-order correlation imaging system. A
CW laser is reflected by a DMD. The synthesized speckle patterns and digital
double-slit are loaded onto the DMD. The reflected light propagates to the
imaging plane by an imaging lens. An 800 $\mu$m diameter iris pinhole is put
right behind the lens to adjust the effective aperture of the lens. The CMOS
sensor at the image plane records the intensity distribution $I(r)$ of the
double-slit image.

An input CW laser is reflected by the DMD, on which the speckle patterns of
synthesized spatial spectrum and a digital object (double-slit pattern) are
loaded. The convoluted double-slit and speckle patterns are then imaged by the
complementary metal oxide semiconductor (CMOS) sensor through an imaging lens.
The effective aperture of the imaging lens can be adjusted by a pinhole behind
the lens. The DMD contains pixelated mirrors with an area of 16 $\mu$m
$\times$ 16 $\mu$m per pixel. The double-slit has a dimension of 4 pixels in
height, 3 pixels in width for each slit, and 5 pixels of center-to-center
separation. The noise pattern consists of 54 $\times$ 98 independent pixels
(each pixel is 4 $\times$ 4 DMD pixels). The diameter of the pinhole is chosen
to be 800 $\mu$m.

##  III Results and Discussion

###  III.1 Superresolution with blue noise speckle patterns

Due to existence of diffraction, the light reflected from each slit diffuses
and can become overlapped with each other. We measured the first-order image
$I(r_{\mathrm{i}})$ and the second-order correlation image with white and blue
noise ($n=4$) speckles illumination, which is constructed from the intensity
fluctuation correlation of measured intensities on each pixel of the image
plane. The intensity distribution of the light field on the image plane can be
expressed as

$\displaystyle\langle
I(r_{\mathrm{i}})\rangle\propto\left\langle\int\mathrm{d}r_{\mathrm{o}}I(r_{\mathrm{o}})|T(r_{\mathrm{o}})|^{2}|h(r_{\mathrm{i}},r_{\mathrm{o}})|^{2}\right\rangle,$
(7)

where $I(r_{\mathrm{o}})$ is the intensity distribution of the speckle pattern
and $r_{\mathrm{o}}$ is the coordinate on the object plane.
$T(r_{\mathrm{o}})$ denotes the object aperture function, and

$h(r_{\mathrm{i}},r_{\mathrm{o}})\propto\operatorname{somb}\left(\frac{2\pi
R|r_{\mathrm{o}}+r_{\mathrm{i}}/m|}{\lambda s_{\mathrm{o}}}\right)$ (8)

is the point spread function (PSF) shih2018introduction . Here $R$ is the
effective radius of the image lens, $m=s_{\mathrm{i}}/s_{\mathrm{o}}$ is the
magnification factor of the imaging system, $s_{\mathrm{o}}$ is the distance
between object plane and lens, $s_{\mathrm{i}}$ is distance between lens and
image plane, $\operatorname{somb}(r)=J_{1}(r)/r$, and $J_{1}(r)$ is the first-
order Bessel function of the first kind.

Since the slit width is smaller than the noise speckle size, we can treat the
two separated slits as two points in the object plane, at $r_{1}$ and $r_{2}$.
The first-order image is then

$\displaystyle\langle I(r_{\mathrm{i}})\rangle\propto
h^{2}(r_{\mathrm{i}},r_{1})+h^{2}(r_{\mathrm{i}},r_{2})\,.$ (9)

The second-order image measured by illumination of white and blue noise
speckle pattern is

$\displaystyle\langle\Delta I_{\mathrm{w/b}}^{2}(r_{\mathrm{i}})\rangle$
$\displaystyle\propto\left\langle\int\mathrm{d}r\mathrm{d}r_{\mathrm{o}}\,\gamma^{(2)}_{\mathrm{w/b}}(r,r_{\mathrm{i}})|T(r_{\mathrm{o}})|^{4}h^{4}(r,r_{\mathrm{o}})\right\rangle\,.$
(10)

In Eq. (10), the image formation depends strongly on the properties of spatial
correlation of illumination light fields with different types of speckle
patterns and cannot be expressed as the overlap of two PSFs from two separated
points $r_{1},r_{2}$. For the blue noise speckles, their significantly
negative cross correlation of the pixelwise intensity fluctuation and the
contribution of the autocorrelation from each pixel results in the
superresolution of the formed image.

Figure 4: second-order image with Rayleigh blue noise, white noise, and first-
order image of the double-slit. The image of the two slits are well separated
in the second-order image with blue noise speckle illumination (a), but are
indistinguishable in the second-order image with white noise speckle
illumination (b), and are completely overlapped in the first-order image
(c).The slight asymmetry of the double-slit image distribution in (a) and (b)
is present because the iris is not exactly oriented to the middle of the two
slits.

The experimental results are shown in Fig. 4(a) for second-order image with
blue noise illumination, Fig. 4(b) for second-order image with white noise
illumination and Fig. 4(c) for first-order image. All three measurements are
performed under the same conditions, i.e., the same laser power, number of
frames, and exposure time of each frame. Data are collected by the CMOS sensor
with exposure time at 30 $\mu$s for each pattern and averaged over 3000
patterns. It can be clearly seen that, due to the finite pinhole size, the
first-order imaging system is not able to resolve the double-slit at all,
while the second-order imaging with white noise cannot distinguish the double-
slit either and the image is blurred. On the other hand, Fig. 4(a) shows a
clearly distinguished slit image.

In Fig. 5, we plot the one-dimensional image of the double-slit when the
effective diameter of the lens is just enough to separate the two slits, i.e.,
at the Rayleigh diffraction limit. The simulation, in solid black lines, shows
that a lens with a minimal effective diameter of $\sim$ 2.1 mm is required. We
also choose the 1D cross-section of vertical pixel 55, from the blue noise
second-order image and the first-order image results from Fig. 4. The data are
fitted to the first-order imaging equation, providing a quantitative
evaluation of the resolving power of the imaging systems. Fitting the first-
order image data results in an effective lens diameter of $\sim 0.8$ mm, which
is well consistent with our experimental setup of the first-order imaging
system. The same analysis shows that, to have an equivalent resolution of the
second-order blue noise imaging system, one needs to use a lens of $\sim 2.5$
mm effective diameter in the first-order imaging system. Thus, using the blue
noise speckle in second-order correlation imaging, the resolution is enhanced
by $\sim 3$ times.

Figure 5: The resolution of the first-order and second-order blue noise
imaging systems, which is characterized by their corresponding effective lens
diameters. The solid black line represents the theoretical Rayleigh limit when
the effective diameter of the lens is just able to resolve the two slits,
which is found to be $2.1$ mm. Red dots are the measured first-order image,
fitted with an effective lens diameter of $0.84$ mm. Blue dots are the
measured second-order blue noise image, fitted to an effective lens diameter
of $2.5$ mm.

##  IV Conclusion

In summary, we proposed a novel method to synthesize blue noise speckle
patterns satisfying Rayleigh statistics by tailoring the amplitudes of the
input light fields. The resulting speckles show a unique feature of negative
spatial correlation between neighboring pixels. The anti-correlation of
spatial intensity fluctuation endows the second-order correlation imaging
system a resolving power that is $\sim 3$ times higher than the first-order
imaging and is equivalent to a $9$th order correlation imaging system using
thermal light illumination. Our method is versatile and compatible with a
broad range of optical setups. It can be used in computational ghost imaging
systems to sharpen the edges of the ghost image and enhance the signal-to-
noise ratio. Our method has potential application in microscopy and biomedical
imaging systems such as laser speckle contrast imaging (LSCI) based on the
second-order correlation imaging mechanism briers1996laser ;
aminfar2019application ; zhang2019laser . With the buttress of blue boise
speckle illumination, LSCI can obtain essential enhancement of resolution.

##  V Conflict of interest

The authors declare no conflict of interest.

##  VI Author contributions

ZL, TP, MSZ, and MOS initiated the project. ZL and TP proposed the approach,
XYN, FY, DYL, XLD and TP did the simulations. ZL, XYN, FY, XPL, XCZ, TP, MSZ
and MOS realized the structured noise generation and imaging system. All
authors prepared the manuscript.

##  VII Acknowledgements

The authors thank R. Nessler and H. Xu for helpful discussion, X. N. and F. Y.
thank A. Svidzinsky for his help during their visit to IQSE, Texas A&M
University. We acknowledge financial support by the Air Force Office of
Scientific Research (Award No. FA9550-20-1-0366 DEF), Office of Naval Research
(Award No. N00014-20-1-2184), Robert A. Welch Foundation (Grant No. A-1261),
National Science Foundation (Grant No. PHY-2013771), and Qatar National
Research Fund (project NPRP 13S-0205-200258).

## References

  * (1) Abbe, E.  Beiträge zur Theorie des Mikroskops und der mikroskopischen Wahrnehmung.  _Archiv für Mikroskopische Anatomie_ 9, 413–468 (1873).
  * (2) Rayleigh, L.  Investigations in optics, with special reference to the spectroscope.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ 8, 261–274 (1879).
  * (3) D’Angelo, M., Chekhova, M. V. & Shih, Y.  Two-photon diffraction and quantum lithography.  _Phys. Rev. Lett._ 87, 013602 (2001).
  * (4) Li, H. _et al._ Optical imaging beyond the diffraction limit via dark states.  _Phys. Rev. A_ 78, 013803 (2008).
  * (5) Boto, A. N. _et al._ Quantum interferometric optical lithography: exploiting entanglement to beat the diffraction limit.  _Phys. Rev. Lett._ 85, 2733 (2000).
  * (6) Gustafsson, M. G.  Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy.  _J. Microsc._ 198, 82–87 (2000).
  * (7) Zeng, X., Al-Amri, M. & Zubairy, M. S.  Nanometer-scale microscopy via graphene plasmons.  _Phys. Rev. B_ 90, 235418 (2014).
  * (8) Classen, A., von Zanthier, J., Scully, M. O. & Agarwal, G. S.  Superresolution via structured illumination quantum correlation microscopy.  _Optica_ 4, 580–587 (2017).
  * (9) Kiffner, M., Evers, J. & Zubairy, M.  Resonant interferometric lithography beyond the diffraction limit.  _Phys. Rev. Lett._ 100, 073602 (2008).
  * (10) Liao, Z., Al-Amri, M. & Zubairy, M. S.  Quantum lithography beyond the diffraction limit via rabi oscillations.  _Phys. Rev. Lett._ 105, 183601 (2010).
  * (11) Hell, S. W. & Wichmann, J.  Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy.  _Opt. Lett._ 19, 780–782 (1994).
  * (12) Hess, S. T., Girirajan, T. P. & Mason, M. D.  Ultra-high resolution imaging by fluorescence photoactivation localization microscopy.  _Biophys. J._ 91, 4258–4272 (2006).
  * (13) Rust, M. J., Bates, M. & Zhuang, X.  Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM).  _Nature Methods_ 3, 793–796 (2006).
  * (14) Betzig, E. _et al._ Imaging intracellular fluorescent proteins at nanometer resolution.  _Science_ 313, 1642–1645 (2006).
  * (15) Hemmer, P. R., Muthukrishnan, A., Scully, M. O. & Zubairy, M. S.  Quantum lithography with classical light.  _Phys. Rev. Lett._ 96, 163603 (2006).
  * (16) Guerrieri, F. _et al._ Sub-Rayleigh imaging via $n$-photon detection.  _Phys. Rev. Lett._ 105, 163602 (2010).
  * (17) Oh, J.-E., Cho, Y.-W., Scarcelli, G. & Kim, Y.-H.  Sub-Rayleigh imaging via speckle illumination.  _Opt. Lett._ 38, 682–684 (2013).
  * (18) Sprigg, J., Peng, T. & Shih, Y.  Super-resolution imaging using the spatial-frequency filtered intensity fluctuation correlation.  _Scientific reports_ 6, 38077 (2016).
  * (19) Chen, X.-H., Kong, F.-H., Fu, Q., Meng, S.-Y. & Wu, L.-A.  Sub-Rayleigh resolution ghost imaging by spatial low-pass filtering.  _Opt. Lett._ 42, 5290–5293 (2017).
  * (20) Zhou, Y., Liu, J., Simon, J. & Shih, Y.  Resolution enhancement of third-order thermal light ghost imaging in the photon counting regime.  _J. Opt. Soc. Am. B_ 29, 377–381 (2012).
  * (21) Oppel, S., Wiegner, R., Agarwal, G. & Von Zanthier, J.  Directional superradiant emission from statistically independent incoherent nonclassical and classical sources.  _Phys. Rev. Lett._ 113, 263606 (2014).
  * (22) Zhang, S., Wang, W., Yu, R. & Yang, X.  High-order correlation of non-Rayleigh speckle fields and its application in super-resolution imaging.  _Laser Phys._ 26, 055007 (2016).
  * (23) Martienssen, W. & Spiller, E.  Coherence and fluctuations in light beams.  _Am. J. Phys._ 32, 919–926 (1964).
  * (24) Shapiro, J. H.  Computational ghost imaging.  _Phys. Rev. A_ 78, 061802 (2008).
  * (25) Goodman, J. W.  _Introduction to Fourier optics_ (Roberts and Company Publishers, 2005).
  * (26) Chen, H., Peng, T. & Shih, Y.  100% correlation of chaotic thermal light.  _Phys. Rev. A_ 88, 023808 (2013).
  * (27) Li, S.-W., Li, F., Peng, T. & Agarwal, G.  Photon statistics of quantum light on scattering from rotating ground glass.  _Phys. Rev. A_ 101, 063806 (2020).
  * (28) Waller, L., Situ, G. & Fleischer, J. W.  Phase-space measurement and coherence synthesis of optical beams.  _Nature Photonics_ 6, 474–479 (2012).
  * (29) Takasaki, K. T. & Fleischer, J. W.  Phase-space measurement for depth-resolved memory-effect imaging.  _Opt. Express_ 22, 31426–31433 (2014).
  * (30) Bromberg, Y. & Cao, H.  Generating non-Rayleigh speckles with tailored intensity statistics.  _Phys. Rev. Lett._ 112, 213904 (2014).
  * (31) Bender, N., Yılmaz, H., Bromberg, Y. & Cao, H.  Customizing speckle intensity statistics.  _Optica_ 5, 595–600 (2018).
  * (32) Kondakci, H. E., Szameit, A., Abouraddy, A. F., Christodoulides, D. N. & Saleh, B. E.  Sub-thermal to super-thermal light statistics from a disordered lattice via deterministic control of excitation symmetry.  _Optica_ 3, 477–482 (2016).
  * (33) Li, X. _et al._ Generation of a super-Rayleigh speckle field via a spatial light modulator.  _Appl. Phys. B_ 122, 82 (2016).
  * (34) Zhou, Y. _et al._ Superbunching pseudothermal light.  _Phys. Rev. A_ 95, 053809 (2017).
  * (35) Ulichney, R. A.  Dithering with blue noise.  _Proc. IEEE_ 76, 56–79 (1988).
  * (36) Bender, N., Yilmaz, H., Bromberg, Y. & Cao, H.  Introducing non-local correlations into laser speckles.  _Opt. Express_ 27, 6057 (2019).
  * (37) Shih, Y.  _An introduction to quantum optics: photon and biphoton physics_ (CRC press, 2018).
  * (38) Briers, J. D. & Webster, S.  Laser speckle contrast analysis (lasca): a nonscanning, full-field technique for monitoring capillary blood flow.  _J. Biomed. Opt._ 1, 174–180 (1996).
  * (39) Aminfar, A., Davoodzadeh, N., Aguilar, G. & Princevac, M.  Application of optical flow algorithms to laser speckle imaging.  _Microvascular Research_ 122, 52–59 (2019).
  * (40) Zhang, R. _et al._ Laser speckle imaging for blood flow based on pixel resolved zero-padding auto-correlation coefficient distribution.  _Opt. Commun._ 439, 38–46 (2019).