\section{Introduction}

Experimental techniques in single-molecule optical trapping and biophysics allow to extract  real-time information  of the state of a nanosystem with exquisite precision~\cite{bustamante2005nonequilibrium,ashkin1986observation,haroche1991trapping}.  Such information is commonly used to infer both thermodynamical and dynamical properties through data-analysis techniques.
Alongside,  as inspired by Maxwell's demon thought experiment, information acquired from a nanosystem  can be delivered into  work by executing feedback-control protocols~\cite{bechhoefer2005feedback,toyabe2010experimental,campisi2017feedback,ciliberto2020autonomous,parrondo1996criticism}.
In parallel to experimental progress, the development of  stochastic thermodynamics (ST) over the last two decades provides a robust  theoretical framework to  describe accurately information-to-work transduction that takes into account  nanoscale   fluctuations~\cite{Sekimoto2010,jarzynski2011equalities,Seifert2012,van2015ensemble,peliti2021stochastic}. Combining  stochastic thermodynamics and feedback-cooling techniques has attracted attention towards refrigerating capabilities of small systems under nonequilibrium conditions~\cite{gieseler2012subkelvin,gieseler2014dynamic}.

An important step to optimize the design of microscopic refrigerators  is to bridge the gap between theoretical proposals and experiments
through the powerful method of all-atoms simulations. Nonequilibrium Molecular Dynamics (MD) studies  provide a suitable platform for the study of heat transfer and fluctuations at the nanoscale~\cite{rajabpour2019thermal,todd2017nonequilibrium, roodbari2022interfacial, muller1997simple, allen2017computer}.  However, little is known yet about the design and performance of information demons at the atomic scale. In particular, are there generic principles that constraint  the forces needed to ensure a prescribed value for the heat transfer between two thermal  baths interacting through nanoscopic objects?
Is it possible to accurately control the net heat transfer between nanoparticles and their respective fluctuations by {\em only} applying non-conservative forces, i.e. forces that do not derive from a potential?

Among the broad class of non-conservative forces, nonreciprocal interactions (i.e. forces that violate Newton's third law ``actio=reactio")
have recently emerged as a topic of lively interest in statistical physics \cite{fruchart2021non,saha2020scalar,you2020nonreciprocity,lavergne2019group,loos2022long}, revealing  nontrivial physical consequences for the dynamical, mechanical and thermodynamic properties of many-body systems. For example, they introduce `odd elasticity' in solids and soft crystals~\cite{braverman2021topological,poncet2022soft} or lead to the formation of travelling waves in binary fluid mixtures~\cite{fruchart2021non,saha2020scalar,you2020nonreciprocity}. In stochastic thermodynamics, recent research has revealed the potential  of nonreciprocal forces in the design of  artificial nano machines  with efficient energetic performance~\cite{Loos2020b}. Inspired by these recent findings, herein we design atomistic MD simulations of trapped nanoparticles immersed in thermal baths at different temperatures that interact through non-conservative, linear forces and that are nonreciprocal. A similar setup was realized experimentally very recently using optical fields \cite{rieser2022tunable}.
Here, we use nonreciprocal interactions to construct a {\em nano refrigerator} which realizes a steady, net heat flow from the cold to the hot bath that does not fulfill Fourier's law for thermal conduction yet  is in  agreement with recent theoretical predictions from ST.  A key advantage of our nano refrigerator design relies on its simplicity as it only requires the usage of nonreciprocal forces acting on each of the nanoparticles. This represents a simplification with respect to previous approaches where heat flows from hot to cold could be achieved using velocity-dependent feedback~\cite{Maxwell1986} or memory registers~\cite{Mandal2012,mandal2013maxwell} as in Maxwell's demons,  or using nonlinear forces in athermal environments~\cite{kanazawa2015minimal}.

Our work establishes theoretical design  principles that ensure a prescribed net heat flux in our MD simulations that could be exported to realistic experimental scenarios with trapped nanoparticles~\cite{bechhoefer2005feedback,proesmans2016brownian,ricci2017optically,midtvedt2021fast}.
We also test  fundamental principles governing the fluctuations of work and the coefficient of performance (COP), some of which follow from recently-established thermodynamic uncertainty relations tested here with realistic atomistic simulations of  nanoparticles~\cite{barato2015thermodynamic,horowitz2020thermodynamic}. These results push forward the synergistic combination of ST  and MD beyond the application of fluctuation theorems in e.g. estimating free energies~~\cite{dellago2014computing,park2003free}.  In particular, our simulations made with  parameters for realistic materials
are a first step towards the engineered design of nanoparticle-based refrigerators powered by thermal fluctuations.

\section{Model and Setup}

\begin{figure}
\centering
\includegraphics[width=0.4\textwidth]{Figures/DemonNRSketch.pdf}
\caption{
	 Sketch of the molecular dynamics (MD) simulation setup. Two Copper nanoparticles (black and yellow spheres) each consisting of 186 atoms, are immersed in two Argon baths at different temperatures $T_\mathrm{C}=100 \mathrm{K}$ and $T_\mathrm{H}=120\mathrm{K}$. The nanoparticles are trapped  with two static three-dimensional harmonic potentials of different strength (see black line for an illustration). A demon-like controller (green) exerts additional nonreciprocal forces to the two nanoparticles as follows. The demon measures the nanoparticles' positions $X_\mathrm{C},X_\mathrm{H}$ and exerts the forces $\kappa_\mathrm{C}X_\mathrm{H}$ to the nanoparticle in the cold bath and $\kappa_\mathrm{H}X_\mathrm{C}$ to the nanoparticle in the hot bath, where in general $\kappa_\mathrm{H}\neq \kappa_\mathrm{C}$ rendering the demon forces nonreciprocal.
	}\label{fig:MD}
\end{figure}

\subsection{Nanodemon setup and MD simulations}

Constructing an atomistic MD simulation that allows us to violate Newton's third law and furthermore realize a nonreciprocal nano refrigerator, requires a highly unconventional setup in nonequilibrium MD studies. Specifically, we simulate two independent Argon baths that are kept thermostatted at different temperatures ($T_\mathrm{C}=100$K and $ T_\mathrm{H}=120$K). The cold and hot bath are in turn separated by a hard wall made of immobile Copper particles (see Figure 1). Within each bath, we immerse a copper-based spherical nanoparticle of radius $1.4 \mathrm{nm}$.
External nonreciprocal forces (sketched as a green demon in Fig.~\ref{fig:MD}) are applied on the center of each nanoparticle through two  forces, $\kappa_\mathrm{C}X_\mathrm{H}$ and $\kappa_\mathrm{H}X_\mathrm{C}$ on the nanoparticles immersed in the cold and hot baths respectively. Here $X_\mathrm{H}$ and $X_\mathrm{C}$ denote respectively the center-of-mass position of the particle in the hot and the particle in the cold bath, respectively.  Such a setup could, in principle, be realized with the help of an external control scheme (e.g., using optical feedback) similar to~\cite{Lavergne2019,Khadka2018}. As we will see shortly, when $\kappa_\mathrm{C}\neq \kappa_\mathrm{H}$, a \textit{nonreciprocal} coupling is introduced which can lead for specific values $\kappa_\mathrm{C}/\kappa_\mathrm{H}$ to a heat flow from the cold to the hot bath. We have further constrained the particle positions by introducing harmonic potentials (with stiffness $\kappa = 1.16 k_\mathrm{B}T_\mathrm{C}/\AA^2$). The traps prevent the particles from hitting the walls, but are, in principle, not needed to construct the nano refrigerator (i.e., we could also set $\kappa=0$), as is evident from our analytical results introduced in the following.

\subsection{Heat transfer}

Our MD setup gives us direct access to thermodynamic quantities allowing for quantitative measurements of the heat transfer between the nanoparticles and their respective solvent baths. Specifically, we determined the
total amount of energy change by extracting both the potential and kinetic energy
of all bath molecules as a function of time which gives the total heat transferred by the copper nanoparticles to both the cold and hot baths, denoted by $\mathrm{d}Q_\mathrm{C}$ and $\mathrm{d}Q_\mathrm{H}$ respectively.
Integrating over the course of the MD simulation yields the $Q_\mathrm{C}$ and $Q_\mathrm{H}$, which directly encodes the stochastic heat dissipated by the nanoparticle into the cold and hot bath respectively. Note that we use the sign convention that
${Q}>0$ when net energy is dissipated from the nanoparticle to the bath and ${Q}<0$ when it is absorbed by the nanoparticle from the bath.
We estimate the heat dissipation rate $\dot{Q}_\mathrm{C}$ and $\dot{Q}_\mathrm{H}$ from the slope of a  linear regression on the cumulative $Q_\mathrm{C}$ and $Q_\mathrm{H}$ over time.

With this protocol in hand, we begin by demonstrating how tuning the relative strength of $\kappa_\mathrm{C}$ and $\kappa_\mathrm{H}$ provides a microscopic mechanism to alter the direction of heat flow.  Figure~\ref{fig:heat} (a) illustrates $Q_\mathrm{C}$ and $Q_\mathrm{H}$ a situation where the effective coupling force on the
hot particle is reduced as
$\kappa_\mathrm{H}\ll \kappa_\mathrm{C}$.
In this case, we observe the canonical situation, where the nanoparticle-duet behaves as a {\em heater}, i.e. heat flows from the hot to the cold bath. On the other hand, by introducing
an effectively enhanced coupling force experienced by the nanoparticle in the hot bath, there is a striking effect where the direction of the heat flow changes as seen in Figure~\ref{fig:heat} (b)  - heat is now pumped from the cold to the hot bath creating a molecular-scale refrigerator.
The preceding results from the MD simulations provide a powerful proof-of-concept on how an nonreciprocal forces applied on two nanoparticles embedded in  a solvent bath, can in principle be used to change both the rate and direction of heat flow.

\begin{figure}[ht]
\centering
\includegraphics[width=0.49\textwidth]{Figures/Figure2_cartoon_v02.pdf}
	\caption{Heat flows $Q$ measured in the MD simulations of the setup sketched in Fig.~\ref{fig:MD}, for two different values of $\kappa_\mathrm{H}$.
	 (a) {\em Heater Type-II}: for $\kappa_\mathrm{H} =0.17\kappa_\mathrm{C}$,  heat flows from the hot bath to the nanoparticle in the right container and from the nanoparticle to the cold bath in the left container.
	(b)  {\em Refrigerator/Heat-pump}:  for $\kappa_\mathrm{H} =4.80\kappa_\mathrm{C}$,  heat flows from the cold bath to the nanoparticle in the left container and from the nanoparticle to the hot bath in the right container. Here,  heat flows in the reverse direction to the temperature gradient, i.e., it is  extracted from the cold bath and released into the hot bath through the two-nanoparticle system.
	\textit{Upper panels} of (a,b) show sketches of the flows of heat and work, respectively.
	\textit{Middle panels} show the cumulative heat as function of time from the MD simulations, dissipated into the cold bath ($Q_\mathrm{C}$, blue line) and from the hot bath ($Q_\mathrm{H}$, red line).  We use the convention $Q>0$ when heat flows from the nanoparticle into the bath and  $Q<0$  when heat flows from the bath to the nanoparticle.
	Black lines show linear fits used to obtain the heat rates. \textit{Lower panels} display the stochastic work exerted on the nanoparticle in the cold bath  ($W_{\rm C}$, blue line),  and on the nanoparticle in the hot bath  ($W_{\rm H}$, red line), and the total work given by their sum  ($W=W_{\rm C}+W_{\rm H}$, green line). The work is obtained from the MD simulations by measuring the center of mass positions of the nanoparticles, and evaluating the work using Eq.~\eqref{eq:stochpower}.
 In both (a, b) we used $\kappa_\mathrm{C}=10\kappa =11.6 k_\mathrm{B}T_\mathrm{C}/\AA ^2 $, $T_\mathrm{C}= 100\mathrm{K}$ and $T_\mathrm{H}=120\mathrm{K}$.}
\label{fig:heat}
\end{figure}

Intuitively, by increasing $\kappa_\mathrm{H}$, the demon tricks the hot particle into \emph{``seeing''} that it was coupled to an even hotter particle, while the cold particle thinks it was coupled to an even colder particle, which results in a heat transfer from the cold to the hot bath. In the following, we will formulate and present a theory that rationalizes this intriguing phenomenon and makes a direct link between thermodynamic observables extracted from the MD simulations and stochastic thermodynamics.

\subsection{Stochastic model}\label{sec:LE}

We employ a mesoscopic stochastic model to describe the  nonequilibrium dynamics of the position and momentum fluctuations of the two-nanoparticle system in their thermal environments.
To this aim, we describe at a coarse-grained level the dynamics of the $x$-components of the positions and velocities of the center of mass of the nanoparticles, $X_\mathrm{C}$ and $X_\mathrm{H}$, by two coupled underdamped Langevin equations,
\begin{align}\label{eq:Le}
m \ddot{X}_\mathrm{C}+\gamma_\mathrm{C} \dot{X}_\mathrm{C} &=  -(\kappa_\mathrm{C} + {\kappa} ) X_\mathrm{C} + \kappa_\mathrm{C} X_\mathrm{H} + \xi_\mathrm{C}
\\ \label{eq:Le2}
m \ddot{X}_\mathrm{H}+\gamma_\mathrm{H} \dot{X}_\mathrm{H} &
=
\kappa_\mathrm{H}
X_\mathrm{C} -({\kappa}+\kappa_\mathrm{H}) X_\mathrm{H} + \xi_\mathrm{H}.
\end{align}
Here, $m$ is the mass of each nanoparticle, and the coefficients $\kappa_\mathrm{C}$, $\kappa_\mathrm{H}$, $\kappa$ have been defined before. Note that the dynamics is independent of the actual distance between both nanoparticles, which has therefore been excluded from the equations of motions \eqref{eq:Le} by a change of variables (see Methods section).
The stochastic forces $\xi_\mathrm{C}$, $\xi_\mathrm{H}$ are independent Gaussian white noises with zero mean $\langle\xi_\mathrm{C}(t)\rangle=\langle\xi_\mathrm{H}(t)\rangle=0$ modeling the thermal noise exerted by the Argon bath that surrounds  each nanoparticle.
Their autocorrelation functions are $\langle \xi_j(t)\xi_l(t') \rangle =2 k_\mathrm{B}T_\mu \gamma_\mu \,\delta_{jl}\delta(t-t')$, where $ j,l\in \{\mathrm{C,H}\}$ are indices denoting the hot or cold bath, $\delta_{jl}$ is Kroneckers' delta, and $k_\mathrm{B}$ is  Boltzmann's constant. Here and in the following, $\langle . \rangle$ denote averages over many realizations of the noise. The averages obtained from the MD simulations are extracted from {single} trajectories of $15$ns long, i.e. exceeding by three orders of magnitude the relaxation times of our system (see below). %
Furthermore, $\gamma_\mathrm{C}$ and $\gamma_\mathrm{H}$ are {\em effective}  coefficients of the friction  that each nanoparticle experiences in its respective environment.
Despite its simplicity, the model \eqref{eq:Le} allows us to infer dynamical and thermodynamic properties of our molecular dynamics simulations, as we describe below.

An important ingredient for the theory are the effective friction coefficients $\gamma_\mathrm{C}$ and $\gamma_\mathrm{H}$, which  emerge from the interaction of the nanoparticles with the baths' particles~\cite{stewart1962transport}. To estimate  $\gamma_\mathrm{C}$ and $\gamma_\mathrm{H}$, we run an equilibrium simulation in the absence of (nonreciprocal) interactions by taking $\kappa_\mathrm{C}=\kappa_\mathrm{H}=0$. In this limit, the position and velocity autocorrelation functions can be solved analytically (see Methods section). By fitting the autocorrelation functions obtained from the MD simulations to the analytical formulas, we extract the estimates $\gamma_\mathrm{C}\simeq \gamma_\mathrm{H}  \sim 3.5 \times 10^{-12}$kg/s. This value is consistent with an independent estimate obtained from a \emph{dragging experiment} simulation where a nanoparticle is pulled with a constant force along  the $x$ axis (see SI for further details), which yields an estimate of $\gamma_\mathrm{C} \approx \gamma_\mathrm{H} \approx 3\times 10^{-12}$kg/s. In the following and for further analyses, we use the estimate $\gamma = 3.5 \times 10^{-12}$kg/s for the friction of the two nanoparticles.

In the ensuing analysis, several timescales are relevant to characterize the dynamics of the particles at the mesoscopic level, following~\eqref{eq:Le} and~\eqref{eq:Le2}.
Firstly, the momentum relaxation time extracted from our data $m/\gamma \approx 5$ps is only one order of magnitude smaller than the relaxation time for the position in the trap $\gamma/{\kappa} \approx 25$ps. %
These results lend credence to the validity of the underdamped description used in our approach, as the  data acquisition timescale $\Delta t\sim 1$ps is below the
momentum relaxation timescale.

\subsection{Stochastic Energetics\label{sec:Energy}}

\begin{figure}[t]
\centering
\includegraphics[width=0.48\textwidth]{Figures/Figure3_heat-epr-work.pdf}\\
	\caption{Comparison between thermodynamic fluxes obtained from MD simulations (symbols) and theoretical predictions from stochastic thermodynamics (lines): 	 The average rate of  heat dissipation rate in each bath (a), the average entropy production rate (b), and the total power applied by the demon to the system (c), as a function of the coupling constant $\kappa_\mathrm{H}$.
   The heat rates  in (a) are measured by accumulating the amount of energy extracted by the thermostats and linearly fitting the results to obtain the slope, see Fig.~\ref{fig:heat}.
The blue shaded area illustrates the parameter regime at which the system works as a refrigerator (given by  the condition~\ref{cond1}).}\label{fig:heat2}
\end{figure}

We now develop and put to the test a theory for the energetics of the nanoparticle system setup in the light of stochastic thermodynamics~\cite{Sekimoto2010}---a framework that enables to describe the fluctuations of heat and work of systems described by e.g. Langevin equations, such as  Eq.~(\ref{eq:Le}).
Applying the framework of stochastic thermodynamics to the model given by  Eq.~(\ref{eq:Le}), and putting forward  mathematical techniques introduced in~\cite{kwon2005structure,bae2021inertial,kwon2011nonequilibrium},  we derive exact closed-form expressions for the expected value of the rate of heat dissipated into the cold and  hot baths' in the steady state, reading as,
\begin{eqnarray}
 \langle\dot{Q}_\mathrm{C}\rangle
 =-k_\mathrm{B}
\frac{ \displaystyle \kappa_\mathrm{C} (\kappa_\mathrm{C} T_\mathrm{H}-\kappa_\mathrm{H} T_\mathrm{C})}{\displaystyle(m/\gamma)  (\kappa_\mathrm{C} + \kappa_\mathrm{H})^2 /2+ \gamma  (2\kappa+\kappa_\mathrm{C}+\kappa_\mathrm{H})} ,
 \label{eq:Qc}
\\
 \langle\dot{Q}_\mathrm{H}\rangle
 =-k_\mathrm{B}
\frac{\displaystyle- \kappa_\mathrm{H} (\kappa_\mathrm{C} T_\mathrm{H}-\kappa_\mathrm{H} T_\mathrm{C})}{\displaystyle(m/\gamma)  (\kappa_\mathrm{C} + \kappa_\mathrm{H})^2 /2+ \gamma  (2\kappa+\kappa_\mathrm{C}+\kappa_\mathrm{H})} . \label{eq:Qh}
\end{eqnarray}
To ensure that the dynamics of the nanoparticles is stable, we further find the necessary condition $\kappa_\mathrm{H} + \kappa_\mathrm{C} > -\kappa$ (see Methods section).
Note that in the above formulas, the steady-state averages $\langle\dot{Q}_\mathrm{C}\rangle= \langle \mathrm{d}{Q}_\mathrm{C}/\mathrm{d}t\rangle$ and $\langle\dot{Q}_\mathrm{H}\rangle= \langle \mathrm{d}{Q}_\mathrm{H}/\mathrm{d}t\rangle$ are obtained using the definitions of stochastic heat used in stochastic thermodynamics  (see Methods section), which are written in terms of the nanoparticles' positions  and velocities and thus not necessarily equal to the ``direct" stochastic heat measured in the MD simulations from the energy fluctuations of the bath molecules.

Notably, Eqs.~(\ref{eq:Qc}-\ref{eq:Qh}) predict a  net heat transfer between the two baths that obeys Fourier's law $\langle\dot{Q}_\mathrm{C,H}\rangle \propto (T_\mathrm{H}-T_\mathrm{C})$
only when the forces exerted by the demon are reciprocal. Moreover,
our theory predicts that a net heat flow can be induced by three mechanisms: first, the existence of a temperature gradient, second, nonreciprocal coupling, or, third, a combination of both.
As we show below, such a net heat flow is a signature of entropy production, with the latter being also accessible from our theory. In particular,  Eqs.~(\ref{eq:Qc},\ref{eq:Qh}) allow for the prediction of a closed-form theoretical expression for the steady-state average rate of entropy production $
\langle \dot{S}_\mathrm{tot}\rangle =
\langle \dot{Q}_\mathrm{H}\rangle /T_\mathrm{H} +\langle \dot{Q}_\mathrm{C} \rangle/T_\mathrm{C}$~\cite{Seifert2012}, yielding
\begin{align}
\langle \dot{S}_\mathrm{tot} \rangle
& =
\frac{k_\mathrm{B}( \kappa_\mathrm{C}T_\mathrm{H}  - \kappa_\mathrm{H} T_\mathrm{C})^2/(T_\mathrm{C}T_\mathrm{H})}
{(m/\gamma) (\kappa_\mathrm{C} +\kappa_\mathrm{H})^2/2 + \gamma(\kappa_\mathrm{C}+\kappa_\mathrm{H}+2\kappa)}
.\label{eq:Stot}
\end{align}
Note that $\langle \dot{S}_\mathrm{tot} \rangle\geq 0$ for any parameter values, in agreement with the second law of stochastic thermodynamics, with equality only for the choice $\kappa_\mathrm{H}/ \kappa_\mathrm{C}= T_\mathrm{H}/T_\mathrm{C}$ for the nonreciprocal coupling constants. Furthermore, from Eqs.~(\ref{eq:Qc}-\ref{eq:Qh}) we can also extract the total power exerted by the demon on the nanoparticles, which is simply given by $\langle \dot{W}\rangle = \langle \dot{Q}_\mathrm{C}+\dot{Q}_\mathrm{H} \rangle$, see below.

To gain further insights, we compare  in Fig.~\ref{fig:heat2}a our theoretical predictions from stochastic thermodynamics (Eqs.~(\ref{eq:Qc}-\ref{eq:Qh}), lines) for the heat transfer evaluated directly over a collection of MD simulations (symbols) ran over a wide range of values for the demon coupling constants $\kappa_\mathrm{H}$ and $\kappa_\mathrm{C}$ spanning three orders of magnitude in $\kappa_\mathrm{H}/\kappa_\mathrm{C}$. Figure ~\ref{fig:heat2}a reveals an excellent semi-quantitative agreement between the direct heat measurement in MD simulations with the stochastic theory over the  parameter regime that we explore.  Notably, we remark that  our predictions  are done   without using any fitting parameter, i.e. we use in Eqs.~(\ref{eq:Qc}-\ref{eq:Qh}) the actual parameter values of the MD simulation and the friction coefficient estimated from the equilibrium fluctuations described above. Importantly, the MD simulation results reveal that the system acts as a heater whenever $\kappa_\mathrm{H}/ \kappa_\mathrm{C}< T_\mathrm{H}/T_\mathrm{C}$ and as a refrigerator (and heat-pump) when
\begin{equation}\label{cond1}
\kappa_\mathrm{H}/ \kappa_\mathrm{C}> T_\mathrm{H}/T_\mathrm{C},
\end{equation}
a feature that is supported by our theory.  This result rationalizes  the findings in Fig.~\ref{fig:heat}, which correspond to $\kappa_\mathrm{H}/ \kappa_\mathrm{C}=0.17$ (Fig.~\ref{fig:heat}a) and $\kappa_\mathrm{H}/ \kappa_\mathrm{C}=4.8$ (Fig.~\ref{fig:heat}b) which correspond respectively to the heater and refrigerator regimes as $T_\mathrm{H}/T_\mathrm{C}=1.2$.

We also compare in Fig.~\ref{fig:heat2}b  the theoretical prediction~\eqref{eq:Stot} for steady-state rate of entropy production (line) with its value  estimated from the MD simulations (symbols), with the latter evaluated by plugging in to $\langle \dot{S}_\mathrm{tot}\rangle =
\langle \dot{Q}_\mathrm{H}\rangle /T_\mathrm{H} +\langle \dot{Q}_\mathrm{C} \rangle/T_\mathrm{C}$, the MD values of the heat flows in the thermostats divided by their respective temperatures.

The estimate for the rate of entropy production obtained from the MD measurements is
in excellent agreement with the  theoretical expression given by~\eqref{eq:Stot}.  Around the  coupling values $\kappa_\mathrm{C} T_\mathrm{H} =\kappa_\mathrm{H} T_\mathrm{C}$ the heat flow and entropy production vanish, which we will refer to as a ``pseudo equilibrium'' point in the following.  In the SI we show that in this case, detailed balance holds, i.e., all probability currents vanish, despite the presence of a temperature gradient in the system together with a nonreciprocal coupling.

\section{Power and Performance}
\begin{figure}[t]
\centering
\includegraphics[width=0.28\textwidth]{Figures/COP_Ref.pdf}
	\caption{ Average coefficient of performance as a function of $\kappa_\mathrm{H}$  in the refrigeration mode: results from MD simulations (symbols) and theoretical prediction from stochastic thermodynamics [\eqref{eq:COP_r}, line].
	The horizontal line indicates the Carnot bound for the coefficient of performance, reached at $\kappa_\mathrm{H}= \kappa_\mathrm{C} T_\mathrm{H}/T_\mathrm{C}$, and given by   $ {\rm{COP}}  = T_\mathrm{C}/ (T_\mathrm{H} -T_\mathrm{C})$.
	The rest of the simulation parameters were set the same as in Fig.~\ref{fig:heat}.
	}
	\label{fig:COP}
\end{figure}
\subsection{Performance of the nano refrigerator}

In any refrigerator, getting the most heat from the temperature source is desirable by doing the least amount of work possible. Therefore, a suitable coefficient of performance (COP) is defined as the ratio of the heat taken from the low-temperature source to work done on the machine. A higher COP indicates a better economic performance of the system. In macroscopic systems, COP is usually in the range of $1-4$~\cite{borgnakke2022fundamentals}.
To quantify the net energetic  performance of the nano refrigerator,  we evaluate its coefficient of performance, defined by the average rate of heat that is extracted from the cold bath
divided by the average total power inputted into the system, COP$=|\langle \dot{Q}_\mathrm{C}\rangle |/\langle \dot{W}\rangle $.
From this definition and upon using  our theoretical  predictions for the heat transfer  given by  Eqs.~(\ref{eq:Qc},\ref{eq:Qh}), we predict that the COP follows
\begin{align}
\rm{COP}
\equiv
\frac{|\langle \dot{Q}_\mathrm{C}\rangle |}{\langle\dot {W}\rangle}
=
\frac{-\langle \dot{Q}_\mathrm{C}\rangle }{\langle\dot{Q}_\mathrm{H}\rangle+\langle\dot{Q}_\mathrm{C}\rangle}
 =
\frac{\kappa_\mathrm{C}}
{\kappa_\mathrm{H}-\kappa_\mathrm{C}} .
\label{eq:COP_r}
\end{align}
In \eqref{eq:COP_r}, the second equality follows from using the first law of thermodynamics for stationary states which in this case reads  $\langle\dot{W}\rangle-\langle \dot{Q}_\mathrm{C}\rangle-\langle \dot{Q}_\mathrm{H}\rangle=0$, whereas the third equality follows from   Eqs.~(\ref{eq:Qc},\ref{eq:Qh}).
Remarkably, our theoretical prediction for the COP   depends solely on the nonreciprocal coupling parameters $\kappa_\mathrm{H}$ and $\kappa_\mathrm{C}$.
Figure~\ref{fig:COP}  reveals a good agreement between the value of the $\mathrm{COP}$ estimated from the MD simulations (symbols) and the prediction from our stochastic theory  \eqref{eq:COP_r} (lines).  %
Interestingly the  agreement between simulation and theory is enhanced especially in far from equilibrium conditions, i.e.  for large  relative strengths of nonreciprocity $\kappa_\mathrm{H}/\kappa_\mathrm{C}\gg 1$. Close to the pseudo equilibrium $\kappa_\mathrm{H}/\kappa_\mathrm{C}=T_\mathrm{H}/T_\mathrm{C}$ point, the noise in the simulation results seems more pronounced. Towards this point, the theory predicts Carnot efficiency, which corresponds to a COP of $T_\mathrm{C}/(T_\mathrm{H}-T_\mathrm{C})=5$ in the present case.

\subsection{Fluctuations of power and performance}

We have shown so far that,  when $\kappa_\mathrm{H}/\kappa_\mathrm{C}>T_\mathrm{H}/T_\mathrm{C}$,  the nanoparticle system behaves like a refrigerator \textit{on average}, i.e. its steady-state average fluxes obey $\langle  \dot{Q}_\mathrm{C}\rangle>0 $, $\langle  \dot{Q}_\mathrm{H}\rangle<0 $ and $\langle  \dot{W}\rangle>0$.  Due to thermal fluctuations, the two-nanoparticle system can give rise to  transient values of the fluxes that do not obey the refrigerator constraints (e.g. transient values  $ \dot{Q}_\mathrm{C}<0 $, $  \dot{Q}_\mathrm{H}>0 $ and $  \dot{W}<0$ in a small time interval) as revealed in Langevin-dynamics models~\cite{rana2016anomalous}. In order to inspect such  fluctuation phenomena it is mandatory to evaluate  quantities such as the  power and  performance of the nano machine  along individual, short time intervals.
As a first approach in this direction, we evaluate the  fluctuations of the stochastic power from the MD simulations using the
positional fluctuations of the center of mass of the nanoparticles. In particular, we evaluate the stochastic power exerted in a small time interval $[t,t+\mathrm{d}t]$ using the expression from stochastic thermodynamics~\cite{Sekimoto2010} associated with the model given by~\eqref{eq:Le}
\begin{equation}\label{eq:stochpower}
\dot W  = \underbrace{\kappa_\mathrm{C} X_\mathrm{H} \circ \dot{X}_\mathrm{C}}_{=\displaystyle \dot{W}_\mathrm{C}} +  \underbrace{\kappa_\mathrm{H} X_\mathrm{C}\circ  \dot{X}_\mathrm{H}}_{=\displaystyle \dot{W}_\mathrm{H}}.
\end{equation}
Here, $\circ$  the Stratonovich product, whereas $ \dot{X}_\mathrm{C} =[X_\mathrm{C} (t+\mathrm{d}t) - X_\mathrm{C} (t)]/\mathrm{d}t$ and $ \dot{X}_\mathrm{H} =[X_\mathrm{H} (t+\mathrm{d}t) - X_\mathrm{H} (t)]/\mathrm{d}t$  are the  time-averaged velocities estimated from the positions of the nanoparticles. We have also introduced  $ \dot{W}_\mathrm{C}$ and $\dot{W}_\mathrm{H}$ as  the fluctuating power exerted on the cold (hot) nanoparticle in the interval $[t,t+\mathrm{d}t]$ respectively.  Note that in~\eqref{eq:stochpower} we take into account only the non-conservative forces exerted on each nanoparticle due to their nonreciprocal coupling.

\begin{figure}
    \centering
    \includegraphics[width=0.48\textwidth]{Figures/Figure5_Distributions-power-cop.pdf}
    \caption{
    (a) Distribution of the stochastic power exerted on the cold nanoparticle: result from the MD simulations (bars) and  analytical prediction given by~\eqref{eq:Pw} (solid blue line). Note that the lower panel in Figure 2 (b) displays the corresponding values of $W_\mathrm{C}$.

    (b) Distribution of the stochastic coefficient of performance defined by~\eqref{def:cop}: empirical density obtained  from the MD simulations
     (see Methods section) and theoretical prediction given by a   power law with exponent $-2$ (red line). For the second panel we took $\kappa_\mathrm{H}/\kappa_\mathrm{C}=7.27$ and fixed the rest of the parameter values the same as in Fig.~\ref{fig:heat2}.
    }
    \label{fig:power+efficiency_fluctuations}
\end{figure}

We now evaluate the stochastic power~\eqref{eq:stochpower} from our MD simulations over time intervals of duration $\mathrm{d}t=1$ps, for the parameter value  $\kappa_\mathrm{H}/\kappa_\mathrm{C}=4.8$ where we obtained a maximum heat  extraction on average from the cold bath. To this aim, we extract the empirical probability density for the  stochastic power  which reveals considerable fluctuations  (gray bars in Fig.~\ref{fig:power+efficiency_fluctuations}a). The distribution estimated from the MD simulations can be described with impressive accuracy using the  closed-form expressions  for the distribution of  $\dot W_\mathrm{C}$ and $\dot W_\mathrm{H}$
\begin{equation}\label{eq:Pw}
P(\dot W_j)=
\frac{1}{\mathcal{Z}_j} \,\exp \left [
 \frac{\beta_j}{\kappa_j}  \dot{W}_j
 -2\frac{\sqrt{\zeta_j \alpha_j}}{|\kappa_j|} |\dot{W}_j| \right],
\end{equation}
which  we derive analytically using the definition of the stochastic power~\eqref{eq:stochpower} and assuming the  effective stochastic model~\eqref{eq:Le}. In~\eqref{eq:Pw}, $j\in\{\mathrm{C,H}\}$ is the particle label,   $\mathcal{Z}_j$ is a normalization constant, and $\alpha_j, \beta_j, \zeta_j$ are functions of $\kappa_\mathrm{C},\kappa_\mathrm{H},T_\mathrm{C},T_\mathrm{H},\gamma,m$ (see Methods section for their explicit expressions).
The power distributions
have exponential tails and are slightly asymmetric. In the refrigerator mode, $P(\dot W_\mathrm{C})$ leans towards positive values, consistent with the net negative work value, while $P(\dot W_\mathrm{H})$ leans towards positive work values.
Furthermore, we find a good agreement between the MD simulation results and \eqref{eq:Pw}  throughout the parameter regime that we explore. The theoretical predictions yield a systematic  overestimation of the variance  which we attribute to the finite timestep $\mathrm{d}t$ in the MD simulations, while in the theoretical calculations we assume it to be infinitesimal (see Supplemental Material for  the values of the variance of the power for different values of $\kappa_\mathrm{H}/\kappa_\mathrm{C}$).

The previous analyses showing that the work and heat in a small time interval are highly fluctuating motivates us to investigate the finite-time fluctuations of the coefficient of performance of the nano machine. To quantify how much the efficiency fluctuates for individual trajectories, we consider the {\em stochastic} coefficient of performance defined as~\cite{joseph2021efficiency}
\begin{align}\label{def:cop}
    \epsilon_r =& \frac{|\dot{Q}_\mathrm{C}|}{\dot{W}}.
\end{align}
Note that since in general  $\langle|\dot{Q}_\mathrm{C}|\rangle/\langle\dot W\rangle\neq \langle|\dot{Q}_\mathrm{C}|/\dot{W}\rangle  $, the ensemble average of $\epsilon_r$ does not coincide  with the average COP given by~\eqref{eq:COP_r}.
Figure~\ref{fig:power+efficiency_fluctuations}b displays the empirical distribution of the stochastic COP obtained  from the MD simulations, which develops a fat tail with values that can exceed significantly  the Carnot value~$\epsilon_\mathrm{C}   =T_\mathrm{C}/ (T_\mathrm{H}-T_\mathrm{C})=5$. Such super-Carnot performance achieved in short time intervals have also been reported in experimental conditions and theoretical models of nanoscopic heat engines ~\cite{polettini2015efficiency,verley2014unlikely,martinez2017colloidal,park2016efficiency} and  nano refrigerators~\cite{rana2016anomalous}; they result from  rare events where the refrigerator works transiently as a heater  reversing the work flux with respect to its average behaviour.
Furthermore we find that the distribution of $\epsilon_r$ extracted form the MD simulations  follows in good approximation a power law
\begin{equation}
    P(\epsilon_r)  \simeq  |\epsilon_r|^{-2},
    \label{eq:pl}
\end{equation}
see red line in Fig.~\ref{fig:power+efficiency_fluctuations}b.
The power-law behavior that we find reinforces the critical  significance of thermal fluctuations for our nano machine system setup. Remarkably, we find that the power law~\eqref{eq:pl} is in good agreement with the numerical results for all values of $\kappa_\mathrm{H}$ and $\kappa_\mathrm{C}$ that we explored, suggesting a universal scaling behavior as predicted by previous theoretical  work within the realm of stochastic efficiency~\cite{polettini2015efficiency}.

\subsection{Uncertainty relations}

The results in previous sections revealed the instrumental role of stochastic thermodynamics to establish design principles for the parameter values of the nano machine to achieve prescribed values of the net power and efficiency.
In the following, we investigate how one can use principles from stochastic thermodynamics --namely the so-called thermodynamic uncertainty relations~\cite{barato2015thermodynamic,hasegawa2019uncertainty,horowitz2020thermodynamic}--  to put fundamental constraints  that regulate the trade-off between dissipation and precision of the nano machine.

A suitable measure to quantify the strength of the power fluctuations is the \textit{uncertainty} of the power, defined by the variance over the squared mean. High values of uncertainty indicate that the dynamics is essentially dominated by  fluctuations.
We have measured the uncertainty of the total power $\dot W = \dot W_\mathrm{C} + \dot W_\mathrm{H}$ in the MD simulations from individual trajectories, by making a statistics over the extracted values of the work. Figure \ref{fig:power_uncertainty} shows the results for the uncertainty of the total power (green symbols), as well as the power uncertainties of $\dot{W}_\mathrm{C}$ and $\dot{W}_\mathrm{H}$ separately (blue and the red symbols). Remarkably, the uncertainties reach extremely high values around the pseudo equilibrium point, and even seem to diverge at $\kappa_\mathrm{H} /\kappa_\mathrm{C} = T_\mathrm{H}/T_\mathrm{C}$.
As we show below, this blow-up can be understood by making use of a recently-developed trade-off relation between the precision of thermodynamic currents and the rate of entropy production.

In the field of stochastic thermodynamics, recently a universally class of results ---often called e thermodynamic uncertainty relation--- governing Markovian nonequilibrium stationary states  were derived, see e.g.~\cite{barato2015thermodynamic,hasegawa2019uncertainty,horowitz2020thermodynamic}.
Such laws connect the uncertainty of a current, quantified by its signal-to-noise ratio, with the total thermodynamic cost, measured by the steady-state rate of entropy production. In particular, the thermodynamic uncertainty relation in Ref.~\cite{pietzonka2017finite} implies  that the uncertainty of the finite-time power fluctuations of stationary Markovian processes is always bounded from below by $2k_\mathrm{B}$ over the mean total entropy production, i.e.,
\begin{equation}\label{eq:TUR}
    \frac{\mathrm{Var}({\dot{W}})}{\langle \dot W \rangle^2} \geq \frac{2 k_{\rm B}}{\langle  \dot S_\mathrm{tot}\rangle}.
\end{equation}
Equation~\eqref{eq:TUR} reveals that there exists  a minimal thermodynamic cost associated with  achieving a certain precision of the power exerted by the controller.
Applying this law to the present case, we find from \eqref{eq:Stot} that $\langle \dot S_\mathrm{tot} \rangle \to 0$ at $\kappa_\mathrm{H} /\kappa_\mathrm{C} \to T_\mathrm{H}/T_\mathrm{C}$, thus, the power uncertainty indeed diverges at the pseudo equilibrium point.
We have complemented the MD simulation results in Fig.~\ref{fig:power_uncertainty} with a black line showing the lower bound $2 k_{\rm B}/\langle  \dot S_\mathrm{tot}\rangle$ according to~\eqref{eq:TUR}, which provides the first test of thermodynamic uncertainty relations in nonequilibrium MD setups.

From the simulation results shown in Fig. \ref{fig:power_uncertainty}, we further detect a minimum in the total power uncertainty in the refrigerator regime. Interestingly, the position of this minimum (that is $\kappa_\mathrm{H}/\kappa_\mathrm{C} \approx 2$) roughly coincides with the  $\kappa_\mathrm{H}/\kappa_\mathrm{C}$ value where the amount of extracted heat is maximal. Thus, there is a regime where the refrigeration is maximal while at the same time, the power to sustain the refrigerator is as precise as possible.

\begin{figure}
    \centering
    \includegraphics[width=0.4\textwidth]{Figures/TUR_w_1.pdf}
    \caption{
    Thermodynamic uncertainty relation for power fluctuations. Relative uncertainty $\text{Var}(\dot{W})/\langle \dot{W}\rangle^2$ of the stochastic power exerted on the cold $\dot{W}_\mathrm{C}$ (blue filled circles) and the hot $\dot{W}_\mathrm{H}$  (red filled squares) nanoparticles, computed using in Eq.~\eqref{eq:stochpower}  the output of our MD simulations. The uncertainty of the total power on the two-nanoparticle setup $\dot{W}=\dot{W}_\mathrm{C}+\dot{W}_\mathrm{H}$ is show by green crosses. The blue and red lines are given by our theoretical predictions for the relative uncertainty of $\dot{W}_\mathrm{C}$ and $\dot{W}_\mathrm{H}$, respectively (see Supplemental for the full derivation).
    The relative uncertainty of the stochastic power lies above  2 divided by   the rate of entropy production [black symbols, MD simulations; black line, theoretical prediction given by \eqref{eq:Stot}], in agreement with the thermodynamic uncertainty relation
    \eqref{eq:TUR}.
    }
    \label{fig:power_uncertainty}
\end{figure}

\section{Discussion}

We have shown with molecular dynamics simulations that it is possible to attain a net heat flow from a cold to a hot bath by connecting two  nanoparticles immersed in fluid containers through nonreciprocal forces. Such nonreciprocal forces could be realized in the laboratory upon using, e.g., feedback traps~\cite{Jun2012} which allow to exert in real time forces  based on measurements of only the position of the center of mass of the nanoparticles. This represents an advantage with respect to traditional Maxwell-demon approaches where the measurement of position and velocities of the bath molecules is required to attain a heat flow from a cold to a hot thermal bath. Notably, our setup  is also advantageous with respect to previous theoretical proposals that require the usage of athermal fluctuations together with nonlinear forces between the nanoparticles~\cite{Kanazawa2013}.

In our simulations we have studied the case of silver nanoparticles immersed in argon, however we expect this effect to be generic for a class of systems whose dynamics can be described by linear underdamped Langevin equations. For such class of systems, we have revealed by using stochastic thermodynamics the necessary conditions (i.e., design principles) that ensure a reverse heat flow (from cold to hot) and provided predictions on key thermodynamic properties such as the net heat transfer and the coefficient of performance. Our theoretical results reveal that the refrigeration effect generated by the nonreciprocal coupling is robust. It can be achieved in a broad parameter range as long as the dynamics is stable (see Appendix).

This work demonstrates fruits of the bridge between molecular dynamics with stochastic thermodynamics, namely the possibility to establish quantitative criteria to control the statistics of thermodynamic fluxes in nanoparticle-based thermal machines. We have developed analytical formulae that describe the average heat fluxes between the nanoparticles, the entropy production rate and the coefficient of performance. Moreover, we have tackled analytically the statistics of the power and the coefficient of performance of the nano machine over  short time intervals within fluctuations play a prominent role. From the theoretical perspective, we expect that stochastic-thermodynamic approaches could shed light in the future on optimal design of nano refrigerator devices taking into account e.g. the effect of delay and/or memory in the exertion of nonreciprocal forces that may play a key role in realistic experimental scenarios.

In the current era when efficient energy conversion is of paramount importance,  optimizing thermal transport in both classical and quantum systems has become a topic of intense theoretical and experimental interest for nanotechnology~\cite{baroni2017,baroni2019}. Furthermore, the possibility of realizing thermal conduits in various biological systems have been proposed to optimize thermal networks and information transfer~\cite{klinman2021,klinman2022}. Transient reverse heat flow in these contexts may confer biological machinery with enhanced functionality for energy harvesting. In this regard, it would be interesting to theoretically investigate whether our theoretical and atomistic setup can be used to generate heat flows in particles immersed in more complex environments including viscoelastic fluids, or even active (e.g. bacterial) baths.

\begin{acknowledgments}
SL acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -- through the project 498288081.
\end{acknowledgments}

\section{Methods}

\subsection{MD simulations}

We employed the LAMMPS for performing all out MD simulations~\cite{plimpton1995fast}. The interatomic force between atoms was accounted by a Lennard-Jones (LJ) potential function,
\begin{equation}
\phi(r_{ij})=4\varepsilon[(\sigma/r_{ij} )^2- (\sigma/r_{ij} )^6 ],
\end{equation} where $r_{ij}$ is the interatomic distance between atom $i$ to atom $j$,  $\varepsilon$ is the depth of the potential well, and $\sigma$  the distance at which the particle--particle potential energy vanishes. The parameter values $\varepsilon$ and $\sigma$ of both Argon--Argon and Copper--Copper interactions are summarized in Table~\ref{tab:1}.
\begin{table}[ht]
\centering
\begin{tabular}[t]{lcc}
\hline
& $\varepsilon (\mathrm{eV})$ & $\sigma (\AA)$ \\
\hline
Ar – Ar&0.0104&3.405\\
Cu – Cu&0.4093&2.338\\
\hline
\end{tabular}
\caption{LJ interaction parameters \cite{zarringhalam2019effects}\label{tab:1}}
\end{table}%
For the interatomic forces between Argon and Copper atoms, we employed the Lorentz-Berthelot mixing rules~\cite{delhommelle2001inadequacy}, i.e.
\begin{eqnarray}
\varepsilon_{12} = \sqrt{\varepsilon_1\: \varepsilon_2} , \hspace{1cm} \sigma_{12}=  (\sigma_2+ \sigma_1)/2 .
\end{eqnarray}

As shown in Fig.~\ref{fig:MD}, our setup consists of two containers of cold and hot (defined by blue and red colours respectively) Argon-Copper mixture, each with volume $11 \rm{nm}^3$ and containing $25280$ Ar fluid atoms with total mass $m _{\rm Ar}= 1.67\times 10^{-21} \rm{kg}$. The  nanoparticles are made each of  $186$ Cu  atoms  with nanoparticle mass $m= 1.96\rm{e}-23 \rm{kg}$ and radius $r = 1.44 \rm{nm}$. To prevent  interactions between the two containers and separate cold and hot baths, one-atom-thick Cu walls are placed between the containers, in the edges and in the middle  of the simulation box along the $X$ direction. Ar atoms of one container do not interact with the Ar atoms inside the other container.
Periodic boundary condition were applied in the $Y$ and $Z$ directions while the $X$ direction is constrained by the walls.
The simulations were carried out for $15-30 \rm{ns}$. Within this time the temperature of fluids at both hot and cold baths were controlled at $T_\mathrm{H}=120\rm{K}$ and $T_\mathrm{C}=100\rm{K}$ using  Nos\'e---Hoover thermostats (NVT)  with a coupling time-constant of $1\rm{ps}$. We have also checked that the reversed heat flow setup is not sensitive to the choice of Nos\'e---Hoover as it is also reproduced using a Langevin-based thermostat as well (data not shown).  A time step of $1 \rm{fs}$ was chosen for the MD simulations and the data was sampled every $1\rm{ps}$.

\subsection{Equations of motion}
The stochastic model that we use is given by two coupled underdamped Langevin equations for the  positions $x_\mathrm{C}$ and $x_\mathrm{H}$ of the center of mass of the nanoparticles along the $x$ axis
\begin{eqnarray} \label{eq:LeStart}
\hspace{-0.5cm}m  \ddot{x}_\mathrm{C}  &=& -  \gamma_\mathrm{C} \dot{x}_\mathrm{C} - \kappa x_\mathrm{C} + \kappa_\mathrm{C} (X_\mathrm{H} -x_\mathrm{C}-R ) + \xi_\mathrm{C}\\
\hspace{-0.5cm}m  \ddot{x}_\mathrm{H}&=& - \gamma_\mathrm{H} \dot{x}_\mathrm{H}   -\kappa(X_\mathrm{H}-R) -\kappa_\mathrm{H} (
x_\mathrm{H} - x_\mathrm{C} - R )+ \xi_\mathrm{H}.\nonumber\\
\end{eqnarray}
The parameters are defined in the Main Text below Eq.~\eqref{eq:Le}.
Introducing the variables,
${x}_\mathrm{H}= {X}_\mathrm{H} - R$ and ${x}_\mathrm{C}= {X}_\mathrm{C}$, we can simplify \eqref{eq:LeStart} to
\begin{align}
m \ddot{X}_\mathrm{C}+\gamma_\mathrm{C} \dot{X}_\mathrm{C} &=  -(\kappa_\mathrm{C} + {\kappa} ) X_\mathrm{C} + \kappa_\mathrm{C} X_\mathrm{H} + \xi_\mathrm{C}
\\
m \ddot{X}_\mathrm{H}+\gamma_\mathrm{H} \dot{X}_\mathrm{H} &
=
\kappa_\mathrm{H}
X_\mathrm{C} -({\kappa}+\kappa_\mathrm{H}) X_\mathrm{H} + \xi_\mathrm{H},
\end{align}
as used throughout in the Main Text.
In all our MD simulations we set the stiffness of the  traps  to the value ${\kappa}=1\,\mathrm{eV/(nm)^2} = 0.16021\,\mathrm{J/m^2}=1.16 k_\mathrm{B}T_\mathrm{C}/\AA ^2$.
We further fix $\kappa_\mathrm{C}=10 \kappa $ and vary $\kappa_\mathrm{H}$.

\subsection{Estimation of the effective friction coefficient}
The noise terms as well as the friction forces appearing in the Langevin equations are not explicitly present in the MD simulations, but they emerge implicitly from the interactions between nanoparticle and surrounding bath particles. %
In our stochastic model, we assume Stokes law, in particular, instantaneous effective friction force that is linearly proportional on the instantaneous velocity, see Eq.~\eqref{eq:LeStart}.
We obtain estimates for the values of the corresponding friction coefficients $\gamma_\mathrm{C,H}$ from MD simulations in the case of no coupling, $\kappa_\mathrm{C,H}=0$, via two distinct routes.
First, we measure the velocity and positional autocorrelation functions in an equilibrium MD simulation, and fit the corresponding analytical expression taken from Ref.~\cite{wang1945theory}:
\begin{eqnarray}
\!\!\frac{\langle {X}_\mu (t){X}_\mu(t+\tau)\rangle}{\langle {X}_\mu^2 (t)\rangle}\!\! &=& \!\!\exp\left(-\frac{\gamma_\mu \tau}{2 m}\right) \left[ \cos( \omega_\mu \tau) + \frac{\gamma\sin(\omega_\mu \tau)}{2 m \omega_\mu}  \right]\nonumber\\
\!\!\frac{\langle {\dot{X}}_\mu (t){\dot{X}}_\mu(t+\tau)\rangle}{\langle {\dot{X}}_\mu^2 (t)\rangle}\!\! &=& \!\! \exp\left(-\frac{\gamma_\mu \tau}{2 m}\right) \left[ \cos( \omega_\mu \tau) -\frac{\gamma \sin(\omega_\mu \tau )}{2 m \omega_\mu} \right],\nonumber\\ \label{acfs}
\end{eqnarray}
with  $\omega_\mu^2 =(\kappa/m)-(\gamma_\mu/m)^2$, and we recall $\mu=\{\rm{C,H}\}$.
Equations~\eqref{acfs} reproduce our numerical estimates for both the position and velocity autocorrelation functions in equilibrium conditions (see the SI for further details). Fitting~\eqref{acfs} to our numerical results by setting $m$, $\kappa$ to their input values in the simulations, we extract the estimates  $\gamma_\mathrm{C}\simeq \gamma_\mathrm{H}  \sim 3.5 \times 10^{-12}$ kg/s for the effective friction coefficient.
 We use this estimate for $\gamma$ and set $\gamma_\mathrm{C} = \gamma_\mathrm{H}$.
Second, we have performed a dragging experiment, pulling the nanoparticle through the bath and measure the resulting velocity (see Appendix for further details), which yields an estimate of $\gamma_\mathrm{C} \approx \gamma_\mathrm{H} \approx 3\times 10^{-12} \text{kg/s}$.
Both routes yield consistent results. We use $\gamma=\gamma_\mathrm{C} = \gamma_\mathrm{H} = 3.5\times 10^{-12} \text{kg/s}$ throughout the paper.

\subsection{Definition of heat and work in Stochastic Thermodynamics}
We address the thermodynamics of the nanoparticle system by applying the framework of stochastic thermodynamics~\cite{Sekimoto2010,Seifert2012} to our stochastic model. For convenience,  we first rewrite Eqs. (1,2) as a two-dimensional Langevin equation
\begin{equation}
    m\ddot{{\bf X}} = -\gamma \dot{{\bf X}} -\nabla E+ {\bf F}_{\rm nr} + \Xi
    \label{LE2}
\end{equation}
where ${\bf X}= ( X_\mathrm{C} , X_\mathrm{H} )^{\mathsf{T}}$, $\Xi=(\xi_\mathrm{C},\xi_\mathrm{H})^{\mathsf{T}}$, and $\nabla = (\partial_{X_\mathrm{C}},\partial_{X_\mathrm{H}})^{\mathsf{T}}$, with ${\mathsf{T}}$ denoting transposition. The term ${\bf F}_{\rm nr}=(\kappa_\mathrm{C} (X_\mathrm{H}-R), \kappa_\mathrm{H} X_\mathrm{C})^{\mathsf{T}}$ is the nonreciprocal force acting on the system, which is non-conservative. On the other hand,  the energy $E$ of the nanoparticles
\begin{equation}
    E = \frac{m}{2}(\dot{X}_\mathrm{C}^2 + \dot{X}_\mathrm{H}^2) + \left(\frac{\kappa+\kappa_\mathrm{C}}{2}\right) X_\mathrm{C}^2 + \left(\frac{\kappa+\kappa_\mathrm{H}}{2}\right)   (X_\mathrm{H}-R)^2
\end{equation}
is a function  of the instantaneous values of the particles' positions and velocities, and thus a quantity that fluctuates during a simulation,  i.e. a stochastic process.
  Following Sekimoto~\cite{Sekimoto2010}, the energy change in a small time interval $[t,t+\rm{d}t]$ can be written as
\begin{equation}\label{def:u}
    \mathrm{d}E = \nabla E({\bf X}) \circ \mathrm{d}{\bf X} + m\ddot{{\bf X}}\circ \mathrm{d}{{\bf X}},
\end{equation}
where $\circ$ denotes the Stratonovich product. Here, $\mathrm{d}{\bf X}= ( \mathrm{d}X_\mathrm{C} , \mathrm{d}X_\mathrm{H} )^T  $ is the stochastic increment of the nanoparticles' positions which follows from the Langevin dynamics~\eqref{LE2}. In $[t,t+\rm{d}t]$ the stochastic work done on the system is given by the non-conservative force times the displacement of the nanoparticles, which we can split as
\begin{equation}\label{def:w}
    \mathrm{d}W = {\bf F_{\rm nr}}\circ  \mathrm{d}{{\bf X}} = \underbrace{\kappa_\mathrm{C} (X_\mathrm{H}-R) \circ  \mathrm{d}X_\mathrm{C}}_{=\displaystyle \mathrm{d}W_\mathrm{C}} +  \underbrace{\kappa_\mathrm{H} X_\mathrm{C}\circ  \mathrm{d}X_\mathrm{H}}_{=\displaystyle \mathrm{d}W_\mathrm{H}},
\end{equation}
where $\mathrm d W_\mathrm{C}$ ($\mathrm d
W_\mathrm{H}$) is the work done on the cold (hot) nanoparticle.
Dividing the stochastic work $\mathrm{d}W$ applied to the particle along a trajectory of length by the trajectory length $\mathrm{d}t$ yields the power $\dot{W}= \mathrm{d}{W}/\mathrm{d}t$.

Similarly, the stochastic heat dissipated in the same time interval
\begin{equation}
    \mathrm{d} Q = \mathrm{d}E - \mathrm{d}W =  \underbrace{(\gamma_\mu \dot{X}_\mathrm{C} - \xi_\mathrm{C} )\circ \mathrm{d}{X}_\mathrm{C}}_{=\displaystyle \mathrm{d}Q_\mathrm{C}} + \underbrace{(\gamma_\mu \dot{X}_\mathrm{H} - \xi_\mathrm{H} )\circ \mathrm{d}{X}_\mathrm{H}}_{=\displaystyle \mathrm{d}Q_\mathrm{H}} ,
    \label{def:heat}
\end{equation}
which ensures that the first law is satisfied for every single trajectory traced by the system.
Note that we use the thermodynamic sign convention: $W>0$ ($W<0$) when work is exerted (extracted) from the system and $Q>0$ ($Q<0$) when heat is dissipated from (absorbed by) the system to (from) its environment.
Dividing the stochastic heat dissipated along a trajectory of length by the trajectory length $\rm{d} t$ yields the stochastic heat dissipation rate $\dot{Q}= {\rm{d}} {Q}/\mathrm{d}t$.

\subsection{Power fluctuations}
The power fluctuations of each nanoparticle predicted by our stochastic model are given in   Eq. \eqref{eq:Pw} as closed-form expressions for $P(\dot W_j)$ with $j\in \{ \text{C,H} \}$. The distributions $P(\dot W_j)$  explicitly depend on $\alpha_j,\beta_j,\zeta_j$, which are in turn functions of the covariances of the positions and velocities of the $j$ nanoparticle  and the $l\neq j$ nanoparticle. Specifically, %
\begin{align}\label{def:alpha}
\alpha_j &= \frac{1}{2 \langle X_l^2 \rangle (1-\psi_j^2)},
~~
\beta_j = \frac{\psi_j}{\sqrt{\langle  X_l^2 \rangle\langle \dot X_j^2 \rangle}(1-\psi^2_j)},
\end{align}
and
\begin{align}
~~\zeta_j = \frac{1}{2 \langle \dot X_j^2 \rangle (1-\psi_j^2)}, ~~
\psi_j=\frac{\langle X_l \dot{X}_j \rangle }{\sqrt{\langle  X_l^2 \rangle\langle \dot X_j^2 \rangle}},
\end{align}
with
\begin{align}
\langle  {X}_l^2 \rangle
&=
k_\mathrm{B}
\frac{
m T_l [\kappa_j^2 (\kappa_j+\kappa_l) + \kappa (\kappa_j^2 + \kappa_l^2)]}
{\kappa(\kappa + \kappa_j + \kappa_l)[m \left[(\kappa_j + \kappa_l)^2 +2 (2\kappa + \kappa_j + \kappa_l)\right]]}
\nonumber \\
&+k_\mathrm{B}
\frac{ 2\gamma^2 [\kappa_l^2 T_j + T_l(\kappa+\kappa_j)(2\kappa+\kappa_j) +\kappa \kappa_l T_l]}
{\kappa(\kappa + \kappa_j + \kappa_l)[m \left[(\kappa_j + \kappa_l)^2 +2 (2\kappa + \kappa_j + \kappa_l)\right]]}
\nonumber \\
&+
k_\mathrm{B}
\frac{
m T_j \kappa_l^2 [2 \kappa + \kappa_j + \kappa_l]
}
{\kappa(\kappa + \kappa_j + \kappa_l)[m \left[(\kappa_j + \kappa_l)^2 +2 (2\kappa + \kappa_j + \kappa_l)\right]]},
\end{align}
\begin{align}
\langle  \dot{X}_j^2   \rangle
=
k_\mathrm{B}\frac{ \kappa_j (\kappa_j T_l-\kappa_l T_j)}{(m/2) (\kappa_j + \kappa_l)^2+ \gamma ^2 (2\kappa+\kappa_j+\kappa_l)}+\frac{k_\mathrm{B} T_j}{m},
\end{align}
and
\begin{align}
\langle  {X}_l{\dot{X}}_j \rangle
=
k_\mathrm{B}
\frac{
2\gamma(\kappa_j T_l - \kappa_l T_j )}
{m (\kappa_j+\kappa_l)^2+ 2(2\kappa+ \kappa_j \kappa_l) \gamma ^2}.
\end{align}
From the power distributions $P(\dot W_j)$, we can further deduce analytical expressions for the variance of the power:
\begin{align}\label{eq:PwVar}
    \text{Var}(\dot W_j)& =
    \frac{1}{\mathcal{Z}_j}\!\left\{ \frac{2\kappa_j^3}{(2\sqrt{\alpha_j \zeta_j}+\beta_j)^3} +
    \frac{2\kappa_j^3}{(2\sqrt{\alpha_j \zeta_j}-\beta_j)^3}

    \right.\nonumber \\&- \left.
    \left[ \frac{\kappa_j^2}{(2\sqrt{\alpha_j \zeta_j}+\beta_j)^2} +
    \frac{2\kappa_j^2}{(2\sqrt{\alpha_j \zeta_j}-\beta_j)^2}\right]^2 \right\} .
\end{align}
See the Supplemental Material for further details and the explicit mathematical derivations.


\section{Decomposition of the $k$-Bonacci Polyominoes}

Let ${\mathcal P}_{n,k}^{(1)}$ and  ${\mathcal P}_{n,k}^{(2)}$ denote the sets of $k$-bonacci polyominoes with $n$ columns, whose last column has height $1$ and 2, respectively. It is clear that ${\mathcal P}_{n,k}={\mathcal P}_{n,k}^{(1)} \cup {\mathcal P}_{n,k}^{(2)}$. If $P\in {\mathcal P}_{n,k}^{(1)}$, then $P$ is a unit square or $P=P'\square$, where $P' \in {\mathcal P}_{n-1,k}={\mathcal P}_{n-1,k}^{(1)} \cup {\mathcal P}_{n-1,k}^{(2)}$ $(n>1)$.  If $P\in {\mathcal P}_{n,k}^{(2)}$, then $P=P'C_j$, where $P' \in {\mathcal P}_{n-j,k}^{(1)}$ (possibly empty) and  $C_j$ is a concatenation of $j$ columns of height 2, for $1\leq j < k$.  Figure~\ref{deco} illustrates the above decomposition.

\begin{figure}[ht]
\centering
\includegraphics[scale=1]{Deco.eps}
\caption{Decomposition of a $k$-bonacci polyomino.}
\label{deco}
\end{figure}

\subsection{Area and Semiperimeter}

In this section we study the distribution of the area and semiperimeter    in ${\mathcal P}_{n,k}$.
Let $P$ be a $k$-bonacci polyomino.  We denote  by ${\texttt{area}}(P)$  the number of cells of $P$ and by
${\texttt{sper}}(P)$ the semiperimeter, that is the half of  the  perimeter of $P$ (since  the perimeter of a $k$-bonacci polyomino is always an even number).

 Define the generating function
$$F_k(x,p,q):=\sum_{n\geq 1} x^n \sum_{P\in {\mathcal P}_{n,k}}p^{{\texttt{sper}}(P)}q^{{\texttt{area}}(P)}, \quad k\geq 2,$$
where $x$ marks the length of the corresponding $k$-bonacci word, i.e., the number of cells in the bottom row of a polyomino.
Analogously, we introduce the generating functions
$$F_{k,j}(x,p,q):=\sum_{n\geq 1} x^n \sum_{P\in {\mathcal P}_{n,k}^{(j)}}p^{{\texttt{sper}}(w)}q^{{\texttt{area}}(P)},  \quad \text{for} \quad j=1, \,  2.$$
It is clear that
\begin{align}\label{GF}
F_k(x,p,q)=F_{k,1}(x,p,q)+F_{k,2}(x,p,q).
\end{align}
In Theorem~\ref{teoareaper1} we give a rational expression for the generating function $F_k(x,p,q)$.
\begin{theorem}\label{teoareaper1}
The generating function for $k$-bonacci polyominoes with respect to
the length of the corresponding word,
semiperimeter and area (marked respectively by $x$, $p$ and $q$)
is
$$F_k(x,p,q)=\frac{p^2 ((q + p q^2 ) x - (p q^3 - p^2 q^3) x^2 - p^k q^{2 k} x^k -
   p^{k+1} q^{2 k+1} x^{k+1})}{1 - (p q + p q^2) x + (p^2 q^3 - p^3 q^3) x^2 +
 p^{k+2} q^{2 k+1} x^{k+1}}.$$
\end{theorem}
\begin{proof}
From the decomposition given in Figure~\ref{deco} we have the functional equations
\begin{align*}
F_{k,1}(x,p,q)&=\underbrace{p^2qx}_{(1)} + \underbrace{pqx(F_{k,1}(x,p,q)+F_{k,2}(x,p,q))}_{(2)+(3)} \\
F_{k,2}(x,p,q)&=\underbrace{\sum_{j=1}^{k-1}p^{i+2}q^{2i}x^j}_{(4)} +\underbrace{\left(\sum_{j=1}^{k-1}p^{i+1}q^{2i}x^j\right) F_{k,1}(x,p,q)}_{(5)}\\
&=\frac{p^2(pq^2x-p^kq^{2k}x^k)}{1-pq^2x} + \frac{p(pq^2x-p^kq^{2k}x^k)}{1-pq^2x}F_{k,1}(x,p,q).
\end{align*}
To make understanding easier, cases (1) to (5) in Figure~\ref{deco} are indicated below their corresponding terms.
Solving the above system of equation and from \eqref{GF} we obtain the desired result.
\end{proof}
For example, the series expansion of the generating function $F_{3}(x,p,q)$ begins with
\begin{multline*}
(p^2 q + p^3 q^2) x + (p^3 q^2 + 2 p^4 q^3 +
    p^4 q^4) x^2 + \bm{(p^4 q^3 + 3 p^5 q^4 + 2 p^5 q^5 +
    p^6 q^5) x^3} \\
    + (p^5 q^4 + 4 p^6 q^5 + 3 p^6 q^6 + 3 p^7 q^6 +  2 p^7 q^7) x^4 + \cdots
    \end{multline*}
Figure~\ref{grafo2} shows the weights of the  $3$-bonacci polyominoes   corresponding to the bold coefficient in the above series.

\begin{figure}[ht]
\centering
\includegraphics[scale=1]{Grafo2.eps}
\caption{Weights for polyominoes in ${\mathcal P}_{3,3}$.}
\label{grafo2}
\end{figure}

\begin{corollary}\label{coroarea}
The generating functions for the total area and semi-perimeter within all of the members in ${\mathcal P}_{n, k}$ are
$$A_k(x)=\frac{3 x - 2 k x^k - 2 (2 - k) x^{k+1} + x^{2 k+1}}{(1 - 2 x + x^{k+1})^2}$$
and
$$P_k(x)=\frac{5 x - 5 x^2 - (2 + k) x^k - (1 - k) x^{k+1} + 4 x^{k+2} - x^{2k+2}}{(1 - 2 x + x^{k+1})^2},$$
respectively.
\end{corollary}
\begin{proof}
From the definition of  $F_k(x,p,q)$ we have the equalities
$$A_k(x)=\left.\frac{\partial F_k(x,1,q)}{\partial q}\right|_{q=1} \quad\text{and}\quad P_k(x)=\left.\frac{\partial F_k(x,p,1)}{\partial p}\right|_{p=1}.$$
Therefore, from Theorem~\ref{teoareaper1} we obtain the desired results.
\end{proof}
In particular, for $k = 2, 3,$ and $4$ we obtain the following generating functions:
\begin{align*}
A_2(x)&=\frac{x (3 + 2 x + x^2)}{(1 -  x -x^2)^2},\\
 A_3(x)&=\frac{x (3 + 6 x + 3 x^2 + 2 x^3 + x^4)}{(1 -  x -x^2-x^3)^2},\\
 A_4(x)&=\frac{x (3 + 6 x + 9 x^2 + 4 x^3 + 3 x^4 + 2 x^5 + x^6)}{(1 -  x -x^2-x^3-x^4)^2},\\
P_2(x)&=\frac{x (5 + x - 2 x^2 - x^3)}{(1 -  x -x^2)^2},\\
 P_3(x)&=\frac{x (5 + 5 x - 3 x^3 - 2 x^4 - x^5)}{(1 -  x -x^2-x^3)^2},\\
 P_4(x)&=\frac{x (5 + 5 x + 5 x^2 - x^3 - 4 x^4 - 3 x^5 - 2 x^6 - x^7)}{(1 -  x -x^2-x^3-x^4)^2}.
   \end{align*}

We will apply the same technique  to obtain other corollaries  throughout the article.

\subsection{Area and Semiperimeter of the Fibonacci polyominoes}

For the Fibonacci polyominoes ($k=2$) we can give some additional results.  Let $t_n(p,q)$ denote the $n$-th coefficient of  $F_2(x,p,q)$, that is,
$$t_n(p,q):= \sum_{P\in {\mathcal P}_{n,2}}p^{{\texttt{sper}}(P)}q^{{\texttt{area}}(P)}.$$

\begin{theorem}\label{fibopoly}
For all $n\geq 3$ we have
$$t_n(p,q)=pqt_{n-1}(p,q) + p^3q^3t_{n-2}(p,q),$$
with the initial values $t_1(p,q)=p^2 q + p^3 q^2$ and $t_2(p,q)=p^3 q^2 + 2 p^4 q^3$. Moreover,  for all $n\geq 1$ we have the combinatorial formula
\begin{align}\label{closedf2b}
t_n(p,q)=\sum_{i=0}^{\lfloor \frac{n+1}{2} \rfloor} \binom{n+1-i}{i}p^{n + i + 1} q^{n + i}.
\end{align}
\end{theorem}
\begin{proof}
Let $P$ be a Fibonacci polyomino with $n$ columns ($n\geq 3$).  If the last column has height 1, then the number of this kind of polyominoes is given by  $qpt_{n-1}(p,q)$. On the other hand, if the last column has height 2, then the previous column has height 1, and these polyominoes are counted by the polynomial  $p^3q^3t_{n-2}(p,q)$.  Hence the total number of Fibonacci polyominoes is given by $pqt_{n-1}(p,q) + p^3q^3t_{n-2}(p,q)$. Finally, we can use mathematical induction and the recurrence relation to prove the combinatorial identity.  This formula can be also proved by means of Zeilberger's creative telescoping method~\cite{ABZ}. Denote the summand on the right side
of the equality in \eqref{closedf2b} by $F(n,i)$, that is
$$F(n,i):= \binom{n+1-i}{i}p^{n + i + 1} q^{n + i}.$$
By the Zeilberger algorithm, $F(n,i)$ satisfies the relation
\begin{align}\label{rel1}
F(n+2,i)-pqF(n+1,i)-p^3q^3F(n,i)=G(n,i+1)-G(n,i),
\end{align}
with the certificate
$$R(n,i)=-\frac{i (2 - i + n) p^2 q^2}{(2 - 2 i + n) (3 - 2 i + n)}.$$
That is, $R(n,i)=G(n,i)/F(n,i)$ is a rational function in both variables.
If $f(n)$ denotes the right sum in the equality \eqref{closedf2b}, then summing both sides of \eqref{rel1} with respect to $i$ yields $f(n+2)-pqf(n+1)-p^3q^3f(n)=0$.  The sequences $f(n)$ and $t_n(p,q)$ satisfy the same recurrence relation and  have the same initial values, therefore these sequences coincides for all positive integers $n$.

\end{proof}

\begin{corollary}
The  total area for all Fibonacci polyominoes in ${\mathcal P}_{n,2}$ is
$$\sum_{i=0}^{\lfloor \frac{n+1}{2} \rfloor} \binom{n+1-i}{i}(n + i)=\frac{1}{5}(6nF_{n+2}+(n+2)F_n),$$
where $F_n$ is the $n$-th Fibonacci number.
\end{corollary}
\begin{proof}
From Theorem~\ref{fibopoly} the total area of the Fibonacci polyominoes is given by
$$\left.\frac{\partial t_n(1,q)}{\partial q}\right|_{q=1} =\sum_{i=0}^{\lfloor \frac{n+1}{2} \rfloor} \binom{n+1-i}{i}(n + i).$$
For the expression in terms of Fibonacci numbers we use the generating function of the total area
$$A_2(x)=\frac{x (3 + 2 x + x^2)}{(1 -  x -x^2)^2}=(3x+2x^2+x^3)\sum_{n\geq 0}c(n+2)x^n,$$
where $c(n)$ is the convolution of the Fibonacci numbers, that is
$$c(n)=\sum_{i=0}^nF_iF_{n-i}=\frac{1}{5}\left((n-1)F_n + 2nF_{n-1}\right).$$
In the last equality we use the identity (32.13) given in~\cite{Koshy}. Therefore,
$$[x^n]A_2(x)=3c(n+1)+2c(n)+c(n-1)=\frac{1}{5}(6nF_{n+2}+(n+2)F_n), \quad  n \geq 1.$$
\end{proof}

The number of Fibonacci polyominoes of area $n$ is related to the
Narayana's cows sequence ~\cite{Allou}. The
Narayana's cows sequence $b_n$ is defined by the recurrence relation  $b_{n}=b_{n-1}+b_{n-3}$ for $n\geq3$, with the initial values  $b_0=0, b_1=1$, and  $b_2=1$. The Narayana's cows sequence can be calculated with the formula
$$b_n=\sum_{i=0}^{\lfloor (n-1)/3 \rfloor} \binom{n-2i-1}{i}.$$

\begin{theorem}
The number of Fibonacci polyominoes of area $n$ is equal to
 the number of Narayana's cows $b_{n+1}$.
\end{theorem}
\begin{proof}
From Theorem~\ref{teoareaper1} the generating function of the number of Fibonacci polyominoes of a fixed area is
$$F_2(1,1,q)=\frac{q (1 + q + q^2)}{1-q-q^3}.$$
On the other hand, the generating function of the Narayana sequence is  $N(x):=\sum_{i\geq 0}b_ix^i=1/(1-x-x^3)$. Then
$N(q)=1+q+qF_2(1,1,q)$. By comparing
the $n$-th coefficient of the generating functions $N(q)$ and $F_2(1,1,q)$ we obtain  the desired result.
\end{proof}

For example,  Figure~\ref{grafo5} shows the Fibonacci polyominoes   of area 5, that is, $b_{6}=6$ polyominoes.

\begin{figure}[ht]
\centering
\includegraphics[scale=1]{Grafo5.eps}
\caption{Fibonacci polyominoes of area 5.}
\label{grafo5}
\end{figure}

\section{Number of Vertices and Edges}

The goal of this section is to enumerate the number of vertices and edges of the $k$-bonacci graphs.
Let $G$ be a $k$-bonacci graph.  We denote  by ${\textsf{ver}}(G)$ and ${\textsf{edg}}(G)$ the number of vertices and edges  of the graph $G$.  Let ${\mathcal G}_{n,k}^{(1)}$ and  ${\mathcal G}_{n,k}^{(2)}$ denote the set of $k$-bonacci graphs associated with the polyominoes in  ${\mathcal P}_{n,k}^{(1)}$ and  ${\mathcal P}_{n,k}^{(2)}$, respectively.
   Define the generating function
$$G_k(x,p,q):=\sum_{n\geq 1} x^n \sum_{G\in {\mathcal G}_{n,k}}p^{{\textsf{edg}}(G)}q^{{\textsf{ver}}(G)}, \quad k\geq 2.$$
Similarly, we have the generating functions
$$G_{k,j}(x,p,q):=\sum_{n\geq 1} x^n \sum_{G\in {\mathcal G}_{n,k}^{(j)}}p^{{\textsf{edg}}(G)}q^{{\textsf{ver}}(G)},  \quad \text{for} \quad j=1, 2.$$
In these trivariate generating functions
the variable $x$ marks the length of the corresponding $k$-bonacci word, i.e., the number of vertices in the bottom row of a graph minus one.

\begin{theorem}
The generating function for $k$-bonacci graphs with respect to
the length of the corresponding word,
the number of edges and vertices (marked respectively by $x$, $p$ and $q$)
is
$$G_k(x,p,q)=\frac{p^2 q^3((p^2 q + p^5 q^3) x - (p^7 q^4 - p^8 q^5) x^2 - p^{5k} q^{3 k} x^k -
   p^{5 k+3} q^{3 k+2} x^{k+1})}{1 - (p^3 q^2 + p^5 q^3) x + (p^8 q^5 - p^9 q^6) x^2 +
 p^{5 k+4} q^{3 k+3} x^{k+1}}.$$

\end{theorem}
\begin{proof}
From the decomposition given in Figure~\ref{deco} we have the functional equations
\begin{align*}
G_{k,1}(x,p,q)&=\underbrace{p^4q^4x}_{(1)} + \underbrace{p^3q^2x(G_{k,1}(x,p,q)+G_{k,2}(x,p,q))}_{(2)+(3)} \\
G_{k,2}(x,p,q)&=\underbrace{\sum_{j=1}^{k-1}p^{3 j + 2 j + 2}q^{3 (j + 1)}x^j}_{(4)} +\underbrace{\left(\sum_{j=1}^{k-1}p^{3 j + 2 j + 1}q^{3 j + 1}x^j\right) G_{k,1}(x,p,q)}_{(5)}\\
&=\frac{p^2 q^3 (p^5 q^3 x - p^{5 k}q^{3 k} x^k)}{1 - p^5 q^3 x} + \frac{p q (p^5 q^3 x - p^{5 k}q^{3 k} x^k)}{1 - p^5 q^3 x}G_{k,1}(x,p,q).
\end{align*}
Solving the above system of equation and from the equality $G_k(x,p,q)=G_{k,1}(x,p,q)+G_{k,2}(x,p,q)$ we obtain the desired result.
\end{proof}
For example, the series expansion of the generating function $G_{3}(x,p,q)$ begins with
\begin{multline*}
x \left(p^7 q^6+p^4 q^4\right)+x^2 \left(p^{12} q^9+2 p^{10}
   q^8+p^7 q^6\right)
   + \bm{x^3 \left(p^{16} q^{12}+2 p^{15} q^{11}+3
   p^{13} q^{10}+p^{10} q^8\right)}\\
   +x^4 \left(2 p^{21} q^{15}+3 p^{19} q^{14}+3 p^{18} q^{13}+4 p^{16} q^{12}+p^{13}
   q^{10}\right)+ \cdots
    \end{multline*}
Figure~\ref{grafo3} shows the weights of the  $3$-bonacci graphs  corresponding to the bold coefficient in the above series.

\begin{figure}[ht]
\centering
\includegraphics[scale=1]{Grafo3.eps}
\caption{Weights for the graphs in ${\mathcal G}_{3,3}$.}
\label{grafo3}
\end{figure}

\begin{corollary}
The generating functions for the total number of vertices and edges  within all of the members in ${\mathcal G}_{n, k}$ are
$$V_k(x)=\frac{10 x - 9 x^2 - 3 (1 + k) x^k - (4 - 3 k) x^{k+1} + 8 x^{k+2} -  2 x^{2k+2}}{(1 - 2 x + x^{k+1})^2}$$
and
$$E_k(x)=\frac{11 x - 5 x^2 - (2 + 5 k) x^k - (9 - 5 k) x^{k+1} + 4 x^{k+2} + 2 x^{2 k+1} - x^{2 k+2}}{(1 - 2 x + x^{k+1})^2},$$
respectively.
\end{corollary}

In particular, for $k = 2, 3$, and $4$ we obtain the following generating functions:
\begin{align*}
V_2(x)&=\frac{2 x (5 + x - 2 x^2 - x^3)}{(1 -  x -x^2)^2},\\
 V_3(x)&=\frac{x (10 + 11 x - 6 x^3 - 4 x^4 - 2 x^5)}{(1 -  x -x^2-x^3)^2},\\
 V_4(x)&=\frac{x (10 + 11 x + 12 x^2 - 2 x^3 - 8 x^4 - 6 x^5 - 4 x^6 - 2 x^7)}{(1 -  x -x^2-x^3-x^4)^2}, \\
E_2(x)&=\frac{x (11 + 5 x - x^3)}{(1 -  x -x^2)^2},\\
 E_3(x)&=\frac{x (11 + 17 x + 6 x^2 + x^3 - x^5)}{(1 -  x -x^2-x^3)^2},\\
 E_4(x)&=\frac{x (11 + 17 x + 23 x^2 + 7 x^3 + 2 x^4 + x^5 - x^7)}{(1 -  x -x^2-x^3-x^4)^2}.
   \end{align*}

Let $v_n(p,q)$ denote the $n$-th coefficient of  $G_2(x,p,q)$, that is,
$$v_n(p,q):= \sum_{G\in {\mathcal G}_{n,2}}p^{{\textsf{edg}}(G)}q^{{\textsf{ver}}(G)}.$$

\begin{theorem}\label{fibopoly2}
For all $n\geq 3$ we have
$$v_n(p,q)=p^3q^2v_{n-1}(p,q) + p^9q^6v_{n-2}(p,q),$$
with the initial values $v_1(p,q)=p^4 q^4 + p^7 q^6$ and $v_2(p,q)=p^7 q^6 + 2 p^10 q^8$. Moreover,  for all $n\geq 1$ we have the combinatorial formula
\begin{align*}
v_n(p,q)=\sum_{i=0}^{\lfloor \frac{n+1}{2} \rfloor} \binom{n+1-i}{i}p^{3 n + 1 + 3 i} q^{2 n + 2 + 2 i}.
\end{align*}
\end{theorem}

\section{Degree Sequences}

In this section we are interested in the degree sequence of a $k$-bonacci graph. Degree sequences have been well studied for Fibonacci cubes~\cite{Klav5} and Fibonacci-run graphs~\cite{OME3}.
Recall that the degree of a vertex of a graph is the number of edges that are incident to the vertex.  A $k$-bonacci graph can have only vertices of degree two, three, or four.  Let $G$ be a $k$-bonacci graph. We denote  by ${\texttt{deg}}_i(G)$ the number of vertices of degree $i$ in the graph  $G$.

We are interested in the  generating function
$$D_k(x, q_2,q_3,q_4):=\sum_{n\geq1}x^n\sum_{G\in {\mathcal G}_{n,k}}q_2^{{\texttt{deg}}_2(G)}q_3^{{\texttt{deg}}_3(G)}q_4^{{\texttt{deg}}_4(G)}, \quad k\geq 2,$$
 where $x$ marks the length of the corresponding $k$-bonacci word, i.e., the number of vertices in the bottom row of a graph minus one.
Analogously, we introduce  the generating functions
$$D_{k,j}(x, q_2,q_3,q_4):=\sum_{n\geq 1} x^n \sum_{G\in {\mathcal G}_{n,k}^{(j)}}q_2^{{\texttt{deg}}_2(G)}q_3^{{\texttt{deg}}_3(G)}q_4^{{\texttt{deg}}_4(G)}, \quad j=1,  2.$$

\begin{theorem}\label{degreegfg}
For all $k\geq 2$, the generating function $D_k(x,q_2,q_3,q_4)$ is given by
$$\frac{q_2^4 \left( (q_3^2 q_4 +q_4) x -(q_3^2
   q_4^2 - 2 q_2 q_3^2 q_4^2 +q_3^4 q_4) x^2 - q_3^{2 k} q_4^k x^k + (q_3^{2 k+2} q_4^k -2 q_2 q_3^{2 k}
   q_4^{k+1}) x^{k+1}\right)}{q_4(1 -(q_3^2 + q_4 q_3^2) x + (q_4 q_3^4 - q_2^2 q_4^2 q_3^2) x^2 + q_2^2 q_4^{k+1} q_3^{2 k}
   x^{k+1})}.$$
\end{theorem}
\begin{proof}
Let $G$ be a Fibonacci graph in ${\mathcal G}_{n,k}$.  If $n=1$, then $G$ contributes to the generating function the term $q_2^4x$. See Figure~\ref{deco2} case (1).  If $n>1$ and $G\in {\mathcal G}_{n,k}^{(1)}$, then this case contributes to the generating function the terms $$q_3^2xD_{k,1}(x,q_2,q_3,q_4) \quad \text{and} \quad q_2q_4xD_{k,2}(x,q_2,q_3,q_4),$$  as seen in Figure~\ref{deco2} cases (2) and (3). Notice that in the case (2) we have only two new vertices of degree 3.  In the case (3) we have a new vertex of degree 2 and another of degree 4. The red vertices denote the vertices of degree 2, the blue vertices denote the vertices of degree 3, and the green vertices denote the vertices of degree 4. For $n>1$ and $G\in {\mathcal G}_{n,k}^{(2)}$, this case contributes the terms   (see Figure~\ref{deco2} cases (4) and (5))
\begin{align*}
\underbrace{q_2^4\sum_{j=1}^{k-1}q_3^{2 (j - 1) + 2}q_4^{j - 1}x^j}_{(4)} +\underbrace{\left(q_2\sum_{j=1}^{k-1}q_3^{2 (j - 1) + 2}q_4^{j}x^j\right) D_{k,1}(x,q_2,q_3,q_4)}_{(5)}.
\end{align*}

\begin{figure}[ht]
\centering
\includegraphics[scale=1]{deco2.eps}
\caption{Decomposition of a $k$-bonacci graph.}
\label{deco2}
\end{figure}
Therefore, we have the functional equations
\begin{align*}
D_{k,1}(x,q_2,q_3,q_4)&=q_2^4x+ q_3^2xD_{k,1}(x,q_2,q_3,q_4) +  q_2q_4xD_{k,2}(x,q_2,q_3,q_4)\\
D_{k,2}(x,q_2,q_3,q_4)&=\frac{q_2^4 (q_3^2q_4 x - q_3^{2k}q_4^k x^k)}{q_4(1-q_3^2q_4x)} +\frac{q_2 (q_3^2q_4 x - q_3^{2k}q_4^k x^k)}{1-q_3^2q_4x}  D_{k,1}(x,q_2,q_3,q_4).
\end{align*}

Solving the above system of equation we obtain the desired result.
\end{proof}
For example, the series expansion of the generating function $D_3(x,q_2,q_3,q_4)$ begins with
\begin{multline*}
\left(q_2^4q_3^2+q_2^4\right) x+\left(2 q_2^5 q_3^2 q_4 + q_2^4 q_3^2
   + q_2^4 q_3^4 q_4 \right) x^2+\\
   \bm{\left(q_2^6 q_3^4 q_4^2 + q_2^6 q_3^2
   q_4^2 +2 q_2^5 q_3^4 q_4^2 + 2 q_2^5 q_3^4 q_4  + q_3^4 q_2^4\right) x^3} + \cdots
    \end{multline*}
Figure~\ref{grafo4} shows the weights of the  $3$-bonacci graphs  corresponding to the bold coefficient in the above series.
\begin{figure}[ht]
\centering
\includegraphics[scale=1]{Grafo4.eps}
\caption{Weights for the graphs in ${\mathcal G}_{3,3}$.}
\label{grafo4}
\end{figure}

\begin{corollary}
The generating function for the total number of vertices of degree 2 within all of the members in ${\mathcal G}_{n, k}$ is
$$D_{k}^{(2)}(x):=\frac{2 (4 x - 7 x^2 - 2 x^k + x^{k+1} + 7 x^{k+2} - x^{2 k+1} - 2 x^{2 k+2})}{(1 - 2 x + x^{k+1})^2}.$$
\end{corollary}

In particular, for $k = 2, 3$, and $4$ we obtain the following generating functions:
\begin{align*}
D_{2}^{(2)}(x)&=\frac{2 x (4 - x - 5 x^2 - 2 x^3)}{(1 -  x -x^2)^2},\\
D_{3}^{(2)}(x)&=\frac{2 x (4 + x - 4 x^2 - 8 x^3 - 5 x^4 - 2 x^5)}{(1 -  x -x^2-x^3)^2},\\
D_{4}^{(2)}(x)&=\frac{2 x (4 - 7 x - 2 x^3 + x^4 + 7 x^5 - x^8 - 2 x^9)}{(1 - 2 x + x^5)^2}.
   \end{align*}
\begin{corollary}
The generating function for the total number of vertices of degree 3 within all of the members in ${\mathcal G}_{n, k}$ is
$$D_{k}^{(3)}(x):=\frac{2 (x + x^2 - k x^k - (1 - k) x^{k+1} - 2 x^{k+2} + x^{2k+2})}{(1 - 2 x + x^{k+1})^2}.$$
\end{corollary}
In particular, for $k = 2, 3$ and $4$ we obtain the following generating functions:
\begin{align*}
D_{2}^{(3)}(x)&=\frac{2 x (1 + x + 2 x^2 + x^3)}{(1 -  x -x^2)^2},\\
D_{3}^{(3)}(x)&=\frac{2 (x + x^2 - 3 x^3 + 2 x^4 - 2 x^5 + x^8)}{(1 - 2 x + x^4 )^2},\\
D_{4}^{(3)}(x)&=\frac{2(x + x^2 - 4 x^4 + 3 x^5 - 2 x^6 + x^{10})}{(1 - 2 x + x^5)^2}.
   \end{align*}

  \begin{corollary}
The generating function for the total number of vertices of degree 4 within all of the members in ${\mathcal G}_{n, k}$ is
$$D_{k}^{(4)}(x):=\frac{3 x^2 + (1-k) x^k - (4 - k) x^{k+1} - 2 x^{k+2} + 2x^{2k+1}}{(1 - 2 x + 2 x^{k+1})^2}.$$
\end{corollary}

In particular, for $k = 2, 3$ and $4$ we obtain the following generating functions:
\begin{align*}
D_{2}^{(4)}(x)&=\frac{2 x^2 (1 + x)}{(1 -  x -x^2)^2},\\
D_{3}^{(4)}(x)&=\frac{x^2 (3 + 4 x + 4 x^2 + 2 x^3)}{(1 - x - x^2 - x^3 )^2},\\
D_{4}^{(4)}(x)&=\frac{3 x^2 - 3 x^4 - 2 x^6 + 2 x^9}{(1 - 2 x + x^5)^2}.
\end{align*}

\subsection{The degree sequence of the Fibonacci graph}
Finally, we consider the particular case of Fibonacci graphs ($k=2$).  Let $d_{n,2}(q)$ denote the $n$-th coefficient of  $D_2(x, q,1,1)$, that is,
$$d_{n,2}(q):= \sum_{G\in {\mathcal G}_{n,2}}q^{{\textsf{deg}}_2(G)}.$$
Similarly, we define the sequences
$$d_{n,3}(q):= \sum_{G\in {\mathcal G}_{n,2}}q^{{\textsf{deg}}_3(G)} \quad \text{and} \quad d_{n,4}(q):= \sum_{G\in {\mathcal G}_{n,2}}q^{{\textsf{deg}}_4(G)}.$$
\begin{theorem}\label{closeformula}
For all $n\geq 3$ we have
$$d_{n,2}(q)=d_{n-1,2}(q) + q^2d_{n-2,2}(q),$$
with the initial values $d_{1,2}(q)=2q^4$ and $d_{2,2}(q)=q^4 + 2 q^4$. Moreover,  for all $n\geq 1$ we have the combinatorial formula
\begin{align}\label{vertdeg2}
d_{n,2}(q)=\sum_{i=1}^{n} \left(\binom{n-1-\lfloor \frac{i}{2}\rfloor}{\lfloor\frac{i-1}{2}\rfloor}+  \binom{n-2-\lfloor \frac{i-1}{2}\rfloor}{\lfloor\frac{i-2}{2}\rfloor}\right)q^{i+3}.
\end{align}
\end{theorem}
\begin{proof}
Define the generating function
$D_2(x,q):=\sum_{n\geq 0} d_{n,2}(q)x^n$.
From the Theorem~\ref{degreegfg} we have the expression
$$D_2(x,q)=\frac{xq^4(2-x+2xq)}{1-x-q^2x^2}.$$
This generating function satisfies
$$D_2(x,q)-xD_2(x,q)-q^2x^2D_2(x,q)=2 q^4 x - q^4 x^2 + 2 q^5 x^2.$$
By comparing the  coefficient  of $x^n$, for  $n\geq 2$,  in the above equality, we obtain the recurrence relation for the sequence $d_{n,2}(q)$. On the other hand, the combinatorial sum \eqref{vertdeg2} is equivalent to
\begin{align}\label{closedf2}
\sum_{i\geq 0} \left(2\binom{n-i-1}{i-1}q^{2i+3}+ \binom{n-i}{i-1}q^{2i+2} + \binom{n-i-1}{i-2}q^{2i+2} \right).
\end{align}
Denote the above summand  by $F(n,i)$, that is
\begin{align*}
F(n,i):=& 2\binom{n-i-1}{i-1}q^{2i+3}+ \binom{n-i}{i-1}q^{2i+2} + \binom{n-i-1}{i-2}q^{2i+2}\\
=&\frac{(n-i-1)!\left(n(2q +1) -4 i q  + 2q-1\right)}{(i-1)!(n-2i+1)!}q^{2i+2}.
\end{align*}
By the Zeilberger algorithm, $F(n,i)$ satisfies the relation
\begin{align}\label{rel2}
F(n+2,i)-F(n+1,i)-q^2F(n,i)=G(n,i+1)-G(n,i),
\end{align}
with the certificate
$$R(n,i)=\frac{(n-i) (i-1)\left((4 i-6)q - n(2q+1)  + 1\right)}{(n+2-2 i) (n+3-2 i) \left(2 n q+n-4 i q+2 q-1\right)}.$$
If $f(n)$ denotes the combinatorial  sum in \eqref{closedf2}, then summing both sides of \eqref{rel2} with respect to $i$ yields $f(n+2)-f(n+1)-q^2f(n)=0$.  The sequences $f(n)$ and $d_{n,2}(q)$ satisfy the same recurrence relation and  have the same initial values, therefore these sequences coincides for all positive integers $n$.

\end{proof}

Notice that the  total vertices of degree 2 of all Fibonacci graphs in ${\mathcal G}_{n,2}$  is given by
$$\sum_{i=1}^{n}(i+3) \left(\binom{n-1-\lfloor \frac{i}{2}\rfloor}{\lfloor\frac{i-1}{2}\rfloor}+  \binom{n-2-\lfloor \frac{i-1}{2}\rfloor}{\lfloor\frac{i-2}{2}\rfloor}\right).$$

From a similar argument as in Theorem~\ref{closeformula} we can prove the following two theorems.

\begin{theorem}
For all $n\geq 3$ we have
$$d_{n,3}(q)=q^2(d_{n-1,3}(q) + d_{n-2,3}(q)),$$
with the initial values $d_{1,3}(q)=q^2+1$ and $d_{2,3}(q)=3q^2$. Moreover,  for all $n\geq 1$ we have the combinatorial formula
$$d_{n,3}(q)=\sum_{i=0}^{n}\left((1+q^2)q^{2(n-1-i)}\binom{n-i-1}{i} + (2q^2-q^4)q^{2(n-1-i)}\binom{n-i-2}{i}\right).$$
\end{theorem}

\begin{theorem}
For all $n\geq 3$ we have
$$d_{n,4}(q)=d_{n-1,4}(q) + q^2d_{n-2,4}(q),$$
with the initial values $d_{1,2}(q)=2q^4$ and $d_{2,2}(q)=q^4 + 2 q^4$. Moreover,  for all $n\geq 1$ we have the combinatorial formula
$$d_{n,4}(q)=\sum_{i=0}^{n+2} \left(\binom{n-1-\lfloor \frac{i-2}{2}\rfloor}{\lfloor\frac{i-3}{2}\rfloor}+  \binom{n-2-\lfloor \frac{i-3}{2}\rfloor}{\lfloor\frac{i-4}{2}\rfloor}\right)q^{i-3}.$$
\end{theorem}

Let $d_n$ denote the total number of vertices of the Fibonacci polyominoes in ${\mathcal G}_{n,2}$. In the following theorems we study the proportion between the sequences $d_{n,i}:=d_{n,i}(1)$ ($i=2, 3, 4$) and $d_n$.   Before,  we need the following result.

\begin{theorem}[Asymptotics of linear recurrences,~\cite{Flajolet2}]\label{FlaTheo}
Assume that a rational generating function $f(x)/g(x)$, with $f(x)$ and $g(x)$ relatively prime and $g(0)\neq 0$, has a unique pole $1/\beta$ of smallest modulus. Then, if the multiplicity of $1/\beta$ is $\nu$, we have
$$[x^n]\frac{f(x)}{g(x)}\sim \nu \frac{(-\beta)^\nu f(1/\beta)}{g^{(\nu)}(1/\beta)}\beta^n n^{\nu-1}.$$

\end{theorem}

\begin{theorem}\label{asymp1}
Among total degree of  vertices of all graphs in ${\mathcal G}_{n,2}$, the proportion of those that are of degree 2 is asymptotically
\[
\lim_{n\to\infty}\frac{d_{n,2}}{d_n} = \frac{7-\sqrt{5}}{22}\approx 0.21654236.
\]
\end{theorem}
\begin{proof}
The generating functions of the sequences $d_n$ and $d_{n,2}$ are rational, therefore we can use the asymptotic analysis for linear recurrences. First, note that the unique pole $1/\beta$ of the rational generating function
$$V_2(x)=\sum_{n \geq 0}d_nx^n=\frac{2 x (5 + x - 2 x^2 - x^3)}{(1 - x - x^2)^2}$$ is $\alpha=\frac{-1+\sqrt{5}}{2}$, with multiplicity 2. Therefore
$d_n\sim \frac{2(3+2\sqrt{5})}{5}\left(\frac{1+\sqrt{5}}{2}\right)^nn$. Similarly, we have $d_{n,2}\sim(\frac{1+\sqrt{5}}{2})^{n+1}n$. From these expressions we obtain the desired result.
\end{proof}

\begin{theorem}
Among total degree of  vertices of all graphs in ${\mathcal G}_{n,2}$, the proportion of those that are of degree 3 is asymptotically
\[
\lim_{n\to\infty}\frac{d_{n,3}}{d_n} = \frac{4+\sqrt{5}}{11}\approx 0.56691527.
\]
\end{theorem}
\begin{proof}
The unique pole $1/\beta$ of the rational generating function
$$\sum_{n \geq 0}d_{n,3}x^n=\frac{2 (x - x^2 + x^3 - 2 x^4 + x^6)}{(1 - 2 x + x^3)^2}$$ is $\alpha=\frac{-1+\sqrt{5}}{2}$, with multiplicity 2. Therefore
$d_{n,2}\sim \frac{2(2+\sqrt{5})}{5}(\frac{1+\sqrt{5}}{2})^nn$. Since \linebreak $d_n\sim \frac{2(3+2\sqrt{5})}{5}(\frac{1+\sqrt{5}}{2})^nn$  we obtain the desired result.
\end{proof}

\begin{theorem}
Among total degree of  vertices of all graphs in ${\mathcal G}_{n,2}$, the proportion of those that are of degree 4 is asymptotically
\[
\lim_{n\to\infty}\frac{d_{n,4}}{d_n} = \frac{7-\sqrt{5}}{22}\approx 0.216542364.
\]
\end{theorem}

\section{Number of Hamiltonian $k$-bonacci graphs}
A \emph{Hamiltonian cycle} is a cycle that visits each vertex exactly once.  Let $G$ be a $k$-bonacci graph.  We define  ${\texttt{Ham}}(G)=1$ if $G$ has a Hamiltonian cycle, and 0 otherwise.  If ${\texttt{Ham}}(G)=1$, we say that $G$ is a Hamiltonian $k$-bonacci graph.
   Define the generating function
$$H_k(x,q):=\sum_{n\geq 1} x^n \sum_{G\in {\mathcal G}_{n,k}}q^{{\texttt{Ham}}(G)}, \quad k\geq 2,$$
 where $x$ marks the length of the corresponding $k$-bonacci word, i.e., the number of vertices in the bottom row of a graph minus one.
Similarly, we have the generating functions
$$H_{k,j}(x,q):=\sum_{n\geq 1} x^n \sum_{G\in {\mathcal G}_{n,k}^{(j)}}q^{{\texttt{Ham}}(G)},  \quad \text{for} \quad j=1, 2.$$

\begin{theorem}
For all $k\geq 2$ we have
$$H_k(x,q)=\frac{x ((1 - x) x (1 - x^{2 \lfloor (k-1)/2\rfloor}) -
   q (1 + x) (1 - 2 x + x^{k+1}) (-2 + x + x^{2 \lfloor k/2\rfloor}))}{(1 - 2 x + x^{k+1}) (1 - x - 2 x^2 + x^3 + x^{2(1 +  \lfloor k/2 \rfloor})}.$$
\end{theorem}
\begin{proof}
From the decomposition given in Figure~\ref{deco} we have the functional equation
\begin{align*}
H_{k,1}(x,q)&=\underbrace{qx}_{(1)} + \underbrace{x(H_{k,1}(x,q)+H_{k,2}(x,q))}_{(2)+(3)}
\end{align*}
The polyomino given in the decomposition $(4)$ corresponds to the grid graph $P_3 \times P_i$, for $1 < i \leq k$. It is known\footnote{~See, for example the discussion on ``Math Stackexchange": \url{https://math.stackexchange.com/questions/1699203/hamilton-paths-cycles-in-grid-graphs}} that the grid graph $P_n \times P_m$ has a hamiltonian cycle if and only if at least one of $m$ or $n$ is even
or $m=n=1$. Therefore, the graph $P_3 \times P_i$ has a hamiltonian cycle if and only if $i$ is even. In this case, the generating function is given by
$$T_k(x,q)=\underbrace{q\sum_{j=1}^{\lfloor k/2 \rfloor}x^{2j-1} + \sum_{j=1}^{\lfloor (k-1)/2 \rfloor}x^{2j}}_{(4)} =\frac{x\left(q + x - x^{2\lfloor\frac{k-1}{2} \rfloor+1} -  qx^{2\lfloor\frac{k}{2} \rfloor} \right)}{1-x^2}.
$$
Therefore, we obtain the functional equation
\begin{align*}
H_{k,2}(x,q)&=\underbrace{T_k(x,q)}_{(4)} +\underbrace{\left(\sum_{j=1}^{\lfloor k/2 \rfloor}x^{2j-1}\right) H_{k,1}(x,q) +\left(\sum_{j=1}^{\lfloor(k-1)/2 \rfloor}x^{2j}\right) H_{k,1}(x,1)}_{(5)}\\
&=T_k(x,q) + \frac{x(1-x^{2\lfloor k/2 \rfloor})}{1-x^2}H_{k,1}(x,q) + \frac{x^2(1-x^{2\lfloor (k-1)/2 \rfloor})}{1-x^2}\frac{x(1-x^k)}{1-2x+x^{k+1}}.
\end{align*}
Notice that $H_{k,1}(x,1)$ is the generating function of the total number of $k$-bonacci words end in \texttt{1}.
Solving the above system of equation  we obtain the desired result.
\end{proof}

\begin{corollary}
The generating function for the total number of Hamiltonian $k$-bonacci graphs is
$$H_{k}(x):=\frac{x(1+x)(2-x-x^{2\lfloor k/2 \rfloor})}{1 - x - 2 x^2 + x^3 + x^{2 \lfloor k/2\rfloor +2}}.$$
\end{corollary}

In particular, for $2 \le k \le 7$
we obtain the following generating functions:
\begin{align*}
H_{2}(x)&=H_{3}(x)=\frac{x (2 + x)}{1 -  x -x^2},\\
H_{4}(x)&=H_{5}(x)=\frac{x (2 + x + x^2 + x^3)}{1 - x - x^2 - x^4},\\
H_{6}(x)&=H_{7}(x)=\frac{x (2 + x + x^2 + x^3 + x^4 + x^5)}{1 - x - x^2 - x^4 - x^5}.
   \end{align*}

Note that every Fibonacci graph admits a Hamiltonian cycle,
just walk along its border.
The grid $P_3 \times P_3$ does not have a Hamiltonian cycle, so
Hamiltonian 3-bonacci graphs cannot have grids $P_3 \times P_3$ as induced subgraphs. Thus, they are precisely 2-bonacci.
In general, Hamiltonian $2k$-bonacci equals
Hamiltonian $(2k+1)$-bonacci graphs, because any Hamiltonian $(2k+1)$-bonacci graph cannot
contain $P_{2k+1} \times P_3$ as
induced subgraph, and there are no $1^{2k}$ factors in the corresponding binary words.