
\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
