# High-resolution Spectroscopic Metallicities of Milky Way Cepheid Standards
and their impact on the Leavitt Law and the Hubble constant

Anupam Bhardwaj INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello
16, 80131 Napoli, Italy  Adam G. Riess Space Telescope Science Institute, 3700
San Martin Drive, Baltimore, MD 21218, USA  Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  Giovanni
Catanzaro  INAF-Osservatorio Astrofisico di Catania, Via S.Sofia 78, 95123,
Catania, Italy  Erasmo Trentin  Leibniz-Institut für Astrophysik Potsdam
(AIP), An der Sternwarte 16, D-14482 Potsdam, Germany  Institut für Physik und
Astronomie, Universität Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, D-14476
Golm (Potsdam), Germany  INAF-Osservatorio Astronomico di Capodimonte, Via
Moiariello 16, 80131 Napoli, Italy  Vincenzo Ripepi  INAF-Osservatorio
Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy  Marina
Rejkuba European Southern Observatory, Karl-Schwarzschild-Straße 2, 85748,
Garching, Germany  Marcella Marconi INAF-Osservatorio Astronomico di
Capodimonte, Via Moiariello 16, 80131 Napoli, Italy  Chow-Choong Ngeow
Graduate Institute of Astronomy, National Central University, 300 Jhongda
Road, 32001 Jhongli, Taiwan  Lucas M. Macri NSF’s NOIRLab, 950 N Cherry Ave,
Tucson, AZ 85719  Martino Romaniello  European Southern Observatory, Karl-
Schwarzschild-Straße 2, 85748, Garching, Germany  Roberto Molinaro  INAF-
Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli,
Italy  Harinder P. Singh Department of Physics and Astrophysics, University of
Delhi, Delhi-110007, India  Shashi M. Kanbur  Department of Physics, State
University of New York, Oswego, NY 13126, USA

(Accepted Sep 5, 2023)

###### Abstract

Milky Way Cepheid variables with accurate Hubble Space Telescope photometry
have been established as standards for primary calibration of the cosmic
distance ladder to achieve a percent-level determination of the Hubble
constant ($H_{0}$). These 75 Cepheid standards are the fundamental sample for
investigation of possible residual systematics in the local $H_{0}$
determination due to metallicity effects on their period-luminosity relations.
We obtained new high-resolution ($R\sim 81,000$), high signal-to-noise
($S/N\sim 50-150$) multi-epoch spectra of 42 out of 75 Cepheid standards using
ESPaDOnS instrument at the 3.6-m Canada-France-Hawaii Telescope. Our
spectroscopic metallicity measurements are in good agreement with the
literature values with systematic differences up to $0.1$ dex due to different
metallicity scales. We homogenized and updated the spectroscopic metallicities
of all 75 Milky Way Cepheid standards and derived their multiwavelength
($GVIJHK_{s}$) period-luminosity-metallicity and period-Wesenheit-metallicity
relations using the latest Gaia parallaxes. The metallicity coefficients of
these empirically calibrated relations exhibit large uncertainties due to low
statistics and a narrow metallicity range ($\Delta\textrm{[Fe/H]}=0.6$ dex).
These metallicity coefficients are up to three times better constrained if we
include Cepheids in the Large Magellanic Cloud and range between $-0.21\pm
0.07$ and $-0.43\pm 0.06$ mag/dex. The updated spectroscopic metallicities of
these Milky Way Cepheid standards were used in the Cepheid-Supernovae distance
ladder formalism to determine $H_{0}=72.9~{}\pm 1.0$ km s-1 Mpc-1, suggesting
little variation ($\sim 0.1$ km s-1 Mpc-1) in the local $H_{0}$ measurements
due to different Cepheid metallicity scales.

††journal: ApJL††software: The IDL Astronomy User’s Library (Landsman, 1993),
Astropy (Astropy Collaboration et al., 2013, 2018, 2022), Libre-ESpRIT (Donati
et al., 1997)††thanks: Marie-Curie Fellow

##  1 Introduction

Classical Cepheids play a vital role in calibrating the first-rung of the
cosmic distance ladder thanks to their well-known Period-Luminosity (PL)
relation (or the “Leavitt Law”, Leavitt & Pickering, 1912), which allows an
accurate and precise measurement of the Hubble constant ($H_{0}$) in the local
universe (Freedman et al., 2001; Riess et al., 2016, 2022a). At present, a
$\sim$1.4% measurement of $H_{0}$ in the late universe based on Cepheids and
Type Ia Supernovae (SNe) is significantly ($\sim$5.3$\sigma$) larger than the
precise Planck inference based on observations of the early universe (Planck
Collaboration et al., 2020; Riess et al., 2022b). This discrepancy between
direct and cosmology-based $H_{0}$ measurements may lead to new physics beyond
the standard model (Di Valentino et al., 2021). Reaching an accurate, percent-
level measurement of $H_{0}$ while investigating all possible sources of
systematic uncertainties in Cepheid-based measurements is now critical
considering this Hubble tension (see Verde et al., 2019, for details).

In the traditional Cepheid-SNe distance ladder from the SH0ES project
(Supernovae and $H_{0}$ for the equation of state of dark energy; Riess et
al., 2011, 2016, 2022a), the main Cepheid-based error budget on $H_{0}$ comes
from three anchors with geometric distances but receives the highest weight
from Gaia Milky Way (MW) parallaxes (Riess et al., 2022a). Since the
calibration of the Galactic Cepheid PL relation is now approaching percent-
level precision (Riess et al., 2022b; Cruz Reyes & Anderson, 2023), it is
crucial to investigate systematic uncertainties in these empirical
calibrations due to metallicity (Riess et al., 2021; Ripepi et al., 2021;
Breuval et al., 2022; Molinaro et al., 2023; Owens et al., 2022).
Specifically, given that the Cepheid-only error budget in $H_{0}$ measurements
is below $1\%$ (Riess et al., 2022a), even modest differences due to past,
heterogeneous metallicities or metallicity scales can challenge the goal of
achieving the ultimate percent-level precision.

High-resolution spectroscopic metallicities for several bright solar
neighbourhood Cepheids are available in the literature that were primarily
used to investigate the abundance gradient in the thin disk and/or the impact
of metallicity on PL relations (e.g., Genovali et al., 2014, 2015; Luck, 2018;
Ripepi et al., 2021; da Silva et al., 2022; Kovtyukh et al., 2022).
Complementing metal-rich Cepheids in the solar vicinity, Trentin et al. (2023)
provided high-resolution spectra for the most metal-poor Galactic Cepheids
that will be useful for an improved quantification of the metallicity
dependence of PL relations. Some of these aforementioned studies also used
multi-epoch spectra for a smaller sample of Cepheids to homogenize the
metallicity scale of a larger sample of literature measurements (e.g.,
Genovali et al., 2014, 2015; Luck, 2018). The literature compilations of a few
hundred Cepheids allow better statistical constraints on the metallicity
dependence of PL relations, but cannot provide a proper account of systematic
uncertainties due to heterogeneous metallicities obtained by employing diverse
methodologies on spectra from different instruments/telescopes in different
spectral ranges. Although there are more than 600 MW Cepheids with high-
resolution spectra (Trentin et al., 2023), only 75 bright nearby Cepheids have
both the high-resolution spectroscopy and accurate Hubble Space Telescope
(HST) photometry (Riess et al., 2021).

Riess et al. (2018a, 2021) established 75 MW Cepheid standards with HST
photometry for calibration of the first-rung of the distance ladder. The HST
photometry for these bright Cepheids was carried out in an observationally
rigorous and specifically designed spatial scanning mode to mitigate
saturation and reduce pixel-to-pixel calibration errors (Riess et al., 2018a).
The HST photometry played a crucial role to establish the set of MW Cepheid
standards calibrated on the same photometric system used to measure their
extragalactic counterparts. The photometric calibration minimized systematic
uncertainties due to different instrumental calibrations and/or ground-to-
space based photometric transformations. However, the spectroscopic data used
by Riess et al. (2021) for the calibration of the $W^{H}_{m}$ Period-
Wesenheit-Metallicity (PWZ) relation in HST filters come from a compilation of
literature measurements. We have therefore started an observational project to
collect high-resolution and high signal-to-noise spectra obtaining homogeneous
spectroscopic metallicities for all these MW Cepheid standards. Our goal is to
derive multiband Period-Luminosity-Metallicity (PLZ) relations and investigate
possible small systematic uncertainties in the local $H_{0}$ determination due
to metallicity dependence of Cepheid PL relations.

##  2 The Spectroscopic and photometric data

Figure 1: High-resolution ESPaDOnS spectra for all 42 MW Cepheid standards in
the range 6000-6050 Å. The star IDs are provided on the right and the vertical
dashed lines mark the positions of atomic lines.

###  2.1 New high-resolution spectroscopic metallicities

We obtained homogeneous high-resolution ($R\sim 81,000$), high signal-to-noise
($S/N\sim 50-150$) optical spectra ($370~{}\textrm{to}~{}1050~{}\textrm{nm}$)
for 42 Cepheids that are observable with the ESPaDOnS instrument at 3.6-m
Canada-France-Hawaii Telescope (CFHT). The observations were carried out in
star-only mode and up to three spectra were obtained for a given target on 13
different nights across six runs between February 2021 and September 2022. A
total of 90 pre-processed spectra (bias subtraction, spectrum extraction, flat
fielding and wavelength calibration) were generated with the Libre-
ESpRIT111https://www.cfht.hawaii.edu/Instruments/Spectroscopy/Espadons/Espadons_esprit.html
data reduction package for reducing echelle spectropolarimetric data at CFHT.

The iron abundances were measured using the methodology discussed in detail in
Ripepi et al. (2021) and Trentin et al. (2023). In brief, the first step was
to measure the atmospheric parameters — effective temperature, surface
gravity, microturbulent velocity and the line broadening parameter. The
effective temperature was measured using the line-depth ratio (LDR) method
(Kovtyukh & Gorlova, 2000). The microturbulent velocities, iron abundances,
and gravities were then determined through an iterative procedure in which it
was assumed that the iron abundances do not depend on the equivalent widths.
The equivalent widths of a sample of 145 Fe i@ lines published by Romaniello
et al. (2008) were used for this purpose. The equivalent widths were converted
to abundances through the WIDTH9 code (Kurucz, 1993) applied to the
corresponding atmospheric model calculated using ATLAS9 (Kurucz, 2005). The
influence of $\log g$ was neglected since neutral iron lines are insensitive
to this parameter. For the surface gravities, an ionization equilibrium was
imposed between Fe i@ and Fe ii@, where the line list for the latter was
provided by Romaniello et al. (2008). The atmospheric parameters were used to
compute stellar atmosphere models and synthetic spectra to analyze our
observed spectra (see Ripepi et al., 2021, for details). For the iron
abundance, we obtained a weighted average of the values derived from the
individual spectra available for each target. Fig. 1 displays random-phase
spectra for each of the 42 MW Cepheids observed with CFHT. The details
regarding the spectral observations, atmospheric and stellar parameters, and
abundance determinations will be presented in Catanzaro et al. (2023, in
prep.) together with new high-resolution near-infrared spectra for some of
these Cepheid standards obtained using the IGRINS instrument at the 8-m Gemini
South Telescope.

Figure 2: Comparison of our [Fe/H]${}_{\textrm{ESP}}$ values with measurements
from Luck (L18, 2018), Genovali et al. (G15, 2015), Ngeow (N12, 2012), and
Acharova et al. (A12, 2012). The mean offsets, standard deviation, and the
number of common stars are listed on the bottom right of each panel. The
dashed and dotted lines represent the mean and $\pm 2.5\sigma$ offsets,
respectively.

###  2.2 Comparison with literature metallicities

We compared our spectroscopic iron abundances ([Fe/H]${}_{\textrm{ESP}}$) with
those available in the literature from different compilations. Riess et al.
(2021) used spectroscopic metallicities for all 75 MW Cepheid standards from
the homogenized compilation by Groenewegen (2018), which were originally
either derived or homogenized by Genovali et al. (2014, 2015) in their
metallicity scale. We note that Genovali et al. (2014) obtained high-
resolution spectra for 42 Cepheids and combined their measurements with other
literature values to derive a common metallicity scale for a sample of 450
Cepheids. We find that 74 Cepheids have [Fe/H] values in Groenewegen (2018) or
Genovali et al. (2014, 2015), but CR Car is not included in these samples.
Riess et al. (2021) adopted a value of $-0.08$ dex for its iron abundance.
Recently, Kovtyukh et al. (2022) derived a [Fe/H] value of $0.14\pm 0.07$ dex
for CR Car, while Trentin et al. (2023) found their measurements
systematically overestimated by $+0.06$ dex. Therefore, we adopted a corrected
[Fe/H] value of $0.08\pm 0.08$ for CR Car. The updated sample222The [Fe/H]
values of 11 Cepheids (V0340 Ara, S Tra, XX Car, V0482 Sco, V0339 Cen, XZ Car,
DR Vel, SS CMa, AG Cru, XY Car, AQ Car) out of 75 in Riess et al. (2021) were
misattributed as small numbers ($<0.03$) due to an inversion between the value
and uncertainty columns in Groenewegen (2018). These were replaced with their
true [Fe/H] values from Genovali et al. (2014, 2015). The impact on $H_{0}$
was an increase by 0.02 km/s/Mpc. of Riess et al. (2021) will be discussed as
Genovali et al. (2015) sample for the subsequent analysis.

Fig. 2 displays a comparison of our [Fe/H]${}_{\textrm{ESP}}$ values for 42
stars with the compilations by Luck (2018), Genovali et al. (2015), Ngeow
(2012), and Acharova et al. (2012). The [Fe/H] values from these studies were
adopted from the compilation of Groenewegen (2018) except for Luck (2018)
measurements. The latter does not provide [Fe/H] values for RY CMa, and
therefore, only 41 stars were considered for comparison. Our [Fe/H]
measurements are in excellent agreement ($\Delta$[Fe/H]$=-0.02$ dex) with
those from Luck (2018) based on multi-phase spectra of Cepheid variables. The
largest difference is noted with the compilation by Ngeow (2012), which is
based on the metallicity scale of Luck & Lambert (2011). This was expected
since Luck & Lambert (2011) found a systematic offset of $+0.07$ dex in their
metallicity measurements with previous literature estimates by their group,
and a larger offset of $+0.11$ dex with independent measurements from
Romaniello et al. (2008). These [Fe/H] offsets were attributed to higher
surface gravity estimates ($\Delta\log g\sim+0.1$ leads to
$\Delta$[Fe/H]$\sim+0.05$ dex) in Luck & Lambert (2011), which likely resulted
from a subjective line editing of a smaller number of Fe ii@ lines during the
aforementioned iterative process.

Luck (2018) reanalysed multi-phase spectra of Cepheids including those in the
sample of Luck & Lambert (2011) and found a small median offset of
$\Delta\textrm{[Fe/H]}=-0.03$ dex with the metallicity scale of Genovali et
al. (2015). Our $\textrm{[Fe/H]}_{\textrm{ESP}}$ values, while being
consistent with Luck (2018), are also systematically smaller by
$\Delta$[Fe/H]$=-0.07$ dex than the measurements from Genovali et al. (2015).
The largest difference of $\Delta$[Fe/H]$=0.35$ dex is noted for GQ Ori.
However, Genovali et al. (2015) mentions that the S/N of spectra for GQ Ori
was not good enough for abundance measurements, which may be contributing to
its relatively larger metallicity estimate as compared to our values and those
in Ngeow (2012) or Acharova et al. (2012). The differences between [Fe/H]
values from different studies are systematically larger for metal-poor stars
which could arise due to small differences in effective temperatures. This is
possible since the number of spectral line pairs is often not enough for LDR
method to determine very accurate temperatures for these lower metallicity
stars (Trentin et al., 2023). Nevertheless, the median uncertainty on our
measurements is $0.09$ dex, and therefore, the offsets with literature
compilations shown in Fig. 2 are within $1\sigma$ of systematic uncertainties
($\sim 0.1$ dex) on spectroscopic metallicities in different
scales/compilations. We also note that the [Fe/H] values of our 42 Cepheids
from these literature compilations are not completely independent. Often the
same literature measurements were homogenized by adopting a fixed zero-point
in different studies, and therefore a detailed investigation into systematic
differences is beyond the scope of this paper.

Figure 3: Histograms of metallicities (in grey) for Cepheids from Genovali et
al. (2015) that were used in Riess et al. (2021). Our homogenized sample of
[Fe/H]${}_{\textrm{ESP}}$ values (in pink) for all 75 Cepheids are also
overplotted. The solid and dashed lines represent the Gaussian fits to our
[Fe/H] values and those from Genovali et al. (2015), respectively. ‘M’ is the
mean value of metallicity for the sample of $N$ stars.

Our homogeneous [Fe/H] values for 41 stars are in excellent agreement with
those provided by Luck (2018) for common Cepheids. Therefore, we adopted their
[Fe/H] values for the remaining 33 stars in our sample. The small mean offset
($-0.02$ dex) was applied to Luck (2018) stars that complement our sample to
construct a homogenized sample of [Fe/H] values for all 75 MW Cepheid
standards. Fig. 3 displays the histograms of corrected [Fe/H] values from
Riess et al. (2021) and our homogenized metallicities
([Fe/H]${}_{\textrm{ESP}}$). It is evident that our new [Fe/H] values for
these Cepheids are on average $0.05$ dex smaller than previously adopted
values from Genovali et al. (2015) sample. However, our average [Fe/H] values
are also in excellent agreement with the mean metallicity of Cepheids in the
local region ($\langle\textrm{[Fe/H]}\rangle\sim 0.05$ dex, Luck, 2018).
Recently, Groenewegen & Lub (2023) also showed that the flux-weighted
gravities from spectral energy distributions are in good agreement with those
determined from high-resolution spectra of Cepheids in Luck (2018) scale, but
do not agree with most other literature measurements.

Table 1: Multiband photometric mean magnitudes of MW Cepheid variables.

ID | $\log P$ | $\overline{\omega}$ | $G$ | $G_{BP}$ | $G_{RP}$ | $V$ | $I$ | $J$ | $H$ | $K_{s}$ | $W^{H}_{m}$ | $E_{BV}$ | [Fe/H] | Ref.f
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| days | mas | mag | dex |
AA Gem | 1.053 | 0.311 | 9.363 | 9.956 | 8.620 | 9.714 | 8.547c | 7.647 | 7.191 | 7.086 | 6.860 | 0.345 | $-0.150$ | ESP
|  | 0.019 | 0.074 | 0.082 | 0.075 | 0.012 | 0.030 | 0.010 | 0.010 | 0.010 | 0.023 | 0.036 | 0.070 |
AD Pup | 1.133 | 0.254 | 9.541 | 10.125 | 8.786 | 9.826 | 8.672 | 7.723d | 7.341 | 7.144 | 7.011 | 0.363 | $-0.130$ | ESP
|  | 0.018 | 0.053 | 0.121 | 0.088 | 0.014 | 0.012 | 0.030 | 0.030 | 0.030 | 0.024 | 0.020 | 0.080 |
AQ Car | 0.990 | 0.361 | 8.586 | 9.070 | 7.937 | 8.864 | 7.896 | 7.164d | 6.811 | 6.573 | 6.373 | 0.168 | $-0.052$ | L18
|  | 0.017 | 0.044 | 0.052 | 0.055 | 0.012 | 0.011 | 0.030 | 0.030 | 0.030 | 0.011 | 0.013 | 0.100 |
AQ Pup | 1.479 | 0.294 | 8.127 | 8.974 | 7.250 | 8.625 | 7.120 | 6.139 | 5.522 | 5.336 | 4.859 | 0.531 | $-0.170$ | ESP
|  | 0.025 | 0.104 | 0.103 | 0.102 | 0.015 | 0.012 | 0.012 | 0.011 | 0.011 | 0.016 | 0.017 | 0.080 |
BK Aur | 0.903 | 0.426 | 9.071 | 9.678 | 8.314 | 9.472 | 8.242c | 7.305 | 6.871 | 6.748 | 6.539 | 0.393 | $0.080$ | ESP
|  | 0.017 | 0.047 | 0.048 | 0.037 | 0.012 | 0.030 | 0.011 | 0.010 | 0.010 | 0.029 | 0.026 | 0.090 |

Note. — The star ID, periods ($\log P$), and the Wesenheit magnitudes
($W^{H}_{m}$) in HST filters are from Riess et al. (2021). The parallaxes
include zero-point corrections from Lindegren et al. (2021). The reddening
values were adopted from Groenewegen (2018). The Gaia photometric magnitudes
were taken from Ripepi et al. (2022) and $VIJHK_{s}$-band magnitudes were
adopted from Bhardwaj et al. (2015), unless specified.

aThe Gaia magnitudes were taken from Gaia Collaboration et al. (2023).

b/cThe $V/I$-band magnitudes were derived using Gaia to Johnson-Cousin
photometric system based on transformation from Pancino et al. (2022).

dNIR ($JHK_{s}$) magnitudes were taken from Groenewegen (2018).

eFor CR Car, NIR magnitudes, reddening, and metallicities were taken from
Trentin et al. (2023).

fThe source of spectroscopic metallicities: CFHT spectra (ESP) and
measurements from Luck (L18, 2018) after correcting for a mean offset of
$-0.02$ dex (see text for details).

(This table is available in its entirety in machine-readable form.)

###  2.3 Multiband photometric mean-magnitudes

We extracted Gaia photometry and astrometry for all 75 MW Cepheid standards.
The intensity-weighted mean magnitudes in $G$, $G_{BP}$, and $G_{RP}$ filters
for 71 Cepheids were taken from Ripepi et al. (2022) based on Gaia light
curves. For the remaining stars, mean magnitudes were taken from Gaia
Collaboration et al. (2023) source catalog. We also adopted Gaia parallaxes
including the zero-point offset corrections from Lindegren et al. (2021) and
the distances from Bailer-Jones et al. (2021). Optical (72/65 stars in $V/I$)
and near-infrared (44 stars in $JHK_{s}$) light curves for Cepheids were
adopted from the homogenized compilation of Bhardwaj et al. (2015). Intensity-
weighted mean magnitudes were derived from their Fourier-fitted light curves.
When the light curves were not available, $JHK_{s}$ magnitudes were taken from
Groenewegen (2018). For the remaining $VI$-band magnitudes, we utilized
photometric transformations provided by Pancino et al. (2022) between Gaia and
the Johnson-Kron-Cousins system. For Cepheids with available $V$ and $I$ band
light curves, we found a mean offset of $-0.015\pm 0.029$ mag and $0.034\pm
0.030$ mag, respectively, between their intensity-weighted mean magnitudes and
photometrically transformed magnitudes. An uncertainty of 0.03 mag was
assigned for all magnitudes for which light curves were not available. The HST
photometry ($W^{H}_{m}$ Wesenheit magnitude) for Cepheids was also adopted
from Riess et al. (2021). Table 1 lists photometric mean magnitudes and
spectroscopic metallicities for all 75 MW Cepheids.

To apply reddening corrections, the color-excess, $E(B-V)$, values were taken
from Groenewegen (2018). We adopted the Fitzpatrick (1999) reddening law
assuming $R_{V}=3.1$, and corrected for extinction using the total-to-
selective absorption ratios - $R_{\lambda}=2.783/1.777/0.812/0.508/0.349$ for
$G$, $I$, $J$, $H$, and $K_{s}$, respectively (Breuval et al., 2022). The
photometric mean magnitudes in the $VIJHK_{s}$ bands were also used to
construct-reddening free Wesenheit magnitudes: $W_{G}=G-1.9(BP-RP)$ from
Ripepi et al. (2019), and
$W_{V}=I-1.387(V-I),~{}W_{JK_{s}}=K-0.735(J-K_{s}),~{}W_{VK_{s}}=K_{s}-0.127(V-K_{s}),$
following Breuval et al. (2022).

##  3 Metallicity dependence of Period-Luminosity relations

###  3.1 Linear regression fits

The photometric properties, parallaxes, and spectroscopic metallicities listed
in Table 1 were used to derive PLZ and PWZ relations for Cepheids at multiple
wavelengths. The extinction-corrected absolute magnitude or the Wesenheit
magnitude ($M_{\lambda_{i}}$) of the $i^{\textrm{th}}$ Cepheid with period
($P_{i}$) and metallicity ([Fe/H]i) at a given wavelength ($\lambda$) is
defined as:

$\displaystyle M_{\lambda_{i}}=\alpha_{\lambda}+\beta_{\lambda}(\log
P_{i}-1)+\gamma_{\lambda}\textrm{[Fe/H]}_{i},$ (1)

Figure 4: Period-Luminosity-Metallicity (left) and Period-Wesenheit-
Metallicity (right) relations for 68 MW Cepheids at multiple bands. The dashed
lines represent the best-fitting linear regression.

where the absolute or Wesenheit magnitudes were obtained using the distances
provided by Bailer-Jones et al. (2021), which include the parallax zero-point
correction by Lindegren et al. (2021). These distances were also corrected for
a parallax over-correction of d$\overline{\omega}=-14~{}\mu$as (Riess et al.,
2021) using the approximation $dr=-r^{2}d\overline{\omega}$ adopted in Breuval
et al. (2022). The $\alpha$, $\beta$, and $\gamma$ represent the zero-point,
slope, and the metallicity coefficient, respectively.

Following Riess et al. (2021), we excluded six Cepheids - DL Cas, RW Cam, SY
Nor, AD Pup, U Aql, and SV Per, which exhibit values greater than 12.5 in the
Gaia astrometric quality flag (astrometric_gof_al). These Cepheids have poor
parallaxes and are all known as members of binary systems: DL Cas (Gieren et
al., 1994), RW Cam and SV Per (Evans, 1995; Riess et al., 2018b), SY Nor
(Kervella et al., 2019), AD Pup (Szabados et al., 2013), U Aql (Evans, 1992;
Gallenne et al., 2019). The faintest Cepheid in our sample, CY Aur, was also
excluded since its mean $G$-band magnitude falls near a sharp inflection point
of larger parallax correction (Riess et al., 2021). For this sample of 68
Cepheids, we fitted a linear regression in the form of equation (1). Given
that the parallax uncertainties for these 68 Cepheids are relatively small
($<12\%$) and they are nearby ($<7$ kpc), we used distances from Bailer-Jones
et al. (2021) instead of parallaxes for the initial analysis to investigate
PLZ/PWZ relations. We applied an iterative outlier clipping method which
removed a Cepheid with the largest residual in each iteration until all stars
were within $\pm 2.5\sigma$ scatter. We adopted the bootstrapping procedure to
create $10^{4}$ random realizations of the fits and determine errors.

Fig. 4 displays PLZ and PWZ relations for Cepheids at multiple wavelengths.
The figure shows that the seven longest-period Cepheids ($\log P>1.42$) in our
sample are brighter than the best-fitting linear regression. This was also
seen for the $H$-band PL relation in Cruz Reyes & Anderson (2023, Figure 9)
and $W^{H}_{m}$ PW relation in Riess et al. (2022b, Figure 3) for these
Cepheid standards. The mean metallicities of these long-period Cepheids is
$+0.07$ dex with a scatter of $0.17$ dex. These Cepheids have
$\textrm{[Fe/H]}>0$ except for AQ Pup ($\textrm{[Fe/H]}=-0.17\pm 0.08$ dex).
The lack of metal-poor stars at the long-period end in our sample can impact
the slopes of the derived PLZ/PWZ relations. Although the non-linearity of
Cepheid PL relations has been extensively investigated previously (e.g.,
Bhardwaj et al., 2016), we only fit a single-slope model due to lower
statistics at the long-period end. It is also evident that one of the Cepheids
RX Cam ($\log P=0.89$) is an outlier in all PWZ relations. Hence this star is
excluded from the subsequent analysis.

Table 2: Multiband MW Cepheid period-luminosity-metallicity and period-
Wesenheit-metallicity relations.

Band | $\alpha$ | $\beta$ | $\gamma$ | $\sigma$ | $N_{\textrm{f}}$ | $\chi^{2}_{\textrm{dof}}$ | $\alpha_{\textrm{LMC}}$ | $\beta_{\textrm{LMC}}$ | $\alpha_{\textrm{MW}}$ | $\gamma_{f}$
---|---|---|---|---|---|---|---|---|---|---
|  PLZ/PWZ relations in the form of eq. (1) |  ABL fits with fixed slope in the form of eq. (2)
$V$ | $-4.01\pm 0.04$ | $-2.26\pm 0.08$ | $-0.21\pm 0.20$ | 0.17 | 63 | 1.05 | $-3.19\pm 0.01$ | $-2.71\pm 0.03$ | $-3.29\pm 0.04$ | $-0.22\pm 0.07$
$G$ | $-4.31\pm 0.05$ | $-2.67\pm 0.09$ | $-0.14\pm 0.24$ | 0.19 | 65 | 0.95 | $-3.33\pm 0.01$ | $-2.82\pm 0.03$ | $-3.52\pm 0.05$ | $-0.40\pm 0.10$
$I$ | $-4.78\pm 0.04$ | $-2.57\pm 0.07$ | $-0.28\pm 0.17$ | 0.13 | 63 | 0.95 | $-3.85\pm 0.01$ | $-2.95\pm 0.02$ | $-3.96\pm 0.04$ | $-0.23\pm 0.07$
$J$ | $-5.36\pm 0.04$ | $-3.00\pm 0.06$ | $-0.30\pm 0.17$ | 0.14 | 65 | 1.02 | $-4.33\pm 0.01$ | $-3.10\pm 0.01$ | $-4.46\pm 0.04$ | $-0.28\pm 0.07$
$H$ | $-5.69\pm 0.03$ | $-3.27\pm 0.06$ | $-0.37\pm 0.18$ | 0.16 | 68 | 0.97 | $-4.68\pm 0.01$ | $-3.16\pm 0.01$ | $-4.80\pm 0.03$ | $-0.27\pm 0.08$
$K_{s}$ | $-5.80\pm 0.04$ | $-3.32\pm 0.06$ | $-0.43\pm 0.18$ | 0.15 | 67 | 1.02 | $-4.71\pm 0.01$ | $-3.22\pm 0.01$ | $-4.86\pm 0.04$ | $-0.33\pm 0.07$
$W_{VI}$ | $-5.81\pm 0.04$ | $-3.08\pm 0.06$ | $-0.36\pm 0.16$ | 0.13 | 65 | 1.03 | $-4.77\pm 0.01$ | $-3.29\pm 0.01$ | $-4.87\pm 0.03$ | $-0.21\pm 0.07$
$W_{G}$ | $-5.94\pm 0.03$ | $-3.21\pm 0.07$ | $-0.33\pm 0.16$ | 0.13 | 64 | 1.04 | $-4.79\pm 0.01$ | $-3.34\pm 0.01$ | $-4.98\pm 0.03$ | $-0.43\pm 0.06$
$W^{H}_{m}$ | $-5.92\pm 0.03$ | $-3.27\pm 0.05$ | $-0.24\pm 0.15$ | 0.12 | 64 | 1.04 | $-4.82\pm 0.02$ | $-3.31\pm 0.04$ | $-4.95\pm 0.03$ | $-0.31\pm 0.07$
$W_{VK_{s}}$ | $-6.02\pm 0.04$ | $-3.40\pm 0.06$ | $-0.46\pm 0.18$ | 0.14 | 65 | 0.99 | $-4.91\pm 0.01$ | $-3.25\pm 0.01$ | $-5.07\pm 0.04$ | $-0.37\pm 0.08$
$W_{JK_{s}}$ | $-6.12\pm 0.04$ | $-3.46\pm 0.06$ | $-0.45\pm 0.19$ | 0.15 | 66 | 0.99 | $-4.99\pm 0.01$ | $-3.32\pm 0.01$ | $-5.16\pm 0.04$ | $-0.37\pm 0.08$

Note. — The zero-point ($\alpha$) at $\log P=1.0$ days, slope ($\beta$),
metallicity coefficient ($\gamma$), dispersion ($\sigma$), and the number of
stars ($N_{f}$) of the final PLZ/PWZ fits of eq. (1) are tabulated in the left
columns. The coefficients of the ABL fits
($\alpha_{\textrm{LMC}},~{}\beta_{MW}$) with a fixed slope for Cepheids in the
LMC are shown in the columns on the right. The metallicity coefficient
($\gamma_{f}$) was determined by comparing the intercepts (at $\log P=0.7$
days) of the PL relations for Cepheids in the MW and LMC as a function of
their mean-metallicities.

Table 2 lists the coefficients of the best-fitting linear regression for each
PLZ/PWZ relation. The slope and the zero-points of these relations are well-
constrained, but the metallicity coefficients exhibit large uncertainties.
Note that the metallicity values in our sample span $-0.21$ to $0.39$ dex,
covering a total range of $\Delta$[Fe/H]$=0.6$ dex. However, the majority of
Cepheids have [Fe/H] values between $-0.1$ and $0.2$ dex as seen in Fig. 3.
Given that the median uncertainty on these spectroscopic measurements is of
the order of $0.1$ dex, this smaller range of metallicities limits the
accurate and precise determination of the metallicity coefficient of the
PLZ/PWZ relations. Nevertheless, the metallicity coefficients at all
wavelengths exhibit a negative sign implying that the metal-rich Cepheids are
brighter. The metallicity coefficients vary between $-0.14$ to $-0.46$
mag/dex, a range that is in agreement with recent empirical studies in the
literature (e.g., Ripepi et al., 2021; Breuval et al., 2022). Our $\gamma$
values for $K_{s}$, $W_{JK_{s}}$, and $W_{JK_{s}}$ PLZ/PWZ relations for MW
Cepheids are in better agreement with the value of $\sim-0.45$ mag/dex
obtained by Ripepi et al. (2021) than that of $\sim-0.32$ mag/dex in Breuval
et al. (2022). This is not surprising since our approach of using individual
metallicities of MW Cepheids is similar to the work of Ripepi et al. (2021).
The negative sign of the metallicity term of PWZ relations has also been
predicted by stellar evolution (Anderson et al., 2016) and pulsation models
(De Somma et al., 2022).

###  3.2 Astrometry-Based Luminosity fits and additional metallicity
constraints with Cepheids in the LMC

The varying precision of Cepheid parallaxes can bias absolute magnitude
determination since the direction-dependent priors constructed from a three-
dimensional model of the Galaxy in Bailer-Jones et al. (2021) can dominate the
distance estimates for faint and distant stars. Therefore, it is more
appropriate to derive the coefficients of the PL relations in parallax space,
which can be done using astrometry-based luminosities (ABL, Arenou & Luri,
1999). Our MW Cepheid sample is small and does not provide strong constraints
on the metallicity coefficient as seen in Table 2. We decided to include
additional constraints from the Cepheids in the LMC, and define ABL as:

ABL $\displaystyle=$ $\displaystyle~{}\overline{\omega}10^{0.2m_{\lambda}-2}$
$\displaystyle=$
$\displaystyle~{}10^{0.2(\alpha_{\lambda,\textrm{MW}}+\beta_{\lambda,\textrm{LMC}}(\log
P_{i}-0.7))},$ (2)

where $m_{\lambda}$ is the extinction-corrected magnitude or the Wesenheit
magnitude at a given wavelength for a Cepheid with period ($P_{i}$). The fixed
LMC slope ($\beta_{\textrm{LMC}}$) was adopted from Breuval et al. (2022) to
constrain the zero-point of the MW Cepheid PL relations
($\alpha_{\textrm{MW}}$). The absolute zero-point was also shifted to $\log
P=0.7$ to compare with $\alpha_{\textrm{LMC}}$ at the same period from Breuval
et al. (2022). The mean metallicity of our homogenized sample is
$\langle\textrm{[Fe/H]}_{\textrm{ESP}}\rangle=+0.04$ dex with a dispersion of
0.09 dex. A mean metallicity of $-0.41\pm 0.02$ was adopted for the LMC
Cepheids (Romaniello et al., 2022). The metallicity coefficient was then
calculated using
$\alpha=\gamma_{f}(\textrm{[Fe/H]}-\langle\textrm{[Fe/H]}_{\textrm{ESP}}\rangle)+\delta$.
Since we have only two intercepts (MW and LMC), measuring the zero-point at
the mean-metallicity of MW Cepheids ensures that the $\delta$ is the same as
$\alpha_{MW}$ in Table 2.

The last column of Table 2 lists the metallicity coefficients of PLZ/PWZ
relations obtained using this approach. The range of $\gamma_{f}$ is similar
to the one obtained using linear regression. However, these coefficients are
determined up to three times more precisely thanks to the additional
constraint on the slope using the Cepheids in the LMC. The metallicity
coefficients of the PLZ/PWZ relations vary between $-0.21$ dex and $-0.43$
mag/dex for different wavelengths and are in good agreement with those from
Breuval et al. (2022). The adopted reddening values for our sample of Cepheids
are different than those adopted in Breuval et al. (2022), which can also
impact the determined zero-points and the metallicity coefficients of PL
relations. The $\gamma_{f}$ values for the most PWZ relations are now smaller
than those from Ripepi et al. (2021) but within $1\sigma$ errors. For this
small sample of Cepheids, it can be noted that the metallicity coefficient
varies depending on the adopted methodology of using individual abundances of
MW Cepheids only or using Galactic and Magellanic Cloud Cepheids together
assuming their mean-metallicities for the quantification of metallicity
coefficient, but appear statistically consistent.

###  3.3 Implications for the Hubble constant

Our new spectroscopic metallicities for MW Cepheid standards were also used to
derive the $W^{H}_{m}$ PWZ relation using their HST photometry from Riess et
al. (2021). From Table 2, the zero-point, slope, and the metallicity
coefficient of $W^{H}_{m}$ PWZ relation are: $-5.92\pm 0.03$ mag, $-3.27\pm
0.05$ mag/dex, and $-0.24\pm 0.15$ mag/dex, which agree within $0.5\sigma$
with the results of free-parameter fits in Riess et al. (2021, their
4-parameter solution in Table 2). The metallicity coefficients differ by
$0.04$ mag/dex, but the uncertainty on this parameter is largely due to low
statistics and limited metallicity range. This $\gamma$ parameter is smaller
(in absolute sense) than $-0.37\pm 0.09$ mag/dex determined by Ripepi et al.
(2021) using a larger sample of MW Cepheids. Comparing the $\gamma$ value
obtained using MW and LMC intercepts for $W^{H}_{m}$ PWZ relation, our value
of $-0.31\pm 0.07$ mag/dex is also in agreement with a value of $-0.28\pm
0.08$ mag/dex determined by Breuval et al. (2022) using a similar approach.

To investigate possible impact of our new spectroscopic metallicities on the
local $H_{0}$ determination, we replaced our [Fe/H] values with those adopted
in the SH0ES project (Riess et al., 2022a). The additional constraints on the
coefficients of the PL relation come from the two independent anchors - the
LMC and NGC 4258. We consider only two anchor variants of $H_{0}$ for our
purpose: (1) MW + LMC + NGC 4258; and (2) MW. These two calibrations will
allow us to explore possible variations in the metallicity coefficient and
corresponding systematics due to differences in the mean metallicities of the
calibrators and the target SNe host galaxies.

Adopting all three anchors and using our [Fe/H] values in the global-matrix
formalism of the SH0ES project (see Section 2.1 of Riess et al., 2022a, b), we
obtained a metallicity coefficient of $\gamma=-0.23\pm 0.04$ mag/dex and an
$H_{0}=72.9\pm 1.0$ km s-1 Mpc-1. The difference in $H_{0}$ is 0.1 km/s/Mpc as
compared to Riess et al. (2022b). In the case of only the MW as an anchor,
these parameters are $\gamma=-0.21\pm 0.05$ mag/dex and an $H_{0}=72.7\pm 1.3$
km s-1 Mpc-1, a difference in $H_{0}$ of 0.2 km/s/Mpc as compared to Riess et
al. (2022b). This suggests that the systematic variations in the metallicities
of Galactic Cepheid standards due to adopted metallicity scale produce changes
in $H_{0}$ that are much smaller than the statistical uncertainty in the local
$H_{0}$ determinations (Riess et al., 2022a).

##  4 Summary

The Milky Way Cepheids photometrically observed with HST play a crucial role
in calibrating the anchor PL relation in the same photometric system as the
target Cepheids in Supernovae host galaxies. We presented new high-resolution,
high signal-to-noise spectra for 42 of these 75 MW Cepheid standards using
ESPaDOnS instrument at CFHT. A homogenized sample of spectroscopic [Fe/H]
values for all 75 Cepheids was obtained for the calibration of Galactic
Cepheid PLZ and PWZ relations for the distance scale applications. Our
spectroscopic metallicities are in good agreement with most literature
measurements in different metallicity scales, but are systematically smaller
by $0.05$ dex than previously adopted metallicities for these 75 Cepheids in
the local determination of $H_{0}$. The mean metallicity of our sample
${\langle\textrm{[Fe/H]}\rangle}=0.04$ dex agrees well with the mean
metallicity of Cepheids in the local region (Luck, 2018).

The Gaia parallaxes together with the literature photometry and new
spectroscopic metallicities were used to derive PLZ and PWZ relations for MW
Cepheids at multiple wavelengths. The low statistics and a relatively limited
range of metallicity ($\Delta\textrm{[Fe/H]}=0.6$ dex) impacted the precision
of the metallicity coefficients of the PLZ/PWZ relations. Adding Cepheids in
the LMC to provide additional constraints, we improve the accuracy and the
precision of metallicity coefficients of the PLZ and PWZ relations. These
metallicity coefficients are in good agreement with recent determinations in
the literature. Moreover, we explored the differences in the metallicity
coefficient of a given PLZ/PWZ relation obtained by fitting individual
metallicities of MW Cepheids or comparing the intercepts of the PL/PW
relations of MW and LMC Cepheids as a function of their mean metallicities.
For most cases, comparing the intercepts of PL/PW relations for Cepheids in
the MW and LMC gives a smaller metallicity coefficient presumably due to a
narrow metallicity range ($\Delta\textrm{[Fe/H]}=0.46$ dex). Since precise
geometric distances to nearby galaxies are limited to the Magellanic Clouds,
it is imperative to increase the metallicity range of MW Cepheids to better
constrain their PLZ/PWZ relations (e.g. in Ripepi et al., 2021; Trentin et
al., 2023).

However, the metallicity coefficient of $W^{H}_{m}$ Wesenheit relation in HST
filters does not show any statistically significant variation with our new
spectroscopic metallicities. Riess et al. (2022a) showed that the $H_{0}$
values do not exhibit any correlation with metallicities because the average
abundances of Cepheids in the anchor and target SNe host galaxies span the
same range. Our new average [Fe/H] value ($+0.04$ dex) for MW Cepheids is
smaller than the previously adopted mean metallicity ($+0.09$ dex). For a
metallicity term of $\sim-0.2$ mag/dex, this would translate into a 0.01 mag
difference in distance moduli measurements when the MW is the sole anchor, or
less than half of this for the combination of all three anchors. Incorporating
these new metallicity measurements in the latest SH0ES distance ladder
formalism (Riess et al., 2022a), a small difference in $H_{0}$ of $0.2$
km/s/Mpc was found for the MW as the only anchor case. If all three anchors
are used as calibrators, which is the baseline solution for the $H_{0}$
determination, the new metallicities result in a little variation of $0.1$
km/s/Mpc in the $H_{0}$ value. Thus MW Cepheid metallicities in different
metallicity scales make a small difference in $H_{0}$ relative to the size of
the Hubble Tension.

We thank the anonymous referee for useful comments that helped improve the
manuscript. This project has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No. 886298. This research was supported by the Munich
Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany´s Excellence Strategy – EXC-2094 – 390783311. CCN thanks the funding
from the National Science and Technology Council (Taiwan) under the contract
109-2112-M-008-014-MY3. Access to the CFHT was made possible by the Institute
of Astronomy and Astrophysics, Academia Sinica. CFHT-ESPaDOnS

## References

  * Acharova et al. (2012) Acharova, I. A., Mishurov, Y. N., & Kovtyukh, V. V. 2012, MNRAS, 420, 1590, doi: 10.1111/j.1365-2966.2011.20161.x
  * Anderson et al. (2016) Anderson, R. I., Saio, H., Ekström, S., Georgy, C., & Meynet, G. 2016, A&A, 591, A8, doi: 10.1051/0004-6361/201528031
  * Arenou & Luri (1999) Arenou, F., & Luri, X. 1999, in Astronomical Society of the Pacific Conference Series, Vol. 167, Harmonizing Cosmic Distance Scales in a Post-HIPPARCOS Era, ed. D. Egret & A. Heck, 13–32, doi: 10.48550/arXiv.astro-ph/9812094
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147, doi: 10.3847/1538-3881/abd806
  * Bhardwaj et al. (2016) Bhardwaj, A., Kanbur, S. M., Macri, L. M., et al. 2016, MNRAS, 457, 1644, doi: 10.1093/mnras/stw040
  * Bhardwaj et al. (2015) Bhardwaj, A., Kanbur, S. M., Singh, H. P., Macri, L. M., & Ngeow, C.-C. 2015, MNRAS, 447, 3342, doi: 10.1093/mnras/stu2678
  * Breuval et al. (2022) Breuval, L., Riess, A. G., Kervella, P., Anderson, R. I., & Romaniello, M. 2022, ApJ, 939, 89, doi: 10.3847/1538-4357/ac97e2
  * Cruz Reyes & Anderson (2023) Cruz Reyes, M., & Anderson, R. I. 2023, A&A, 672, A85, doi: 10.1051/0004-6361/202244775
  * da Silva et al. (2022) da Silva, R., Crestani, J., Bono, G., et al. 2022, A&A, 661, A104, doi: 10.1051/0004-6361/202142957
  * De Somma et al. (2022) De Somma, G., Marconi, M., Molinaro, R., et al. 2022, ApJS, 262, 25, doi: 10.3847/1538-4365/ac7f3b
  * Di Valentino et al. (2021) Di Valentino, E., Mena, O., Pan, S., et al. 2021, Classical and Quantum Gravity, 38, 153001, doi: 10.1088/1361-6382/ac086d
  * Donati et al. (1997) Donati, J. F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658, doi: 10.1093/mnras/291.4.658
  * Evans (1992) Evans, N. R. 1992, ApJ, 389, 657, doi: 10.1086/171238
  * Evans (1995) —. 1995, ApJ, 445, 393, doi: 10.1086/175704
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63, doi: 10.1086/316293
  * Freedman et al. (2001) Freedman, W. L., Madore, B. F., Gibson, B. K., et al. 2001, ApJ, 553, 47, doi: 10.1086/320638
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940
  * Gallenne et al. (2019) Gallenne, A., Kervella, P., Borgniet, S., et al. 2019, A&A, 622, A164, doi: 10.1051/0004-6361/201834614
  * Genovali et al. (2014) Genovali, K., Lemasle, B., Bono, G., et al. 2014, A&A, 566, A37, doi: 10.1051/0004-6361/201323198
  * Genovali et al. (2015) Genovali, K., Lemasle, B., da Silva, R., et al. 2015, A&A, 580, A17, doi: 10.1051/0004-6361/201525894
  * Gieren et al. (1994) Gieren, W. P., Welch, D. L., Mermilliod, J.-C., Matthews, J. M., & Hertling, G. 1994, AJ, 107, 2093, doi: 10.1086/117019
  * Groenewegen & Lub (2023) Groenewegen, M., & Lub, J. 2023, arXiv e-prints, arXiv:2307.07559, doi: 10.48550/arXiv.2307.07559
  * Groenewegen (2018) Groenewegen, M. A. T. 2018, A&A, 619, A8, doi: 10.1051/0004-6361/201833478
  * Kervella et al. (2019) Kervella, P., Gallenne, A., Evans, N. R., et al. 2019, A&A, 623, A117, doi: 10.1051/0004-6361/201834211
  * Kovtyukh et al. (2022) Kovtyukh, V., Lemasle, B., Bono, G., et al. 2022, MNRAS, 510, 1894, doi: 10.1093/mnras/stab3530
  * Kovtyukh & Gorlova (2000) Kovtyukh, V. V., & Gorlova, N. I. 2000, A&A, 358, 587
  * Kurucz (1993) Kurucz, R. L. 1993, SYNTHE spectrum synthesis programs and line data
  * Kurucz (2005) —. 2005, Memorie della Societa Astronomica Italiana Supplementi, 8, 14
  * Landsman (1993) Landsman, W. B. 1993, in ASP Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 246
  * Leavitt & Pickering (1912) Leavitt, H. S., & Pickering, E. C. 1912, Harvard College Observatory Circular, 173, 1
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2, doi: 10.1051/0004-6361/202039709
  * Luck (2018) Luck, R. E. 2018, AJ, 156, 171, doi: 10.3847/1538-3881/aadcac
  * Luck & Lambert (2011) Luck, R. E., & Lambert, D. L. 2011, AJ, 142, 136, doi: 10.1088/0004-6256/142/4/136
  * Molinaro et al. (2023) Molinaro, R., Ripepi, V., Marconi, M., et al. 2023, MNRAS, 520, 4154, doi: 10.1093/mnras/stad440
  * Ngeow (2012) Ngeow, C.-C. 2012, ApJ, 747, 50, doi: 10.1088/0004-637X/747/1/50
  * Owens et al. (2022) Owens, K. A., Freedman, W. L., Madore, B. F., & Lee, A. J. 2022, ApJ, 927, 8, doi: 10.3847/1538-4357/ac479e
  * Pancino et al. (2022) Pancino, E., Marrese, P. M., Marinoni, S., et al. 2022, A&A, 664, A109, doi: 10.1051/0004-6361/202243939
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Riess et al. (2021) Riess, A. G., Casertano, S., Yuan, W., et al. 2021, ApJL, 908, L6, doi: 10.3847/2041-8213/abdbaf
  * Riess et al. (2011) Riess, A. G., Macri, L., Casertano, S., et al. 2011, ApJ, 730, 119, doi: 10.1088/0004-637X/730/2/119
  * Riess et al. (2016) Riess, A. G., Macri, L. M., Hoffmann, S. L., et al. 2016, ApJ, 826, 56, doi: 10.3847/0004-637X/826/1/56
  * Riess et al. (2018a) Riess, A. G., Casertano, S., Yuan, W., et al. 2018a, ApJ, 855, 136, doi: 10.3847/1538-4357/aaadb7
  * Riess et al. (2018b) —. 2018b, ApJ, 861, 126, doi: 10.3847/1538-4357/aac82e
  * Riess et al. (2022a) Riess, A. G., Yuan, W., Macri, L. M., et al. 2022a, ApJL, 934, L7, doi: 10.3847/2041-8213/ac5c5b
  * Riess et al. (2022b) Riess, A. G., Breuval, L., Yuan, W., et al. 2022b, ApJ, 938, 36, doi: 10.3847/1538-4357/ac8f24
  * Ripepi et al. (2019) Ripepi, V., Molinaro, R., Musella, I., et al. 2019, A&A, 625, A14, doi: 10.1051/0004-6361/201834506
  * Ripepi et al. (2021) Ripepi, V., Catanzaro, G., Molinaro, R., et al. 2021, MNRAS, 508, 4047, doi: 10.1093/mnras/stab2460
  * Ripepi et al. (2022) Ripepi, V., Clementini, G., Molinaro, R., et al. 2022, arXiv e-prints, arXiv:2206.06212, doi: 10.48550/arXiv.2206.06212
  * Romaniello et al. (2008) Romaniello, M., Primas, F., Mottini, M., et al. 2008, A&A, 488, 731, doi: 10.1051/0004-6361:20065661
  * Romaniello et al. (2022) Romaniello, M., Riess, A., Mancino, S., et al. 2022, A&A, 662, C1, doi: 10.1051/0004-6361/202142441e
  * Szabados et al. (2013) Szabados, L., Derekas, A., Kiss, L. L., et al. 2013, MNRAS, 430, 2018, doi: 10.1093/mnras/stt027
  * Trentin et al. (2023) Trentin, E., Ripepi, V., Catanzaro, G., et al. 2023, MNRAS, 519, 2331, doi: 10.1093/mnras/stac2459
  * Verde et al. (2019) Verde, L., Treu, T., & Riess, A. G. 2019, Nature Astronomy, 3, 891, doi: 10.1038/s41550-019-0902-0


adieresis=ä,germandbls=ß

# In-Network Processing Acoustic Data for Anomaly Detection in Smart Factory

Huanzhuo Wu1, Yunbin Shen1, Xun Xiao3, Artur Hecker3 and Frank H.P. Fitzek12
<EMAIL_ADDRESS><EMAIL_ADDRESS>dresden.de<EMAIL_ADDRESS>This is a preprint of the
work [1], that has been accepted for publication in the proceedings of the
2021 IEEE Global Communications Conference. 1Deutsche Telekom Chair of
Communication Networks - Technische Universität Dresden, Germany  2Centre for
Tactile Internet with Human-in-the-Loop (CeTI)  3Munich Research Center,
Huawei Technologies, Muenchen, Germany

###### Abstract

Modern manufacturing is now deeply integrating new technologies such as 5G,
Internet-of-things (IoT), and cloud/edge computing to shape manufacturing to a
new level – Smart Factory. Autonomic anomaly detection (e.g., malfunctioning
machines and hazard situations) in a factory hall is on the list and expects
to be realized with massive IoT sensor deployments. In this paper, we consider
acoustic data-based anomaly detection, which is widely used in factories
because sound information reflects richer internal states while videos cannot;
besides, the capital investment of an audio system is more economically
friendly. However, a unique challenge of using audio data is that sounds are
mixed when collecting thus source data separation is inevitable. A traditional
way transfers audio data all to a centralized point for separation.
Nevertheless, such a centralized manner (i.e., data transferring and then
analyzing) may delay prompt reactions to critical anomalies. We demonstrate
that this job can be transformed into an in-network processing scheme and thus
further accelerated. Specifically, we propose a progressive processing scheme
where data separation jobs are distributed as microservices on intermediate
nodes in parallel with data forwarding. Therefore, collected audio data can be
separated $43.75\%$ faster with even less total computing resources. This
solution is comprehensively evaluated with numerical simulations, compared
with benchmark solutions, and results justify its advantages.

###### Index Terms:

In-network Computing, IoT, Smart Factory, Anomaly Detection

##  I Introduction

The contemporary industry is in the dawn of the 4th revolution towards full
digitization and intelligence. Deep integration of emerging technologies such
as 5G, Internet-of-things (IoT), artificial intelligence (AI), and cloud
computing [2] is happening right now and brings manufacturing to a new level –
Smart Factory. A critical operational task in manufacturing is anomaly
detection for machines malfunctioning on the production line. It prevents
machines and/or products from any serious damage and economic losses [3]. In
the past, this detection was manually done by on-site inspectors, who were
replaced by remote inspectors monitoring on screens. In the future, such
inspection would be fully autonomous by analyzing data from IoT sensors and
necessary reactions will be triggered without human intervention.

For anomaly detection, acoustic data can reflect the internal states of
machines that are not visible through videos [4]. For a long time, acoustic
data are used by experienced workers, who can directly judge whether or not a
machine works properly by hearing. In addition, audio devices are much cheaper
than professional cameras thus more friendly to capital investments.
Therefore, acoustic data-based anomaly detection will still play an important
role in future smart factory operations.

Video data is naturally separated when captured with cameras, but this is not
the case for acoustic data because sounds interfere with each other by nature.
Thus, anomaly detection based on acoustic data is more challenging as original
signal data have to be restored first. A natural idea is to first transfer all
data to a centralized node; when all data are received, a sort of Blind Source
Separation (BSS) [5] algorithm is applied to separate mixed data. BSS
candidates include Independent Component Analysis (ICA)-based methods [6, 7,
8] or neural network-based methods [9, 10]. However, forwarding and then
analyzing could delay critical decision-making actions due to $i$) possibly
long waiting time of transferring the data, and $ii$) possibly long execution
time of running the algorithm on a single node. Clearly, the realization of
autonomous anomaly detection requires a better solution.

In this paper, we tackle this problem from another angle: instead of sending
all data and then waiting for the separation result, we are thinking if the
whole task can be accelerated by distributing the data separation job on
intermediate forwarding nodes. In other words, we try to transform the
centralized manner into an in-network processing manner to speed up the entire
job. The key idea is sketched as follows: a new _lightweight processing logic_
is proposed and deployed on every intermediate node running with local spare
compute resources as microservices; every node best-effort computes a temporal
result, which is a _solution matrix_ that is roughly estimated to restore the
original data; this temporal result (i.e., the solution matrix) will be
progressively optimized along the forwarding path. Specifically, a modified
ICA algorithm is proposed so that progressive improvements on every node can
be maximized in order to fit such a distributed processing scheme; by doing
so, at the final destination (i.e., the last hop), an optimal solution matrix
shall be ready with sufficient precision approximately. In summary, our key
contributions can be summarized as follows:

  1. 1.

We propose an in-network processing solution for acoustic data-based anomaly
detection, which is demonstrated as an example how audio data separation can
be accelerated up to $43.75\%$ by utilizing intermediate computing resources
in networks;

  2. 2.

We design a specific processing logic for intermediate nodes with a modified
ICA algorithm, making the acoustic data restoration to be distributedly
executable as microservices and yield a faster convergence speed;

  3. 3.

We conduct comprehensive simulations and numerical results which justify the
effectiveness of our proposed scheme.

To the best of our knowledge, technically, this is the first work that studies
how to transform a BSS algorithm into an in-network processing scheme,
overcoming a key constraint where traditional BSS can be mainly executed on a
centralized node.

The rest of the paper is organized as follows. In Section II, a literature
review is provided and the main differences of our solution are highlighted;
in Section III, we present full details of our solution. After that, in
Section IV, comprehensive numerical results will be presented, and Section V
concludes this paper.

##  II Related Work

In-network computing/processing – a joint consideration of computing and
communication – raises increasing research interests because network elements
now have much more powerful computing capabilities. In-network processing
empowered with microservices provides a new possibility and flexibility to
support emerging applications that require low latency, high bandwidth, and/or
high reliability [11].

Existing in-network processing studies focus more on how to embed computing
tasks into a network, deriving the processing results closer to the users
[12]. For example, the work in [13] decomposes an image detection job as a
service function chain. Since the processing can be simply done at closer
nodes, it reduces the latency by more than $25\%$. Similarly, in [14] a
computer vision task, a real-time line-following service, is decomposed and
deployed in a network with a series of programmable network devices. These
works consider only the processing location and network transport issues. On
the other hand, the processing logic does not have to be changed or modified.

In contrast, not every task can be deployed in a network straightforwardly. An
exception is running a BSS algorithm for mixture data separation. This is
because the algorithm cannot be simply split into sub-tasks, running on
multiple nodes in parallel or a chain. Differently, a joint task is solved by
distributed nodes coordinating with each other. Therefore, our problem
considers how a BSS algorithm can be run on distributed nodes. This goal is
clearly more sophisticated than pure task decomposition.

When it comes to a BSS problem, many candidate options are available. One
school is machine learning (ML) based on neural networks (NNs), such as Y-Net
[9], Conv-TasNet [10]. However, for our problem, ML-based solutions are less
interesting because $i$) it is hard to obtain enough representative and
labeled training data, $ii$) training an NN model is time-consuming and
resource-intensive, and $iii$) once deployed on nodes, NN models are
inflexible to be updated. Additionally, to maximize their performance, ML-
based solutions require special hardware (e.g., GPUs), which barely appears on
a network device. Another school is ICA algorithms, working directly with
input data and require only average hardware capability. Typical ICA
algorithms are FastICA [6], InfoMax [7], and CdICA [8]. They are free from the
constraints of ML-based solutions, so more feasible to fit an in-network
processing scheme.

However, existing work only provides centralized ICA algorithms, which cannot
be trivially transplanted to an in-network processing framework. The main
reason is: they require all data to calculate a final result (i.e., the
solution matrix mentioned before). Therefore, simply executing an ICA
algorithm on every node equals repeating the same procedure by multiple nodes.
This does not improve the ultimate convergence.

Some recent works realized this issue and looked for a distributed version. We
note that Adaptive Extraction-based ICA (AeICA) [15] has the potential to be
performed in networks but is quite sensitive to initial parameter
configurations, which heavily depend on the prior knowledge of input data.
According to reported results, the processing time would rather increase in
some cases. Its heterogeneous performances fail to fulfill our goal. We thus
only pick it as one of our benchmarks.

In summary, the main differences of this work are that: we study a non-trivial
in-network processing problem where the task – data separation – cannot be
simply decomposed to multiple equivalent sub-tasks; secondly, we fill the gap
by converting a traditional ICA algorithm to fit in a distributed environment
and overcome the deficiency.

##  III Our Solution

###  III-A Main Idea

According to our observations, the performance of a traditional ICA algorithm
is mainly constrained by the following two factors. The first one is the time
spent on transferring all data. Before all data are received, the algorithm
suspends, with unnecessary waiting time even with perfect network conditions.
The second one is the execution time to wait for a full convergence when
running on a centralized node. As we will see in Section III-C1, the main
procedure of an ICA algorithm is to gradient-descent optimize a solution
matrix in iterations. Marginal gains usually become much smaller when
approaching the convergence point. In other words, more time is spent at later
stages but the improvement is little.

Motivated by these observations, our main idea is to transform this sequential
processing manner into an in-network processing manner so that data processing
jobs can start early on intermediate nodes in parallel with data forwarding.
The achieved acceleration depends on the availability of spare compute
resources along the forwarding path. Obviously, the more intermediate compute
resources available are, the more accelerations intermediate nodes can
contribute. With this idea, new strategies are introduced as our design
principle.

We first introduce a _greedy strategy_ , wherein once a node finds that its
improvement from gradient-descent gets slower, it stops its local iterations
and hands over the temporal result to the next hop. However, simply relaying
the temporal results does not bring any acceleration. This leads to the next
strategy.

We further introduce a _growing strategy_ on the size of the input dataset to
make every node only cache a subset data from the whole dataset but the size
of the subset data _progressively_ increases on the following nodes. With more
input data, this thus guarantees that later nodes can better improve the
result.

Together with the greed strategy, the synthetic effect is that every node
takes a temporal solution from its last hop as an input (i.e., a better
starting point) and consumes more data to derive an even better solution for
the following nodes. This can mitigate the problem of spending too much
processing time for little marginal improvements at later iteration stages
when running in a centralized manner. Clearly, all these new strategies can be
easily realized as microservices on distributed nodes.

Before we start to introduce technical details, basic assumptions are
clarified here:

  * •

We focus on the algorithmic part and its distributed design in this work. Thus
we assume that possible node failures are handled by mechanisms at the network
resource layer (i.e., a secondary concern here);

  * •

We assume that the network resource allocation happens at earlier stages by
management layer; our starting point begins with considering the actual
service directly (i.e., in-network processing);

  * •

We assume that an intermediate node is not as powerful as a Commercial off-
the-shelf (COTS) server machine and a single network device cannot afford the
whole processing task. However, collectively, the sum of all distributed
compute resources is enough.

###  III-B Scenario and System Architecture

Figure 1: An in-network processing architecture.

An example scenario considered here is illustrated in Fig. 1. Specifically,
every IoT sensor collects sound data from the machine “source” it attaches to.
One audio sound signal interferes in the air with sounds from other machines
(as noise). Every IoT sensor constantly sends collected data to a wireless
access point (AP). The AP node can be either a WiFi AP or a cellular base
station of a (non-public) network. This AP keeps forwarding the data (from all
sensors) to the backend “Remote Agent” over a forwarding path. This path can
be either dynamically determined based on a routing protocol or statically
configured. The forwarding path consists of some intermediate “Node $k$” and
“Node $k+1$”.

Given this system, there are $n$ working machines indexed by $i$. The original
data denoted by $s_{i}$ is generated from the $i$-th working machine. As said,
the original signal data $s_{i}$ will be distorted to imperfect data $x_{i}$
due to mutual interference. With $m$ time slots, collected data at the AP
aggregate to a data matrix $X$. Accordingly, the original counterpart of $X$
is thus a source matrix $S$. We formulate this mutual interference effect as
the operation in Eq. 1:

$X=A\times S$ (1)

where distorting the original source data $S$ to $X$ is modeled by a mixing
matrix $A$. Mathematically, a BSS problem is an inverse operation of Eq. 1:
$\hat{S}=A^{-1}\times X=W\times X$, where original data $\hat{S}$ are
estimated by deriving the solution matrix $W$, which will be applied on the
input data $X$ for data restoration.

Our in-network processing scheme introduces a new _processing logic_ as a
microservice running on every intermediate node (gear icons in Fig. 1). As
introduced in the previous section, while forwarding the collected data $X$,
node $k$ will start processing with a subset data sampled from $X$. This aims
to roughly but quickly calculate a temporal solution matrix $W_{k}$, then
passes $W_{k}$ to node $k+1$ (recall the greedy strategy). Node $k+1$ will
evolve the temporal solution to a better solution matrix $W_{k+1}$ with a
larger size of subset data (recall the growing strategy). We will see that
such a progressive manner yields a faster converging speed without sacrificing
the precision on the final solution matrix $W$.

Next, we detail the new processing logic, namely progressive ICA (pICA), which
functions as the key component calculating the solution matrix $W$. Besides,
all notations are summarized in Table I.

TABLE I: Notation Notation |  Meaning
---|---
$s_{i}$ |  Original source data from the $i$-th machine.
$S$ |  Original data matrix from all machines.
$x_{i}$ |  Data collected by the $i$-th IoT sensor.
$X$ |  Data matrix received by AP.
$m$ |  Total time span.
$A$ |  Mixing matrix mimic mutual interference.
$W$ |  Solution matrix separating $S$ out of $X$.
$W_{k}$ |  Solution matrix on the $k$-th network node.
$\hat{s}_{i}$ |  The $i$-th separated source signal.
$\hat{S}$ |  Separated source matrix consisting all $\hat{s}_{i}$.
$Tol$ |  Threshold precision of convergence.
$\mu_{k}$ |  Sampling step for $X$ on node $k$.
$\alpha_{k}$ |  Parameter to denser the sampling from $X$ on node $k$.
$g_{k}$ |  Gradient of Newton’s iteration on node $k$.
$\hat{h}_{g}$ |  Threshold of gradient to Newton’s iteration on node $k$.
$\Delta_{l}$ |  Difference between two consecutive iterations

###  III-C Intermediate Node Processing Logic

####  III-C1 Preliminary

A general BSS problem blindly recovers original data from mixture data without
prior knowledge. One of the popular methods is ICA [5] estimating the solution
matrix $W$ mentioned above. If original data is $i$) statistically independent
and $ii$) non-Gaussian distributed, conventional ICA-based algorithms (e.g.,
FastICA [6]) maximize the non-gaussianity of observed data to separate the
mixtures. This can be done with an approximation of Newton’s Iteration
$\mathcal{N}$:

$(W_{l},\Delta_{l})=\mathcal{N}(X,W_{l-1}),~{}l\in\mathbb{N},$ (2)

where $X$ is the given mixture data, $W_{l}$ is the estimation of separation
matrix in the $l$-th iteration. $\Delta_{l}$ is the difference between two
consecutive iterations, defined by:

$\Delta_{l}=W_{l}\times W_{l-1}^{T}-I,~{}l\in\mathbb{N},$ (3)

where $I$ is an identity matrix. Newton’s Iterations keep updating $W_{l}$
until it converges. It sets a desired tolerance $Tol$ indicating one of the
stoppage criteria (e.g., $\Delta_{l}\leq Tol$) to exit the iterations. The
original ICA above requires all data available on a single point, without
considering the issue of running in a distributed environment.

Based on the classical version, ICA variants are proposed. For example, AeICA
is introduced in [15], where a very basic growing size of subsets of data is
used to calculate the solution matrix $W$ in iterations. However, the
performance of AeICA is heterogeneous case-by-case, because the parameter
controlling the subset data sampling distance is chosen based on an ideal
assumption where the distribution of the data features is uniform. This
assumption does not hold in reality because the prior knowledge on $X$ is
usually unknown. This means that the sampled subset data are not always
representative, which leads to inconsistent separation results.

####  III-C2 Our progressive ICA (pICA)

pICA runs as an in-network microservice. On an intermediate node $k$, $i$) it
keeps forwarding the data to the next node $k+1$; and $ii$) it starts with the
temporal result (i.e., a solution matrix $W_{k-1}$) provided by node $k-1$ and
further improves the result for the next node $k+1$.

Different to AeICA, pICA makes two substantial modifications. The first one is
a new sampling strategy to sample subset data, which eliminates the dependence
on the prior knowledge on input data $X$; and the second one is a new set of
stoppage criteria where node $k$ uses to judge whether or not local iterations
should continue.

Now we introduce the first modification – the new sampling strategy. Our pICA
introduces a controlling parameter $\alpha_{k}$, instead of assuming any prior
knowledge as in AeICA. Every node will dynamically adjust $\alpha_{k}$ in
order to control the variety of the sampled subset data based on the outcome
from the last hop. Specifically, a node $k$ can tune the value of $\alpha_{k}$
so that the sampling step $\mu_{k}$ becomes $\alpha_{k}$ times smaller than
the step value used in the last hop:

$\mu_{k}\leftarrow\frac{\mu_{k-1}}{\alpha_{k}}~{},$ (4)

The decreasing value of $\mu_{k}$ leads to the size of sampled data
sequentially increasing on consecutive nodes. For example, $\alpha_{k}=2$
means on every hop the sampled data are doubled than that of the previous hop,
since the sampling step is halved .

With the sampled subset data, denoted as ${}_{\mu_{k}}X$, instead of re-
computing a solution matrix $W_{k}$, node $k$ continues with the temporal
result $W_{k-1}$ provided by the last hop:

$(W_{k},\Delta_{k})\leftarrow\mathcal{N}({}_{\mu_{k}}X,W_{k-1})~{},$ (5)

where stoppage criteria of a local Newton’s Iteration, which is the second
modification to AeICA, are as follows.

The first criterion is the local convergence tolerance $Tol$, which
characterizes the required precision of solution matrix $W_{k}$. Obviously, if
$W_{k}$ is seen to arrive at the required precision (i.e., $\Delta_{k}\leq
Tol$), the iteration on node $k$ can stop. Note that this does not mean
$W_{k}$ is globally optimal because it is calculated based on a local subset
of data sampled with $\mu_{k}$. On node $k+1$, the size of the subset of data
will increase (e.g., $\alpha_{k+1}=\alpha_{k}\times 2$), so $W_{k}$ will be
further improved if possible.

Another criterion is to indicate whether or not the local marginal gain
becomes too small (recall our greedy strategy). Node $k$ continues only if the
gradient value $g_{k}$ is still large enough in iterations. If the current
$g_{k}$ appears too small ($g_{k}\leq\hat{h}_{g}$), further iterations will
not improve the solution matrix $W_{k}$ significantly anymore. Recall the
existing AeICA, its stoppage criterion however only employs a sole convergence
tolerance based on $\mu_{k}$ without measuring the marginal gain. This leads
AeICA to waste execution time on pursuing little improvements on one node.

After node $k$ exits its iteration, it will relay its solution matrix $W_{k}$
together with the stepping parameter $\mu_{k}$ to node $k+1$. If a node sees
$\mu_{k}$ diminishes to smaller than $1$, which means that all data have been
used on the previous node. Then this node knows that it is the last stop of
the progressive procedure. The action is to run a classical ICA until the
precision of the final solution matrix $W$ satisfies the predefined $Tol$.
Note that most of the jobs for optimizing $W$ have been done on previous
nodes, thus the last step only requires little extra effort. As a final
output, the solution matrix $W$ is applied on the input data $X$ to estimate
original data $\hat{S}$, which will be eventually used for anomaly detection.
The processing logic of node $k$ (i.e., the microservice’s template) is
summarized in Algorithm 1111The algorithm is open sourced at
https://github.com/Huanzhuo/pICA..

input : $X\in\mathbb{R}^{n\times m}$, $W_{k-1}\in\mathbb{R}^{n\times n}$,
$\mu_{k-1}$, $\alpha_{k}$.

output : $W_{k}\in\mathbb{R}^{n\times n}$, $\mu_{k}$, $\alpha_{k+1}$ or
$\hat{S}$

1 Update stepping parameter $\mu_{k}\leftarrow\frac{\mu_{k-1}}{\alpha_{k}}$;

2 Sample subset of sensing data ${}_{\mu_{k}}X$;

3 while  _True_ do

4 $W_{k}\leftarrow\mathcal{N}({}_{\mu_{k}}X,W_{k-1})$ ;

5 if  _$\mu_{k} >1$_ then // Intermediate node of computing.

6 if  _$\Delta_{k} <Tol$_ then // Reach tolerance.

7 $\alpha_{k+1}=\alpha_{k}\times 2$; break;

8

9 else if  _$g_{k} <\hat{h}_{g}$_ then // Small gradient variation.

10 $\alpha_{k+1}=max(2,\frac{\alpha_{k}}{2})$; break;

11

12 $k++$;

13

14 else // Last node of computing.

15 if  _$\Delta_{k} <Tol$_ then

16 $\hat{S}=W_{k}\times X$; break;

17

18 end if

19

20 end if

21

22 end while

Algorithm 1 pICA algorithm on node $k$.

##  IV Numerical Results

###  IV-A Simulation Setup

####  IV-A1 Dataset

We pick a published data set from [16], called MIMII222MIMII: Malfunctioning
Industrial Machine Investigation and Inspection for evaluation. It collects
normal and anomalous operating sound data of $n=4$ types of machines
(including valves, pumps, fans, and slide rails). Every segment is a
$10$-second audio (single-channel and sample rate is $16$kHz). The size of one
data source $s_{i}$ is $m=160k$ (sample rate$\times$duration). Since we have
$4$ types of data sources, the original data $S$ is a $4\times 160k$ matrix.

A $4\times 4$ mixing matrix $A$ generated from a standard distribution will be
applied to the original data matrix $S$ according to Eq. 1 to simulate the
effect of mutual interference. At the AP node, it will receive the data matrix
$X$ with the same size as $S$.

Note that this is just $10$s audio data of $4$ types of machines. In reality,
there will be much more machine types and much longer audio data, time of data
transferring and after that processing them may delay critical anomaly
detection.

####  IV-A2 Scenarios

Our simulation considers five network configurations: $k$ intermediate nodes
($k=0,3,7,10$ and $15$, respectively), representing a low to a high level of
resource availability. For each given number of intermediate nodes $k$, we run
our simulation $50$ times to exhibit the randomness of the mixing matrix $A$.

We implemented the proposed pICA in Python, which can be directly deployed as
microservices. Multi-node network configurations are virtualized on a COTS
server with an i7-6700T CPU with 16GB RAM using Ubuntu 18.04 LTS.

For comparison, we choose FastICA and AeICA discussed in Section II. Note that
FastICA can only be executed on a single node. When $k=0$ (i.e., no
intermediate node available), actually all candidates (FastICA, AeICA, and
pICA) run in a centralized manner.

####  IV-A3 Measured Metrics

The first metric is total _processing time_ $t_{p}$. We use a python module
time (with $1$ microsecond ($\mu s$) precision). $t_{p}$ covers the entire
data processing time in our network system, i.e., from the moment the data
enters the network to the end of the estimated original data $\hat{S}$
derived.

The second metric is _precision_. A Source-to-Distortion Ratio (SDR) metric
from [17] is used to quantify the precision of the estimated original data
$\hat{S}$. Its definition is given in Eq. 6:

$\lx@glossaries@gls@link{acronym}{sdr}{\leavevmode
SDR}=10\cdot\log_{10}\frac{\|s_{truth}\|^{2}}{\|e_{interf}+e_{noise}+e_{artif}\|^{2}}~{},$
(6)

where $s_{truth}$ is the known truth from the picked dataset, $e_{interf}$,
$e_{noise}$, and $e_{artif}$ are the respective errors of interference, noise,
and artifacts errors, which can be evaluated with the output of any separation
algorithm (like our $\hat{S}$) with an open-source BSS Eval Toolbox [18]. SDR
is most widely used metric nowadays because different types of errors are
comprehensively considered.

###  IV-B Processing Time and Precision

Figure 2: Processing time comparison, $\mu_{0}=4130,\alpha_{0}=2$.

Fig. 2 provides a comparison of the processing time of selected algorithms. We
can observe that pICA gains much more speed-ups with increasing numbers of
intermediate nodes $k$. The processing time of pICA decreases from ca. $80ms$
to ca. $45ms$ ($43.75\%$ faster). Specifically, when $k=0$, non-surprisingly,
the three algorithms show similar performances because no acceleration with
intermediate nodes. When $k=3$, pICA starts outperforming the other two. The
gap becomes larger when $k=7$ increasing to $k=15$. This proves that pICA can
accelerate more with intermediate resources compared with others.
Additionally, such gains tend to be deterministic (i.e., smaller variances)
while the other two candidates do not show good stability.

Figure 3: Processing time with different initial sampling steps,
$k=15,\alpha_{0}=2$.

Fig. 3 further presents the impact of how the sampling step distance $\mu_{0}$
on processing time. The result shows that our pICA even prefers starting with
a small subset of data (i.e., larger values of $\mu_{0}$). For example, when
$\mu_{0}=1$ (meaning with all data), pICA does not yield any acceleration at
all; however, with $\mu_{0}=3$, the processing time declines from $65ms$ to
$53ms$; with $\mu_{0}>21$, the processing time constantly reduces to around
$45$ms. On the other hand, since AeICA is very sensitive to the parameter
$\mu_{0}$, which depends on the prior knowledge of the input data, its
performance becomes worse. Thus, it cannot always benefit from available
intermediate nodes. This justifies our proposed growing strategy.

Figure 4: Processing precision comparison, $\mu_{0}=511,\alpha_{0}=2$.

Fig. 4 compares the achieved precision (SDR) with different methods. It proves
that our pICA does not compromise its precision for acceleration but yields an
equivalent precision as the other two methods to restore the original data. It
again justifies the benefits of the progressive and greedy strategies when
pICA is executed hop-by-hop.

###  IV-C Performance of Intermediate Nodes

Figure 5: Intermediate execution costs of pICA, $k=7$ and
$\mu_{0}=500,\alpha_{0}=2$.

At last, we reveal the intermediate procedure by showing the execution costs
of pICA (with $k=7$ and $\mu_{0}=500$) on individual nodes. The result shows
that $60\%$ of processing tasks are finished on intermediate nodes and only
$40\%$ jobs left to the last node. Meanwhile, the precision of the estimated
original data $\hat{S}$ improves quickly ($>60\%$ after node $4$’s processing)
with the size of the subset data growing. It reveals that with small amounts
of data, a solution matrix with relatively high precision can be derived. This
again justifies the effectiveness of our strategy of introducing the only
subset of data on every hop.

##  V Conclusion

This paper provides a novel in-network processing solution suitable with
microservices for acoustic data separation for anomaly detection: progressive
ICA. It is built with in-network nodes and results prove that it indeed
accelerates the processing time of mixture data separation by up to $43.75\%$.
This gives a new representative use case demonstrating the advantage of in-
network processing and how spare resources can be better utilized for critical
tasks.

## Acknowledgment

This work was funded by the Federal Ministry of Education and Research of
Germany (Software Campus Net-BliSS 01IS17044), and by the German Research
Foundation (Deutsche Forschungsgemeinschaft) as part of Germany’s Excellence
Strategy (EXC 2050/1 – Project ID 390696704) Cluster of Excellence “Centre for
Tactile Internet with Human-in-the-Loop” (CeTI) of Technische Universität
Dresden.

## References

  * [1] H. Wu, Y. Shen, X. Xiao, A. Hecker, and F. H. Fitzek, “In-Network processing acoustic data for anomaly detection in smart factory,” in _2021 IEEE Global Communications Conference: IoT and Sensor Networks (Globecom2021 IoTSN)_.
  * [2] M. Wollschlaeger, T. Sauter, and J. Jasperneite, “The future of industrial communication: Automation networks in the era of the internet of things and industry 4.0,” _IEEE Industrial Electronics Magazine_ , vol. 11, no. 1, pp. 17–27, 2017.
  * [3] F. Civerchia, S. Bocchino, C. Salvadori, E. Rossi, L. Maggiani, and M. Petracca, “Industrial internet of things monitoring solution for advanced predictive maintenance applications,” _Journal of Industrial Information Integration_ , vol. 7, pp. 4–12, 2017.
  * [4] D. P. Ellis, “Prediction-driven computational auditory scene analysis,” Ph.D. dissertation, Columbia University, 1996.
  * [5] P. Comon and C. Jutten, _Handbook of Blind Source Separation: Independent component analysis and applications_. Academic press, 2010.
  * [6] A. Hyvarinen, “Fast and robust fixed-point algorithms for independent component analysis,” _IEEE transactions on Neural Networks_ , vol. 10, no. 3, pp. 626–634, 1999.
  * [7] D. T. Pham and P. Garat, “Blind separation of mixture of independent sources through a quasi-maximum likelihood approach,” _IEEE transactions on Signal Processing_ , vol. 45, no. 7, pp. 1712–1725, 1997.
  * [8] H. Wu, Y. Shen, J. Zhang, I. A. Tsokalo, H. Salah, and F. H. P. Fitzek, “Component-dependent independent component analysis for time-sensitive applications,” in _2020 IEEE International Conference on Communications (ICC)_. Dublin, Ireland: IEEE, 2020.
  * [9] H. Wu, J. He, M. Tömösközi, and F. H. P. Fitzek, “Y-net: A dual path model for high accuracy blind source separation,” in _2020 IEEE Globecom Workshops (GC Wkshps)_ , 2020, pp. 1–6.
  * [10] Y. Luo and N. Mesgarani, “Conv-tasnet: Surpassing ideal time–frequency magnitude masking for speech separation,” _IEEE/ACM transactions on audio, speech, and language processing_ , vol. 27, no. 8, pp. 1256–1266, 2019\.
  * [11] I. Kunze, K. Wehrle, and D. Trossen, “Use Cases for In-Network Computing,” Internet Engineering Task Force, Internet-Draft draft-kunze-coin-industrial-use-cases-04, Nov. 2020, work in Progress (Accessed 2021-01-28). [Online]. Available: https://datatracker.ietf.org/doc/html/draft-kunze-coin-industrial-use-cases-04
  * [12] ——, “Transport Protocol Issues of In-Network Computing Systems,” Internet Engineering Task Force, Internet-Draft draft-kunze-coinrg-transport-issues-03, Nov. 2020, work in Progress (Accessed 2021-01-28). [Online]. Available: https://datatracker.ietf.org/doc/html/draft-kunze-coinrg-transport-issues-03
  * [13] Z. Xiang, P. Seeling, and F. H. P. Fitzek, “You only look once, but compute twice: Service function chaining for low-latency object detection in softwarized networks,” _Applied Sciences_ , vol. 11, no. 5, 2021. [Online]. Available: https://www.mdpi.com/2076-3417/11/5/2177
  * [14] R. Glebke, J. Krude, I. Kunze, J. Rüth, F. Senger, and K. Wehrle, “Towards executing computer vision functionality on programmable network devices,” in _Proceedings of the 1st ACM CoNEXT Workshop on Emerging in-Network Computing Paradigms_ , 2019, pp. 15–20.
  * [15] H. Wu, Y. Shen, J. Zhang, H. Salah, I. A. Tsokalo, and F. H. Fitzek, “Adaptive Extraction-Based independent component analysis for Time-Sensitive applications,” in _2020 IEEE Global Communications Conference: Selected Areas in Communications: Internet of Things and Smart Connected Communities (Globecom2020 SAC IoTSCC)_ , Taipei, Taiwan, Dec. 2020.
  * [16] H. Purohit, R. Tanabe, T. Ichige, T. Endo, Y. Nikaido, K. Suefusa, and Y. Kawaguchi, “MIMII Dataset: Sound dataset for malfunctioning industrial machine investigation and inspection,” in _Proceedings of the Detection and Classification of Acoustic Scenes and Events 2019 Workshop (DCASE2019)_ , November 2019, pp. 209–213.
  * [17] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” _IEEE transactions on audio, speech, and language processing_ , vol. 14, no. 4, pp. 1462–1469, 2006.
  * [18] F.-R. Stöter, A. Liutkus, and N. Ito, “The 2018 signal separation evaluation campaign,” in _International Conference on Latent Variable Analysis and Signal Separation_. Springer, 2018, pp. 293–305.