\section{Introduction}\label{intro}

56 Ursae Majoris (also known as 56~UMa, HD\,98839, or HIP\,55560) is a binary system that contains a G8-type giant and a faint lower-mass companion \citep{Keenan89}. The giant was identified as a mild barium (Ba) star by \cite{Keenan77}. Ba stars are chemically peculiar stars that show an overabundance of carbon and elements formed by the slow-neutron-capture (s-) process of nucleosynthesis \citep{BidelmanKeenan51}, which mostly operates in the interiors of Asymptotic Giant Branch (AGB) stars \citep[e.g.][]{Lugaro03,Cristallo09,Karakas10}. Ba stars form when an AGB star transfers enriched material to a companion in a binary system \citep{McClure80}. According to this scenario, the donors in these systems evolved off the AGB long ago and are now dim white dwarfs (WDs).

Even though the Ba-star classification of the giant star in 56~UMa has been debated \citep[e.g.][]{Keenan89, Lu91, Griffin1996, Karinkuzhi18}, and its high luminosity and mass made it an outlier among Ba giants \citep{Jorissen19}, the nature of its companion has not been extensively discussed. \cite{Bohm-Vitense84} noticed that the International Ultraviolet Explorer (IUE) spectrum of 56~UMa showed two weak emission lines at 147.4 and 166.8~nm, likely due to the N IV] and O III] semi-forbidden lines. They proposed that these lines arise from an extended region of high-temperature, low-density gas coming from a small, hot source in the system, and suggested that this source might be a massive WD. However, it is interesting to note that the same lines arise in the hot gas surrounding accreting neutron stars \citep[NSs;][]{Anderson1994,CastroSegura2022}.

\begin{figure*}[t]
\centering
\includegraphics[width=0.9\textwidth]{HD98839.png}
\caption{Fits to the RV curve (left) and the \textit{Hipparcos} and \textit{Gaia} proper motions (middle and right) of 56~UMa. The black line is the orbit with the maximum likelihood, and 40 additional well-fitting solutions are included, colour-coded as a function of companion mass.}\label{fits}
\end{figure*}

At the location of 56~UMa, there is also a strong detection in the GALEX (Galaxy Evolution Explorer) far-UV band (153~nm), at the level of $7.66 \times 10^{-5}$~Jy (GALEX-DR5; \citealt{Bianchi2011}). This is comparable with the $6.09 \times 10^{-6}$~Jy measured at 212~nm (XMM-SUSS5.0; \citealt{Page2012}) for the closest (120~pc) NS to Earth, RX J1856.5-3754 \citep{Walter2010}. 56~UMa is located at 170~pc \citep{GaiaEDR3}. Finally, in a study of coronal X-ray emission in late-type giants, \cite{Gondoin1999} reports an X-ray luminosity of $L_X = 1.07 \times 10^{29}$~erg~s$^{-1}$ for 56~UMa. This could be ascribed to coronal emission or to X-ray emission associated with a NS. Neutron stars are observed to have masses between 1.2 and 2.35~M$_{\odot}$\ \citep{Romani22} and typically emit at radio and X-ray wavelengths when they are young and rapidly rotating. As they age, their rotation slows down and their emission fades, unless they are fuelled by mass transfer from a companion. Neutron stars not powered by accretion are difficult to detect and, therefore, very rarely observed \citep{Mazeh2022}.

Here, we revisit the stellar and orbital parameters of 56~UMa (Sect. \ref{orbit}), the stellar properties of the giant component (Sect. \ref{m1}), and its heavy-element abundances (Sect. \ref{abundances}). We show possible evidence of a supernova explosion in the system (Sect. \ref{sec:cavity}), and we discuss different evolutionary scenarios compatible with these observations (Sect. \ref{sec:discussion}), including the possibility that the faint component of 56~UMa is a NS.

\section{Orbital parameters and companion mass}\label{orbit}

In order to obtain the orbital parameters of 56~UMa, we fit a Keplerian orbital model to the available radial-velocity (RV) and astrometric data. We used CAMCOR (Coravel
at the Cambridge 36-inch telescope; \citep{Griffin08}) and HERMES (High-Efficiency and high-Resolution Mercator Echelle Spectrograph; \citealt{Raskin11,Jorissen19}) RV data, the positions and proper motions provided by the \textit{Hipparcos} catalogue \citep{ESA1997, HippCat97} and the \textit{Gaia} Early Data Release 3 \citep[EDR3;][]{GaiaEDR3summary, GaiaEDR3}, the individual astrometric measurements from the re-reduction of the \textit{Hipparcos} intermediate astrometric data \citep{vanLeeuwen07} and the \textit{Hipparcos}-\textit{Gaia} Catalogue of Accelerations (HGCA; \citealt{Brandt18, HGCA-EDR3}). We used the code \textsc{orvara}, developed by \cite{Brandt21}, to exploit the combinations of these datasets. \textsc{orvara} fits a single Keplerian model to the datasets described above, employing a parallel-tempering Markov chain Monte Carlo minimisation routine \citep[\textsc{ptmcmc};][]{PTMCMC}. We fit ten parameters: the six Keplerian orbital elements, the masses of the two components and a RV jitter per instrument to be added to the RV uncertainties. We assumed uninformative priors for the orbital elements and the faint companion mass, but we adopted a Gaussian prior of $4.3\pm 0.6$~M$_{\odot}$\ \citep{Jorissen19} for the primary mass (we used a 3$\sigma$ uncertainty to account for systematic errors in the mass determination) and another Gaussian prior for the parallax that corresponded to the EDR3 value. The results are presented in Table \ref{table:MCMCresults}.

\renewcommand{\arraystretch}{1.2}
\begin{table}[t]
\begin{center}
\caption{Stellar and orbital parameters of the 56~UMa system}\label{table:MCMCresults}
\vspace{1mm}
\begin{small}
\begin{tabular}{l c}
\hline
{\bf Parameter} & {\bf Median $\pm$ 1$\sigma$}\\
\hline
Temperature, $T_{\rm eff}$ [K] & 4917 $\pm$ 34\\
Surface gravity, log \textit{g} & 2.3 $\pm$ 0.6\\
Metallicity, [Fe/H] & $-0.05$\\
Microturbulence, $\xi$ [km/s] & 1.56\\
Primary mass, $M_{\rm 1}$ [$M_{\odot}$] & 4.3 $\pm$ 0.2 \\
\hline
Period, $P$ [days] & 16911$^{+438}_{-401}$\\
Eccentricity, $e$ & 0.562$^{+0.012}_{-0.012}$\\
Semi-major axis, $a$ [AU] & 22.9$^{+1.0}_{-1.1}$\\
Argument of periastron, $\omega_{\rm 1}$ [$^{\circ}$] & 286$^{+2.3}_{-2.3}$\\
Time of periastron, $T_{0}$ [HJD] & 2,468,401$^{+432}_{-385}$\\
Parallax, $\varpi$ [mas] & 5.86$^{+0.03}_{-0.04}$ \\
Ascending node, $\Omega$ [$^{\circ}$] & 60$^{+3}_{-3}$\\
Inclination, $i$ [$^{\circ}$] & 68$^{+3.6}_{-3.4}$ \\
Secondary mass, $M_{\rm 2}$ [$M_{\odot}$] & 1.31$^{+0.11}_{-0.12}$\\
Center-of-mass velocity [km/s] & $0.13\pm 0.01$\\
Center-of-mass $\mu_{\alpha*}$ [mas/yr] & $-37.32\pm 0.01$\\
Center-of-mass $\mu_{\delta}$ [mas/yr] & $-12.18\pm 0.01$\\
\hline
\end{tabular}
\end{small}
\end{center}
\end{table}

\begin{figure*}[t]
\centering
\includegraphics[width=0.7\textwidth]{hrd.png}
\caption{Location of 56~UMa in the HRD. On the left, we included STAREVOL solar-metallicity tracks, and on the left, we compared STAREVOL (solid lines), Geneva (dot-dashed lines) and Padova (dashed lines) tracks for 3, 4 and 5~M$_\odot$ stars.}\label{hrd}
\end{figure*}

Figure \ref{fits} shows the fit to the RV data and to the \textit{Hipparcos} and \textit{Gaia} proper motions, and Fig. \ref{fig:corner} shows a corner plot with some derived parameters. The system is in an orbit of about 46 years and has a significant eccentricity of 0.56. Our joint RV and astrometric fit yields orbital parameters that are compatible with the results published by \cite{Jorissen19}. Additionally, for the first time, we constrained the inclination to $67.5 \pm 4^{\circ}$ and the secondary mass to $1.31 \pm 0.12$~M$_{\odot}$. This mass estimate, which is compatible with the value estimated by \cite{Kervella19} from the proper motion anomaly only, raises three issues: (i) the companion mass is surprisingly high for WDs and close to the Chandrasekhar mass \citep{Chandrasekhar1939}; (ii) in the context of Ba stars, this WD would be significantly more massive than any other Ba-star companion (see \citealt{Escorza23}); and (iii) the AGB progenitor of such a WD would be at the upper end of intermediate-mass stars, where the s-process is predicted to be fairly inefficient \citep{Kappeler11}. Alternatively, this high mass is compatible with the mass of a NS.

\section{Constraining the primary mass}\label{m1}

Given the strong correlation between the primary and secondary masses (see Fig. \ref{fig:corner}) and the fact that we used a prior for the primary mass, we re-evaluated this value. To do this, we co-added 27 high-resolution (R $= 86\,000$) HERMES \citep{Raskin11} spectra to first recalculate the stellar parameters of the giant.
The signal-to-noise ratio (S/N) of the co-added spectrum was about 500, and the stellar parameters were derived using the Brussels Automatic Code for Characterizing High accUracy Spectra (BACCHUS) in a semi-automated mode \citep{Masseron2016}. BACCHUS uses interpolated MARCS model atmospheres \citep{Gustafsson2008} and the 1D local thermodynamical equilibrium (LTE) spectrum-synthesis code TURBOSPECTRUM \citep{Plez2012}. It computes abundances using equivalent widths or spectral synthesis, allowing excitation and ionization equilibria to be checked for, thereby constraining $T_{\rm eff}$ and $\log$ g. It also calculates the microturbulence velocity ($\xi$) by ensuring consistency between the Fe abundances derived from lines of various equivalent widths. We adopted solar abundances from \cite{Asplund2009} and the line lists assembled in the framework of the \textit{Gaia}-ESO survey \citep{Heiter2015, Heiter2020}.

Our analysis confirms the stellar parameters derived by \cite{Jorissen19} and the updated results are listed in Table \ref{table:MCMCresults}. Our results are also in agreement with those previously derived by other authors (Table~\ref{tab:atmospheric_parameters}). A luminosity of $L = 332$~L$_\odot$ was also derived by \cite{Jorissen19} by integrating the spectral energy distribution of the star and using the \textit{Gaia} Data Release (DR) 2 parallax ($5.87\pm0.19$~mas; \citealt{GaiaDR2}), following \cite{Escorza17}. The \textit{Gaia} DR2 parallax is in agreement with its re-evaluation in DR3 ($5.88\pm0.093$~mas, \citealt{GaiaEDR3}). With the latter as an input, our Markov chain Monte Carlo (MCMC) optimisation converges to $5.86\pm0.04$~mas (Table \ref{table:MCMCresults}), so the luminosity quoted above remains valid.

Once $T_{\rm eff}$ and $L$ are confirmed, the evolutionary mass of the giant can be revisited. \cite{Jorissen19} obtained it by comparing its location in the Hertzsprung-Russell diagram (HRD) with solar-metallicity tracks from STAREVOL \citep{Siess00, Siess06}. However, evolutionary tracks differ from one another depending on the prescriptions used for overshooting, rotation, mass loss, and so on. The right panel of Fig.~\ref{hrd} shows a comparison between STAREVOL, Geneva \citep{Ekstrom12} and Padova \citep{Girardi00} models of solar metallicity ($Z_{\odot} = 0.0134, 0.014$, and $0.019$, respectively) and demonstrates that a range from 4 to 5~M$_{\odot}$\ is suitable for the red giant, independently of the used models. That range for the primary mass implies a range from 1.1 to 1.5 M$_{\odot}$\ for the companion.

\begin{figure*}[t]
\centering
\includegraphics[width=0.8\textwidth]{abundances.png}
\caption{56~UMa compared to the CEMP stars from \cite{Karinkuzhi-2021}. Left: [La/Fe]-[Eu/Fe] plane, where the dashed green line corresponds to a pure solar r-process, and the continuous blue line to s-process predictions. Right: Histograms of the abundance distance $d$ (Appendix \ref{sec:d}) of 56~UMa and the reference stars. The vertical dashed line at $d$ = 0.6 marks the threshold adopted between CEMP-s and CEMP-rs stars.
}\label{Fig:abundances}
\end{figure*}

\section{Chemical abundance profile of the red giant}\label{abundances}

To investigate the origin of the chemical enrichment of the giant, we derived its chemical abundances, including the abundances for 14 heavy elements, some of which are predominantly produced by the s-process, some by the r-process, and some by both. The atomic lines we used, including hyperfine (HF) splitting when available, are listed in Table~\ref{Tab:linelist}.

The abundances of four s-process elements (Y, Zr, La, and Ce) were published by \cite{Jorissen19}, but we present a new analysis that uses the co-added spectrum with a higher S/N. The measured abundances and the non-LTE corrections, when available, are listed in Table \ref{Tab:abundances}. A non-LTE correction of 0.19 dex, corresponding to the atmospheric parameters of the \ion{Sr}{i} line, was adopted from \cite{Bergemann2012}. Table \ref{Tab:uncertainties} lists the sensitivity of the abundances to variations in the atmospheric parameters (determined following \citealt{Karinkuzhi18}). Figure \ref{Fig:HD98839_synth} shows the fit to a few heavy-element lines.

56~UMa has a (moderate) enhancement in s-process elements, consistent with its mild Ba-star classification, but also an unexpected r-process enhancement. A small group of carbon-enriched metal-poor (CEMP) stars are known to have such a hybrid (s+r) profile \citep{Karinkuzhi-2021}, but contrarily to these CEMP-rs stars, 56~UMa has a solar metallicity and a solar C/O ratio. The left panel of Fig. \ref{Fig:abundances} shows the location of 56~UMa in the [La/Fe]-[Eu/Fe] plane compared to the CEMP stars studied by \cite{Karinkuzhi-2021}. Additionally, we computed an abundance distance, $d$, from the solar r-process abundance pattern, which is described in Appendix \ref{sec:d}. The r-process was adopted as reference because, although it shows some variations \citep{Goriely97b,Roederer11}, it appears less sensitive to its operating conditions (metallicity, radiative or convective environment, mixing of protons into the C-rich intershell, etc.) than the s-process. The abundance distance, $d$, allowed us to discern whether the star abundance profile is pure r ($d < 0.25$), predominantly s ($d > 0.6$), or hybrid s+r ($0.25 < d < 0.6$) as defined by \cite{Karinkuzhi-2021}. 56~UMa, with $d=0.54$, falls among the CEMP-rs objects (right panel of Fig.~\ref{Fig:abundances}). \cite{Karinkuzhi-2021} showed that AGB nucleosynthesis may account for such a mixed s+r profile, but only at low metallicity. To our knowledge, 56~UMa is the first object at solar metallicity that shows an enrichment of both s- and r-process origin.

\section{A neutral hydrogen cavity}\label{sec:cavity}
We also investigated the morphology of the interstellar gas in the vicinity of 56~UMa. Figure \ref{cavity} shows an integrated $\lambda$ 21 cm map of neutral atomic hydrogen in Galactic coordinates, centred on 56~UMa for H\textsc{i} with velocities between $\pm$~25~km~s$^{-1}$ (the center-of-mass velocity of 56~UMa is 0.13 km/s; Table~\ref{table:MCMCresults}). The data are from the HI4PI survey \citep{Bekhti16}, which has an angular resolution of $\Theta_{\rm FWHM}=16^{\prime}2$, a sensitivity of $\sigma = 43$~mK, and full spatial sampling of 5$^{\prime}$ in both Galactic longitude and latitude.

While there is complex structure in the area, the image shows a clear absence of emission, a cavity, surrounding 56~UMa, although the star is not exactly at its centre. \cite{Schmelz2022} propose that this cavity is associated with the remnant of a SN that exploded about $10^5$~years ago. Its rim is steepest in the lower-right quadrant, where we have the clearest view of the material that may have piled up during the snow-plough phase. The minimum column density is found at Galactic coordinates $l, b = 170.67, 62.42$ is 6.5 $\times$ 10$^{18}$ cm$^{-2}$, and there is also about 1.2M$_{\odot}$\ \citep{Schmelz2022} of low-velocity gas in the cavity that could be associated with mass lost by the giant after the explosion. The clear identification of this cavity is only possible because it is within the Local Chimney \citep{Lallement03}, a low-density extension of the Local Bubble \citep{Zucker22} that stretches from the solar neighbourhood to the Galactic halo.

The unseen companion of 56~UMa has traditionally been identified as a WD, but if it were a NS, it could be the remains of the SN that evacuated the MI cavity and blasted MI itself outwards at 120~km/s \citep{Schmelz2022}. Then, one might expect to find evidence for high energy emission, as was indeed found by ROSAT as diffuse \nicefrac{1}{4}~keV X-rays in the lower-left half of the cavity, the area of the lowest H\textsc{i} emission \citep{Herbstmeier95,Schmelz2022}.

\begin{figure}[t]
\centering
\includegraphics[width=0.4\textwidth]{Cavity_arrow.pdf}
\caption{Galactic-coordinate map of H\textsc{i} at $\lambda$ 21 cm integrated between $\pm$ 25 km s$^{-1}$ and centred on 56 UMa. The arrow corresponds to the motion of the star over the last $10^5$~yr \citep{Schmelz2022}, and was derived from the orbital solution (Table~\ref{table:MCMCresults}). Colour bar units are kelvins.}\label{cavity}
\end{figure}

\section{Summary, discussion and conclusions}\label{sec:discussion}
We have presented the orbital, stellar, and chemical properties of the wide binary system 56~UMa and determined the mass of the unseen component to be 1.31$\pm$0.12~M$_{\odot}$, which is compatible with both a WD and a NS. In the remainder of this section, we discuss the possible configurations of the system along with their pros and cons.

\begin{enumerate}
\item The faint companion could be a 1.3~M$_{\odot}$\ WD, implying an AGB progenitor of a rather high mass, of the order of 6~M$_{\odot}$\ \citep{Marigo22}. However, massive AGB stars are not expected to be efficient enough s-process producers \citep{Kappeler11} to account for the observed barium overabundance. Additionally, the lack of rubidium overabundance, which is considered the major signature of super-AGB nucleosynthesis \citep{GarciaHernandez06, Doherty17}, is also an argument against this scenario. A new exploration of the nucleosynthesis parameter space could add to this discussion. Finally, this scenario would require the observed H\textsc{i} cavity to have been created by an AGB (super-)wind or by a last-thermal pulse ejection event that produced a shock wave that cleared out the accumulated material around the AGB star. However, we would then need to understand why not all AGB stars are surrounded by interstellar medium cavities. As already mentioned, this cavity was discovered only because it is in a low-density environment, and such favourable conditions probably do not occur for each AGB star. An alternative option is that the cavity is not associated with any of the components of the 56~UMa system, given its off-centre location.

\item The companion could be a close pair formed by two low-mass stars \citep[e.g.][]{vandenheuvel2020}. If they were both main-sequence stars, mass transfer could not explain the peculiar abundances of the giant. Therefore, at least one companion must be an evolved star, a WD. Assuming we can extrapolate the stability criterion discussed by \cite{Tokovinin14} to this more evolved system, the ratio between the two periods must be higher than 4.7. If 56~UMa were a triple system, the inner orbit would have to be shorter than 3600~days. Most Ba stars are in orbits below that threshold \citep{Jorissen19, Escorza19}, so that period should give enough space for an AGB star to evolve. A problem in this scenario is that this WD must still originate from a star more massive than the initial mass of the current primary, so we cannot circumvent the difficulties listed for the previous scenario.

\item The faint companion could be a NS. This scenario can explain the cavity better than the previous ones; however, the SN rate in the Milky Way is estimated to be one to two every 100~yr \citep[e.g.][]{Rozwadowska21}, so the probability of such an event occurring within 170~pc of Earth and within $10^{5}$~yr is very low. Even if we spotted such an event, the fact that the system is bound, yet very eccentric, together with the current location of the giant on the HRD require dynamical and evolutionary fine-tuning. We estimate that for this scenario to work, it would require an asymmetric explosion, with a small associated kick opposite to the orbital motion, in order to match the mass budget imposed by the evolutionary timescales (see Appendix \ref{sec:kick}). \cite{Schmelz2022} suggest that the SN exploded about $10^{5}$~yrs ago, and the Padova isochrones
\citep{Bressan2012} predict that a solar-metallicity 8~M$_{\odot}$\ star would reach carbon ignition at an age of 40 Myr. At that age, a 4.3~M$_{\odot}$\ star (the current mass of the giant) would still be on the main sequence, which is incompatible with its current status. A possible way to solve this tension is to assume that the initial masses of the two components were almost equal, and that the SN blast removed a couple of solar masses from the current red giant envelope. \cite{Leigh16} show, however, that red giant stars must be nearly completely stripped for their luminosity to change significantly. Hence, 56~UMa should have kept a luminosity typical of a 7-8 M$_{\odot}$ star, which it did not (Fig. \ref{hrd}). Additionally, numerical simulations of SN explosions in binaries \citep{Hirai2014} indicate that there is no mixing between the SN ejecta and the companion material, so the SN would leave no chemical imprint on the surviving companion. Thus, the chemical peculiarities detected on the giant remain unexplained in this scenario too.
\end{enumerate}

We conclude that none of the potential scenarios can explain all observables, and we cannot yet pin down the nature of the faint companion in 56~UMa. A crucial ingredient is whether or not the H\textsc{i} cavity is a SN remnant associated with 56~UMa. Even if it is not, the massive chemically peculiar giant in combination with the massive unseen companion is difficult to reconcile with stellar evolution.  We stress that the chemistry in 56~UMa, bearing imprints from both the s- and r-process of nucleosynthesis, is the first such star at solar metallicity. This system provides a strong test for binary evolution.\\

\begin{acknowledgements}
We thank N. Chamel, L. Siess and S. Goriely for the useful discussions. The HERMES spectrograph is supported by the Fund for Scientific Research of Flanders (FWO), the Research Council of KU Leuven, Belgium, the Fonds National de la Recherche Scientifique (F.R.S.-FNRS), the Royal Observatory of Belgium, the Observatoire de Genève, and the Thüringer Landessternwarte Tautenburg. This work makes use of SIMBAD, operated at CDS, Strasbourg, and of data from the European Space Agency (ESA), processed by the \textit{Gaia} Data Processing and Analysis Consortium (DPAC).
\end{acknowledgements}

\bibliographystyle{aa}


\section{Transition rates}\label{aptr}
As stated in the main text, collective effects from ordered structures have
been investigated for two models for $q=2$ and $q=3$.  They can be grouped in a common expression for $q=2$, in which engine B is precisely set by $\alpha=0$.
The system dynamics  is governed  by the master equation $ {\dot p}_j=\sum_{\nu=1}^2\sum_{i\neq j}(\omega^{(\nu)}_{ji}p_i-\omega^{(\nu)}_{ij}p_j) $, where
transition rates  $i$ to  $j$ are  given by
  \begin{equation}
  \omega^{(1)}_{ji}=\Gamma e^{-\frac{\beta_1}{2}\{\mp\epsilon(1+\alpha)(1-\frac{2N_{i\uparrow(\downarrow)}\mp 1}{N}) \mp F\}} \quad {\rm and} \quad
  \omega^{(2)}_{ji}=\Gamma e^{-\frac{\beta_2}{2}\{\mp\epsilon(1+\alpha)(1-\frac{2N_{i\uparrow(\downarrow)} \mp 1}{N}) \pm F\}},
  \end{equation}
  where signs $\mp\epsilon(1+\alpha) $ attempt for the transition from state $i$ to $j$ in which $N_{i\uparrow}\rightarrow N_{i\uparrow}\pm 1$ and $N_{i\downarrow}\rightarrow N_{i\downarrow}\mp 1$, respectively, with associate driving contribution $\mp F$ and $\pm F$ due to the contact  with the cold ($\nu=1$) and hot ($\nu=2$) thermal baths, respectively.
  For a given $N$, the dynamics is generated from Gillespie algorithm
  \cite{gillespie1977exact} in which each above transition is chosen with probability
   $\omega^{(\nu)}_{ji}/{\tilde w}$, where
  $\tilde w=\Gamma \left(N_{i\uparrow} e^{-\frac{\beta_1}{2}\{\epsilon(1+\alpha)(1-\frac{2N_{i\uparrow}+1}{N}) + F\}}+N_{i\uparrow} e^{-\frac{\beta_2}{2}\{\epsilon(1+\alpha)(1-\frac{2N_{i\uparrow}+1}{N}) - F\}}+N_{i\downarrow} e^{-\frac{\beta_1}{2}\{\epsilon(1+\alpha)(1-\frac{2N_{i\downarrow}+1}{N}) - F\}}+N_{i\downarrow} e^{-\frac{\beta_2}{2}\{\epsilon(1+\alpha)(1-\frac{2N_{i\downarrow}+1}{N}) + F\}}\right)$.

The thermodynamic limit $ N\rightarrow \infty$ is obtained by replacing
$N_{i\uparrow(\downarrow)}/N\rightarrow p_{\uparrow(\downarrow)}$, in which the time evolution of probability $p_\beta$ [$\beta\in(\downarrow,\uparrow)$]  is ruled by
the master equation
${\dot p}_\beta=\sum_{\nu=1}^2J_{\beta \beta'}^{(\nu)}$, where
$J_{\beta \beta'}^{(\nu)}=\omega^{(\nu)}_{\beta \beta'}p_{\beta'}-\omega^{(\nu)}_{\beta'\beta}p_\beta$ with transition rates listed below:
\begin{equation}
  \omega^{(1)}_{\uparrow\downarrow}=\Gamma e^{-\frac{\beta_1}{2}\{-\epsilon(1+\alpha)(1-2p_{\uparrow}) - F\}} \quad {\rm and} \quad
  \omega^{(2)}_{\uparrow\downarrow}=\Gamma e^{-\frac{\beta_2}{2}\{-\epsilon(1+\alpha)(1-2p_{\uparrow}) + F\}},
  \end{equation}
 for $(\beta',\beta)\equiv(\downarrow,\uparrow)$ and
 \begin{equation}
  \omega^{(1)}_{\downarrow\uparrow}=\Gamma e^{-\frac{\beta_1}{2}\{\epsilon(1+\alpha)(1-2p_{\uparrow}) + F\}}\quad {\rm and} \quad
  \omega^{(2)}_{\downarrow\uparrow}=\Gamma e^{-\frac{\beta_2}{2}\{\epsilon(1+\alpha)(1-2p_{\uparrow}) - F\}},
  \end{equation}
for $(\beta,\beta')\equiv(\downarrow,\uparrow)$   respectively.

  Transition rates are evaluated
 in a similar fashion for $q=3$, in which energy differences can be grouped
in three distinct sets.
Starting with engine A, they are identical to $q=2$
for transitions of type $\uparrow(\downarrow)\rightarrow \downarrow(\uparrow)$,
whereas
the energy difference reads $\epsilon(N_{ik}-\alpha N_{i\ell})/N$ for transitions like $0\rightarrow \uparrow(\downarrow)$, where $k=\uparrow(\downarrow)$ and $\ell= \downarrow(\uparrow)$
 (likewise for all other ways round of transitions).
Likewise, for engine B,  a given transition $N_{i\ell}\rightarrow N_{i\ell}- 1$ and
  $N_{ik}\rightarrow N_{ik}+ 1$ (where $k,\ell \in (\uparrow,0,\downarrow)$)
  has energy difference given by
  $\epsilon (N_{ik}-N_{i\ell}+1)/N$
  \cite{herpich,herpich2}. Numerical simulations are performed
  as before, but now there are $2\,q\,(q-1) = 12$ distinct transitions. As for $q=2$,
 the limit $N\rightarrow \infty$ is promptly obtained and described by the master equation $ {\dot p}_\beta=\sum_{\nu=1}^2\sum_{\beta'\neq \beta} J_{\beta \beta'}^{(\nu)}$ [$\beta\in(\downarrow,0,\uparrow)$], but now
 there are six distinct transitions. Below, we list some of them:
 \begin{equation}
  \omega^{(1)}_{\uparrow\downarrow}=\Gamma e^{-\frac{\beta_1}{2}\{-\epsilon(1+\alpha)(p_\downarrow-p_\uparrow) + F\}}, \qquad
  \omega^{(1)}_{\uparrow 0}=\Gamma e^{-\frac{\beta_1}{2}\{\epsilon(p_\uparrow-\alpha p_\downarrow)- F\}} \quad {\rm and} \quad
  \omega^{(1)}_{0\downarrow}=\Gamma e^{-\frac{\beta_1}{2}\{\epsilon(\alpha p_\uparrow- p_\downarrow)- F\}},
  \end{equation}
for engine A and
\begin{equation}
  \omega^{(1)}_{\uparrow\downarrow}=\Gamma e^{-\frac{\beta_1}{2}\{\epsilon(p_\uparrow-p_\downarrow) + F\}}, \qquad
  \omega^{(1)}_{\uparrow 0}=\Gamma e^{-\frac{\beta_1}{2}\{\epsilon(p_\uparrow-p_0)- F\}}
\quad {\rm and} \quad
  \omega^{(1)}_{0\downarrow}=\Gamma e^{-\frac{\beta_1}{2}\{\epsilon(p_0- p_\downarrow)- F\}},
  \end{equation}
 for engine B, where in all cases  the limit $N_{i\beta}/N \rightarrow p_\beta$ was taken. The  remaining  transitions  are calculated analogously and    corresponding
$\omega^{(2)}_{ij}$'s  are promptly obtained from $\omega^{(1)}_{ij}$'s just by replacing $F\rightarrow -F$. Finally,
thermodynamic quantities expressions in the limit $N\rightarrow \infty$
are similar to those presented in
the main text (for finite  $N$), where the power ${\cal P}$ and heat $\langle\dot{Q}_\nu\rangle$  per unity are given by
\begin{equation}
\label{work1}
{\cal P}=-\sum_{(\nu,\gamma)} F_\gamma^{(\nu)}\sum_{(\beta,\beta')}d_{\beta \beta'}^{(\nu)}J_{\beta \beta'}^{(\nu)} \qquad{\rm and} \qquad
   \langle \dot{Q}_\nu\rangle=\sum_{(\beta,\beta')}\left(E_\beta-E_{\beta'}+\sum_{\gamma}F_\gamma^{(\nu)}d_{\beta \beta'}^{(\nu)}\right)J_{\beta \beta'}^{(\nu)},
\end{equation}
where the energy difference $E_\beta-E_{\beta'}$ was evaluated previously for each transition.

\section{Phenomenological description for the probability distribution in the  regime of strong collective effects} \label{apc}
Despite the nonlinear shape of the master equation for $N\rightarrow \infty$, it
is possible  to get some
insight about the probability distribution
in the regime of strong collective
effects,  as described as follows. Starting
with $q=2$, the  ordered phase has degeneracy 2 and characterized by the predominance of a spin of type ($p^{st}_\downarrow>>p^{st}_\uparrow$ or $p^{st}_\uparrow>>p^{st}_\downarrow$). By focusing on the
former solution,  steady probabilities are approximately
given by
\begin{equation}
p^{st}_\downarrow \approx \frac{\omega^{(1)}_{\downarrow\uparrow}+\omega^{(2)}_{\downarrow\uparrow}}{\omega^{(1)}_{\downarrow\uparrow}+\omega^{(2)}_{\downarrow\uparrow}+\omega^{(1)}_{\uparrow\downarrow}+\omega^{(2)}_{\uparrow\downarrow}},
\label{theor}
\end{equation}
where $p^{st}_\uparrow=1-p^{st}_\downarrow$ and,
at this level of approximation, the probabilities
 $p^{st}_\uparrow$ and $p^{st}_\downarrow$ appearing in transition rates have been neglected.
 Taking into account that $\beta_2<\beta_1$, $F>0$ and
 $-\epsilon>>F$ (strong collective effects), it follows that $\omega^{(2)}_{\downarrow\uparrow}>\omega^{(1)}_{\downarrow\uparrow}\gg\omega^{(1)}_{\uparrow\downarrow}>\omega^{(2)}_{\uparrow\downarrow}$ and  Eq. (\ref{theor}) is approximately
 given by $p^{st}_\downarrow \approx 1-(\omega^{(1)}_{\uparrow\downarrow}+\omega^{(2)}_{\uparrow\downarrow})/(\omega^{(1)}_{\downarrow\uparrow}+\omega^{(2)}_{\downarrow\uparrow})$. Since $\omega^{(2)}_{\downarrow\uparrow}>\omega^{(1)}_{\downarrow\uparrow}$
 and $\omega^{(1)}_{\uparrow\downarrow}>\omega^{(2)}_{\uparrow\downarrow}$,  the ratio $\omega^{(1)}_{\uparrow\downarrow}/\omega^{(2)}_{\downarrow\uparrow}$ is dominant
  and finally $p^{st}_\downarrow$  acquires the simpler form
$p^{st}_{\downarrow}\approx 1-e^{\frac{1}{2}\{(\beta_1+\beta_2)(1+\alpha)\epsilon\}}$$e^{\frac{(\beta_1-\beta_2)F}{2}}$, where $\epsilon<0$. Note that it has the general form $p^{st}_\downarrow=1-\gamma e^{\delta x}$ described in the main text, where  $\delta=(\beta_1+\beta_2)(1+\alpha)/2$ and $\gamma=\exp(\beta_1-\beta_2)F/2$  for $x=\epsilon$ and $\delta=(\beta_1-\beta_2)/2$ and $\gamma=\exp\{(\beta_1+\beta_1)(1+\alpha)\epsilon/2\}$
for $x=F$.

By inserting above expressions for $p^{st}_\downarrow$ and $p^{st}_\uparrow$ into expressions from Eq. (\ref{work1}), one arrives at the following expressions for ${\cal P}_{ph}\equiv {\cal P}$ and $\langle \dot{Q}_2\rangle_{ph}\equiv \langle \dot{Q}_2\rangle$ described in the main text.
All similar procedures can be extended for $q=3$ just by replacing $\uparrow\rightarrow 0$ and neglecting the other transition rates. In this case,
$p^{st}_\downarrow \approx 1-\omega^{(1)}_{0\downarrow}/\omega^{(2)}_{\downarrow0}=1-e^{\frac{1}{2}\{(\beta_1+\beta_2)\epsilon+(\beta_1-\beta_2)F\}}$.
Fig. \ref{appr} exemplifies   the validity of  approximate expression for $p^{st}_{\downarrow}$ for distinct sets of temperatures
$\beta_1,\beta_2$ and $F$ for $q=3$.
\begin{figure}[h!]
          \includegraphics[scale=0.5]{Fig1sp.pdf}
               \caption{Monolog plot of $1-p^{st}_\downarrow$ versus $\epsilon$ for
               distinct sets of $\beta_1,\beta_2$ and $F$. Continuous lines and symbols
               denote exact results for $N\rightarrow \infty$ and those from expression $1-p^{st}_{\downarrow}=e^{\frac{(\beta_1+\beta_2)\epsilon+(\beta_1-\beta_2)F}{2}}$, respectively. }
       \label{appr}
    \end{figure}

\section{Thermodynamics of q=2 engines}\label{apt}
The main features of the thermodynamics of $q=2$ engines are
summarized in Figs. \ref{fig0} and   \ref{fig4}. The former shows, for finite $N$ and
$N\rightarrow \infty$, the power per unity, efficiency and reliability of the two-state
the model discussed in the main text:
\begin{figure}[h!]
    \centering
   \includegraphics[scale=0.23]{Fig2sp.pdf}
\caption{Left: Schematics of the driving operation for $q=2$ engines. For $\beta_1=2,\beta_2=F=1$ and distinct $\alpha$'s, panels $(a)$-$(b)$  depict the efficiency $\eta$ and the power per unity ${\cal P} $ versus interaction strength $\epsilon$, respectively.   Circles show the optimal $\epsilon_{ME}$'s in which $\eta=\eta_c$. In both cases, symbols correspond to the phenomenological description. Insets: plot of $p^{st}_\downarrow$ and monolog plot  $p_\uparrow$  versus $\epsilon$, respectively. For the same  parameters as in the top panels, the bottom panels show ${\cal P}=\langle{\mathtt P}\rangle/N$ and $\eta$  for finite $N$'s and $\alpha=1$. }
 \label{fig0}
    \end{figure}

Fig. \ref{fig4}  extends the power and efficiency heat maps for $\alpha=1$. Note the similar behavior to $q=3$, but in this case ideal heat engine operations are obtained for all values
of $\alpha$ at $\epsilon_{ME}$ ($F$ fixed) or $F_{ME}$ ($\epsilon$ fixed). They satisfy the conditions ${\cal P}=\langle {\dot Q_2}\rangle =0$
and according to the two-state model, ideal efficiencies (dashed lines) fulfill the condition $M=-\tanh \left[\frac{\beta_2}{2} (F+(\alpha +1) M \epsilon )\right]$
at $\epsilon_{ME}/F_{ME}$ when $F/\epsilon$ is held fixed.
    \begin{figure}[h!]
          \includegraphics[scale=0.35]{Fig3sp.png}
               \caption{For the same parameters from Fig. \ref{fig0}, the heat maps for the ${\cal P}$ and $\eta$. Heat engine, pump, and dud regimes are described
               by symbols $HE,P$, and $D$, respectively. Continuous and dotted
               lines denote maximizations of power for $F$ and $\epsilon$ held fixed,
               respectively. The dashed line corresponds to the crossover from heat engine to pump regimes, in which $\eta=\eta_c$.}
       \label{fig4}
    \end{figure}

\section{Heat maps for q=3 and N=2 engines}\label{apn2}
This section portrays two important aspects discussed in the main text: the reliability of numerical simulations for finite $N$ and the fact a minimal setup of $N=2$ interacting unities already captures the essential ingredients introduced here.  For such a minimal setup, collective effects are measured by a dominant fraction of quantity $-M=p^{st}_{\downarrow\downarrow}-p^{st}_{\uparrow\uparrow}$.

In both cases, results will be shown for $\beta_1=1$ and $\beta_2=0.4$, in which the engine performance is somewhat superior to $\beta_1=2$
and $\beta_2=1$.
The former aspect is depicted in Fig.  \ref{n2exactvsgill} by comparing thermodynamic quantities evaluated from numerical simulations (Gillespie algorithm) and those from  exact  steady probabilities via largest eigenvalue of the master equation.
\begin{figure}[h!]
    \centering
    \includegraphics[scale=0.32]{Fig4sp.pdf}
 \caption{Left and right panels show the  power per unity ${\cal P}=\langle{\mathtt P}\rangle/N$  and efficiency for $N=2$.
    Continuous lines and symbols
    show the comparison between the exact solution, via exact probability
    distribution, and numerical simulations from Gillespie algorithm, respectively.}
    \label{n2exactvsgill}
\end{figure}

 Fig. \ref{fig5n} extends the heat maps for the minimal setup of $N=2$, showing that despite the substantial reduced performance, all characteristics from collective effects are already present in the minimum
 $N=2$ setup.
\begin{figure}[h!]
\includegraphics[scale=0.35]{Fig5sp.png}
\caption{For $q=3$ and $N=2$, top panels depict power per  unity ${\cal P}=\langle{\mathtt P}\rangle/N$ and efficiency heat maps, respectively.  Heat engine  and dud regimes are described by symbols $HE$ and $D$, respectively. For a better visualization  the pump regime has been omitted. Parameters: $\alpha=1$, $\beta_1=1$ and $\beta_2=0.4$.}
        \label{fig5n}
    \end{figure}

\section{Linear regime}\label{lin}
Additional insights regarding the influence of the collective effects on the efficiency and power are obtained from a linear analysis, valid near equilibrium ($\beta_1-\beta_2,F<<1$). By resorting to the linear stochastic thermodynamics ideas \cite{callen,PhysRevLett.116.220601,proesmans2015onsager,proesmans16b,fiorek}, one introduces the following
thermodynamic forces  $f_1=\beta_1-\beta_2$ and $f_2=\beta_1F$, in such a way that
the entropy production $\langle {\dot \sigma} \rangle$ is expressed
in the bilinear form
\begin{equation}
\langle {\dot \sigma} \rangle =\left(\beta_1 - \beta_2\right) \langle \dot{Q}_2\rangle+\beta_1\mathcal{P}= J_1f_1+J_2f_2,
\end{equation}
where $J_1=\langle \dot{Q}_2\rangle$ and $J_2=\mathcal{P}/F$ denote the thermodynamic fluxes. Near equilibrium, ie. in the linear regime, these fluxes can be expressed in terms of the Onsager coefficients, $J_1=L_{11}f_1+L_{12}f_2$ and  $J_2=L_{21}f_1+L_{22}f_2$, which satisfy the conditions $L_{11},L_{22}\ge 0$
and $L_{12}=L_{21}$. From the above, the efficiency $\hat{\eta}$ promptly reads
\begin{equation}
    \hat{\eta} =-\frac{\beta_1\mathcal{P}}{\left(\beta_1 - \beta_2\right) \langle \dot{Q}_2\rangle}=-\frac{L_{21}f_2f_1+L_{22}f_2^2}{L_{11}f_1^2+L_{12}f_1f_2},
    \label{efff}
\end{equation}
from which $\eta=\hat{\eta} \eta_c$ follows immediately. As previously, the heat engine  (${\cal P}<0$) or pump (${\cal P}>0$)  regimes  impose boundaries to optimization with respect to $f_2$, whose absolute value must lie in the interval $0\le \vert f_2 \vert \le \vert f_m \vert$, where $f_m=-L_{21}f_1/L_{22}$ the so-called stopping force for which ${\cal P}=0$. As previously, the optimization can be performed to ensure maximum power ${\cal P}_{\textrm{\tiny{MP}}}$  (with efficiency ${\eta}_{\textrm{\tiny{MP}}}$)
or maximum efficiency ${\eta}_{\textrm{\tiny{ME}}}$   (with power ${\cal P}_{\textrm{\tiny{ME}}}$), by adjusting
the output force $f_2$  to optimal values $f_{2,\textrm{\tiny{MP}}}$ and $f_{2,\textrm{\tiny{ME}}}$, respectively. These optimal output forces can be expressed in terms of the Onsager coefficients as
\begin{equation}
  f_{2,\textrm{\tiny{ME}}}=\frac{L_{11} }{L_{12} }\left(-1+ \sqrt{1-\frac{L_{12}^2 }{L_{11} L_{22}}}\right)f_1,
  \label{eq:x2meta}
\end{equation}
and
\begin{equation}
  f_{2,\textrm{\tiny{MP}}}=-\frac{1}{2}\frac{L_{12}}{L_{22}}f_1,
  \label{x2mp}
\end{equation}
respectively, where the property $L_{21}=L_{12}$ has been considered. By inserting  $f_{2,\textrm{\tiny{ME}}}$ or  $f_{2,\textrm{\tiny{MP}}}$ into  expressions for $\hat{\eta}$, we obtain $\hat{\eta}_{\textrm{\tiny{ME}}}$ and the efficiency at maximum power $\hat{\eta}_{\textrm{\tiny{MP}}}$  given by,
\begin{equation}
  \hat{\eta}_{\textrm{\tiny{ME}}}=-1+\frac{2L_{11}L_{22}}{L_{12}^2}\left(1-\sqrt{1-\frac{L_{12}^2}{L_{11} L_{22}}}\right),
  \label{etame}
\end{equation}
and
\begin{equation}
  \hat{\eta}_{\textrm{\tiny{MP}}}=\frac{L_{12}^2}{4 L_{11}L_{22}-2L_{12}^2},
\label{etamp}
\end{equation}
and similarly the expressions for ${\cal P}_{\textrm{\tiny{MP}}}$ and ${\cal P}_{\textrm{\tiny{ME}}}$ are obtained.
 In fact, these quantities are not independent of each other, instead they are related as
\begin{equation}
  \hat{\eta}_{\textrm{\tiny{MP}}}=\frac{\hat{\eta}_{\textrm{\tiny{ME}}}}{1+\hat{\eta}^2_{\textrm{\tiny{ME}}}} \qquad {\rm and} \qquad \frac{{\cal P}_{\textrm{\tiny{ME}}}}{{\cal P}_{\textrm{\tiny{MP}}}}=1-\hat{\eta}^2_{\textrm{\tiny{ME}}},
  \label{opti}
\end{equation}
where the symmetry between crossed Onsager coefficients $L_{12}=L_{21}$ was taken into account.
It is convenient to introduce the \emph{coupling parameter}
$\kappa=L_{12}/\sqrt{L_{11}L_{22}}$ \cite{cleuren2015universality,TF9656101897}, in such a way that maximum efficiencies  $\hat{\eta}_{\textrm{\tiny{MP}}}$ and  $\hat{\eta}_{\textrm{\tiny{ME}}}$ are solely expressed in terms of $\kappa$  given by
\begin{equation}
  \hat{\eta}_{\textrm{\tiny{ME}}}=-1+\frac{2}{\kappa^2}\left(1-\sqrt{1-\kappa^2}\right),
  \label{etame1}
\end{equation}
and
\begin{equation}
    \hat{\eta}_{\textrm{\tiny{MP}}} = \frac{1}{2}\frac{\kappa^2}{2-\kappa^2},
    \label{etamp1}
\end{equation}
respectively. Since $\langle {\dot \sigma} \rangle \ge 0$, it follows that $\kappa$ must be constrained in the interval  $-1 \le \kappa \le 1$, implying that both $\hat{\eta}_{\textrm{\tiny{MP}}}$ and $\hat{\eta}_{\textrm{\tiny{ME}}}$ are confined to $0 \le \hat{\eta}_{\textrm{\tiny{MP}}} \le 1/2$ and $0\le\hat{\eta}_{\textrm{\tiny{ME}}}\le 1$, respectively. Notice that $\kappa = \pm 1$ corresponds to the determinant of the ($2 \times 2$) Onsager Matrix being equal to zero. And this, in turn, implies proportionality between the two thermodynamic fluxes, $J_1 \propto J_2$, for all forces $f_1$ and $f_2$.  Fig. \ref{ons} shows that all signatures about collective effects are also captured by the linear regime, describing very well the system behavior near the equilibrium regime (panels $(b)$ and $(c)$). Remarkably, the increase of efficiencies toward the ideal case as $\epsilon$ and driving force goes up described in the main text is understood  from the interplay among $L_{ij}$'s in which the coupling $|\kappa|$ monotonously increases toward  the ideal limit $\kappa \rightarrow -1$ (panel $(a)$ and inset). Also, ${\hat \eta}_{ME}$ and ${\hat \eta}_{MP}$ follow Eqs. (\ref{etame1})-(\ref{etamp1}) (d).
\begin{figure}[h!]
\centering
\includegraphics[scale=0.3]{Fig6sp.png}
\caption{For engine A, $q=3$ and $\alpha=1$, the depiction of system thermodynamics near the equilibrium regime. Panel $(a)$ and inset show the  Onsager coefficients $L_{ij}$'s and coupling $\kappa$ versus interaction parameter $\epsilon$. In $(b)$ and $(c)$ shows $\hat{\eta}$  and ${\cal P}$ versus $f_2$  for $f_1=9.10^{-3}$. Continuous lines and symbols correspond to exact results and those obtained from Eq. (\ref{efff}), respectively. Each heat engine is delimited by vertical lines yielding at $f_m$. Panel $(d)$ shows the behavior of maximum efficiency $\hat{\eta}_{\textrm{\tiny{ME}}}$ and efficiency at maximum power $\hat{\eta}_{\textrm{\tiny{MP}}}$ versus $\kappa$, where continuous lines follow Eqs. (\ref{etame1}) and (\ref{etamp1}), respectively. The symbol $\bullet$ denotes the crossover
to the independent regime yielding at $\epsilon_c\approx-1.5$. }
\label{ons}
\end{figure}

 \section{Linear stability of disordered phase solution for models A and B  for q=3}\label{apb}
For completeness, we provide additional information about the crossover between collective and independent regimes for engines A  (for $\alpha=1$) and B for $q=3$ which manifests through continuous and discontinuous phase transitions, respectively.
They can be analyzed
in a similar way to their equilibrium counterparts, by means of
order parameters $M$ and $\phi={3}(p_{max}-1)/2$ ($p_{max}={\rm max}\{p^{st}_\downarrow,p^{st}_0,p^{st}_\uparrow\}$), respectively, the former one
characterized by the classical exponent $\beta=1/2$ \cite{fiorejcp2013,challa}.
However, contrasting to the equilibrium Potts model,   nonequilibrium
ingredients shift the phase transition for engine B, from continuous
to discontinuous, as exemplified in panel $(b)$ of Fig. \ref{app1}.
Alternatively, both phase transitions can be
investigated under a linear expansion of master equation given by ${\dot p}_m=\sum_nA_{mn}p_n$,  where $A$ is the Jacobian matrix with elements $A_{mn}=\partial (\omega^{(1)}_{mn}+\omega^{(2)}_{mn})/\partial p_n|_{p_n=p^*}$ evaluated at  fixed points $\sum_nA_{mn}p_n^*=0$ for $n=0,\uparrow$
 (where $p_\downarrow^*=1-p_\uparrow^*-p_0^*$). In particular, the solution  $p_{n}^*$ is linearly stable if  real parts
 of eigenvalues of Jacobian matrix are negative \cite{strogatz2018nonlinear}.
 In both cases, the independent regime is  characterized by equal population $p_{\downarrow}^*=p_0^*=p_\uparrow^*=1/3$ for $\epsilon \ge\epsilon_c$, in which associate eigenvalues can be written as $\lambda_\pm=\lambda_0\pm \lambda_1$,
 with $\lambda_0$, in both cases, given by
 \begin{equation}
 \lambda_0=-(3+\beta_1 \epsilon) \cosh \left(\frac{\beta_1 \
F}{2}\right)-(3+\beta_2 \epsilon) \cosh \left(\frac{\beta_2 F}{2}\right),
 \end{equation}
whereas $\lambda_1$'s read
\begin{align}
 \lambda_1=\Big[&6+\epsilon ^2 \left(\beta_1^2+\beta_2^2\right)+\left(\beta_1^2 \epsilon ^2-3\right) \cosh
   (\beta_1 F)-3 \cosh (\beta_2 F)+\beta_2 \epsilon ^2 \left(4 \beta_1 \cosh \left(\frac{\beta_1 F}{2}\right)
   \cosh \left(\frac{\beta_2 F}{2}\right)+\beta_2 \cosh (\beta_2 F)\right)+\nonumber \\&+12 \sinh
   \left(\frac{\beta_1 F}{2}\right) \sinh \left(\frac{\beta_2 F}{2}\right)\Big]^{1/2},
   \end{align}
 for engine A and
  \begin{equation}
     \lambda_1=i \sqrt{3}\left[\sinh \left(\frac{\beta_1 F}{2}\right)-\sinh \left(\frac{\beta_2 F}{2}\right)\right],
 \end{equation}
  for engine B,  respectively. Since $\lambda_1$ is imaginary for engine B, the linear stability of disordered solution is  granted provided $\lambda_0<0$. Conversely, for engine A,
due to the fact that $\beta_1$ and $\beta_2$ are always positive, $\lambda_-$
is always negative. Conversely, $\lambda_+$ is always positive (negative) for sufficiently $\epsilon<0$ ($\epsilon\ge 0$), respectively,  in which the order-disorder phase transition  corresponds to  a transcritical bifurcation when $\lambda_+=0$.

Fig. \ref{app1} depicts the phase diagrams
$\Delta \beta=\beta_1-\beta_2$ versus $\epsilon_c$ for different $F$
    obtained from the linear analysis. In particular, for $F=0$, both $\lambda_+$'s and
 $\lambda_-$'s read $-6$ and $-2[3+\epsilon(\beta_1+\beta_2)]$  (engine A) and $-[6+\epsilon(\beta_1+\beta_2)]$ (engine B), respectively, consistent to
phase transitions yielding at $\epsilon_c=-3/(\beta_1+\beta_2)$ and $\epsilon_c=-6/(\beta_1+\beta_2)$. The crossover from collective to independent regime smoothly changes with the driving and it is more sensitive to the difference of temperatures. Note the excellent agreement between $\epsilon_c$'s
obtained from linear
analysis and those from order parameter behaviors (bottom panels for $F=1$).
     \begin{figure}[h!]
    \centering
              \includegraphics[scale=0.3]{Fig7sp.pdf}
  \caption{For $q=3$ and  distinct $F$'s, left and right panels depict the phase diagrams $\Delta \beta=\beta_1-\beta_2$ versus $\epsilon_c$ for  engines A and B, respectively,
               obtained from the linear stability of the disordered phase. For the sake of comparison, the bottom panels
               show (for $F=1$) the  location of phase transitions from order-parameter behaviors. From left to right $\beta_1-\beta_2=0,1/3,1,3$ and $9$. In all cases, we set $\beta_2=1$.}
        \label{app1}
    \end{figure}
\twocolumngrid


\section{Introduction}
Virtual characters are present everywhere in computer assisted media; populating virtual worlds, embodying players, telling stories. Producing qualitative motion for them is crucial to achieve engaging animation and raise interest in the viewer. Character animation software provide tools and interfaces which help artists do so. Many of such tools have been developed over the years, striving to find balance between precision, ease of use, and controllability. \cite{lasseterPrinciplesTraditionalAnimation1987, catmullProblemsComputerassistedAnimation1978}

The human eye perceives a lot of information through motion: it is a complex behavior through which character convey intent, emotions or personality. A same gesture can be performed in an wide variety of ways, even by a same character. Editing an existing animation to change the way it is performed, but not the gesture itself, is then an interesting research topic. Many techniques have been proposed to tackle the problem of animation style edition, but most of them act from the highest level, and lack the fine-grained control users look for.

This paper presents a new framework for character motion edition using pose metrics as a middle-ground between raw bone manipulation and high-level "style" labels.
We introduce these intermediate-level descriptors, then propose a neural-networks-based edition pipeline which allows animators to edit a character's motion by manipulating the metrics.

\section{Related work}

Virtual characters are traditionally represented by a mesh on which various textures and materials are applied. Its motion is driven by an underlying skeleton on which the mesh is fitted ("skinned"). The role of an animator is to design the character's motion by editing the skeleton's poses at specific time stamps, called key-frames. Instead of key-framing a new animation from the ground up this way, it is also possible to capture the motion of real-life actors through dedicated motion capture (MoCap) setups. MoCap data is a good source of realistic motion, but the resulting  animation still often needs to be cleaned up and edited to fit the animator's vision.

Previous research on character animation investigates tools, algorithms and interfaces to help make authoring and editing motion a simpler process.
Increased availability and quality of natural motion acquired by MoCap systems has led to raising interest in data-oriented methods. Meaningful information can be extracted from ground truth data and put to use to provide animators with more efficient tooling.

\subsection{Animation edition}
Whether originally key-framed by an an animator or recorded in a MoCap setting, animation data is set to change for any number of reasons, ranging from tweaking details to fit an artist's vision, to re-purposing it to be used on a different character. This paper builds upon a large body of work whose common goal is to provide animators with tools facilitating the edition process. The common animation edition pipeline can be split in two phases: posing and interpolation control. During the first, animators manipulate the skeleton to define key poses at specified key frames. In the second, they manipulate the interpolation process which is used to fill the gaps between key-frames.

\paragraph*{Posing.}
In most posing techniques the skeleton is conceptualized as a set of kinematic chains \cite{magnenat-thalmannJointDependentLocalDeformations1988}. With forward kinematics (FK), the modification of a joint's position or rotation trickles down, modifying joints down the chain. The opposite technique, inverse kinematics (IK) \cite{wangCCDInverseKinematics, aristidouFABRIKFastIterative2011}, find a suitable parameterization of the chain given a target end-effector, allowing the user to pose the character much like a puppet. These two methods are at the core of most modern animation software and lead to precise results, but can be tedious to use from an artist's perspective. Other techniques and interfaces have then been studied to alleviate this problem.

Some describe alternative frameworks in which the character's pose is derived from a line drawn by the user, much like a paintbrush stroke \cite{guayLineActionIntuitive2013, choiSketchiMoSketchbasedMotion2016}. This method has also been extended to virtual reality, tracking the user's hand in 3 dimensions and interpreting their motion as strokes \cite{garciaSpatialMotionDoodles2019}. On the data-oriented side, Wu \etal \cite{wuNaturalCharacterPosing2011} use learned clustering techniques to extract from a motion database a pose which most closely match the user's input. Victor \etal \cite{victorLearningbasedPoseEdition2021} train deep neural networks to solve the IK problem, learning skeleton constraints from existing data to lift some of the burden of respecting them from the user. Finally, some work has been dedicated to designing physical devices to facilitate the posing process \cite{yoshizakiActuatedPhysicalPuppet2011, jacobsonTangibleModularInput2014}.

\paragraph*{Interpolations.} Once key poses have been set, different methods exist to control how the software will interpolate between them to produce motion. Animators usually provide the interpolation factor used for each joint's position and rotation values over the elapsed time between key frames by manipulating animation curves through control points \cite{catmullClassLocalInterpolating1974}.
In the same veins, optimization-based methods have been proposed to improve the curve edition process \cite{koyamaOptiMoOptimizationGuidedMotion2018, cicconeTangentspaceOptimizationInteractive2019}. In the family of data-oriented solutions, Chai \etal \cite{chaiConstraintbasedMotionOptimization2007} propose to learn a low-dimensional space of motion, and solve the same kind of optimization problems within it. Recently, Harvey \etal \cite{harveyRobustMotionInbetweening2020} trained deep neural networks to learn to produce inbetweening motion from key frames.

Creating high-quality, appealing motion is a complex task. Animating a character according to fit one's vision, while depicting the desired action and conveying the right emotion requires an advanced knowledge of anatomy and physics rules. All the previously mentioned work allow animators to edit animation, but they all focus on controlling the character's body at bone-level. The process, while facilitated, is still a complex one. In order to provide animators and the greater public access to fast an expressive animation tools, higher-level concepts and tools are needed.

\subsection{Motion style manipulation and transfer}

Stemming from the need of methods to reuse existing animation, another body of work has taken interest on editing animation by manipulating the motion style. The notion of \emph{style} is most commonly used as opposed to the notion of \emph{content}. Given an animation sequence, the terms are used as a way to separate \emph{what} the character is doing (the content) from \emph{how} it is doing it (the style). For example, several animations sequence could represent a character walking in various ways: fast, slow, making shorter or largers strides, swinging its arms... In this case they all share the same \emph{content}, walking, but the various \emph{styles} set them apart.

Being able to generate stylized animation, and to modify the style of an existing one, is then a desirable goal. Such control would allow the user to easily experiment with different styles, and it would be possible to generate a wide array of variations from an existing sequence.

Early solutions to the style edition problem use hand-made labels of style and content to categorize existing motion, and interpolate between categories to produce new motion. Rose \etal \cite{roseVerbsAdverbsMultidimensional1998} use such a method, coupled with an interface using a notion of verbs and adverbs to parameterize the generated motion. Brand \etal \cite{brandStyleMachines2000} describe a learned state machine, in which motion sequences are organized by shared style characteristics. More recently, Rodriguez \etal \cite{rodriguezParameterizedAnimatedActivities2019} stray away from the usual style paradigm and propose to edit animation trough high level trajectories parameterization.

Some other methods leverage the work produced in related fields, such as dance theory, to categorize motion style. Durupinar \etal \cite{durupinarPERFORMPerceptualApproach2016} build upon Laban Movement Analysis \cite{labanMasteryMovement1971} to attempt to formulate a link between personality and body motion. Aristidou \etal \cite{aristidouEmotionControlUnstructured2017} map the Laban features to a system of emotions coordinates, which can be interpolated in to edit and existing animation style.

Drawing inspiration from the success achieved by style transfer methods on images, Holden \etal \cite{holdenDeepLearningFramework2016} train a deep convolutional neural network to extract "style" features from of an animation and apply them to another. In the same vein, other types of neural networks have been used, such as variational autoencoders \cite{duStylisticLocomotionModeling2019}, adversarial learning \cite{wangAdversarialLearningModeling2020} or generative flow \cite{wenAutoregressiveStylizedMotion2021}.

On the opposite end of the spectrum from the posing methods discussed above, these style edition approaches share the common drawback of being \emph{too} high-level. Because the notion of animation style is still loosely defined (and this shortcoming is discussed further in the next section), tools using styles as control parameters could be met with incomprehension on the user's end when reaching broader audiences.

\section{Problem formulation}

The previously defined usage of \emph{style} as a separation from an animation's \emph{content} is commonly accepted in the computer animation literature. The limits of what the term itself entails are however not clearly defined.

Some associate an animation's style to emotions experienced by the character, such as "joy", "anger" or "sadness" \cite{fouratiEmilyaEmotionalBody2014}. Some link it to existing motion analysis systems, such as Laban's \cite{labanMasteryMovement1971, durupinarPERFORMPerceptualApproach2016}, or to psychological concepts such as the arousal/valence scale \cite{ekmanHeadBodyCues1967, randhavaneIdentifyingEmotionsWalking2019}. Some others use arbitrary categories, such as "gorilla", "childish", "zombie", or "neutral" \cite{xiaRealtimeStyleTransfer2015}.

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.8\columnwidth]{figures/style_categories.pdf}
    \caption{Examples of walk animations in exaggerated styles from existing motion databases \cite{holdenLearningMotionManifolds2015, xiaRealtimeStyleTransfer2015}. (1) "zombie"; (2) "old man"; (3) "childish". Distinguishing each style with no context is not straight-forward.}
    \label{fig:exaggerated_styles}
\end{figure}

We argue that this categorization is too subjective to be considered as an expressive parameter to provide to users. Various factors can introduce bias in style categories, leading edition tools based on them to fail to meet the user's expectation.

The bias lies in how the data used to calibrate the style categories is recorded. It is usually not captured in the wild, but rather played by actors \cite{xiaRealtimeStyleTransfer2015}. The source sequence are thus often overplayed, exaggerated representations of the actor's and director's vision of the categories (see Fig. \ref{fig:exaggerated_styles}). In the most grounded acquisition procedures, emotions states are elicited in a subject by exposing him to a pre-determined scenario. The bias is however not totally avoided, possibly induced by the design of the experiment \cite{mccambridgeSystematicReviewHawthorne2014, henrichWeirdestPeopleWorld2010}, or even by the divergent interpretations of a common emotion in different cultures \cite{elfenbeinUniversalityCulturalSpecificity2002}.

The main goal of animation edition research is to provide tools to animators so that they can transform their mental representation of a character motion into tangible animation data. If low-level methods are fine-grained but tedious to manipulate, and if high-level ones are likely biased and take too much of the control out of the user's hands, finding a middle ground is an important open problem. Efficient tooling, allowing easy manipulation for the naïve user and controllable enough for the professional one, are yet to be developed.

The parameterized edition paradigm, as seen in \cite{rodriguezParameterizedAnimatedActivities2019}, offer interesting perspective in working around the bias problem while still allowing higher-level control. Following this line of thought, this paper proposes an alternative solution, introducing a method to manipulate style from higher-level, but still objective, controls: pose metrics.

\section{Method}

We propose a new character animation edition paradigm in which the animator manipulates intermediate-level parameters called pose metrics, leveraging neural networks to learn and automatically account for complex morphological constraints. We first introduce this new concept, then describe the neural networks pipeline we use in our experiments. We finally describe the method we use to edit whole animation sequences according to a user's input.

In the subsequent sections the term \emph{skeleton} is used to define the virtual armature of a character composed of a hierarchy of joints and bones \cite{magnenat-thalmannJointDependentLocalDeformations1988}. We call \emph{pose} the data vector describing the concatenated position and/or orientation of each joint in the skeleton at a single time frame.  An \emph{animation sequence} is then a sequence of key poses, which represent the character's motion.

In equations the notation is as follow: we denote $A$ an animation sequence containing $N$ frames, with $p = A_n$ the pose at frame $n$. In a skeleton composed of $J$ joints, the pose can be indexed to access a specific joint's data: $p_{head}$.

\subsection{Pose metrics}
 Being able to edit an animation sequence by manipulating higher-level parameters is a desirable capacity, but the literature is polarized between low level controls (operating on crude pose data such as bone position or orientation using forward and inverse kinematics) and high level ones, using \emph{style} labels relating to emotions or behavior (i.e. "angry", "childish", etc. )  \cite{wangAdversarialLearningModeling2020,duStylisticLocomotionModeling2019,abermanUnpairedMotionStyle2020}.
These style controllers lack control in the sense that they produce a stylized animation at once based on a generalized style definition, which users might not share.

As a compromise, filling the gap between these poles, we introduce \emph{pose metrics}: higher-level, objective, style agnostic values, that can be directly computed from a single pose and do not rely on any assumed universal human perceptive behavior.

The vision for these metrics is to be as flexible as possible, so as not to restrict the user. All characters are different, and the concept must make as little assumption as possible on the way they move. That way the user is responsible for designing his own metrics, mentally link them to his own definition of a particular style, and edit the animation through them.

For the sake of an example, imagine an artist editing an animation sequence of a character walking, expressing cold rage. In his mind, the character should be taking short, fast steps, his arms tense, stretched toward the floor. With previous style edition methods, he might be able to hit an "angry" button, maybe use a slider to determine how angry the result should be. Unfortunately, it might be that the system was calibrated with another definition of angry in mind: the character now walks faster, the arms bent and swinging, with the fists close to the torso. Using metrics would allow the user to circumvent this problem. Creating a few metrics such as the angle formed by the elbows, distance between feet, or between hands and the floor, and editing the animation through them, we avoid relying on assumed style characteristics.

The only requirement for a new metric is to be computable from a single pose, and to provide a function to do so. Given a function calculating such a user-provided metric, the proposed associated framework is extendable enough to accommodate it and edit animation sequences accordingly.

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.8\columnwidth]{figures/example_metrics.pdf}
    \caption{Visual depiction of the three metrics used in this paper's examples.}
    \label{fig:metrics_example}
\end{figure}

\paragraph*{Experimental pose metrics.} As previously explained pose metrics are intended to be provided by the user. In order to demonstrate the application of our system, the following experiments use three sample metrics applicable to a human skeleton: \metric{spine\_flexion}, \metric{shoulders\_openness} and \metric{legs\_spread}.

All three represent the angle formed at the intersection of two vectors $\vec{u}$ and $\vec{v}$ obtained from the pose. Formally the angle at the intersection is computed following Eq. \ref{eq:vector_angle}.

\begin{equation}
    \label{eq:vector_angle}
    \theta(\vec{u}, \Vec{v}) = \arccos \frac{
    \Vec{u} \cdot \Vec{v}
    }
    {
        \lVert \Vec{u} \rVert \cdot \lVert \Vec{v} \rVert
    }
\end{equation}
The \metric{spine\_flexion} metric represents the angle between $\vec{y}$ the world "up" vector and the spine axis. It is computed as shown in eq. \ref{eq:spine_flexion}. It could be used to bend the character forward when in a hurry, or backward when relaxed.

\begin{equation}
    \label{eq:spine_flexion}
    SF(p) = \theta(p_{neck} - p_{pelvis}, \vec{y})
\end{equation}

Next the metric \metric{shoulders\_openness} (Eq. \ref{eq:shoulders_openness}) describes how open the character's torso/shoulders area is. Example usages might be to raise the character shoulders to denote anxiety or shyness, and to open up to denote a more relaxed or proud behavior.

\begin{equation}
    \label{eq:shoulders_openness}
    SO(p) = \theta(
        p_{spine1} - p_{rshoulder},
        p_{lshoulder} - p_{spine1}
    )
\end{equation}

Finally, \metric{legs\_spread} (Eq.\ref{eq:legs_spread}) denotes how close the character's knees are to each other. Possible usages include restricting walking strides, or tweaking whether the character appears in balance.

\begin{equation}
    \label{eq:legs_spread}
    LS(p) = \theta(
        p_{pelvis} - p_{rknee},
        p_{lknee} - p_{pelvis}
    )
\end{equation}

In this paper's equations, the value $m$ of a metric $M$ for a pose $p$ is denoted as the function $m = M(p)$.

\subsection{Learned metric pipeline}

This section introduces a pipeline through which an animation sequence can be edited using pose metrics. It is designed with modularity in mind, as we want users to be able to decide on which metric should be used at any time, as well as extending it with custom metrics if need be. We first describe the global layout of the pipeline, then each of its pieces separately.
We also explain how the pipeline can be extended by the user, and set up to satisfy any combination of metrics. Finally, we propose a procedure to edit whole animation clips, by taking advantage of our architecture to propagate changes made to a selected frame to the rest of the sequence.

\paragraph*{Pipeline layout. }
Our pipeline (Fig. \ref{fig:pose_pipeline}) relies on several small neural networks. We opt for these networks because of their low memory footprint and fast run time. Using smaller independent networks also allows us to easily configure the pipeline by adding, removing or modifying modules.

The outer part of the pipeline is an encoder-decoder model pair whose purpose is to learn a latent pose space, in which a pose representation is more compact and easier to manipulate. As the decoder can only map latent codes to the real pose space, there is no risk of accidentally generating a broken pose when averaging and interpolating samples. The edition itself is done by the metrics network, which operates only on the latent space. An instance of the metric network is trained to modify a latent pose given a \emph{single} metric value. The reduced size of the network and the fact that we re-use the already trained encoder-decoder, mean that training times are fairly short. This in turn further facilitates the possibility for users to design their own metrics, train a network, and incorporate it in the pipeline.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.9\columnwidth]{figures/pipeline_pose_maths.pdf}
    \caption{Overview of the proposed pipeline with its pose metric module. The process modifies an input pose $p$ according to an arbitrary number of user-provided target metric values $\tilde{m}_{0...i}$. An encoder-decoder pair ($E$ and $D$) maps poses to and from a latent pose space, on which the other modules operate. Each metric is processed by a dedicated network $M_{0..i}$. Their resulting latent samples $\hat{z}_{0..i}$ are averaged, then mapped back to the pose space.}
    \label{fig:pose_pipeline}
\end{figure}

\paragraph*{Encoder-Decoder.} The encoder $E$ and decoder $D$ networks are trained together to map a pose $p$ to a latent representation $z$, as shown in Eq. \ref{eq:autoencoder}. In practice they are trained to minimize the squared reconstruction error between $p$ and its reconstructed version $\hat{p}$. Both networks are simple 1-layer perceptrons, using 512 units in the hidden layer and rectified linear unit (ReLU) activation functions. In our experiment, the latent size dimension is set to 64.

\begin{equation}
    \label{eq:autoencoder}
    \begin{split}
    z &= E(p) \\
    \hat{p} &= D(z) \\
    \mathcal{L_{\textrm{ED}}} &= \frac{1}{J}\sum_{j=0}^{J}{\big(p_j - \hat{p}_j \big)^2}
    \end{split}
\end{equation}

\paragraph*{Metric network(s)}. A metric network $N_m$ learns to produce a modified latent pose $\hat{z}_m$, taking as input a starting one $z$ concatenated with a target metric value $\targetnotation{m}$ (see Eq. \ref{eq:metric_network_output}).

\begin{equation}
    \label{eq:metric_network_output}
    \hat{z}_m = N_m \big(z \oplus \targetnotation{m}\big)
\end{equation}

To reflect this objective during training, we sample two poses from the same original animation sequence: an input $p = A_t$ and a target $\targetnotation{p} = A_{t+n}$. Sampling from a neighboring frame prevents the target pose from being too far from the starting one and ensures the network learns to \emph{edit} a pose rather than to generate an unrelated one. In our experiments, $n$ is sampled randomly for each iteration in $n \in [-10, 10]$. As we target an interactive edition context, aiming for smaller changes in metric is reasonable.
The input target metric is computed from $\targetnotation{p}$ and fed to the network along with $z = E(p)$. The resulting modified latent pose $\hat{z}$ is mapped back to the pose space using the decoder, i.e. $\hat{p} = D(\hat{z})$. Finally the networks are trained to minimize the mean squared error between the target pose and the modified one, as shown in the loss function in Eq. \ref{eq:metric_loss_fn}.

\begin{equation}
    \label{eq:metric_loss_fn}
    \mathcal{L} = \frac{1}{J} \sum_{j=0}^{J}{\big( \targetnotation{p_j} - \hat{p}_j \big)^2}
\end{equation}

Keeping with the idea of smaller networks we opt for a single layer perceptron with 126 hidden units and ReLU activation. The input layer's number of units is equal to the size of the auto-encoder's latent dimension plus one (to account for the metric value). The output is the same size as the latent dimension. We empirically find that this architecture suits our needs, larger/deeper networks yielding no perceptible value. If a users deemed a different one more suitable for one of his metrics, the only requirements would be on the size of the input and output layers.

\paragraph*{Data}. We use a pose dataset of about 1.5 million poses as described by Victor \etal \cite{victorLearningbasedPoseEdition2021}. The dataset is obtained by merging several pre-existing animation databases\cite{fouratiEmilyaEmotionalBody2014, carnegiemellonuniversityCMUGraphicsLab2003, holdenDeepLearningFramework2016}.
Each animation clip is retargeted to a common skeleton \cite{holdenLearningMotionManifolds2015}. The poses data is composed of the positions of each joint relative to the projection of the pelvis joint on the floor. When used to train the networks, the dataset is normalized by subtracting the mean and dividing by the standard deviation of each component value.

\paragraph*{Training parameters}. The encoder-decoder and metrics networks are trained separately, but using the same hyper-parameters: using the Adam optimizer \cite{kingmaAdamMethodStochastic2015} with a learning rate of $0.0001$, and batches of 1024 pose pairs until convergence, which takes around 15 minutes on an NVIDIA Quadro T1000.

\paragraph*{Pipeline modularity}
In cases where it is desirable to edit multiple metrics at once, it is possible to combine the results of a set of metric networks. As the output of the network is a latent representation of a pose we can combine them. In \cite{victorLearningbasedPoseEdition2021}, the authors ran multiple networks sequentially. We instead suggest averaging multiple latent poses, so that we don't have to pick an order for the modules.

As explicited in Eq. \ref{eq:metrics_average}, the source pose $p$ is first mapped to the latent space, then each specialized metric network $M{0..i}$ produces modified latent result. These intermediary results are all averaged, and finally mapped back to pose space by the decoder.

\begin{equation}
    \label{eq:metrics_average}
    \hat{p} = D\Big(\frac{1}{n}\sum_{i=0}^{n}{M_i\big(E(p)\big)}\Big)
\end{equation}

This method operating on single key-frames, can be used conjointly with other posing solutions, and integrated with interpolation methods, to offer as much control to the user as possible. Using a shared latent space also allow us to incorporate other modules from the literature into the pipeline. For example, the IK solver from \cite{victorLearningbasedPoseEdition2021} operates on a similar latent space, and, if trained on the pipeline's one, its output can be averaged with the metric network ones.

\subsection{Animation edition}

The pose pipeline described above only operates on a single pose at once. In a real-world scenario however animators typically edit a specific key frame, expecting their changes to be gradually applied to the neighboring ones. To meet these expectations we propose a method to edit an entire animation sequence through a single set of target metrics.

Our approach (Fig. \ref{fig:pipeline_anim}) resembles the traditional animation curve paradigm, using a user provided weight curve $W$ to control how much each key pose will be impacted by a modification. We once again leverage the good interpolation properties of the latent space.

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.9\columnwidth]{figures/pipeline_anim.pdf}
    \caption{Visual overview of the proposed method to echo modifications made to a key pose $A_n$ of an animation clip. The blend operation occurs between latent poses (source $z_n$ and modified by the metric module $\overline{z}_n$), and use a factor $W(n)$ sampled from the user-provided curve.}
    \label{fig:pipeline_anim}
\end{figure}

To edit an animation $A$, the user selects a key frame $A_n$ and provides an arbitrary number of target metrics values $\tilde{m}_{0..i}$ as well as a weight curve $W$. For each frame $A_t$ in the sequence, a latent representation $z_t$ is encoded. Each of them is then processed by the pose metric module, resulting in an intermediary latent pose $\overline{z}_t$. The final modified latent pose $\hat{z}_t$ is given by blending the between $z_t$ and $\overline{z}_t$ by the $W(t)$ factor (Eq. \ref{eq:interpolation}).

\begin{equation}
    \label{eq:interpolation}
    \hat{z}_t = \big(1 - W(t)\big) \cdot z_t + W(t) \cdot \overline{z}_t
\end{equation}

In our experiments $W$ is a hat function peaking at $1$ at the selected frame, and reaching $0$ at the clips extremities. This weight function is however an external parameter: experiments with a sinusoidal yield comparable results, and the curve could be edited by animators to fit their expectations, much like existing animation curves. We also note that each one of the clip's poses is processed independently from the others. The edition can thus happen in parallel, maintaining real time interactivity.

\section{Results}

This section presents visual animation edition results of the pipeline in different scenarios. More examples can be found in the accompanying video.

\subsection{Pose edition}

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.9\columnwidth]{figures/metric_edition.pdf}
    \caption{Incremental modification of a pose by providing a target $spine\_flexion$ input. On each step, the user's input is a $0.1$ increase over the computed metric value.}
    \label{fig:spine_flexion_pose}
\end{figure}

In Fig. \ref{fig:spine_flexion_pose}, a single pose is modified by a user by providing gradually increasing target \metric{spine\_flexion} metric values. The results respect the desired changes as shown by the measured metric values below each pose. While it could be argued that changing the orientation angle of the pelvis bone would produce similar results, a closer look to the resulting poses illustrates the interesting properties of the latent pose space. At every iteration, the skeleton's arms are slightly raised, and its knees bent, so as to compensate the loss of balance induced by the flexion of the spine.

\subsection{Modular pipeline}

\begin{figure}[h]
    \centering
    \includegraphics[width=0.6\columnwidth]{figures/multiple_metrics_and_ik.pdf}
    \caption{Sample results of the pipeline using multiple modules. In (1), \metric{shoulders\_openness} (so) and \metric{spine\_flexion} are used simultaneously. In (2) we incorporate the latent inverse kinematic solver from \cite{victorLearningbasedPoseEdition2021} into the pipeline, using it along \metric{spine\_flexion}.}
    \label{fig:multiple_metrics_results}
\end{figure}

Fig. \ref{fig:multiple_metrics_results} presents some results obtained by averaging the latent outputs from multiple modules. The first example show the result of modifying the \metric{spine\_flexion} and \metric{shoulders\_openness} metrics at the same time. Both modifications can be seen on the skeleton: its spine is straightened, and its shoulders lowered. The rest of the pose is also slightly edited to accommodate for those change: the hands are lowered, and knees bent a little further.

The second one shows the result of using the \metric{spine\_flexion} metric module along with the neural IK method from \cite{victorLearningbasedPoseEdition2021}. Averaging the latent outputs leads our method to find a compromise between each provided constraint. Here, the skeleton's back is straightened, and the IK targets for the hands are reached for.

\subsection{Animation edition}

\begin{figure}[b]
    \centering
    \includegraphics[width=0.9\columnwidth]{figures/anim_edition_results.pdf}
    \caption{Editing an animation clip according to the user input. The output is produced by the pipeline given a selected frame $A_n$, a target value for the \metric{legs\_spread} metric, and a weight curve. }
    \label{fig:animation_edition_result}
\end{figure}

Fig. \ref{fig:animation_edition_result} showcases the result of an animation $A$ being edited using the pipeline described earlier. The user selects a frame $A_n$, provides a target \metric{leg\_spread} metric value, and and interpolation factor curve, a hat function peaking at one on the $A_n$ and reaching down to 0 on both sides at $A_{n-3}$ and $A_{n+3}$. The remaining frames use a factor of 0.

Each frame pose is processed separately by the network, and the result is given by interpolating the latent representations of the input and output poses, using the frame-specific factor.

As a result, the output animation $\hat{A}$ is modified, with the impact being most important on the selected frame and gradually lowering on each side. At frame $n$ the output pose is modified the most, and the reduced legs angle is very visible. At frames $n-2$ and $n+2$ this impact is barely visible, while further frames are not modified at all, as expected.

Interpolating between original and modified poses allow the whole sequence to keep its temporal coherence even though the pipeline processes one frame at a time, with no knowledge of the past and future frame context. We also take advantage of the properties of latent spaces: any latent interpolation result is mapped back to the real pose space, thus avoiding artifacts such as skeleton interpenetration.

\section{Conclusion}
In order to avoid the pitfalls of categorizing animation styles, we introduce the concept of pose metrics, values computed on a single pose that user can define themselves to fit their \emph{own} perception of styles. We propose a pipeline to edit an animation by manipulating these metrics, which uses a modular set of tiny neural networks, running in real time. A user can extend the pipeline by designing a custom metric and training an associated metric network in a short 15 minutes.

The main shortcoming of this approach is that it is limited to a single standardized humanoid skeleton, with which the networks are trained. This is a common problem with data-driven methods: in order to adapt to a new skeleton, the network would need to be retrained. This could be difficult if access to existing animation data for the new skeleton is limited, but recent progress in animation retargeting \cite{abermanSkeletonawareNetworksDeep2020} could help by allowing more accessible data augmentation. The modularity of metrics helps here as well, as different skeleton will require a new set of metrics.

Future work will focus on grounding the proposed pipeline: securing the process by automatically detecting overlapping metrics, and conducting user studies to evaluate the usability of the method.

\bibliographystyle{ACM-Reference-Format}