#  Amplifying Dark Sector Particle Non-Resonant Production
in Scattering Dominance Regime

Mingxuan Du<EMAIL_ADDRESS>School of Physics and State Key Laboratory
of Nuclear Physics and Technology, Peking University, Beijing 100871, China
Center for High Energy Physics, Peking University, Beijing 100871, China  Jia
Liu<EMAIL_ADDRESS>School of Physics and State Key Laboratory of Nuclear
Physics and Technology, Peking University, Beijing 100871, China  Center for
High Energy Physics, Peking University, Beijing 100871, China  Xiao-Ping Wang
<EMAIL_ADDRESS>School of Physics, Beihang University, Beijing
100083, China  Beijing Key Laboratory of Advanced Nuclear Materials and
Physics, Beihang University, Beijing 100191, China  Tianhao Wu
<EMAIL_ADDRESS>School of Physics and State Key Laboratory of Nuclear
Physics and Technology, Peking University, Beijing 100871, China

###### Abstract

We investigate the enhancement of dark sector particle production within the
scattering dominant regime. These particles typically exhibit a slight mixing
with Standard Model particles through various portals, allowing for their
generation through in-medium oscillation from Standard Model particle sources.
Our analysis reveals that in the scattering dominance regime, with a
significantly smaller scattering mean free path $\lambda_{\rm sca}$ compared
to the absorption mean free path $\lambda_{\rm abs}$, the non-resonant
production of sterile states can experience an enhancement by a factor of
$\lambda_{\rm abs}/\lambda_{\rm sca}$. This phenomenon is demonstrated within
the context of kinetic mixing dark photon production at a reactor, precisely
satisfying this condition. By incorporating this collisional enhancement, we
find that the current sensitivity to the mixing parameter $\epsilon$ for dark
photons in the TEXONO experiment can be significantly improved across a range
spanning from tens of eV to MeV. This advancement establishes the most
stringent laboratory constraint within this mass spectrum for the dark photon.
Sterile neutrino production, however, does not exhibit such enhancement,
either due to the failure to meet the scattering dominance criterion or the
neutrino damping in resonant production.

##  I Introduction

The dark photon $A^{\prime}$ is a well-motivated new physics particle, which
is the gauge boson corresponding to an extra $U(1)^{\prime}$ gauge symmetry
and can couple to the Standard Model (SM) through a kinetic mixing term [1, 2,
3, 4, 5, 6, 7]. After normalization of kinetic terms, it couples to the
electric charge of the SM particle with a mixing strength $\epsilon$, serving
as a bridge between the SM and the dark sector [8, 9, 10, 11, 12] or a dark
matter candidate as well [13, 14, 15, 16]. Its mass can be originated from the
Stueckelberg or the Higgs mechanisms [17, 18, 19, 20, 21, 22, 23, 24], thus
there are only two phenomenology parameters, mass $m_{A^{\prime}}$ and
interaction strength $\epsilon$ in the model. Recent experiments have
extensively studied the parameter space from particle physics, astrophysics
and cosmology, see Ref. [25, 26] for a recent review on the dark photon and
the various constraints.

Among various terrestrial laboratories [27, 28], nuclear reactor experiments
such as TEXONO [29] and JUNO-TAO [30], can provide the competitive limits on
the dark photon in the mass range of ($10$ eV $\sim$ MeV). At the core of a
nuclear reactor, numerous photons with MeV energy are produced, which can
potentially transform into dark photons with masses up to MeV if allowed by
the kinematics. If the dark photon mass is well below the effective photon
mass $m_{\gamma}$ (which is around 10 eV in the reactor [29]), the conversion
probability will suffer the suppression $\sim
m_{A^{\prime}}^{4}/m_{\gamma}^{4}$ [31, 32]. Consequently, the nuclear reactor
experiments are the suitable place to explore the high mass range $10$ eV
$\sim$ MeV.

Estimating dark photon production in nuclear reactor cores is performed in the
interaction basis, where the dark photon is sterile to the SM electromagnetic
current. Consequently, the dark photon can be generated only through the
oscillation between the SM photon and the dark photon when photons are
produced in the reactor. This scenario requires considering the interaction
between the SM photon and reactor cores, which directly affects the
oscillation process. For photons in reactor cores, two important effects are
addressed: absorption, which terminates the oscillation, and scattering with
the medium, which destroys the old photon and generate a new one with a softer
energy. We emphasize that the scattering rate is significantly larger than the
absorption rate for MeV reactor photons [33], resulting in multiple
scatterings before absorption occurs. Previous studies only considered
absorption when calculating the conversion of photons to dark photons in the
reactor [29, 34]. However, this study will also include the scattering effects
and examine how they modify the results.

##  II The dark photon and photon oscillation scheme

In the kinetic-mixing dark photon model [5], the Hamiltonian for the
transverse modes of the photon and dark photon in the medium in the
interaction basis is given by [35],

$\displaystyle
H=\frac{1}{2E_{\gamma}}\left(\begin{array}[]{cc}\epsilon^{2}m_{A^{\prime}}^{2}+m_{\gamma}^{2}&-\epsilon
m_{A^{\prime}}^{2}\\\ -\epsilon
m_{A^{\prime}}^{2}&m_{A^{\prime}}^{2}\end{array}\right),$ (3)

where the dark photon $A^{\prime}$ has mass $m_{A^{\prime}}$, $E_{\gamma}$ is
the (dark) photon energy, and $m_{\gamma}=\sqrt{4\pi\alpha n_{e}/m_{e}}$ is
the effective mass of photon induced by the medium, which is about
$m_{\gamma}=20$ eV [29] in the reactor core. Our analysis focuses only on the
non-resonance region, specifically when $m_{A^{\prime}}\gg m_{\gamma}$ [29],
as exploring the resonance region requires a detailed understanding of the
reactor core structure and materials.

We do not include the longitudinal photon in our calculation because its
energy, equal to the plasma frequency (approximately 20 eV), is below the
threshold of the neutrino detector used in this study. Additionally, the
longitudinal plasmon is not a propagation mode, so the physical picture of
Compton scattering during propagation is not valid. Therefore, we only
consider the transverse mode of the photon in our study.

##  III The absorption and scattering in the medium

We will examine matter effects of photon oscillations to dark photons in a
medium, including absorption and scattering. When dealing with MeV photons in
the reactor, we will consider the photoelectric process for absorption and
Compton scattering for scattering, as they are the dominant processes [33].

The damping processes (absorption and scattering) terminate photon propagation
and oscillation to the dark photon. The conversion probability with damping
effect is given by [36]

$P_{\gamma\to
A^{\prime}}(\ell,\Gamma,E_{\gamma})=\epsilon^{2}m_{A^{\prime}}^{4}\frac{1+e^{-\Gamma\ell}-2\cos\frac{\Delta
m^{2}\ell}{2E_{\gamma}}e^{-\frac{\Gamma\ell}{2}}}{\left(\Delta
m^{2}\right)^{2}+(E_{\gamma}\Gamma)^{2}},$ (4)

where the damping rate is $\Gamma=\Gamma_{\rm abs}+\Gamma_{\rm sca}$. Here
$\Gamma_{\rm abs}=n_{A}\sigma_{{pe}}$ represents the absorption rate for the
photoelectric process, and $\Gamma_{\rm sca}=n_{A}\sigma_{C}$ is the Compton
scattering rate. Assuming a pure thorium ($Z=90$) reactor core [37], $n_{A}$
is the corresponding number density of thorium atoms, and $\sigma_{pe}$ and
$\sigma_{C}$ denote the cross-section for photoelectric process and Compton
scattering in pure thorium [38]. By neglecting the scattering rate and
consider a constant absorption rate $\Gamma_{\rm abs}\simeq 1/(10\rm{cm})$ for
MeV photons in reactor cores [29], for $\ell\to\infty$, Eq. (4) recovers to
[39, 40, 35, 36]:

$P_{\gamma\to
A^{\prime}}^{\rm{a}}=\epsilon^{2}\frac{m_{A^{\prime}}^{4}}{(\Delta
m^{2})^{2}+E_{\gamma}^{2}\Gamma_{\rm abs}^{2}},$ (5)

which is the conversion probability of a photon oscillating to a dark photon
at infinity solely under the influence of absorption.

Note that $\Gamma_{\rm sca}\gg\Gamma_{\rm abs}$ for MeV photons in reactor
cores [33, 38]. Therefore, the scattering process $\gamma+e\to\gamma+e$ will
significantly impact the oscillation by absorbing and emitting photons,
altering their energy and direction. It is essential to incorporate these
newly released photons into the simulation and recalculate the oscillation
probability. The reactor core produces isotropic photons [41], and we assume a
homogeneous core [37] which is much larger ($O(\rm{m})$) than the absorption
length ($O(10\leavevmode\nobreak\ \rm{cm})$) of MeV photons [29]. These
assumptions lead us to expect isotropic distribution of resulting reactor dark
photons. Thus, our analysis focuses solely on the resulting dark photon
spectrum. Notably, a photon becoming a dark photon retains its original energy
during oscillation. Therefore, our study considers only energy changes in
Compton scattering as the direction change does not affect the dark photon
energy spectrum.

##  IV Simulation process

In our simulation, the absorption and scattering effects in the reactor cores
are explored by tracking the individual photon. We consider a photon generated
from the reactor core with initial energy $E_{1}$. This photon first travels
through a mean free path of a Compton scattering process, $\lambda_{C}$, in
pure thorium medium [38]

$\lambda_{C}(E_{1})=1/(n_{A}\sigma_{C}(E_{1})).$ (6)

For a propagation distance smaller than $\lambda_{C}$, no Compton scattering
is assumed to take place. Then the probability of the photon oscillating to a
dark photon with absorption effect in this distance is given by $P_{\gamma\to
A^{\prime},\ 1}(\lambda_{C}(E_{1}),\Gamma_{\rm abs}(E_{1}),E_{1})$ as Eq. (4),
where the subscript $1$ records the step of the propagation. Considering loss
of photon due to the medium absorption and the oscillation to dark photon, the
resulting survival probability for photon is

$\displaystyle P^{s}_{1}=\left[1-P_{\gamma\to A^{\prime},\
1}\right]e^{-\lambda_{C}\Gamma_{\rm abs}}.$ (7)

After traveling the mean free path distance $\lambda_{C}$, the photon is
assumed to undergo next Compton scattering immediately, where the old photon
with energy $E_{1}$ is absorbed by the medium and a new photon with lower
energy $E_{2}<E_{1}$ is generated. The energy $E_{2}$ of the new photon is
randomly chosen based on the final state photon energy distribution in the
laboratory frame for Compton scattering.

Figure 1: The flow chart of our simulations for reactor photons oscillating to
dark photons with Compton scattering included.

In the next step, the new photon travels through the medium over another mean
free path distance $\lambda_{C}(E_{2})$ and undergoes the processes of
propagation, absorption, oscillation, and scattering, see Fig. 1 for our
simulation flow chart. We loop these processes until the survival probability
of the final photon falls below 0.1, denoted as $\Pi_{i=0}^{n}P^{s}_{i}$ (with
initial value $P^{s}_{0}=1$), where $i$ is the step index and $n$ is the
number of loops when the survival probability falls below 0.1. Here we choose
0.1 as the threshold to reduce the computing power requirement, since a
smaller number will not significantly change the production of dark photons.

Finally, after considering absorption and scattering effects, we simulated the
total conversion probability from a photon to a dark photon in the medium,
which can be expressed as

$P_{\gamma\to A^{\prime}}=\sum_{i=1}^{n}\left(P^{a}_{\gamma\to A^{\prime},\
i}\times\Pi_{j=0}^{i-1}P^{s}_{j}\right).$ (8)

It sums up the probabilities of generating dark photons in each step along the
photon propagation path. Furthermore, to suppress the fluctuations in our
results, we simulated $1000$ photons for each initial energy. In Fig. 2, when
counting the scattering effect, the conversion probability is significantly
enhanced for $m_{A^{\prime}}>20$ eV.

Figure 2: The rescaled conversion probability for a photon with energy 1 MeV
to dark photon in reactor cores under different scenarion. The solid curve
considers both the absorption and Compton scattering effects. The dashed curve
only considers the absorption effect.

##  V Analytic estimation

To comprehend the simulation results, we aim to derive an analytical
explanation for the observed enhancement in the total conversion probability.
Eq. (4) contains two lengths for the MeV photon: interaction length
$\lambda=1/\Gamma\sim 2$ cm and oscillation length $\ell_{\rm osc}=4\pi
E_{\gamma}/\Delta m^{2}\sim 10^{-10}(m_{A^{\prime}}/\rm{MeV})^{-2}$ cm.
Therefore, $\ell_{\rm osc}$ leads to a highly oscillating pattern with
propagation length $\ell$. To simplify, we disregard the highly oscillating
linear term $\cos{(l/l_{\rm osc})}$ in Eq. (4), and obtain an average
estimation of the single conversion probability. Then the average probability
of a photon not undergoing any scattering and converting to a non-resonant
dark photon, is

$P_{0}(\ell)=\frac{\epsilon^{2}m_{A^{\prime}}^{4}\left(1+e^{-\Gamma\ell}\right)}{\left(\Delta
m^{2}\right)^{2}+(E_{\gamma}\Gamma)^{2}}\simeq\frac{\epsilon^{2}m_{A^{\prime}}^{4}\left(1+e^{-\Gamma\ell}\right)}{\left(\Delta
m^{2}\right)^{2}}.$ (9)

We drop the $E_{\gamma}\Gamma$ term in the denominator since $\Delta m^{2}\sim
m_{A^{\prime}}^{2}\gg E_{\gamma}\Gamma\sim(3\leavevmode\nobreak\ \rm{eV})^{2}$
for non-resonant dark photon production at the reactor, resulting in an
energy-independent conversion probability.

For a rough estimate, assuming no energy loss during medium scattering, the
probability of a photon undergoing $n$ times of scatterings at locations
$0<\ell_{1}<\ell_{2},...<\ell_{n}<\ell$ and converting to a non-resonant dark
photon is (see a detailed derivation in Appendix A)

$\displaystyle P_{n}(\ell)$ $\displaystyle=$
$\displaystyle\int_{0}^{\ell}d\ell_{1}\int_{\ell_{1}}^{\ell}d\ell_{2}...\int_{\ell_{n-1}}^{\ell}d\ell_{n}e^{-\Gamma\ell_{n}}\Gamma_{\rm
sca}^{n}P_{0}(\ell-\ell_{n})$ (10) $\displaystyle=$
$\displaystyle\frac{\epsilon^{2}m_{A^{\prime}}^{4}}{(\Delta
m^{2})^{2}}\frac{\Gamma_{\rm
sca}^{n}}{\Gamma^{n}}\left(1+e^{-x}\frac{x^{n}}{n!}-\frac{\Gamma(n,x)}{\Gamma(n)}\right),$

where $x=\Gamma\ell$, while $\Gamma(n,x)$ and $\Gamma(n)$ are the incomplete
and complete Gamma functions, respectively. The damping effect on the active
state at the $n$-th scattering is described by $e^{-\Gamma\ell_{n}}$,
$\Gamma_{\rm sca}^{n}$ represents the $n$ times of scatterings, and the last
term $P_{0}(\ell-\ell_{n})$ gives the conversion probability from the new
photon generated by the $n$-th scattering.

The total conversion probability, summing up the multi-scattering between
photons and medium, leads to

$\displaystyle P^{\rm tot}$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}P_{i}=\frac{\epsilon^{2}m_{A^{\prime}}^{4}}{(\Delta
m^{2})^{2}}\left(\frac{\Gamma}{\Gamma_{\rm abs}}+B(\ell)\right),$ (11)
$\displaystyle B(\ell)$ $\displaystyle\equiv$ $\displaystyle e^{-\Gamma_{\rm
abs}\ell}-\sum_{n=1}^{\infty}\frac{\Gamma_{\rm
sca}^{n}}{\Gamma^{n}}\frac{\Gamma(n,\Gamma\ell)}{\Gamma(n)}.$

For $\ell\to 0$, $P^{\rm tot}(0)=P_{0}(0)$ thus the short propagation distance
suppresses the multiple scattering effect. When the active state propagates
long distance in the medium, $\ell\gtrsim 1/\Gamma_{\rm abs}$, one has

$P^{\rm tot}(\ell)\approx\frac{\Gamma}{\Gamma_{\rm abs}}P_{0}(\infty).$ (12)

In the scattering dominant region where $\Gamma\sim\Gamma_{\rm
sca}\gg\Gamma_{\rm abs}$, the non-resonant dark photon production is enhanced
by the factor $\Gamma_{\rm sca}/\Gamma_{\rm abs}\sim\lambda_{\rm
abs}/\lambda_{\rm sca}$ compared to the case without considering scattering.
For MeV photons in reactor, this enhancement factor is about
$10\leavevmode\nobreak\ \rm{cm}/2\leavevmode\nobreak\ \rm{cm}\sim 5$, which in
the ballpark agrees with the factor of $10$ in our full numeric simulation for
the non-resonant region $m_{A^{\prime}}>50$ eV in Fig. 2.

Comparing $\gamma-A^{\prime}$ mixing with sterile neutrino mixing is
intriguing. Unlike non-resonant dark photons in reactors, the production of
sterile neutrinos in supernovae lacks a similar enhancement. Two mechanisms
drive sterile neutrino production in supernovae: adiabatic MSW resonance and
collisional production [42, 43, 44, 45]. The dominant mechanism, adiabatic MSW
resonance, loses efficacy due to a high scattering rate, significantly
reducing neutrino mean free path [44]. Conversely, collisional production is
less influenced by scattering rate, as the damping factor of $\Gamma^{2}$ in
conversion probability is counterbalanced by scattering rate in the numerator
and enhanced resonant length [39]. Thus, elevated scattering rates generally
suppress sterile neutrino production in supernovae. In contrast, for reactor
dark photons in the non-resonant regime, the mass term $(\Delta m^{2})^{2}$
outweighs the significance of the scattering damping term $E^{2}\Gamma^{2}$ in
the conversion probability. Therefore, increased scattering rates linearly
amplify dark photon production at the reactor.

##  VI The reactor experiment search of the dark photon

Nuclear reactor cores produce numerous MeV photons, making them ideal for dark
photon production. Assuming isotropic generation of prompt photons, the
approximate spectrum for photon flux at $E_{\gamma}>0.2$ MeV is [41]:

$\frac{dF_{\gamma}}{dE_{\gamma}}=\frac{0.58\times
10^{21}}{\mathrm{s}\,\mathrm{MeV}}\left(\frac{P}{\mathrm{GW}}\right)\exp\left[-\frac{1.1E_{\gamma}}{\mathrm{MeV}}\right],$
(13)

where $P$ represents the reactor thermal power. The dark photon production
rate from photon oscillation can be expressed as

$\frac{dF_{A^{\prime}}}{dE_{A^{\prime}}}=\int
dE_{\gamma}\frac{dF_{\gamma}}{dE_{\gamma}}\frac{dP_{\gamma\to
A^{\prime}}}{dE_{A^{\prime}}}(E_{\gamma}),$ (14)

where ${dP_{\gamma\to A^{\prime}}}/{dE_{A^{\prime}}}$ represents the
differential conversion probability in our simulation. Fig. 3 illustrates the
energy spectrum for dark photons produced by a 2.9 GW nuclear reactor with
$m_{A^{\prime}}=0.1$ MeV. The solid line accounts for absorption and
scattering effects, while the dashed line only considers absorption. Since the
scattering softens the photon energy and $E_{\gamma}=E_{A^{\prime}}$ in the
oscillation, more low-energy $A^{\prime}$ are produced with inclusion of
scatterings, resulting in a larger enhancement for the lower energy dark
photon in the solid line compared to the dashed line.

Figure 3: The spectra of dark photon with $m_{A^{\prime}}=0.1$ MeV produced in
the reactor with 2.9 GW thermal power under different scenarios. The solid
curve considers both the absorption and Compton scattering. The dashed curve
only considers the absorption effect.

The dark photons, with negligible attenuation due to oscillation back into
photons ($\epsilon^{2}\ll 1$), can be detected through Compton scattering of
electrons once they reach the detector by measuring the electron recoil
energy. We analyzed experimental data [46, 47, 48] from high-energy threshold
($3\sim 8$ MeV) CsI(Tl) scintillating detectors [46] and low-energy threshold
($0.3\sim 12.3$ keV) Germanium (Ge) ionization detectors [47] at TEXONO [49].
Then in Fig. 4, we set $95\%$ C.L. constraints on dark photons with a mass
lighter than $1$ MeV obtained from the CsI (red) and PCGe (blue) detectors in
TEXONO experiment, see Appendix B for more details. For dark photon heavier
than $1$ MeV, it can quickly decay to a pair of electrons thus does not apply
to the scenario described here. The solid curves represent the constraints
considering both absorption and scattering effects of the medium on photons,
while the dashed curves only consider the absorption effect.

Figure 4: TEXONO $95\%$ C.L. constraints on dark photons lighter than $1$ MeV
from CsI (red) and PCGe (blue) detectors. The solid curves are the constraints
incorporating both absorption and scattering effects on photons by the medium,
and the dashed curves incorporate only the absorption. Existing constraints
from terrestrial experiments are shown as dark gray regions [50, 51, 52, 25].
Astrophysical and cosmological constraints are shown as light gray regions
[39, 35, 53, 54, 55].

##  VII Results and discussion

The inclusion of scatterings increases the conversion probability of
oscillation into the dark photon by a factor of 10 in Fig. 2. The number of
signals $N_{\text{sig}}$ in the detector is proportional to $\epsilon^{4}$, so
the $\epsilon^{2}$ limit will be enhanced about $\sqrt{10}\sim 3$ times.
Additionally, scattering softens the energy of the photons, reducing
$N_{\text{sig}}$ for high-energy threshold detectors (e.g., CsI detector with
threshold of $3\sim 8$ MeV) and increasing $N_{\text{sig}}$ for low-energy
threshold detectors (e.g., PCGe detector with threshold of $0.3\sim 12.3$
keV). Consequently, the $\epsilon^{2}$ limit is enhanced by approximately 2
times for CsI detectors limits, and 6 times for PCGe detectors limits. This
highlights the significance of considering the scattering effect in dark
photon constraints. Consequently, TEXONO provides the most stringent
constraints on dark photons among terrestrial experiments [50, 51, 52, 25]
Within the range of dark photon masses spanning from 10 eV to MeV.

The photon energy spectrum in Eq. (13) only includes photons with energies
$E_{\gamma}>0.2$ MeV. Thus, when $m_{A^{\prime}}>0.2$ MeV, low-energy photons
cannot pass through the reactor to reach the detector by oscillating to heavy
dark photons, resulting in a decrease in the signal rate of low electron
recoil energy. Consequently, the PCGe detector with a low-energy threshold of
$0.3\sim 12.3$ keV has a weakened sensitivity in the mass range of $0.2\sim 1$
MeV. In contrast, the CsI detector with a high-energy threshold of $3\sim 8$
MeV, remains constant in the $0.2\sim 1$ MeV mass range.

##  VIII Conclusion

In this study, we have calculated the detection flux for the conversion
probability of $\gamma\to A^{\prime}$ in neutrino experiments, such as TEXONO.
Unlike conventional methods that solely consider the absorption process, we
also take into account the multiple scattering effect resulting from Compton
scattering in the medium. Our findings indicate that this multiple Compton
scattering effect significantly enhances the sensitivity to dark photons
across all probed mass region. This kind of enhancement does not exhibit in
the sterile neutrino production, either due to the failure to meet the
scattering dominance criterion or the neutrino damping in resonant production.
In the future, the new reactor experiment, such as JUNO [56] and NEOS [57] can
provide better sensitivity for dark photon with the inclusion of the multiple
scattering effect.

###### Acknowledgements.

We thank Zuowei Liu and Shuailiang Ge for useful discussions. The work of JL
is supported by the National Science Foundation of China under Grant No.
12075005, 12235001. The work of X.P.W. is supported by the National Science
Foundation of China under Grant No. 12005009 and the Fundamental Research
Funds for the Central Universities.

## References

  * [1] P. Fayet, “Effects of the Spin 1 Partner of the Goldstino (Gravitino) on Neutral Current Phenomenology,” Phys. Lett. B 95 (1980) 285–289.
  * [2] P. Fayet, “On the Search for a New Spin 1 Boson,” Nucl. Phys. B 187 (1981) 184–204.
  * [3] L. B. Okun, “LIMITS OF ELECTRODYNAMICS: PARAPHOTONS?” Sov. Phys. JETP 56 (1982) 502.
  * [4] P. Galison and A. Manohar, “TWO Z’s OR NOT TWO Z’s?” Phys. Lett. B 136 (1984) 279–283.
  * [5] B. Holdom, “Two U(1)’s and Epsilon Charge Shifts,” Phys. Lett. B 166 (1986) 196–198.
  * [6] P. Fayet, “Extra U(1)’s and New Forces,” Nucl. Phys. B 347 (1990) 743–768.
  * [7] R. Foot and X.-G. He, “Comment on Z Z-prime mixing in extended gauge theories,” Phys. Lett. B 267 (1991) 509–512.
  * [8] M. Pospelov, A. Ritz, and M. B. Voloshin, “Secluded WIMP Dark Matter,” Phys. Lett. B 662 (2008) 53–61 [arXiv:0711.4866].
  * [9] N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, and N. Weiner, “A Theory of Dark Matter,” Phys. Rev. D 79 (2009) 015014 [arXiv:0810.0713].
  * [10] R. Essig et al. in Community Summer Study 2013: Snowmass on the Mississippi.  2013\.  arXiv:1311.0029.
  * [11] J. Alexander et al. 2016\.  arXiv:1608.08632.
  * [12] M. Battaglieri et al. in U.S. Cosmic Visions: New Ideas in Dark Matter.  2017\.  arXiv:1707.04591.
  * [13] J. Redondo and M. Postma, “Massive hidden photons as lukewarm dark matter,” JCAP 02 (2009) 005 [arXiv:0811.0326].
  * [14] A. E. Nelson and J. Scholtz, “Dark Light, Dark Matter and the Misalignment Mechanism,” Phys. Rev. D 84 (2011) 103501 [arXiv:1105.2812].
  * [15] P. Arias, et al., “WISPy Cold Dark Matter,” JCAP 06 (2012) 013 [arXiv:1201.5902].
  * [16] P. W. Graham, J. Mardon, and S. Rajendran, “Vector Dark Matter from Inflationary Fluctuations,” Phys. Rev. D 93 (2016) 103520 [arXiv:1504.02102].
  * [17] B. Kors and P. Nath, “Aspects of the Stueckelberg extension,” JHEP 07 (2005) 069 [hep-ph/0503208].
  * [18] D. Feldman, Z. Liu, and P. Nath, “Probing a very narrow Z-prime boson with CDF and D0 data,” Phys. Rev. Lett. 97 (2006) 021801 [hep-ph/0603039].
  * [19] D. Feldman, Z. Liu, and P. Nath, “The Stueckelberg $Z$ Prime at the LHC: Discovery Potential, Signature Spaces and Model Discrimination,” JHEP 11 (2006) 007 [hep-ph/0606294].
  * [20] D. Feldman, Z. Liu, and P. Nath, “The Stueckelberg Z-prime Extension with Kinetic Mixing and Milli-Charged Dark Matter From the Hidden Sector,” Phys. Rev. D 75 (2007) 115001 [hep-ph/0702123].
  * [21] D. Feldman, Z. Liu, P. Nath, and B. D. Nelson, “Explaining PAMELA and WMAP data through Coannihilations in Extended SUGRA with Collider Implications,” Phys. Rev. D 80 (2009) 075001 [arXiv:0907.5392].
  * [22] M. Du, Z. Liu, and V. Q. Tran, “Enhanced Long-Lived Dark Photon Signals at the LHC,” JHEP 05 (2020) 055 [arXiv:1912.00422].
  * [23] M. Du, R. Fang, Z. Liu, and V. Q. Tran, “Enhanced long-lived dark photon signals at lifetime frontier detectors,” Phys. Rev. D 105 (2022) 055012 [arXiv:2111.15503].
  * [24] M. Redi and A. Tesi, “Dark photon Dark Matter without Stueckelberg mass,” JHEP 10 (2022) 167 [arXiv:2204.14274].
  * [25] M. Fabbrichesi, E. Gabrielli, and G. Lanfranchi, “The Dark Photon.” arXiv:2005.01515.
  * [26] A. Caputo, A. J. Millar, C. A. J. O’Hare, and E. Vitagliano, “Dark photon limits: A handbook,” Phys. Rev. D 104 (2021) 095029 [arXiv:2105.04565].
  * [27] P. Ilten, Y. Soreq, M. Williams, and W. Xue, “Serendipity in dark photon searches,” JHEP 06 (2018) 004 [arXiv:1801.04847].
  * [28] M. Bauer, P. Foldenauer, and J. Jaeckel, “Hunting All the Hidden Photons,” JHEP 07 (2018) 094 [arXiv:1803.05466].
  * [29] M. Danilov, S. Demidov, and D. Gorbunov, “Constraints on hidden photons produced in nuclear reactors,” Phys. Rev. Lett. 122 (2019) 041801 [arXiv:1804.10777].
  * [30] M. Smirnov, G. Yang, J. Liao, Z. Hu, and J. Ling, “Light dark bosons in the JUNO-TAO neutrino detector,” Phys. Rev. D 104 (2021) 116024 [arXiv:2109.04276].
  * [31] T. Mizumoto, et al., “Experimental search for solar hidden photons in the eV energy range using kinetic mixing with photons,” JCAP 07 (2013) 013 [arXiv:1302.1000].
  * [32] J.-F. Fortin and K. Sinha, “Photon-dark photon conversions in extreme background electromagnetic fields,” JCAP 11 (2019) 020 [arXiv:1904.08968].
  * [33] Particle Data Group Collaboration, “Review of Particle Physics,” PTEP 2022 (2022) 083C01.
  * [34] S. H. Seo and Y. D. Kim, “Dark Photon Search at Yemilab, Korea,” JHEP 04 (2021) 135 [arXiv:2009.11155].
  * [35] J. Redondo and G. Raffelt, “Solar constraints on hidden photons re-visited,” JCAP 08 (2013) 034 [arXiv:1305.2920].
  * [36] J. Redondo, “Atlas of solar hidden photon emission,” JCAP 07 (2015) 024 [arXiv:1501.07292].
  * [37] J. B. Dent, et al., “New Directions for Axion Searches via Scattering at Reactor Neutrino Experiments,” Phys. Rev. Lett. 124 (2020) 211804 [arXiv:1912.05733].
  * [38] M. Berger, J. Hubbell, S. Seltzer, J. Coursey, and D. Zucker, “XCOM: Photon Cross Section Database (version 1.2),” 1999.
  * [39] J. Redondo, “Helioscope Bounds on Hidden Sector Photons,” JCAP 07 (2008) 008 [arXiv:0801.1527].
  * [40] H. An, M. Pospelov, and J. Pradler, “New stellar constraints on dark photons,” Phys. Lett. B 725 (2013) 190–195 [arXiv:1302.3884].
  * [41] H. S. H. Bechteler, H. Fasissner and Yogeshwar, IKP Annual report.  1984\.
  * [42] L. Stodolsky, “On the Treatment of Neutrino Oscillations in a Thermal Environment,” Phys. Rev. D 36 (1987) 2273.
  * [43] G. G. Raffelt and S. Zhou, “Supernova bound on keV-mass sterile neutrinos reexamined,” Phys. Rev. D 83 (2011) 093014 [arXiv:1102.5124].
  * [44] C. A. Argüelles, V. Brdar, and J. Kopp, “Production of keV Sterile Neutrinos in Supernovae: New Constraints and Gamma Ray Observables,” Phys. Rev. D 99 (2019) 043012 [arXiv:1605.00654].
  * [45] A. Ray and Y.-Z. Qian, “Evolution of Tau-Neutrino Lepton Number in Protoneutron Stars due to Active-Sterile Neutrino Mixing.” arXiv:2306.08209.
  * [46] TEXONO Collaboration, “Measurement of Nu(e)-bar -Electron Scattering Cross-Section with a CsI(Tl) Scintillating Crystal Array at the Kuo-Sheng Nuclear Power Reactor,” Phys. Rev. D 81 (2010) 072001 [arXiv:0911.1597].
  * [47] TEXONO Collaboration, “Characterization and Performance of Germanium Detectors with sub-keV Sensitivities for Neutrino and Dark Matter Experiments,” Nucl. Instrum. Meth. A 836 (2016) 67–82 [arXiv:1411.4802].
  * [48] TEXONO Collaboration, “Constraints on millicharged particles with low threshold germanium detectors at Kuo-Sheng Reactor Neutrino Laboratory,” Phys. Rev. D 99 (2019) 032009 [arXiv:1808.02719].
  * [49] TEXONO Collaboration, “A CsI(Tl) scintillating crystal detector for the studies of low-energy neutrino interactions,” Nucl. Instrum. Meth. A 459 (2001) 93 [hep-ex/0001001]. [Erratum: Nucl.Instrum.Meth.A 485, 821 (2002)].
  * [50] J. Jaeckel and S. Roy, “Spectroscopy as a test of Coulomb’s law: A Probe of the hidden sector,” Phys. Rev. D 82 (2010) 125020 [arXiv:1008.3536].
  * [51] M. Schwarz, et al., “Results from the Solar Hidden Photon Search (SHIPS),” JCAP 08 (2015) 011 [arXiv:1502.04490].
  * [52] D. Banerjee et al., “Dark matter search in missing energy events with NA64,” Phys. Rev. Lett. 123 (2019) 121801 [arXiv:1906.00176].
  * [53] N. Vinyoles, et al., “New axion and hidden photon constraints from a solar data global fit,” JCAP 10 (2015) 015 [arXiv:1501.01639].
  * [54] J. H. Chang, R. Essig, and S. D. McDermott, “Revisiting Supernova 1987A Constraints on Dark Photons,” JHEP 01 (2017) 107 [arXiv:1611.03864].
  * [55] D. K. Hong, C. S. Shin, and S. Yun, “Cooling of young neutron stars and dark gauge bosons,” Phys. Rev. D 103 (2021) 123031 [arXiv:2012.05427].
  * [56] JUNO Collaboration, “Neutrino Physics with JUNO,” J. Phys. G 43 (2016) 030401 [arXiv:1507.05613].
  * [57] NEOS Collaboration, “Development and Mass Production of a Mixture of LAB- and DIN-based Gadolinium-loaded Liquid Scintillator for the NEOS Short-baseline Neutrino Experiment,” J. Radioanal. Nucl. Chem. 310 (2016) 311–316 [arXiv:1511.05551].

##  Appendix A Derivation of conversion probability with multiple scattering

In this section, we provide the derivation of Eq. (10) in the main text, for
the probability of a photon undergo $n$ times of scatterings and converting to
a non-resonant dark photon.

Taking into account the damping effect of the medium (absorption and
scattering), the scattering probability of a photon at differential location
of $\ell\sim\ell+d\ell$ is

$dP_{\rm sca}=e^{-\Gamma\ell}(1-e^{-\Gamma_{\rm sca}d\ell})\simeq
e^{-\Gamma\ell}\Gamma_{\rm sca}d\ell.$ (S1)

Thus, assuming no energy loss during medium scattering for a rough estimate,
the probability of a photon being scattered once by the medium at
$0<l_{1}<\ell$ and then converting into a dark photon is

$P_{1}(\ell)=\int_{0}^{\ell}dl_{1}e^{-\Gamma l_{1}}\Gamma_{\rm sca}P_{0}(\ell-
l_{1}).$ (S2)

Similarly, the conversion probability with $n$ scatterings can be obtained
using the following iterative relationship

$P_{n}(\ell)=\int_{0}^{\ell}dl_{n}e^{-\Gamma l_{n}}\Gamma_{\rm
sca}P_{n-1}(\ell-l_{n}).$ (S3)

Here, the photon firstly scatters once with the medium at $l_{n}$, then
scatters $n-1$ times along the remaining propagation path $\ell-l_{n}$ and
converts into a dark photon. By explicitly expressing this iterative
relationship, we have

$\displaystyle P_{n}(\ell)$ $\displaystyle=$
$\displaystyle\int_{0}^{\ell}dl_{n}e^{-\Gamma l_{n}}\Gamma_{\rm
sca}\int_{0}^{\ell-l_{n}}dl_{n-1}e^{-\Gamma l_{n-1}}\Gamma_{\rm
sca}...\int_{0}^{\ell-\sum_{i=2}^{n}l_{i}}dl_{1}e^{-\Gamma l_{1}}\Gamma_{\rm
sca}P_{0}(\ell-\sum_{i=1}^{n}l_{i})$ (S4) $\displaystyle=$
$\displaystyle\int_{0}^{\ell}dl_{n}\int_{0}^{\ell-
l_{n}}dl_{n-1}...\int_{0}^{\ell-\sum_{i=2}^{n}l_{i}}dl_{1}e^{-\Gamma\sum_{i=1}^{n}l_{i}}\Gamma_{\rm
sca}^{n}P_{0}(\ell-\sum_{i=1}^{n}l_{i})$ $\displaystyle=$
$\displaystyle\int_{0}^{\ell}d\tilde{\ell}_{1}\int_{0}^{\ell-\tilde{\ell}_{1}}d\tilde{\ell}_{2}...\int_{0}^{\ell-\sum_{i=1}^{n-1}\tilde{\ell}_{i}}d\tilde{\ell}_{n}e^{-\Gamma\sum_{i=1}^{n}\tilde{\ell}_{i}}\Gamma_{\rm
sca}^{n}P_{0}(\ell-\sum_{i=1}^{n}\tilde{\ell}_{i})$ $\displaystyle=$
$\displaystyle\int_{0}^{\ell}d\ell_{1}\int_{\ell_{1}}^{\ell}d\ell_{2}...\int_{\ell_{n-1}}^{\ell}d\ell_{n}e^{-\Gamma\ell_{n}}\Gamma_{\rm
sca}^{n}P_{0}(\ell-\ell_{n}).$

In the last two steps, we employ the following variable transformations

$\tilde{\ell}_{i}=l_{n+1-i},\quad\ell_{i}=\sum_{j=1}^{i}\tilde{\ell}_{j}.$
(S5)

The last line of Eq. (S4) corresponds to the probability of a photon
undergoing $n$ times of scatterings at locations
$0<\ell_{1}<\ell_{2},...<\ell_{n}<\ell$ and converting to a non-resonant dark
photon. Thus, we arrive at Eq. (10) in the main text.

##  Appendix B The reactor experiment search of the dark photon

In this section, we briefly describe the two TEXONO detectors used in our
analysis: CsI(Tl) scintillating detectors and Germanium (Ge) ionization
detectors. We provide a detailed description on how to utilize their data to
establish the dark photon limit.

The TEXONO neutrino experiment is located 28 m from the reactor core of the
Kuo-Sheng Nuclear Power station with thermal power of 2.9 GW [49]. It uses
CsI(Tl) scintillating detectors [49] and Germanium (Ge) ionization detectors
[47] to search for neutrinos through $\bar{\nu}_{e}-e$ elastic scattering.

The CsI(Tl) scintillating crystal array measures $40\times 31.2\times
37\leavevmode\nobreak\ \text{cm}^{3}$. For CsI(Tl) detectors, the energy
thresholds for the electron recoil energy are $E_{e}^{\min}=3$ MeV and
$E_{e}^{\max}=8$ MeV [46], respectively. We use TEXONO II-V data with a
duration of 160 days. In the typical energy range of CsI(Tl), the Compton
scattering is the dominant process, so the probability of generating a Compton
scattering signal for each dark photon is approximated by the conversion
probability of the dark photon [29]:

$P_{\text{sig}}=\epsilon^{2}\frac{m_{A^{\prime}}^{4}}{(\Delta m^{2})^{2}}.$
(S6)

Thus, the number of signal events is obtained via

$\displaystyle N_{\text{sig}}=ST\int dE_{A^{\prime}}\frac{1}{4\pi
D^{2}}\frac{dF_{A^{\prime}}}{dE_{A^{\prime}}}P_{\text{sig}}\epsilon_{acc}(E_{A^{\prime}}),$
(S7)

where $S$ represents the cross-sectional area of the detector facing the
reactor core, $T$ is the data collection time, and $D$ is the detector-to-
reactor core distance. The differential flux for $A^{\prime}$ is $1/(4\pi
D^{2})dF_{A^{\prime}}/dE_{A^{\prime}}$, and after multiplying with the
conversion probability $P_{\rm sig}$, $A^{\prime}$ converts back into $\gamma$
without an energy change ($E_{A^{\prime}}=E_{\gamma}$) [29]. The detector
acceptance factor $\epsilon_{acc}$ reads

$\displaystyle\epsilon_{acc}(E_{\gamma})=\frac{1}{\sigma_{C}(E_{\gamma})}\int^{E_{e}^{\max}}_{E_{e}^{\min}}dE_{e}\frac{d\sigma_{C}(E_{\gamma})}{dE_{e}},$
(S8)

where $E_{\gamma}$ is the initial photon energy in Compton scattering. The
$95\%$ C.L. upper limit on the signal events from the dark photon is $1174.1$
from TEXONO data [46]. Therefore, we require $N_{\rm sig}<1174.1$ to set
limits in our scenario.

TEXONO also includes an n-type PCGe detector with a of mass is 500 g and a low
electron recoil energy threshold of $0.3\sim 12.3$ keV. The expected events
rate $N_{\text{sig}}/T$ at PCGe detector is derived from Eq. (S7) with the
cross-sectional area $S$ replaced by $N_{e}\sigma_{C}(E_{\gamma})$, where
$N_{e}$ is the electron number in the detector. This is because the size of
PCGe detector is comparable to the Compton scattering interaction length (a
few cm), and then we can no longer assume that every converted photon is
scattered in the PCGe detector. We use the reactor ON/OFF data with exposure
of $124.2/70.3$ kg-days [47, 48]. Our expected event rates take the same bins
as the TEXONO PCGe data: 60 bins with bin width 0.2 keV for each in the
electron recoil energy from 0.3 to 12.3 keV [48]. At last, we use the
minimal-$\chi^{2}$ fit method to determine the $95\%$ C.L. upper limit in our
scenario.


* [53] White, H.  A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct Test for Heteroscedasticity.  Econometrica, 48:817–838, 1980.

Online Appendix for
“Double robust inference for continuous updating GMM”


# Real time enhancement of operator’s ergonomics in physical human - robot
collaboration scenarios using a multi-stereo camera system

Gerasimos Arvanitis∗, Nikos Piperigkos†, Christos Anagnostopoulos†, Aris S.
Lalos†, Konstantinos Moustakas∗ ∗Department of Electrical and Computer
Engineering, University of Patras, Patra, Greece
†Industrial Systems Institute, ATHENA Research Center, Patra, Greece
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In collaborative tasks where humans work alongside machines, the robot’s
movements and behaviour can have a significant impact on the operator’s
safety, health, and comfort. To address this issue, we present a multi-stereo
camera system that continuously monitors the operator’s posture while they
work with the robot. This system uses a novel distributed fusion approach to
assess the operator’s posture in real-time and to help avoid uncomfortable or
unsafe positions. The system adjusts the robot’s movements and informs the
operator of any incorrect or potentially harmful postures, reducing the risk
of accidents, strain, and musculoskeletal disorders. The analysis is
personalized, taking into account the unique anthropometric characteristics of
each operator, to ensure optimal ergonomics. The results of our experiments
show that the proposed approach leads to improved human body postures and
offers a promising solution for enhancing the ergonomics of operators in
collaborative tasks.

###### Index Terms:

operator’s ergonomics, anthropometrics, human-robot collaboration, multi-
stereo camera system.

##  I Introduction

In manufacturing, when operators are physically interacting with a robot,
trying to accomplish a collaborative task, their posture is inevitably
influenced by the robot’s movement and trajectory [1, 2]. Work-related
Musculoskeletal Disorders (WMSD) is a major health problem in developed
countries [3] and there is a worldwide interest to reduce the conditions and
risk factors [4] that may cause this problem [5], decreasing in this way the
substantial costs for therapies and negative impacts on the life quality of
operators [6]. To overcome this challenging situation, adapted systems are
required, able to be adjusted according to the anthropometrics and special
characteristics of each operator, improving their physical situation and
ergonomics [7]. In literature, there are a lot of works [8] that have used
depth sensors (Kinect) or low-cost RGB devices to compute the Rapid Upper Limb
Assessment (RULA) [9, 10] score or to detect awkward postures [11]. To achieve
this, open-source landmark estimation libraries [12, 13] are used for
calculating the body joint angles. Additionally, a lot of research has been
done to improve physical ergonomics [14] in human-robot collaboration [15, 16]
by generating task assignments [17, 18], finding the optimal trade-off motions
[19], providing wearable feedback [20], estimating the physical load required
in a daily job [21], or providing an automatic assessment of ergonomics [22].

In this study, we operate under the premise that collaborative tasks between
robots and humans in industrial environments are typically comprised of pre-
defined and known actions for the sake of safety. One factor that can impact
the robot’s movement is the operator’s body characteristics, such as height.
This factor does not alter the primary action of the robot (e.g., delivering a
tool), but can be utilized to select the optimal adjustment of the robot’s
position and configuration for optimized ergonomics assessment of the operator
(e.g., placing the tool in a comfortable location without causing strain).
This adjustment takes place only at the start of the robot’s operation, as
dynamic adjustments are deemed inappropriate for safety reasons. The
adaptation of the robot’s trajectory facilitates the improvement of the
operator’s ergonomics and it is performed based on their unique anthropometric
characteristics. The key contributions of this work can be summarized as
follows:

  * •

A graph Laplacian-based scheme that provides a 3D landmark optimization,
allowing for the best use of available captured information at any time,
regardless of the camera source.

  * •

The development of a manufacturing simulator from scratch, emulating robotic
movements and operator’s actions in human-robot collaboration scenarios.

  * •

A framework that increases the operators’ situational awareness regarding
uncomfortable work postures that may be harmful in long-term exposure. This is
achieved through a personalized anthropometric estimation component that
provides a comprehensive ergonomic analysis of the operators’ actions during
their collaboration with the robot.

  * •

An offline tool, suitable for the ergonomics experts to analyse and evaluate
the operator’s action and help them to set personalized rehabilitation
strategies or recommendations to avoid muscle injuries.

The rest of this paper is organized as follows: In Section 2, we discuss in
detail each step of the proposed method. Section 3 presents the experimental
results and in Section 4 we draw the conclusions.

##  II General Architecture & Methodology

The architecture of the proposed multi-stereo camera system is illustrated in
Fig. 1. In our implementation, the system consists of three stereo cameras
placed in strategic locations to provide continuous monitoring of the
operator’s pose from multiple angles, increasing the possibility of
continuously having a good capture of the operator’s pose (i.e., a view in
front of the operator), at least from one of the cameras, while he/she is
freely moving in different directions.

Figure 1: Proposed concept architecture.

We use each RGB camera of each stereo set to capture the operator’s actions.
To extract 2D posture landmarks, we utilize the MediaPipe framework [23],
which runs locally on the Jetson device where the stereo camera is connected.
The extracted landmarks are then sent to a PC server for 3D landmark
estimation using a Direct Linear Transform (DTL) triangulation approach [24]
and optimization. Next, we estimate the operator’s height and calculate their
current physical ergonomic state based on the RULA score111https://ergo-
plus.com/rula-assessment-tool-guide. For seamless integration of the system
components and algorithms, we utilize the Robotic Operating System (ROS) as
our middleware. The system runs on isolated nodes that communicate using a
publish-subscribe model, and we have developed a ROS node for each camera.
Each node captures a frame of the camera at a frequency of 10Hz, extracts and
processes the landmarks, and publishes them to a relevant ROS topic. The
decentralized approach of ROS allows us to connect the cameras to different
physical machines, such as the Jetson TX2 embedded device, and all components
must be connected to the same local network.

Algorithm 1 summarizes the main steps of the proposed framework. More details
will be presented in the next sections.

/* A) Calibration steps (It is performed only once) */

1 for  _$i=1:M$_ do

2 Single camera $C_{i}$ calibration;

3 end for

4for  _$j=1:L$_ do

5 Stereo camera $S_{j}$ calibration;

6 end for

Data: Frames’ sequency from the $M$ cameras

Result: Improve operators’ physical ergonomics

/* B) 3D landmarks estimation */

7 for  _$i=1:M$_ do

/* It runs locally to the Jetson devices */

8 2D pose landmark estimation from $C_{i}$;

9 end for

10Send the 2D pose landmarks to the central server via ROS;

11 for  _$j=1:L$_ do

12 3D pose landmark estimation from $S_{j}$ via DLT

13 end for

/* C) 3D landmarks optimization */

14 Follow the steps 1-5 of section II.B;

/* D) Anthropometric and ergonomic analysis */

15 Operator’s height estimation;

16 Real-time RULA score analysis;

Algorithm 1 Algorithmic steps of the framework

Figure 2: Examples of the simulator in different configurations (height of the
operator), corresponding to 3 different anthropometric classes.

###  II-A 3D Landmarks Calculation

1) Single Camera Calibration. To calibrate each camera individually, we begin
by using the same checkerboard pattern positioned in a location that is
visible to all cameras at the same time. This step is crucial for performing
stereo calibration later. We also ensure that all the used cameras have the
same resolution analysis and capture images at the same frame-per-second (fps)
ratio to ensure that all captured frames are synchronized [25].

2) Stereo Calibration. Let’s assume that each camera of the integrated system
can be defined as $C_{i}\ (\forall\ i=1,\cdots,6)$ with
$\textbf{c}_{i}\in\mathbb{R}^{3\times 1}$ 3D coordinates. The single-camera
calibration of camera $C_{i}$ returns the rotation matrix
$\textbf{R}_{i}\in\mathbb{R}^{3\times 3}$, the translation
$\textbf{T}_{i}\in\mathbb{R}^{3\times 1}$ matrix and the distortion
coefficients. However, the R and T matrices alone are not enough to
triangulate the 3D points. For the stereo-calibration procedure, we need two
sets of cameras. For the sake of simplicity, let’s define these two cameras as
$C_{1}$ and $C_{2}$. Notice here that the same approach can be used for the
stereo calibration of all other cameras. The rotation matrix R and translation
vector T shows how to go from the $C_{2}$ coordinate system to the $C_{1}$. We
also obtain world coordinates to $C_{2}$ rotation and translation, by
calculating: $\textbf{R}_{2}=\textbf{R}\textbf{R}_{1}$,
$\textbf{T}_{2}=\textbf{R}\textbf{T}_{1}+\textbf{T}$. Next, we choose the
$C_{1}$ position as world space origin ($x=0,y=0,z=0$). Thus, the world origin
to $C_{1}$ rotation is equal to the identity matrix and translation is a zeros
vector. Then, the R matrix and the T vector that have been estimated from the
previous stereo-calibration step becomes the rotation and translation from
world origin to $C_{2}$. Basically, this means that the 3D triangulated points
will be in a space with respect to the coordinate system of $C_{1}$’s lens. In
other words, we simply overlap world coordinates with the coordinates of the
$C_{1}$ camera. This means $\textbf{R}_{1}=\text{eye}(3)$,
$\textbf{T}_{1}=\text{zeros}(3)$ and $\textbf{R}_{2}=\textbf{R}$,
$\textbf{T}_{2}=\textbf{T}$. Therefore, all triangulated 3D points are
measured from the $C_{1}$ camera position in the world.

3) Direct Linear Transform. We assume the existence of a 3D point p in real
space with coordinates given as $\textbf{p}=[x,y,z]$. This point can be
observed from both 2 cameras, which have pixel coordinates
$\textbf{G}_{i}=[u_{i},v_{i},1]$ for $C_{i},\forall\ i=1,2$. Using the camera
projection matrix
$\textbf{K}_{1}=[\textbf{R}_{1}\textbf{T}_{1}]\in\mathbb{R}^{3\times 4}$,
$\textbf{G}_{1}$ can be written as: ${\rm\textbf{G}}_{1}=a\textbf{K}_{1}$. In
a triangulation problem, we do not know the coordinates of p. But we can
determine the pixel coordinates and the projection matrix through camera
calibration. Since $\textbf{G}_{1}$ and $\textbf{K}_{1}\textbf{p}$ are
parallel vectors, the cross product of these should be zero. The $j$ row
vector of $\textbf{K}_{1}$ can be written as $\textbf{k}_{j}\forall\
j=1,...,4$. This gives us:

$\begin{aligned} \noindent\begin{bmatrix}u_{1}\\\ v_{1}\\\
1\end{bmatrix}\times\begin{bmatrix}{\textbf{k}}_{1}{\textbf{p}}\\\
{\textbf{k}}_{2}{\textbf{p}}\\\
{\textbf{k}}_{3}{\textbf{p}}\end{bmatrix}=\begin{bmatrix}v_{1}{\textbf{k}}_{3}{\textbf{p}}-{\textbf{k}}_{2}{\textbf{p}}\\\
{\textbf{k}}_{3}{\textbf{p}}-u_{1}{\textbf{k}}_{3}{\textbf{p}}\\\
u_{1}{\textbf{k}}_{2}{\textbf{p}}-v_{1}{\textbf{k}}_{1}{\textbf{p}}\end{bmatrix}=\\\
\begin{bmatrix}v_{1}{\textbf{k}}_{3}-{\textbf{k}}_{2}\\\
{\textbf{k}}_{3}-u_{1}{\textbf{k}}_{3}\\\
u_{1}{\textbf{k}}_{2}-v_{1}{\textbf{k}}_{1}\end{bmatrix}{\textbf{p}}=\begin{bmatrix}0\\\
0\\\ 0\end{bmatrix}\end{aligned}$ (1)

This gives us an equation of the form $\textbf{A}\textbf{x}=0$. But the third
row is a linear combination of the first two rows giving 2 systems of
equations, which is not enough to solve the 3 unknowns in p. Since we have two
cameras, we can extend the matrix to have more rows. In fact, we simply add on
more rows for any number of views. This gives us the equation:

$\begin{bmatrix}v_{1}{\textbf{k}}_{3}-{\textbf{k}}_{2}\\\
{\textbf{k}}_{3}-u_{1}{\textbf{k}}_{3}\\\
u_{1}{\textbf{k}}_{2}-v_{1}{\textbf{k}}_{1}\\\
v_{2}{\textbf{k}}_{3}-{\textbf{k}}_{2}\\\
{\textbf{k}}_{3}-u_{2}{\textbf{k}}_{3}\\\
u_{2}{\textbf{k}}_{2}-v_{2}{\textbf{k}}_{1}\end{bmatrix}{\textbf{p}}=0$ (2)

DLT is a method for calculating a matrix equation of the form
$\textbf{A}\textbf{p}=0$. In the real world, there can be some noise, so we
write the equation as $\textbf{A}\textbf{p}=\textbf{n}$, and we solve for p
such that n is minimized. The first step is to determine the SVD decomposition
of A.

$\textbf{A}{\textbf{p}}=\textbf{U}\textbf{S}\textbf{V}^{T}{\textbf{p}}$ (3)

Our goal is to minimize n for some p. This can be done by taking the dot
product:

${\textbf{n}}^{T}{\textbf{n}}=({\textbf{p}}^{T}\textbf{V}\textbf{S}\textbf{U}^{T})\cdot(\textbf{U}\textbf{S}\textbf{V}^{T}{\textbf{p}})={\textbf{p}}^{T}\textbf{V}\textbf{S}^{2}\textbf{V}^{T}{\textbf{p}}$
(4)

U and V are orthonormal matrices and S is a diagonal matrix. The entries on
the diagonal of S are decreasing, so that the last entry on the diagonal is
the minimum value. Since our goal was to minimize $\textbf{n}^{T}\textbf{n}$,
this tells us that it is equivalent to choosing the smallest value of
$\textbf{S}^{2}$ by selecting the corresponding column vector of
$\textbf{V}^{T}$ as p.

###  II-B 3D Landmark Optimization Algorithm

As a result of varying camera viewpoints and operator movements, we expect
that landmarks detected from different cameras will not be the same and may
differ in accuracy. Therefore, it is necessary to design a fusion approach
that integrates the multi-camera system to provide the best possible landmarks
configuration from all the involved cameras.

Figure 3: Indicative example of graph topology for cameras and landmarks.

For this task, we create an undirected graph of $N+M$ points
$\hat{\textbf{p}}=[\textbf{p}\ \textbf{c}]^{T}\in\mathbb{R}^{(N+M)\times 3}$,
consisting of the $N$ detected landmarks as edges p, and the $M$ cameras as
nodes c. Fig 3 presents a simplified example of a graph with $M=4$ and $N=4$.
In order to take advantage of this graph topology formulation, we apply the
Graph Laplacian Processing (GLP) technique [26] for re-estimating landmarks in
an optimal manner, following the next steps:

1) Encode the spatial relationship of cameras $C$ and landmarks p via GLP.

2) Construct the binary Laplacian matrix
$\bm{L}\in\mathbb{R}^{(N+M)\times(N+M)}$ of connectivity graph, indicating
which landmark is related to which camera. To mention here that there are not
connections between two cameras and between two landmarks.

3) Use as anchor points $\textbf{a}=[\bar{\textbf{p}}\
\textbf{c}]^{T}\in\mathbb{R}^{(N+M)\times 3}$ the known location of cameras c
and the average position of landmarks
$\bar{\textbf{p}}^{i}=\sum_{i=1}^{M}\textbf{p}^{i}\in\mathbb{R}^{N\times 3}$
as estimated by each one of the cameras.

4) Relative measurements and anchors are used to formulate differential
coordinates $\bm{\delta}\in\mathbb{R}^{(N+M)\times 3}$ [27, 28], where the
$i^{th}$ row of $\bm{\delta}_{i}$ is equal to:

$\begin{aligned} \noindent\bm{\delta}_{i}=[\bm{\delta}_{xi}\ \bm{\delta}_{yi}\
\bm{\delta}_{zi}]=\hat{\textbf{p}}_{i}-\frac{1}{|\Psi_{i}|}\sum_{j\in\Psi_{i}}\hat{\textbf{p}}_{j}=\bm{L}\hat{\textbf{p}}_{i},\\\
\ \ \forall\ i=1,...,N+M\end{aligned}$ (5)

where $\bm{l_{i}}$ is the $i$-row of Laplacian matrix, and $|\Psi_{i}|$ is the
number of immediate neighbors of $\hat{\textbf{p}}_{i}$.

5) The optimized positions of point $\hat{\textbf{p}}$ are given by the
minimization of the following cost function:

$\operatorname*{argmin}\limits_{\hat{\textbf{p}}}{\left\lVert\bm{\tilde{L}}\hat{\textbf{p}}-\bm{b}\right\rVert_{2}^{2}}$
(6)

where $\bm{\tilde{L}}=[\bm{L}\
\mathbb{I}_{N+M}]^{T}\in\mathbb{R}^{2(N+M)\times(N+M)}$, $\mathbb{I}_{N+M}$ is
the identity matrix and $\bm{b}=[\bm{\delta}\
\textbf{a}]^{T}\in\mathbb{R}^{2(N+M)\times 3}$. Finally, the re-estimated
positions of landmarks are derived from linear least-squares minimization:

$\hat{\textbf{p}}=\left(\bm{L^{T}L}\right)^{-1}\bm{L^{T}b}$ (7)

To note that during the simulation horizon the topology of landmarks and
cameras remains unchanged. As such, the term
$\left(\bm{L^{T}L}\right)^{-1}\bm{L^{T}}$ is constant and can be computed only
once at the start of simulation, in order to be used for the successive
frames. In that way, we reduce the computational time required for solving (7)
since the costly operation of matrix inversion is performed only at the
beginning of the process.

##  III Experimental Analysis and Results

###  III-A Experimental Setup and Implementations

The operating system of the server is UBUNTU 20.04 which is compatible with
ROS 1. All the algorithms are written in Python 3.10. There are no special
requirements regarding the computational power of the server since the
implementation can easily run to low-cost devices (e.g., Jetson).

1) Operator’s Height Estimation. For this use case two steps are followed: i)
Estimation of the operator’s height; ii) Adaptation of the robot’s position.
Once the operator’s height has been estimated, the robot’s control system
receives the information and configures the robot’s movement parameters
according to the selected scenario. This adaptation ensures that the robot’s
movements are ergonomically comfortable for the interacting user. 3 different
classes $c$ are used, based on the operators’ height ($c_{1}<1.68,1.68\leq
c_{2}\leq 1.82,c_{3}>1.82$) (Fig. 2). Each class leads to an adaptable robot
response, ensuring that the user can work in an optimal ergonomic state.

2) Real-Time RULA Estimation. We perform a real-time joint angle estimation to
calculate the RULA score (ranging from 1 to 7, used by ergonomists to evaluate
work activities involving upper-body motion [29]). Based on this value, the
appropriate warning messages are sent to the operators, informing them if
their working pose is ergonomically correct or not (Fig. 4), so as to correct
their posture in real-time.

Figure 4:  Safe and wrong pose identification.

3) Simulator. To evaluate the effectiveness of our framework, we have
developed an in-house simulation framework from scratch (Fig. 5). The
simulator emulates a robotic movement and the corresponding operator’s actions
in a digital workspace created in Unity3D [30]. Fig. 7 presents an example of
a user’s action before and after the robot’s adaptation. 3 stereo cameras are
involved, located in 3 different areas of the working environment, monitoring
the operator’s actions while performs a collaborative tasks with a KUKA robot.
The primary objective of the manufacturing simulator is to create a digital
twin, which provides a safe environment for evaluating the developed framework
and algorithms. Additionally, the simulator can be also used for further
ergonomic analysis in order to optimize the workstation of a real industrial
environment, so as to obtain a trajectory of the robot as much as possible
adapted to the human. Furthermore, the simulator serves as a validation
framework since the ground truth information of the 3D landmarks’ locations is
not known in the real case.

Figure 5: Screenshot from the environment of the simulator.

4) Ergonomics Evaluation Tool. The 3D landmarks are stored and are accessible
offline (i.e., 3D view, pause, zoom in/out) by ergonomic experts that can
evaluate the poses during a specific action, and suggest how it could be
improved. Fig. 6 illustrates an example showing how the joint angles are
strained during an action. A heatmap visualizes the level of strength (deep
blue is the lowest, deep red is the highest).

Figure 6:  Ergonomics tool for the evaluation of the joint angles stress. TABLE I:  RMSE of landmarks using multicamera system and GLP. | $\textbf{p}_{1}$ | $\textbf{p}_{2}$ | $\textbf{p}_{3}$ | $\textbf{p}_{4}$ | $\textbf{p}_{5}$ | $\textbf{p}_{6}$ | $\textbf{p}_{7}$ | $\textbf{p}_{8}$ | $\textbf{p}_{9}$ | $\textbf{p}_{10}$ | $\textbf{p}_{11}$ | $\textbf{p}_{12}$
---|---|---|---|---|---|---|---|---|---|---|---|---
$S_{1}$ | 0.34 | 0.38 | 0.19 | 0.24 | 0.15 | 0.18 | 0.24 | 0.27 | 0.27 | 0.27 | 0.22 | 0.35
$S_{2}$ | 0.36 | 0.40 | 0.19 | 0.28 | 0.15 | 0.21 | 0.27 | 0.29 | 0.28 | 0.32 | 0.24 | 0.38
$S_{3}$ | 0.56 | 0.57 | 0.34 | 0.35 | 0.32 | 0.30 | 0.41 | 0.42 | 0.43 | 0.40 | 0.31 | 0.43
fusion | 0.38 | 0.39 | 0.20 | 0.24 | 0.16 | 0.16 | 0.27 | 0.27 | 0.30 | 0.29 | 0.21 | 0.35
Figure 7: Screenshot from the environment of the simulator before and after
robot’s adaptation.

###  III-B Experimental Results

Table I summarizes the main outcomes (RMSE values of the 12 landmarks p, per
each stereo camera $S_{i}\ \forall\ i\in(1-3)$) from a testing scenario
running in the simulator. As we can see, in most of the cases, the proposed
GF-based fusion approach achieved results very close to the optimal, while in
two cases ($\textbf{p}_{6}$ and $\textbf{p}_{11}$), landmarks were detected
slightly better than the optimal camera. To mention here that the evaluation
of the fusion algorithm’s performance can be achieved only in a simulated
environment since in a real-case scenario the ground truth locations of the
landmarks, as well as the accuracy of the cameras, can not be known.

Fig 8 shows a plot with the mean RULA score of different users with different
heights (1.5-2 m.) that perform a very specific task before and after the
robot’s adjustment. As we can see, only operators with 1.75 m. height has a
low RULA score before the adaptation, while it is improved for all of the
users after the robot’s adaptation.

Figure 8:  RULA score for operators with different heights.

Fig 9 shows the mean score per different body areas for all operators. The
improvement is more apparent for the cases of ”Neck” and ”Lower arms”. Fig 10
shows how the distribution of the angles per each joint decreases after the
robot’s adaptation.

Figure 9:  Mean score per different body areas. Figure 10:  Distribution of
angles per joint.

##  IV Conclusions & Future Work

In this work, we introduced a multi-stereo camera pipeline aimed at enhancing
the operator’s physical ergonomics by optimizing the estimation of 3D
landmarks. This is achieved in two ways: 1) by adjusting the movement
parameters of the collaborative machine to the unique anthropometrics of the
operator and 2) by modifying the operator’s motion and actions through real-
time visual warnings or expert recommendations, based on RULA score
estimation. To facilitate the evaluation of our solution, we also developed a
simulator. Future plans include implementation in a real industrial setting
and a survey to gauge the effectiveness of the method through real user
feedback and additionally the evaluation in more complex collaborative tasks.

## References

  * [1] Ilias El Makrini, Kelly Merckaert, Joris De Winter, Dirk Lefeber, and Bram Vanderborght,  “Task allocation for improved ergonomics in human-robot collaborative assembly,”  Interaction Studies, vol. 20, no. 1, pp. 102–133, 2019.
  * [2] Efstratios Kontopoulos, Gerasimos Arvanitis, Alessandro Zanella, Panagiotis Mitzias, Evangelia I. Zacharaki, Pavlos Kosmides, Nikos Piperigkos, Konstantinos Moustakas, and Aris S. Lalos,  “Semantic data integration for monitoring operators’ ergonomics in an automotive manufacturing setting,”  in The Semantic Web: ESWC 2022 Satellite Events, Paul Groth, Anisa Rula, Jodi Schneider, Ilaria Tiddi, Elena Simperl, Panos Alexopoulos, Rinke Hoekstra, Mehwish Alam, Anastasia Dimou, and Minna Tamper, Eds., Cham, 2022, pp. 303–306, Springer International Publishing.
  * [3] Pauline Maurice, Vincent Padois, Yvan Measson, and Philippe Bidaud,  “Human-oriented design of collaborative robots,”  International Journal of Industrial Ergonomics, vol. 57, pp. 88–102, 2017.
  * [4] Soheil Gholami, Marta Lorenzini, Elena De Momi, and Arash Ajoudani,  “Quantitative physical ergonomics assessment of teleoperation interfaces,”  IEEE Transactions on Human-Machine Systems, vol. 52, no. 2, pp. 169–180, 2022.
  * [5] André Cardoso, Ana Colim, Estela Bicho, Ana Cristina Braga, Marino Menozzi, and Pedro Arezes,  “Ergonomics and human factors as a requirement to implement safer collaborative robotic workstations: A literature review,”  Safety, vol. 7, no. 4, 2021.
  * [6] Marta Lorenzini, Wansoo Kim, Elena De Momi, and Arash Ajoudani,  “A new overloading fatigue model for ergonomic risk assessment with application to human-robot collaboration,”  in 2019 International Conference on Robotics and Automation (ICRA), 2019, pp. 1962–1968.
  * [7] Sujee Lee, Li Liu, Robert Radwin, and Jingshan Li,  “Machine learning in manufacturing ergonomics: Recent advances, challenges, and opportunities,”  IEEE Robotics and Automation Letters, vol. 6, no. 3, pp. 5745–5752, 2021.
  * [8] Ernesto Rocha-Ibarra, Marvella-Izamar Oros-Flores, Dora-Luz Almanza-Ojeda, Gabriel-Armando Lugo-Bustillo, Andres Rosales-Castellanos, Mario-Alberto Ibarra-Manzano, and Juan Carlos Gomez,  “Kinect validation of ergonomics in human pick and place activities through lateral automatic posture detection,”  IEEE Access, vol. 9, pp. 109067–109079, 2021.
  * [9] Pierre Plantard, Hubert P. H. Shum, Anne-Sophie Le Pierres, and F. Multon,  “Validation of an ergonomic assessment method using kinect data in real workplace conditions.,”  Applied ergonomics, vol. 65, pp. 562–569, 2017.
  * [10] Darius Nahavandi and Mohammed Hossny,  “Skeleton-free rula ergonomic assessment using kinect sensors,”  Int. Dec. Tech., vol. 11, no. 3, pp. 275–284, jan 2017.
  * [11] Vito Modesto Manghisi, Antonio Emmanuele Uva, Michele Fiorentino, Vitoantonio Bevilacqua, Gianpaolo Francesco Trotta, and Giuseppe Monno,  “Real time rula assessment using kinect v2 sensor,”  Applied Ergonomics, vol. 65, pp. 481–491, 2017.
  * [12] Z. Cao, G. Hidalgo, T. Simon, S. Wei, and Y. Sheikh,  “Openpose: Realtime multi-person 2d pose estimation using part affinity fields,”  IEEE Transactions on Pattern Analysis Machine Intelligence, vol. 43, no. 01, pp. 172–186, jan 2021.
  * [13] Woojoo Kim, Jaeho Sung, Daniel Saakes, Chunxi Huang, and Shuping Xiong,  “Ergonomic postural assessment using a new open-source human pose estimation technology (openpose),”  International Journal of Industrial Ergonomics, vol. 84, pp. 103164, 2021.
  * [14] Ilias El Makrini, Glenn Mathijssen, Sten Verhaegen, Tom Verstraten, and Bram Vanderborght,  “A virtual element-based postural optimization method for improved ergonomics during human-robot collaboration,”  IEEE Transactions on Automation Science and Engineering, vol. 19, no. 3, pp. 1772–1783, 2022.
  * [15] Marianna Ciccarelli, Alessandra Papetti, Cecilia Scoccia, Giacomo Menchi, Leonardo Mostarda, Giacomo Palmieri, and Michele Germani,  “A system to improve the physical ergonomics in human-robot collaboration,”  Procedia Computer Science, vol. 200, pp. 689–698, 2022,  3rd International Conference on Industry 4.0 and Smart Manufacturing.
  * [16] Amir Yazdani, Roya Sabbagh Novin, Andrew Merryweather, and Tucker Hermans,  “Ergonomically intelligent physical human-robot interaction: Postural estimation, assessment, and optimization,” 2021.
  * [17] Margaret Pearce, Bilge Mutlu, Julie Shah, and Robert Radwin,  “Optimizing makespan and ergonomics in integrating collaborative robots into manufacturing processes,”  IEEE Transactions on Automation Science and Engineering, vol. 15, no. 4, pp. 1772–1784, 2018.
  * [18] Baptiste Busch, Guilherme Maeda, Yoan Mollard, Marie Demangeat, and Manuel Lopes,  “Postural optimization for an ergonomic human-robot interaction,”  in 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2017, pp. 2778–2785.
  * [19] Waldez Gomes, Pauline Maurice, Eloïse Dalin, Jean-Baptiste Mouret, and Serena Ivaldi,  “Multi-objective trajectory optimization to improve ergonomics in human motion,”  IEEE Robotics and Automation Letters, vol. 7, no. 1, pp. 342–349, 2022.
  * [20] Wansoo Kim, Marta Lorenzini, Kağan Kapıcıoğlu, and Arash Ajoudani,  “Ergotac: A tactile feedback interface for improving human ergonomics in workplaces,”  IEEE Robotics and Automation Letters, vol. 3, no. 4, pp. 4179–4186, 2018.
  * [21] Marta Lorenzini, Wansoo Kim, and Arash Ajoudani,  “An online multi-index approach to human ergonomics assessment in the workplace,”  IEEE Transactions on Human-Machine Systems, vol. 52, no. 5, pp. 812–823, 2022.
  * [22] Adrien Malaisé, Pauline Maurice, Francis Colas, and Serena Ivaldi,  “Activity recognition for ergonomics assessment of industrial tasks with automatic feature selection,”  IEEE Robotics and Automation Letters, vol. 4, no. 2, pp. 1132–1139, 2019.
  * [23] Charles Andrew Q. Bugarin, Juan Miguel M. Lopez, Scud Gabriel M. Pineda, Ma. Franzeska C. Sambrano, and Pocholo James M. Loresco,  “Machine vision-based fall detection system using mediapipe pose with iot monitoring and alarm,”  in 2022 IEEE 10th Region 10 Humanitarian Technology Conference (R10-HTC), 2022, pp. 269–274.
  * [24] Y.I Abdel-Aziz, H.M. Karara, and Michael Hauck,  “Direct linear transformation from comparator coordinates into object space coordinates in close-range photogrammetry*,”  Photogrammetric Engineering and Remote Sensing, vol. 81, no. 2, pp. 103–107, 2015.
  * [25] Vaishnav Chunduru, Mrinalkanti Roy, Dasari Romit N. S, and Rajeevlochana G. Chittawadigi,  “Hand tracking in 3d space using mediapipe and pnp method for intuitive control of virtual globe,”  in 2021 IEEE 9th Region 10 Humanitarian Technology Conference (R10-HTC), 2021, pp. 1–6.
  * [26] Nikos Piperigkos, Aris S. Lalos, Kostas Berberidis, and Christos Anagnostopoulos,  “Cooperative multi-modal localization in connected and autonomous vehicles,”  in 2020 IEEE 3rd Connected and Automated Vehicles Symposium (CAVS), 2020, pp. 1–5.
  * [27] Olga Sorkine-Hornung,  “Laplacian mesh processing,”  in Eurographics, 2005.
  * [28] Gerasimos Arvanitis, Aris S. Lalos, and Konstantinos Moustakas,  “Fast spatio-temporal compression of dynamic 3d meshes,”  in 2021 IEEE 23rd International Workshop on Multimedia Signal Processing (MMSP), 2021, pp. 1–6.
  * [29] Lynn McAtamney and E. Nigel Corlett,  “Rula: a survey method for the investigation of work-related upper limb disorders,”  Applied Ergonomics, vol. 24, no. 2, pp. 91–99, 1993.
  * [30] Pavlou M., Laskos D., Zacharaki E., Risvas K., and Moustakas K.,  “Xrsise: An xr training system for interactive simulation and ergonomics assessment,”  Frontiers in Virtual Reality, vol. 2, pp. 17, 2021.