# Particle acceleration controlled by ambient density in the southwestern rim
of RCW 86

Hiromasa Suzuki Department of Physics, Faculty of Science and Engineering,
Konan University, 8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan
Satoru Katsuda  Graduate School of Science and Engineering, Saitama
University, 255 Shimo-Ohkubo, Sakura, Saitama 338-8570, Japan  Takaaki Tanaka
Department of Physics, Faculty of Science and Engineering, Konan University,
8-9-1 Okamoto, Higashinada, Kobe, Hyogo 658-8501, Japan  Nobuaki Sasaki
Graduate School of Science and Engineering, Saitama University, 255 Shimo-
Ohkubo, Sakura, Saitama 338-8570, Japan  Tsuyoshi Inoue Department of Physics,
Faculty of Science and Engineering, Konan University, 8-9-1 Okamoto,
Higashinada, Kobe, Hyogo 658-8501, Japan  Federico Fraschetti Center for
Astrophysics $|$ Harvard & Smithsonian, Cambridge, MA, 02138, USA  Department
of Planetary Sciences, Lunar & Planetary Laboratory, University of Arizona,
Tucson, AZ 85721, USA

###### Abstract

Particle acceleration physics at supernova remnant (SNR) shocks is one of the
most intriguing problems in astrophysics. SNR RCW 86 provides a suitable
environment for understanding the particle acceleration physics because one
can extract the information of both accelerated particles and acceleration
environment at the same regions through the bright X-ray emission. In this
work, we study X-ray proper motions and spectral properties of the
southwestern region of RCW 86. The proper motion velocities are found to be
$\sim 300$–2000 km s-1 at a distance of 2.8 kpc. We find two inward-moving
filaments, which are more likely reflected shocks rather than reverse shocks.
Based on the X-ray spectroscopy, we evaluate thermal parameters such as the
ambient density and temperature, and non-thermal parameters such as the power-
law flux and index. From the flux decrease in time of several non-thermal
filaments, we estimate the magnetic field amplitudes to be $\sim 30$–100
$\mu$G. Gathering the physical parameters, we then investigate parameter
correlations. We find that the synchrotron emission from thermal-dominated
filaments is correlated with the ambient density $n_{\rm e}$ as $\text{(power-
law flux)}\propto n_{\rm e}^{1.0\pm 0.2}$ and $\text{(power-law index)}\propto
n_{\rm e}^{0.38\pm 0.10}$, not or only weakly with the shock velocity and
shock obliquity. As an interpretation, we propose a shock-cloud interaction
scenario, where locally enhanced magnetic turbulence levels have a great
influence on local acceleration conditions.

acceleration of particles — shock waves — ISM: supernova remnants — X-rays:
ISM — ISM: individual objects (RCW 86)

††journal: ApJ††facilities: Chandra ††software: HEAsoft (v6.20; HEASARC 2014),
CIAO (v4.11; Fruscione et al. 2006), mkacispback (v2021-07-15; Suzuki et al.
2021)

##  1 Introduction

Particle acceleration physics at supernova remnant (SNR) shocks is one of the
principal problems in astrophysics as a promising mechanism to produce cosmic
rays below the “knee” energy ($\approx 3\times 10^{15}$ eV). X-ray and gamma-
ray studies of SNRs have revealed several aspects of the particle acceleration
physics in SNRs. Völk et al. (2005) and Vink (2006) suggested that magnetic
field amplification is very effective in young ($<2$ kyr) SNRs. High
turbulence levels of magnetic fields in some parts of young SNRs were found
with the highest level close to the Bohm limit at $\sim 2$ kyr (Tsuji et al.,
2021). Reynolds et al. (2021) studied the energy amount of accelerated
electrons and magnetic field in young SNRs based on several physical
parameters such as radio luminosity, plasma density, and shock velocity,
finding large variations in them among objects controlled by unknown factors.
Suzuki et al. (2022) also found that the maximum energies of accelerated
protons differ by more than one order of magnitude among objects at similar
ages.

As the essential part of the mechanism to accelerate particles up to the knee
energy, the enhancement of magnetic-field strength and turbulence has been
attracting particular interest (e.g., Bell 2004; Bamba et al. 2005; Vink et
al. 2006; Amato & Blasi 2006; Uchiyama et al. 2007). As a possible cause of
such a field enhancement, shock-cloud interactions are thought to be important
(Giacalone & Jokipii, 2007; Inoue et al., 2012; Fraschetti, 2013). The
magnetic turbulence is expected to be amplified around dense clumps and
enhanced synchrotron X-rays have been in fact observed (e.g., Sano et al.
2013, 2015, 2017).

SNR RCW 86 provides a suitable environment for understanding the particle
acceleration microphysics. One can extract the information of both accelerated
particles and acceleration environment at the same regions because the bright
X-ray emission exhibits both thermal and non-thermal components. RCW 86 is
believed to be the remnant of the oldest historical supernova of A.D. 185
(Stephenson & Green, 2002; Green & Stephenson, 2003). RCW 86 is located at
($l,b$) = (315.4, $-$2.5) and has a radio shell with a radius of $\sim
21^{\prime}$, which is almost completely surrounded by Balmer-dominated
filaments (Smith, 1997). The distance is estimated to be 2.8 kpc (Rosado et
al., 1996). RCW 86 is thought to have evolved in a low-density cavity region
and is currently interacting with dense materials (Williams et al., 2011). The
highly irregular morphology of the SNR shell indicates that RCW 86 is
currently expanding in a very inhomogeneous ambient medium. Such inhomogeneity
yields a broad range of shock velocities and magnetic turbulence around the
shocks that will affect the efficiency of particle acceleration. The north-
eastern (NE) corner is thought to be expanding in a rather low-density medium
with large velocities (Yamaguchi et al., 2016) and emits hard non-thermal
X-rays (Bamba et al., 2000; Vink et al., 2006). On the other hand, the radio-
brightest south-western (SW) corner is likely interacting with a dense cloud
that modifies the shock structure, thereby reducing or enhancing the non-
thermal X-ray emission (Rho et al., 2002; Sano et al., 2017, 2019).

In this work, we focus on the SW region. We investigate the shock velocities,
spectral features around shocks, and particle acceleration environments. The
observation logs and data reduction process are described in Section 2. Our
analysis procedure and results are presented in Section 3. The shock
structures and acceleration environments are discussed in Section 4.

##  2 Observations and Data Reduction

We use all the four existing Chandra observations of the RCW 86 SW region
listed in Table 1. The baseline for the proper motion study is $\approx 12$
year, which consists of the first-epoch observation in 2001 (OBSID 1993) and
the second-epoch ones in 2013 (OBSIDs 13748, 15610, and 15611). The total
exposure time is 177 ksec.

The RCW 86 SW region is observed with the Advanced CCD Imaging Spectrometer
(ACIS; Garmire 1997) S2, S3, and I3 in 2001 and with S2, S3, and I2 in 2013.
All the data were taken in the FAINT mode. We process the raw data following
the standard data reduction method (chandra_repro). We used CIAO (v4.11;
Fruscione et al. 2006) and calibration database 4.8.3 for the data reduction.

Table 1: Chandra observation logs of the RCW 86 SW region  OBSID | R.A. (2000.) | Decl. (2000.) | Roll angle | Date | Exposure (ksec) | PI
---|---|---|---|---|---|---
1993 | 220.19279 |  $-$62.66287 | 80$\fdg$2 | 2001 Feb 01 | 92 | S. Reynolds
13748 | 220.11392 |  $-$62.71971 | 70$\fdg$7 | 2013 Feb 14 | 36 | S. Katsuda
15610 | 220.11389 |  $-$62.71975 | 70$\fdg$7 | 2013 Feb 17 | 26 | S. Katsuda
15611 | 220.11392 |  $-$62.71971 | 70$\fdg$7 | 2013 Feb 12 | 23 | S. Katsuda

##  3 Analysis and Results

We perform proper motion study, spectroscopy, and filament-width measurement
for the RCW 86 SW region. The procedures and results are presented in this
section. In our analysis, we use HEAsoft (v6.20; HEASARC 2014), XSPEC
(v12.9.1; Arnaud 1996), and AtomDB 3.0.9. Throughout the paper, uncertainties
in the text, figures, and tables indicate $1\sigma$ confidence intervals.

The wide-band (0.5–7.0 keV) and hard band (2.0–7.0 keV) images extracted from
the observation in 2001 are presented in Figures 1 and 2, respectively, with
an indication of the analysis regions. Note that only SW5, SW6, SW7, SW9, and
SW10 are prominent in the hard X-ray image. We also make an image showing the
difference between the exposure-corrected 1.0–5.0 keV fluxes in 2001 and 2013
(Figure 3). For some filaments such as SW6, SW7, and SW9, their motions are
visible in this image. These 12 analysis regions are selected to enclose all
the bright and sharp filament structures seen in the field of view. The flux-
profile extraction directions are defined by eye to approximately match the
directions perpendicular to the filament structures.

Figure 1: Exposure-corrected 0.5–7.0 keV image obtained with Chandra in 2001.
The white and green lines indicate the analysis regions for proper motion and
spectral studies, respectively. The arrows beside the regions indicate the
moving directions and velocities (proportional to the arrow lengths) of the
filaments as described in Section 3.1.2. Their directions also correspond to
the positive directions of the positions in Figures 17 and 18 in the Appendix
A.  Figure 2: Exposure-corrected 2.0–7.0 keV image in 2001. The white
rectangular regions are used for the non-thermal spectral variation study. The
arrows beside the regions indicate the moving directions and velocities
(proportional to the arrow lengths) of the filaments as described in Section
3.1.2. The arrows also correspond to the positive directions of the angular
positions in Figure 9.  Figure 3: The image difference of the exposure-
corrected 1.0–5.0 keV images between 2001 and 2013. The image is smoothed with
a Gaussian kernel of $\sigma=0\farcs 8$. The white and black pixels show
increased and decreased fluxes from 2001 to 2013, respectively. The cyan
regions and arrows indicate the analysis regions and proper motion directions
and velocities of the filaments. Note that a triangular region at the bottom-
left corner is covered only by the observation in 2001.

###  3.1 Proper motion measurements

####  3.1.1 Aspect Correction

In order to perform the proper motion study with the best position accuracy
available, we apply the aspect correction to the three observations in 2013.
The six point-like sources indicated in Figure 4 are selected and used for the
aspect correction because they are bright and close to the nominal points in
all the observations. Before the aspect correction, the positions of these
sources differ among the observations typically by $\sim 0\farcs 3$. We first
run the wavdetect tool for each observation to determine the central positions
of these sources, which are summarized in Table 2. Then we run the wcs_match
tool to find the best-fit transformation matrices to correct the coordinate
systems of the three observations in 2013 to match those of the observation in
2001111All the six point sources are used in all the calculations without
being excluded in the wcs_match process.. The transformation matrix only uses
two-dimensional translation without rotation and scaling because the number of
available point-like sources is relatively small. Finally, the wcs_update tool
applies these corrections to the observations in 2013.

After the correction, for each of the three observations in 2013, the
resultant correction accuracy is evaluated based on the standard deviation of
the position offsets of the six sources with respect to those in 2001. The
correction accuracies are obtained as $0\farcs 18$, $0\farcs 24$, and $0\farcs
26$, for OBSIDs 13748, 15610, and 15611, respectively. We merge the three
observations in 2013 for use in the following analysis. The aspect correction
accuracy for these merged data is estimated by taking an exposure-weighted
average of those for the three observations, to be $0\farcs 22$, which is
converted to $0\farcs 018$ yr-1. This is considered as the systematic
uncertainty associated with the position accuracy.

Figure 4: Exposure-corrected 0.5–7.0 keV image in 2001 with the six point-like source locations used for aspect correction indicated with white circles.  Table 2: Positions of the six point-like sources  | Position: 1993 | Diff.: 1993–15610aaPosition differences between two observations in the 0.5–7.0 keV energy range. | Diff.: 1993–13748aaPosition differences between two observations in the 0.5–7.0 keV energy range. | Diff.: 1993–15611aaPosition differences between two observations in the 0.5–7.0 keV energy range.
---|---|---|---|---
Source | R.A. (2000.) | Decl. (2000.) | R.A. (′′) | Decl. (′′) | R.A. (′′) | Decl. (′′) | R.A. (′′) | Decl. (′′)
PS1 | 220.1181474 |  $-$62.7402549 |  $-$0.17 (0.33) |  $-$0.04 (0.17) |  $-$0.57 (0.29) | 0.35 (0.16) | 0.67 (0.43) | 0.71 (0.17)
PS2 | 220.0325486 |  $-$62.7092583 | 0.29 (0.32) |  $-$0.45 (0.10) | 0.29 (0.35) | 0.03 (0.10) | 0.09 (0.33) | 0.26 (0.10)
PS3 | 220.0828321 |  $-$62.7526891 |  $-$0.13 (0.26) |  $-$0.31 (0.17) |  $-$0.57 (0.22) | 0.22 (0.18) |  $-$0.39 (0.47) | 0.46 (0.12)
PS4 | 220.1688151 |  $-$62.7111788 |  $-$0.26 (0.21) |  $-$0.08 (0.12) |  $-$0.40 (0.21) | 0.45 (0.12) |  $-$0.20 (0.21) | 0.64 (0.18)
PS5 | 220.1801661 |  $-$62.7627948 | 0.64 (0.57) |  $-$0.71 (0.26) | 0.00 (0.36) |  $-$0.24 (0.17) | 0.13 (0.50) |  $-$0.61 (0.24)
PS6 | 220.1572650 |  $-$62.7548930 |  $-$1.15 (0.46) |  $-$0.44 (0.15) |  $-$0.84 (0.32) | 0.12 (0.16) |  $-$0.62 (0.36) |  $-$0.14 (0.21)

####  3.1.2 Proper motions

Proper motion velocities of the filaments indicated in Figure 1 are measured
as follows. Flux profiles are extracted from the vignetting-corrected images
in 2001 and 2013, which are presented in Figures 17 and 18. To evaluate their
proper motion velocities, a $\chi^{2}$-test is used as below. We artificially
shift the profile in 2013 by $\Delta x$ and calculate $\chi^{2}(\Delta x)$,
which is defined as

$\chi^{2}(\Delta x)=\sum_{i}\frac{(f_{i}-g(\Delta x)_{i})^{2}}{\Delta
f_{i}^{2}+\Delta g(\Delta x)_{i}^{2}},$ (1)

where $f_{i}$ and $\Delta f_{i}$ indicate the flux and error of the bin number
$i$ in 2001, and $g_{i}$ and $\Delta g_{i}$ indicate those of the shifted
profile in 2013. This calculation is repeated with various values of $\Delta
x$ to plot $\chi^{2}$ as a function of $\Delta x$. The minimum $\chi^{2}$
value ($\chi^{2}_{\rm min}$) and corresponding profile shift ($\Delta x_{\rm
min}$) are determined by fitting the $\chi^{2}$–$\Delta x$ plot with a
parabola function. An example of the $\chi^{2}$–$\Delta x$ plot and parabola
fitting is presented in Figure 5. The best-fit $\Delta x_{\rm min}$ is
converted to the proper motion velocity. The profile shift is not limited to
an integer multiple of the bin width. We re-bin the shifted profile $g(\Delta
x)$ with the same bin arrangement as $f$ with an assumption of a uniform
probability distribution inside each bin. Then, the profile-shift ranges which
give $\chi^{2}(\Delta x)=\chi^{2}_{\rm min}+1$ are calculated from the best-
fit parabola functions. These ranges are considered to be 1$\sigma$ confidence
ranges of the proper motion velocities.

Figure 5: Example of the $\chi^{2}$ test for proper motion measurement for the
SW1 region. The $\chi^{2}$ values (black solid line) and best-fit parabola
model (red line) are shown. The black dotted lines indicate the angular range
for the parabola fitting.  Figure 6: Proper motion estimates for individual
filaments in units of arcsec yr-1 (on the left) and km s-1 (on the right). A
distance of 2.8 kpc is assumed.

The resultant proper motion estimates are summarized in Table 3, and shown in
Figures 17 and 18. Note that the filaments SW7 and SW10 are moving inward,
toward the SNR center. These inward movements are confirmed by spectral
softening toward presumable downstream in SW7 and SW10 as described in Section
3.3. We also note that the projection effect, i.e., the difference between the
measured radial velocities and actual three-dimensional velocities, will be
small ($\lesssim 15\%$) considering the positions of the filaments with
respect to the apparent SNR radius.222We assume the SNR radius of $\approx
22.6^{\prime}$ and radial distances of the filaments from the explosion center
of $\approx 19.5\text{--}22.6^{\prime}$.

Table 3: Proper motion velocities  Name | Velocity (′′ yr-1)aaMinus velocities indicate movements toward the SNR center. | Velocity (km s-1)aaMinus velocities indicate movements toward the SNR center.
---|---|---
SW1 | $0.051\pm 0.020$ | $660\pm 250$
SW2 | $0.028\pm 0.005$ | $360\pm 70$
SW3 | $0.050\pm 0.015$ | $650\pm 190$
SW4 | $0.065\pm 0.007$ | $840\pm 90$
SW5 | $0.180\pm 0.026$ | $2330\pm 340$
SW6 | $0.072\pm 0.026$ | $930\pm 340$
SW7 | $-0.118\pm 0.029$ | $-1530\pm 370$
SW8 | $0.105\pm 0.021$ | $1350\pm 270$
SW9 | $0.165\pm 0.018$ | $2130\pm 230$
SW10 | $-0.166\pm 0.038$ | $-2140\pm 500$
SW11 | $0.093\pm 0.024$ | $1200\pm 310$
SW12 | $0.109\pm 0.025$ | $1410\pm 320$

###  3.2 Filament Widths and Their Energy Dependence

Figure 7: Flux profiles of the non-thermal filaments. The extraction energy
ranges are 0.5–2.0 keV (black) and 2.0–7.0 keV (red). Increasing positions
correspond to the directions of the arrows shown in Figure 1. The displayed
flux ranges are different for different panels. The best-fit models for the
two profiles are overplotted with black and red solid lines, respectively.
Radial ranges used for the fitting are indicated with black dotted lines.
Figure 8: Upstream width parameter $w_{\rm u}$ and filament widths in FWHM of
the non-thermal filaments in the 0.5–2.0 keV and 2.0–7.0 keV energy ranges.

In order to investigate what determines the maximum energies of accelerated
electrons seen in X-rays, the filament widths and their energy dependence are
measured. Here, we focus on the non-thermal-dominated filaments in 0.5–7.0
keV, SW5, SW6, SW7, and SW10. We merge the images in 2001 and 2013 and use the
merged one in this section. We extract the flux radial profiles from the two
energy ranges, 0.5–2.0 keV and 2.0–7.0 keV. To model the profiles, we use a
function

$\begin{cases}F_{\rm u}(x)&={C}\,\exp\left[-\frac{x-x_{0}}{w_{\rm
u}}\right]\quad{(x\geq x_{0})}\\\ F_{\rm
d}(x)&={C}\,\exp\left[\frac{x-x_{0}}{w_{\rm
d}}\right]\quad{(x<x_{0})},\end{cases}$ (2)

where the $F_{\rm u}(x)$ and $F_{\rm d}(x)$ are presumable upstream and
downstream fluxes as functions of angular position $x$, respectively. The
parameters $C$, $x_{0}$, $w_{\rm u}$ and $w_{\rm d}$ are normalization
parameter, angular position of the flux peak, and parameters to determine the
widths in both regions, respectively.

The flux profiles in the soft and hard energy bands are presented in Figure 7.
We fit the profile models to the data using the radial ranges where the
filaments are bright compared to the background emission. We show the upstream
width parameter $w_{\rm u}$ determined for the two energy ranges in the left
panel of Figure 8. Based on the derived parameters, we also calculate their
filament widths in Full Width Half Maximum (FWHM), which are presented in the
right panel of Figure 8. Most cases show significantly narrower widths in the
higher energies. Thus, their maximum energies will be determined by
synchrotron cooling, as in the case of Tycho’s SNR (Tran et al., 2015) and SN
1006 (Ressler et al., 2014), not by alternative processes such as the damping
of downstream magnetic field (e.g., Pohl et al. 2005), which predicts filament
widths independent of photon energy.

We confirm that the energy dependence of the point spread function (PSF) at
the positions of these filaments is $\lesssim 0\farcs 5$ by comparing the PSFs
at two representative energies, 1.0 and 4.0 keV (with the tool psfmap; Allen
et al. 2004). This small difference is because these filaments are located
within $\sim 5^{\prime}$ from the on-axis direction. Thus, the effect of the
energy dependent PSF on the filament widths will be negligible. Also, we check
the energy dependence of the widths of the thermal-dominated filaments to
examine possible systematics. We find that the widths in the 0.5–0.7 and
0.7–1.2 keV energy ranges of the thermal-dominated filaments are consistent
with each other. 333For SW1, as an example, the filament widths in the lower
and higher energy ranges are $42\pm 6^{\prime\prime}$ and $49\pm
3^{\prime\prime}$ in FWHM, respectively.

###  3.3 Spectral Variations in the Downstream Regions of the Non-thermal
Filaments

Here we examine the spatial variation of the spectral shapes downstream the
shocks for the non-thermal-dominated filaments. In this section, we focus on
the filaments which are non-thermal-dominated in the 1.0–7.0 keV band, SW5,
SW6, SW7, SW9, and SW10. 444We find that the non-thermal parameters are
constrained well if a spectrum is non-thermal-dominated above $\sim 1$ keV.
Spectral extraction regions are indicated in Figure 2. We merge the spectra in
2001 and 2013 to increase the statistics. Their proper motions are less than
30% of the widths of the extraction regions. For the spectral modeling here,
we use the model, Abs. (source emission) + (sky background) + PB. The source
emission is assumed to be powerlaw \+ vpshock for SW9 and powerlaw for the
others. For the vpshock model, only the normalization is treated as a free
parameter whereas the other parameters are fixed to the best-fit values
determined in Section 3.4.2.

Resultant power-law indices are presented in Figure 9. The SW7, SW9, and SW10
regions show spectral hardening in their shock downstream regions toward the
shock fronts. This is another evidence that the synchrotron emission dominates
the non-thermal component for these filaments (e.g., Katsuda et al. 2010;
Kishishita et al. 2013). We confirm that the spectra extracted from all the
regions are well explained with the powerlaw model and thus thermal
contamination will be negligible.

Figure 9: Variation of power-law indices in the downstream regions toward the
filaments for SW5, SW6, SW7, SW9, and SW10. Increasing positions correspond to
the directions of the arrows beside the regions in Figure 2. For each panel,
the cross with the largest position is the data at the filament.

###  3.4 Spectroscopy

In order to extract thermal and non-thermal properties from individual
filaments, we model here their energy spectra.

####  3.4.1 Background Estimation

Regarding the background (sky and particle-induced background), we estimate
the sky background from a region outside the SNR shell (shown in Figure 1),
and model the particle-induced background using the tool mkacispback (Suzuki
et al., 2021). We simultaneously model the spectra extracted from the
background region for both 2001 and 2013 observations. The spectral model for
the background region is FE + Abs. (MWH + CXB) + PB, where the FE, Abs., MWH,
CXB, and PB indicate the foreground emission (Yoshino et al., 2009; Kuntz &
Snowden, 2000), interstellar absorption, Milky Way halo (transabsorption
emission) (Masui et al., 2009; Yoshino et al., 2009; Kuntz & Snowden, 2000),
cosmic X-ray background (Snowden et al., 1992; Kushino et al., 2002; Hickox &
Markevitch, 2006), and particle-induced background model suited for our
background region.

The FE component is described by the apec model with the fixed temperature of
0.1 keV and metal abundances of solar values. The emission measure, which is
defined as (4$\pi$D2)-1$\int$nenHdV cm-5, where D, ne. and nH stand for
distance, electron and hydrogen number densities, respectively, is treated as
a free parameter. The absorption column density of the Abs. component (tbabs
model) is fixed to $N_{\rm H}=6.4\times 10^{21}$ cm-2 (HI4PI Collaboration et
al., 2016). The MWH component is also described with the apec model with the
fixed metal abundances of solar values and free temperature and emission
measure. The CXB component is described by the powerlaw model with the fixed
spectral index of 1.4 and the normalization corresponding to the flux of
$6.38\times 10^{-8}$ erg cm-2 s-1 str-1 (Kushino et al., 2002). As for the
particle-induced background, the acispback model predicts a lower flux than
the observation in the $\sim 2$–7 keV band only in 2001 (Suzuki et al., 2021).
Thus, we apply an additional powerlaw model to acispback only for the spectra
in 2001.

The spectral fitting results are shown in Figure 10. The best-fit spectral
parameters are summarized in Table 4. The spectral parameters are consistent
with Kuntz & Snowden (2000) and Yoshino et al. (2009). Note that the visible
flux decrease from 2001 to 2013 at the energies below 1 keV is due to the
increased contamination on the sensor surface.

Table 4: Best-fit spectral parameters for the background region

Model | Parameter | Value
---|---|---
Abs. |  $N_{\rm H}$ ($10^{22}$ cm-2) | 0.64 (fixed)
FE |  $kT$ (keV) | 0.1 (fixed)
| Abundance (solar) | 1 (fixed)
| EMaaEmission measure, 10-10(4$\pi$D2)-1$\int$nenHdV cm-5, where D, ne. and nH stand for distance (cm), electron and hydrogen number densities (cm-3), respectively. | $(9.7\pm 1.2)\times 10^{-1}$
MWH |  $kT$ (keV) | $0.18\pm 0.03$
| Abundance (solar) | 1 (fixed)
| EMaaEmission measure, 10-10(4$\pi$D2)-1$\int$nenHdV cm-5, where D, ne. and nH stand for distance (cm), electron and hydrogen number densities (cm-3), respectively. | $4.5\pm 0.3$
CXB | Photon index | 1.4 (fixed)
| NormalizationbbEmission measure in the same units as that in Table 4 |  $4.32\times 10^{-6}$ (fixed)

bbEmission measure in the same units as that in Table 4

Normalization of the power-law model in units of cm-2 s-1 keV-1 at 1 keV.

Figure 10: Energy spectra and best-fit models for the background region. The
black and red crosses represent the data in 2001 and 2013, respectively. The
model spectra only for the data in 2001 are overplotted with lines. The black
solid line represents the entire model. The black dotted, dashed, and dash-
dotted lines show the FE, MWH, and CXB spectral models, respectively. The
orange dotted and solid lines represent the additional powerlaw model to the
PB component and the entire PB spectral model, respectively.

####  3.4.2 Spectral Modeling for Filaments

Figure 11: Spectral modeling results for individual filaments. Black and red
crosses are the observations in 2001 and 2013, respectively. In each panel,
the upper and lower panels show the flux and data-to-model ratio,
respectively. The best-fit model spectra only for the data in 2001 are
overplotted with lines. The black solid lines represent the entire models. The
green, blue, and magenta solid lines represent the thermal (vpshock), non-
thermal (powerlaw), and Gaussian components, respectively. The orange dotted
and solid lines show the additional powerlaw models to the PB component and
the entire PB models, respectively. The black dotted lines show the sky
background components.

As the source emission, we apply the vpshock (non-equilibrium ionizing plasma)
model to the thermal, and the powerlaw model to the non-thermal components
following Tsubone et al. (2017). In addition, a Gaussian with the fixed energy
centroid of 6.4 keV and width of zero is added for all the regions (Yamaguchi
et al., 2011; Tsubone et al., 2017). The spectral model is described as Abs.
(source emission) + (sky background) + PB. The source emission is powerlaw for
SW5, SW6, SW7, and SW10 because thermal components are found to be
unnecessary. The model vpshock \+ powerlaw is applied for the other regions.
We note that the regions SW9 and SW12 are intermediate cases between the non-
thermal- and thermal-dominated ones. The absorption column density of the Abs.
component (tbabs model) is treated as a free parameter. For vpshock, the
temperature, ionization timescale ($n_{\rm e}t$), emission measure, and the
metal abundances of O, Ne, Mg, Si, and S are treated as free parameters. The
other metal abundances are fixed to solar values. All the parameters of
vpshock are tied between 2001 and 2013. The powerlaw normalization and index
are treated as free parameters and are basically tied between 2001 and 2013.
Only for SW5, SW6, SW7, and SW10, both parameters are untied between 2001 and
2013. The sky background component is the same as that determined in Section
3.4.1. Their overall normalizations are scaled by the area ratio between the
source regions and background region. The PB component for the spectra in 2001
is composed of acipback plus powerlaw models, whereas that for 2013 is a
simple acispback model, which is the same treatment as that for the background
estimation.

The spectra and best-fit models are presented in Figure 11. The resultant
spectral parameters are summarized in Table 5. The power-law fluxes of SW5,
SW6, SW7, and SW10 all show significant decrease from 2001 to 2013. Similarly,
the power-law indices increase (spectra soften) from 2001 to 2013. The
absorption column densities, temperatures, metal abundances, and ionization
timescales are roughly consistent with those presented in Tsubone et al.
(2017). From the plasma densities inferred from the emission measures ($\sim
1$ cm-3) and ionization timescales ($5\text{--}13\times 10^{10}$ s cm-3), we
confirm555With an assumption that our analysis regions are cuboids. that the
elapsed times after the shock heating are comparable to the remnant age ($\sim
2000$ yr). The SW3 region may show smaller elapsed time of $\sim 600$ yr. This
is consistent with the fact that this region is more distant from the remnant
center than the other filaments and will be less affected by the projection
effect, which mixes the newly heated and older regions. We also search for
possible variations of the non-thermal parameters in the thermal-dominated
cases where the spectral model fits the data relatively well (e.g., SW9 and
SW12), and find consistent parameters in 2001 and 2013. 666For example, the
power-law fluxes of SW9 in the 0.5–7.0 keV band in 2001 and 2013 are
$3.6\text{--}4.5\times 10^{-13}$ erg and $4.4\text{--}5.5\times 10^{-13}$ erg,
respectively.

Table 5: Best-fit spectral-model parameters for the filament regions  Name | Abs. | powerlaw | vpshock | Gaussian | C-stat (d.o.f.)
---|---|---|---|---|---
|  $N_{\rm H}$ (1022 cm-2) |  $\log_{10}$(Flux (erg cm-2 s-1)) aaUnabsorbed flux in the 0.5–7.0 keV energy range. For those with two values, the values in 2001 (2013) are shown. | Photon index |  $kT$ (keV) |  $n_{\rm e}t$ (1010 cm-3 s-1) | EM bbEmission measure in the same units as that in Table 4 | O (solar) | Ne (solar) | Mg (solar) | Si (solar) |  ${\rm Norm}$ccTotal flux of the Gaussian model in units of 10-7 photons cm-2 s-1. The energy centroid of the Gaussian component is fixed to 6.4 keV. |
SW1 | $0.23\pm 0.02$ | $-12.21\pm 0.02$ | $3.14\pm 0.12$ |  | $4.9\pm 0.9$ | $6.9\pm 1.0$ | $0.74\pm 0.05$ | $1.43\pm 0.09$ | $1.31\pm 0.09$ | $1.26\pm 0.14$ | $9.3\pm 4.5$ | 1603.64 (1288)
SW2 | $0.34\pm 0.02$ | $-11.89\pm 0.02$ | $4.31\pm 0.12$ |  | $12.5\pm 1.5$ | $25.8\pm 2.9$ | $1.02\pm 0.04$ | $2.11\pm 0.07$ | $1.44\pm 0.06$ | $1.62\pm 0.09$ | $2.4\pm 2.0$ | 1798.70 (1288)
SW3 | $0.44\pm 0.02$ | $\cdots$ | $\cdots$ |  | $0.97\pm 0.05$ | $4.2\pm 0.9$ | $0.78\pm 0.07$ | $1.23\pm 0.10$ | $0.79\pm 0.09$ | $1.00\pm 0.22$ | $<3.5$ | 1373.35 (1290)
SW4 | $0.43\pm 0.02$ | $-11.89\pm 0.02$ | $4.18\pm 0.12$ |  | $7.7\pm 1.0$ | $25.2\pm 4.4$ | $1.59\pm 0.08$ | $2.49\pm 0.12$ | $1.45\pm 0.08$ | $2.28\pm 0.16$ | $3.3\pm 2.3$ | 1468.98 (1287)
SW5 | $0.26\pm 0.01$ |  $-12.32\pm 0.01\,(-12.34\pm 0.01$) |  $2.59\pm 0.07\,(2.81\pm 0.06$) |  | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $1.2\pm 1.5$ | 1405.35 (1293)
SW6 | $0.33\pm 0.01$ |  $-12.04\pm 0.01\,(-12.10\pm 0.01$) |  $2.46\pm 0.05\,(2.66\pm 0.05$) |  | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $2.4\pm 2.1$ | 1315.36 (1293)
SW7 | $0.32\pm 0.02$ |  $-12.28\pm 0.01\,(-12.38\pm 0.01$) |  $2.82\pm 0.08\,(3.23\pm 0.08$) |  | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $<1.8$ | 1231.83 (1293)
SW8 | $0.27\pm 0.01$ | $-12.00\pm 0.02$ | $4.10\pm 0.11$ |  | $7.7\pm 0.09$ | $9.15\pm 0.07$ | $1.15\pm 0.06$ | $2.35\pm 0.09$ | $1.54\pm 0.07$ | $1.51\pm 0.11$ | $2.7\pm 1.8$ | 1799.19 (1287)
SW9 | $0.40\pm 0.06$ | $-12.34\pm 0.01$ | $2.66\pm 0.08$ |  | $9.5\pm 3.6$ | $2.6\pm 1.1$ | $1.52\pm 0.18$ | $2.67\pm 0.41$ | $1.85\pm 0.36$ | $1.90\pm 0.86$ | $1.4\pm 1.2$ | 1211.74 (1288)
SW10 | $0.25\pm 0.01$ |  $-11.93\pm 0.01\,(-11.97\pm 0.01$) |  $2.57\pm 0.05\,(2.69\pm 0.04$) |  | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $4.9\pm 2.5$ | 1468.25 (1293)
SW11 | $0.33\pm 0.03$ | $-12.52\pm 0.02$ | $2.98\pm 0.12$ |  | $7.3\pm 1.7$ | $3.0\pm 0.6$ | $1.00\pm 0.08$ | $2.16\pm 0.16$ | $1.36\pm 0.13$ | $1.95\pm 0.25$ | $3.9\pm 1.7$ | 1336.28 (1287)
SW12 | $0.28\pm 0.05$ | $-12.49\pm 0.02$ | $2.57\pm 0.09$ |  | $4.5\pm 1.2$ | $1.7\pm 0.7$ | $1.01\pm 0.09$ | $2.49\pm 0.39$ | $1.43\pm 0.26$ | $2.18\pm 0.66$ | $<2.0$ | 1316.22 (1288)

We also use the spectral model by Zirakashvili & Aharonian (2007) (hereafter
the ZA07 model) for the non-thermal component to discuss the acceleration
efficiency. The ZA07 model is described as

$\frac{dN}{dE}\propto\left(\frac{E}{E_{0}}\right)^{-2}\left[1+0.38\left(\frac{E}{E_{0}}\right)^{1/2}\right]^{11/4}\exp\left[-\left(\frac{E}{E_{0}}\right)^{1/2}\right],$
(3)

where $E$ and $E_{0}$ indicate photon energy and spectral turnover energy,
respectively. This model describes the synchrotron emission spectrum with
several assumptions: loss-limited maximum energy, arbitrary energy dependence
of the diffusion coefficient, shock compression ratio of four, and upstream-
to-downstream magnetic field ratio of $\sqrt{11}$. We replace the powerlaw
model in the spectral model set used above with the ZA07 model. The free
parameters for spectral modeling are the same as those in the case with the
powerlaw model except for the normalization and turnover $E_{0}$ of the ZA07
model. The results are summarized in Table 6. One can see that the non-
thermal-dominated filaments, SW5, SW6, SW7, and SW10, exhibit decreasing
$E_{0}$ in time, which is the same tendency as our results with the powerlaw
model.

Table 6: Best-fit parameters of the ZA07 model  Name |  $E_{0}$ (keV)aaValues are tied between the data in 2001 and 2013 except for those for SW5, SW6, SW7, and SW10, for which the values in 2001 (2013) are shown.
---|---
SW1 | $0.18\pm 0.03$
SW2 | $0.041\pm 0.004$
SW3 | $\cdots$
SW4 | $0.051\pm 0.005$
SW5 | $0.42\pm 0.06\,(0.24\pm 0.02)$
SW6 | $2.35\pm 0.76\,(1.92\pm 0.64)$
SW7 | $0.26\pm 0.03\,(0.13\pm 0.01)$
SW8 | $0.08\pm 0.01$
SW9 | $1.72\pm 0.72$
SW10 | $0.45\pm 0.04\,(0.33\pm 0.02)$
SW11 | $0.24\pm 0.05$
SW12 | $2.52\pm 1.69$

##  4 Discussion

###  4.1 Magnetic field amplitude

We have found that the non-thermal-dominated filaments SW5, SW6, SW7, and SW10
showed flux decrease from 2001 to 2013. This is seen in Figure 12. The
decrease rates are from $\sim 0.4\%\text{ yr}^{-1}$ (SW5) to $\sim 2\%\text{
yr}^{-1}$ (SW7). Given that the non-thermal emissions are dominated by the
synchrotron radiation, it is reasonable to attribute such flux decreases to
the synchrotron cooling. The cooling timescale is described as

$t_{\rm syn}=50\,{\rm yr}\,\left(\frac{E}{1~{}{\rm
keV}}\right)^{-0.5}\left(\frac{B}{100~{}\mu{\rm G}}\right)^{-1.5},$ (4)

where $E$ and $B$ represent the synchrotron photon energy and magnetic field
strength, respectively. Thus, the flux decreases roughly require the field
strengths of $\sim 30$–100 $\mu$G. Such field strengths are similar to the
previous estimates (24 $\mu$G by Vink et al. 2006 based on filament widths; 35
$\mu$G by Helder et al. 2012 based on filament widths; 14–20 $\mu$G by Yuan et
al. 2014 based on the broadband spectrum).

Figure 12: Power-law fluxes of the individual filaments (top) and their
2013-to-2001 ratios (bottom; only for the non-thermal dominated filaments).
The flux is calculated for the 0.5–7.0 keV energy range.

The magnetic field strength at a filament can also be estimated based on the
filament width seen in X-rays. As presented in Vink et al. (2006), the
magnetic field strength is simply determined by the filament width $w$ at the
target distance $d$ as

$\displaystyle B$ $\displaystyle=54~{}\mu{\rm G}\,(w/0.1~{}{\rm pc})^{-2/3}$
(5) $\displaystyle\approx 54~{}\mu{\rm
G}\,(w/7.3^{\prime\prime})^{-2/3}(d/2.8~{}{\rm kpc})^{-2/3},$

if the energy range for the width estimation is near the roll-off energy of
the synchrotron emission. We thus obtain magnetic field strengths of $\sim
20$–50 $\mu$G for SW5, SW6, SW7, and SW10 based on Figure 8. These estimates
are also consistent with the previous estimates (Vink et al., 2006; Helder et
al., 2012; Yuan et al., 2014).

###  4.2 Magnetic-field turbulence level

Based on the shock velocity ($v_{\rm sh}$) and turnover energy of the ZA07
model ($E_{0}$), one can estimate the gyrofactor ($\eta=(B/\delta B)^{2}$,
where $B$ and $\delta B$ are the magnetic field amplitude and fluctuation,
respectively). Zirakashvili & Aharonian (2007) derived $E_{0}$ as a function
of $\eta$ and $v_{\rm sh}$ as

$E_{0}=0.24~{}{\rm keV}\,\eta^{-1}(v_{\rm sh}/2000~{}{\rm km~{}s^{-1}})^{2}.$
(6)

Figure 13 exhibits the plots of $E_{0}$ against the shock velocity, for which
the proper motion velocity is simply substituted. Compared to the
“$\eta=\text{const.}$” lines, SW3, SW4, SW5, SW7, SW8, SW10, SW11, and SW12
are consistent with $\eta=1$–4, whereas SW1, SW2, SW6, and SW9 have $\eta<1$.
Such small values of $\eta$ indicate the existence of highly amplified
magnetic turbulence. This is similar to the situations of SN 1006 (Ressler et
al., 2014) and RX J1713.7$-$3946 (Tsuji et al., 2019). Dickel et al. (2001)
found a low degree of polarization in the RCW 86 SW region and proposed a
Faraday depolarization scenario. The high level of magnetic turbulence
suggested from our results can be an alternative for this low degree of
polarization.

Given that RCW 86 is thought to have begun to interact with the cavity wall
recently (Williams et al., 2011), such small $\eta$ values might indicate that
the shocks have been decelerated just recently, whereas the accelerated
particles remain as they were before the interaction. This would be reasonable
if the synchrotron cooling timescales we evaluated ($\gtrsim 50$–200 yr) are
longer than the deceleration timescales of the shocks. We basically believe
that the projection effect on the proper motions will be small ($\lesssim
15\%$; see Section 3.1.2), but it still possibly contributes to such small
$\eta$ values because of the complicated shock structure. It is also possible
that these small $\eta$ values are due to the uncertainties in the distance to
RCW 86. A factor of two larger distance would result in $\eta\gtrsim 1$ for
most of the filaments. Note that larger distances would, on the other hand,
make the magnetic-field estimates based on the filament widths be in tension
with those from the flux decrease (Section 4.1).

Figure 13: Plots of the cutoff energy of the ZA07 model ($E_{0}$) against the
shock velocity (proper motion velocity is substituted). A distance of 2.8 kpc
is assumed. The black solid, dashed, and dotted lines indicate the ZA07 model
with different gyrofactors ($\eta$).

###  4.3 Nature of the inward-moving filaments

We have found that the SW7 and SW10 filaments are moving inward. The spectral
softening toward the downstream regions is also confirmed (Section 3.3), which
supports their inward movements. Here we discuss their nature: reverse or
reflected shocks.

Assuming that they are reverse shocks, we can estimate their shock velocities
in the unshocked-ejecta frame with an assumption of the explosion center
coordinates. The explosion center is substituted with the geometrical center,
and is assumed to be $(l,b)=(315\fdg 4068,-2\fdg 3151)$ by fitting the radio
shell with a circle by eye. The ejecta velocities in the laboratory frame
($v_{\rm ej}$) and the estimated shock velocities ($v_{\rm rev}$) are
summarized in Table 7. These velocities would be too large considering their
values of $E_{0}$, i.e., they would require much larger values of the
gyrofactor $\eta$ compared to the other filaments (Figure 13). Given that SW7
and SW10 are associated with rather hard X-ray emission, it would be
unreasonable if they were outliers in Figure 13 with particularly large $\eta$
values.777Note that the large $\eta$ values of $\sim 20$ suggested in the
reverse-shock scenario for SW7 and SW10 are not unreasonable themselves,
because similar values of $\eta$ have been found in other SNRs with similar
ages (Tsuji et al., 2021). Besides, the presumable reverse-shock velocities of
$\sim 10000$ km s-1 are even larger than the shock velocities inferred for the
inward-moving filaments observed in Cassiopeia A (e.g., Sato et al. 2018).

On the other hand, if we assume that the inward-moving filaments are reflected
shocks, we cannot simply estimate their shock velocities without certain
assumptions (e.g., Truelove & McKee 1999). In this case, the reflected shocks
should be moving in the shocked ejecta. If the ejecta would have been
decelerated significantly (e.g., by the cavity wall as suggested by Williams
et al. 2011 or dense clouds as suggested by Sano et al. 2019), the shock
velocities with respect to the shocked ejecta would be as small as those of
the other filaments (Figure 13). Thus, we propose the reflected-shock scenario
as the most likely nature of the two inward-moving filaments.

Table 7: Shock velocities in the reverse shock scenario Name | Radius (′)aaFilament position with respect to the geometrical center $(l,b)=(315\fdg 4068,-2\fdg 3151)$ |  $v_{\rm ej}$ (km s-1)bbFast-moving ejecta velocity in the laboratory frame |  $v_{\rm rev}$ (km s-1)ccShock velocity in the ejecta frame
---|---|---|---
SW5 | 22.8826 | 9668 | $7300\pm 700$
SW6 | 17.6610 | 7462 | $6500\pm 700$
SW7 | 18.4815 | 7809 | $9300\pm 900$
SW10 | 18.1032 | 7649 | $9800\pm 1000$

###  4.4 Maximum energies of accelerated protons

As we have obtained several parameters related to particle acceleration, the
age $t$, shock velocity $v_{\rm sh}$, magnetic field strength $B$, and
gyrofactor $\eta$, we derive here the maximum energies of particles. Those of
protons are of particular interest. With an assumption that the proton maximum
energy ($E_{\rm max,p}$) is limited by the acceleration time and that the
diffusion coefficient is time-invariant, $E_{\rm max,p}$ can be written as
(e.g., Reynolds 2008; Yamazaki et al. 2014)

$E_{\rm max,p}\approx 200~{}{\rm TeV}\,\eta^{-1}\left(\frac{v_{\rm
sh}}{2000\,\text{km
s${}^{-1}$}}\right)^{2}\left(\frac{B}{50\,\text{$\mu$G}}\right)\left(\frac{t}{2\,\text{kyr}}\right).$
(7)

If we substitute the parameters for SW7, $\eta\approx 1$, $v_{\rm sh}\approx
1600$ km s-1, and $B\sim 100~{}\mu$G (based on the non-thermal flux decrease),
we obtain $E_{\rm max,p}\sim 210$ TeV. Similarly, for SW5, SW6, and SW10, we
obtain $E_{\rm max,p}\sim 140$, 130, and 240 TeV, respectively. These
estimates are consistent with the maximum energies evaluated from the gamma-
ray spectrum assuming hadronic gamma-rays, $\approx 10\text{--}200$ TeV (Yuan
et al., 2014; Zeng et al., 2019; Suzuki et al., 2020, 2022)888We note that our
estimates are also consistent with the leptonic scenario, where the proton
maximum energy of $\gtrsim 20$ TeV is indicated.. Thus, our parameter
estimates such as the magnetic field strength and turbulence level would be
reasonable.

###  4.5 Acceleration parameters and their dependence on environment

Figure 14: Plots of the power-law flux (PL flux) and index (PL index) as
functions of the shock velocity (proper motion velocity is substituted),
emission measure ($\propto\text{density }(n_{\rm e})^{2}$), and electron
temperature ($kT$). The PL flux is calculated for 0.5–7.0 keV. The EM is shown
in the same units as Table 4. The 95% upper limits of the EMs of the non-
thermal-dominated filaments are shown with the blue arrows in panels (b). The
best-fit power-law functions for the plots of the PL flux and index over the
emission measure are also shown with red lines.  Figure 15: Plots of the
emission measure ($\propto n_{\rm e}^{2}$) as a function of the shock velocity
($v_{\rm sh}$; proper motion velocity is substituted). A distance of 2.8 kpc
is assumed. The EM is shown in the same units as Table 4. The solid lines
indicate two conditions, $n_{\rm e}v_{\rm sh}=\text{const.}$ (green) and
$n_{\rm e}v_{\rm sh}^{2}=\text{const.}$ (blue).

####  4.5.1 Parameter correlations

Based on the proper motions and spectral parameters we have obtained, we
discuss here their parameter correlations to see how the environmental
parameters affect the acceleration processes. Figure 14 shows the power-law
parameters versus shock velocity, thermal emission measure, and electron
temperature. Focusing on the filaments with both thermal and non-thermal
parameters obtained (black crosses), the power-law parameters clearly depend
on the emission measure ($\propto n_{\rm e}^{2}V$, where $n_{\rm e}$ and $V$
are plasma density and volume, respectively) and do not or only weakly on the
shock velocity and post-shock electron temperature. If we evaluate the
correlation between the downstream plasma density and power-law parameters
with a power-law function, we obtain $\text{(power-law flux)}\propto n_{\rm
e}^{1.0\pm 0.2}$ and $\text{(power-law index)}\propto n_{\rm e}^{0.38\pm
0.10}$. Note that the non-thermal-dominated filaments, SW5, SW6, SW7, and SW10
are not included in these evaluations. Such correlations were indeed implied
by Tsubone et al. (2017) using Suzaku by investigating the whole remnant, but
are more clarified in this work with Chandra. An increase in the non-thermal
flux associated with high-density regions was also observed in Cassiopeia A
(Sato et al. 2018; Fraschetti et al. 2018)

As seen in Figure 15, the emission measure is negatively correlated with the
shock velocity $v_{\rm sh}$. This correlation appears to follow the function
$n_{\rm e}v_{\rm sh}=\text{const.}$ A simple assumption of the constant ram
pressure $n_{\rm e}v_{\rm sh}^{2}=\text{const.}$ seems to be inapplicable.
This may reflect the inhomogeneity in the ambient density suggested by the
complicated structure of the remnant and discussed in previous works (Williams
et al., 2011; Tsubone et al., 2017). We note that additional uncertainties of
the shock velocity are expected in the moving direction and projection effect,
which would be less than a factor of $\sim$ two.

Certain correlations among magnetic field strength $B$, density $n_{\rm e}$,
and shock velocity $v_{\rm sh}$ are expected due to the Bell instability
(Bell, 2004), and were in fact confirmed from observations of several young
SNRs (Völk et al., 2005; Vink, 2006; Helder et al., 2012; Vink, 2017):
$B^{2}\propto n_{\rm e}v_{\rm sh}^{3}$ or $B^{2}\propto n_{\rm e}v_{\rm
sh}^{2}$. If we plot the filament width ($\propto B^{-3/2}$ by Eq. 5) over the
shock velocity for SW5, SW6, SW7, and SW10, we obtain Figure 16. As can be
seen, the data do not follow the relations expected above.

Figure 16: Plots of the filament width ($w$) against the shock velocity
($v_{\rm sh}$; proper motion velocity is substituted). A distance of 2.8 kpc
is assumed. The solid lines indicate two functions, $w\propto v_{\rm
sh}^{-3/2}$ (green) and $w\propto v_{\rm sh}^{-9/4}$ (blue).

In order to compare the parameter correlations of the non-thermal-dominated
filaments to those of the thermal-dominated ones, we evaluate the upper limits
of the EMs for the non-thermal-dominated filaments (SW5, SW6, SW7, SW9, and
SW10). We repeat the spectral modeling described in Section 3.4.2 with an
additional ionizing plasma model with various electron temperatures (0.3–0.8
keV) and ionization timescales ($10^{10\text{--}11}$ s cm-3) with the metal
abundances fixed to solar, and derive the upper limit of the EM. 999The two
parameter ranges are selected based on the values determined for the thermal-
dominated filaments (Table 5). The resulting upper limits of the EMs versus
the power-law parameters are shown in Figure 14 (b). The non-thermal-dominated
filaments located at the inner regions, SW6, SW7, and SW10, do not follow the
same tendency as that of the other filaments (shown with the red solid lines).

We have also investigated the dependence of the non-thermal parameters on the
shock obliquity (angle between the shock normal and background magnetic
field). The background field direction is estimated from the starlight
polarization at similar distances to RCW 86 (see Appendix B for more details).
The field direction around RCW 86 is found to be nearly parallel to the
Galactic plane. We find no correlation between the non-thermal parameters and
the shock obliquity as can be seen in Figure 20.

####  4.5.2 Scenario to explain the measurements: shock-cloud interaction

We propose here a shock-cloud interaction scenario to explain the parameter
correlations we have obtained above. As the RCW 86 SW region is interacting
with dense atomic and molecular clouds (Sano et al., 2017, 2019), X-ray
emission is expected to trace such an interaction. The shock-cloud interaction
will slow down the shocks and will damp the magnetic turbulence inside dense
clumps but amplify the magnetic turbulence around them (Inoue et al., 2012;
Fraschetti, 2013). If our analysis has resolved such clumpy structures, we
expect harder synchrotron X-rays for lower-density regions due to higher
maximum acceleration energies, which matches the measured trend (Figure 14
(b)). We note that our spectral extraction regions have a size of $\sim
0.2$–0.4 pc at a distance of 2.8 kpc, and the shock crossing time is estimated
as $\sim 100$–1000 yr. Since the cooling timescale derived for the non-
thermal-dominated filaments of $\sim 50$–200 yr is similar to the shock
crossing timescale, it will be reasonable to assume that the spatial extent of
clumps are similar to our region sizes. The fact that the power-law flux is
proportional to the density (Figure 14 (b-1)) can be understood as more
enhanced non-thermal particle density in higher-density regions. The shock-
cloud interaction leads highly spatially inhomogeneous magnetic fields, and so
the power-law parameters show large scatters even at similar shock velocities
(Figure 14 (a)). However, interestingly, the plasma density seems to remain as
a good tracer of the modification of the acceleration conditions even after
the interaction (Figure 14 (b)).

The reason for the apparent inconsistency of our estimates of the magnetic
field strengths and shock velocities with the Bell instability (Figure 16) is
understandable if the magnetic turbulence induced by the shock-cloud
interaction determines local acceleration activity. Note that, according to a
comparison with other SNRs, the base magnetic field amplification level over
the whole remnant seems to be controlled by the Bell instability (Vink, 2006;
Helder et al., 2012). Thus, it is suggested that the local magnetic turbulence
levels $\delta B/B$ are of great importance to understand the local
acceleration conditions.

In Figure 14 (b), the tendency of the inner non-thermal-dominated filaments
(SW6, SW7, and SW10) differs from that of the outer filaments. If we assume
that the outer filaments are newly interacting with dense gas and the inner
ones are still in the wind-blown bubble, our results may suggest that the
cavity region has different acceleration conditions from those of the
interacting regions, which is naturally expected.

As a conclusion, our findings suggest that the acceleration physics at the
outer filaments of the RCW 86 SW region is governed by the ambient density,
not or only weakly by the shock velocity and shock obliquity. We find that the
shock-cloud interaction scenario can explain the measurements consistently,
although not yet in a quantitative manner. Radio observations with high
angular resolutions (e.g., Atacama Large Millimeter/submillimeter Array: ALMA)
will greatly help test our scenario.

##  5 Conclusion

In this work, we studied the X-ray proper motions and spectral properties of
the RCW 86 SW region. The proper motion velocities were found to be $\sim
300$–2000 km s-1 at a distance of 2.8 kpc. We found two inward-moving
filaments. They were found to be non-thermal dominated and the spectral
softening toward downstream were seen, which confirmed their inward movements.
It is likely that they are reflected shocks rather than reverse shocks. Based
on the X-ray spectroscopy, we evaluated thermal parameters such as the ambient
density and temperature, and non-thermal parameters such as the power-law flux
and index. Also, based on the flux decrease of several non-thermal filaments,
we were able to estimate the magnetic field amplitudes of $\sim 30$–100
$\mu$G.

Gathering the proper motion and X-ray properties, we then studied the
parameter correlations. We found that, at the outer thermal-dominated
filaments, the non-thermal parameters were correlated with the ambient density
as $\text{(power-law flux)}\propto n_{\rm e}^{1.0\pm 0.2}$ and $\text{(power-
law index)}\propto n_{\rm e}^{0.38\pm 0.10}$, not or only weakly with the
shock velocity and shock obliquity. These indicate harder and fainter
synchrotron emission for lower-density regions. As an interpretation of the
measured physical parameters, we propose the shock-cloud interaction scenario,
where the locally enhanced magnetic turbulence levels ($\delta B/B$) have a
great influence on the local acceleration conditions. The inner non-thermal-
dominated filaments showed a different tendency from that of the outer
filaments, which is understandable if the inner ones are still in the wind-
blown bubble and have different acceleration conditions.

We appreciate helpful suggestions provided by the anonymous referee, which has
improved the paper significantly. We are grateful to K. Kawabata for providing
the software to extract the starlight polarization data. This work was
partially supported by JSPS grant Nos. 19J11069 and 21J00031 (HS), 20H00174,
21H01121 (SK), 19H01936, 21H04493 (TT), and 20H01944 (TI). FF was supported,
in part, by NASA through Chandra Theory Award Number $TM0-21001X$, issued by
the Chandra X-ray Observatory Center, which is operated by the Smithsonian
Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060.

##  Appendix A Flux profiles across the filaments

The X-ray flux profiles of individual filaments in 2001 and 2013 are shown in
Figures 17 and 18. The radial ranges indicated with the dotted lines are used
to calculate the proper motion velocities.

Figure 17: Flux profiles of the regions SW1–SW6. The black and red crosses
represent the data taken in 2001and 2013, respectively. The positive
directions of positions correspond to the directions of the arrows shown in
Figure 1. The displayed flux ranges are different for different panels. The
vertical dashed lines represent the ranges used for proper motion measurement.
Figure 18: Flux profiles of the regions SW7–SW12. Same convention for the
plots is used as in Figure 17.

##  Appendix B Other parameter correlations

In order to investigate the dependence of the non-thermal parameters on the
shock obliquity, we first estimate the magnetic field directions as follows.
We use the starlight polarization data compiled by Heiles (2000). This
database includes the polarization properties of $\sim 10000$ stars with their
positions and distances. We extract the data of the stars in a 20∘ $\times$
20∘ square region centered on $(l,b)=(315\fdg 4015,-2\fdg 31664)$ at distances
of 2–4 kpc. The resultant magnetic field directions around RCW 86, which
directly correspond to the starlight polarization directions, are presented in
Figure 19. The directions are found to be nearly parallel to the Galactic
plane. We then calculate the shock obliquity of each filament assuming that
the filaments have been moving straight against the explosion center estimated
in Section 4.3 (case (a)). The plots of the power-law parameters over the
shock obliquity are presented in Figure 20. No clear correlations are found.
Instead of assuming the simple filament motion against the explosion center,
we have also assumed that the shock normal directly corresponds to the moving
direction we defined in Figure 1 (case (b)). This assumption also results in
similar plots (Figure 20), without significant correlations.

Figure 19: Magnetic field directions estimated at the stars around RCW 86
region. $\Delta l^{\parallel}$, $\Delta b^{\parallel}$, and $P$ indicate the
Galactic longitude and latitude with respect to the coordinates
$(l,b)=(315\fdg 4015,-2\fdg 31664)$, and polarization degree, respectively.
Left and right panels show the estimates for distances of 2–3 kpc and 3–4 kpc,
respectively.  Figure 20: Plots of the power-law flux and index against the
shock obliquity. As a proxy of the shock obliquity, we assume two cases: we
use the angle made by the Galactic plane and vector connecting the filament
position and geometric center of RCW 86 (case (a)), and the angle made by the
Galactic plane and moving direction of the filament defined in Figure 1 (case
(b)). For the shock obliquity, a typical error of 1 degree is assigned to all
the regions.

## References

  * Allen et al. (2004) Allen, C., Jerius, D. H., & Gaetz, T. J. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5165, X-Ray and Gamma-Ray Instrumentation for Astronomy XIII, ed. K. A. Flanagan & O. H. W. Siegmund, 423–432, doi: 10.1117/12.509411
  * Amato & Blasi (2006) Amato, E., & Blasi, P. 2006, MNRAS, 371, 1251, doi: 10.1111/j.1365-2966.2006.10739.x
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Bamba et al. (2000) Bamba, A., Koyama, K., & Tomida, H. 2000, PASJ, 52, 1157, doi: 10.1093/pasj/52.6.1157
  * Bamba et al. (2005) Bamba, A., Yamazaki, R., Yoshida, T., Terasawa, T., & Koyama, K. 2005, ApJ, 621, 793, doi: 10.1086/427620
  * Bell (2004) Bell, A. R. 2004, MNRAS, 353, 550, doi: 10.1111/j.1365-2966.2004.08097.x
  * Dickel et al. (2001) Dickel, J. R., Strom, R. G., & Milne, D. K. 2001, ApJ, 546, 447, doi: 10.1086/318259
  * Fraschetti (2013) Fraschetti, F. 2013, ApJ, 770, 84, doi: 10.1088/0004-637X/770/2/84
  * Fraschetti et al. (2018) Fraschetti, F., Katsuda, S., Sato, T., Jokipii, J. R., & Giacalone, J. 2018, Phys. Rev. Lett., 120, 251101, doi: 10.1103/PhysRevLett.120.251101
  * Fruscione et al. (2006) Fruscione, A., McDowell, J. C., Allen, G. E., et al. 2006, in Observatory Operations: Strategies, Processes, and Systems, ed. D. R. Silva & R. E. Doxsey, Vol. 6270, International Society for Optics and Photonics (SPIE), 586 – 597, doi: 10.1117/12.671760
  * Garmire (1997) Garmire, G. P. 1997, in American Astronomical Society Meeting Abstracts, Vol. 190, American Astronomical Society Meeting Abstracts #190, 34.04
  * Giacalone & Jokipii (2007) Giacalone, J., & Jokipii, J. R. 2007, ApJ, 663, L41, doi: 10.1086/519994
  * Green & Stephenson (2003) Green, D. A., & Stephenson, F. R. 2003, Historical Supernovae, ed. K. Weiler, Vol. 598, 7–19, doi: 10.1007/3-540-45863-8_2
  * HEASARC (2014) HEASARC. 2014, HEAsoft: Unified Release of FTOOLS and XANADU.  http://ascl.net/1408.004
  * Heiles (2000) Heiles, C. 2000, AJ, 119, 923, doi: 10.1086/301236
  * Helder et al. (2012) Helder, E. A., Vink, J., Bykov, A. M., et al. 2012, Space Sci. Rev., 173, 369, doi: 10.1007/s11214-012-9919-8
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116, doi: 10.1051/0004-6361/201629178
  * Hickox & Markevitch (2006) Hickox, R. C., & Markevitch, M. 2006, ApJ, 645, 95, doi: 10.1086/504070
  * Inoue et al. (2012) Inoue, T., Yamazaki, R., Inutsuka, S.-i., & Fukui, Y. 2012, ApJ, 744, 71, doi: 10.1088/0004-637X/744/1/71
  * Katsuda et al. (2010) Katsuda, S., Petre, R., Mori, K., et al. 2010, ApJ, 723, 383, doi: 10.1088/0004-637X/723/1/383
  * Kishishita et al. (2013) Kishishita, T., Hiraga, J., & Uchiyama, Y. 2013, A&A, 551, A132, doi: 10.1051/0004-6361/201220525
  * Kuntz & Snowden (2000) Kuntz, K. D., & Snowden, S. L. 2000, ApJ, 543, 195, doi: 10.1086/317071
  * Kushino et al. (2002) Kushino, A., Ishisaki, Y., Morita, U., et al. 2002, PASJ, 54, 327, doi: 10.1093/pasj/54.3.327
  * Masui et al. (2009) Masui, K., Mitsuda, K., Yamasaki, N. Y., et al. 2009, PASJ, 61, S115, doi: 10.1093/pasj/61.sp1.S115
  * Pohl et al. (2005) Pohl, M., Yan, H., & Lazarian, A. 2005, ApJ, 626, L101, doi: 10.1086/431902
  * Ressler et al. (2014) Ressler, S. M., Katsuda, S., Reynolds, S. P., et al. 2014, ApJ, 790, 85, doi: 10.1088/0004-637X/790/2/85
  * Reynolds (2008) Reynolds, S. P. 2008, ARA&A, 46, 89, doi: 10.1146/annurev.astro.46.060407.145237
  * Reynolds et al. (2021) Reynolds, S. P., Williams, B. J., Borkowski, K. J., & Long, K. S. 2021, ApJ, 917, 55, doi: 10.3847/1538-4357/ac0ced
  * Rho et al. (2002) Rho, J., Dyer, K. K., Borkowski, K. J., & Reynolds, S. P. 2002, ApJ, 581, 1116, doi: 10.1086/344248
  * Rosado et al. (1996) Rosado, M., Ambrocio-Cruz, P., Le Coarer, E., & Marcelin, M. 1996, A&A, 315, 243
  * Sano et al. (2013) Sano, H., Tanaka, T., Torii, K., et al. 2013, ApJ, 778, 59, doi: 10.1088/0004-637X/778/1/59
  * Sano et al. (2015) Sano, H., Fukuda, T., Yoshiike, S., et al. 2015, ApJ, 799, 175, doi: 10.1088/0004-637X/799/2/175
  * Sano et al. (2017) Sano, H., Reynoso, E. M., Mitsuishi, I., et al. 2017, Journal of High Energy Astrophysics, 15, 1, doi: 10.1016/j.jheap.2017.04.002
  * Sano et al. (2019) Sano, H., Rowell, G., Reynoso, E. M., et al. 2019, ApJ, 876, 37, doi: 10.3847/1538-4357/ab108f
  * Sato et al. (2018) Sato, T., Katsuda, S., Morii, M., et al. 2018, ApJ, 853, 46, doi: 10.3847/1538-4357/aaa021
  * Smith (1997) Smith, R. C. 1997, AJ, 114, 2664, doi: 10.1086/118676
  * Snowden et al. (1992) Snowden, S. L., Plucinsky, P. P., Briel, U., Hasinger, G., & Pfeffermann, E. 1992, ApJ, 393, 819, doi: 10.1086/171549
  * Stephenson & Green (2002) Stephenson, F. R., & Green, D. A. 2002, Historical supernovae and their remnants, 5
  * Suzuki et al. (2020) Suzuki, H., Bamba, A., Yamazaki, R., & Ohira, Y. 2020, PASJ, 72, 72, doi: 10.1093/pasj/psaa061
  * Suzuki et al. (2022) —. 2022, ApJ, 924, 45, doi: 10.3847/1538-4357/ac33b5
  * Suzuki et al. (2021) Suzuki, H., Plucinsky, P. P., Gaetz, T. J., & Bamba, A. 2021, A&A, 655, A116, doi: 10.1051/0004-6361/202141458
  * Tran et al. (2015) Tran, A., Williams, B. J., Petre, R., Ressler, S. M., & Reynolds, S. P. 2015, ApJ, 812, 101, doi: 10.1088/0004-637X/812/2/101
  * Truelove & McKee (1999) Truelove, J. K., & McKee, C. F. 1999, ApJS, 120, 299, doi: 10.1086/313176
  * Tsubone et al. (2017) Tsubone, Y., Sawada, M., Bamba, A., Katsuda, S., & Vink, J. 2017, ApJ, 835, 34, doi: 10.3847/1538-4357/835/1/34
  * Tsuji et al. (2019) Tsuji, N., Uchiyama, Y., Aharonian, F., et al. 2019, ApJ, 877, 96, doi: 10.3847/1538-4357/ab1b29
  * Tsuji et al. (2021) Tsuji, N., Uchiyama, Y., Khangulyan, D., & Aharonian, F. 2021, ApJ, 907, 117, doi: 10.3847/1538-4357/abce65
  * Uchiyama et al. (2007) Uchiyama, Y., Aharonian, F. A., Tanaka, T., Takahashi, T., & Maeda, Y. 2007, Nature, 449, 576, doi: 10.1038/nature06210
  * Vink (2006) Vink, J. 2006, in ESA Special Publication, Vol. 604, The X-ray Universe 2005, ed. A. Wilson, 319.  https://arxiv.org/abs/astro-ph/0601131
  * Vink (2017) Vink, J. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 2063, doi: 10.1007/978-3-319-21846-5_92
  * Vink et al. (2006) Vink, J., Bleeker, J., van der Heyden, K., et al. 2006, ApJ, 648, L33, doi: 10.1086/507628
  * Völk et al. (2005) Völk, H. J., Berezhko, E. G., & Ksenofontov, L. T. 2005, A&A, 433, 229, doi: 10.1051/0004-6361:20042015
  * Williams et al. (2011) Williams, B. J., Blair, W. P., Blondin, J. M., et al. 2011, ApJ, 741, 96, doi: 10.1088/0004-637X/741/2/96
  * Yamaguchi et al. (2016) Yamaguchi, H., Katsuda, S., Castro, D., et al. 2016, ApJ, 820, L3, doi: 10.3847/2041-8205/820/1/L3
  * Yamaguchi et al. (2011) Yamaguchi, H., Koyama, K., & Uchida, H. 2011, PASJ, 63, S837, doi: 10.1093/pasj/63.sp3.S837
  * Yamazaki et al. (2014) Yamazaki, R., Ohira, Y., Sawada, M., & Bamba, A. 2014, Research in Astronomy and Astrophysics, 14, 165, doi: 10.1088/1674-4527/14/2/005
  * Yoshino et al. (2009) Yoshino, T., Mitsuda, K., Yamasaki, N. Y., et al. 2009, PASJ, 61, 805, doi: 10.1093/pasj/61.4.805
  * Yuan et al. (2014) Yuan, Q., Huang, X., Liu, S., & Zhang, B. 2014, ApJ, 785, L22, doi: 10.1088/2041-8205/785/2/L22
  * Zeng et al. (2019) Zeng, H., Xin, Y., & Liu, S. 2019, ApJ, 874, 50, doi: 10.3847/1538-4357/aaf392
  * Zirakashvili & Aharonian (2007) Zirakashvili, V. N., & Aharonian, F. 2007, A&A, 465, 695, doi: 10.1051/0004-6361:20066494


Topology-independent GEVD-based
distributed adaptive node-specific signal estimation
in ad-hoc wireless acoustic sensor networks

This research was carried out at the ESAT Laboratory of KU Leuven, in the frame of Research Council KU Leuven C14-21-0075 “A holistic approach to the design of integrated and distributed digital signal processing algorithms for audio and speech communication devices”, and was supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 956369: `Service-Oriented Ubiquitous Network-Driven Sound — SOUNDS'. The scientific responsibility is assumed by its authors. This paper reflects only the authors' views and the Union is not liable for any use that may be made of the contained information.

Paul Didier, Toon van Waterschoot, Marc Moonen
KU Leuven, Department of Electrical Engineering (ESAT),

STADIUS Center for Dynamical Systems, Signal Processing and Data Analytics, Belgium

Email: paul.didier, toon.vanwaterschoot<EMAIL_ADDRESS>
A low-rank approximation-based version of the
topology-independent distributed adaptive node-specific signal estimation (TI-DANSE) algorithm is introduced, using a generalized eigenvalue decomposition (GEVD) for application in ad-hoc wireless acoustic sensor networks. This TI-GEVD-DANSE algorithm as well as the original TI-DANSE algorithm exhibit a non-strict convergence, which can lead to numerical instability over time, particularly in scenarios where the estimation of accurate spatial covariance matrices is challenging. An adaptive filter coefficient normalization strategy is proposed to mitigate this issue and enable the stable performance of . The method is validated in numerical simulations including dynamic acoustic scenarios, demonstrating the importance of the additional normalization.

wireless acoustic sensor networks, distributed signal estimation, topology-independent, low-rank approximation

§ INTRODUCTION

In recent years, the ever-increasing ubiquity of multi-microphone devices capable of exchanging and processing acoustic signals has motivated the development of distributed audio signal processing algorithms.
As opposed to traditional localized microphone arrays, distributed systems do not rely on a fusion center; they instead leverage the computing capacities of each device (i.e., each node) in wasn. Distributed systems are typically able to use signals spanning a large acoustic area while maintaining a high degree of flexibility in their physical design [1, 2].

The tasks of distributed algorithms may be categorized depending on the estimated quantity [3].
Here, the focus is set on signal estimation for applications that require the retrieval of entire (possibly multichannel) signals of interest which may be node-specific and non-stationary, e.g., for noise reduction in speech enhancement tasks.
In a fully connected wasn, one may let all nodes transmit all their local sensor signals to all other nodes, effectively corresponding to a centralized case. For obvious reasons, this strategy suffers from an inefficient usage of communication bandwidth. Instead, it has been shown that the nodes can exchange fused versions of their local sensor signals while retaining the same performance as if they were transmitting all of their sensor signals. This is a core idea in the danse algorithm [4], which is considered in this paper.

The *danse algorithm is iterative and converges towards the centralized lmmse optimum. It operates in fully connected wasn, where every node can communicate with every other node. However, in many practical applications, there is no guarantee that the network topology will be fully connected or be static over time (e.g., due to link failures).
A solution to this is provided by the topology-independent (TI) danse algorithm (tidanse) [5], which allows a new tree to be pruned from the ad-hoc topology at any algorithm iteration while retaining convergence, thus being robust to dynamic topologies.

The *danse and *tidanse algorithms rely at their core on the computation of a mwf.
A significant performance improvement can be obtained when the number of latent target sources, or an estimate thereof, is known a priori – a reasonable assumption in many applications. The rank of the desired signal scm can be set equal to that number via a gevd. It has been shown that such gevd-based low-rank approximation of the mwf (gevdmwf) indeed outperforms the unconstrained mwf, particularly in challenging signal-to-noise ratio conditions [6]. Although the concept has been applied to danse in [7], leading to the gevddansegevddanse algorithm, its potential for  has remained unexplored.

The contribution of this paper is two-fold. First, we introduce the tigevddanse algorithm, which can operate in any topology and incorporates the advantages of a . Second, we address the observed non-strict convergence of the  algorithm, also observed in the original  algorithm, which allows the lmmse optimum to be reached even though the filters themselves diverge.
To this end, we propose a filter coefficient normalization strategy which stabilizes the behavior of ti_or_gevddanse even when signal statistics are estimated on the fly.

This paper is organized as follows. The problem statement is given in Section <ref>, defining the signal model and the centralized mwf solution. The  algorithm is reviewed in Section <ref>. The  algorithm is presented in Section <ref>. The proposed normalization strategy addressing the stability of  is presented in Section <ref>. Numerical experiments are provided in Section <ref>. Finally, conclusions are given in Section <ref>.

§ PROBLEM STATEMENT

We consider a wasn consisting of $K$ nodes, where node $k$ has $M_k$ sensors ($k\in\mathcal{K}:=\{1,...,K\}$) such that the total number of sensors is $M=\sum_{k\in\mathcal{K}}M_k$. The signals are assumed to be complex-valued to allow representation of, e.g., processing in a particular bin of a filter bank. The acoustic scenario includes $S$ desired sources whose latent signals at time $t$ are grouped in $\hat{\mathbf{s}}[t]\in\mathbb{C}^{S\times1}$, and $S_n$ noise sources with latent signals $\hat{\mathbf{n}}[t]\in\mathbb{C}^{S_n\times1}$. The sensor signals at all nodes are stacked in $\mathbf{y}[t]\in\mathbb{C}^{M\times 1}$, which is modeled as:

\begin{equation}\label{eq:signalmodel}
\mathbf{y}[t] = \mathbf{s}[t] + \mathbf{n}[t]
= \mathbf{A}\hat{\mathbf{s}}[t] + \mathbf{B}\hat{\mathbf{n}}[t] + \mathbf{q}[t],
\end{equation}

where $\mathbf{n}[t]$ is the noise component, $\mathbf{A}\in\mathbb{C}^{M\times S}$ and $\mathbf{B}\in\mathbb{C}^{M\times S_n}$ are the steering matrices of the desired signal and the noise, respectively, and $\mathbf{q}[t]$ is uncorrelated (thermal) noise. Time indices $[t]$ are omitted from here on for conciseness.
Each vector in (<ref>) can be partitioned as
$\mathbf{y} = [\mathbf{y}_1^\T\dots \mathbf{y}_K^\T]^\T$, where
$\cdot^\T$ denotes the transpose operator and $\mathbf{y}_k\in\mathbb{C}^{M_k\times 1}$ are the sensor signals of node $k$. Each node strives to estimate its own $J$-channel desired signal $\mathbf{d}_k = \mathbf{E}_{kk}^\T\mathbf{s}_k$ where $\mathbf{E}_{kk}$ is an $M_k\times J$ selection matrix.
A node-specific lmmse estimation problem is considered, where the optimal filter matrix $\mathbf{W}_k\in\mathbb{C}^{M\times J}$ used to estimate $\mathbf{d}_k$ from $\mathbf{y}$ is defined as:

\begin{equation}\label{eq:mmseCentr}
\mathbf{W}_k = \underset{\mathbf{W}}{\mathrm{arg\,min}}\:
\mathbb{E}\left\{
\left\|
\mathbf{d}_k - \mathbf{W}^\Her\mathbf{y}
\right\|^2_2
\right\},
\end{equation}

where $\mathbb{E}\{\cdot\}$, $\|\cdot\|_2$, and $\cdot^\Her$ denote the expected value operator, the Euclidean norm, and the Hermitian operator, respectively.
The solution of (<ref>) is the well-known mwf:

\begin{equation}
\mathbf{W}_k=\left(\mathbf{R}_{\mathbf{yy}}\right)^{-1}\mathbf{R}_{\mathbf{ss}}\mathbf{E}_{k},
\end{equation}

with the scm $\mathbf{R}_{\mathbf{yy}}=\mathbb{E}\{\mathbf{yy}^\Her\}$ and $\mathbf{R}_{\mathbf{ss}}=\mathbb{E}\{\mathbf{ss}^\Her\}$ and $\mathbf{E}_{k}$ the $M\times J$ selection matrix extracting $\mathbf{d}_k$ from $\mathbf{s}$. In practice, under the assumption that noise and desired signal components are uncorrelated, $\mathbf{R}_{\mathbf{ss}}$ may be estimated as $\mathbf{R}_{\mathbf{yy}} - \mathbf{R}_{\mathbf{nn}}$, where $\mathbf{R}_{\mathbf{nn}}=\mathbb{E}\{\mathbf{nn}^\Her\}$.
In real-world applications, scm must be estimated in an online fashion.
This can be achieved via, e.g., exponential averaging with a forgetting factor $0<\beta<1$ as:

\begin{equation}\label{eq:scmsOnline}
\begin{split}
\mathbf{R}_{\mathbf{y}\mathbf{y}}[t]
\beta
\mathbf{R}_{\mathbf{y}\mathbf{y}}[t-1]
+ (1 - \beta)
\mathbf{y}[t]\mathbf{y}^\Her[t],\\
\mathbf{R}_{\mathbf{n}\mathbf{n}}[t]
\beta
\mathbf{R}_{\mathbf{n}\mathbf{n}}[t-1]
+ (1 - \beta)
\mathbf{n}[t]\mathbf{n}^\Her[t],
\end{split}
\end{equation}

where the same symbols are used for estimated and true quantities for simplicity of notation.
In practice, $\mathbf{n}[t]$ may not be directly available and can be extracted via, e.g., an activity detector exploiting the ON-OFF structure of speech-like target signals [4, 8].
The desired signal estimate is finally obtained as $\hat{\mathbf{d}}_k = \mathbf{W}_k^\Her\mathbf{y}$.

In practice, making $\mathbf{y}$ available requires either a fusion center or the exchange of $M_k$-dimensional signals between the nodes of a fully connected wasn.

§ *TIDANSE

The danse algorithm [4] may be used to reach the centralized solution of (<ref>) while substantially reducing the communication bandwidth requirements. However, it requires a static and fully connected wasn topology, which is rarely the case in practice.
In an ad-hoc, possibly time-varying wasn topology, one may instead employ the  algorithm [5] which is reviewed in this section.

The  algorithm is iterative with iteration index $i$, i.e., it asymptotically converges towards the centralized *lmmse solution. At any iteration, the wasn topology can be pruned to a new tree using, e.g., Prim's algorithm. One node is chosen as the updating node, a natural choice being the root of the tree. The index of the updating node is set to cycle through all $k\in\mathcal{K}$ in a round-robin fashion, such that all nodes have updated once after $K$ iterations.
Each node $k$ defines its so-called fusion matrix $\mathbf{P}_k^i\in\mathbb{C}^{M_k\times J}$ and uses it to compute a $J$-channel fused version of its local sensor signals $\mathbf{z}_k^i = \mathbf{P}_k^{i\Her}\mathbf{y}_k$. All fused signals in the wasn are summed up via a sequence of partial in-network signal summations from the leaf nodes towards the root node. Once the full $J$-dimensional in-network sum $\boldsymbol{\eta}^i = \sum_{k\in\mathcal{K}}\mathbf{z}_k^i$ is built at the root node, it is flooded back towards the leaf nodes. Through this sequence of operations, each node has access to $D_k = M_k+J$ signals grouped in an observation vector:

\begin{equation}\label{eq:yTildeTI}
\tilde{\mathbf{y}}_k^i = [\mathbf{y}_k^\T\:|\: \boldsymbol{\eta}_{-k}^{i\T}]^\T
\:\:\:\text{where}\:\:\:
\boldsymbol{\eta}_{-k}^{i} = \boldsymbol{\eta}^{i} - \mathbf{z}_k^i = \sum_{q\in\mathcal{K}\backslash\{k\}}\mathbf{z}_q^i.
\end{equation}

The objective of node $k$ is then to solve its own node-specific *lmmse problem to estimate $\mathbf{d}_k$ from $\tilde{\mathbf{y}}_k^i$:

\begin{equation}\label{eq:mmseTIDANSE}
\tilde{\mathbf{W}}_k^{i+1} = \underset{\mathbf{W}}{\mathrm{arg\,min}}\:
\mathbb{E}\left\{
\left\|
\mathbf{d}_k - \mathbf{W}^\Her\tilde{\mathbf{y}}_k^i
\right\|^2_2
\right\},
\end{equation}

which, as in the centralized case, is solved by an mwf as:

\begin{equation}\label{eq:mwf_danse}
\tilde{\mathbf{W}}_k^{i+1} = (
\tilde{\mathbf{R}}_{\mathbf{y}_k\mathbf{y}_k}^i
\tilde{\mathbf{R}}_{\mathbf{s}_k\mathbf{s}_k}^i\tilde{\mathbf{E}}_{k},
\end{equation}

with the scm $\tilde{\mathbf{R}}_{\mathbf{y}_k\mathbf{y}_k}^i = \mathbb{E}\{\tilde{\mathbf{y}}_k^i\tilde{\mathbf{y}}_k^{i\Her}\}$ and $\tilde{\mathbf{R}}_{\mathbf{s}_k\mathbf{s}_k}^i = \mathbb{E}\{\tilde{\mathbf{s}}_k^i\tilde{\mathbf{s}}_k^{i\Her}\}$, $\tilde{\mathbf{y}}_k^i = \tilde{\mathbf{s}}_k^i + \tilde{\mathbf{n}}_k^i$, and $\tilde{\mathbf{E}}_{k} = [\mathbf{E}_{kk}^\T \:|\: \mathbf{0}_{J\times J}]^\T$.
A partitioning of $\tilde{\mathbf{W}}_k^i$ is defined here as $[\mathbf{W}_{kk}^{i\T} \:|\: \mathbf{G}_{k}^{i\T}]^\T$, where $\mathbf{W}_{kk}^{i}\in\mathbb{C}^{M_k\times J}$ is applied to the $\mathbf{y}_k$ while $\mathbf{G}_{k}^{i}\in\mathbb{C}^{J\times J}$ is applied to $\boldsymbol{\eta}_{-k}^i$.

In [5], the convergence and optimality of  is proved when the fusion rule is $\mathbf{P}_k^i = \mathbf{W}_{kk}^i(\mathbf{G}_{k}^i)^{-1}$, which completes the algorithm definition. It can be noted that the $(\mathbf{G}_k^i)^{-1}$ term serves to decouple the contributions of each individual node in $\boldsymbol{\eta}^i$. The  definition of $\mathbf{P}_k^i$ is a fundamental difference with that used in the danse algorithm, where the $(\mathbf{G}_k^i)^{-1}$ term is omitted [4].
As $i\to\infty$, the  algorithm converges to the centralized solution to (<ref>) while reducing the amount of information exchanged between nodes with respect to the centralized case, as well as to danse. The algorithm assumes that the nodes update sequentially, in a round-robin fashion. Simultaneous or asynchronous node-updating strategies are not considered here.

The mwf solving (<ref>) depends on $\mathbf{R}_\mathbf{ss}$, which is a rank-$S$ matrix since $\mathbf{R}_\mathbf{ss} = \mathbf{A}\mathbb{E}\{\hat{\mathbf{s}}\hat{\mathbf{s}}^\Her\}\mathbf{A}^\Her$ (cf. (<ref>)). However, in practice, $\mathbf{R}_\mathbf{ss}$ is not directly available and must be estimated, for instance as $\mathbf{R}_\mathbf{ss} = \mathbf{R}_\mathbf{yy} - \mathbf{R}_\mathbf{nn}$, often resulting in
a rank greater than $S$. In many practical applications such as speech enhancement, $S$ is known a priori or can be well-estimated. Constraining the rank of $\mathbf{R}_\mathbf{ss}$ can then be done via a gevd of the pencil $\{\mathbf{R}_\mathbf{yy}, \mathbf{R}_\mathbf{nn}\}$. This low-rank approximation has been shown to yield more robust performance in scenarios where accurate estimation of $\mathbf{R}_\mathbf{ss}$ is challenging.

Since the  algorithm is also based on the signal model from (<ref>), $\tilde{\mathbf{R}}_{\mathbf{s}_k\mathbf{s}_k}^i$ is also rank-$S$. This can be seen from the definitions below (<ref>), rewriting:

\begin{gather}
\tilde{\mathbf{R}}_{\mathbf{s}_k\mathbf{s}_k}^i
= \mathbb{E}\left\{
\tilde{\mathbf{s}}_k^i\tilde{\mathbf{s}}_k^{i\Her}
\right\}
= \mathbb{E}\left\{
\mathbf{C}_k^{i\Her}\mathbf{s}\mathbf{s}^\Her\mathbf{C}_k^{i}
\right\}
= \mathbf{C}_k^{i\Her} \mathbf{R}_{\mathbf{ss}} \mathbf{C}_k^{i},\\
\text{where:}\:\:\:
\mathbf{C}_k^{i\Her} = \begin{bmatrix}
\mathbf{0} \:\:\:\: \cdots \:\:\:\:\: \mathbf{0} \:\:\:\: \mathbf{I}_{M_k} \:\:\: \mathbf{0} \:\:\:\: \cdots \:\:\:\: \mathbf{0}\\
\mathbf{P}_1^{i\Her} \: \cdots \: \mathbf{P}_{k-1}^{i\Her} \:\: \mathbf{0} \:\:\: \mathbf{P}_{k+1}^{i\Her} \: \cdots \: \mathbf{P}_{K}^{i\Her}
\end{bmatrix}.
\end{gather}

Since $\mathbf{R}_{\mathbf{ss}}$ is rank-$S$, so is $\tilde{\mathbf{R}}_{\mathbf{s}_k\mathbf{s}_k}^i$. This low-rank property can be guaranteed by first making use of a gevd on the pencil $\{
\tilde{\mathbf{R}}_{\mathbf{y}_k\mathbf{y}_k}^i,
\tilde{\mathbf{R}}_{\mathbf{n}_k\mathbf{n}_k}^i
\}$ to rewrite the scm as:

\begin{equation}\label{eq:GEVDoutcome}
\tilde{\mathbf{R}}_{\mathbf{y}_k\mathbf{y}_k}^i
= \tilde{\mathbf{Q}}_k^i \mathbf{\Sigma}_k^i \tilde{\mathbf{Q}}_k^{i\Her}
\:\:\:\text{and}\:\:\:
\tilde{\mathbf{R}}_{\mathbf{n}_k\mathbf{n}_k}^i
= \tilde{\mathbf{Q}}_k^i \tilde{\mathbf{Q}}_k^{i\Her},
\end{equation}

where $\mathbf{\Sigma}_k^i = \mathrm{diag}\{
\tilde{\sigma}_{k1}^i,...,\tilde{\sigma}_{kD_k}^i
\}$ contains the generalized eigenvalues ordered from largest to smallest and $\tilde{\mathbf{Q}}_k^i$ contains the corresponding generalized eigenvectors.
The rank can be constrained to $R$ (ideally equal to $S$, when $S$ is known a priori) by setting the $D_k - R$ smallest eigenvalues in $\tilde{\mathbf{R}}_{\mathbf{y}_k\mathbf{y}_k}^i$ to zero, obtaining an estimate of $\tilde{\mathbf{R}}_{\mathbf{s}_k\mathbf{s}_k}$ as:

\begin{equation}\label{eq:RssGEVDest}
\hat{\tilde{\mathbf{R}}}_{\mathbf{s}_k\mathbf{s}_k}
= \hat{\tilde{\mathbf{R}}}_{\mathbf{y}_k\mathbf{y}_k}^i - \tilde{\mathbf{R}}_{\mathbf{n}_k\mathbf{n}_k}^i
\tilde{\mathbf{Q}}_k^i \mathbf{\Delta}_k^i \tilde{\mathbf{Q}}_k^{i\Her},
\end{equation}

where $\mathbf{\Delta}_k^i = \mathrm{diag}\{
\tilde{\sigma}_{k1}^i-1, \dots, \tilde{\sigma}_{kR}^i-1, 0, \dots, 0
\}$.
This results in the  algorithm, which solution to (<ref>) is obtained by substituting (<ref>) into (<ref>):

\begin{equation}\label{eq:filterUpdateTIGEVDDANSE}
\tilde{\mathbf{W}}_k^{i+1}
\end{equation}

where $\tilde{\mathbf{\Lambda}}_k^i = \mathrm{diag}\{
1-1/\tilde{\sigma}_{k1}^i \dots 1-1/\tilde{\sigma}_{kR}^i, 0\dots 0
\}$. The rest of the algorithm remains unchanged with respect to Section <ref>. In addition to enabling enhanced and more robust signal estimation performance, it is observed that this rank-$R$ approximation preserves convergence of *tigevddanse even in cases where $J$ or $R$ underestimates the number of latent desired sources $S$. This is not the case of *tidanse which relies on the assumption $S\leq J$. Although a convergence proof is not provided in this paper, the observed property is illustrated via the experimental results showcased in Section <ref>.

§ IMPROVED ROBUSTNESS

§.§ Non-Strict Convergence of

It is observed that the  algorithm exhibits a non-strict convergence due to the formulation of its fusion matrices.
As derived in Section <ref>, $\tilde{\mathbf{y}}_k^i = \mathbf{C}_k^{i\Her}\mathbf{y}$, implying that the desired signal estimate can be written as $\hat{\mathbf{d}}_k^{i+1} = \tilde{\mathbf{W}}_k^{i+1,\Her}\tilde{\mathbf{y}}_k^i
= \tilde{\mathbf{W}}_k^{i+1,\Her}\mathbf{C}_k^{i\Her}\mathbf{y}$.
Therefore, the $M \times J$ matrix $\mathbf{W}_k^{i+1} = \mathbf{C}_k^{i}\tilde{\mathbf{W}}_k^{i+1}$ is the network-wide version of the  filter matrix, with structure:

\begin{equation}\label{eq:netWideTIDANSEmat}
\mathbf{W}_k^{i+1}
\left[
\cdots\:|\:
\right]^\T\\
\mathbf{W}_{kk}^{i+1}\\
\left[
\cdots\:|\:
\right]^\T
\end{bmatrix}.
\end{equation}

Inspecting the terms $\mathbf{P}_q^{i}\mathbf{G}_k^{i+1} = \mathbf{W}_{qq}^{i}(\mathbf{G}_q^{i})^{-1}\mathbf{G}_k^{i+1}$ in (<ref>) reveals that convergence of $\mathbf{W}_k^i$ (and thus $\tilde{\mathbf{W}}_k^i$) does not require the filter coefficients in the $\{\mathbf{G}_k^i\}_{k\in\mathcal{K}}$ matrices to strictly converge. In fact, only the strict convergence of $(\mathbf{G}_q^i)^{-1}\mathbf{G}_k^{i+1}\:\forall\:(k,q)\in\mathcal{K}\times\mathcal{K}\backslash\{k\}$ is necessary. In this state, the  algorithm allows the elements of the $\{\mathbf{G}_k^i\}_{k\in\mathcal{K}}$ matrices to grow infinitely large or small as $i$ increases. Such behavior is indeed observable in practice and can lead to numerical overflow or significant precision errors.

§.§ Normalization strategy

In this section, we address the non-strict convergence of .
Let us define a normalization factor $\gamma^i\in\mathbb{C}$ and suppose that every node simultaneously start normalizing their $\mathbf{G}_k^i$ matrix at iteration $i$. The normalized fused signal $\bar{\mathbf{z}}_k^i$ can be related to its non-normalized counterpart via:

\begin{equation}\label{eq:z_norm}
\bar{\mathbf{z}}_k^{i} =
\mathbf{W}_{kk}^{i}(
\mathbf{G}_k^{i}
\gamma^{i}
\mathbf{y}_k
\gamma^{i,\ast}
\mathbf{P}_k^{i\Her}
\mathbf{y}_k
\gamma^{i,\ast}
\mathbf{z}_k^i.
\end{equation}

where $\cdot^\ast$ denotes the complex conjugate.
Following the same notation logic, $\bar{\boldsymbol{\eta}}^i = \sum_{k\in\mathcal{K}} \bar{\mathbf{z}}_k^{i} = \gamma^{i,\ast}\boldsymbol{\eta}^i$, meaning that the normalized version of $\tilde{\mathbf{y}}_k^i$ can be expressed as $\bar{\tilde{\mathbf{y}}}_k^i = \mathbf{N}_k^i\tilde{\mathbf{y}}_k^i$ where $\mathbf{N}_k^i = \mathrm{blkdiag}\{\mathbf{I}_{M_k}, \gamma^i\mathbf{I}_{J}\}$.
It follows that:

\begin{equation}
\bar{\tilde{\mathbf{R}}}_{\mathbf{y}_k\mathbf{y}_k}^i
= \mathbb{E}\left\{\bar{\tilde{\mathbf{y}}}_k^i\bar{\tilde{\mathbf{y}}}_k^{i\Her}\right\}
= \mathbf{N}_k^i\tilde{\mathbf{R}}_{\mathbf{y}_k\mathbf{y}_k}^i\mathbf{N}_k^{i\Her},
\end{equation}

and likewise for $\bar{\tilde{\mathbf{R}}}_{\mathbf{n}_k\mathbf{n}_k}^i$. Substituting in (<ref>) gives:

\begin{equation}
\bar{\tilde{\mathbf{R}}}_{\mathbf{y}_k\mathbf{y}_k}^i
= \bar{\tilde{\mathbf{Q}}}_k^i \mathbf{\Sigma}_k^i \bar{\tilde{\mathbf{Q}}}_k^{i\Her}
\:\:\:\text{and}\:\:\:
\bar{\tilde{\mathbf{R}}}_{\mathbf{n}_k\mathbf{n}_k}^i
= \bar{\tilde{\mathbf{Q}}}_k^i \bar{\tilde{\mathbf{Q}}}_k^{i\Her},
\end{equation}

where $\bar{\tilde{\mathbf{Q}}}_k^i = \mathbf{N}_k^i\tilde{\mathbf{Q}}_k^i$. Consequently, (<ref>) gives:

\begin{equation}
\bar{\tilde{\mathbf{W}}}_k^{i+1}
\tilde{\mathbf{W}}_k^{i+1},
\end{equation}

where use is made of the fact that the desired signal only includes contributions from local sensors, i.e., $\mathbf{N}_k^{i\Her}
\tilde{\mathbf{E}}_{k} = \tilde{\mathbf{E}}_{k}$.
It can be shown that the normalization does not alter the network-wide filters of (<ref>). This can be seen from:

\begin{equation}\label{eq:d_withNorm}
\bar{\hat{\mathbf{d}}}_k^{i+1} = \bar{\tilde{\mathbf{W}}}_k^{i+1,\Her}\bar{\tilde{\mathbf{y}}}_k^i
\tilde{\mathbf{W}}_k^{i+1,\Her}
\mathbf{N}_k^i
\tilde{\mathbf{y}}_k^i
= \hat{\mathbf{d}}_k^{i+1}.
\end{equation}

This normalization procedure can be incorporated in , resulting in the following algorithm:

[1]
Initialize $u=0$, $\gamma^{1} = 1$, $r\in\mathcal{K}$.
Form tree topology rooted at node $u$.
At all $k\in\mathcal{K}$, compute $\bar{\mathbf{z}}_k^{i} = \gamma^{i,\ast} \bar{\mathbf{P}}_k^{i\Her} \mathbf{y}_k$.
At node $u$, compute the in-network sum $\boldsymbol{\eta}^{i}$ and flood it back through the wasn.
Build observation vector $\bar{\tilde{\mathbf{y}}}_k^i$ as in (<ref>).
Compute scm $\bar{\tilde{\mathbf{R}}}_{\mathbf{y}_k\mathbf{y}_k}^{i}$ and $\bar{\tilde{\mathbf{R}}}_{\mathbf{n}_k\mathbf{n}_k}^{i}$.
Compute $\bar{\tilde{\mathbf{W}}}_k^{i+1}=[\bar{\mathbf{W}}_{kk}^{i+1,\T}\:|\:\bar{\mathbf{G}}_{k}^{i+1,\T}]^\T$ via (<ref>).
$k\neq u$
Compute $\bar{\tilde{\mathbf{W}}}_k^{i+1} = \left(\mathbf{N}_k^i\right)^{-\Her} \bar{\tilde{\mathbf{W}}}_k^{i}$.
Update fusion matrix as $\bar{\mathbf{P}}_k^{i+1} = \bar{\mathbf{W}}_{kk}^{i+1}(\bar{\mathbf{G}}_{k}^{i+1})^{-1}$.
At reference node $r$, compute $\gamma^{i+1} = \|\bar{\mathbf{G}}_r^{i+1}\|_F$ and flood it through the *wasn.
At all $k\in\mathcal{K}$, compute $\hat{\mathbf{d}}_k^{i+1} =
\bar{\tilde{\mathbf{W}}}_k^{i+1,\Her}\bar{\tilde{\mathbf{y}}}_k^{i+1}$.
$u \gets (u + 1) \mod K$.

Remark 1:
The choice for $\gamma^i$ is to use $\|\bar{\mathbf{G}}_r^{i+1}\|_F$ for a fixed reference node $r$ (step 15). It should be noted that the algorithm can in principle select any node $r$ as reference, as long as the normalization is equal for all nodes.

Remark 2:
Step 12 ensures consistency across iterations for non-updating nodes $k\neq u$. Indeed, the $\gamma^{i,\ast}$ factor from step 4 is carried along in the scm and thus impacts the filter of the updating node. Conversely, a non-updating node will keep its filters from the previous iteration, meaning that the normalization effect must be accounted for via $(\mathbf{N}_k^i)^{-\Her}$.

Remark 3:  The communication bandwidth increase generated by the broadcasting of $\gamma^{i+1}$ (step 15) is negligible in comparison to the unnormalized  algorithm, as it represents the exchange of a single scalar at most at every iteration. Depending on the severity of the non-strictly convergent behavior of the $\{\bar{\mathbf{G}}_k^i\}_{k\in\mathcal{K}}$ coefficients, it may even be sufficient to update (and thus broadcast) $\gamma$ less frequently.

When estimating scm in an online fashion based on (<ref>), the normalization must be accounted for by setting:

\begin{equation}\label{eq:scmsOnlineWithNorm}
\begin{split}
\bar{\tilde{\mathbf{R}}}_{\mathbf{y}_k\mathbf{y}_k}^i
\beta
\mathbf{N}_k^i \bar{\tilde{\mathbf{R}}}_{\mathbf{y}_k\mathbf{y}_k}^{i-1} \mathbf{N}_k^{i\Her}
+ (1 - \beta)
\bar{\tilde{\mathbf{y}}}_k^i\bar{\tilde{\mathbf{y}}}_k^{i\Her},\\
\bar{\tilde{\mathbf{R}}}_{\mathbf{n}_k\mathbf{n}_k}^i
\beta
\mathbf{N}_k^i \bar{\tilde{\mathbf{R}}}_{\mathbf{n}_k\mathbf{n}_k}^{i-1} \mathbf{N}_k^{i\Her}
+ (1 - \beta)
\bar{\tilde{\mathbf{y}}}_k^i\bar{\tilde{\mathbf{y}}}_k^{i\Her}.
\end{split}
\end{equation}

where the index $i$ here represents both the iteration index and the time index, for simplicity.

§ NUMERICAL EXPERIMENTS

The performance of  with and without normalization is assessed via simulations in an acoustic environment composed of an ad-hoc non-fully connected wasn and 6 localized sources, $S{=}3$ of which are considered as targets and the $S_n{=}3$ others as noise.
The number of nodes is fixed but the specific wasn topology does not have to be, as it is simply assumed that $\boldsymbol{\eta}^i$ is available at all nodes at any $i$.
Although the gevd-rank $R$ can be chosen independently from the number of exchanged channels $J$, we here set $R{=}J$ for simplicity (other cases are discussed in, e.g., [7]).
The uncorrelated (thermal) noise at node $k$ is set to have a power equal to 10% of the power of the combined target source signals as observed by the first sensor of the node. Signal samples and steering matrices entries are drawn from the uniform distribution over .

§.§ Batch-mode simulations

Batch-mode simulations without normalization ($\gamma^i = 1\:\forall\:i$) are first performed with $K{=}5$ and $M_k{=}4$ for all $k$ to demonstrate the convergence of  towards the centralized . The number of observations is set to $N{=}15000$ samples. The scm are estimated as sample means. The results are averaged over 3 runs with different random steering matrices.
In order to isolate scm estimation errors, $\mathbf{R}_\mathbf{nn}$ and $\bar{\tilde{\mathbf{R}}}^i_{\mathbf{n}_k\mathbf{n}_k}$ are estimated using oracle knowledge of the noise-only signals.
The performance of the  algorithm is shown for different values of $R{=}J$ in Fig. <ref>, where two mse metrics are used. The first, $\mathrm{MSE}_{W_k}^i$, is defined between the network-wide expansion of the  solution of (<ref>) and the centralized  solution of (<ref>). The second, $\mathrm{MSE}_{d_k}^i$, is defined between the true desired signal and its estimate $\hat{\mathbf{d}}_k^i$:

\begin{align}
\mathrm{MSE}_{W_k}^i &= \frac{1}{MJ}\|
\bar{\mathbf{W}}_k^i - \mathbf{W}_k
\|_F^2,\\
\mathrm{MSE}_{d_k}^i &= \frac{1}{JN}\sum_{n=0}^{N-1}
\|
\hat{\mathbf{d}}_k^i[n] - \mathbf{d}_k[n]
\|_F^2.
\end{align}

Batch processing results. The $x$-axis represents the iteration index $i$. Top: $\mathrm{MSE}_{d_k}^i$ averaged over nodes, where the horizontal dotted lines represent the corresponding centralized values. Bottom: $\mathrm{MSE}_{W_k}^i$ averaged over nodes.

The convergence of the  algorithm towards the centralized  is clearly visible from the average $\mathrm{MSE}_{W_k}^i$, which would decrease to machine precision for $i{\to}\infty$. Larger values of $R{=}J$ show a faster convergence since more distinct signals are available at each nodes. It is important to notice that convergence is preserved even in cases where $R{=}J{<}S$, which is not the case for the original *tidanse [5]. The average $\mathrm{MSE}_{d_k}^i$ shows that  also matches the target signal estimation performance of the  with corresponding rank, and that a rank $R$ closer (or equal) to $S$ yields better performance.

§.§ Online-mode simulations

The behavior of  is assessed in an online processing scenario with and without normalization, with $K{=}10$ and $M_k{=}15$ for all $k$. The scm estimation strategy defined in (<ref>) is used with $\beta{=}0.7$ and all scm entries randomly initialized.
At every iteration, a new frame of $B=500$ samples is drawn from the uniform distribution over  for each sound source.
A dynamic scenario is created where the entries of steering matrices $\mathbf{A}$ and $\mathbf{B}$ defined in (<ref>) have, at each time frame, a 0.05 probability to be changed by drawing entries from the same uniform distribution over , thus changing the relative positioning of sensors and sources. After an acoustic scenario change, the probability is set to 0 for 30 frames before coming back to 0.05.

The quantity $\mathrm{MSE}_{d_k}^i$ averaged over all nodes is shown in Fig. <ref> for  with and without normalization and for the centralized , all with $J{=}R{=}S{=}3$.
To highlight the effect of $\gamma^i$, the average over all nodes of $\|\bar{\mathbf{G}}_k^{i}\|_F$ is also shown for  with and without normalization. The reference node index $r$ is kept equal to 1 through the entire simulation.

Online processing results. The $x$-axis represents the index $i$. Top: $\mathrm{MSE}_{d_k}^i$ averaged over nodes. Bottom: $\|\bar{\mathbf{G}}_k^i\|_F$ averaged over nodes.

The results show that, as long as  is stable without normalization, the normalization has no impact on the performance as per (<ref>), even in a dynamic scenario. It can be noticed that a change in the steering matrices triggers a new adaptation phase where the algorithm must  the scm to be able to  towards the centralized  (visible as peaks in the average $\mathrm{MSE}_{d_k}^i$).

The importance of normalization is clearly visible on the lower plot. Without normalization, the non-strict convergence of  leads to a decreasing $\|\bar{\mathbf{G}}_k^i\|_F = \|\mathbf{G}_k^i\|_F$. Steeper decreases occur when the acoustic scenario is changed, as the scm must be re-estimated. This behavior results in numerical instability over time as can be seen from $i=488$ in the upper plot. With normalization, however, the average norm of $\bar{\mathbf{G}}_k^i$ remains close to 1 even if the acoustic scenario changes, which allows  to indefinitely perform without numerical overflow or critical precision errors.

§ CONCLUSION

In this paper, we have introduced a gevd-based  algorithm () which provides an enhanced and robust performance in scenarios where the accurate estimation of the desired-signal scm is challenging. An adaptive normalization procedure has been included, ensuring the stability of the estimated filter coefficients through time, even in online processing with time-varying acoustic scenarios.

[1] A. Bertrand, “Applications and trends in wireless acoustic sensor networks: A signal processing perspective,” Proc. 18th IEEE Symp. Commun. Veh. Tech., 2011, pp.1–6.
[2] A. Boukerche and P. Sun, “Design of Algorithms and Protocols for Underwater Acoustic Wireless Sensor Networks,” ACM Comput. Surv., vol. 53, no. 6, art. 134, pp 1–34, Nov. 2021.
[3] S. X. Wu, H.-T. Wai, L. Li and A. Scaglione, "A Review of Distributed Algorithms for Principal Component Analysis," in Proc. of the IEEE, vol. 106, no. 8, pp. 1321–1340, Aug. 2018.
[4] A. Bertrand and M. Moonen, “Distributed adaptive node-specific signal estimation in fully connected sensor networks–Part I: Sequential node updating,” IEEE Trans. Signal Process., vol. 58, no. 10, pp. 5277–5291, Oct. 2010.
[5] J. Szurley, A. Bertrand, and M. Moonen, “Topology-independent distributed adaptive node-specific signal estimation in wireless sensor networks,” IEEE Trans. Signal Inf. Process. Netw., vol. 3, no. 1, pp. 130–144, Mar. 2017.
[6] R. Serizel, M. Moonen, B. Van Dijk, and J. Wouters, “Low-rank approximation based multichannel Wiener filter algorithms for noise reduction with application in cochlear implants,” IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 22, no. 4, pp.785-799, Feb. 2014.
[7] A. Hassani, A. Bertrand, and M. Moonen, “GEVD-based low-rank approximation for distributed adaptive node-specific signal estimation in wireless sensor networks,” IEEE Trans. Signal Process., vol. 64, no. 10, pp. 2557–2572, May 2016.
[8] Y. Zhao, J. K. Nielsen, J. Chen, and M. G. Christensen, “Model-based distributed node clustering and multi-speaker speech presence probability estimation in wireless acoustic sensor networks,” J. Acoust. Soc. Amer., vol. 147, no. 6, pp. 4189–4201, 2020.