\section{Introduction}
\label{sec:intro}

Thermoelectric effects in solid state devices allow to convert heat exchanged with the environment into an electric current.
Nanoscale conductors can this way work as on-chip converters of waste heat into power at low temperatures~\cite{benenti:2017}. To this end, a mechanism that breaks particle-hole symmetry is needed. Among the desired properties of a good thermoelectric engine, the conductor should be electrically isolated from the thermal source, such that the absorbed heat and the generated charge currents are well separated as is the case, e.g., in a thermocouple. Despite the fact that superconductors are good thermal insulators and (obviously) good electrical conductors~\cite{mazza:2015}, they are rarely considered as components of thermoelectric generators due to their intrinsic electron-hole symmetry~\cite{ginzburg1991thermoelectric}.

Recent approaches to this problem exploit tunnel junctions between unequal superconducting electrodes using nonlinearities~\cite{marchegiani_nonlinear_2020,marchegiani_phase_2020,germanese_bipolar_2022} or single-electron transistors~\cite{donald,kamp_phase_2019,bauer_phase_2021}, hybrid normal-superconductor junctions including quantum dot energy filters~\cite{Cao2015,sanchez_cooling_2018,Hussein2019,Kirsanov2019,Tan2021,verma:2022}, spin-dependent scattering~\cite{machon_nonlocal_2013,giazotto_proposal_2014,ozaeta:2014,hwang:2016,hwang:2016b,kolenda:2016,marchegiani:2016,savander_thermoelectric_2020,keidel_ondemand_2020,Blasi2020a,blasi_nonlocal_2021}, photon-assisted tunneling~\cite{hofer:2016prb}, or interference phenomena~\cite{jacquod_coherent_2010,kalenkov_large_2017}. The gapped density of states has also been successfully used for cooling~\cite{Nahum1994Dec,Leivo1996Apr} and for power sources~\cite{Marin-Suarez:2022,pekola:2022}.

Mesoscopic versions of the thermocouple geometry are based on three terminal configurations, where one terminal acts as the thermal reservoir and the other two support the charge current. The coupling to the heat source must be such that no particle currents are injected in the conductor, at least on average. Different mechanisms have been proposed depending on the way carriers couple to the heat source. This coupling can be mediated by microscopic interactions, e.g, charge correlations~\cite{hotspots,sothmann:2012,thierschmann:2015,hartmann:2015,roche:2015,thierschmann_thermal_2015} and electron-boson interactions~\cite{rutten:2009,entin:2010,krause:2011,sothmann:2012epl,ruokola:2012,jiang:2012,bergenfeldt:2014,arrachea:2014,bosisio_nanowire_2016,yamamoto:2016,jiang:2017,dorsch:2020,dorsch_characterization_2021}, or due to electronic relaxation in hot probes~\cite{jordan:2013,mazza:2014,donsa:2014,jiang_enhancing_2014,jaliel:2019},
to mention a few quantum dot systems.
Of particular interest are heat engines based on capacitively coupled  dots~\cite{molenkamp:1995,chan:2002,hubel:2007,mcclure:2007}, where the separation of charge and heat currents is explicit: current flows in one conductor, that we call passive, with the other one being connected to the heat source via charge fluctuations~\cite{hotspots,sothmann:2012,thierschmann:2015,hartmann:2015,roche:2015,strasberg:2013,zhang:2016,dare:2017,walldorf:2017,strasberg:2018,mayrhofer:2021}. Reversely the charge current (in response to a voltage bias or to temperature differences) can be used to operate the system as a refrigerator~\cite{sanchez:2013,zhang:2015,koski:2015,sanchez_correlation_2017,erdman_absorption_2018,dare:2019}.

\begin{figure}[b]
\includegraphics[width=\linewidth]{fig1.pdf}
\caption{\label{fig:scheme}Hybrid normal-superconductor heat engine based on coupled quantum dots ($p$, for passive and $a$, for active). Heat injected from the hot reservoir (H) induces the transfer of either (a) a Cooper pair in the superconductor (S) via an Andreev reflection in the normal lead (N), or (b) a quasiparticle transferred between S and N. In both cases, an electric current is generated in the passive subsystem.
}
\end{figure}

Similar coupling schemes are relevant for mesoscopic Coulomb drag configurations~\cite{mortensen_coulomb_2001,khrapai:2006,gasser:2009,shinkai:2009,Moldoveanu2009Jul,sanchez:2010,hussein:2012,hussein:2015,bischoff:2015,kaasbjerg:2016,keller:2016,lim:2018,sierra:2019,takada:2021} which generate a current in the passive circuit by coupling it to another one (the drive) that is voltage biased. Differently from extended samples, that rely on momentum exchange~\cite{Narozhny2016May}, the mesoscopic drag is based on the exchange of energy. This has been emphasized in proposals of thermal drag currents where only heat flows in the drive system~\cite{whitney:2015,lu:2016,bhandari:2018,ben-abdallah:2019,sanchez:2019,berdanier:2019,Lu:2021,Xi:2021,wang:2022,idrisov:2022} and of absorption refrigerators~\cite{venturelli:2013}.
In a recent work~\cite{mojtaba}, we found that the interplay of charge fluctuations and Andreev reflection processes gives rise to a drag current when the passive system contains a superconducting electrode. These processes compete with a second mechanism due to single quasiparticle tunneling. Here we explore how these mechanisms can make the system work as a nonlocal thermoelectric heat engine. For this purpose, we consider a hybrid three-terminal configuration consisting of two Coulomb-coupled quantum dots ($a$, for active and $p$, for passive), as sketched in Fig.~\ref{fig:scheme}. The active dot is coupled to a hot normal terminal (H). The passive one is connected to one normal (N) and one superconducting (S) terminals, see Refs.~\cite{graber:2004,deacon:2010,deFranceschi:2010} for discussions of related experimental implementations. We analyze the regimes where the generated current is mediated by the two mentioned mechanisms [Andreev reflection, cf. Fig.~\ref{fig:scheme}(a), and correlated quasiparticles filtered by the superconducting density of states, cf. Fig.~\ref{fig:scheme}(b)]. They are shown to give opposite contributions, which  in the intermediate regime where the two processes coexist reduces the thermoelectric performance of the engine (in terms of the generated power and efficiency).

Additionally it is interesting to explore the properties of the heat currents, which may make the device work as a refrigerator (when heat is extracted from the coldest reservoir) or as a heat pump (when heat flows into the hottest one). Multitask operations in three terminal conductors have been recently identified~\cite{manzano:2020} where combinations of two or more such operations are met, see also Refs.~\cite{entin:2015,hammam_2022}.

The remaining of the manuscript is organized as follows. The system and the model are described in Sec.~\ref{sec:model}. The limiting regimes where Andreev and quasiparticle transport dominate are discussed in Sec.~\ref{sec:andreev} and \ref{sec:quasipart} respectively, by using a master equation approach~\cite{Schaller2014,Strasberg2022} that provides a physical understanding of the involved mechanisms. The two regimes require different master equations which will be discussed separately. The intermediate regime is explored numerically in Sec.~\ref{sec:keldysh} by invoking a non-equilibrium Green's functions technique~\cite{Haug2008}. Finally, conclusions are presented in Sec.~\ref{sec:conclusions}.

\section{Hybrid coupled quantum dot system}
\label{sec:model}

\subsection{Hamiltonian}

We aim at the simplest description of the main mechanisms involved in the passive current generation. Spin is essential for pairing in the passive system. In the active dot, only charge fluctuations are required in order to generate a passive current~\cite{sierra:2019}.  Therefore, we ignore the spin degree of freedom and intradot Coulomb interactions in the active system for simplicity. This provides a good description of transport (up to spin degeneracy prefactors) as long as double occupancy of the active dot is negligible and in the absence of interdot spin-spin interactions. Note however that intradot Coulomb interactions in the passive dot do affect the pairing processes and need to be taken into account.

Our system is hence described with the Hamiltonian
\begin{align}\label{eq_H}
{\cal H}={\cal H}_{l}
+{\cal H}_{{\rm {dqd}}}+{\cal H}_{t}.
\end{align}
The first term is the Hamiltonian for the reservoirs,
\begin{gather}
\begin{aligned}
{\cal H}_{l}&= \sum_{k,\sigma}\varepsilon_{k,S}\hat{c}_{k,S,\sigma}^{\dagger}\hat{c}_{k,S,\sigma}+
\sum_k\Delta(\hat{c}_{k,S,\uparrow}^{\dagger}\hat{c}_{k,S,\downarrow}^{\dagger}+\rm{h.c.})\\
&+\sum_{k}\varepsilon_{k,H}\hat{c}_{k,H}^{\dagger}\hat{c}_{k,H}+
\sum_{k,\sigma}\varepsilon_{k,N}\hat{c}_{k,N,\sigma}^{\dagger}\hat{c}_{k,N,\sigma},
\label{eq:Hleads}
\end{aligned}
\end{gather}
where $\hat{c}_{k,\beta,\sigma}$ is the annihilation operator
for electrons with energy $\varepsilon_{k}$, momentum $k$
and spin $\sigma{=}\{\uparrow,\downarrow\}$ in terminal $\beta$=H,N,S (note that for $\beta$=H we drop the spin index),
and $\Delta$ is the order parameter in the superconducting electrode. Each of the metallic reservoirs is in local thermodynamic equilibrium with electrochemical potential $\mu_{\beta}$ and temperature $T_{\beta}$. The superconducting chemical potential $\mu_S$ sets the common Fermi energy, which we take from now on as the reference energy, $\mu_S=0$.
The passive subsystem is held at a lower temperature
($T\equiv T_N=T_S$) than that of the active subsystem ($T<T_H$).

The second term in Eq.~\eqref{eq_H} accounts for the
double quantum dot,
\begin{equation}\label{eq_dqd}
{\cal H}_{{\rm dqd}}=\underset{\alpha}{\sum}\varepsilon_{\alpha}\hat{n}_{\alpha}+U_{p}\hat{n}_{p,\uparrow}\hat{n}_{p,\downarrow}+U_{ap}\hat{n}_{a}\hat{n}_{p},
\end{equation}
where the number operators are defined as $\hat{n}_{p}=\sum_{\sigma}\hat{n}_{p,\sigma}=\sum_{\sigma}\hat{d}_{p,\sigma}^{\dagger}\hat{d}_{p,\sigma}$ for the passive dot,
and $\hat{n}_{a}=\hat{d}_{a}^{\dagger}\hat{d}_{a}$ for the active one. Here, $\hat{d}_{p,\sigma}$ and $\hat{d}_{a}$ denote  electron annihilation operators and $\varepsilon_{\alpha}$ are the dot energy levels. The interdot charging energy is $U_{ap}$, and $U_p$ is the (intradot) charging energy of the passive dot.
The coupling between the dots and the leads is given by the term:
\begin{gather}
\begin{aligned}
\label{eq_Ht}
{\cal H}_{t} &= \sum_{k,\sigma}\left(t_{N}\hat{d}_{p,\sigma}^{\dagger}\hat{c}_{k,N,\sigma}+t_{S}\hat{d}_{p,\sigma}^{\dagger}\hat{c}_{k,S,\sigma}+\rm{h.c.}\right)\\
&+ \sum_{k}\left(t_{H}\hat{d}_{a}^{\dagger}\hat{c}_{k,H}+\rm{h.c.}\right)
,
\end{aligned}
\end{gather}
where the $t_\beta$ are the dot-lead tunnel couplings. In the following, we consider the
wide-band approximation, in which case the tunnel hybridization strength is given by $\Gamma_{\beta}=2\pi|t_{\beta}|^{2}\rho_{0}^{\beta}$,
where $\rho_{0}^{\beta}$ is the corresponding electrode's density
of states in its normal state.

\subsection{Power and efficiency}

Let $J_\beta$ be the heat flux and $I_\beta$ the charge current in terminal $\beta$. While charge conservation ensures $I\equiv I_N=-I_S$ and $I_H=0$, energy conservation in the conductor is only expressed as $J_N+J_S=P-J_H$, where \begin{equation}
\label{eq:power}
P=-VI
\end{equation}
is the power dissipated (Joule heating) by a charge flowing in favor of a voltage bias $V=(\mu_N-\mu_S)/e$. In that case, $P<0$.

In the absence of additional forces ($V=0$), the difference $T_H-T>0$ causes heat to be transferred from the active to the passive dot via the Coulomb coupling $U_{ap}$~\cite{hotspots}. A nonlocal thermoelectric engine is able to convert this injected heat flow into a charge current generated in the passive
subsystem.
Then, a finite voltage
$V$ can be applied to counteract $I$.
In the range of $V$ such that $P$ is positive, work can be done against an external load. The ratio
between the output power in the passive circuit
and the heat absorbed from the active one,
\begin{equation}\label{eq_eff}
    \eta = \frac{P}{J_H}
\end{equation}
is the thermodynamic efficiency of the engine, which is limited
by the Carnot bound, $\eta\le \eta_C=1-T/T_H$~\cite{benenti:2017}. For voltages larger than the stall voltage (a non-local analogue of the thermovoltage, where the thermoelectric current is compensated by electrons flowing in favor of the bias), power is dissipated and the engine stops working.
It is then necessary to calculate $P$, $J_H$ and $\eta$ to fully
describe the heat engine properties of our device.

Sign changes of the heat currents define the other thermodynamic operations. When positive, $J_N$ and $J_S$ will characterise the cooling power of the system working as a refrigerator of the corresponding reservoir. A heat pump occurs when $J_H<0$.

\section{Andreev regime}
\label{sec:andreev}

\begin{figure}[t]
\includegraphics[width=\linewidth]{fig2.pdf}
\caption{\label{fig:andreevtransit}(a) Eigenstates of the proximitized quantum dot. The arrows represent the different transitions due to tunneling in the passive (black) and active systems (purple arrows). Note that the latest cannot change the parity of the passive system.
(b) Separation of the even and odd states, $\delta E_{\pm,n}{\equiv} E_{\pm,n}-E_{\sigma,n}$, as a function of the pairing, $\Gamma_S$, with $U_p=2U_{ap}=\unit[1.2]{meV}$ and $\varepsilon_a=0$.}
\end{figure}

When the superconducting gap $\Delta$ is the largest energy scale of the problem, one can replace the superconductor and its coupling to the passive dot with a pairing term in Eq.~\eqref{eq_dqd}:
\begin{equation}
    {\cal H}_{\rm dqd} \to {\cal H}_{\rm dqd} + \Gamma_{S}(\hat{d}_{p,\uparrow}^{\dagger}\hat{d}_{p,\downarrow}^{\dagger}+{\rm {h.c.})},
\end{equation}
see e.g., Refs.~\cite{rozhkov:2000,Trocha:2014,Nigg:2015} for microscopic justifications of this approximation.
Then, ${\cal H}_{\rm dqd}$ can be exactly diagonalized
yielding the eigenenergies
\begin{align}
    \label{eq:Esig}
     E_{\sigma,n}&=\varepsilon_{p}+n(\varepsilon_{a}+U_{ap}) \\
    \label{eq:Epm}
     E_{\pm,n}&=n\varepsilon_{a}+A_{\mp,n},
\end{align}
with the associated eigenstates
\begin{align}
    \{\left|\sigma,n\right\rangle& \,,
    \left|\pm,n\right\rangle  =\mathcal{N}_{\pm,n}^{-1}(A_{\pm,n}\left|0,n\right\rangle -\Gamma_{S}\left|2,n\right\rangle )\},
    \label{eq:evenst}
\end{align}
represented in Fig.~\ref{fig:andreevtransit}(a). Here, $n=0,1$ keeps record of the charge number in the active dot and we have defined
\begin{equation}
A_{\pm,n}=\tilde{\varepsilon}_{n}\pm\sqrt{\left(\tilde{\varepsilon}_{n}\right)^{2}+\Gamma_{S}^{2}},
\end{equation}
with $\tilde{\varepsilon}_{n}=\varepsilon_{p}+U_{p}/2+nU_{ap}$.
The states in the proximitized dot are expressed
in terms of the charge basis $\{0,\sigma,2\}$ and fall into
either the odd ($\left|\sigma\right\rangle$)
or the even ($\left|\pm\right\rangle$) charge sector.
The latter are important because they involve
coherent superpositions of states with 0 or 2 electrons,
and consequently contribute to the transfer of Cooper
pairs through the dot.

For small lead-dot couplings, a master equation
approach correctly describes the system dynamics in terms of sequential transitions between
the different states, as sketched in Fig.~\ref{fig:andreevtransit}(a). We get the stationary occupation of the different states, $P_\lambda$, by solving the set of equations
\begin{equation}
\label{eq:rateeq}
\sum_{\alpha\beta\kappa}\left(W_{\lambda\kappa}^{\alpha\beta}P_\kappa-W_{\kappa\lambda}^{\alpha\beta}P_\lambda\right)=0.
\end{equation}
The rates
$W_{\lambda\kappa}^{\alpha\beta}=\Gamma_{\lambda\kappa}^{\alpha\beta}+\gamma_{\lambda\kappa}^{\alpha\beta}$
for the transition $|\kappa\rangle\to|\lambda\rangle$
due to the tunneling of an electron from (to) lead $\beta$
in (out of) dot $\alpha$ are given respectively by
\begin{align}
    \Gamma_{\lambda\kappa}^{\alpha\beta}&=\Gamma_{\beta}|\langle\lambda|\hat{\delta}_{\alpha}^{\dagger}|\kappa\rangle|^{2}f_{\beta}\left(E_{\lambda}-E_{\kappa}\right)\\
    \gamma_{\lambda\kappa}^{\alpha\beta}&=\Gamma_{\beta}\bigl|\bigl\langle\lambda\bigl|\hat{\delta}_{\alpha}\bigr|\kappa\bigr\rangle\bigr|^{2}[1-f_{\beta}\left(E_{\kappa}-E_{\lambda}\right)],
\end{align}
where $f_\beta(E)=\{1+\exp[(E{-}\mu_\beta)/k_\text{B}T_\beta]\}^{-1}$
is the Fermi function. The operator $\hat\delta_\alpha$ annihilates an electron in dot $\alpha$ i.e., $\hat\delta_a=\hat d_a$ and $\hat \delta_p=\sum_\sigma\hat d_{p\sigma}$. These rates allow us to
determine both the charge current,
\begin{equation}\label{eq_I}
    I=e\underset{\lambda,\kappa}{\sum}\left(\gamma_{\lambda\kappa}^{p N}-\Gamma_{\lambda\kappa}^{p N}\right)P_\kappa,
\end{equation}
and the heat flux out of terminal $\beta$,
\begin{equation}\label{eq_J}
    J_{\beta} =\underset{\lambda,\kappa}{\sum}\left(E_{\lambda}-E_{\kappa}-\mu_\beta\right)\left(\gamma_{\kappa\lambda}^{\alpha \beta}P_\lambda-\Gamma_{\lambda\kappa}^{\alpha\beta}P_\kappa\right).
\end{equation}

\begin{figure}[t]
    \includegraphics[width=\linewidth]{fig3.pdf}
    \caption{Charge current in the passive subsystem as a function of active and passive dot levels for (a)
    $\Gamma_S=\unit[0.01]{meV}$
    and (b) $\Gamma_S=\unit[0.5]{meV}$. Parameters (in meV): $k_B T_H=1.5$, $k_B T_C=1$, $\Gamma_H=\Gamma_N=0.01$,
    $U_{ap}=0.6$, $U_p=1.2$ and $V=0$. }
    \label{fig_curAndreev}
\end{figure}

The transitions between odd and even parity states can be both due to an electron or to a hole tunneling process. The relative rate at which one or the other process contributes more depends on the ratio $A_{\pm,n}/\Gamma_S$ that dictates the asymmetry of the even superpositions in Eq.~\eqref{eq:evenst}. It depends on the energy of the passive dot and, in particular, on the occupation of the active one, $n$.
In fact, the gap between even and odd
states is itself a function of $n$, as shown in Fig.~\ref{fig:andreevtransit}(b).
Hence, fluctuations of the charge in the active dot change the contribution of electron- or hole-like processes in the passive dot in a dynamical way. This results in a rectification effect that leads to a finite $I$ at $V=0$.

In Fig.~\ref{fig_curAndreev} we plot the thermoelectric current in the passive circuit as a function of the energies of the two dot levels, which can be shifted with external gate voltages~\cite{vanderwiel:2002}.
A finite current is generated in the passive circuit due to Coulomb interactions with the hotter active circuit.
Similarly to Ref.~\cite{mojtaba},
the passive system response occurs around the center of the stability diagram, where fluctuations of the charge of both dots are enhanced. However, while the drag current in Ref.~\cite{mojtaba} is driven by nonequilibrium fluctuations created by a dc bias applied across the active dot, the generated current here is driven by purely thermal means.

\begin{figure}
    \includegraphics[width=\linewidth]{fig4.pdf}
    \caption{(a-b) Zero bias charge and (c-d) heat current from H for different couplings to the superconductor: $\Gamma_S=0.01$, $0.2$, $0.5$, $0.7$, and $\Gamma_S=\unit[0.9]{meV}$. The different panels show the dependence of the currents as a function of both the active (a), (c) and passive (b), (d) dot levels, along the dashed lines in Fig.~\ref{fig_curAndreev}(a) and (b). The current in panel (b) is divided by 10 to fit with the axis of panel (a). Other parameters are as in Fig.~\ref{fig_curAndreev}}
    \label{fig_heatAndreev}
\end{figure}

For small coupling to the S electrode,
the current changes sign by tuning $\varepsilon_p$, see Fig.~\ref{fig_curAndreev}(a). The current direction is determined by the character (electron- or hole-like)
of the dominant transitions, which is strongly dependent on the position of the passive dot energy.
The generated current can be seen at a much wider range of gate voltages when $\Gamma_S$ is further increased, see Fig.~\ref{fig_curAndreev}(b).
An important aspect is that the generated current changes sign in this case also when tuning the level position of the active dot, keeping $\varepsilon_p$ fixed, as can be seen in Fig.~\ref{fig_heatAndreev}(a). This way, both the magnitude and the sign of the passive current can be controlled by using the parameters of the active system as an external knob.
Linecuts of the currents along the dotted lines in Fig.~\ref{fig_curAndreev} are plotted in Fig.~\ref{fig_heatAndreev}(a) and (b) for clarity.

Let us now investigate the heat currents.
To make a comparison, we show in Figs.~\ref{fig_heatAndreev}(c) and (d) the heat current $J_H=-J_N$  together with the charge current
through the passive dot. While $I$ has the characteristic
sign changes discussed above, $J_\beta$ does not change sign, as expected for heat between a hot and a cold reservoir in the absence of external work done on the system.
This energy is transferred through the
interdot repulsion $U_{ap}$ into the passive subsystem in the region where the charge states of the two quantum dots fluctuate strongly.
We recall that in the infinite gap approximation
the superconducting reservoir is only treated
as a source of pair correlations for the $p$ dot
and therefore we have two heat fluxes only.
In the next sections, we will relax this approximation
and study also the heat that flows in the superconductor.

Furthermore, the charge and heat currents in Fig.~\ref{fig_heatAndreev} show a nonmonotonic behaviour of the maximal charge current with the coupling to the superconductor: for small values of $\Gamma_S$, the generated current increases with the coupling to S, as expected. However, as it becomes comparable to and larger than the charging energy, both $I$ and $J_H$ get reduced. This is understood in terms of the increased splitting between the even and odd parity states, see Fig.~\ref{fig:andreevtransit}(b). For large $\Gamma_S$, the fluctuations between the lowest energy $|+,n\rangle$ state and the $|\sigma,n\rangle$ states (due to either the tunneling of an electron or a hole) are hence suppressed, as soon as $E_{\sigma,n}{-}E_{+,n}>k_\text{B}T$, in a similar way as the charging energy suppresses single-electron transport~\cite{kouwenhoven:1997}. These fluctuations are necessary both for the charge and the heat currents, resulting in a backaction on the thermal transport between H and N. Hence, pairing cooperates with $U_{ap}$ in blocking the current. Differently from Coulomb blockade, this pairing blockade effect is quantum coherent.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{fig5.pdf}
    \caption{Thermoelectric performance of the heat engine in the Andreev regime for different values of $\varepsilon_a$. Power and efficiency are computed for voltages ranging between zero and the stall voltage, where $P=\eta=0$. Other parameters are as in Fig.~\ref{fig_curAndreev}(a).}
    \label{fig_effAndreev}
\end{figure}

We are now in a position to assess the thermoelectric efficiency $\eta$. Positive power $P$ is produced when the current flows against the applied voltage $V$. This occurs between $V=0$ and the stall voltage where the bias compensates the thermoelectric current (making $I=0$).
Thus, in our efficiency calculations we select the maximum negative current, which occurs around $\varepsilon_p=-(U_p+U_{ap})/2+U_{ap}/2$, as shown in Fig.~\ref{fig_heatAndreev}.
By tuning $V$, we calculate both $P$ and $\eta$, which both vanish at $V=0$ and at the corresponding stall voltage where $I=0$, see Eqs.~\eqref{eq:power} and \eqref{eq_eff}, as shown in Fig.~\ref{fig_effAndreev}.
We observe that power and efficiency are optimized for different values of $\varepsilon_a$. When the gates are such that current is maximal (around $\varepsilon_a=\unit[-0.6]{meV}$), the power is enhanced. However, the efficiency increases when both $\varepsilon_a$ and $\varepsilon_a+U_{ap}$ go well below $\mu_H$, at the expense of reducing the fluctuations and hence the generated power.
In order to have a sizable performance, the value of $\Gamma_S$ needs to be optimized as well. We need it to be of the order of $A_{\pm,n}$ in order to maximize the effect of the fluctuations of the pairing wavefunction.
In the limits where either  $\Gamma_S\gg\tilde{\varepsilon}_n$ or the opposite, the asymmetry between electron- and hole-like processes disappears and thereby the thermoelectric effect vanishes.

\section{Quasiparticle regime}
\label{sec:quasipart}

The other contribution to transport is due to quasiparticles. They are expected to dominate in different configurations: either when the passive dot level is close to the superconducting gap, when $\Gamma_S\rightarrow0$, or at high temperatures.
In these cases, tunneling of Cooper pairs
is strongly suppressed. In order to discriminate its contribution, we consider a simple model where the effect of the superconductor is only due to its density of states (DoS). The spin of the electron does not play a crucial role, so for simplicity we assume the limit $U_p\ggk_\text{B}T$ where a charge basis captures the main effect.

We are then left with four states $|n_p, n_a\rangle$, with 0 or 1 electron in the passive and active dots.
The hopping rates in the passive subsystem are
\begin{align}
\Gamma_{1n,0n}^{p,\beta}{=}\Gamma_{\beta}\nu_{\beta n}f_{\beta n}\ \ \text{and}\ \
\gamma_{0n,1n}^{p,\beta}{=}\Gamma_{\beta}\nu_{\beta n}(1{-}f_{\beta n}),
\end{align}
where $\nu_{Nn}=1$, and $\nu_{Sn}$ is the normalized superconductor DoS.
Notably, we take into account the dependence on the active
dot occupation $n=0,1$ via the Coulomb repulsion strength
$U_{ap}$. This is better seen in the Fermi functions as
$f_{\beta n}=f_{\beta}(nU_{ap}+\varepsilon_{p})$
for $\beta$=N,S.
We model the superconductor DoS with the Dynes
form~\cite{dynes:1978},
\begin{equation}
    \nu_{Sn} = \left|
    {\Re}\,\frac{\varepsilon_p+nU_{ap}+i\chi}{\sqrt{(\varepsilon_p+nU_{ap}+i\chi)^2-\Delta^2}}
    \right|,
\end{equation}
where $\chi$ takes into account quasiparticle occupation in the gap.
The rates in the active dot read
\begin{align}
    \Gamma_{n1,n0}^{a,H}=\Gamma_{H}f_{H n}\quad\text{and}\quad
\gamma_{n0,n1}^{a,H}=\Gamma_{H}(1-f_{H n}),
\end{align}
where $f_{H n}=f_{H}(nU_{ap}+\varepsilon_{a})$.
We can find the stationary occupation probabilities $P_{n_pn_a}$ by solving the system of equations~\eqref{eq:rateeq}. With these, we again calculate charge and heat currents using Eqs.~\eqref{eq_I} and~\eqref{eq_J}.

\begin{figure}
    \includegraphics[width=\linewidth]{fig6.pdf}
    \caption{Quasiparticle transport enabling sequences at zero bias. A finite current is induced by a temperature difference $T_H>T_N=T_S$ whose sign depends on the passive dot energy. Charge flows (a) from $N$ to $S$ for   $\varepsilon_p-\mu_N\in(\Delta{-}U_{ap},\Delta)$, and (b) in the opposite direction when $\varepsilon_p-\mu_N\in(-\Delta{-}U_{ap},-\Delta)$. The order of the tunneling events is indicated.}
    \label{fig_procQP}
\end{figure}

\begin{figure}
    \includegraphics[width=\linewidth]{fig7.pdf}
    \caption{Zero bias (a) charge and (b) heat currents in the different terminals in the quasiparticle regime. They are plotted as functions of the passive dot energy.
    Parameters: $\Delta=\unit[0.2]{meV}$, $k_B T_H=\unit[0.3]{meV} $, $k_B T=\unit[0.2]{meV}$,
    $U_{ap}=\unit[0.1]{meV}$, $\varepsilon_a=-U_{ap}/2$, $\Gamma_H=\Gamma_N=\Gamma_S=\unit[0.01]{meV}$, and $\chi=10^{-6}$.
    }
    \label{fig_heatQP}
\end{figure}

It was shown in Ref.~\cite{hotspots} (and later confirmed experimentally~\cite{thierschmann:2015}) that a similar configuration (based there on all-normal reservoirs) leads to a non-local thermoelectric response provided the tunneling rates in the passive system are energy-dependent and left-right asymmetric. In our hybrid configuration, this property is enabled by the superconducting DoS: tunneling is suppressed at energies $\varepsilon_p+nU_{ap}$ within the gap, and enhanced at energies close to the coherence peaks, as sketched in Fig.~\ref{fig_procQP} for two different level positions. It can be shown that in the case where the tunneling couplings $\Gamma_\beta$ are energy independent, a zero-bias current \begin{equation}
\label{eq:IqpV0}
I(V{=}0)\propto(\nu_{S0}-\nu_{S1})(f_{N0}-f_{N1})
\end{equation}
is generated~\cite{mojtaba}.
This current is shown in Fig.~\ref{fig_heatQP}(a). As expected, $I=0$ when both energies $\varepsilon_p+nU_{ap}$ lie in the gap [region with $\varepsilon_p-\mu_N\in(-\Delta,\Delta{-}U_{ap})$, marked as C].
The conditions for transport appear in regions of width given by the Coulomb repulsion $U_{ap}$, marked as B [with $\varepsilon_p-\mu_N\in(-\Delta{-}U_{ap},-\Delta)$] and D [with $\varepsilon_p-\mu_N\in(\Delta{-}U_{ap},\Delta)$] in Fig.~\ref{fig_heatQP}(a). There, electrons/holes around the Fermi energy in N (such that $\nu_{S0}=0$ at those energies) can exchange energy with the passive dot and be transferred over/below the gap (when $\nu_{S1}\neq0$). In this configuration, the charge current is antisymmetric around $\varepsilon_p=-U_{ap}/2$. Note that the change of sign around the particle-hole symmetry point is opposite to the one obtained in the Andreev regime discussed in Sec.~\ref{sec:andreev}. Also differently from that case, the charge current does not change sign with the position of the active dot, due to the lack of coherence in the passive dot.

\begin{figure}[t]
    \includegraphics[width=.9\linewidth]{fig8.pdf}
    \caption{Map of the different thermal operations as a function of the passive dot energy and the applied voltage. Different colors mark different thermodynamic operations, as indicated (white means no useful operation). They are labeled with E, when the system works as a heat engine, and with R$_\alpha$ and P$_\alpha$, when it works as a refrigerator or a heat pump in terminal $\alpha$. Parameters: $k_\text{B}T=\Delta$, $k_\text{B}T_H=1.1\Delta$, $\varepsilon_a=-0.3\Delta$, $\chi=10^{-4}$. Other parameters as in Fig.~\ref{fig_heatQP}.

    }
    \label{fig_multitask}
\end{figure}

The same transport windows are evident in the heat transport, as shown in Fig.~\ref{fig_heatQP}(b). In region C, where transport in terminal S is avoided by the gap, all the heat injected from H is absorbed by the normal contact. Out of this region, the heat current sign gives insights of the relevant mechanism for the generated current. In region B electrons below the gap in S enter the passive dot when the active one is empty and tunnel out below the chemical potential of N once the active dot has been occupied. This way, $I<0$, the superconductor is heated up ($J_S<0$) and the normal contact is refrigerated ($J_N>0$), as sketched in Fig.~\ref{fig_procQP}(b). In the opposite region D, electrons tunnel from over $\mu_N$ into the passive dot, with an empty active dot, and tunnel out over the superconducting gap after gaining an energy $U_{ap}$ from the interaction with an electron having occupied the passive dot, see Fig.~\ref{fig_procQP}(a). In this case, while the charge flows in the opposite direction, $I>0$, one still finds $J_S<0$ and $J_N>0$, as expected. Noticeably, we find that close to $\varepsilon_p-\mu_N\approx\pm(\Delta+U_{ap})$ (in the frontiers of A-B and D-E regions), the peaks of the superconducting DOS enable an additional change of sign in the charge current that cools the superconductor. In all these cases (except for region C), the system works as a hybrid thermal engine~\cite{manzano:2020} capable of simultaneously generate power and cool (either terminal N or S) from a single resource (heat in reservoir H). This possibility is further explored in Fig.~\ref{fig_multitask} in the presence of a finite voltage. This way, we find an additional hybrid operation where high enough voltages are able to reverse the heat current $J_H$, while $J_S>0$ (violet region in Fig.~\ref{fig_multitask}, labelled as R$_S$P$_H$).  Hence, the system works both as a refrigerator of S and a heat pump into the hottest reservoir, H. Tuning the voltage facilitates the cooling of the N or S terminals (light blue regions in Fig.~\ref{fig_multitask}), very much like a usual quantum dot Peltier refrigerator~\cite{edwards:1993,prance:2009}. Here these operations are affected by the coupling to the other dot. The configuration map in Fig.~\ref{fig_multitask} is repeated by inverting over the $V=0$ and $\varepsilon_p-\mu_N=-U_{ap}/2$.

\begin{figure}[t]
    \includegraphics[width=\linewidth]{fig9.pdf}
    \caption{Thermoelectric efficiency in the quasiparticle regime, for different interaction $U_{ap}$. Each point is optimized with respect to the value of $\varepsilon_p$. The inset zooms the curve for $U_{ap}=0.8\Delta$ in, for clarity. Other parameters as in Fig.~\ref{fig_heatQP}.
    }
    \label{fig_effPQP}
\end{figure}

The asymmetric filtering of some electronic tunneling events, here facilitated by the superconducting DoS, is expected to enable optimal heat to power conversion in terms of power and efficiency~\cite{hotspots}.
Consider for instance region D, and the scheme in Fig.~\ref{fig_procQP}(a): the gap avoids the electron tunneling from N to subsequently tunnel to S without having exchanged an energy $U_{ap}$ with the active system. Once this is done, the coherence peak makes it more favorable to tunnel to S than back to N. This way, most of the cycles for which an amount of energy $U_{ap}$ is transferred from the active to the passive system result in an electron being transferred from N to S.
This is confirmed in Fig.~\ref{fig_effPQP}, where power increases by two orders of magnitude with respect to the Andreev regime, with the efficiency being doubled. Furthermore, for low $T$, high efficiencies approaching $\eta_C$ are obtained at low power (but still one order of magnitude larger) by taking advantage of the BCS peak. The maximal efficiency is hence limited by the sharpness of the coherence peak, which is affected by the Dynes parameter, $\chi$. In this case, the thermoelectric performance improves when increasing $U_{ap}$, such that electrons gain enough energy from H to overcome the gap. Increasing temperature smears the effect of the spectral features and hence affects the performance.

\section{Intermediate regime}
\label{sec:keldysh}

In this section we explore the thermoelectric performance of the engine in the intermediate regime where the temperature is low and the superconducting energy gap $\Delta$ is comparable to other energy scales. Then it is important to take into account the contributions from both Andreev and quasiparticle transport mechanisms in the calculations. This is possible by employing the non-equilibrium Green's functions (NEGF) technique~\cite{mojtaba}, which furthermore takes the finite linewidth of the quantum dot states (which is neglected in the sequential tunneling master equations discussed in previous sections) into account.

Within this formalism, the different currents through terminal $\beta$=N,S,H, can be calculated by using
\begin{align}
I_{\beta}&=  \frac{e}{2\hbar}\int\frac{d\omega}{2\pi}\text{Tr}\left(\hat\sigma_z\cal{I}_{\beta}\right)\\
J_{\beta}&= \frac{1}{2h}\int\frac{d\omega}{2\pi}\text{Tr}\left[\left(\omega I_{2}-eV_{\beta}\hat{\sigma}_{z}\right){\cal I}_{\beta}\right],
\end{align}
for charge and heat, respectively, with
\begin{gather}
\label{eq:currop}
\begin{aligned}
\cal{I}_{\beta}= ~
G_{\alpha}^{R}\left(\omega\right)\Sigma_{\alpha,\beta}^{<}\left(\omega\right)+G_{\alpha}^{<}\left(\omega\right)\Sigma_{\alpha,\beta}^{A}\left(\omega\right
 \\
 &
 -\Sigma_{\alpha,\beta}^{R}\left(\omega\right)G_{\alpha}^{<}\left(\omega\right)
 -\Sigma_{\alpha,\beta}^{<}\left(\omega\right)G_{\alpha}^{A}\left(\omega\right)
\end{aligned}
\end{gather}
Here the terminal index $\beta$ fixes the index $\alpha=p,a$ of the quantum dot coupled to it, and
$I_2$ and $\hat\sigma_i$ are the $2{\times}2$ identity and Pauli matrices, respectively. The retarded and lesser selfenergies are given by $\Sigma_{a,H}^{R}\left(\omega\right)=-i\Gamma_{a,H}$, $\Sigma_{p,N}^{R}\left(\omega\right)=-i\Gamma_{p,N}I_{2},$ and
$\Sigma_{p,S}^{R}\left(\omega\right)=-i\Gamma_{S}b\left(\omega\right)(I_{2}-\Delta\omega^{-1}\hat\sigma_{x}),$ where
$b\left(\omega\right)\equiv\left|\omega\right|(\omega^{2}-\Delta^{2})^{-1/2}\theta\left(\left|\omega\right|-\Delta\right)-i\omega(\Delta^{2}-\omega^{2})^{-1/2}\theta\left(\Delta-\left|\omega\right|\right)$, with
the Heaviside step function, $\theta\left(E\right)$, and $\Sigma_{\alpha,\beta}^<(\omega)=-2i{\rm Im}(\Sigma_{\alpha,\beta}^R)f_\beta(\omega)$. The interacting retarded and lesser components of the Keldysh Green's function of dot $\alpha$, $G_{\alpha}^{R/<}$, are calculated self-consistently by using the Dyson and the Keldysh equations
\begin{align}
\label{eq:Gar_interacting}
G_{\alpha}^{R}&=\left\{ \left[g_{\alpha}^{R}\right]^{-1}-\Sigma_{\alpha,{\rm {int}}}^{R}\right\} ^{-1}\\%\end{equation}
G_{\alpha}^{<}&=G_{\alpha}^{R}\left(\Sigma_{\alpha,{\rm {leads}}}^{<}+\Sigma_{\alpha,{\rm {int}}}^{<}\right)G_{\alpha}^{A},\label{eq:Glpinteracting}
\end{align}
where $g_{\alpha}^{R}$ is the non-interacting retarded Green's function of dot $\alpha$, and $\Sigma_{\alpha,{\rm {int}}}$ is the corresponding interaction selfenergy of the dot $\alpha$ due to its Coulomb interaction with other charges in the dots. We refer the interested reader to the supplementary material in Ref.~\cite{mojtaba} for the explicit expressions of these self energies and Green's functions and also for the details of the numerical procedure used to solve the above equations selfconsistently. Here we chose the parameters such that the system falls in the intermediate regime where both Andreev and quasiparticle mechanisms contribute. In particular, their relative contribution is controlled by the coupling $\Gamma_S$.

\begin{figure}
    \includegraphics[width=\linewidth]{fig10.pdf}
    \caption{(a)-(d) Charge $I$ and heat currents $J_N$, $J_S$ and $J_H$ as functions of the active and passive dot energies for $\Gamma_S=3\Gamma_N$. (e)-(h) Line cuts of left panels at $\varepsilon_a=-(U_{ap}+U_p)/2=\unit[-0.15]{meV}$ (indicated by dashed lines in left panels) and for different values of $\Gamma_S=0.5\Gamma_N$, $\Gamma_N$, $3\Gamma_N$ and $8\Gamma_N$. Other parameters are $\Delta=\unit[0.2]{meV}$, $k_\text{B}T_H=2k_\text{B}T=\unit[0.2]{meV}$, $\Gamma_H=\Gamma_N=\unit[0.01]{meV}$, $U_p=1.6U_{ap}=\unit[0.16]{meV}$ and $V=0$.
    }
    \label{fig_NEGF-EaEp}
\end{figure}

As we learned in the previous sections, the contribution of the two different mechanisms (Cooper pair and quasiparticle dominated transport) depends strongly on the position of the quantum dot levels. The Andreev (pairing) induced current is finite around the particle-hole symmetric point, see Fig.~\ref{fig_curAndreev}. Differently, the quasiparticle contribution occurs close to the gap borders, see Fig.~\ref{fig_heatQP}. Hence, their competition results in a rich behaviour of the current generated at $V=0$. This is shown in Fig.~\ref{fig_NEGF-EaEp}(a): in the center of the stability diagram the Andreev contribution dominates, with the current sign being as predicted by the rate equation description in Sec.~\ref{sec:andreev}. As the passive dot level deviates from this region, the quasiparticle contribution starts to dominate and the current changes sign, in agreement with the results in Sec.~\ref{sec:quasipart}.
Importantly, for fixed values of $\varepsilon_p$, the charge current changes sign as a function of $\varepsilon_a$. This effect was found in the Andreev regime, cf. Figs.~\ref{fig_curAndreev}(b) and \ref{fig_heatAndreev}(a), and is hence a signature of quantum coherence persisting in the intermediate regime. This is interesting since it allows us to control the generated current by only using the active dot parameters: $T_H$ generates the current and $\varepsilon_a$ controls its direction.

The different contributions are also visible in the heat currents, plotted in Figs.~\ref{fig_NEGF-EaEp}(b)-(d). As expected, the heat current injected from H is always positive (for $V=0$). It is restricted to the region where charge fluctuations in both dots are present~\cite{koski:2015}, see Fig.~\ref{fig_NEGF-EaEp}(d), emphasizing the role of electron-electron interaction in the heat exchange between the systems. The latter occurs as well in the passive system currents, with two particular aspects to point out: heat always flows into the superconductor ($J_S<0$) but is slightly suppressed by the gap for $\varepsilon_p$ close to the particle-hole symmetry point. Most interestingly, $J_N$ changes sign at the crossover between the Andreev and the quasiparticle regimes, leading to cooling, see Fig.~\ref{fig_NEGF-EaEp}(b). Hence, we recover the ER$_N$ operation shown in Fig.~\ref{fig_multitask} as soon as a small voltage $V$ is applied.

The relative contribution of the two mechanisms depends strongly on the parameters. Crucially,  it depends on the coupling to the superconductor, because the Andreev reflection is suppressed for low $\Gamma_S$. This is shown in Figs.~\ref{fig_NEGF-EaEp}(e)-(h), where cuts of the left-column panels at $\varepsilon_a=-(U_{ap}+U_p)/2$ (marked by dotted lines) are plotted for different $\Gamma_S$. Remarkably, Fig.~\ref{fig_NEGF-EaEp}(e) shows how the Andreev contribution to the charge current dominates for large $\Gamma_S$ and totally disappears for sufficiently small couplings. Note that the quasiparticle contribution gets correspondingly reduced but can still be visible for large $\Gamma_S$ by the effect of the gap border. In this case, the sharp features in Fig.~\ref{fig_heatQP} due to the gap are smeared here due to the finite level-width of the passive dot. We observe a different behaviour in Fig.~\ref{fig_NEGF-EaEp}(f), where the contribution of the strong gap features to the cooling of N are favored by increasing the coupling to the superconductor.
While $J_S$ gets strongly reduced for low $\Gamma_S$, cf. Fig.~\ref{fig_NEGF-EaEp}(g), the injected current from the hot bath is barely unaffected, except for large couplings, see Fig.~\ref{fig_NEGF-EaEp}(h).
The latter can be attributed to the pairing blockade effect introduced in Sec.~\ref{sec:andreev}. Note however a stronger nonmonotonic behaviour of $J_S$ in Fig.~\ref{fig_NEGF-EaEp}(g), that is caused by the increasing role  at high $\Gamma_S$ of pair tunneling, which does not carry heat.

\begin{figure}[t]
    \includegraphics[width=\linewidth]{fig11.pdf}
    \caption{(a) Charge current, $I$, and (b) heat current, $J_N$, as functions of the interdot Coulomb interaction strength $U_{ap}$ and the passive dot energy for $\Gamma_S=3\Gamma_N$. $\varepsilon_a=\unit[-0.15]{meV}$ and other parameters are as in Fig.~\ref{fig_NEGF-EaEp}.   }
    \label{fig_NEGF-UapEp}
\end{figure}

The interplay of the two processes furthermore depends on the interaction between the two quantum dots, $U_{ap}$. This is shown in Fig.~\ref{fig_NEGF-UapEp}. The quasiparticle contribution to the charge current is maximal when $U_{ap}\lesssim\Delta$ i.e., when quasiparticles gain/lose enough energy from the interaction with the active dot to overcome the gap, see Fig.~\ref{fig_NEGF-UapEp}(a). For the same reason, the cooling of terminal N, due to quasiparticles, is suppressed for $U_{ap}>\Delta$, see Fig.~\ref{fig_NEGF-UapEp}(b). The Andreev contribution is on the contrary reduced by the interaction, as it reduces the hybridization of the even parity states (essential for the current generation via pairing) when it becomes large compared to $\Gamma_S$, see Eq.~\eqref{eq:evenst}. For large $U_{ap}$, the Andreev contribution is washed out by quasiparticles, as expected.

\begin{figure}[t]
    \includegraphics[width=\linewidth]{fig12.pdf}
    \caption{(a) Power and (b) normal terminal heat current in the low voltage regime as functions of voltage and the passive dot energy. Panels (c)-(e) show line-cuts of the power and the different heat currents for different $\Gamma_S$ and two different voltages indicated by dotted lines in (a) and (b). Same parameters as in Fig.~\ref{fig_NEGF-EaEp}, with $\varepsilon_a=\unit[-0.15]{meV}$.}
    \label{fig_NEGF-power-VbEp}
\end{figure}

The generated current leads to useful power in the presence of a finite bias voltage. This affects the properties of the heat currents as well. In Fig.~\ref{fig_NEGF-power-VbEp} we show the generated power and the heat current in N. Figures~\ref{fig_NEGF-power-VbEp}(a) and (b) show the low voltage behaviour, with line-cuts for fixed $V$ plotted in Figs.~\ref{fig_NEGF-power-VbEp}(c)-(f). In the low (but positive) voltage regime, finite power is generated in two regions as we tune $\varepsilon_p$, which correspond to the quasiparticle and Andreev dominated regimes. The quasiparticle dominated region (around $\varepsilon_p\approx\unit[-0.4]{meV}$) results in larger power and larger stall voltages for this intermediate configuration. The cooling of terminal N discussed above coexists with the power production, though the limiting voltage for having $J_N>0$ is smaller than the stall voltage for $P>0$, compare Figs.~\ref{fig_NEGF-power-VbEp}(a) and (b). This is consistent with the behaviour shown in Fig.~\ref{fig_multitask}, where the E and ER$_N$ are adjacent for quasiparticle induced transport. For the configuration chosen here, the efficiency is reduced to $\eta\approx0.02\eta_C$, which we attribute to the finite width of the passive level and to the fact that the competition of the two mechanisms unavoidably reduces the generated power, without necessarily reducing the transferred heat $J_H$.

The line-cuts in Figs.~\ref{fig_NEGF-power-VbEp}(c)-(f) show $P$ and $J_\beta$, for different $\Gamma_S$. The same behaviour that we observed for the current is shown here for the power: as we decrease the coupling to S, the quasiparticle contribution dominates and give $P>0$ even close to the particle-hole symmetric point. The injected current from H is little affected. However, we observe that in the region with Cooper pair contributions ($\varepsilon_p$ around the particle-hole symmetry point), $J_H$ decreases with $\Gamma_S$, while it increases as quasiparticles take over (for larger $|\varepsilon_p|$), confirming the coherent nature of the reduction of $J_H$ in the gap, in agreement with the pairing backaction effect identified in Sec.~\ref{sec:andreev}. Note however that here we are in the regime $k_\text{B}T>\Gamma_S$.

Decreasing $\Gamma_S$ suppresses the ER$_N$ regions, as seen by comparing Figs.~\ref{fig_NEGF-power-VbEp}(c) and (e) for $\varepsilon_p\approx\unit[-0.4]{meV}$, as expected for a quasiparticle-only effect. For even more negative $\varepsilon_p$, we find cooling of the superconductor induced by $V$, for small $\Gamma_S$, see Figs.~\ref{fig_NEGF-power-VbEp}(d) and (f). However, the operation ER$_S$ is not found, which we attribute to the effect of the finite linewidth of the quantum dot states. We speculate that $P>0$ and $J_S>0$ may coexist for even smaller $\Gamma_N$ and $\Gamma_S$. Unfortunately this will also considerably reduce the heat flows.

\begin{figure}
    \includegraphics[width=\linewidth]{fig13.pdf}
    \caption{Plots of (a) power, (b) $J_N$, (c) $J_S$ and (d) $J_H$, as functions of the bias voltage $V$ and the passive dot energy $\varepsilon_p$. Arrows in panel (a) indicate the regions where $P>0$, shown in Fig.~\ref{fig_NEGF-power-VbEp}(a) for clarity. The same configuration as in Fig.~\ref{fig_NEGF-power-VbEp} is considered, for a larger range of $V$. Note that the regions of $P>0$ are almost invisible in panel (a) for being restricted to very low voltages.
    }
    \label{fig_NEGF-VbEp}
\end{figure}

The currents in the presence of larger voltages are shown in Fig.~\ref{fig_NEGF-VbEp}. There, the charge current flowing in favor of the bias leads to power dissipation, $P<0$. The arrows in Fig.~\ref{fig_NEGF-VbEp}(a) mark the regions where power is produced, as zoomed in in Fig.~\ref{fig_NEGF-power-VbEp}(a). Dissipation leads to heating terminal N, see Fig.~\ref{fig_NEGF-power-VbEp}(b), except for large and positive $\varepsilon_p$, where electrons tunnel from N into the passive dot over $\mu_N$. Voltage also allows for the heat extracted from the superconductor to increase, see Fig.~\ref{fig_NEGF-VbEp}(c), again due to quasiparticles (we remind that Cooper pairs do not lead to heat flows, as discussed in Sec.~\ref{sec:andreev}). Note the cooling power $J_S$ is maximal for voltages of the order of the gap. This process is related to refrigerators based on NIS tunnel junctions~\cite{giazotto:2006}, here mediated by the fluctuations in a quantum dot, see also Ref.~\cite{sanchez_correlation_2017}. Additionally, the dissipated Joule power reverts the heat flow at terminal H close to the gap borders, see Fig.~\ref{fig_NEGF-VbEp}(d) for $\varepsilon_p\approx\unit[-0.4]{meV}$  and $\varepsilon_p\gtrsim\unit[0.2]{meV}$, which makes the system work as a heat pump and recovers the R$_S$P$_H$ operation discussed in Sec.~\ref{sec:quasipart}. Remarkably, for high enough $V$, we observe a subgap region (around $\varepsilon_p\approx\unit[-0.1]{meV}$) where $J_H<0$, i.e., induced by Andreev processes and their correlation with the charge fluctuations of the active dot.

\section{Conclusions}
\label{sec:conclusions}

We have shown how a hybrid conductor composed by a quantum dot coupled to a normal and a superconductor contact can work as a quantum heat engine, when  capacitively coupled to a second quantum dot in contact with a hot terminal. The injection of heat is mediated by the interaction of electrons in different dots, whose charge fluctuates. The asymmetry of the normal-superconductor contact induces a zero-bias current via two competing mechanisms: on one hand, the superconducting density of states acts as an efficient energy filter for quasiparticle tunneling. On the other hand, the hybridization of even parity states due to pairing depends on the occupation of the hot dot, hence changing the relative contribution of empty and doubly occupied states to the superposition and affecting the electron- or holelike character of the electrons tunnelling from N. This way, the symmetry of a wavefunction is used as the piston of an autonomous quantum heat engine.

We have analyzed the limiting cases where these two contributions operate separately by using simple models described by quantum master equations. More realistic configurations where both effects coexist are treated numerically with NEGF techniques. The relative contribution of the two mechanisms can be controlled with gate voltages acting on the quantum dot levels. The contribution of Andreev processes dominates for states close to the chemical potential of the superconductor, while quasiparticle processes are enhanced for transitions via states aligned with the gap. The opposite contributions of the two mechanisms lead to changes of sign in the generated zero-bias current as the passive dot level is tuned which are not present in all-normal configurations~\cite{dare:2017,walldorf:2017}. The quantum coherent character of this mechanism manifests in an additional change of sign of the generated current as the active dot level is tuned. In this way, the coupled quantum dot acts not only as the heat source but also as an external knob of the response. We also identify a backaction of pairing on the active system in the form of a coherent pairing-Coulomb blockade of the injected heat current.
The gap and the coherence peak in the superconducting DoS are beneficial for the conversion of heat into power, which reaches high efficiencies in the quasiparticle regime. Furthermore, it allows the system to work as a hybrid engine where several thermodynamic operations (power generation, cooling and/or heat pumping) coexist. Most of these features survive in the intermediate regime, where the finite linewidth of the quantum dot states introduces a limitation to the performance of the device. The robustness of the Andreev induced currents in this regime makes this system appealing for the experimental realization of quantum themodynamic engines.

\acknowledgements
We acknowledge P. Burset and G. Steffensen for useful comments. Work funded by Iran Science Elites Federation (ISEF) and by Spanish State Research Agency through Grants Nos.\ PID2019-110125GB-I00, PID2020-117347GB-I00, PID2020-117671GB-I00 and RYC2016-20778 and through the Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence in R\&D (Nos.\ MDM2017-0711 and CEX2018-000805-M) funded by MCIN/AEI/10.13039/501100011033.


\section{Introduction}

Recent years have witnessed a revival of random channel access theory and techniques, fostered by new applications such as grant-free access schemes in \ac{MMA} scenarios \cite{Chen2021:Massive}.
This renewed interest has brought back the problem, known for decades (e.g., \cite{wolf1981:coding,Gallager1985:Perspective}), of framing random access under an information-theoretic setting able to capture not only collisions and intermittent users' activity but also the physical communication process.
Several research efforts have been made into this direction over the years, among which we mention \cite{Mathys1990:Class,Bar-David1993:Forward,Luo2012:New,Chen2017:Capacity,Recep2021:Random}.
Moreover, many recent works in this area have considered an ``unsourced'' model \cite{Polyanskiy2017:Perspective}, in which all users have the same codebook, the receiver is only tasked with producing an unordered list of messages, and the main performance metric is the per-user error probability \cite{Ordentlich2017:low_complexity,Ngo2021:Random_user_activity,vem2019:unsourced,Fengler2021:SPARCs, Liva2021:Coding}.

Coded random access has recently emerged as a new approach to reliably support \ac{MMA} applications, featuring a combination of packet-based coding and \ac{SIC} and having a strong connection with codes on sparse graphs.
Instances of coded random access are \ac{CSA} \cite{paolini2015:csa} and \ac{IRSA} \cite{liva2011:irsa}; in the latter one, packet repetitions are used as a form of packet erasure coding.

Notably, so far the above-mentioned efforts towards a new framework for random access have addressed coded random access only marginally, which motivates us to explore this direction of investigation.
Therefore, in this paper we propose a new approach to \ac{IRSA}, capturing information-theoretic concepts such as codebook construction and information rate.
The specific choice of the users' codebooks, based on \cite{Bar-David1993:Forward}, and the adopted adder channel mode allows incorporating in the analysis some known key quantities, such as the \ac{IRSA} asymptotic load threshold over the \ac{MPR} channel.
This allows us obtaining achievability and converse results on the average sum rate in the regime in which both the users' population size and the number of slots per frame grow unbounded, their ratio remaining constant.

\section{IRSA over MPR Channel}\label{sec:irsa_mpr}

\subsection{Irregular Repetition Slotted ALOHA}
\ac{IRSA} is a grant-free and uncoordinated access protocol for \ac{MMA}, strongly related to codes on sparse graphs \cite{liva2011:irsa}.
The protocol is grant-free as users perform transmissions with no prior handshake procedure with the receiver.
It is uncoordinated as users cannot cooperate with each other.
The time is slotted and framed, and users are slot- and frame-synchronous. The number of slots per frame is $N_{\mathrm{s}}$, the total number of users is $K$, the number of simultaneously contending users over the same frame, unknown to the receiver, is $K_{\mathrm{c}} \leq K$.

At the beginning of each frame some users are contending and some others are idle.
Each non-idle user contends for transmission of one packet over the current frame.
The user samples, independently of the other users, a discrete random variable $D \in \{2,\dots,\bar{d}\}$ with probability distribution $\Pr(D=d) = \Lambda_d$ and \ac{PGF}
$\Lambda(x) = \sum_{d} \Lambda_d x^d$.
The distribution $\Lambda$ is the same for all users and is hereafter referred to as the \ac{IRSA} distribution.
The user then draws, uniformly at random and without replacement, $d$ integers between $0$ and $N_{\mathrm{s}}-1$ and transmits $d$ replicas of its packet in the corresponding slots of the frame.
The status of the frame after transmissions can be described by means of a bipartite graph with $K$ user nodes, one per user, and $N_{\mathrm{s}}$ slot nodes, one per slot.
An edge connects user node $u_i$, $i\in\{0,\dots,K-1\}$, to slot node $s_j$, $j\in\{0,\dots,N_{\mathrm{s}}-1\}$, if user $i$ performed transmission of a packet replica in slot $j$.

\begin{example}
An example with $K = 11$ users ($u_0$, $\dots$, $u_{10}$) and $N_{\mathrm{s}} = 13$ slots ($s_0$, $\dots$, $s_{12}$) is depicted in Fig.~\ref{fig:frame}.
Users are represented by circles and slots by squares.
There are $K_{\mathrm{c}} = 7$ contending users, represented by colored circles; the number of transmitted packet replicas may change from user to user.
Blank circles represent users that are idle in the current frame.
\end{example}

\begin{figure}[t]
\centering
         \resizebox{0.50\textwidth}{!}{%
     		\input{./Figures/graph}
     	}%
       \caption{Graphical representation of \ac{IRSA} access with $K=11$ users ($K_{\mathrm{c}}=7$ contending) and $N_{\mathrm{s}}=13$ slots. Light-blue circles are contending users and blank circles idle users. Assuming $T=2$, light-blue squares are resolvable slots, light-red squares collision slots, and blank squares empty slots.}
    \label{fig:frame}
\end{figure}

The efficiency of \ac{IRSA}, denoted by $\eta$, is defined as $\eta = [\Lambda'(1)]^{-1}$ where $\Lambda'(1)$ represents the expected number of replicas transmitted by a contending user.
If users become contending independently of each other with probability $\pi$ at the beginning of each frame, the average load is
$\mathsf{G} = \pi K/N_{\mathrm{s}}\,\mathrm{[packets/slot]}
$.
The average number of packet replicas per slot is $\mathsf{G} / \eta$.

\subsection{\ac{SIC}-Based \ac{IRSA} Receiver over the \ac{MPR} Channel}

The $T$-\ac{MPR} channel is a packet-based channel model representing a simple extension of the collision channel \cite{Ghez1988:Stability}.
The $T$-\ac{MPR} channel model is based on three assumptions that may be summarized as follows:

\emph{Assumption~1}: The receiver can always discriminate between a slot with no packet arrivals, a slot with at most $T$ arrivals, and a slot with more than $T$ arrivals;

\emph{Assumption~2}: If there are more than $T$ arrivals in a slot, a condition detected by the receiver, then none of these packets can be correctly received;

\emph{Assumption~3}: If there are at most $T$ arrivals in a slot, a condition detected by the receiver, then all of these packets are correctly received with zero error probability.

Over the $T$-\ac{MPR} channel model, the \ac{IRSA} receiver attempts recovery of all packets using an iterative \ac{SIC}-based procedure.
At the beginning, the receiver marks all slots as empty, resolvable (at most $T$ arrivals), or collision (more than $T$ arrivals) slots.
Then, each iteration comprises two steps and may be described as follows, under the assumption that each replica carries information about the number and the positions of the other replicas:

1. In every resolvable slot, all packets are correctly received.

2. For each such packet, the interference of every replica of the packet is subtracted from the corresponding slot, possibly leading to new resolvable slots.

\begin{example}
With reference again to Fig.~\ref{fig:frame}, assume $T=2$.
Since slots $s_0$, $s_2$, $s_5$, and $s_7$, are resolvable, packet replicas transmitted by users $u_1$, $u_5$, $u_6$, and $u_{10}$ in these slots are correctly received at the first iteration.
Interference cancellation in the first iteration makes slots $s_4$ and $s_{12}$ resolvable.
In the second iteration, at least one packet replica transmitted by users $u_0$, $u_3$, and $u_8$ is correctly received in these slots.
Interference cancellation in the second iteration makes all slots empty and the procedure terminates successfully.
\end{example}

The \ac{SIC}-based procedure has been analyzed in the asymptotic setting $N_{\mathrm{s}} \rightarrow \infty$, $K \rightarrow \infty$, and $a=K/N_{\mathrm{s}}$ constant \cite{Ghanbarinejad2013:Irregular,stefanovic2018:multipacket}, under the following further assumption:

\emph{Assumption~4}: Interference subtraction across slots is ideal.

The analysis resembles density evolution of LDPC codes over the erasure channel and shows the existence of an asymptotic load threshold $\mathsf{G}^{*} = \mathsf{G}^{*}(\Lambda,T)$.
Let $N_{\mathrm{s}} \rightarrow \infty$, $K \rightarrow \infty$, and $a=K/N_{\mathrm{s}}$ be constant.
In this asymptotic setting, let $p_{\ell}$ be the probability that the generic edge is connected to a slot node not yet resolvable at the end of iteration $\ell$.
Then, the evolution of $p_{\ell}$ through \ac{SIC} iterations is governed by the recursion
\begin{align}
  p_{\ell} &= 1 - \exp\left( - \mathsf{G} \Lambda'(p_{\ell-1}) \right) \sum_{k = 0}^{T - 1} \frac{ 1 }{k!} \left( \mathsf{G} \Lambda'(p_{\ell-1}) \right)^{k}
\end{align}
with starting point $p_1 = 1-\exp(-\frac{\mathsf{G}}{\eta})\sum_{k=0}^{T-1} \frac{1}{k!} (\frac{\mathsf{G}}{\eta})^k$.
The load threshold is defined as
\begin{align}\label{eq:MPR_threshold}
    \mathsf{G}^{*} = \sup \{ \mathsf{G} > 0 : p_{\ell} \rightarrow 0 \text{ as } \ell \rightarrow \infty\} \, .
\end{align}
\begin{remark}\label{remark:plp}
The quantity $p_{\ell}^d$ represents the probability that the packet of a user, that transmitted $d$ packet replicas, has not yet been correctly received after $\ell$ \ac{SIC} iterations.
Thus the condition $p_{\ell} \rightarrow 0$, which holds for all $\mathsf{G} \leq \mathsf{G}^{*}(\Lambda,T)$, corresponds to a vanishing packet loss probability in the asymptotic setting.
\end{remark}

\section{A New Perspective on \ac{IRSA}} \label{sec:abc}

The description of \ac{IRSA} provided in Section~\ref{sec:irsa_mpr} is the usually adopted one and is typical of collision resolution approaches: it focuses on bursty arrivals and traffic, while it ignores any information-theoretic quantity.
In this section we formulate a new description in order to fill this gap.
We address the channel model at symbol level, the multiple access code, and the decoding function.
We propose a framework for the $T$-\ac{MPR} channel, including all these aspects, which, on one hand, represents a simple way to realize Assumptions 1-4 above and, on the other hand, provides a connection between random access method and encoding functions and sum rates.

\subsubsection{Channel model} The channel model we employ is the noiseless adder channel \cite{Dyachkov1981:Coding}.
The channel has $K$ inputs $X_1$, $\dots$, $X_{K}$ and one output $Y$, where $X_i \in \{0,1\} \subset \mathbb{R}$ for all $i \in \{1,\dots,K\}$, $Y \in \{0,\dots,K\} \subset \mathbb{R}$, and
\begin{align}
    Y = \sum_{i=1}^{K} X_i
\end{align}
the sum being over the reals.
Due to the binary input alphabets we refer to this channel as the \ac{BAC}.

\subsubsection{Codebook construction and encoding} To define the encoding function of each user we resort on the codebook construction proposed by \ac{BPR} in \cite{Bar-David1993:Forward}, where a coding technique to achieve zero error probability over the ``$K$-choose-$T$'' \ac{BAC} is proposed.
In this channel model there are $K$ users, of which exactly $T$ are contending.
The $T$ binary codewords are transmitted over a \ac{BAC}; the decoder aims at decoding the $T$ messages, with $T$ known a priori, and at associating each message to one of the $K$ users.

The \ac{BPR} construction starts from the parity-check matrix of a $T$-error-correcting binary \ac{BCH} code of length $M=2^m-1$, in the form
\begin{align}\label{eq:H}
    \mathbf{H} = \left[
    \begin{array}{ccccc}
        1 & \alpha & \alpha^2 & \dots & \alpha^{2^m-2} \\
        1 & \alpha^3 & \alpha^6 & \dots & \alpha^{3(2^m-2)} \\
        \vdots & \vdots & \vdots & \vdots & \vdots \\
        1 & \alpha^{2T-1} & \alpha^{2(2T-1)} & \dots & \alpha^{(2T-1)(2^m-2)}
    \end{array}
    \right]
\end{align}
with $mT$ binary rows (as the binary image of $\alpha$, a primitive element of $\mathrm{GF}(2^m)$, has length $m$) and $M$ columns.
The columns are partitioned into $K$ disjoint subsets, each of size $\frac{M}{K}$ (assuming $K$ divides $M$) and each subset is employed as the codebook of the corresponding user.\footnote{The assumption that $K$ divides $N_{\mathrm{s}}$ is made throughout the paper. In case $M$ is not a multiple of $K$, one can simply re-define $M$ as the largest multiple of $K$ that is less than $2^m-1$, and assign $M/K$ columns of $\mathbf{H}$ to each user.}
All users have the same information rate and the same codeword length $mT$.
Note that no codeword can be shared by any two users, that all codebooks are nonlinear, and that no codebook contains the all-$0$ codeword.
Note also that, although we stick to \ac{BCH} codes as in \cite{Bar-David1993:Forward}, the codebook construction may in principle rely on other error correcting code families.

The \ac{BPR} construction achieves zero error probability over the $K$-choose-$T$ \ac{BAC}.
In fact, since by construction any $2T$ columns of the parity-check matrix of a $T$-error-correcting linear block code are linearly independent, any two different $T$-tuples of codewords cannot have the same sum.
Hence, the set of $T$ messages can be decoded with no error at the receiver and its message can be uniquely associated with a user.

Next we address the proposed \ac{IRSA} encoder.
We refer to the described coding scheme as ``\ac{IRSA}-\ac{BPR}''.
\ac{BPR} $(m,T,K)$ codebooks are assigned to the $K$ users, as described above.
At the beginning of a frame, each user has an idle message $W=0$ (meaning nothing to communicate) with probability $1-\pi$ and one message of information $W \in \{1,\dots,\frac{M}{K}\}$ with probability $\pi$, independently from user to user.
Users with a message $W \neq 0$ are contending ones and the message of each contending user is uniformly distributed in $\{1,\dots,\frac{M}{K}\}$.
An idle user generates an idle sequence $X^N(0)$ of $N = (1+mT) N_{\mathrm{s}}$ symbols, all equal to $0 \in \mathbb{R}$.
On the other hand, a contending user generates a binary codeword $X^N(W)$ of length $N$ symbols, in the following way.
The message $W$ is encoded into a codeword of length $mT$ belonging to the user's \ac{BPR} codebook and one extra-symbol equal to $1 \in \mathbb{R}$ is appended to the \ac{BPR} codeword.
Then, the length-$N$ codeword is constructed as the concatenation of $d$ ``\ac{BPR} blocks'', each containing the \ac{BPR} codeword with the extra symbol $1$, and of $N_{\mathrm{s}}-d$ ``idle blocks'', each of length $1+mT$ and containing the symbol $0 \in \mathbb{R}$ only.
Importantly, to reflect the grant-free and uncoordinated nature of the access scheme, the number of \ac{BPR} blocks in the codeword, $d$, and the positions of these $d$ blocks are functions solely of the random message $W$.
A way to generate the repetition degree $d$ and the positions of the $d$ \ac{BPR} blocks out of $W$ is discussed in the Appendix.

\begin{remark}
Each of the $N_{\mathrm{s}}$ blocks composing a codeword $X^N(W)$ corresponds to one slot in the description provided in Section~\ref{sec:irsa_mpr}.
As such, the codeword of a contending user corresponds to the entire frame.
This is consistent with the message $W$ being transmitted using the channel $N$ times.
The information rate of a contending user is $R=r/N_{\mathrm{s}}$ where
\begin{align}\label{eq:rate_bpr}
    r = \frac{1}{1+ m T} \log_2 \frac{M}{K}
\end{align}
is the \ac{BPR} encoder rate with the extra symbol $1$.
The average sum rate, $R_{\mathrm{sum}}$, is given by $R_{\mathrm{sum}} = \mathbb{E}[K_{\mathrm{c}}]\, R$.
We note also that the encoder may be seen as the concatenation of two component encoders: a nonlinear \ac{BPR} encoder and an \ac{IRSA} ``random access'' encoder.
\end{remark}

\begin{remark}
The IRSA-BPR construction provides different codebooks to different users, thus allowing the receiver to identify the users that are sending messages.
A slightly simplified IRSA coding scheme may be also defined for the unsourced setting, where each user can use the entire codebook without any prior assignment, and the receiver is only able to produce an unordered list of messages.  The information rate for this unsourced construction can be easily derived, leading to an encoder rate expression as in (1), but without the term $\log_2K$, which appears as the cost of identifying the users that send messages. However, we should note that in the unsourced model the function of message association to users is left to upper layer protocols that spend at least $\log_2K$ bits  in the overhead of each message to accomplish this task.
\end{remark}

\subsubsection{Iterative Decoding} The vector $Y^N$ at the output of the \ac{BAC} is a sequence of $N_{\mathrm{s}}$ blocks, each of length $1+mT$ symbols.
We say a contending user is active in any such block when a \ac{BPR} block is present in that position in the user's codeword.
Decoding is performed iteratively; in every iteration all non-zero blocks of $Y^N$ are processed and, for each of them, \ac{BPR} decoding is run if the number of active users, known owing to the extra-symbol $1$, does not exceed $T$.
Each message $W$ decoded from a block can be associated with no error to a specific user; $W$ is re-encoded and the obtained length-$N$ codeword is subtracted from $Y^N$.
No action is taken in a block if the number of active users in it exceeds $T$.
Decoding terminates when $Y^N = 0^N$ (success) or when $Y^N \neq 0^N$ but no further subtraction can be performed (premature stop).

Denoting the iterative decoding function by $\phi(\cdot)$, we have
\begin{align}
\phi(Y^N) = \{ (W_{1},\dots,W_{q}), (i_1,\dots,i_q), \mathsf{F} \}
\end{align}
where $(W_1,\dots,W_q)$ is an ordered message list, $(i_1,\dots,i_q)$ is the ordered list of the corresponding user's indexes, $\mathsf{F}=0$ flags a success ($q=K_{\mathrm{c}}$, complete message list), and $\mathsf{F}=1$ flags a premature stop ($q<K_{\mathrm{c}}$, incomplete message list).

For any random access coding scheme we distinguish between the decoding error probability (probability that the obtained message list is different from the actual one and/or some incorrect user-message association is done) and the per-user error probability (probability that the message of a user is not in the decoded message list and/or that the user is not associated with its message).

\section{Results}\label{sec:results}

We provide both achievability and converse results on the expected sum rate.
The achievability one applies to the \ac{IRSA}-\ac{BPR} scheme, while the converse has a more general validity.
We always assume that, in each frame, each user is contending with probability $\pi$, independently of the other users.

\begin{theorem}\label{th:achievability}
Consider an \ac{IRSA}-\ac{BPR} scheme over the \ac{BAC}, with distribution $\Lambda$ and parameters $T$ and $m = \frac{1}{\epsilon} \log_2 K$, $0<\epsilon<1$.
If
\begin{align}\label{eq:achievability}
  R_{\mathrm{sum}} \leq \frac{1-\epsilon}{T} \, \mathsf{G}^{*}(\Lambda,T)
\end{align}
then the per user error probability under iterative decoding tends to zero in the limit as $N_{\mathrm{s}} \rightarrow \infty$ and constant $a=K/N_{\mathrm{s}}$, where $\mathsf{G}^{*}(\Lambda,T)$ is the load threshold of the \ac{IRSA} distribution $\Lambda$ over the $T$-\ac{MPR} channel model, defined in \eqref{eq:MPR_threshold}.
\end{theorem}
\begin{IEEEproof}
In the considered setting, Assumptions 1-4 reviewed in Section~\ref{sec:irsa_mpr} are all satisfied.
In fact, due to the addition of the extra symbol to each \ac{BPR} codeword, the receiver has constantly perfect knowledge of the number of \ac{BPR} blocks in each slot (Assumption~1) and takes no action in any slot with more than $T$ such blocks (Assumption~2).
Moreover, up to $T$ \ac{BPR} blocks in the same slot can always be decoded with zero error probability (Assumption~3).
Finally, over the \ac{BAC} interference subtraction is perfect since the channel is noiseless and no channel estimation is required. As such, iterative decoding for \ac{IRSA}-\ac{BPR} over the \ac{BAC} has a one-to-one correspondence with \ac{SIC} over the \ac{MPR} channel (Section~\ref{sec:irsa_mpr}) and the per-user error probability coincides with the packet loss probability mentioned in Remark~\ref{remark:plp}.
Next, the average sum rate is
\begingroup
\allowdisplaybreaks
\begin{align}
    R_{\mathrm{sum}} &= \mathbb{E}[K_{\mathrm{c}}] \frac{r}{N_{\mathrm{s}}} = \frac{\mathsf{G}}{1+ \frac{T}{\epsilon} \log_2 K} \log_2 \frac{K^{\frac{1}{\epsilon}} -1}{K} \label{eq:rsum_exact}
\end{align}
\endgroup
where \eqref{eq:rsum_exact} follows from \eqref{eq:rate_bpr} and $m = \frac{1}{\epsilon} \log_2 K$.
As $N_{\mathrm{s}} \rightarrow \infty$ for constant $a$, the per-user error probability vanishes for all $\mathsf{G} \leq \mathsf{G}^{*}(\Lambda,T)$ under \ac{SIC} decoding.
Since
\begin{align}\label{eq:limit_K_eps}
\frac{\log_2 \frac{K^{\frac{1}{\epsilon}} -1}{K} }{1+ \frac{T}{\epsilon} \log_2 K} \rightarrow \frac{1-\epsilon}{T}
\end{align}
as $K = a N_{\mathrm{s}} \rightarrow \infty$, we obtain the statement.
\end{IEEEproof}

\begin{remark}
We note that the slot size (number of channel uses per slot), and then the ``packet'' size, is $\Theta(\log_2 K)$.
We may interpret this as a consequence of the decoder task to perform association between decoded messages and contending users.
\end{remark}

We now address a converse result. We start by proving the following useful lemma.

\begin{lemma}\label{lemma:loeve}
Let $C$ be any positive random variable.
For constant $a = K/N_{\mathrm{s}}$, if
$
\lim_{\substack{N_{\mathrm{s}} \rightarrow \infty}} \Pr(K_{\mathrm{c}} > C) = 0
$
then
\begin{align}\label{eq:G_ineq_conv}
    \mathsf{G} \leq \lim_{N_{\mathrm{s}} \rightarrow \infty} \frac{\mathbb{E}[C]}{N_{\mathrm{s}}} \, .
\end{align}
\end{lemma}
\begin{IEEEproof}
The statement can be proved by applying the simple lower bound on the probability of the upper tail provided in Lo{\'e}ve \cite{loeve1977:Probability}, stating that for a real random variable $X$ with support $\mathcal{X}$ and a nonnegative and nondecreasing function $g(\cdot) : \mathbb{R} \mapsto \mathbb{R}$ we have
\begin{align}\label{eq:loeve}
    \Pr(X > a) \geq \frac{\mathbb{E}[g(X)]-g(a)}{\sup g(x)}
\end{align}
where $\sup g(x) = \sup \{g(x) : x \in \mathcal{X} \}$.
Applying \eqref{eq:loeve} with $X = K_{\mathrm{c}} - C$, $g(x)=\max\{0,x\}$, and $a=0$, we obtain
\begingroup
\allowdisplaybreaks
\begin{align}
\Pr(K_{\mathrm{c}} > C) &\geq \frac{\mathbb{E}[\max\{0,K_{\mathrm{c}}-C\}]}{\sup g(k_{\mathrm{c}}-c)} \\
&\geq \frac{\max\{0,\mathbb{E}[K_{\mathrm{c}}]-\mathbb{E}[C]\}}{\sup g(k_{\mathrm{c}}-c)} \label{eq:jensen}
\end{align}
\endgroup
where \eqref{eq:jensen} follows from Jensen's inequality since $g(x)$ is convex.
By hypothesis, for any $\epsilon>0$ there exists $N_{\mathrm{s}}(\epsilon)$ such that $\Pr(K_{\mathrm{c}} > C) \leq \epsilon$ for all $N_{\mathrm{s}} > N_{\mathrm{s}}(\epsilon)$.
It follows that, for all $N_{\mathrm{s}} > N_{\mathrm{s}}(\epsilon)$, we also have
\begin{align}
\frac{\max\{0,\mathbb{E}[K_{\mathrm{c}}]-\mathbb{E}[C]\}}{\sup g(k_{\mathrm{c}}-c)} \leq \epsilon
\end{align}
and then
\begingroup
\allowdisplaybreaks
\begin{align}
\max\{0,\mathbb{E}[K_{\mathrm{c}}]-\mathbb{E}[C]\} &\leq \epsilon \, \sup g(k_{\mathrm{c}}-c) \\
&\leq \epsilon\, K \label{eq:eps_K}
\end{align}
\endgroup
where \eqref{eq:eps_K} is due to $k_{\mathrm{c}} \leq K$ and $c > 0$. Inequality \eqref{eq:eps_K} implies
$\pi K -\mathbb{E}[C] \leq \epsilon \, K$. Dividing both sides by $N_{\mathrm{s}}$ we obtain
$\mathsf{G} \leq \frac{\mathbb{E}[C]}{N_{\mathrm{s}}} + \epsilon \, a$.
Taking the limit $N_{\mathrm{s}} \rightarrow \infty$, $\epsilon$ becomes arbitrarily small; noting that $a$ is constant (so $\epsilon \, a$ vanishes) yields the statement.
\end{IEEEproof}

\medskip
We can use Lemma~\ref{lemma:loeve} to obtain a converse result valid for any framed and slotted access scheme over the \ac{BAC} (even cooperative or grant-based) with \ac{BPR} codebooks assigned to users and the slot size matching the \ac{BPR} codeword size.
We keep appending the extra symbol $1$ to the \ac{BPR} codeword, although not strictly necessary for the proof.

\begin{theorem}\label{th:converse}
Let \ac{BPR} codebooks be assigned to the $K$ users, with $m = \frac{1}{\epsilon} \log_2 K$, and the channel model be the \ac{BAC}.
Assume each contending user generates one \ac{BPR} codeword of length $n=1+mT$ out of its message and transmits replicas of it in any number of slots, where the slot size is also $n$, according to any access protocol.
Length-$n$ all-zero words are transmitted by the user in the other slots.
Moreover, let $l_i$ be the number of \ac{BPR} codewords added in slot $i$ and
\begin{align}\label{eq:ci}
C_i = \left\{
\begin{array}{ll}
   l_i  & \text{if } l_i \leq T \\
   T  & \text{if } l_i > T .
\end{array}
\right.
\end{align}
Let $P_e$ be the error probability under optimum decoding.
If $P_e \rightarrow 0$ as $N_{\mathrm{s}} \rightarrow \infty$ for constant $a = K / N_{\mathrm{s}}$, then
\begin{align}\label{eq:converse}
R_{\mathrm{sum}} \leq \frac{1-\epsilon}{T} \, \mathbb{G}
\end{align}
where $\mathbb{G}$ is the sup of the set of all $\mathsf{G}$ fulfilling \eqref{eq:G_ineq_conv} and $C = \sum_{i=0}^{N_{\mathrm{s}}-1} C_i$ with $C_i$ defined in \eqref{eq:ci}.
\end{theorem}
\begin{IEEEproof}
An optimum, idealized decoder is the one aware of which users are active in each slot and formulating the decoding problem over the \ac{BAC} as the solution of a linear system of equations in which the unknowns are the $K_{\mathrm{c}}$ \ac{BPR} codewords.
The number of linearly independent equations available in slot $i$ is $C_i$ defined in \eqref{eq:ci}, so the total number of such equations is $C = \sum_{i=0}^{N_{\mathrm{s}}-1} C_i$.
If the error probability of this optimum decoder vanishes as $N_{\mathrm{s}} \rightarrow \infty$ and constant $a = K/N_{\mathrm{s}}$, then the probability that the number of unknown exceeds the number of equations, $\Pr(K_{\mathrm{c}} > C)$, also vanishes.
Lemma~\ref{lemma:loeve}, with \eqref{eq:rsum_exact} and \eqref{eq:limit_K_eps}, then leads to the statement.
\end{IEEEproof}

\begin{remark}
Comparing \eqref{eq:achievability} and \eqref{eq:converse}, we can note how their structure is analogous.
As mentioned above, however, while Theorem~\ref{th:achievability} applies to \ac{IRSA}-\ac{BPR} under \ac{SIC} decoding, Theorem~\ref{th:converse} has a broader scope.
A more subtle difference is that Theorem~\ref{th:achievability} considers the per-user error probability under iterative decoding while Theorem~\ref{th:converse} the decoding error probability under an optimum decoder.
\end{remark}

\begin{example}
Consider the following three slotted and framed access schemes: (i) A coordinated scheme where each contending user fills one slot with a \ac{BPR} block and the other slots with all-zero blocks; the slot index is scheduled by the receiver to maximize the sum rate. (ii) An uncoordinated scheme where all contending users employ \ac{IRSA}-\ac{BPR}. (iii) A mixed scheme in which $\nu K$ users have high priority and are served in a coordinated fashion, as in scheme (i), while the other $(1-\nu) K$ users employ \ac{IRSA}-\ac{BPR}; $\gamma N_{\mathrm{s}}$ slots are reserved to high-priority users while the other $(1-\gamma) N_{\mathrm{s}}$ slots are available to low-priority ones, but the high-priority slots still empty after scheduled transmissions can be used by low-priority users as per decoder feedback.
Specifying \eqref{eq:G_ineq_conv} for the three schemes, we obtain $\lim_{N_{\mathrm{s}} \rightarrow \infty} \mathbb{E}[C]/N_{\mathrm{s}} = \min \{\mathsf{G},T\}$, hence $\mathbb{G}=T$, for scheme (i). Proof sketch: $K_{\mathrm{c}}/N_{\mathrm{s}} \rightarrow \mathsf{G}$ in probability as $N_{\mathrm{s}} \rightarrow \infty$; $f(x)=\min\{x,T\}$ is continuous so $C/N_{\mathrm{s}}=f(K_{\mathrm{c}}/N_{\mathrm{s}}) \rightarrow f(\mathsf{G})$ in probability; $C/N_{\mathrm{s}}$ is bounded then we can conclude that $\mathbb{E}[C/N_{\mathrm{s}}] \rightarrow f(\mathsf{G})$. For scheme (iii) we obtain (proof omitted)
\begin{align*}
    \mathsf{G} &\leq T - (1 - \beta)
    \exp\! \bigg(\!\! - \frac{\mathsf{G} (1-\nu)}{\eta (1-\beta)} \bigg) \!\sum_{t=0}^{T-1} \frac{T-t}{t!} \bigg( \frac{\mathsf{G} (1-\nu)}{\eta (1 - \beta)} \bigg)^{\!t}
\end{align*}
where $\beta = \min \{ \frac{\nu \mathsf{G}}{T}, \gamma \}$; this expression holds also for scheme (ii) with $\nu=\beta=0$. In both schemes (ii) and (iii), $\mathbb{G}$ must be computed numerically.
The obtained values of \eqref{eq:converse} are plotted in Fig.~\ref{fig:converse} versus $T$ for $\epsilon=0.99$.
Two values of $\eta$ are considered, namely, $\eta=1/2$ and $\eta=1/3$; concerning the mixed scheme, the plot is for $\nu=\gamma=0.2$.
The region above each curve is not achievable by the corresponding scheme.
The value of $\mathbb{G}$ for scheme (ii) also appears in \cite{stefanovic2018:multipacket}, but it was proved only under \ac{SIC}, while it is shown here that it holds for any decoder.
\end{example}

\begin{figure}
    \begin{center}
    \includegraphics[width=0.98\columnwidth]{./Figures/Figure2.pdf}
    \end{center}
    \vspace{-5mm}
    \caption{Converse \eqref{eq:converse} versus $T$ for $\epsilon=0.99$. Solid: coordinated scheme; dashed: \ac{BPR}-\ac{IRSA} scheme; dot-dashed: mixed scheme. \ac{IRSA} efficiency $\eta=1/2$ ($\mathsf{o}$) and $\eta=1/3$ ($\mathsf{x}$).}
    \label{fig:converse}
\end{figure}

\begin{comment}
\begin{figure}[t]
        \begin{center}
        \resizebox{0.5\textwidth}{!}{%
    		\input{ISIT_2022/Figures/converse}
    	}%
    	\end{center}
    \caption{Converse \eqref{eq:converse} versus $T$ for $\epsilon=0.99$. Solid: coordinated scheme; dashed: \ac{BPR}-\ac{IRSA} scheme; dot-dashed: mixed scheme. \ac{IRSA} efficiency $\eta=1/2$ ($\mathsf{o}$) and $\eta=1/3$ ($\mathsf{x}$).}
    \label{fig:converse}
\end{figure}
\end{comment}

\section{Conclusion}\label{sec:conclusions}

In this paper a new framework for \ac{IRSA} has been proposed, in which rigorous results on the average sum rate can be established.
A number of possible research directions, we believe, arise from this work.
Among them, the development of converse results capturing the \ac{IRSA} distribution (and not just its efficiency), extension to \ac{CSA} schemes, extensions to channel models incorporating noise and fading and to different settings, e.g., to asynchronous schemes or to a finite number of users generating an increasing traffic as the frame size grows.

\section*{Appendix}
This appendix discusses generation of the repetition rate $d$ and of a $d$-tuple of different integers in $\{0,\dots,N_{\mathrm{s}}-1\}$ as a function of a message $W$ uniformly distributed in $\{1,\dots,\frac{M}{K}\}$.
For $\frac{M}{K} \Lambda_d$ an integer for all $d\in\{2,\dots,\bar{d}\}$, an instance of the random variable $D$ with \ac{PGF} $\Lambda(x)$, mentioned in Section~\ref{sec:irsa_mpr}, can be generated as
\begin{align}\label{eq:d_from_W}
    d = \min \bigg\{ i \in \{2,\dots,\bar{d}\} : \frac{M}{K} \sum_{h=2}^{i} \Lambda_h \geq W \bigg\} \, .
\end{align}
Once the repetition degree $d$ has been obtained from \eqref{eq:d_from_W}, the positions of the $d$ \ac{BPR} blocks can be generated by resorting on the combinadic representation of an integer (e.g., \cite{Knuth2011:Art}).
Accordingly, for any integer $0 \leq u < \binom{N_{\mathrm{s}}}{d}$, there exists a unique $d$-tuple of integers $(u_1,\dots,u_d)$, with $0 \leq u_1 < \dots < u_{d}$, such that
$u = \binom{u_1}{1} + \binom{u_2}{2} + \dots + \binom{u_d}{d}$. %
In particular, for $h = d, \dots, 1$, $u_h$ is computed as
\begin{align}\label{eq:u_h}
    u_h = \max \bigg\{ i \in \mathbb{N} : \binom{i}{h} \leq u - \sum_{j=h+1}^{d} \binom{u_j}{j} \bigg\} \, .
\end{align}
The positions of the $d$ blocks can be computed by first letting $v = W-\frac{M}{K}\sum_{i=2}^{d-1} \Lambda_i-1$ (note that $v \in \{0,\dots,\frac{M}{K}\Lambda_d-1\}$), then by defining $u = v \mod \binom{N_{\mathrm{s}}}{d}$, and finally by returning $(u_1, \dots, u_d)$ through \eqref{eq:u_h}.
The desired uniform probability distribution over the $\binom{N_{\mathrm{s}}}{d}$ combinations is obtained when $\binom{N_{\mathrm{s}}}{d}$ divides $\frac{M}{K} \Lambda_d$ for all $d\in\{2,\dots,\bar{d}\}$.

\section*{Acknowledgement}
This work has been carried out in the framework of the CNIT National Laboratory WiLab and the WiLab-Huawei Joint Innovation Center.
The authors would like to thank Dr. Lei Wang for useful discussions.