# Dynamical pattern formation without self-attraction in quorum-sensing active
matter:
the interplay between nonreciprocity and motility

Yu Duan  Max Planck Institute for Dynamics and Self-Organization (MPI-DS),
37077 Göttingen, Germany  Jaime Agudo-Canalejo  Max Planck Institute for
Dynamics and Self-Organization (MPI-DS), 37077 Göttingen, Germany  Ramin
Golestanian<EMAIL_ADDRESS>Max Planck Institute for Dynamics and
Self-Organization (MPI-DS), 37077 Göttingen, Germany  Rudolf Peierls Centre
for Theoretical Physics, University of Oxford, Oxford OX1 3PU, United Kingdom
Benoît Mahault<EMAIL_ADDRESS>Max Planck Institute for Dynamics and
Self-Organization (MPI-DS), 37077 Göttingen, Germany

###### Abstract

We study a minimal model involving two species of particles interacting via
quorum-sensing rules. Combining simulations of the microscopic model and
linear stability analysis of the associated coarse-grained field theory, we
identify a mechanism for dynamical pattern formation that does not rely on the
standard route of intra-species effective attractive interactions. Instead,
our results reveal a highly dynamical phase of chasing bands induced only by
the combined effects of self-propulsion and nonreciprocity in the inter-
species couplings. Turning on self-attraction, we find that the system may
phase separate into a macroscopic domain of such chaotic chasing bands
coexisting with a dilute gas. We show that the chaotic dynamics of bands at
the interfaces of this phase-separated phase results in anomalously slow
coarsening.

Active systems are driven out-of-equilibrium at the level of their microscopic
constituents. Since activity may arise in various forms, including particle
motility [1], local generation of nonconservative forces [2, 3] and torques
[4, 5], growth [6], or sustained chemical reactions [7], the active matter
field is rapidly expanding in multiple directions [8]. Recently, the
generation of interactions breaking action-reaction symmetry [9, 10, 11, 12,
13] as a paradigm for activity has received increasing attention. Such
nonreciprocal interactions occur in a variety of contexts, including particles
interacting through a nonequilibrium medium [14, 15, 16, 17, 18, 19, 20, 21],
or via social forces [22, 23, 24].

Due to connections built with non-Hermitian physics [25, 26], nonreciprocity
has been argued to constitute a generic route for the emergence of steady
states breaking time reversal symmetry (TRS) [18, 10, 11, 27, 12, 28, 29, 30].
Examples include rotating chiral phases in flocking systems involving multiple
species [12], or traveling patterns in phase-separating mixtures [10, 11]. The
latter are described by the nonreciprocal Cahn-Hilliard (NRCH) model [10, 11],
which generalizes the Cahn-Hilliard theory of phase separation [31] by
including a nonequilibrium chemical potential with antisymmetric couplings
between species. Importantly, the NRCH model predicts that preconditions for
the emergence of TRS broken phases are the presence of intra-species
attraction that drives demixing, as well as chasing inter-species interactions
[10, 11]. In addition to dynamical patterns, the similarity between the NRCH
model and reaction-diffusion equations [32, 33] has led to additional
connections such as the presence of Turing-like instabilities [27, 34].

A model mechanism for self-organization in motile active matter are quorum-
sensing interactions. Quorum-sensing mediated by chemical signals is for
example known to drive swarming and pattern formation in bacterial populations
[35, 36, 37]. Collective aggregation was moreover shown to arise in
suspensions of light-activated colloids whose motility is locally regulated by
their density via feedback control loops [38, 23, 39]. Often, quorum-sensing
interactions are modelled as a direct response of particle motility to local
density variations [40, 41, 42]. Despite their simplicity, minimal quorum-
sensing models are able to capture complex collective behaviors [23, 37]. A
notable example is the motility induced phase separation (MIPS), whose origin
is rooted in the effective attraction induced by the self-inhibition of
motility [40]. In addition, nonreciprocal couplings naturally arise when
multiple species exhibit quorum-sensing with asymmetric responses [17, 43],
thus without resorting to explicit nonreciprocal pairwise forces [44, 45, 46].
For strong nonreciprocity and in the phase separation regime, related models
were shown to lead to TRS broken phases in qualitative agreement with the
picture provided by the NRCH model [16, 17, 43].

In this Letter, we study a binary model of quorum-sensing self-propelled
particles in two dimensions. Our simulations reveal the presence of thin
traveling bands even in the absence of intra-species couplings. Instead, and
as confirmed by the associated coarse-grained field theory, they arise from a
new mechanism relying only on self-propulsion and chasing interactions. In
large systems, these bands self-organize into a spatio-temporal chaotic phase
which we find to be involved in phase-separated configurations at moderate
self-attraction. We show that this phase-separated phase coarsens with a
scaling exponent distinct from that of standard MIPS [47, 48, 49, 50].

Figure 1:  (a-d) Representative simulation snapshots for strong (a,b) and
vanishing (c,d) motility self-inhibition in the presence of chasing
interactions ($\chi<0$). The top (bottom) row shows the particles color-coded
by their species (polarity orientation), while the thick arrows indicate the
travelling direction of the patterns.  The zoom in (a) illustrates how
disordered configurations appear in grey. The top captions give a graphical
representation of the interaction rules with arrows (bars) denoting motility
activation (inhibition). Parameters: $\bar{\rho}_{A,B}=\rho_{0}=80$,
$\mathrm{Pe}=10$ and $L_{x}=L_{y}=40$. (e) ${\rm Var}(v)$ (see definition in
the text) as function of $\chi$ and ${\rm Pe}$ in the vanishing self-
inhibition regime, triangles and crosses indicate points where traveling bands
and the homogeneous phase are observed. The green shaded region marks the
domain of existence of SOI. (f) Typical linear stability diagram for $\chi<0$
and $\bar{\rho}_{A,B}=\rho_{0}$.

## Description of the model.—

We consider a dynamics where the position $\mathbf{r}_{i}$ and orientation
$\hat{\mathbf{u}}_{i}=(\cos\theta_{i},\sin\theta_{i})$ of particle $i$ from
species $S\in\\{A,B\\}$ evolve at time $t$ according to

$\dot{\mathbf{r}}_{i,S}=v_{S}\left[\tilde{{\rho}}_{A},\tilde{\rho}_{B}\right]\hat{\mathbf{u}}_{i,S},\qquad\dot{\theta}_{i,S}=\sqrt{2D_{r}}\xi_{i}(t).$
(1)

$\xi_{i}$ is a Gaussian white noise with zero mean and unit variance, while
$D_{r}$ denotes the corresponding rotational diffusivity assumed equal for the
two species. Due to quorum-sensing, the self-propulsion velocity $v_{S}$ in
Eq. (1) depends on the coarse-grained density fields
$\tilde{\rho}_{S}(\mathbf{r},t)$ measured over a finite interaction scale $R$
via a short-ranged kernel $\hat{w}(r)$ whose expression is given in Appendix.
The linear stability analysis performed below reveals that, approximating
$\tilde{\rho}_{S}(\mathbf{r},t)\simeq\rho_{S}(\mathbf{r},t)$ with
$\rho_{S}(\mathbf{r},t)\equiv\sum_{i}\delta(\mathbf{r}-\mathbf{r}_{i,S}(t))$,
the emergence of patterns from the dynamics described by (1) is controlled by
the dimensionless couplings
$\eta_{SS^{\prime}}(\rho_{A},\rho_{B})\equiv\rho_{S}\partial\ln(v_{S})/\partial\rho_{S^{\prime}}$.
For $\eta_{SS^{\prime}}<0$, a particle from species $S$ moves slower in higher
$S^{\prime}$ density regions, leading to an effective attraction via motility
inhibition. In turn, ${\eta}_{SS^{\prime}}>0$ makes $S$ particles spend less
time in regions of high $S^{\prime}$ density, such that motility activation
amounts to an effective repulsion. The sign of $\eta_{SS}$ thus determines
whether the effective interaction between particles of same species $S$ is
attractive (repulsive), as a result of the self-inhibition (-activation) of
their motilities. For multiple species, nonreciprocity arises whenever
$\rho_{S^{\prime}}\eta_{SS^{\prime}}\neq\rho_{S}\eta_{S^{\prime}S}$, such that
the dynamics (1) cannot be coarse-grained to an effective equilibrium field
theory [43]. In particular, effective chasing interactions between species $A$
and $B$ are achieved when ${\eta}_{AB}{\eta}_{BA}<0$.

For simplicity, we consider
$v_{S}(\rho_{A},\rho_{B})={v}_{0}L_{SA}(\rho_{A})L_{SB}(\rho_{B})$, where
$L_{SS^{\prime}}(x)>0$ is a logistic function such that $v_{S}$ varies
monotonously with $\rho_{A}$ and $\rho_{B}$. The values of the couplings
$\eta_{SS^{\prime}}$ can then be varied changing either the values of the
densities $\rho_{A,B}$, or the shape of the curves $L_{SS^{\prime}}$. For
convenience, we define a reference value $\rho_{0}$ such that
$v_{S}(\rho_{0},\rho_{0})=v_{0}$ while
$\eta_{SS^{\prime}}^{0}\equiv\eta_{SS^{\prime}}(\rho_{0},\rho_{0})$ is the
value at which $|L^{\prime}_{SS^{\prime}}|$ reaches its maximum. Rescaling
space and time, we set $R=1$ and ${v}_{0}=1$ so that the remaining control
parameters of the microscopic model are the mean particle densities
$\\{\bar{\rho}_{S}\\}$, the nominal couplings $\\{\eta_{SS^{\prime}}^{0}\\}$,
and $D_{r}$. We moreover define the Péclet number $\mathrm{Pe}\equiv
v_{B}(\bar{\rho}_{A},\bar{\rho}_{B})/D_{r}$ as a measure of the self-
propulsion strength, while the parameter
$\chi\equiv\eta_{AB}(\bar{\rho}_{A},\bar{\rho}_{B})\eta_{BA}(\bar{\rho}_{A},\bar{\rho}_{B})$
is used to quantify nonreciprocity. All simulations are performed in periodic
domains of size $L_{x}\times L_{y}$, with total particle number ranging from
$N=10^{4}$ to $10^{7}$.  Additional details about the microscopic model are
given in Appendix.

Traveling patterns induced by chasing interactions.— To start with, we fix the
mean particle densities $\bar{\rho}_{S}=\rho_{0}$, such that the couplings
$\eta_{SS^{\prime}}$ evaluated at $\bar{\rho}_{A,B}$ are given by
$\eta^{0}_{SS^{\prime}}$ and ${\rm Pe}=v_{0}/D_{r}$. For ${\rm Pe}\ll 1$ or
$|\chi|\ll 1$, our observations are in line with predictions from the NRCH
model [10, 11, 27]. Namely, if either of $\eta_{AA}^{0}$ or $\eta_{BB}^{0}$
are sufficiently negative, systems initialized in the homogeneous state are
unstable. For $\chi>0$, this instability leads to the phase separation of one
of the two species, or demixing. For negative $\chi$ and small systems, on the
other hand, nonreciprocity gives rise to the formation of system spanning
traveling bands rich in either of the two species and chasing each other (Fig.
1(a) and Supplemental Movie (SMov) 1). As this TRS broken phase mainly relies
on the presence of two phase-separated domains whose cohesion is maintained by
self-inhibition, TRS can be restored when one or both species see their self-
coupling vanish, as shown in Fig. 1(b).

Strikingly, for $\chi<0$ and $\rm Pe$ and $|\chi|$ large enough we find
traveling patterns arising even at weak or vanishing self-attraction. These
patterns take the form of thin chasing bands which for small $|\chi|$ and
moderate system size self-organize into a superposition of regular
arrangements traveling along different directions (Fig. 1(c) and SMov 1). Such
structure is however rather fragile such that increasing $|\chi|$ (Fig. 1(d)),
or in presence of moderate self-inhibition of motilities, it is destabilized
and replaced by a chaotic chasing bands (CCB) phase (Fig. 2(b)).

Contrary to the regime of strong motility self-inhibition, the cohesion of the
chasing bands here is ensured by the combined effects of $\rm Pe$ and
$\chi<0$. For weak or vanishing self-attraction, traveling bands can only
exist if particles inside them move coherently over large distances, as
evidenced by the local polarization shown in the bottom row of Figs. 1(c,d).
This feature, which is absent in strongly self-attracting mixtures (bottom row
of Figs. 1(a,b)), remarkably arises despite the absence of explicit aligning
interactions between particle velocities (Eq. (1)). The residual particle flux
to the outside of the bands induced by rotational noise is then balanced by
the effect of nonreciprocal inter-species couplings. To understand this, let
us consider the case of Fig. 1(c,d) where species $A$ inhibits the motility of
$B$ particles and is activated by them. A cluster of $A$ particles traveling
in the gas locally slows down $B$ particles, so that they aggregate at its
rear. In the case of vanishing self-attraction, $B$ particles can then follow
the $A$ cluster only when they move in the same direction. Similarly, the
motility activation of $A$ particles provoked by a $B$ cluster leads them to
reside at its front so long as they move in the same direction. This emergent
noise rectification mechanism relying on chasing interactions allows the $AB$
cluster pair to continuously recruit particles from the gas. To confirm this
picture, we define the system- and time-averaged speed variance ${\rm
Var}(v)\equiv\langle\langle|\dot{\bf r}_{i}|^{2}\rangle_{i}-\langle|\dot{\bf
r}_{i}|\rangle_{i}^{2}\rangle_{t}$ which is nonzero only in the presence of
motility induced patterns. Scanning the $(-\chi,{\rm Pe})$ plane at
$\eta_{AA}^{0}=\eta_{BB}^{0}=0$, Fig. 1(e) shows that the CCB phase — whose
onset is characterized by an abrupt growth of ${\rm Var}(v)$ — is found when
both self-propulsion and nonreciprocity are sufficiently strong.

Mesoscopic description & linear stability analysis.— To get a theoretical
understanding of the onset of CCB, we derived the coarse-grained description
of the microscopic model.  The full derivation uses standard coarse-graining
techniques and is detailed in the Appendix. It leads to a pair of equations
for the noise-averaged particle density and polarity. Linearizing these
equations around their solution with homogeneous densities
$(\bar{\rho}_{A},\bar{\rho}_{B})$ and vanishing polarities, we find that its
stability is determined by five parameters: $\eta_{AA}$, $\eta_{BB}$, $\chi$,
$\mathrm{Pe}$ and $\sigma_{v}\equiv{v}_{A}/{v}_{B}$, where $v_{S}$ and
$\eta_{SS^{\prime}}$ are evaluated at $(\bar{\rho}_{A},\bar{\rho}_{B})$.
Taking $\sigma_{v}=1$ (see [51] for a discussion of the general case leading
to similar results), the eigenvalue controlling the stability is
$\lambda(q)=-\tfrac{q}{2}\gamma+\tfrac{q}{2}\sqrt{\gamma^{2}-2\mu+2\sqrt{\alpha}}$,
where $q$ denotes the wavenumber of the perturbation while
$\gamma(q)\equiv(\mathrm{Pe}\,q)^{-1}+\mathrm{Pe}\,q/16$, $\mu\equiv
1+\tfrac{1}{2}({\eta}_{AA}+{\eta}_{BB})$, and
$\alpha\equiv\chi+\tfrac{1}{4}({\eta}_{AA}-{\eta}_{BB})^{2}$.

For $\alpha>0$, $\lambda$ is real as $q\to 0$ and a long-wave stationary
instability (LSI) arises when $\mu<\sqrt{\alpha}$, i.e. $\mu<0$ or
$(1+{\eta}_{AA})(1+{\eta}_{BB})<\chi$. On the other hand, $\alpha=0$
corresponds to an exceptional point [10] such that for $\alpha<0$, the
imaginary part of $\lambda(q)$ is always nonzero, while its real part
$\Re(\lambda)$ becomes positive iff [51]

$2\mu\gamma^{2}(q)+\alpha=2\mu\left(\frac{1}{\mathrm{Pe}\,q}+\frac{\mathrm{Pe}\,q}{16}\right)^{2}+\alpha<0.$
(2)

For $\alpha<0$ and $\mu<0$ the condition (2) is always true down to $q=0$,
giving rise to a long-wave oscillatory instability (LOI). The emergence of LSI
and LOI are thus both mainly controlled by the self-interaction couplings,
such that these instabilities may arise at arbitrary small $|\chi|$ and their
range of existence is insensitive to $\mathrm{Pe}$.

Equation (2) moreover shows that for $\alpha<0$ and $\mu>0$, although
$\Re(\lambda)$ is negative as $q\to 0$, it may turn positive at finite $q$
when $\gamma^{2}(q)<-\alpha/2\mu$. As such a scenario can only happen when
$\lambda$ is complex, it is associated to a short-wave oscillatory instability
(SOI). This is in contrast with the NRCH model [10, 11, 27] or reaction-
diffusion systems 111SOI is forbidden in standard binary reaction-diffusion
systems [32], but can occur under certain conditions, such as systems on
curved surfaces [64] or in the presence of differential flow [65] for which
SOI is absent in binary mixtures. In the limit $\eta_{AA}=\eta_{BB}=0$ of
vanishing self-couplings, the condition for SOI reduces to
$\chi<-\tfrac{1}{2}$ (dotted line in Fig. 1(e)), while the lowest unstable
wavenumber is given by ${\rm
Pe^{2}}q_{c}^{2}=-32(1+4\chi+8\sqrt{(1/2+\chi)\chi})$. Using the interaction
range $R=1$ as a natural lower cutoff for the admissible values of
$q_{c}^{-1}$, we thus get a line in the $(\chi,{\rm Pe})$ plane that defines
the boundary of the instability region shown in Fig. 1(e), and which
qualitatively agrees with the numerical simulations of the microscopic model.
We furthermore summarize the linear stability results for $\chi<0$ in the
diagram of Fig. 1(f) drawn in the $(\eta_{AA},\eta_{BB})$ plane, which
confirms that SOI typically arises when both species present weak or vanishing
self-couplings, while its range of existence increases with ${\rm Pe}$.

CCB-gas phase coexistence.— As we now show, the range of existence of CCB
actually extends beyond that of SOI, since it may correspond to one of the
phases involved in the phase-separation configurations following LOI. We now
fix $\eta^{0}_{AA}=\eta^{0}_{BB}=-2$, $\eta^{0}_{AB}=-\eta^{0}_{BA}=0.5$,
$D_{r}=0.1$, and perform simulations of the microscopic model (1) scanning the
composition plane $(\bar{\rho}_{A},\bar{\rho}_{B})$.

Figure 2: (a) Snapshots from particle-based simulations in the composition
plane in square domains of size $L_{x}=40$. The regions of static phase
separation are indicated by the dashed lines as guides to the eye.  Solid
squares and dots respectively mark the compositions used to analyze the
profiles of Figs. 3(a-c) and coarsening in Fig. 3(e). Inset: the corresponding
linear stability diagram. (b–d) Simulations in square domains of size
$L_{x}=160$, each with the same composition as its top right inset where
$L_{x}=40$.

The resulting phase diagram shown in Fig. 2(a) exhibits four distinctive lobes
encircled by black dashed lines where static phase separation takes place
[10], in qualitative agreement with the location of LSI in the linear
stability diagram ( blue regions in the inset of Fig. 2(a)).  Along the
diagonal $\bar{\rho}_{A}=\bar{\rho}_{B}$ marked by a white dashed line in Fig.
2(a), the homogeneous phase becomes unstable below densities $\approx
2\rho_{0}$ and is superseded by the CCB phase (Fig. 2(b) and SMov 2)
previously described. Further decreasing the densities, we observe system
spanning traveling bands for $\bar{\rho}_{B}\gtrsim\bar{\rho}_{A}$ and
rotating spirals for $\bar{\rho}_{B}\lesssim\bar{\rho}_{A}$.  Examples of
these patterns are shown in the insets of Figs. 2(c,d) and SMov 2. Increasing
the system size, the CCB phase remains qualitatively unchanged,  as can be
seen comparing the main panel and inset of Fig. 2(b). On the other hand, both
traveling bands and spirals found at lower densities are replaced by a
macroscopic CCB domain coexisting with a dilute homogeneous gas (Fig. 2(c,d)).
Accordingly, the locations of phase-separated CCB domains qualitatively
correspond to regions of the composition plane where LOI is found (maked in
magenta in the inset of Fig. 2(a)). The system-wide CCB phase, on the other
hand, borders the uniform high-density phase, in line with the location of SOI
highlighted in green in the inset of Fig. 2(a).

Figure 3: (a–c) $y$\- and time-averaged density (a), $x$-component of the
polarity (b) and speed (c) profiles for both species at various compositions
along the tie line with $L_{x}=6L_{y}=960$. (d) CCB-gas interface for the
configuration corresponding to the red curves in (a–c). (e) Typical domain
size $\zeta$ evaluated from the first moment of the structure factor [51] as
function of time for randomly initialized systems in the static phase
separation (SPS), CCB-gas phase coexistence and pure CCB regimes. For each
curve in (a–c,e), the corresponding composition is marked by a symbol of the
same color in Fig. 2(a).

To further characterize the CCB-gas coexistence phase, we performed
simulations within this regime in a large rectangular domain. This way, the
CCB domain connects with itself along the shortest dimension $y$, leading to a
well-defined interface (Fig. 3(d)).  The $y$ and time-averaged density
profiles shown in Fig. 3(a) indeed highlight two bulk phases of nearly
constant densities. The values taken by $\langle\rho_{A,B}\rangle_{y,t}$ far
away from the interfaces thus define a pair of points in the composition plane
which can be joined by a tie line. Consistently with the usual picture of
phase separation, shifting $(\bar{\rho}_{A},\bar{\rho}_{B})$ along this tie
line changes the relative proportions of the gas and CCB phases while the
corresponding bulk densities are left unchanged.

TRS is however obviously broken at mesoscopic scales in the CCB-gas
coexistence phase, since chasing bands are constantly created in the dense CCB
domain and expelled in the gas, where they quickly dissolve (Fig. 3(d) and
SMov 3). The outward mass flux from the chasing bands must then be balanced by
the diffusive flux from the resulting excess particles in the gas, thereby
maintaining the cohesion of the CCB domain. The nontrivial structure of the
CCB-gas interface is also highlighted by the presence of local polarization
pointing towards the dilute regions (Fig. 3(b)), at odds with the usual MIPS
phenomenology [53, 54]. We moreover note from Figs. 3(a,c) that $A$ particles
accumulate at interfaces where their self-propulsion speed is lowest, but move
on average faster in the dense CCB domain than in the dilute gas.

Remarkably, the nonequilibrium nature of the CCB-gas phase coexistence also
emerges over macroscopic scales, as revealed by the anomalous coarsening
behavior shown in Fig. 3(f). Although in the regimes of pure CCB and static
phase separation the coarsening is arrested or follows the Lifshitz–Slyozov
$t^{1/3}$ law [55] (green and blue symbols in Fig. 3(f)), the late-time
coarsening of phase-separated CCB domains is characterized by an exponent
$\approx 0.13$ (magenta symbols in Fig. 3(f)), significantly smaller than the
$\tfrac{1}{3}$ value expected in passive systems [56] and for the coarsening
of dense MIPS domains [47, 49, 48, 50, 57]. We rationalize this result by
noting that larger CCB domains generate more bands —and thus expel more mass
into the surrounding gas— than smaller domains, which may naturally slow down
coarsening. 222We note that for an underdamped dynamics the mapping of MIPS to
an equilibrium coarse-grained model at large scales may not hold [66].

Using a minimal model including self-propulsion and nonreciprocity, we have
shown how the combination of these two sources of activity leads to a chaotic
chasing band phase. Since it is involved in phase-separated configurations at
large nonreciprocity, this phase is moreover found in a large portion of the
phase diagram. Additionally, although the corresponding phase-separated
domains are globally static, they still defy an equilibrium mapping as
evidenced by the observed abnormal coarsening behavior.  As they allow for the
design of programmable quorum-sensing motility responses, both genetically
engineered E. coli [58] and light-controlled microswimmers [59, 38, 60, 61]
offer promising experimental platforms to observe this new type of self-
organized behavior.

###### Acknowledgements.

We thank Xiaqing Shi and Hugues Chaté for their critical reading of our
manuscript. This work has received support from the Max Planck School Matter
to Life and the MaxSynBio Consortium, which are jointly funded by the Federal
Ministry of Education and Research (BMBF) of Germany, and the Max Planck
Society.
Appendix on simulations of the microscopic model.— Here, we give additional
details on the microscopic model (1) and discuss the underlying assumptions
behind its formulation. The Langevin equations (1) were integrated using the
Euler-Maruyama scheme with time resolution $\Delta t=0.1$. The local density
fields $\tilde{\rho}_{S}({\bf r},t)$ entering the expression of the quorum-
sensing interaction in Eq. (1) were computed as
$\tilde{\rho}_{S}(\mathbf{r},t)=\sum_{i}\hat{w}(|\mathbf{r}-\mathbf{r}_{i,S}(t)|),$
with the linear weight function

$\hat{w}(r)=\begin{cases}\frac{3}{\pi
R^{2}}(1-\frac{r}{R})\quad&\mathrm{if}\quad r<R,\\\ 0\quad&\mathrm{if}\quad
r\geq R,\end{cases}$

that satisfies the normalization condition $\int_{\mathbb{R}^{2}}{\rm
d}\mathbf{r}\,\hat{w}(|\mathbf{r}|)=1$. As mentioned in the main text, the
dependency of the particles self-propulsion speed with their density is
modeled as
$v_{S}(\tilde{\rho}_{A},\tilde{\rho}_{B})={v}_{0}L_{SA}(\tilde{\rho}_{A})L_{SB}(\tilde{\rho}_{B})$
with

$L_{SS^{\prime}}(x)\equiv
1+\kappa\tanh\left(\frac{\eta_{SS^{\prime}}^{0}}{\kappa}\frac{x-\rho_{0}}{\rho_{0}}\right).$

Therefore, the reference density value $\rho_{0}$ satisfies
$L_{SS^{\prime}}(\rho_{0})=1$ and
$\partial_{x}L_{SS^{\prime}}(\rho_{0})=\eta_{SS^{\prime}}^{0}$. Here, $\kappa$
controls the lower and upper limits of $L_{SS^{\prime}}(x)$. To prevent
jamming of particles at high densities, we set $\kappa=0.7$ so that the
minimum speed of species $S$ is given by $\approx 0.09\,v_{0}$.

To keep the microscopic model minimal, we have assumed a direct dependency of
the particle motility in the local density fields. We expect this
approximation to be well-verified for synthetic colloids controlled by optical
feedback loops, which allow for a fast motility response to dynamical
variations of the local density [38, 23]. In the case of microorganisms
communicating via chemical signals, the limit of instantaneous quorum-sensing
response requires that the timescales associated with diffusion of chemicals
and the internal gene regulatory network are fast as compared to the spatial
dynamics of particles. Whether this assumption is verified depends on the
system of interest, but the limit of fast quorum-sensing response should in
general qualitatively capture the relevant dynamics [37].

Another simplification used in the microscopic model (1) consists in
neglecting the short-range repulsive interactions between particles. This
assumption is well-justified so long as the range of quorum-sensing
interactions $R$ is much larger than the typical particles size $\sigma$. In
synthetic systems, $R$ can be tuned arbitrarily, while $R=10\sigma$ was used
in Ref. [23]. In the case of quorum-sensing regulated by chemicals, the scale
$R$ is set by the typical distance a signaling molecule can diffuse before
being degraded. Estimates for the acyl-homoserine lactone molecules that
mediate the quorum-sensing interactions of bacteria in Refs. [36, 37] can
reach a few millimeters [62], thus several orders of magnitude larger then the
typical size of a bacterium.

Appendix on the coarse-grained equations.— To derive the field theory
describing the binary mixture of self-propelled particles with quorum-sensing
interactions, we follow the coarse-graining framework developed for a single-
species [63]. As a starting point, we consider the many-body probability
distribution $P(\mathbf{X},t)$ with
$\mathbf{X}\equiv\\{\mathbf{r}_{1},\mathbf{r}_{2},\dots,\mathbf{r}_{N},\theta_{1},\theta_{2},\dots,\
\theta_{N}\\}$, and where $N$ denotes the total number of particles of the two
species. From the Langevin formulation (1), we get the statistically
equivalent Fokker-Planck equation

$\partial_{t}P=-\sum_{i=1}^{N}\left[\nabla_{\mathbf{r}_{i}}\cdot\left(v_{s_{i}}\hat{\mathbf{u}}(\theta_{i})P\right)-D_{r}\partial_{\theta_{i}}^{2}P\right],$
(A1)

where $\hat{\mathbf{u}}(\theta)=(\cos\theta,\sin\theta)$ and $s_{i}=A$ or $B$
denotes the species of particle $i$. Without loss of generality, we set
$s_{i}=A$ for $i=1,2,\dots,N_{A}$ and $s_{i}=B$ for
$i=N_{A}+1,N_{A}+2,\dots,N$ with $N_{A}$ the total number of $A$ particles.
The one-particle probability densities $f_{S}(\mathbf{r},\theta,t)$ are
obtained from $P$ by integrating over all degrees of freedom except that of
one particle of species $S$. Hence, we have
$f_{A}(\mathbf{r}_{1},\theta_{1},t)\equiv
N_{A}\left(\prod_{i=2}^{N}\int\mathrm{d}\mathbf{r}_{i}\int_{0}^{2\pi}\mathrm{d}\theta_{i}\right)P(\mathbf{X},t)$
while $f_{B}$ is defined in a similar way. Integrating Eq. (A1) over the
relevant coordinates, we thus determine the dynamics of the $A$ and $B$
single-particle distributions. For $f_{A}$ (the generalization to the species
$B$ being straightforward), we obtain

$\displaystyle\partial_{t}f_{A}(\mathbf{r}_{1},\theta_{1},t)$
$\displaystyle=-\nabla_{\bf{r}_{1}}\cdot\left[{\cal
G}_{A}(\mathbf{r}_{1},\theta_{1},t)\hat{\bf{u}}(\theta_{1})f_{A}(\mathbf{r}_{1},\theta_{1},t)\right]$
$\displaystyle+D_{r}\partial_{\theta_{1}}^{2}f_{A}(\mathbf{r}_{1},\theta_{1},t),$
(A2)

where ${\cal
G}_{A}(\mathbf{r}_{1},\theta_{1},t)\equiv\left(\prod_{k=2}^{N}\int\mathrm{d}\mathbf{r}_{k}\right)v_{A}({\bf
X}_{r})g_{A}({\bf X}_{r},\theta_{1},t)$, and we have used the shorthand
notation $\mathbf{X}_{r}\equiv\\{\mathbf{r}_{1},\dots,\mathbf{r}_{N}\\}$. The
function $g_{A}$ is defined by
$N_{A}\prod_{k=2}^{N}\int_{0}^{2\pi}\mathrm{d}\theta_{k}P({\bf X},t)\equiv
f_{A}(\mathbf{r}_{1},\theta_{1},t)g_{A}({\bf X}_{r},\theta_{1},t)$. In
general, the quorum-sensing interaction term ${\cal G}_{A}$ in (Description of
the model.—) depends on the full many-body distribution. However, we now use
the fact that the particle speed is formally a function of
$\tilde{\rho}_{A,B}({\bf X}_{r})$, which in the mean field limit can be
expressed in terms of the single particle densities:

$\displaystyle\tilde{\rho}_{S}({\bf X}_{r})$
$\displaystyle=\sum_{j}\hat{w}(|\mathbf{r}_{1}-\mathbf{r}_{j}(t)|)\delta_{s_{j},S}$
$\displaystyle\simeq\int\mathrm{d}\mathbf{r}^{\prime}\hat{w}(|\mathbf{r}_{1}-\mathbf{r}^{\prime}|)\phi_{S}(\mathbf{r}^{\prime},t)\equiv\tilde{\phi}_{S}({\bf
r}_{1},t),$

where the $S$ particle density is formally defined as $\phi_{S}({\bf
r},t)\equiv\int_{0}^{2\pi}{\rm d}\theta f_{S}(\mathbf{r},\theta,t)$. We expect
this approximation to be reasonably valid in sufficiently dense systems. In
the mean field limit considered here, the function $v_{A}$ in the expression
of ${\cal G}_{A}$ thus only depends on the position variable ${\bf r}_{1}$.
Using that by definition
$\left(\prod_{k=2}^{N}\int\mathrm{d}\mathbf{r}_{k}\right)g_{A}({\bf
X}_{r},\theta_{1},t)=1$, we therefore obtain after dropping indices on the
$\bf r$ and $\theta$ variables:

$\partial_{t}f_{S}=-\nabla\cdot\left[v_{S}(\tilde{\phi}_{A},\tilde{\phi}_{B})\hat{\bf{u}}(\theta)f_{S}\right]+D_{r}\partial_{\theta}^{2}f_{S}.$
(A3)

To simplify Eq. (A3) further, we use the local approximation
$\tilde{\phi}_{S}(\mathbf{r})\approx\phi_{S}(\mathbf{r})$, such that
$v_{S}(\tilde{\phi}_{A},\tilde{\phi}_{B})\approx v_{S}(\phi_{A},\phi_{B})$.

Due to angular diffusion, we expect the dynamics of the system to be well
captured over long timescales by that of the low order orientational moments
of the distributions $f_{A}$ and $f_{B}$. Therefore, we expand the
distribution $f_{S}$ in angular Fourier modes:
$f_{S}(\mathbf{r},\theta,t)=(2\pi)^{-1}\sum_{k=-\infty}^{\infty}f_{k,S}(\mathbf{r},t)\exp(-ik\theta)$.
It is straightforward to check that the first three modes of $f_{S}$
correspond to the complex representation of the density $\phi_{S}({\bf r},t)$,
polarity $\mathbf{p}_{S}(\mathbf{r},t)$ and nematic order
$\mathbf{Q}_{S}(\mathbf{r},t)$ fields. For the derivation below, it is
convenient to work with complex notations for which $\hat{\bf
u}(\theta)\leftrightarrow e^{i\theta}$. We therefore define the complex
gradient $\triangledown\equiv\partial_{x}+i\partial_{y}$, and obtain from Eq.
(A3)

$\displaystyle\partial_{t}f_{k,S}$
$\displaystyle=-\frac{1}{2}\triangledown^{*}\left[v_{S}(\phi_{A},\phi_{B})f_{k+1,S}\right]$
$\displaystyle-\frac{1}{2}\triangledown\left[v_{S}(\phi_{A},\phi_{B})f_{k-1,S}\right]-D_{r}k^{2}f_{k,S},$
(A4)

where star denotes complex conjugate. The equation for the $k^{\rm th}$
angular mode of $f_{S}$ contains a linear damping term $-D_{r}k^{2}f_{k,S}$.
Considering the long time and large scale limits, we thus enslave the dynamics
of the high order modes to that of the slow ones. As the densities are the
only conserved fields, it is customary in this context to retain only them as
hydrodynamic fields, while enslaving higher order modes [40]. Motivated by the
presence of local polarization in the CCB phase, here we instead retain both
$\phi_{S}$ and the polarity fields $\mathbf{p}_{S}(\mathbf{r},t)$. Namely,
neglecting $f_{k,S}$ for $k\geq 3$, we get from (Description of the model.—)
closed equations for $\phi_{S}$, $f_{1,S}$ and $f_{2,S}$. Setting
$\partial_{t}f_{2,S}=0$, we solve the equation for $f_{2,S}$ which leads to
$f_{2,S}=-(8D_{r})^{-1}\triangledown(v_{S}f_{1,S})$. Replacing this expression
in the equation for $f_{1,S}$, we get after going back to vector notations

$\displaystyle\partial_{t}\phi_{S}=$
$\displaystyle-\nabla\cdot(v_{S}\mathbf{p}_{S}),$ (A5a)
$\displaystyle\partial_{t}\mathbf{p}_{S}=$
$\displaystyle-\frac{1}{2}\nabla(v_{S}\phi_{S})-D_{r}\mathbf{p}_{S}+\frac{v_{S}}{16D_{r}}\Delta\left(v_{S}\mathbf{p}_{S}\right)$
$\displaystyle+(8D_{r})^{-1}\left[\nabla(v_{S}\mathbf{p}_{S})\right]_{\rm
ST}\cdot\nabla v_{S},$ (A5b)

where $\left[{\bf A}\right]_{\rm ST}\equiv\tfrac{1}{2}\left[{\bf A}+{\bf
A}^{T}-{\bf I}{\rm Tr}({\bf A})\right]$ is the symmetric and traceless part of
the tensor $\bf A$. Performing the linear stability analysis of Eqs. (A5)
around their homogeneous disordered solution $\phi_{S}=\bar{\rho}_{S}$ and
${\bf p}_{S}={\bf 0}$, we obtain the results presented in the main text.

## References

  * Romanczuk  _et al._ [2012] P. Romanczuk, M. Bär, W. Ebeling, B. Lindner, and L. Schimansky-Geier, Active brownian particles, Eur. Phys. J. Spec. Top. 202, 1 (2012).
  * Prost  _et al._ [2015] J. Prost, F. Jülicher, and J.-F. Joanny, Active gel physics, Nat. Phys. 11, 111 (2015).
  * Doostmohammadi  _et al._ [2018] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Active nematics, Nat. Commun. 9, 3246 (2018).
  * Uchida and Golestanian [2010] N. Uchida and R. Golestanian, Synchronization and collective dynamics in a carpet of microfluidic rotors, Phys. Rev. Lett. 104, 178103 (2010).
  * Fürthauer  _et al._ [2012] S. Fürthauer, M. Strempel, S. W. Grill, and F. Jülicher, Active chiral fluids, Eur. Phys. J. E 35, 89 (2012).
  * Dell’Arciprete  _et al._ [2018] D. Dell’Arciprete, M. L. Blow, A. T. Brown, F. D. C. Farrell, J. S. Lintuvuori, A. F. McVey, D. Marenduzzo, and W. C. K. Poon, A growing bacterial colony in two dimensions as an active nematic, Nat. Commun. 9, 4190 (2018).
  * Golestanian [2022] R. Golestanian, Phoretic Active Matter, in  _Active Matter and Nonequilibrium Statistical Physics: Lecture Notes of the Les Houches Summer School: Volume 112, September 2018_ (Oxford University Press, 2022).
  * Gompper et al. [2020] G. Gompper et al., The 2020 motile active matter roadmap, J. Phys. Condens. Matter 32, 193001 (2020).
  * Ivlev  _et al._ [2015] A. V. Ivlev, J. Bartnick, M. Heinen, C.-R. Du, V. Nosenko, and H. Löwen, Statistical mechanics where newton’s third law is broken, Phys. Rev. X 5, 011035 (2015).
  * Saha  _et al._ [2020] S. Saha, J. Agudo-Canalejo, and R. Golestanian, Scalar active mixtures: The nonreciprocal Cahn-Hilliard model, Phys. Rev. X 10, 041009 (2020).
  * You  _et al._ [2020] Z. You, A. Baskaran, and M. C. Marchetti, Nonreciprocity as a generic route to traveling states, Proc. Natl. Acad. Sci. USA. 117, 19767 (2020).
  * Fruchart  _et al._ [2021] M. Fruchart, R. Hanai, P. B. Littlewood, and V. Vitelli, Non-reciprocal phase transitions, Nature 592, 363 (2021).
  * Gelbwaser-Klimovsky  _et al._ [2021] D. Gelbwaser-Klimovsky, N. Graham, M. Kardar, and M. Krüger, Near field propulsion forces from nonreciprocal media, Phys. Rev. Lett. 126, 170401 (2021).
  * Dzubiella  _et al._ [2003] J. Dzubiella, H. Löwen, and C. N. Likos, Depletion forces in nonequilibrium, Phys. Rev. Lett. 91, 248301 (2003).
  * Soto and Golestanian [2014] R. Soto and R. Golestanian, Self-assembly of catalytically active colloidal molecules: Tailoring activity through surface chemistry, Phys. Rev. Lett. 112, 068301 (2014).
  * Wysocki  _et al._ [2016] A. Wysocki, R. G. Winkler, and G. Gompper, Propagating interfaces in mixtures of active and passive Brownian particles, New J. Phys. 18, 123030 (2016).
  * Wittkowski  _et al._ [2017] R. Wittkowski, J. Stenhammar, and M. E. Cates, Nonequilibrium dynamics of mixtures of active and passive colloidal particles, New J. Phys. 19, 105003 (2017).
  * Agudo-Canalejo and Golestanian [2019] J. Agudo-Canalejo and R. Golestanian, Active phase separation in mixtures of chemically interacting particles, Phys. Rev. Lett. 123, 018101 (2019).
  * Nasouri and Golestanian [2020] B. Nasouri and R. Golestanian, Exact phoretic interaction of two chemically active particles, Phys. Rev. Lett. 124, 168003 (2020).
  * Banerjee  _et al._ [2022] J. P. Banerjee, R. Mandal, D. S. Banerjee, S. Thutupalli, and M. Rao, Unjamming and emergent nonreciprocity in active ploughing through a compressible viscoelastic fluid, Nat. Commun. 13, 1 (2022).
  * Gupta  _et al._ [2022] R. K. Gupta, R. Kant, H. Soni, A. K. Sood, and S. Ramaswamy, Active nonreciprocal attraction between motile particles in an elastic medium, Phys. Rev. E 105, 064602 (2022).
  * Chen  _et al._ [2017] Q.-s. Chen, A. Patelli, H. Chaté, Y.-q. Ma, and X.-q. Shi, Fore-aft asymmetric flocking, Phys. Rev. E 96, 020601(R) (2017).
  * Bäuerle  _et al._ [2018] T. Bäuerle, A. Fischer, T. Speck, and C. Bechinger, Self-organization of active particles by quorum sensing rules, Nat. Commun. 9, 1 (2018).
  * Gómez-Nava  _et al._ [2022] L. Gómez-Nava, R. Bon, and F. Peruani, Intermittent collective motion in sheep results from alternating the role of leader and follower, Nat. Phys. 18, 1494 (2022).
  * El-Ganainy  _et al._ [2018] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Non-Hermitian physics and PT symmetry, Nat. Phys. 14, 11 (2018).
  * Shankar  _et al._ [2022] S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti, and V. Vitelli, Topological active matter, Nat. Rev. Phys. 4, 380 (2022).
  * Frohoff-Hülsmann  _et al._ [2021] T. Frohoff-Hülsmann, J. Wrembel, and U. Thiele, Suppression of coarsening and emergence of oscillatory behavior in a Cahn-Hilliard model with nonvariational coupling, Phys. Rev. E 103, 042602 (2021).
  * Ouazan-Reboul  _et al._ [2021] V. Ouazan-Reboul, J. Agudo-Canalejo, and R. Golestanian, Non-equilibrium phase separation in mixtures of catalytically active particles: size dispersity and screening effects, Eur. Phys. J. E 44, 113 (2021).
  * Loos  _et al._ [2023] S. A. M. Loos, S. H. Klapp, and T. Martynec, Long-range order and directional defect propagation in the nonreciprocal xy model with vision cone interactions, Phys. Rev. Lett. 130, 198301 (2023).
  * Saha and Golestanian [2022] S. Saha and R. Golestanian, Effervescent waves in a binary mixture with non-reciprocal couplings, arXiv preprint arXiv:2208.14985 (2022).
  * Cahn and Hilliard [2004] J. W. Cahn and J. E. Hilliard, Free Energy of a Nonuniform System. I. Interfacial Free Energy, J. Chem. Phys. 28, 258 (2004).
  * Diego  _et al._ [2018] X. Diego, L. Marcon, P. Müller, and J. Sharpe, Key features of turing systems are determined purely by network topology, Phys. Rev. X 8, 021071 (2018).
  * Haas and Goldstein [2021] P. A. Haas and R. E. Goldstein, Turing’s diffusive threshold in random reaction-diffusion systems, Phys. Rev. Lett. 126, 238101 (2021).
  * Frohoff-Hülsmann  _et al._ [2023] T. Frohoff-Hülsmann, U. Thiele, and L. M. Pismen, Non-reciprocity induces resonances in a two-field Cahn–Hilliard model, Philos. Trans. Royal Soc. A 381, 20220087 (2023).
  * Daniels  _et al._ [2004] R. Daniels, J. Vanderleyden, and J. Michiels, Quorum sensing and swarming migration in bacteria, FEMS Microbiol. Rev. 28, 261 (2004).
  * Liu  _et al._ [2011] C. Liu, X. Fu, L. Liu, X. Ren, C. K. Chau, S. Li, L. Xiang, H. Zeng, G. Chen, L.-H. Tang, P. Lenz, X. Cui, W. Huang, T. Hwa, and J.-D. Huang, Sequential establishment of stripe patterns in an expanding cell population, Science 334, 238 (2011).
  * Curatolo  _et al._ [2020] A. Curatolo, N. Zhou, Y. Zhao, C. Liu, A. Daerr, J. Tailleur, and J. Huang, Cooperative pattern formation in multi-component bacterial systems through reciprocal motility regulation, Nat. Phys. 16, 1152 (2020).
  * Gomez-Solano  _et al._ [2017] J. R. Gomez-Solano, S. Samin, C. Lozano, P. Ruedas-Batuecas, R. van Roij, and C. Bechinger, Tuning the motility and directionality of self-propelled colloids, Scientific Reports 7, 14891 (2017).
  * Lavergne  _et al._ [2019] F. A. Lavergne, H. Wendehenne, T. Bäuerle, and C. Bechinger, Group formation and cohesion of active particles with visual perception–dependent motility, Science 364, 70 (2019).
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, Motility-induced phase separation, Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
  * Solon  _et al._ [2018] A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri, and J. Tailleur, Generalized thermodynamics of motility-induced phase separation: phase equilibria, laplace pressure, and change of ensembles, New J. Phys. 20, 075001 (2018).
  * Fischer  _et al._ [2020] A. Fischer, F. Schmid, and T. Speck, Quorum-sensing active particles with discontinuous motility, Phys. Rev. E 101, 012601 (2020).
  * Dinelli  _et al._ [2022] A. Dinelli, J. O’Byrne, A. Curatolo, Y. Zhao, P. Sollich, and J. Tailleur, Non-reciprocity across scales in active mixtures, arXiv preprint arXiv:2203.07757 (2022).
  * Chiu and Omar [2023] Y.-J. Chiu and A. K. Omar, Phase coexistence implications of violating newton’s third law, J. Chem. Phys. 158 (2023).
  * Mandal  _et al._ [2022] R. Mandal, S. S. Jaramillo, and P. Sollich, Robustness of travelling states in generic non-reciprocal mixtures, arXiv preprint arXiv:2212.05618 (2022).
  * Kreienkamp and Klapp [2022] K. L. Kreienkamp and S. H. L. Klapp, Clustering and flocking of repulsive chiral active particles with non-reciprocal couplings, New J. Phys. 24, 123009 (2022).
  * Stenhammar  _et al._ [2014] J. Stenhammar, D. Marenduzzo, R. J. Allen, and M. E. Cates, Phase behaviour of active Brownian particles: the role of dimensionality, Soft Matter 10, 1489 (2014).
  * Caporusso  _et al._ [2020] C. B. Caporusso, P. Digregorio, D. Levis, L. F. Cugliandolo, and G. Gonnella, Motility-induced microphase and macrophase separation in a two-dimensional active Brownian particle system, Phys. Rev. Lett. 125, 178004 (2020).
  * Shi  _et al._ [2020] X.-q. Shi, G. Fausti, H. Chaté, C. Nardini, and A. Solon, Self-organized critical coexistence phase in repulsive active particles, Phys. Rev. Lett. 125, 168001 (2020).
  * Caporusso  _et al._ [2023] C. B. Caporusso, L. F. Cugliandolo, P. Digregorio, G. Gonnella, D. Levis, and A. Suma, Dynamics of Motility-Induced Clusters: Coarsening beyond Ostwald Ripening, Phys. Rev. Lett. 131, 068201 (2023).
  * [51] See Supplemental Material at [URL will be inserted by publisher] for further details.
  * Note [1] SOI is forbidden in standard binary reaction-diffusion systems [32], but can occur under certain conditions, such as systems on curved surfaces [64] or in the presence of differential flow [65].
  * Omar  _et al._ [2020] A. K. Omar, Z.-G. Wang, and J. F. Brady, Microscopic origins of the swim pressure and the anomalous surface tension of active matter, Phys. Rev. E 101, 012604 (2020).
  * Mahault  _et al._ [2023] B. Mahault, P. Godara, and R. Golestanian, Emergent organization and polarization due to active fluctuations, Phys. Rev. Res. 5, L022012 (2023).
  * Lifshitz and Slyozov [1961] I. Lifshitz and V. Slyozov, The kinetics of precipitation from supersaturated solid solutions, J. Phys. Chem. Solids 19, 35 (1961).
  * Bray [2002] A. J. Bray, Theory of phase-ordering kinetics, Adv. Phys. 51, 481 (2002).
  * Note [2] We note that for an underdamped dynamics the mapping of MIPS to an equilibrium coarse-grained model at large scales may not hold [66].
  * Weiss  _et al._ [2008] L. E. Weiss, J. P. Badalamenti, L. J. Weaver, A. R. Tascone, P. S. Weiss, T. L. Richard, and P. C. Cirino, Engineering motility as a phenotypic response to luxi/r-dependent quorum sensing in Escherichia coli, Biotechnol. Bioeng. 100, 1251 (2008).
  * Bregulla  _et al._ [2014] A. P. Bregulla, H. Yang, and F. Cichos, Stochastic localization of microswimmers by photon nudging, ACS Nano 8, 6542 (2014).
  * Massana-Cid  _et al._ [2022] H. Massana-Cid, C. Maggi, G. Frangipane, and R. Di Leonardo, Rectification and confinement of photokinetic bacteria in an optical feedback loop, Nat. Commun. 13, 2740 (2022).
  * Yang  _et al._ [2021] S. Yang, M. Huang, Y. Zhao, and H. P. Zhang, Controlling cell motion and microscale flow with polarized light fields, Phys. Rev. Lett. 126, 058001 (2021).
  * Marenda  _et al._ [2016] M. Marenda, M. Zanardo, A. Trovato, F. Seno, and A. Squartini, Modeling quorum sensing trade-offs between bacterial cell density and system extension from open boundaries, Sci. Rep. 6, 39142 (2016).
  * te Vrugt  _et al._ [2023] M. te Vrugt, J. Bickmann, and R. Wittkowski, How to derive a predictive field theory for active brownian particles: a step-by-step tutorial, J. Phys. Condens. Matter 35, 313001 (2023).
  * Nishide and Ishihara [2022] R. Nishide and S. Ishihara, Pattern propagation driven by surface curvature, Phys. Rev. Lett. 128, 224101 (2022).
  * Rovinsky and Menzinger [1992] A. B. Rovinsky and M. Menzinger, Chemical instability induced by a differential flow, Phys. Rev. Lett. 69, 1193 (1992).
  * Mandal  _et al._ [2019] S. Mandal, B. Liebchen, and H. Löwen, Motility-induced temperature difference in coexisting phases, Phys. Rev. Lett. 123, 228001 (2019).


# Variance-reduced sampling importance resampling

Yao Xiao1, Kang Fu2, and Kun Li3
1 School of Mathematics and Statistics, Beijing Institute of Technology,
Beijing, China.
2 School of Mathematics and Statistics, Central China Normal University,
Wuhan, China.
3 School of Statistics and Mathematics, Zhongnan University of Economics and
Law, Wuhan, China Corresponding author<EMAIL_ADDRESS>
###### Abstract

The sampling importance resampling method is widely utilized in various
fields, such as numerical integration and statistical simulation. In this
paper, two modified methods are presented by incorporating two variance
reduction techniques commonly used in Monte Carlo simulation, namely
antithetic sampling and Latin hypercube sampling, into the process of sampling
importance resampling method respectively. Theoretical evidence is provided to
demonstrate that the proposed methods significantly reduce estimation errors
compared to the original approach. Furthermore, the effectiveness and
advantages of the proposed methods are validated through both numerical
studies and real data analysis.

Keywords: Statistical simulation; sampling importance resampling; antithetic
sampling; Latin hypercube sampling; variance reduction.

##  1 Introduction

The Monte Carlo (MC) method is a widespread simulation technique to validate
the efficiency and effect of the statistical theory and method (Ross, 2022).
The basic idea of the MC method is to simulate random sampling using computer
software and then analyze statistical problems by collecting sample values
generated through simulation. It is well-known that the estimators obtained
using the MC method often have a large variance. As a result, a series of
variance reduction techniques for the MC method have been proposed, such as
antithetic variable method, control variable method, Latin hypercube sampling
method, importance sampling method, and so on (Givens and Hoeting, 2013;
Robert and Casella, 2005). Among them, importance sampling has gained much
popularity (Tokdar and Kass, 2010). The basic idea behind importance sampling
is to introduce a proposal distribution and sample from it, thereby giving
greater “importance” to the parts of the original distribution that have a
significant contribution. Furthermore, some modified methods have also been
proposed, such as adaptive importance sampling, sequential importance
sampling, and sampling importance resampling.

In Bayesian statistics, one often needs to sample from a target distribution
whose regularization constant is unknown or is too complex to compute. The
sampling importance resampling (SIR) method proposed by Rubin (1987) is one of
the efficient methods to deal with this sampling problem. The process of this
method is to generate random samples from the introduced proposal distribution
to construct a sample pool, and then randomly resample from the sample pool
with a certain importance weight, hence the final obtained samples are
asymptotically i.i.d from the target distribution. In this method, the
proposal distribution is usually given and easy to implement. The importance
weight measures the difference between the proposal distribution and the
target distribution, which can make samples from the target distribution
appeal with higher frequencies in the resampling process.

In order to improve the performance and applicability of the SIR method, some
modified methods have been proposed. For example, Givens and Raftery (1996)
proposed the adaptive SIR method in which the sampling process is repeated
many cycles, and for each cycle, the proposal distribution is updated through
the kernel density estimator. Skare et al. (2003) proposed a simple adjustment
of the sample probabilities for the SIR algorithm which results in a faster
convergence. Perez et al. (2005) proposed two quasi-Monte Carlo versions of
the SIR method, that is, the MC points (simple random samples) involved in the
SIR method are replaced with the quasi-Monte Carlo points (such as low
discrepancy sequence and number-theoretic net (Fang and Wang, 1994)), which
can lead to more representative target distribution samples than SIR method.
And the convergence of the quasi-random SIR method has been proved by
Vandewoestyne and Cools (2010). Recently, Ning and Tao (2020) embedded the
randomized quasi-Monte Carlo technique into the SIR method, which introduces
randomness into the quasi-random SIR method. Following the idea of SIR method,
Yi et al. (2023) proposed the global likelihood sampler that is independent of
the proposal distributions. In this work, in order to reduce variance and to
improve the sampling effect of the SIR method, two variance reduction
techniques are embedded into the resampling process of SIR method. The
modified sampling methods can get more accurate estimators in statistical
inference problems compared with the original method.

The rest of this paper is organized as follows. In Section 2, the original SIR
method is introduced briefly. Then the sampling process and properties of the
proposed methods are demonstrated in Section 3. The numerical simulations and
a real data case are used to validate the effect of the proposed methods in
Section 4 and Section 5. Section 6 concludes this paper.

##  2 Sampling importance resampling

Suppose that the density function of a random variable $\bm{X}$ is
$f(\bm{x})$, then the expectation of the random function $h(\bm{X})$ can be
denoted as $\mu:=E[h(\bm{X})]=\int h(\bm{x})f(\bm{x})d\bm{x}$. If
$\bm{x}_{1},\bm{x}_{2},\cdots,\bm{x}_{n}$ are i.i.d. samples drawn from
$f(\bm{x})$, according to the law of large number, the expectation can be
approximated by the sample mean, i.e.,

$\hat{\mu}_{MC}=\frac{1}{n}\sum_{i=1}^{n}h(\bm{x}_{i})\rightarrow\int
h(\bm{x})f(\bm{x})d\bm{x}=\mu,\quad n\rightarrow\infty.$ (2.1)

Thus $\hat{\mu}_{MC}$ is usually deemed as the MC estimator of $\mu$ when
$\bm{x}_{i},i=1,\cdots,n$ are MC samples (i.e., simple random samples).

In the process of MC estimator, a crucial issue is to generate samples
$\bm{x}_{1},\bm{x}_{2},\cdots,\bm{x}_{n}$ from the density function
$f(\bm{x})$. This can be achieved using methods such as the inverse-
transformation method. This sampling method usually requires the target
distribution to have an explicit density function. However, in statistics,
particularly in Bayesian inference, the target density function is often
expressed as $f(\bm{x})=c\cdot p(\bm{x})$ where $c$ is an unknown constant or
with an extremely complex expression, in such cases, directly sampling from
the target distribution can be prohibited.

Thus the sampling importance resampling (SIR) method is one of the efficient
methods that can deal with the above issue. The basic process of this method
can be divided into the following three steps:

  * •

Step 1 (Sampling): Generating $N$ random samples
$\bm{x}_{1},\bm{x}_{2},\cdots,\bm{x}_{N}$ from the proposal distribution
$g(\bm{x})$;

  * •

Step 2 (Weighting): Calculating the importance weight:

$\omega_{i}=\frac{p(\bm{x}_{i})/g(\bm{x}_{i})}{\sum_{j=1}^{N}p(\bm{x}_{j})/g(\bm{x}_{j})},i=1,\cdots,N;$

  * •

Step 3 (Resampling): Sampling $n$ samples from the sample pool
$\\{\bm{x}_{i},i=1,\cdots,N\\}$ with the corresponding weight
$\omega_{i},i=1,\cdots,N$, denoting as
$\\{\bm{x}_{1}^{*},\cdots,\bm{x}_{n}^{*}\\}$.

According to Givens and Hoeting (2013), the samples
$\\{\bm{x}_{1}^{*},\cdots,\bm{x}_{n}^{*}\\}$ generating through the SIR method
asymptotically converge to target distribution $f(\bm{x})$ when
$N\rightarrow\infty$ and $n$ is fixed. Hence, once the samples
$\\{\bm{x}_{1}^{*},\cdots,\bm{x}_{n}^{*}\\}$ are obtained, the MC integration
in (2.1) can be calculated through the SIR samples instead of the MC samples,
and the corresponding estimator of integral $\mu$ can be denoted as
$\hat{\mu}_{SIR}$.

##  3 Two modified methods

In nature, step 3 of the SIR method illustrated in Section 2 is to generate
$n$ random samples from the following multinomial distribution
$M(N;\omega_{1},\omega_{2},\cdots,\omega_{N})$:

$\left(\begin{matrix}\bm{x}_{1}&\bm{x}_{2}&\cdots&\bm{x}_{N}\\\
\omega_{1}&\omega_{2}&\cdots&\omega_{N}\end{matrix}\right).$ (3.1)

For this discrete multinomial distribution, simple random sampling method can
be implemented to generate the samples via the inverse-transformation
operator. That is, firstly generating MC samples from the uniform
distribution, then the samples from this multinomial distribution can be found
via transforming the MC samples with the inverse cumulative distribution
function. However, the simple random sampling above will result in a large
variance (Ross, 2022; Robert and Casella, 2005; Fang and Wang, 1994).

###  3.1 SIR with antithetic samples

In this work, we first propose a modified SIR method to reduce the variance
with the help of the antithetic variable method. The antithetic variable
method is one of the widespread variance reduction techniques (Robert and
Casella, 2005). The basic idea of this method is to generate another estimator
$\hat{\mu}_{2}$ that is identically distributional but negatively
correlational with the given MC estimator $\hat{\mu}_{1}$ of the integral
$\mu$. Then by averaging these two estimators, the antithetic estimator is
obtained. That is, the key issue of the antithetic variable method is to
generate antithetic samples corresponding to the MC samples to estimate the
integral $\mu$.

Through embedding antithetic sampling into the process of the SIR method, a
new sampling method, named as “antithetic sampling importance resampling”
method (denoted as “Anti-SIR” for simplicity), is proposed. The first two
steps of the Anti-SIR method are the same as the SIR method. Step 3 of the new
method is to generate the antithetic samples from the sample pool
$\\{\bm{x}_{i},i=1,\cdots,N\\}$. Specifically,

  * •

Step 3’ (Antithetic resampling): Sampling $n$ antithetic samples
$\\{\bm{z}_{1},\bm{z}_{2},\cdots,\bm{z}_{n}\\}$ from the sample pool
$\\{\bm{x}_{i},i=1,\cdots,N\\}$ by using the antithetic variable method.

When the resampling size $n$ is even, $\bm{z}_{i}$ and $\bm{z}_{n/2+i}$ are
the corresponding antithetic samples for $i=1,\cdots,n/2$; When $n$ is odd,
one can first sample $n-1$ antithetic samples and then randomly sample the
last one.

####  3.1.1 The sampling details

In order to demonstrate the process of Anti-SIR method in detail, the size of
antithetic samples is assumed to be $2n$ in the following analysis. The key
issue of Anti-SIR method is to generate samples
$\\{\bm{z}_{1},\cdots,\bm{z}_{n},\widetilde{\bm{z}}_{1},\cdots,\widetilde{\bm{z}}_{n}\\}$
from the sample pool $\\{\bm{x}_{i},i=1,\cdots,N\\}$, where
$\widetilde{\bm{z}}_{i}$ is the antithetic sample correspond to $\bm{z}_{i}$
for $i=1,\cdots,n$. Thus, firstly $n$ random numbers $u_{1},\cdots,u_{n}$
following the uniform distribution U$[0,1]$ are generated. Through the
inverse-transformation method, we can get

$\bm{z}_{i}=\begin{cases}\bm{x}_{1},&\text{if }u_{i}\leq\omega_{1},\\\
\bm{x}_{2},&\text{if }\omega_{1}<u_{i}\leq\omega_{1}+\omega_{2},\\\ \cdots\\\
\bm{x}_{N},&\text{if }\sum_{k=1}^{N-1}\omega_{k}<u_{i}\leq 1,\end{cases}\quad
i=1,\cdots,n.$

Then let $\widetilde{u}_{i}=1-u_{i},i=1,\cdots,n$, the antithetic samples
corresponding to these $n$ samples $\bm{z}_{1},\cdots,\bm{z}_{n}$ can be
generated as follows:

$\widetilde{\bm{z}}_{i}=\begin{cases}\bm{x}_{1},&\text{if
}\widetilde{u}_{i}\leq\omega_{1},\\\ \bm{x}_{2},&\text{if
}\omega_{1}<\widetilde{u}_{i}\leq\omega_{1}+\omega_{2},\\\ \cdots\\\
\bm{x}_{N},&\text{if }\sum_{k=1}^{N-1}\omega_{k}<\widetilde{u}_{i}\leq
1,\end{cases}\quad i=1,\cdots,n.$

Following the framework of SIR method, the obtained samples
$\\{\bm{z}_{1},\cdots,\bm{z}_{n},\widetilde{\bm{z}}_{1},\cdots,\widetilde{\bm{z}}_{n}\\}$
via Anti-SIR method are distributed from the target distribution $f(\bm{x})$
asymptotically.

####  3.1.2 Theoretical properties

In this subsection, the statistical properties of $2n$ samples
$\\{\bm{z}_{1},\cdots,\bm{z}_{n},\widetilde{\bm{z}}_{1},\cdots,\widetilde{\bm{z}}_{n}\\}$
generating by Anti-SIR method are discussed.

######  Theorem 3.1.

When using SIR and Anti-SIR methods to estimate the integral $\mu=\int
h(\bm{x})f(\bm{x})d\bm{x}$ respectively, suppose the sample size to be $2n$,
the estimator corresponding to SIR method can be written as
$\hat{\mu}_{SIR}=\frac{1}{2n}\sum_{i=1}^{2n}h(\bm{x}_{i}^{*})$ and the
estimator corresponding to Anti-SIR method can be written as
$\hat{\mu}_{ASIR}=\frac{1}{2n}\sum_{i=1}^{n}[h(\bm{z}_{i})+h(\tilde{\bm{z}}_{i})]$,
we have

  * (i)

$E(\hat{\mu}_{ASIR})=E(\hat{\mu}_{SIR})$;

  * (ii)

$Var(\hat{\mu}_{ASIR})\leq Var(\hat{\mu}_{SIR})$.

Theorem 3.1 shows that Anti-SIR method has the same bias as SIR method and
both of them are unbiased estimators of integral $\mu$. In addition, by
embedding the antithetic sampling method into the original SIR method, the
variance of the new estimator $\hat{\mu}_{ASIR}$ can be reduced.

###  3.2 Latin hypercube sampling (LHS)-based SIR

It is well-known that Latin hypercube sampling (LHS) is another popular
variance reduction technique proposed by McKay et al. (1979), which is widely
used in computer experiments (Santner et al., 2018) and MC simulations (Ross,
2022). The basic idea of LHS method is to divide the domain along with each
dimension into sub-intervals with equal marginal probability, and sample once
from each stratum. In this way, the integration estimator via LHS can possess
lower variance than that via simple random sampling.

Hence, in this subsection, another variance-reduced SIR method based on Latin
hypercube sampling is proposed, denoted as “LHS-SIR” for convenience. Similar
to Anti-SIR method, the only difference between this modified method with the
original SIR method is the resampling operation. Specifically,

  * •

Step 3” (Latin hypercube resampling): Generating $n$ Latin hypercube samples
$\\{\bm{y}_{1},\bm{y}_{2},\cdots,\bm{y}_{n}\\}$ from the sample pool
$\\{\bm{x}_{i},i=1,\cdots,N\\}$ by using the LHS method.

Through this resampling process, it can be observed that the variance for the
integral estimator via the original SIR method can be reduced.

####  3.2.1 The sampling details

The pivotal element in the proposed LHS-SIR method is to generate samples from
the multinomial distribution (3.1) via the Latin hypercube sampling technique.
The detailed procedures are demonstrated as following three steps:

  * 1.

Dividing the interval $[0,1]$ into $n$ strata with equidistance, i.e.,

$\left[0,\frac{1}{n}\right],\left[\frac{1}{n},\frac{2}{n}\right],\cdots,\left[\frac{n-1}{n},1\right];$

  * 2.

Generating a random number $v_{i}$ in each sub-interval, i.e.,

$v_{i}\sim U\left[\frac{i-1}{n},\frac{i}{n}\right],i=1,\cdots,n;$

  * 3.

Finding $\bm{y}_{i}$ corresponding to each $v_{i}$ via the inverse-
transformation method, i.e.,

$\bm{y}_{i}=\begin{cases}\bm{x}_{1},&\text{if }v_{i}\leq\omega_{1},\\\
\bm{x}_{2},&\text{if }\omega_{1}<v_{i}\leq\omega_{1}+\omega_{2},\\\ \cdots\\\
\bm{x}_{N},&\text{if }\sum_{k=1}^{N-1}\omega_{k}<v_{i}\leq 1,\end{cases}\quad
i=1,\cdots,n.$

Hence based on the framework of SIR method, the obtained samples
$\\{\bm{y}_{i},i=1,\cdots,n\\}$ via LHS-SIR method asymptotically converge to
the target distribution $f(\bm{x})$.

####  3.2.2 Theoretical properties

In this subsection, the statistical properties of $n$ samples
$\\{\bm{y}_{1},\bm{y}_{2},\cdots,\bm{y}_{n}\\}$ generated by LHS-SIR method
are discussed. Similar to Anti-SIR method, LHS-SIR method also reduces the
variance of the estimation by generating the samples with a negative
correlation. Then, we have the following result.

######  Theorem 3.2.

When using SIR and LHS-SIR methods to estimate the integral $\mu=\int
h(\bm{x})f(\bm{x})d\bm{x}$ respectively, suppose the sample size to be $n$
here, the estimator corresponding to SIR method can be written as
$\hat{\mu}_{SIR}=\frac{1}{n}\sum_{i=1}^{n}h(\bm{x}_{i}^{*})$ and the estimator
corresponding to LHS-SIR method can be written as
$\hat{\mu}_{LSIR}=\frac{1}{n}\sum_{i=1}^{n}h(\bm{y}_{i})$, we have

  * (i)

$E(\hat{\mu}_{LSIR})=E(\hat{\mu}_{SIR})=\mu$;

  * (ii)

$Var(\hat{\mu}_{LSIR})\leq Var(\hat{\mu}_{SIR})$.

The results of Theorem 3.2 are similar to that of Theorem 3.1. By embedding
the Latin hypercube resampling method into the original SIR method, the
variance of the estimation can be reduced. In the next sections, we will
evaluate the performances of these methods by extensive numerical studies.

##  4 Numerical studies

In this section, some numerical simulations with different dimensions are
conducted to validate the performance of the modified SIR method.

Suppose $\bm{X}=(x_{1},x_{2},\cdots,x_{d})^{\top}$ is the $d$-dimensional
random variable and the density function is $f(\bm{x}),\bm{x}\in\mathcal{X}$,
where $\mathcal{X}$ is the corresponding sample space. Here, we want to
estimate the expectation of random variable $\bm{X}$:

$\mathrm{E}(\bm{X})=\int_{\mathcal{X}}\bm{x}f(\bm{x})\mathrm{d}\bm{x},$ (4.1)

by the corresponding samples mean with $n$ samples from target distribution
$f(\bm{x})$ generated through Anti-SIR and SIR methods. In order to alleviate
the randomness in the sampling process, the operation is replicated $K$ times
and the obtained samples each time are denoted as
$\bm{x}_{i}^{k},i=1,2,\cdots,n,k=1,2,\cdots,K$. Then $\mathrm{E}(\bm{X})$ can
be estimated by

$\hat{\mathrm{E}}_{k}(\bm{X})=\dfrac{1}{n}\sum_{i=1}^{n}\bm{x}_{i}^{k},k=1,2,\cdots,K,$

and the mean square error (MSE) of these estimations is given by

$\text{MSE}=\dfrac{1}{K}\sum_{k=1}^{K}\left[\hat{\mathrm{E}}_{k}(\bm{X})-\overline{\hat{\mathrm{E}}_{k}(\bm{X})}\right]^{2},$
(4.2)

where
$\overline{\hat{\mathrm{E}}_{k}(\bm{X})}=\dfrac{1}{K}\sum_{k=1}^{K}\hat{\mathrm{E}}_{k}(\bm{X})$.
If the true value of expectation (4.1) is known in advance, it could be used
to replace the term $\overline{\hat{\mathrm{E}}_{k}(\bm{X})}$ in (4.2). The
smaller the MSE value, the more accurate the estimation, indicating the better
the performance of the sampling method.

Example 1: Univariate distributions.

In this example, several cases in dimension 1 are considered. The target
distributions in this example are Beta distribution, normal distribution,
$t$-distribution, and $F$-distribution; and the proposal distributions are
uniform distribution, logistic distribution, Cauchy distribution, and inverse-
Gaussian distribution. The sample size chosen for all the cases is $N=20000$
and the resampling size is $n=1000$. All the simulations are repeated $K=1000$
times. The results obtained by SIR, Anti-SIR, and LHS-SIR methods are listed
in Table 1. It can be shown clearly that both the proposed Anti-SIR and LHS-
SIR methods perform better than the original SIR method in all cases. These
two modified methods can reduce the mean squared errors for the SIR method
with the help of antithetic resampling and Latin hypercube resampling.
However, in this example, neither of these two modified methods is
consistently better than the other.

Table 1: MSEs in the estimation of the mean $f(x)$ | $g(x)$ | SIR | Anti-SIR | LHS-SIR
---|---|---|---|---
Beta(2, 3) | U(0, 1) | $4.190\times 10^{-5}$ | $4.007\times 10^{-5}$ | $3.990\times 10^{-5}$
Beta(0.9, 0.9) | U(0, 1) | $1.037\times 10^{-4}$ | $9.613\times 10^{-5}$ | $9.023\times 10^{-5}$
N(0, 1) | Log(0, 1) | $1.144\times 10^{-3}$ | $1.070\times 10^{-3}$ | $1.053\times 10^{-3}$
N(0, 1) | C(0, 1) | $1.136\times 10^{-3}$ | $9.898\times 10^{-4}$ | $1.069\times 10^{-3}$
$t$(2, 0, 1) | C(0, 1) | $2.377\times 10^{-2}$ | $1.478\times 10^{-2}$ | $1.545\times 10^{-2}$
$F$(10, 6) | IG(1, 1) | $3.512\times 10^{-3}$ | $2.963\times 10^{-3}$ | $3.091\times 10^{-3}$

Example 2: Multivariate distribution.

In this example, the target distribution is the Kotz-type distribution which
is a generalization of the multivariate normal distribution (Kotz, 1975; Fang
and Zhang, 1990). Suppose the $d$ dimensional random variable $\bm{X}$ is
following the Kotz-type distribution with the density function:

$f(\bm{x})=C_{p}|\bm{\Sigma}|^{-1/2}\left[(\bm{x}-\bm{\beta})^{\top}\bm{\Sigma}^{-1}(\bm{x}-\bm{\beta})\right]^{M-1}\exp\left\\{-r\left[(\bm{x}-\bm{\beta})^{\top}\bm{\Sigma}^{-1}(\bm{x}-\bm{\beta})\right]^{s}\right\\},$

where $r>0,s>0,2M+d>2$, $\bm{\Sigma}$ is the $d\times d$ positive-definite
dispersion matrix, $\bm{\beta}=(\beta_{1},\cdots,\beta_{d})^{\top}$ is the
mean vector, and $C_{p}$ is the normalizing constant. Here, we consider a
4-dimensional Kotz-type distribution with the parameters:
$r=0.5,s=2,M=3,\bm{\beta}=(0,0,0,0)^{\top}$, and

$\Sigma=\begin{pmatrix}5.3&0.0&0.0&-0.2\\\ 0.0&4.0&-0.4&0.3\\\
0.0&-0.4&6.8&0.0\\\ -0.2&0.3&0.0&9.0\end{pmatrix}.$

In order to generate samples from the Kotz-type distribution through the SIR
and Anti-SIR methods, the multivariate normal distribution with mean vector
$\bm{\beta}_{g}=\bm{\beta}$ and covariance matrix
$\bm{\Sigma}_{g}=\bm{\Sigma}$ is taken as the proposal distribution. We set
$N=2000,n=400,K=1000$. The samples drawn by the SIR and Anti-SIR methods are
used to estimate the expectation of the Kotz-type distribution. Table 2 lists
the MSEs of the estimated components of the mean and the overall MSEs (OMSE).
It can be shown that for the MSEs of each component and the overall MSE, both
Anti-SIR and LHS-SIR methods are superior to SIR method. Hence, the proposed
methods in our work are effective sampling methods in this example. In this
case, LHS-SIR method performs a little better than Anti-SIR method.

Table 2: MSEs and OMSE in the estimation of the mean vector Methods | MSE ($\beta_{1}$) | MSE ($\beta_{2}$) | MSE ($\beta_{3}$) | MSE ($\beta_{4}$) | OMSE
---|---|---|---|---|---
SIR | 0.008897 | 0.006810 | 0.01227 | 0.01561 | 0.04359
Anti-SIR | 0.008950 | 0.006707 | 0.01170 | 0.01526 | 0.04262
LHS-SIR | 0.006330 | 0.004795 | 0.009020 | 0.01149 | 0.03163

##  5 Real data analysis

In this section, a real data analysis case is used to illustrate the
performance in Bayesian inference for our proposed method. Here, the data on
the number of coal-mining disasters per year between 1851 and 1962 is plotted
in Figure 1. It can be shown that the rate of accidents per year appears to
decrease around 1900, hence a change-point model is constructed to analyze
these data (Carlin et al., 1992).

Figure 1: The number of coal-mining disasters per year between 1851 and 1962

Let $X_{j}$ be the number of coal-mining disasters in year $j$, where we
denote $j=1$ in the year 1851 and index each year thereafter, i.e., $j=112$ in
the year 1962. Hence suppose the change-point occurs after the $\theta$year in
the series, $\theta\in\\{1,\cdots,111\\}$, then the change-point model can be
expressed as

$X_{1},\cdots,X_{\theta}\sim\mathrm{Poisson}(\lambda_{1}),$

and

$X_{\theta+1},\cdots,X_{112}\sim\mathrm{Poisson}(\lambda_{2}),$

where $\theta,\lambda_{1}$ and $\lambda_{2}$ are the unknown parameters and
the estimator of $\theta$ is our primary interest.

In order to estimate the parameters
$\Theta=\\{\theta,\lambda_{1},\lambda_{2}\\}$ for this model, suppose the
prior distribution for $\Theta$ is $\pi(\Theta)$, and the likelihood function
is:

$\ell(\Theta|\bm{x})=\prod_{i=1}^{\theta}\dfrac{\lambda_{1}^{x_{i}}}{x_{i}!}\mathrm{e}^{-\lambda_{1}}\prod_{j=\theta+1}^{112}\dfrac{\lambda_{2}^{x_{j}}}{x_{j}!}\mathrm{e}^{-\lambda_{2}}.$
(5.1)

Then with prior and the likelihood (5.1), the posterior for $\Theta$ can be
written as $f(\Theta|\bm{x})\propto\ell(\Theta|\bm{x})\pi(\Theta)$. In this
work, the SIR, Anti-SIR, and LHS-SIR methods are used to draw samples from
this posterior distribution $f(\Theta|\bm{x})$. That is, the posterior is the
target density. Let the prior distribution $\pi(\Theta)$ sever as the proposal
distribution, i.e., sampling $\Theta_{1},\cdots,\Theta_{N}$ i.i.d. from
$\pi(\Theta)$, then the $i$th unstandardized weight equals
$\ell(\Theta_{i}|\bm{x})$, $i=1,\cdots,N$. Then the i.i.d resampling,
antithetic resampling, and Latin hypercube resampling samples from
$\Theta_{1},\cdots,\Theta_{N}$ can be employed to estimate the unknown
parameters respectively.

Here, we set two different priors for a Bayesian analysis of this model:

  * •

Case 1: $\theta$ is from the discrete uniform distribution on
$\\{1,2,\cdots,111\\}$, $\lambda_{i}|a_{i}\sim\mathrm{Gamma}(3,a_{i})$ and
$a_{i}\sim\mathrm{Gamma}(10,10)$ for $i=1,2$;

  * •

Case 2: $\theta$ is from the discrete uniform distribution on
$\\{1,2,\cdots,111\\}$, $\lambda_{2}=\alpha\lambda_{1}$,
$\lambda_{1}|a\sim\mathrm{Gamma}(3,a)$, $a\sim\mathrm{Gamma}(10,10)$ and
$\log\alpha\sim\mathrm{U}(\log 1/8,\log 2)$.

where “Gamma” represents the Gamma distribution and “U” is the uniform
distribution. The sample size chosen for these two cases is $N=5000$ and the
resample size is $n=2000$. Repeating the sampling processes $K=1000$ times,
the results about the estimated parameters are listed in Table 3. It can be
shown that in both cases, the estimated parameters through SIR, Anti-SIR, and
LHS-SIR methods are consistent on average. However, the standard deviations
obtained by Anti-SIR and LHS-SIR methods are a little less than those
corresponding to SIR method. That is, these two proposed methods scan provide
more robust results than that of SIR method in statistical inference. In
addition, LHS-SIR method produces better performance than Anti-SIR method
except for the estimation of parameter $\lambda_{1}$ in case 1.

Table 3: The estimated parameters of the change-point model

Cases | Methods | Parameters
---|---|---
$\theta$ | $\lambda_{1}$ | $\lambda_{2}$
1 | SIR | 39.81 (0.8862) | 3.122 (0.1103) | 0.9569 (0.04513)
Anti-SIR | 39.81 (0.8831) | 3.121 (0.1101) | 0.9569 (0.04503)
| LHS-SIR | 39.81 (0.8828) | 3.122 (0.1103) | 0.9569 ( 0.04494)
2 | SIR | 39.97 (0.5880) | 3.116 (0.07230) | 0.9243 (0.02926)
Anti-SIR | 39.97 (0.5876) | 3.116 (0.07198) | 0.9241 (0.02915)
| LHS-SIR | 39.98 (0.5850) | 3.116 (0.07173) | 0.9242 ( 0.02904)

  * ∗ The results are listed as the form “mean (standard deviations)”.

##  6 Conclusion

In this paper, two modified SIR methods, namely Anti-SIR and LHS-SIR methods,
are proposed based on the variance reduction techniques. The core step is to
replace the simple random resampling from the sample pool
$\\{\bm{x}_{1},\cdots,\bm{x}_{N}\\}$ with the antithetic resampling and Latin
hypercube resampling respectively. The variance reduction effects for these
two proposed methods are demonstrated in both theoretical and numerical
aspects. It can be observed that the new proposed methods can produce lower
estimation errors than the original SIR method. In addition, from numerical
results, LHS-SIR method performs a little better than Anti-SIR method,
especially in the high-dimensional case.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China
[grant numbers: 12371261 and 12171187] and the Humanities and Social Science
Foundation of Ministry of Education of China [grant number: 23YJC910004].

## Declarations

The authors have no relevant financial or non-financial interests to disclose.

## References

  * Carlin et al. (1992) Carlin, B. P., A. E. Gelfand, and A. F. Smith (1992), Hierarchical Bayesian analysis of changepoint problems, Journal of the Royal Statistical Society: Series C (Applied Statistics), 41(2), 389–405, doi:10.2307/2347570.
  * Fang and Wang (1994) Fang, K.-T., and Y. Wang (1994), Number-theoretic Methods in Statistics, CRC Press, London.
  * Fang and Zhang (1990) Fang, K.-T., and Y.-T. Zhang (1990), Generalized Multivariate Analysis, Springer-Verlag, Berlin.
  * Givens and Hoeting (2013) Givens, G. H., and J. A. Hoeting (2013), Computational Statistics, John Wiley & Sons, New Jersey.
  * Givens and Raftery (1996) Givens, G. H., and A. E. Raftery (1996), Local adaptive importance sampling for multivariate densities with strong nonlinear relationships, Journal of the American Statistical Association, 91(433), 132–141, doi:10.1080/01621459.1996.10476670.
  * Kotz (1975) Kotz, S. (1975), Multivariate distributions at a cross road, in Modern Course on Statistical Distributions in Scientific Work, edited by G. P. Patil, S. Kotz, and J. K. Ord, pp. 247–270, Springer, Dordrecht, the Netherlands.
  * McKay et al. (1979) McKay, M. D., R. J. Beckman, and W. J. Conover (1979), A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics, 21(2), 239–245, doi:10.1080/00401706.2000.10485979.
  * Ning and Tao (2020) Ning, J., and H. Tao (2020), Randomized quasi-random sampling/importance resampling, Communications in Statistics-Simulation and Computation, 49(12), 3367–3379, doi:10.1080/03610918.2018.1547398.
  * Perez et al. (2005) Perez, C., J. Martin, M. Rufo, and C. Rojano (2005), Quasi-random sampling importance resampling, Communications in Statistics-Simulation and Computation, 34(1), 97–112, doi:10.1081/SAC-200047112.
  * Robert and Casella (2005) Robert, C. P., and G. Casella (2005), Monte Carlo Statistical Methods, Springer, New York.
  * Ross (2022) Ross, S. M. (2022), Simulation (6th Edition), Academic Press, New Jersey.
  * Rubin (1987) Rubin, D. B. (1987), The calculation of posterior distributions by data augmentation: Comment: A noniterative sampling/importance resampling alternative to the data augmentation algorithm for creating a few imputations when fractions of missing information are modest: The sir algorithm, Journal of the American Statistical Association, 82(398), 543, doi:10.1080/01621459.1987.10478461.
  * Santner et al. (2018) Santner, T. J., B. J. Williams, and W. I. Notz (2018), The Design and Analysis of Computer Experiments (Second Edition), Springer, New York.
  * Skare et al. (2003) Skare, O., E. Bolviken, and L. Holden (2003), Improved sampling-importance resampling and reduced bias importance sampling, Scandinavian Journal of Statistics, 30(4), 719–737, doi:10.1111/1467-9469.00360.
  * Tokdar and Kass (2010) Tokdar, S. T., and R. E. Kass (2010), Importance sampling: a review, Wiley Interdisciplinary Reviews: Computational Statistics, 2(1), 54–60, doi:10.1002/wics.56.
  * Vandewoestyne and Cools (2010) Vandewoestyne, B., and R. Cools (2010), On the convergence of quasi-random sampling/importance resampling, Mathematics and Computers in Simulation, 81(3), 490–505, doi:10.1016/j.matcom.2009.09.004.
  * Yi et al. (2023) Yi, S.-Y., Z. Liu, M.-Q. Liu, and Y.-D. Zhou (2023), Global likelihood sampler for multimodal distributions, Journal of Computational and Graphical Statistics, 32(3), 927–937, doi:10.1080/10618600.2023.2165499.

##  Appendix A The proof of the theorems

###  A.1 The proof of Theorem 1

In order to prove the theorem 1, a significant lemma is introduced.

######  Lemma A.1 (Corollary in Ross (2022)).

If $l(x_{1},\cdots,x_{m})$ is a monotone function of each of its arguments,
then, for a set $U_{1},\cdots,U_{m}$ of independent random numbers drawn from
the uniform distribution,

$Cov(l(U_{1},\cdots,U_{m}),l(1-U_{1},\cdots,1-U_{m}))\leq 0.$

This lemma indicates that the samples obtained by the antithetic sampling
method are negatively correlated. This property plays a key role in reducing
the variance of samples by the Anti-SIR method. Formally, Theorem 1 is proved
as follows.

Proof: (i) For SIR method, according to page 183 of Givens and Hoeting (2013),

$E(\hat{\mu}_{SIR})=E\left\\{\frac{1}{2n}\sum_{i=1}^{2n}h(\bm{x}_{i}^{*})\right\\}=E(h(\bm{x}_{1}^{*})).$

For Anti-SIR method, antithetic samples
$\\{\bm{z}_{1},\cdots,\bm{z}_{n},\tilde{\bm{z}}_{1},\cdots,\tilde{\bm{z}}_{n}\\}$,
$\bm{z}_{i}$ and $\tilde{\bm{z}_{i}}$ are correlated, but pairs
$(\bm{z}_{i},\tilde{\bm{z}_{i}})$, $i=1,\cdots,n,$ are independent, hence

$E(\hat{\mu}_{ASIR})=E\left\\{\frac{1}{2n}\sum_{i=1}^{n}[h(\bm{z}_{i})+h(\bm{\tilde{}}{z}_{i})]\right\\}=\frac{1}{2}E[h(\bm{z}_{1})+h(\tilde{\bm{z}}_{1})]=E(h(\bm{z}_{1}))=E(\hat{\mu}_{SIR}).$

(ii) From (i), we can have

$Var(\hat{\mu}_{SIR})=Var\left\\{\frac{1}{2n}\sum_{i=1}^{2n}h(\bm{x}_{i}^{*})\right\\}=\frac{1}{2n}Var(h(\bm{x}_{1}^{*})).$

Then,

$\begin{split}Var(\hat{\mu}_{ASIR})&=Var\left\\{\frac{1}{2n}\sum_{i=1}^{n}[h(\bm{z}_{i})+h(\bm{\tilde{}}{z}_{i})]\right\\}=\frac{1}{4n}Var[h(\bm{z}_{1})+h(\bm{\tilde{}}{z}_{1})]\\\
&=\frac{1}{4n}[Var(h(\bm{z}_{1}))+Var(h(\bm{\tilde{}}{z}_{1}))+2Cov(h(\bm{z}_{1}),h(\bm{\tilde{}}{z}_{1}))]\\\
&=\frac{1}{2n}Var(h(\bm{z}_{1})+\frac{1}{2n}Cov(h(\bm{z}_{1}),h(\bm{\tilde{}}{z}_{1}))\\\
&=Var(\hat{\mu}_{SIR})+\frac{1}{2n}Cov(h(\bm{z}_{1}),h(\bm{\tilde{}}{z}_{1}).\end{split}$

According to Lemma A.1, because the inverse cumulative distribution function
of multinomial distribution is monotone,
$Cov(h(\bm{z}_{1}),h(\bm{\tilde{}}{z}_{1}))\leq 0$. Hence,
$Var(\hat{\mu}_{ASIR})\leq Var(\hat{\mu}_{SIR})$. The proof is finished.

###  A.2 The proof of Theorem 2

Proof: (i) Because LH samples $\\{\bm{y}_{i},i=1,\cdots,2\\}$ are negatively
correlated each other, we have

$E(\hat{\mu}_{LSIR})=E\left(\frac{1}{n}\sum_{i=1}^{n}h(\bm{y}_{i})\right)=E(h(\bm{y}_{1}))=E(\hat{\mu}_{SIR}).$

(ii) From (i), we can get that

$\begin{split}Var(\hat{\mu}_{LSIR})&=Var\left(\frac{1}{n}\sum_{i=1}^{n}h(\bm{y}_{i})\right)\\\
&=\frac{1}{n^{2}}\left[\sum_{i}Var(h(\bm{y}_{i}))+2\sum_{i,j}Cov(h(\bm{y}_{i}),h(\bm{y}_{j}))\right]\\\
&=\frac{1}{n}Var(h(\bm{y}_{1}))+\frac{2}{n^{2}}\sum_{i,j}Cov(h(\bm{y}_{i}),h(\bm{y}_{j}))\\\
&=Var(\hat{\mu}_{SIR})+\frac{2}{n^{2}}\sum_{i,j}Cov(h(\bm{y}_{i}),h(\bm{y}_{j})).\end{split}$

According to the results in McKay et al. (1979), we have
$Cov(h(\bm{y}_{i}),h(\bm{y}_{j}))\leq 0$ for any $i,j\in\\{1,\cdots,n\\}$.
Hence, we can obtain $Var(\hat{\mu}_{LSIR})\leq Var(\hat{\mu}_{SIR})$. The
proof is finished.


$\mathbf{R}(\bm{\theta},n)\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{r}\left(\bm{\theta},n\right)$
exist and are continuous in ${\bm{\theta}}\in\bm{\Theta}$. Moreover, we
require that for all $\mathbf{u}\in{\rm I\\!R}^{p}$ with
$||\mathbf{u}||_{2}=1$, there exists a sequence $(\delta_{n})\in\mathfrak{D}$
such that

$\scalebox{0.98}{$\sqrt{n}\left(1+\frac{1}{H}\right)^{-\nicefrac{{1}}{{2}}}\mathbf{u}^{T}\bm{\Sigma}(\bm{\theta}_{0},n)^{-\nicefrac{{1}}{{2}}}\Big{[}\mathbf{B}(\bm{\theta}_{0},n)-\mathbf{B}\left({\bm{\theta}}^{(\mathbf{a},\mathbf{r})},n\right)\Big{]}\left(\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}\right)=\mathcal{O}_{\rm
p}(\delta_{n})$},$ (G.2)

where

$\mathbf{B}(\bm{\theta},n)\equiv\mathbf{I}+\mathbf{A}(\bm{\theta})+\mathbf{L}(n)+\mathbf{R}(\bm{\theta},n),$
$\mathbf{B}\big{(}{\bm{\theta}}^{(\mathbf{a},\mathbf{r})},n\big{)}\equiv\mathbf{I}+\mathbf{A}({\bm{\theta}}^{(\mathbf{a})})+\mathbf{L}(n)+\mathbf{R}\left({\bm{\theta}}^{(\mathbf{r})},n\right).$

while $\mathbf{A}({\bm{\theta}}^{(\mathbf{a})})$ and
$\mathbf{R}({\bm{\theta}}^{(\mathbf{r})})$ are defined using notation in
(F.1).

Assumption D5 allows us to quantify how “far” the matrices
$\mathbf{B}(\bm{\theta}_{0},n)$ and
$\mathbf{B}\big{(}\hat{\bm{\theta}}_{(n,H)},n\big{)}$ are from each other. In
the case where $p\to c<\infty$ as $n\to\infty$, the consistency of
$\hat{\bm{\theta}}_{(n,H)}$ and the continuity of
$\mathbf{B}(\bm{\theta})\equiv\displaystyle\lim_{n\to\infty}\mathbf{B}(\bm{\theta},n)$
in $\bm{\theta}_{0}$ would be sufficient so that Assumption D5 would not be
needed to establish the asymptotic distribution of
$\hat{\bm{\theta}}_{(n,H)}$. However, when $p\to\infty$ as $n\to\infty$ the
requirement in (G.2) can be strong and difficult to verify for a specific
model and auxiliary estimator. Having stated this, the following proposition
defines the distribution of the IB-estimator.

######  Proposition 4:

Under Assumptions A, B, C∗ and D5, for all $\mathbf{u}\in{\rm I\\!R}^{p}$ such
that $||\mathbf{u}||_{2}=1$, there exist a sample size
$n^{\ast}\in\mathbb{N}^{\ast}$ and a sequence
$\left(\delta_{n}\right)\in\mathfrak{D}$ such that for all
$n\in\mathbb{N}^{\ast}$ satisfying $n\geq n^{\ast}$ we have

$\scalebox{0.99}{$\sqrt{n}\left(1+\frac{1}{H}\right)^{-\nicefrac{{1}}{{2}}}\mathbf{u}^{T}\bm{\Sigma}(\bm{\theta}_{0},n)^{-\nicefrac{{1}}{{2}}}\mathbf{B}(\bm{\theta}_{0},n)\left(\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}\right)\overset{d}{=}Z+\delta_{n}\mathcal{O}_{\rm
p}\left(\max\left(1,\frac{p^{2}}{\sqrt{H}}\right)\right)$},$

where $Z\sim\mathcal{N}(0,1)$.

Proof: Without loss of generality we can assume that the sample sizes
$n^{*}\in\mathbb{N}^{*}$ of Assumption C∗ and D5 are the same. Let
$n\in\mathbb{N}^{*}$ be such that $n\geq n^{*}$. By the definition of the IB-
estimator we have

$\hat{\bm{\pi}}\left({\bm{\theta}_{0}},n,{\bm{\omega}_{0}}\right)=\frac{1}{H}\sum_{h=1}^{H}\hat{\bm{\pi}}\left({\hat{\bm{\theta}}_{(n,H)}},n,\bm{\omega}_{h}\right).$
(G.3)

Since the auxiliary estimator may be expressed, for any
$({\bm{\theta}},n,\bm{\omega})\in\bm{\Theta}\times\mathbb{N}^{*}\times\bm{\Omega}$,
as

$\hat{\bm{\pi}}\left({\bm{\theta}},n,{\bm{\omega}}\right)=\bm{\theta}+\mathbf{a}\left({\bm{\theta}}\right)+\mathbf{c}(n)+\mathbf{L}_{n}\bm{\theta}+\mathbf{r}\left({\bm{\theta}},n\right)+\mathbf{v}\left({\bm{\theta}},n,{\bm{\omega}}\right).$
(G.4)

we deduce from (G.3), by using the mean value theorem on
$\mathbf{a}({\bm{\theta}})$ and $\mathbf{r}({\bm{\theta}},n)$, that

$\mathbf{B}\left({\bm{\theta}}^{(\mathbf{a},\mathbf{r})},n\right)\bm{\Delta}=\mathbf{v}\left({\bm{\theta}_{0}},n,{\bm{\omega}_{0}}\right)-\frac{1}{H}\sum_{h=1}^{H}\mathbf{v}\left({\hat{\bm{\theta}}_{(n,H)}},n,\bm{\omega}_{h}\right),$
(G.5)

where
$\mathbf{B}\left({\bm{\theta}}^{(\mathbf{a},\mathbf{r})},n\right)\equiv\mathbf{I}+\mathbf{A}({\bm{\theta}}^{(\mathbf{a})})+\mathbf{L}(n)+\mathbf{R}\left({\bm{\theta}}^{(\mathbf{r})},n\right)$
and $\bm{\Delta}\equiv\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}$. Applying the
mean value theorem on
$\mathbf{v}\left({\hat{\bm{\theta}}_{(n,H)}},n,{\bm{\omega}_{h}}\right)$ and
using the notation defined in (F.1), we obtain for all $h=1,\dots,H$

$\mathbf{v}\left({\hat{\bm{\theta}}_{(n,H)}},n,{\bm{\omega}_{h}}\right)=\mathbf{v}\left({\bm{\theta}_{0}},n,{\bm{\omega}_{h}}\right)+\mathbf{V}\left({\bm{\theta}}^{(\mathbf{v})},n,\bm{\omega}_{h}\right)\bm{\Delta}.$

Hence, let $\mathbf{u}\in{\rm I\\!R}^{p}$ be such that $||\mathbf{u}||_{2}=1$,
then by Assumption C∗ we have,

$\sqrt{n}\mathbf{u}^{T}\bm{\Sigma}(\bm{\theta}_{0},n)^{-\nicefrac{{1}}{{2}}}\mathbf{v}\left({\bm{\theta}_{0}},n,{\bm{\omega}_{0}}\right)=Z(\bm{\theta}_{0},n,\bm{\omega}_{0},\mathbf{u})+\delta_{n}^{(0)}W(\bm{\theta}_{0},n,\bm{\omega}_{0},\mathbf{u})$

and for all $h=1,\dots,H$

$\sqrt{n}\mathbf{u}^{T}\bm{\Sigma}(\bm{\theta}_{0},n)^{-\nicefrac{{1}}{{2}}}\mathbf{v}\left({\bm{\theta}_{0}},n,{\bm{\omega}_{h}}\right)=Z(\bm{\theta}_{0},n,\bm{\omega}_{h},\mathbf{u})+\delta_{n}^{(h)}W(\bm{\theta}_{0},n,\bm{\omega}_{h},\mathbf{u})$

where
$(\delta_{n}^{(0)}),(\delta_{n}^{(1)}),\dots,(\delta_{n}^{(H)})\in\mathfrak{D}$
and

$\displaystyle
Z(\bm{\theta}_{0},n,\bm{\omega}_{0},\mathbf{u})\sim\mathcal{N}(0,1),\;\;\;\;\;\;Z(\bm{\theta}_{0},n,\bm{\omega}_{h},\mathbf{u})\sim\mathcal{N}(0,1),$
$\displaystyle
W(\bm{\theta}_{0},n,\bm{\omega}_{0},\mathbf{u})=\mathcal{O}_{\rm
p}(1),\;\;\;\;\;\;W(\bm{\theta}_{0},n,\bm{\omega}_{h},\mathbf{u})=\mathcal{O}_{\rm
p}(1).$

Without loss of generality, we can suppose that
$(\delta_{n})\equiv(\delta_{n}^{(0)})=(\delta_{n}^{(1)})=\dots=(\delta_{n}^{(H)})$
as one may simply modify $W(\bm{\theta}_{0},n,\bm{\omega}_{0},\mathbf{u})$ and
$W(\bm{\theta}_{0},n,\bm{\omega}_{h},\mathbf{u})$ for all $h=1,\dots,H$. For
simplicity, we write
$\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}\equiv\sqrt{n}\mathbf{u}^{T}\bm{\Sigma}(\bm{\theta}_{0},n)^{-\nicefrac{{1}}{{2}}}$,
$Z_{0}\equiv Z(\bm{\theta}_{0},n,\bm{\omega}_{0},\mathbf{u})$ and $Z_{h}\equiv
Z(\bm{\theta}_{0},n,\bm{\omega}_{h},\mathbf{u})$. Multiplying the right hand
side of (G.5) by $\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}$, we have

$\displaystyle\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}\left(\mathbf{v}\left({\bm{\theta}},n,{\bm{\omega}_{0}}\right)-\frac{1}{H}\sum_{h=1}^{H}\mathbf{v}\left({\hat{\bm{\theta}}_{(n,H)}},n,\bm{\omega}_{h}\right)\right)$
(G.6) $\displaystyle=Z_{0}+\mathcal{O}_{\rm
p}(\delta_{n})-\dfrac{1}{H}\sum_{h=1}^{H}\big{(}Z_{h}+\mathcal{O}_{\rm
p}(\delta_{n})\big{)}-\dfrac{1}{H}\sum_{h=1}^{H}\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}\mathbf{V}\left({\bm{\theta}}^{(\mathbf{v})},n,\bm{\omega}_{h}\right)\bm{\Delta}$

The first three terms of (G.6) are considered separately in the following
computation

$\displaystyle Z_{0}+\mathcal{O}_{\rm
p}(\delta_{n})-\dfrac{1}{H}\sum_{h=1}^{H}\big{(}Z_{h}+\mathcal{O}_{\rm
p}(\delta_{n})\big{)}$
$\displaystyle=Z_{0}-\dfrac{1}{H}\sum_{h=1}^{H}Z_{h}+\mathcal{O}_{\rm
p}(\delta_{n})-\dfrac{1}{H}\sum_{h=1}^{H}\mathcal{O}_{\rm p}(\delta_{n})$
(G.7)
$\displaystyle\overset{d}{=}\left(1+\frac{1}{H}\right)^{\nicefrac{{1}}{{2}}}Z+\mathcal{O}_{\rm
p}(\delta_{n})-\mathcal{O}_{\rm p}\left(\frac{\delta_{n}}{\sqrt{H}}\right)$
$\displaystyle=\left(1+\frac{1}{H}\right)^{\nicefrac{{1}}{{2}}}Z+\mathcal{O}_{\rm
p}(\delta_{n}),$

where $Z\sim\mathcal{N}\left(0,1\right)$. We now consider the last term of
(G.6). Since
$\bm{\Sigma}({\bm{\theta}},n)^{-\nicefrac{{1}}{{2}}}=\mathcal{O}(1)$ and
$\mathbf{V}\left({\bm{\theta}}^{(\mathbf{v})},n,\bm{\omega}_{h}\right)=\mathcal{O}_{p}\left(n^{-\nicefrac{{1}}{{2}}}\right)$
elementwise by Assumption C∗ and $||\mathbf{u}||_{2}=1$ we have for all
$h=1,\dots,H$

$\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}\mathbf{V}\left({\bm{\theta}}^{(\mathbf{v})},n,\bm{\omega}_{h}\right)\bm{\Delta}=\mathcal{O}_{\rm
p}(p)\mathbf{1}_{p}^{T}\bm{\Delta},$

where $\mathbf{1}_{p}\in{\rm I\\!R}^{p}$ is a vector of ones. Since
Assumptions A, B, C∗ and D5 imply that $\hat{\bm{\theta}}_{(n,H)}$ is
consistent estimator of $\bm{\theta}_{0}$ (see Figure 4), there exists a
sequence $(\delta^{\prime}_{n})\in\mathfrak{D}$ such that
$\bm{\Delta}=\mathcal{O}_{\rm p}(\delta^{\prime}_{n})$ elementwise for all
$n\in\mathbb{N}^{*}$. Without loss of generality, we may assume that
$(\delta^{\prime}_{n})=(\delta_{n})$. Therefore, we have

$\displaystyle\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}\mathbf{V}\left({\bm{\theta}}^{(\mathbf{v})},n,\bm{\omega}_{h}\right)\bm{\Delta}$
$\displaystyle=\mathcal{O}_{\rm
p}(p)\mathbf{1}_{p}^{T}\bm{\Delta}=\mathcal{O}_{\rm
p}(p)\mathbf{1}_{p}^{T}\mathcal{O}_{\rm p}(\delta_{n})\mathbf{1}_{p}$
$\displaystyle=\mathcal{O}_{\rm
p}(p\delta_{n})\mathbf{1}_{p}^{T}\mathbf{1}_{p}=\mathcal{O}_{\rm
p}(p^{2}\delta_{n}),$

which leads to the computation for the last term of (G.6),

$\dfrac{1}{H}\sum_{h=1}^{H}\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}\mathbf{V}\left({\bm{\theta}}^{(\mathbf{v})},n,\bm{\omega}_{h}\right)\bm{\Delta}=\dfrac{1}{H}\sum_{h=1}^{H}\mathcal{O}_{\rm
p}(p^{2}\delta_{n})=\mathcal{O}_{\rm
p}\left(\delta_{n}\frac{p^{2}}{\sqrt{H}}\right).$ (G.8)

Combining (G.5), (G.6) and (G.8), we obtain

$\displaystyle\left(1+\frac{1}{H}\right)^{-\nicefrac{{1}}{{2}}}\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}\mathbf{B}\left({\bm{\theta}}^{(\mathbf{a},\mathbf{r})},n\right)\bm{\Delta}$
$\displaystyle=$
$\displaystyle\left(1+\frac{1}{H}\right)^{-\nicefrac{{1}}{{2}}}\mathfrak{X}_{\bm{\theta}_{0},n,\mathbf{u}}\left(\mathbf{v}\left({\bm{\theta}},n,{\bm{\omega}_{0}}\right)-\frac{1}{H}\sum_{h=1}^{H}\mathbf{v}\left({\hat{\bm{\theta}}_{(n,H)}},n,\bm{\omega}_{h}\right)\right)$
$\displaystyle\overset{d}{=}\left(1+\frac{1}{H}\right)^{-\nicefrac{{1}}{{2}}}\left(\left(1+\frac{1}{H}\right)^{\nicefrac{{1}}{{2}}}Z+\mathcal{O}_{\rm
p}(\delta_{n})+\mathcal{O}_{\rm
p}\left(\delta_{n}\frac{p^{2}}{\sqrt{H}}\right)\right)$
$\displaystyle\overset{d}{=}Z+\mathcal{O}_{\rm
p}\left(\delta_{n}\sqrt{\frac{H}{H+1}}\right)+\mathcal{O}_{\rm
p}\left(\delta_{n}\frac{p^{2}}{\sqrt{H+1}}\right)$
$\displaystyle=Z+\delta_{n}\mathcal{O}_{\rm
p}\left(1\right)+\delta_{n}\mathcal{O}_{\rm
p}\left(\frac{p^{2}}{\sqrt{H}}\right)=Z+\delta_{n}\mathcal{O}_{\rm
p}\left(\max\left(1,\frac{p^{2}}{\sqrt{H}}\right)\right).$

The proof is completed using (G.2) of Assumption D5. ∎

######  Remark D:

From the proof of Proposition 4, we have that

$\scalebox{0.99}{$\sqrt{n}\left(1+\frac{1}{H}\right)^{-\nicefrac{{1}}{{2}}}\mathbf{u}^{T}\bm{\Sigma}(\bm{\theta}_{0},n)^{-\nicefrac{{1}}{{2}}}\mathbf{B}(\bm{\theta}^{(\mathbf{a},\mathbf{r})},n)\left(\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}\right)\overset{d}{=}Z+\delta_{n}\mathcal{O}_{\rm
p}\left(\max\left(1,\frac{p^{2}}{\sqrt{H}}\right)\right)$},$ (G.9)

without using (G.2) in Assumption D5. As a result, Assumption D5 is not
strictly necessary for the Gaussian approximation of the distribution of
$\hat{\bm{\theta}}_{(n,H)}$. However, while this assumption is strong, it is
quite convenient to deliver a reasonable approximation of
$\mathbf{B}(\bm{\theta}^{(\mathbf{a},\mathbf{r})},n)$ in (G.9) when $p$ is
allowed to diverge, since its applicability is quite limited in practice.

##  Appendix H Proofs of the Main Results

The results and assumptions presented so far can now be combined to deliver
the three theorems of this article. While Theorem 1 and 2 were stated in the
main text, Theorem 3, which is stated below, presents the asymptotic
distribution of the IB-estimator. Theorem 3 is based on Assumption D∗∗, which
is the union of Assumption D∗ and Assumption D5. Assumption D∗∗ is provided
below for clarity of exposition.

######  Assumption D∗∗:

The bias function $\mathbf{d}\left(\bm{\theta},n\right)$ is such that:

  1. 1.

The asymptotic bias function $\mathbf{a}(\bm{\theta})$ can be written as

$\mathbf{a}({\bm{\theta}})=\mathbf{M}{\bm{\theta}}+\mathbf{s},$

where $\mathbf{M}\in{\rm I\\!R}^{p\times p}$ with $||\mathbf{M}||_{F}<1$ and
$\mathbf{s}\in{\rm I\\!R}^{p}$.

  2. 2.

There exists a sample size $n^{*}\in\mathbb{N}^{*}$ such that for all
$n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$ the matrix
$(\mathbf{M}+\mathbf{L}(n)+\mathbf{I})^{-1}$ exists.

  3. 3.

There exist real $\beta,\gamma>0$ such that for all
$\bm{\theta}\in\bm{\Theta}$ and any $j,l=1,\dots,p$, we have

$\displaystyle\mathbf{L}_{j,l}(n)=\mathcal{O}(n^{-\beta}),\;\;\;\mathbf{r}_{j}\left(\bm{\theta},n\right)=\mathcal{O}(n^{-\gamma}),\;\;\;\lim_{n\to\infty}\;\frac{p^{\nicefrac{{3}}{{2}}}}{n^{\beta}}=0\;\;\;\text{and}\;\;\;\lim_{n\to\infty}\;\frac{p^{2}}{n^{\gamma}}=0.$

  4. 4.

Defining $c_{n}\equiv\displaystyle{\max_{j=1,\dots,p}}\mathbf{c}_{j}(n)$ for
all $n\in\mathbb{N}^{*}$, we require that the sequence
$\left\\{c_{n}\right\\}_{n\in\mathbb{N}^{*}}$ is such that

$\lim_{n\to\infty}\;p^{\nicefrac{{1}}{{2}}}c_{n}=0.$

  5. 5.

There exists a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ such that for all
$n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$ the Jacobian matrices
$\mathbf{A}(\bm{\theta})\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{a}\left(\bm{\theta}\right)$
and
$\mathbf{R}(\bm{\theta},n)\equiv\frac{\partial}{\partial\,{\bm{\theta}}^{T}}\mathbf{r}\left(\bm{\theta},n\right)$
exist and are continuous in ${\bm{\theta}}\in\bm{\Theta}$. Moreover, we
require that for all $\mathbf{u}\in{\rm I\\!R}^{p}$ with
$||\mathbf{u}||_{2}=1$, there exists a sequence $(\delta_{n})\in\mathfrak{D}$
such that

$\scalebox{0.99}{$\sqrt{n}\left(1+\frac{1}{H}\right)^{-\nicefrac{{1}}{{2}}}\mathbf{u}^{T}\bm{\Sigma}(\bm{\theta}_{0},n)^{-\nicefrac{{1}}{{2}}}\Big{[}\mathbf{B}(\bm{\theta}_{0},n)-\mathbf{B}\left({\bm{\theta}}^{(\mathbf{a},\mathbf{r})},n\right)\Big{]}\left(\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}\right)=\mathcal{O}_{\rm
p}(\delta_{n})$},$

where

$\mathbf{B}(\bm{\theta},n)\equiv\mathbf{I}+\mathbf{A}(\bm{\theta})+\mathbf{L}(n)+\mathbf{R}(\bm{\theta},n)$

and

$\mathbf{B}\big{(}{\bm{\theta}}^{(\mathbf{a},\mathbf{r})},n\big{)}\equiv\mathbf{I}+\mathbf{A}({\bm{\theta}}^{(\mathbf{a})})+\mathbf{L}(n)+\mathbf{R}\left({\bm{\theta}}^{(\mathbf{r})},n\right).$

Using this new assumption, we now restate Theorem 1 and 2, and state Theorem
3, followed by their proof.

######  Theorem 1:

Under Assumptions A, B, C and D, for all $H\in\mathbb{N}^{\ast}$,

  1. 1.

There exist a $n^{*}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}$
with $n\geq n^{*}$,
$\left\\{\hat{\bm{\theta}}_{(n,H)}\right\\}=\widehat{\bm{\Theta}}_{(n,H)}$,
i.e. the set $\widehat{\bm{\Theta}}_{(n,H)}$ is a singleton.

  2. 2.

There exist a $n^{*}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}$
with $n\geq n^{*}$, the sequence
$\left\\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\\}_{k\in\mathbb{N}}$ has the
following limit

$\lim_{k\to\infty}\;\hat{\bm{\theta}}_{(n,H)}^{(k)}=\hat{\bm{\theta}}_{(n,H)}.$

Moreover, there exists a real $\epsilon\in(0,\,1)$ such that for any
$k\in\mathbb{N}^{\ast}$

$\left\lVert\hat{\bm{\theta}}_{(n,H)}^{(k)}-\hat{\bm{\theta}}_{(n,H)}\right\rVert_{2}=\mathcal{O}_{\rm
p}({p}^{\nicefrac{{1}}{{2}}}\,\epsilon^{k}).$

  3. 3.

$\hat{\bm{\theta}}_{(n,H)}$ is a consistent estimator of $\bm{\theta}_{0}$, in
that for all $\varepsilon>0$ and all $\delta>0$, there exists a sample size
$n^{\ast}\in\mathbb{N}^{\ast}$ such that for all $n\in\mathbb{N}^{\ast}$
satisfying $n\geq n^{\ast}$ we have:

$\Pr\left(||\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}||_{2}\geq\varepsilon\right)\leq\delta.$

######  Theorem 2:

Under Assumptions A, B, C and D∗, and for all $H\in\mathbb{N}^{\ast}$, the IB-
estimator $\hat{\bm{\theta}}_{(n,H)}$ is consistent and PT-unibiased.
Therefore, the IB-estimator satisfies:

  1. 1.

$\hat{\bm{\theta}}_{(n,H)}$ is a consistent estimator of $\bm{\theta}_{0}$, in
that for all $\varepsilon>0$ and all $\delta>0$, there exists a sample size
$n^{\ast}\in\mathbb{N}^{\ast}$ such for all $n\in\mathbb{N}^{\ast}$ satisfying
$n\geq n^{\ast}$ we have:

$\Pr\left(||\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}||_{2}\geq\varepsilon\right)\leq\delta.$

  2. 2.

$\hat{\bm{\theta}}_{(n,H)}$ is PT-unbiased, that is there exists a
$n^{*}\in\mathbb{N}^{*}$ such that for all $n\in\mathbb{N}^{*}$ satisfying
$n\geq n^{*}$, we have
$\big{\lVert}\mathbb{E}[\hat{\bm{\theta}}]-\bm{\theta}_{0}\big{\rVert}_{2}=0$.

######  Theorem 3:

Under Assumptions A, B, C∗ and D∗∗, and for all $H\in\mathbb{N}^{\ast}$, the
IB-estimator $\hat{\bm{\theta}}_{(n,H)}$ satisfies:

  1. 1.

$\hat{\bm{\theta}}_{(n,H)}$ is a consistent estimator of $\bm{\theta}_{0}$, in
that for all $\varepsilon>0$ and all $\delta>0$, there exists a sample size
$n^{\ast}\in\mathbb{N}^{\ast}$ such for all $n\in\mathbb{N}^{\ast}$ satisfying
$n\geq n^{\ast}$ we have:

$\Pr\left(||\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}||_{2}\geq\varepsilon\right)\leq\delta.$

  2. 2.

$\hat{\bm{\theta}}_{(n,H)}$ is PT-unbiased, that is there exists a
$n^{*}\in\mathbb{N}^{*}$ such that for all $n\in\mathbb{N}^{*}$ satisfying
$n\geq n^{*}$, we have
$\big{\lVert}\mathbb{E}[\hat{\bm{\theta}}]-\bm{\theta}_{0}\big{\rVert}_{2}=0$.

  3. 3.

For all $\mathbf{u}\in{\rm I\\!R}^{p}$ such that $||\mathbf{u}||_{2}=1$, there
exist a sample size $n^{\ast}\in\mathbb{N}^{\ast}$ and a sequence
$\left(\delta_{n}\right)\in\mathfrak{D}$ such that for all for all
$n\in\mathbb{N}^{*}$ satisfying $n\geq n^{*}$ we have

$\scalebox{0.99}{$\sqrt{n}\left(1+\frac{1}{H}\right)^{-\nicefrac{{1}}{{2}}}\mathbf{u}^{T}\bm{\Sigma}(\bm{\theta}_{0},n)^{-\nicefrac{{1}}{{2}}}\mathbf{B}(\bm{\theta}_{0},n)\left(\hat{\bm{\theta}}_{(n,H)}-\bm{\theta}_{0}\right)\overset{d}{=}Z+\delta_{n}\mathcal{O}_{\rm
p}\left(\max\left(1,\frac{p^{2}}{\sqrt{H}}\right)\right)$},$

where $Z\sim\mathcal{N}(0,1)$.

Proof: The figure below delivers a proof of Theorems 1, 2 and 3.

Proposition 3 PT-Unbiased Corollary 1 Consistency Proposition 1 IB convergence
Proposition 2 Consistency Proposition 4 Asym. normality
D∗DD∗∗AABD3ABC1D4A1BCD1A2BCD2ABC∗D5 Assumptions used in Theorem 1 Assumptions
used in Theorem 2 Assumptions used in Theorem 3

Figure 4: Illustration of the implication links ($\longrightarrow$) of the
different assumptions used in Appendices D to G and the sufficient conditions
A, B, C, C∗, D, D∗ and D∗∗, needed for Theorems 1 to 3. The double arrows
($\Longrightarrow$) are used to denote implication of boxes (i.e. a set of
assumptions) while simple arrows are used to represent implications between
assumptions. The color boxes are used to highlight the assumptions used in
Theorems 1 to 3.

∎

######  Remark E:

Other results can be derived from Figure 4. For example, one can state a
theorem ensuring that the IB-estimator $\hat{\bm{\theta}}_{(n,H)}$ is unbiased
and asymptotically normal (and therefore consistent) without the guarantee
that the IB sequence
$\left\\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\\}_{k\in\mathbb{N}}$ is
convergent. This can be of interest if another method is used to compute the
IB-estimator. Another case is when one is only interested in the convergence
of the IB sequence
$\left\\{\hat{\bm{\theta}}_{(n,H)}^{(k)}\right\\}_{k\in\mathbb{N}}$ and the
consistency of the IB-estimator $\hat{\bm{\theta}}_{(n,H)}$. Assumption A is
stronger than necessary for such a result. The reason we use this assumption
to prove Theorem 1 is guided by the rationale that we want a restricted number
of general assumptions for the main Theorems 1, 2 and 3.