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

Athlete's career trajectory is a popular topic of discussion in the
media nowadays. Questions regarding whether a player has reached their
peak, is past their prime, or is good enough to remain in their
respective professional league are often seen in media outlets such as
news articles, television debate shows, and podcasts. The average
performance of players by age throughout their careers is visually
represented by an \emph{aging curve}. This graph typically consists of a
horizontal axis representing a time variable (usually age or season) and
a vertical axis showing a performance metric at each time point in a
player's career.

One significant challenge associated with the study of aging curves in
sports is \emph{survival bias}, as pointed out by Lichtman (2009),
Turtoro (2019), Judge (2020a), and Schuckers et al. (2021). In
particular, the aging effects are not often determined from a full
population of athletes in a given league. That is, only players that are
good enough to remain are observed; whereas those might be involved, but
do not actually participate or not talented enough to compete, are being
completely disregarded. This very likely results in an overestimation in
the construction of aging curves.

As such, player survivorship and dropout can be viewed in the context of
missing data. There are different cases of player absence from
professional sport at different points in their careers. At the
beginning, teams may elect to assign their young prospects to their
minor/development league affiliates for several years of nurture. Many
of those players would end up receiving a call-up to join the senior
squad, when the team believes they are ready. During the middle of one's
career, a nonappearance could occur due to various reasons. Injury is
unavoidable in sports, and this could cost a player at least one year of
their playing time. Personal reasons such as contract situation and more
recently, concerns regarding a global pandemic, could also lead to
athletes sitting out a season. Later on, a player might head for
retirement because they cannot perform at a level like they used to.

The primary aim of this paper is to apply missing data techniques to the
estimation of aging curves. In doing so, we focus on baseball and pose
one research question: What would the aging curve look like if players
competed in every season from a fixed range of age? In other words, what
would have happened if a player who was forced to retire from their
league at a certain age had played a full career? The manuscript
continues with a review of existing literature on aging curves in
baseball and other sports in Section \ref{sec:lit}. Next, we describe
our data and methods used to perform our analyses in Section
\ref{sec:meth}. After that, our approach is implemented through
simulation and analyses of real baseball data in Sections \ref{sec:sim}
and \ref{sec:app}. Finally, in Section \ref{sec:discuss}, we conclude
with a discussion of the results, limitations, and directions for future
work.

\hypertarget{sec:lit}{%
\section{Literature Review}\label{sec:lit}}

To date, we find a considerable amount of previous work related to aging
curves and career trajectory of athletes. This body of work consists of
several key themes, a wide array of statistical methods, and
applications in many sports besides baseball such as basketball, hockey,
and track and field.

A typical notion in the baseball aging curves literature is the
assumption of a quadratic form for modeling the relationship between
performance and age. Morris (1983) looks at Ty Cobb's batting average
trajectory using parametric empirical Bayes and uses shrinkage methods
to obtain a parabolic curve for Cobb's career performance. Albert (1992)
proposes a quadratic random effects log-linear model for smoothing a
batter's home run rates throughout their career. A nonparametric method
is implemented to estimate the age effect on performance in baseball,
hockey, and golf by Berry et al. (1999). However, Albert (1999) weighs
in on this nonparametric approach and questions the assumptions that the
peak age and periods of growth and decline are the same for all players.
Albert (1999) ultimately prefers a second-degree polynomial function for
estimating age effect in baseball, which is a parametric model.
Continuing his series of work on aging trajectories, Albert (2002)
proposes a Bayesian exchangeable model for modeling hitting performance.
This approach combines quadratic regression estimates and assumes
similar careers for players born in the same decade. Fair (2008) and
Bradbury (2009) both use a fixed-effects regression to examine age
effects in the MLB, also assuming a quadratic aging curve form.

In addition to baseball, studies on aging curves have also been
conducted for other sports. Early on, Moore (1975) looks at the
association between age and running speed in track and field and
produces aging curves for different running distances using an
exponential model. Fair (1994) and Fair (2007) study the age effects in
track and field, swimming, chess, and running, in addition to their
latter work in baseball, as mentioned earlier. In triathlon, Villaroel
et al. (2011) assume a quadratic relationship between performance and
age, as many have previously considered. As for basketball, Page et al.
(2013) use a Gaussian process regression in a hierarchical Bayesian
framework to model age effect in the NBA. Additionally, Lailvaux et al.
(2014) use NBA and WNBA data to investigate and test for potential sex
differences in the aging of comparable performance indicators. Vaci et
al. (2019) apply Bayesian cognitive latent variable modeling to explore
aging and career performance in the NBA, accounting for player position
and activity levels. In tennis, Kovalchik (2014) studies age and
performance trends in men's tennis using change point analysis.

Another convention in the aging curve modeling literature is the
assumption of discrete observations. Specifically, most researchers use
regression modeling and consider a data measurement for each season
played throughout a player's career. In contrast to previous approaches,
Wakim \& Jin (2014) take a different route and consider functional data
analysis as the primary tool for modeling MLB and NBA aging curves. This
is a continuous framework which treats the entire career performance of
an athlete as a smooth function.

A subset of the literature on aging and performance in sports provides
answers to the question: At what age do athletes peak? Schulz \& Curnow
(1988) look at the age of peak performance for track and field,
swimming, baseball, tennis, and golf. A follow-up study to this work was
done by Schulz et al. (1994), where the authors focus on baseball and
find that the average peak age for baseball players is between 27 and
30, considering several performance measures. Later findings on baseball
peak age also show consistency with the results in Schulz et al. (1994).
Fair (2008) determines the peak-performance age in baseball to be 28,
whereas Bradbury (2009) determines that baseball hitters and pitchers
reach their career apex at about 29 years old. In soccer, Dendir (2016)
determines that the peak age for footballers in the top leagues falls
within the range of 25 to 27.

Last and most importantly, the idea of player survivorship is only
mentioned in a small number of articles. To our knowledge, not many
researchers have incorporated missing data methods into the estimation
of aging curves to account for missing but observable athletes. Schulz
et al. (1994) and Schell (2005) note the selection bias problem with
estimating performance averages by age in baseball, as better players
tend to have longer career longevity. Schall \& Smith (2000) predict
survival probabilities of baseball players using a logit model, and
examine the link between first-year performance and career length.
Lichtman (2009) studies different aging curves for different eras and
groups of players after correcting for survival bias, and shows that
survival bias results in an overestimation of the age effects.
Alternatively, Judge (2020a) concludes that survivorship bias leads to
an underestimation, not overestimation, of the aging curves. In
analyzing NHL player aging, Brander et al. (2014) apply their quadratic
and cubic fixed-effects regression models to predict performance for
unobserved players in the data.

Perhaps the most closely related approach to our work is that by
Schuckers et al. (2021), which considers different regression and
imputation frameworks for estimating the aging curves in the National
Hockey League (NHL). First, they investigate different regression
approaches including spline, quadratic, quantile, and a delta plus
method, which is an extension to the delta method previously studied by
Lichtman (2009), Turtoro (2019), and Judge (2020b). This paper also
proposes an imputation approach for aging curve estimation, and
ultimately concludes that the estimation becomes stronger when
accounting for unobserved data, which addresses a major shortcoming in
the estimation of aging curves. However, it appears that the aging
curves are constructed without taking into account the variability as a
result of imputing missing data. This could be improved by applying
multiple imputation rather than considering only one imputed dataset. As
pointed out by Gelman \& Hill (2006) (Chapter 25), conducting only a
single imputation essentially assumes that the filled-in values
correctly estimate the true values of the missing observations. Yet,
there is uncertainty associated with the missingness, and multiple
imputation can incorporate the missing data uncertainty and provide
estimates for the different sources of variability.

\hypertarget{sec:meth}{%
\section{Methods}\label{sec:meth}}

\hypertarget{sec:data}{%
\subsection{Data Collection}\label{sec:data}}

In the forthcoming analyses, we rely on one primary source of publicly
available baseball data: the Lahman baseball database (Lahman, 1996 --
2021). Created and maintained by Sean Lahman, this database contains
pitching, hitting, and fielding information for Major League Baseball
players and teams dating back to 1871. The data are available in many
different formats, and the \texttt{Lahman} package in \texttt{R}
(Friendly et al., 2021; R Core Team, 2022) is utilized for our
investigation.

Due to our specific purpose of examining the aging curves for baseball
offensive players, the following datasets from the \texttt{Lahman}
library are considered: \texttt{Batting}, which provides
season-by-season batting statistics for baseball players; and
\texttt{People}, to obtain the date of birth of each player and
calculate their age for each season played. In each table, an athlete is
identified with their own \texttt{playerID}, hence this attribute is
used as a joining key to merge the two tables together. A player's age
for a season is determined as their age on June 30, and the formula
suggested by Marchi et al. (2018) for age adjustment based on one's
birth month is applied.

Throughout this paper, we consider on-base plus slugging (OPS), which
combines a hitter's ability to reach base on a swing and power-hitting,
as the baseball offensive performance measure. We scale the OPS for all
players and then apply an arcsine transformation to ensure a reasonable
range for the OPS values when conducting simulation and imputation. We
also assume a fixed length for a player's career, ranging from age 21 to
39. In terms of sample restriction, we observe all player-seasons with
at least 100 plate appearances, which means a season is determined as
missing if one's plate appearances is below that threshold.

\hypertarget{multiple-imputation}{%
\subsection{Multiple Imputation}\label{multiple-imputation}}

Multiple imputation (Rubin, 1987) is a popular statistical procedure for
addressing the presence of incomplete data. The goal of this approach is
to replace the missing data with plausible values to create multiple
completed datasets. These datasets can each be analyzed and results are
combined across the imputed versions. Multiple imputation consists of
three steps. First, based on an appropriate imputation model, \(m\)
copies of the dataset are created by filling the missing values. Next,
\(m\) analyses are performed on each of the \(m\) completed datasets.
Finally, the results from each of the \(m\) datasets are pooled together
to create a combined estimate and standard errors are estimated that
account for the between and within imputation variability. This last
step can be accomplished using asymptotic theory with Rubin's combining
rules (Little \& Rubin, 1987), which are as follows.

Let \(Q\) be a parameter of interest and \(\widehat Q_i\) where
\(i=1,2,\dots,m\) are estimates of \(Q\) obtained from \(m\) imputed
datasets, with sampling variance \(U\) estimated by \(\widehat U_i\).
Then the point estimate for \(Q\) is the average of the \(m\) estimates
\[
\overline Q = \frac{1}{m} \sum_{i=1}^m \widehat Q_i
\,.
\] The variance for \(\overline Q\) is defined as \[
T=\overline U + \left(1 + \frac{1}{m}\right)B
\,,
\] where \[
\overline U = \frac{1}{m} \sum_{i=1}^m \widehat U_i
\] and \[
B=\frac{1}{m-1} \sum_{i=1}^m (\widehat Q_i - \overline Q)^2
\] are the estimated within and between variances, respectively.

Inferences for \(Q\) are based on the approximation \[
\frac{Q - \overline Q}{\sqrt{T}} \sim t_\nu
\,,
\] where \(t_\nu\) is the Student's \(t\)-distribution with
\(\displaystyle \nu = (m-1)\left(1+\frac{1}{r}\right)^2\) degrees of
freedom, with
\(\displaystyle r=\left(1+\frac{1}{m}\right)\frac{B}{\overline U}\)
representing the relative increase in variance due to missing data.

Accordingly, a \(100(1-\alpha)\)\% Wald confidence interval for \(Q\) is
computed as \[
\overline Q \ \pm \ t_{\nu,1-\alpha/2}\sqrt{T}
\,,
\] where \(t_{\nu,1-\alpha/2}\) is the \(1-\alpha/2\) quantile of
\(t_\nu\).

It is important to understand the reasons behind the missingness when
applying multiple imputation to handle incomplete data. Based on the
degree of bias attributable to the missingness, three types of missing
data are defined: missing completely at random (MCAR), missing at random
(MAR), and missing not at random (MNAR) (Rubin, 1976). MCAR occurs when
a missing observation is statistically independent of both the observed
and unobserved data. In the case of MAR, the missingness is associated
with the observed but not with the unobserved data. When data are MNAR,
there is a link between the missingness and the unobserved values in the
dataset.

Among the tools for performing multiple imputation, multivariate
imputations by chained equation (MICE) (van Buuren \&
Groothuis-Oudshoorn, 1999) is a flexible, robust, and widely used
method. This algorithm imputes missing data via an iterative series of
conditional models. In each iteration, each incomplete variable is
filled in by a separate model of all the remaining variables. The
iterations continue until apparent convergence is reached.

In this paper, we implement the MICE framework in \texttt{R} via the
popular \texttt{mice} package (van Buuren \& Groothuis-Oudshoorn, 2011).
Moreover, we focus on multilevel multiple imputation, due to the
hierarchical structure of our data. Specifically, we consider multiple
imputation by a two-level normal linear mixed model with heterogeneous
within-group variance (Kasim \& Raudenbush, 1998). In context, our data
consist of baseball seasons (ages) which are nested within the class
variable, player; and the season-by-season performance is considered to
be correlated for each athlete. The described imputation model can be
specified as the \texttt{2l.norm} method available in the \texttt{mice}
library.

\hypertarget{sec:sim}{%
\section{Simulation}\label{sec:sim}}

In this simulation, we demonstrate our aging curve estimation approach
with multiple imputation, and evaluate how different types of player
dropouts affect the curve. There are three steps to our simulation.
First, we fit a model for the performance-age relationship and utilize
its output to generate fictional careers for baseball players. Next, we
generate missing data by dropping players from the full dataset based on
different criteria, and examine how the missingness affects the original
aging curve obtained from fully observed data. Finally, we apply
multiple imputation to obtain completed datasets and assess how close
the imputed aging curves are to the true curve.

\hypertarget{generating-player-careers}{%
\subsection{Generating Player Careers}\label{generating-player-careers}}

We fit a mixed-effects model using the player data described in Section
\ref{sec:data}. Our goal is to obtain the variance components of the
fitted model to simulate baseball player careers. The model of
consideration is of the form \[
\displaylines{
Y_{pq} = (\beta_0 + b_{0p}) + \beta_1X_q + \beta_2X_q^2 + \beta_3X_q^3 + \epsilon_{pq}
\cr
b_{0p} \sim N(0, \tau^2)
\cr
\epsilon_{pq} \sim N(0, \sigma^2).
}
\] In detail, this model relates the performance metric \(Y_{pq}\) (in
our case, transformed OPS) for player \(p\) at age (season) \(q\) to a
baseline level via the fixed effect \(\beta_0\). The only covariate
\(X\) in the model is age, which is assumed to have a cubic relationship
with the response variable, transformed OPS. Another component is the
observational-level error \(\epsilon_{pq}\) with variance \(\sigma^2\)
for player \(p\) at age \(q\). We also introduce the random effects
\(b_{0p}\), which represents the deviation from the grand mean
\(\beta_0\) for player \(p\), allowing a fixed amount of shift to the
performance prediction for each player. In addition, to incorporate the
variability in production across the season \(q\), a random effect
parameter \(\tau^2\) is included. Our modeling approach is implemented
using the \texttt{lme4} package in \texttt{R} (Bates et al., 2015). We
utilize the sources of variance from the fitted model to simulate 1000
careers for baseball players of age 21 to 39.

\hypertarget{sec:drop}{%
\subsection{Generating Missing Data}\label{sec:drop}}

After obtaining reasonable simulated careers for baseball players, we
create different types of dropouts and examine how they cause deviations
from the fully observed aging curve. We consider the following cases of
players' retirement from the league:

\begin{enumerate}
\def(\arabic{enumi}){(\arabic{enumi})}
\tightlist
\item
  Dropout players with 4-year OPS average below a certain threshold, say
  0.55.
\item
  Dropout players with OPS average from age 21 (start of career) to 25
  of less than 0.55.
\item
  25\% of the players randomly retire at age 30.
\end{enumerate}

For the first two scenarios, the missingness mechanism is MAR, since
players get removed due to low previously observed performance. Dropout
case (3) falls under MCAR, since athletes are selected at random to
retire without considering past or future offensive production.

Figure \ref{fig:drop-compare} displays the average OPS aging curves for
all baseball players obtained from the original data with no missingness
and data with only the surviving players associated with the dropout
mechanisms mentioned above. These are smoothed curves obtained from
loess model fits, and we use mean absolute error (MAE) to evaluate the
discrepancy between the dropout and true aging curves. It is clear that
randomly removing players have minimal effect on the aging curve, as the
curve obtained from (3) and the original curve essentially overlap (MAE
\(= 7.42 \times 10^{-4}\)). On the other hand, a positive shift from the
fully observed curve occurs for the remaining two cases of dropout based
on OPS average (MAE \(= 0.031\) for (1) and MAE \(=0.019\) for (2)).
This means the aging curves with only the surviving players are
overestimated in the presence of missing data due to past performance.
More specifically, the player performance drops off faster as they age
than when it is estimated with only complete case analysis.

\begin{figure}

{\centering \includegraphics{drop-compare-1}

}

\caption{Comparison of the average OPS aging curve constructed with the fully observed data and different cases of dropouts (obtained only for the surviving players and without imputation).}\label{fig:drop-compare}
\end{figure}

\hypertarget{sec:imp}{%
\subsection{Imputation}\label{sec:imp}}

Next, we implement the multiple imputation with a hierarchical structure
procedure described in Section \ref{sec:meth} to the cases of dropout
that shifts the aging effect on performance. We perform \(m=5\)
imputations with each made up of 30 iterations, and apply Rubin's rules
for combining the imputation estimates. The following results are
illustrated for dropout mechanism (2), where players with a low OPS
average at the start of their careers (ages 21--25) are forced out of
the league.

Figure \ref{fig:drop-imp} (left) shows smoothed fitting loess aging
curves for all 5 imputations and a combined version of them, in addition
to the curves constructed with fully observed and only surviving players
data. The 95\% confidence interval for the mean OPS at each age point in
the combined curve obtained from Rubin's rules is further illustrated in
Figure \ref{fig:drop-imp} (right). It appears that the combined imputed
curve follows the same shape as the true, known curve. Moreover,
imputation seems to capture the rate of change for the beginning and end
career periods quite well, whereas the middle of career looks to be
slightly underestimated. The resulting MAE of \(0.0039\) confirms that
there is little deviation of the combined curve from the true one.

Additionally, we perform diagnostics to assess the plausibility of the
imputations, and also examine whether the algorithm converges. We first
check for distributional discrepancy by comparing the distributions of
the fully observed and imputed data. Figure \ref{fig:diag} (left)
presents the density curves of the OPS values for each imputed dataset
and the fully simulated data. It is obvious that the imputation
distributions are well-matched with the observed data. To confirm
convergence of the MICE algorithm, we inspect trace plots for the mean
and standard deviation of the imputed OPS values. As shown in Figure
\ref{fig:diag} (right), convergence is apparent, since no definite trend
is revealed and the imputation chains are intermingled with one another.

\begin{figure}

{\centering \includegraphics{drop-imp-1}

}

\caption{At left, comparison of the average OPS aging curve constructed with the fully observed data, only surviving players, and imputation. At right, combined imputed curve with 95\% confidence intervals obtained from Rubin's rules. Results shown here are for dropout case of players having OPS average from age 21 to 25 below 0.55.}\label{fig:drop-imp}
\end{figure}

\begin{figure}

{\centering \includegraphics{diag-1} \includegraphics{diag-2}

}

\caption{At left, kernel density estimates for the fully observed and imputed OPS values. At right, trace plots for the mean and standard deviation of the imputed OPS values against the iteration number for the imputed data. Results shown here are for dropout case of players having OPS average from age 21 to 25 below 0.55.}\label{fig:diag}
\end{figure}

\hypertarget{sec:app}{%
\section{Application: MLB Data}\label{sec:app}}

Lastly, we apply the multilevel multiple imputation model to estimate
the average OPS aging curve for MLB players. For this investigation,
besides the data pre-processing tasks mentioned in Section
\ref{sec:data}, our sample is limited to all players who made their
major league debut no sooner than 1985, resulting in a total of 2323
players. To perform imputation, we pass in the parameters similar to our
simulation study (\(m=5\) with 30 iterations for each imputation).

Figure \ref{fig:mlb-imp} shows the OPS aging curves for MLB players
estimated with and without imputation. The plot illustrates a similar
result as from simulation, as the combined imputed curve is lower than
the curve obtained when ignoring the missing data. It is clear that the
aging effect is overestimated without considering the unobserved
player-seasons that are observable. In other words, the actual
performance declines with age more rapidly than estimates based on only
the observed data.

\begin{figure}

{\centering \includegraphics{mlb-imp-1}

}

\caption{Comparison of the average OPS aging curve constructed with only observed players and imputation for MLB data.}\label{fig:mlb-imp}
\end{figure}

\hypertarget{sec:discuss}{%
\section{Discussion}\label{sec:discuss}}

The concept of survivorship bias is frequently seen in professional
sports, and our paper approaches the topic of aging curves and player
dropout in baseball as a missing data problem. We utilize multiple
imputation with a multilevel structure to improve estimates for the
baseball aging curves. Through simulation, we highlight that ignoring
the missing seasons leads to an overestimation of the age effect on
baseball offensive performance. With imputation, we achieve a better
aging curve which shows that players actually decline faster as they get
older than previously estimated.

There are many limitations in our study which leave room for improvement
in future work. In our current imputation model, age is the only
predictor for estimating performance. It is possible to include more
covariates in the algorithm and determine whether a better aging curve
estimate is achieved. In particular, we could factor in other baseball
offensive statistics (e.g.~home run rate, strikeout rate, WOBA, walk
rate,\ldots) in building an imputation model for OPS.

Furthermore, the aging curve estimation problem can be investigated in a
completely different statistical setting. As noted in Section
\ref{sec:lit}, rather than considering discrete observations, another
way of studying aging curves is through a continuous approach, assuming
a smooth curve for career performance. As pointed out by Wakim \& Jin
(2014), methods such as functional data analysis (FDA) and principal
components analysis through conditional expectation (PACE) possess many
modeling advantages, in regard to flexibility and robustness. There
exists a number of proposed multiple imputation algorithms for
functional data (Ciarleglio et al., 2021; He et al., 2011; Rao \&
Reimherr, 2021), which all can be applied in future studies on aging
curves in sports.

\hypertarget{references}{%
\section*{References}\label{references}}
\addcontentsline{toc}{section}{References}

\hypertarget{refs}{}
\begin{CSLReferences}{1}{0}
\leavevmode\vadjust pre{\hypertarget{ref-Albert2002smoothing}{}}%
Albert, J. (2002). \emph{Smoothing career trajectories of baseball
hitters}. Unpublished manuscript, Bowling Green State University.
\url{https://bayesball.github.io/papers/career_trajectory.pdf}

\leavevmode\vadjust pre{\hypertarget{ref-Albert1999comment}{}}%
Albert, J. (1999). Bridging different eras in sports: comment.
\emph{Journal of the American Statistical Association}, \emph{94}(447),
677. \url{https://doi.org/10.2307/2669974}

\leavevmode\vadjust pre{\hypertarget{ref-Albert1992bayesian}{}}%
Albert, J. (1992). A bayesian analysis of a poisson random effects model
for home run hitters. \emph{The American Statistician}, \emph{46}(4),
246. \url{https://doi.org/10.2307/2685306}

\leavevmode\vadjust pre{\hypertarget{ref-Bates2015lme4}{}}%
Bates, D., Mächler, M., Bolker, B., \& Walker, S. (2015). Fitting linear
mixed-effects models using {lme4}. \emph{Journal of Statistical
Software}, \emph{67}(1), 1--48.
\url{https://doi.org/10.18637/jss.v067.i01}

\leavevmode\vadjust pre{\hypertarget{ref-Berry1999bridging}{}}%
Berry, S. M., Reese, C. S., \& Larkey, P. D. (1999). Bridging different
eras in sports. \emph{Journal of the American Statistical Association},
\emph{94}(447), 661--676.
\url{https://doi.org/10.1080/01621459.1999.10474163}

\leavevmode\vadjust pre{\hypertarget{ref-Bradbury2009peak}{}}%
Bradbury, J. C. (2009). Peak athletic performance and ageing: Evidence
from baseball. \emph{Journal of Sports Sciences}, \emph{27}(6),
599--610. \url{https://doi.org/10.1080/02640410802691348}

\leavevmode\vadjust pre{\hypertarget{ref-Brander2014estimating}{}}%
Brander, J. A., Egan, E. J., \& Yeung, L. (2014). Estimating the effects
of age on {NHL} player performance. \emph{Journal of Quantitative
Analysis in Sports}, \emph{10}(2).
\url{https://doi.org/10.1515/jqas-2013-0085}

\leavevmode\vadjust pre{\hypertarget{ref-Ciarleglio2021elucidating}{}}%
Ciarleglio, A., Petkova, E., \& Harel, O. (2021). Elucidating age and
sex-dependent association between frontal {EEG} asymmetry and
depression: An application of multiple imputation in functional
regression. \emph{Journal of the American Statistical Association},
\emph{117}(537), 12--26.
\url{https://doi.org/10.1080/01621459.2021.1942011}

\leavevmode\vadjust pre{\hypertarget{ref-Dendir2016soccer}{}}%
Dendir, S. (2016). When do soccer players peak? A note. \emph{Journal of
Sports Analytics}, \emph{2}(2), 89--105.
\url{https://doi.org/10.3233/jsa-160021}

\leavevmode\vadjust pre{\hypertarget{ref-Fair2007estimated}{}}%
Fair, R. C. (2007). Estimated age effects in athletic events and chess.
\emph{Experimental Aging Research}, \emph{33}(1), 37--57.
\url{https://doi.org/10.1080/03610730601006305}

\leavevmode\vadjust pre{\hypertarget{ref-Fair2008estimated}{}}%
Fair, R. C. (2008). Estimated age effects in baseball. \emph{Journal of
Quantitative Analysis in Sports}, \emph{4}(1).
\url{https://doi.org/10.2202/1559-0410.1074}

\leavevmode\vadjust pre{\hypertarget{ref-Fair1994fast}{}}%
Fair, R. C. (1994). {How Fast Do Old Men Slow Down?} \emph{The Review of
Economics and Statistics}, \emph{76}(1), 103--118.
\url{https://ideas.repec.org/a/tpr/restat/v76y1994i1p103-18.html}

\leavevmode\vadjust pre{\hypertarget{ref-Friendly2021Lahman}{}}%
Friendly, M., Dalzell, C., Monkman, M., \& Murphy, D. (2021).
\emph{{Lahman: Sean 'Lahman' Baseball Database}}.
\url{https://CRAN.R-project.org/package=Lahman}

\leavevmode\vadjust pre{\hypertarget{ref-Gelman2006data}{}}%
Gelman, A., \& Hill, J. (2006). \emph{Data analysis using regression and
multilevel/hierarchical models}. Cambridge University Press.
\url{https://doi.org/10.1017/cbo9780511790942}

\leavevmode\vadjust pre{\hypertarget{ref-He2011functional}{}}%
He, Y., Yucel, R., \& Raghunathan, T. E. (2011). A functional multiple
imputation approach to incomplete longitudinal data. \emph{Statistics in
Medicine}, \emph{30}(10), 1137--1156.
\url{https://doi.org/10.1002/sim.4201}

\leavevmode\vadjust pre{\hypertarget{ref-Judge2020approach}{}}%
Judge, J. (2020a). \emph{An approach to survivor bias in baseball}.
BaseballProspectus.com.
\url{https://www.baseballprospectus.com/news/article/59491/an-approach-to-survivor-bias-in-baseball}

\leavevmode\vadjust pre{\hypertarget{ref-Judge2020delta}{}}%
Judge, J. (2020b). \emph{{The Delta Method, Revisited: Rethinking Aging
Curves}}. BaseballProspectus.com.
\url{https://www.baseballprospectus.com/news/article/59972/the-delta-method-revisited}

\leavevmode\vadjust pre{\hypertarget{ref-Kasim1998application}{}}%
Kasim, R. M., \& Raudenbush, S. W. (1998). Application of gibbs sampling
to nested variance components models with heterogeneous within-group
variance. \emph{Journal of Educational and Behavioral Statistics},
\emph{23}(2), 93--116. \url{https://doi.org/10.3102/10769986023002093}

\leavevmode\vadjust pre{\hypertarget{ref-Kovalchik2014older}{}}%
Kovalchik, S. A. (2014). The older they rise the younger they fall: Age
and performance trends in men's professional tennis from 1991 to 2012.
\emph{Journal of Quantitative Analysis in Sports}, \emph{10}(2).
\url{https://doi.org/10.1515/jqas-2013-0091}

\leavevmode\vadjust pre{\hypertarget{ref-Lahman2021baseball}{}}%
Lahman, S. (1996 -- 2021). \emph{Lahman{'}s baseball database}.
SeanLahman.com.
\url{https://www.seanlahman.com/baseball-archive/statistics/}

\leavevmode\vadjust pre{\hypertarget{ref-Lailvaux2014trait}{}}%
Lailvaux, S. P., Wilson, R., \& Kasumovic, M. M. (2014). Trait
compensation and sex-specific aging of performance in male and female
professional basketball players. \emph{Evolution}, \emph{68}(5),
1523--1532. \url{https://doi.org/10.1111/evo.12375}

\leavevmode\vadjust pre{\hypertarget{ref-Lichtman2009baseball}{}}%
Lichtman, M. (2009). \emph{How do baseball players age? (Part 2)}. The
Hardball Times.
\url{https://tht.fangraphs.com/how-do-baseball-players-age-part-2}

\leavevmode\vadjust pre{\hypertarget{ref-Little1987statistical}{}}%
Little, R. J. A., \& Rubin, D. B. (1987). \emph{Statistical analysis
with missing data}. Wiley.

\leavevmode\vadjust pre{\hypertarget{ref-marchi2018analyzing}{}}%
Marchi, M., Albert, J., \& Baumer, B. S. (2018). \emph{Analyzing
baseball data with {R}} (2nd ed., p. 83). Boca Raton, FL: Chapman;
Hall/CRC Press. \url{https://doi.org/10.1201/9781351107099}

\leavevmode\vadjust pre{\hypertarget{ref-Moore1975study}{}}%
Moore, D. H. (1975). A study of age group track and field records to
relate age and running speed. \emph{Nature}, \emph{253}(5489), 264--265.
\url{https://doi.org/10.1038/253264a0}

\leavevmode\vadjust pre{\hypertarget{ref-Morris1983parametric}{}}%
Morris, C. N. (1983). Parametric empirical bayes inference: Theory and
applications. \emph{Journal of the American Statistical Association},
\emph{78}(381), 47--55.
\url{https://doi.org/10.1080/01621459.1983.10477920}

\leavevmode\vadjust pre{\hypertarget{ref-Page2013effect}{}}%
Page, G. L., Barney, B. J., \& McGuire, A. T. (2013). Effect of
position, usage rate, and per game minutes played on {NBA} player
production curves. \emph{Journal of Quantitative Analysis in Sports},
\emph{0}(0), 1--9. \url{https://doi.org/10.1515/jqas-2012-0023}

\leavevmode\vadjust pre{\hypertarget{ref-R2022language}{}}%
R Core Team. (2022). \emph{R: A language and environment for statistical
computing}. R Foundation for Statistical Computing.
\url{https://www.R-project.org/}

\leavevmode\vadjust pre{\hypertarget{ref-Rao2021modern}{}}%
Rao, A. R., \& Reimherr, M. (2021). Modern multiple imputation with
functional data. \emph{Stat}, \emph{10}(1).
\url{https://doi.org/10.1002/sta4.331}

\leavevmode\vadjust pre{\hypertarget{ref-Rubin1987multiple}{}}%
Rubin, D. B. (Ed.). (1987). \emph{Multiple imputation for nonresponse in
surveys}. John Wiley {\&} Sons, Inc.
\url{https://doi.org/10.1002/9780470316696}

\leavevmode\vadjust pre{\hypertarget{ref-Rubin1976inference}{}}%
Rubin, D. B. (1976). Inference and missing data. \emph{Biometrika},
\emph{63}(3), 581--592. \url{https://doi.org/10.1093/biomet/63.3.581}

\leavevmode\vadjust pre{\hypertarget{ref-Schall2000career}{}}%
Schall, T., \& Smith, G. (2000). Career trajectories in baseball.
\emph{{CHANCE}}, \emph{13}(4), 35--38.
\url{https://doi.org/10.1080/09332480.2000.10542233}

\leavevmode\vadjust pre{\hypertarget{ref-Schell2005baseball}{}}%
Schell, M. J. (2005). Calling it a career: Examining player aging. In
\emph{Baseball's all-time best sluggers: Adjusted batting performance
from strikeouts to home runs} (pp. 45--57). Princeton University Press.
\url{https://www.jstor.org/stable/j.ctt19705ks}

\leavevmode\vadjust pre{\hypertarget{ref-Schuckers2021observed}{}}%
Schuckers, M., Lopez, M., \& Macdonald, B. (2021). \emph{What does not
get observed can be used to make age curves stronger: Estimating player
age curves using regression and imputation}. arXiv.
\url{https://doi.org/10.48550/arXiv.2110.14017}

\leavevmode\vadjust pre{\hypertarget{ref-Schulz1988peak}{}}%
Schulz, R., \& Curnow, C. (1988). Peak performance and age among
superathletes: Track and field, swimming, baseball, tennis, and golf.
\emph{Journal of Gerontology}, \emph{43}(5), 113--120.
\url{https://doi.org/10.1093/geronj/43.5.p113}

\leavevmode\vadjust pre{\hypertarget{ref-Schulz1994relationship}{}}%
Schulz, R., Musa, D., Staszewski, J., \& Siegler, R. S. (1994). The
relationship between age and major league baseball performance:
Implications for development. \emph{Psychology and Aging}, \emph{9}(2),
274--286. \url{https://doi.org/10.1037/0882-7974.9.2.274}

\leavevmode\vadjust pre{\hypertarget{ref-Turtoro2019flexible}{}}%
Turtoro, C. (2019). \emph{Flexible aging in the NHL using GAM}. RPubs.
\url{https://rpubs.com/cjtdevil/nhl_aging}

\leavevmode\vadjust pre{\hypertarget{ref-Vaci2019large}{}}%
Vaci, N., Cocić, D., Gula, B., \& Bilalić, M. (2019). Large data and
bayesian modeling{\textemdash}aging curves of {NBA} players.
\emph{Behavior Research Methods}, \emph{51}(4), 1544--1564.
\url{https://doi.org/10.3758/s13428-018-1183-8}

\leavevmode\vadjust pre{\hypertarget{ref-Vanbuuren1999flexible}{}}%
van Buuren, S., \& Groothuis-Oudshoorn, C. G. M. (1999). \emph{Flexible
multivariate imputation by MICE} (Vol. PG/VGZ/99.054). TNO Prevention;
Health.

\leavevmode\vadjust pre{\hypertarget{ref-vanBuuren2011mice}{}}%
van Buuren, S., \& Groothuis-Oudshoorn, K. (2011). {mice}: Multivariate
imputation by chained equations in {R}. \emph{Journal of Statistical
Software}, \emph{45}(3), 1--67.
\url{https://doi.org/10.18637/jss.v045.i03}

\leavevmode\vadjust pre{\hypertarget{ref-Villaroel2011elite}{}}%
Villaroel, C., Mora, R., \& Parra, G. C. G. (2011). Elite triathlete
performance related to age. \emph{Journal of Human Sport and Exercise},
\emph{6}(2 (Suppl.)), 363--373.
\url{https://doi.org/10.4100/jhse.2011.62.16}

\leavevmode\vadjust pre{\hypertarget{ref-Wakim2014functional}{}}%
Wakim, A., \& Jin, J. (2014). \emph{Functional data analysis of aging
curves in sports}. arXiv. \url{https://doi.org/10.48550/arXiv.1403.7548}

\end{CSLReferences}

\end{document}


\section{Introduction}

	The fast-growing field of quantum thermodynamics brings together
	methods and tools from a variety of research areas ranging
	from quantum information,
	open quantum systems, quantum optics, non equilibrium thermodynamics,
	theory of fluctuations, estimation theory
	and so on \cite{Vinjanampathy_2016,Kosloff_2013,Millen_2016,quan2007quantum,Threelevelmaser,alicki2014quantum,rivas2020strong,goold2016role,PhysRevE.103.032130,Albash_2012,RevModPhys.91.025001,Esposito_2010,RevModPhys.81.1,Parrondo,Popescu2006,Alicki1979TheQO, AJM2008,PRXQuantum.3.030315,RevModPhys.93.035008,esposito2015quantum, VTJ2020,  Maffei2021, Rubio2021, LandiPRA2022}. The quantum extensions of the
	concepts of heat, work, and entropy have in turn led to generalizations
	of the classical heat cycles. The so-called
	quantum thermal machines exploit new thermodynamic
	resources such as quantum entanglement, coherence,
	quantum interactions and quantum statistics \cite{PhysRevA.99.062103,Alet2021,zhang2007four,albayrak2013entangled,Diaz_de_la_Cruz_2014,Pekola2015, Myers_BOF2020,
	WatanabePRL2020}. For instance,
	quantum Otto engine (QOE) based on various platforms has
	been widely studied for its possible quantum
	advantages---both in its quasi-static formulations
	as well as the ones based on time-dependent constraints
	\cite{PhysRevE.90.062134,PhysRevE.101.012116,das2020quantum,chand2021finite,
		lee2020finite,turkpencce2019coupled,thomas2014friction,
    geva1992quantum,feldmann2000performance,ccakmak2017irreversible,Shastri2022,PhysRevB.101.054513,Papadatos_2021,e21111131}.
	A QOE offers conceptual simplicity by virtue of
	 a clear separation of heat and work steps in its heat cycle.
	The quantum working medium used in these models may be
	taken in the form of spins, quantum harmonic oscillator,
	interacting systems and so on \cite{e22070755,lin2003performance,Rezek2006,
		zhang2008entangled,hubner2014spin,
		azimi2014quantum,insinga2016thermodynamical,mehta2017quantum}.
	Further, theoretical models have
	motivated experimental realizations which promise
	to be a boost for future applications in devices \cite{peterson2019experimental, MyersAVS2022}.

	One of the prominent analytical tools
	guiding the theoretical developments is
	the notion of majorisation \cite{Marshallmajorisationbook,TSagava,
		Bhatia1996MatrixA,PhysRevA.95.012110,Joe1990MajorisationAD,
		2020,vedral2015,10.5555/2011326.2011331} and its generalizations \cite{PhysRevResearch.2.033455,2016,RUCH1980222}.
	The concept, originally perhaps from matrix analysis,
	finds  wide applications in various areas of science, mathematics, economics, social sciences, including modern applications to
	entanglement theory and thermodynamic resource theories.
	The majorisation partial order was developed to quantify the notion
	of disorder, in a relative sense, when comparing probability distributions. Transformations between pure
	bipartite states by means of local operations and classical
	communication can be determined in terms of majorisation of the
	Schmidt coefficients of the states \cite{PhysRevLett.83.436,
		PhysRevA.91.052120,PhysRevLett.83.3566,Horodecki2003}. Majorisation has been
	shown to determine
	the possibility of state transformations in the resource theories
	of entanglement, coherence and purity.
	It also provides the first
	complete set of necessary and sufficient conditions for arbitrary
	quantum state transformations under thermodynamic processes
	\cite{majorisationcomplete,2015,2013}, which rigorously accounts for quantum
	coherence among other quantum mechanical properties.

	A majorisation relation may be defined in one of the equivalent ways, as follows.
	Suppose $x^{\downarrow}=(x_{1},...,x_{n}$) and $y^{\downarrow}=(y_{1},...,y_{n}$) are two real $n$-dimensional vectors, where $x_{}^{\downarrow}$ and $ y_{}^{\downarrow}$ indicate
	that the elements are taken in the descending order.
	Then, the vector $x$ is said to be majorised by the vector $y$,
	denoted as $x\prec y$, if,
	for each $m = 1,\dots, n$, we have
	\begin{equation}
		\sum_{k=m}^{n}x_{k}^{} \geq \sum_{k=m}^{n}y_{k}^{},
	\end{equation}
	with equality holding for $m=1$. The notion of majorisation can
	be readily applied to compare how two probability distributions0 deviate
	from a uniform distribution $u = (1/n,\dots,1/n)$.
	Thus, $x\prec y$ implies that the distribution $y$ is more ordered
	than $x$. An important consequence of this relation is
	the inequality: $\sum_{k=1}^{n} f(x_k) \geq \sum_{k=1}^{n} f(y_k)$, where
	$f$ is any continuous, real-valued concave function.
	For example, the relation $x\prec y$ implies that
	the corresponding
	Shannon entropies are related as: $S(x) \geq S(y)$, where
	$S(x) = - \sum_{k=1}^{n} x_k \ln x_k$.
	Further, there are many equivalent characterizations of the majorisation
	relation. For instance, $x$ is majorised by $y$ only when $y$
	can be obtained from $x$ by the action of a bistochastic matrix
	\cite{Marshallmajorisationbook}.

	In the present work, we characterize the
	operation of a QOE through the notion of majorisation.
	We show that a spin-based quantum working substance
	provides a natural platform by which the majorisation conditions
	characterize the operation of
	a QOE. Thus, majorisation provides
	sufficient criteria for the operation of a spin-based Otto engine.
	In fact, the analysis can be extended to a model of two spins
	coupled via Heisenberg exchange interactions. Further, majorisation
	provides insight into the local thermodynamics of individual
	spins in the coupled model. Using the majorisation conditions,
	we also validate an upper bound for Otto efficiency in
	the coupled case, which is
	tighter than the Carnot value.

	The paper is organized as follows. In Section II,
	we describe the quantum Otto cycle and its various
	stages based on a quantum
	working substance.
	In Section III, we express the work output in terms
	of the relative entropy between the two  equilibrium
	distributions corresponding to hot and cold reservoirs,
	and we show that a greater value of the
	Shannon entropy of the system at the hot reservoir
	(as compared to the cold reservoir) does not ensure that
	net work may be extracted in the Otto cycle.
	In Section IV, we show how the majorisation
	relation, $P \prec P'$, between the hot and cold
	reservoir equilibrium distributions lead to
	positive work condition for the QOE.
	This is shown in Section IV.A for a single spin-$s$.
	In Section V, we analyze the coupled spins model,
	showing our main results for a special case of $(1/2,1)$
	system. In Section VI, local work by individual
	spins is analyzed based on global conditions.
	Lastly, Section VII shows proves the enhancement
	in Otto efficiency based on majorisation conditions.
	We end our paper by summarizing our main results
	in Section VIII. The derivations of various results are
	presented in Appendix.

	\section{Quantum Otto Engine (QOE)}

	The classical Otto cycle is a textbook example of
	a four-step heat cycle in which a classical
	working medium, in the form of a gas or air-fuel mixture,
	undergoes two adiabatic and two isochoric steps \cite{Zemansky}.
	A QOE is based on
	a quantum generalization of the classical cycle,
	in which the quantum working substance
	(referred to below as the system)
	undergoes two quantum adiabatic steps and two isochoric
	steps. Here, we are interested in a quasi-static Otto cycle
	in which each of the steps can take an arbitrarily long time.
	A quantum adiabatic process, either in the
	compression or expansion stage to be explained
	below, is performed by
	varying an externally controllable parameter.
	Secondly, for such a process,
	the quantum adiabatic theorem \cite{Adiabatic_Fock} is assumed to hold so
	the process does not cause any transitions between
	the energy levels, thus preserving their occupation probabilities.
	The remaining two steps are
	the isochoric heating and cooling processes which involve thermal
	interaction of the system with the hot or cold reservoir.
	Here, the system gets enough time to reach thermal
	equilibrium with the corresponding reservoir.
	The heat cycle is described in a more quantitative
	detail as follows.

	\begin{figure}[h]
		\centering
		{\includegraphics[width=0.45\linewidth]{Otto5.jpeg}\label{fig:fig1}}
		\caption{Schematic of a quantum Otto cycle consisting of the stages
			$1\to 2\to 3 \to 4 \to 1$.
			Steps $1\to 2$
			and $3\to 4$ respectively denote quantum adiabatic expansion
			and compression processes, which preserve occupation probability  and involve work $W_{1\to 2}^{} <0$
			and $W_{3\to 4}^{} > 0$.
			The net work extracted in one cycle is $W = W_{1\to 2}^{} + W_{3\to 4}^{} < 0$. In step $4 \to 1$, heat
			$Q_1 >0$ enters the working substance from a hot reservoir,
			while in step $2 \to 3$, heat
			$Q_2 < 0$ is rejected to the cold reservoir.
			Conservation of energy implies,
			$Q_1 + Q_2 + W =0$.
		}
		\label{fig:fig3}
	\end{figure}
	\par\noindent
	\textbf{Stage 1.}
	Consider an $n$-level quantum system
	with Hamiltonian $H(B_{1})$ whose eigenvalues
	can be arranged (in descending order) as:
	$\varepsilon^{\downarrow} =
	(\varepsilon_{n}, \dots, \varepsilon_{1})$.
	The system is in thermal equilibrium
	with the hot reservoir at
	temperature $T_1$. The canonical occupation probabilities
	for different energy levels, $P_{k}=e^{-\varepsilon_{k}/T_{1}}/{\sum_{k} e^{-\varepsilon_{k}/T_{1}}}$, are arranged as
	$P^{\downarrow}=(P_{1},\dots,P_{n}$),
	where we have set the Boltzmann constant $k_{\rm B}$ equal to unity.
	The energy eigenstates are represented by
	$\{\ket{\psi_{k}} \vert k=1,...n\} $.
	Thus, the density matrix representing the thermal
	state of the system is given by:
	\begin{equation} \label{eqn6}
		\rho=\sum_{k=1}^{n}P_{k}\ket{\psi_{k}}\bra{\psi_{k}},
	\end{equation}
	%
	\textbf{Stage 2.} The system is detached from the
	hot reservoir and undergoes a quantum adiabatic process
	in which the external field strength is {\it lowered} from $B_{1}$ to
	$B_{2}$. Here, the quantum adiabatic theorem
	ensures that no transitions are induced between
	the energy levels in the change from
	$\varepsilon_{k}$ to $\varepsilon_{k}^{\prime}$. Suppose that the energies after the first adiabatic process
	are given by:
	${\varepsilon^{\prime}}^{\downarrow} =
	(\varepsilon_{n}^{\prime}, \dots, \varepsilon_{1}^{\prime})$,
	where we assume no level-crossing as the Hamiltonian changes from $H(B_1)$ to $H(B_2)$.
	\par\noindent
	\textbf{Stage 3.} The system is brought in thermal contact with the cold reservoir at temperature $T_{2} (<T_{1})$. The energy eigenvalues remain at $\varepsilon_{k}^{\prime}$
	while the occupation probabilities change from $P_{k}$ to $P_{k}^{\prime}=e^{-\varepsilon_{k}^{\prime}/T_{2}}
	/{\sum_{k} e^{-\varepsilon_{k}^{\prime}/T_{2}}}$, which are ordered as:
	$P^{{\prime}{\downarrow}} = (P_{1}^{\prime},\dots, P_{n}^{\prime}$).
	Thus, the density matrix of the system at the end of Stage-3
	is given by:
	\begin{equation} \label{eqn66}
		\rho^{\prime}=\sum_{k=1}^{n}P_{k}^{\prime}\ket{\psi_{k}}\bra{\psi_{k}}.
	\end{equation}
	\par\noindent
	\textbf{Stage 4.} The system is detached from the cold reservoir and the field strength is changed back
	to $B_{1}$.
	The occupation
	probabilities $\{P_{k}^{\prime}\}$ remain unchanged, while the energy levels change back
	from $\{\varepsilon_{k}^{\prime}\}$ to
	$\{\varepsilon_{k}\}$.

	Finally, the system is attached to the hot reservoir again whereby
	the initial state ($\rho$) is recovered, thus
	completing one heat  cycle. Note that only heat is exchanged between the system and the reservoir
	during an isochoric process,
	which is given by  the difference between the final and initial
	mean energies of the system in that process. Thus, in
	Stage-1 and Stage-3, the heat exchanged is given respectively as:
	\begin{align} \label{eqn7}
		Q_{1}&=\sum_{k=1}^n\varepsilon_{k}(P_{k}-P_{k}^{\prime}),~~~~~
		Q_{2} =\sum_{k=1}^n\varepsilon_{k}^{\prime}(P_{k}^{\prime}-P_{k}).
	\end{align}
	On the other hand, only work is performed during
	the adiabatic branches of the
	quantum Otto cycle. Let $W$ be the net work performed in one cycle. Applying
	the law of conservation of energy to the cyclic process,
	we have: $Q_1 + Q_2 + W =0$.
	The operation of a heat engine requires that
	heat is absorbed (rejected)
	by the system at the hot (cold) reservoir,
	while net work is extracted from the system by
	the end of the cycle.
	These conditions can be satisfied by choosing the sign convention: $Q_1 >0$, $Q_2 <0$
	and $W <0$. The net
	work performed by the QOE can then be written as
	\begin{equation} \label{eqn8}
		\vert W \vert = Q_{1}+Q_{2} = \sum_{k=1}^{n}(\varepsilon_{k}-
		\varepsilon_{k}^{\prime})(P_{k}-P_{k}^{\prime}).
	\end{equation}
	We denote $\vert W \vert \geq 0$ as the positive work condition (PWC)
	of our engine. The efficiency of the QOE
	is defined as $\eta=\vert W \vert/Q_{1} = 1+ Q_2/Q_1$.

	\section{Relative entropy and QOE}

	In this section, we cast the thermodynamic
	quantities for a QOE in terms of the relative
	entropy which is defined as
	$D(x||y) \equiv \sum_k x_{k} \left(
	\ln x_{k} - \ln y_{k} \right) \geq 0$.
	Also known as the Kullback-Leibler divergence\cite{10.1214/aoms/1177729694},
	this quantity is a measure of the 'distance'
	between two discrete probability distributions,
	and vanishes only when
	the distributions $x$ and $y$ are identical.

	Now, the expressions for canonical probabilities
	may be inverted as:
	$\varepsilon_{k} \equiv -T_{1}\ln{(P_{k}
		\sum_{j} e^{-\varepsilon_{j}/T_{1}})}$ and
	$\varepsilon_{k}^{\prime} \equiv -T_{2}\ln{(P_{k}^{\prime}
		\sum_{j} e^{-\varepsilon_{j}^{\prime}/T_{2}})}$.
	Substituting these expressions in
	Eq. (\ref{eqn7}),  and after some algebra
	(Appendix A), the heat exchanged with each
	the reservoir is expressed as:
	\begin{align}
		Q_1 & = T_1 (S_1-S_2) - T_1 D(P^{\prime}||P),
		\label{q1D}
		\\
		%
		Q_2 & =  -T_2 (S_1-S_2) - T_2 D(P||P^{\prime}),
		\label{q2D}
	\end{align}
	where $S_1 =-\sum_{k} P_{k} \ln P_{k}$ and $S_2 =-\sum_{k} P_{k}' \ln P_{k}'$
	are the Shannon entropies of the system in equilibrium
	with hot and cold reservoirs, respectively. as we have let $k_B=1$ Shannon entropy is equal to canonical entropy.

	Thus, it can be seen that $D(P^{\prime}||P)$
	is equal to the entropy generated in the
	hot isochoric step.
	$D(P||P^{\prime})$ has a similar meaning
	for the cold isochoric step.
	The net work extracted in an Otto cycle is given by:
	\begin{equation}\label{mj6}
		|W|=(T_{1}-T_{2}) (S_{1}-S_{2})-T_{1}D(P^{\prime}||P)-T_{2}D(P||P^{\prime}).
	\end{equation}
	%
	Using, Eq. (\ref{q1D}) and Eq. (\ref{q2D}) the total entropy generated
	in the heat cycle,
	$\Delta_{\rm tot} S = -Q_2/T_2 - Q_1/T_1$,
	can be expressed as:
	\begin{equation}
		\Delta_{\rm tot} S = D(P||P^{\prime}) + D(P^{\prime}||P).
		\label{dstot}
	\end{equation}
	The total entropy generated in a quantum Otto cycle
	is thus equal to the symmetric sum of the relative
	entropies. This quantity is also known as the
	symmetrized divergence and is distinguished by
	the fact that it serves as a metric in the
	space of probability distributions.
	Finally, the positivity of the total entropy
	generated proves the consistency of the QOE
	with the second law of thermodynamics, and hence
	its efficiency is bounded by the Carnot value:
	$\eta \leq 1-T_2/T_1$.

	Returning to Eq. (\ref{mj6}), the positivity
	of the relative entropy indicates that
	for $T_1 >T_2$,
	$S_1 >S_2$ is a necessary,
	but not a sufficient condition for
	$|W| \geq 0$. The necessity of the condition $S_1>S_2$
	can be reasoned due to the fact that
	heat is absorbed by the system at the hot reservoir,
	while heat is rejected by the system at the cold reservoir,
	and  the intermediate, quantum adiabatic processes do not
	alter the entropy of the system.
	These considerations suggest that more general conditions
	are desirable to characterize
	the probability distributions, which not only
	ensure $S_1>S_2$, but also the PWC or $|W| \geq 0$.
	In this paper, we show
	that the majorisation relation ($P \prec  P'$) provides
	sufficient conditions for the operation of
	a spins-based quantum Otto cycle as a heat engine.

	\section{Majorisation and QOE}

	As mentioned earlier, the majorisation relation $P\prec P^{\prime}$
	implies the following set of inequalities:
	\begin{equation}
		\sum_{k=m}^{n}P_{k}^{} \geq \sum_{k=m}^{n}P_{k}^{\prime},
		\quad (m =1,\dots,n)
		\label{majn2}
	\end{equation}
	with the equality holding for $m=1$ owing to the normalization
	property of each distribution.

	Specifically, we obtain from the above inequalities, for
	$m=n$
	\begin{equation} \label{eqna}
		P_{n}\geq P_{n}^{\prime},
	\end{equation}
	and, for $m=2$, along with normalization
	\begin{equation} \label{majn1}
		P_{1}^{\prime} \geq P_1.
	\end{equation}
	These inequalities may be combined as:
	$ {P_{1}^{\prime}}/{P_{n}^{\prime}} \geq {P_{1}}/{P_{n}}$,
	to yield the condition:
	\begin{equation} \label{majcond1}
		\frac{\varepsilon_{n}^{\prime}- \varepsilon_{1}^{\prime}}{T_{2}}
		\geq \frac{\varepsilon_{n}^{}- \varepsilon_{1}^{}}{T_{1}}.
	\end{equation}
	For $T_1>T_2$, the above inequality will yield a nontrivial condition,
	provided that $\varepsilon_{n} - \varepsilon_{1}>\varepsilon_{n}^{\prime}-
	\varepsilon_{1}^{\prime}$. In other words, we must assume that the range of
	the energy spectrum shrinks during the first quantum adiabatic process. Apart from that, the condition (\ref{majcond1})
	is derived for a generic, non-degenerate spectrum.

	Now, an important question arises regarding the circumstances
	under which the majorisation inequalities, Eq. (\ref{majn2}), hold.
	Naturally, this is dependent on the form of Hamiltonian
	(or the energy spectrum which enters the expressions for
	the canonical probabilities).
	In the following, we show for a spin system,
	a set of necessary
	and sufficient conditions to satisfy the majorisation relation.

	\subsection{QOE with a single spin-$s$}

	Suppose the system is in the form of
	a quantum spin of magnitude $s$. The energy
	eigenvalues in Stage-1 are: $\varepsilon_k = 2(k -s -1) B_1$, where
	$k =1, \dots, 2s+1$. Explicitly, we have
	\begin{equation*}
		\varepsilon_{1}=-2sB_{1}, \; \varepsilon_{2}=-2(s-1)B_{1},
		\dots,
		\; \varepsilon_{2s}^{}=2(s-1)B_{1}, \; \varepsilon_{2s+1}^{}=2sB_{1}.
	\end{equation*}
	After the first quantum adiabatic step,
	the energy spectrum is given by:  $\varepsilon_{k}^{\prime} = 2(k -s -1) B_2$, where $k =1, \dots, 2s+1$ and $B_2 < B_1$.

	Now, for this system,
	Eq. (\ref{majcond1})
	simplifies to the following condition:
	\begin{equation}
		\frac{B_{2}}{T_{2}}\geq \frac{B_{1}}{T_{1}}.
		\label{bt12}
	\end{equation}
	The above condition was first derived in Ref. \cite{kieu2004second}
	for a two-level quantum system (equivalent to $s=1/2$).
	Here, we see
	it as a consequence of the majorisation relation
	between the canonical distributions corresponding to
	hot and cold reservoirs. In fact, the above condition
	is necessary and sufficient to satisfy all the majorisation
	inequalities, (\ref{majn2}) \cite{2021}.

	Next, using the definitions in Eq. (\ref{eqn7}),
	the  heat exchanged between the system
	and each reservoir is calculated to be:
	\begin{align} \label{eqn52}
		Q_{1}&=2B_{1}X, \quad Q_{2}=-2B_{2}X,
	\end{align}
	where

		$X =\sum_{k=2}^{2s+1}(k-1)\left(P_{k}-P_{k}^{\prime}\right)$.

	Thus, the magnitude of the work performed in one cycle is:
	\begin{equation} \label{eqn53}
		\vert W \vert=2(B_{1}-B_{2})X.
	\end{equation}
	Now, assuming the relation
	$P\prec P^{\prime}$, and for $n=2s+1$,
	we add up all the inequalities
	(\ref{majn2}) corresponding to $m=2,\dots,n$.
	The resulting inequality can be rewritten in the form
	$X \geq 0$. In other words, we obtain that
	$P\prec P^{\prime}$ implies
	$|W| \geq 0$, provided $B_1>B_2$.

	Thus, we may say that for a spin-$s$ system,
	the majorisation relation $P\prec P^{\prime}$
	is a sufficient and necessary condition for the
	operation of QOE.

	As an extreme case scenario, we may have the
	conditions:   $P_{k} \geq P_{k}^{\prime}$,
	for $k = 2, \dots, n$.

	The normalization property then ensures $P_{1} \leq P_{1}^{\prime}$.
	It is clear that the above inequalities
	satisfy the majorisation conditions (\ref{majn2}),
	implying that $P \prec P^{\prime}$.
	Thus, the above
	extreme case, applied to the case of
	a spin, also leads to PWC (see \cite{2021} for details).

	\section{QOE with two coupled spins}

	Next, we consider a system of two coupled spins, a spin-$1/2$ particle
	interacting with an arbitrary spin-$s$, via 1-d isotropic,  Heisenberg
	exchange interaction. The Hamiltonian of the working substance is
	\begin{equation} \label{eqn5}
		H=2B(s_{z}^{(1)}\otimes I^{(2)}+I^{(1)} \otimes s_{z}^{(2)})+
		8J(s_{x}^{(1)}\otimes s_{x}^{(2)}+s_{y}^{(1)}\otimes s_{y}^{(2)}
		+s_{z}^{(1)}\otimes s_{z}^{(2)}),
	\end{equation}
	%
	where $J\geq 0$ is the coupling strength parameter. Here,
	$s^{(1)} \equiv \{s_{x}^{(1)},s_{y}^{(1)},s_{z}^{(1)}\}$ and $s^{(2)}
	\equiv \{s_{x}^{(2)},s_{y}^{(2)},s_{z}^{(2)}\}$ are the spin-1/2 and
	spin-$s$ operators, respectively, and $I$ denotes the identity operator.
	We set $\hbar=1$,
	Bohr magneton $\mu_{b}=1$, and
	assume there is no orbital angular momentum so that the
	gyromagnetic ratio $\gamma$ is the same for both spins, $\gamma=2$.
	The total number of levels of the bipartite system
	is $n=2(2s+1)$.
	The energy spectrum is displayed in Fig. 1 of SM. Note that
	the energy eigenvalues contain a constant term $4 s J$
	which can be adjusted, for convenience,
	by an overall shift of the energy spectrum.
	Further, note that only the field parameter
	$B$ is varied cyclically while the
	coupling parameter $J$ is held fixed.

	QOE of the above kind was first studied with two spin-1/2 particles \cite{PhysRevE.83.031135} where, amongst other things, an enhancement in Otto efficiency was reported as a result of
	coupling between the spins.
	Further, the model was extended incorporating
	the above Hamiltonian \cite{PhysRevE.92.022142}. As the energy spectrum
	becomes more complex, numerical results were used to gain insights
	into the performance of QOE. In \cite{2021}, a heuristics-based approach was used to analyze the general case of
	spin-$s_1$ coupled with
	spin-$s_2$. Thus, sufficient criteria for work extraction
	were inferred using the extreme case scenarios. It was also
	argued with numerical results that majorisation leads to
	a more robust characterization of QOE than the extreme case
	scenario. Motivated by these findings, in this paper,
	we develop a characterization of the QOE in terms
	of the majorisation relation.

	In the following, we show how the majorisation conditions
	also lead to sufficient criteria for QOE based
	on the above coupled-spins model.
	In order to illustrate our main results, we treat the case
	of a spin-1/2 particle coupled to a spin-1, for which $n=6$.
	The results for the more general case of $(1/2,s)$
	system are reported in SM.

	\subsection{The coupled $(1/2,1)$ system}

	By introducing a constant shift of $4sJ \equiv 4J$
	in the energy eigenvalues at the hot reservoir,
	these are given by:
	$\varepsilon_{1}=-3B_{1}, \; \varepsilon_{2}=-B_{1}-12J,
	\; \varepsilon_{3}=-B_{1},
	\; \varepsilon_{4}=B_{1}-12J,\;\varepsilon_{5}=B_{1},
	\;\varepsilon_{6}=3B_{1}$.
	The corresponding eigenstates are:
	$\ket{\psi_{1}} =\ket{12^\prime}, ~\ket{\psi_{2}}
	=(\sqrt{2}\ket{02^\prime}-\ket{11^\prime})/\sqrt{3}, ~\ket{\psi_{3}}
	=(\ket{02^\prime}+\sqrt{2}\ket{11^\prime})/\sqrt{3},~
	\ket{\psi_{4}}
	=(\ket{01^\prime}-\sqrt{2}\ket{10^\prime})/\sqrt{3}, 	~\ket{\psi_{5}}
	=(\sqrt{2}\ket{01^\prime}+\ket{10^\prime})/\sqrt{3}, 	~\ket{\psi_{6}}=\ket{00^\prime}$,
	where $\ket{0} \equiv (1,0)^T$ and $\ket{1} \equiv (0,1)^T$ are eigen-kets for the bare Hamiltonian of spin-1/2.
	Similarly, $\ket{0^{\prime}}\equiv (1,0,0)^T,~\ket{1^{\prime}} \equiv (0,1,0)^T$ and $\ket{2^{\prime}} \equiv (0,0,1)^T$ are the eigen-kets
	for the bare Hamiltonian of spin-1.
	The density matrix in the initial state
	of the coupled system is:
	$\rho=\sum_{k=1}^{6}P_{k}\ket{\psi_{k}}\bra{\psi_{k}}$.
	Similarly at the cold reservoir, the
	density matrix is
	$\rho^{\prime}=\sum_{k=1}^{6}P_{k}^{\prime}\ket{\psi_{k}}
	\bra{\psi_{k}}$.

	Now, by inspection, for $B_{1}>6J$,
	the energy levels are ordered as:
	$\varepsilon^{\downarrow} =
	(\varepsilon_{6}, \dots,
	\varepsilon_{1})$, and
	therefore,
	$P^{\downarrow} = (P_{1},\dots, P_{6}$).
	Similarly, for $B_{2}>6J$, we have
	${\varepsilon^{\prime}}^{\downarrow} =
	(\varepsilon_{6}^{\prime}, \dots, \varepsilon_{1}^{\prime})$,
	as well as $P^{{\prime}{\downarrow}}=(P_{1}^{\prime},\dots, P_{6}^{\prime}$).
	Then, the heat $Q_{1} = \sum_{k=1}^{6} \varepsilon_{k}(P_{k}-P_{k}^{\prime})$ can be expressed as:
	\begin{equation} \label{eqn56}
		Q_{1}=2B_{1}{\cal X}-12 J {\cal Y},
	\end{equation}
	where
	\begin{align}
		{\cal X}&=3(P_{6}-P_{6}^{\prime})+2(P_{5}-P_{5}^{\prime})+2(P_{4}-P_{4}^{\prime}
		)+(P_{3}-P_{3}^{\prime})+(P_{2}-P_{2}^{\prime}) \label{calx},\\
		{\cal Y} &=(P_{2}-P_{2}^{\prime})+(P_{4}-P_{4}^{\prime}).
		\label{caly}
	\end{align}
	Similarly, we can evaluate:
	$Q_{2}=\sum_{k=1}^{6} \varepsilon_{k}'(P_{k}-P_{k}^{\prime}) = 2B_{2}{\cal X}-12 J {\cal Y}$.
	Thus, the work extracted in one cycle is
	\begin{equation} \label{eqn58}
		|W|=2(B_{1}-B_{2}){\cal X}.
	\end{equation}
	For $B_{1}>B_{2}$, PWC requires ${\cal X} \geq 0$.
	Now, under the relation $P\prec P^{\prime}$,
	the following set of inequalities must hold:
	\begin{align}
		P_{6}&\geq P_{6}^{\prime} \label{eqn60}\\
		P_{5}+P_{6}&\geq P_{6}^{\prime}+P_{5}^{\prime} \label{eqn61}\\
		P_{4}+P_{5}+P_{6}&\geq P_{6}^{\prime}+P_{5}^{\prime}+P_{4}^{\prime} \label{eqn62}\\
		P_{3}+P_{4}+P_{5}+P_{6}&\geq
		P_{6}^{\prime}+P_{5}^{\prime}+P_{4}^{\prime}+P_{3}^{\prime} \label{eqn63}\\
		P_{2}+P_{3}+P_{4}+P_{5}+P_{6}&\geq
		P_{6}^{\prime}+P_{5}^{\prime}+P_{4}^{\prime}+P_{3}^{\prime}+P_{2}^{\prime}. \label{eqn64}
	\end{align}
	Further, as we have seen above,  Eq. (\ref{eqn64}) implies:
	\begin{equation} \label{eqn66}
		P_{1}\leq P_{1}^{\prime}.
	\end{equation}
	Eqs. (\ref{eqn60}) and (\ref{eqn66}) can be combined
	to yield the condition $B_2/T_2 \geq B_1/T_1$, which
	is the same condition as for a QOE based on a single spin.

	Then, upon adding Eqs. (\ref{eqn60}), (\ref{eqn62}) and (\ref{eqn64}),
	we obtain
	%
	\begin{equation}
		P_{2}+P_{3}+2P_{4}+2P_{5}+3P_{6}\geq
		3P_{6}^{\prime}+2P_{5}^{\prime}+2P_{4}^{\prime}+P_{3}^{\prime}+P_{2}^{\prime},
	\end{equation}
	which can be rearranged as the inequality: ${\cal X} \geq0$.
	In this manner, we see that
	the majorisation relation ($P \prec P'$) directly
	implies the positive work condition for the QOE based on the coupled $(1/2,1)$ system, provided $B_1 > B_2$. The proof can
	be straightforwardly generalized to the case of a $(1/2,s)$
	system, as discussed in Section II of SM.

	Now, in order to ascertain the conditions under
	which the majorisation inequalities themselves hold,
	for the given Hamiltonian of the system,
	we employ numerical evidence.
	As Fig. \ref{fig_maj21} shows,
	$B_2/T_2 \geq B_1/T_1$ is a necessary, but not a sufficient
	condition in the case of the coupled system.
	The majorisation relation may be verified only
	for a limited range of $J$ values (for given $T_1,T_2,B_1,B_2$).
	Thus, condition (\ref{eqn64}),
	$\sum_{k=2}^{6} (P_{k}-P_{k}^{\prime})\geq 0$,
	is violated beyond a certain range of $J$.
	Now, it is difficult to estimate this range analytically for,
	say, arbitrary reservoir temperatures. Numerically, it is seen that
	as the temperatures of the reservoirs are raised,
	the range of validity of the majorisation relation (in terms of
	the range of $J$)
	also broadens. Below, we infer  a {\it sufficient} criterion for majorisation, in terms of the permissible range of $J$ values,
	which is followed with a good accuracy at lower temperatures.

	Upon combining Eqs. (\ref{eqn62}) and (\ref{eqn66}), we arrive
	at the condition:	$0 \leq J_{}\leq {\Phi}/{3}$,
	where $\Phi=({B_{2}}/{T_{2}}-{B_{1}}/{T_{1}})
	({1}/{T_{2}}-{1}/{T_{1}})^{-1}$.
	The same condition is obtained upon combining
	Eqs. (\ref{eqn63}) and (\ref{eqn66}). However,
	the validity of Eq. (\ref{eqn66}) requires that
	\begin{equation} \label{eqn69}
		0 \leq J\leq\frac{\Phi}{6} +\frac{1}{12}\left(\frac{1}{T_2}-\frac{1}{T_1}\right)^{-1}\ln\left(\frac{1+e^{-2B_1/T_1}}{1+e^{-2B_2/T_2}}\right)\equiv J_c.
	\end{equation}
	Thus, we may infer that the above range for $J$ is the
	strictest one for which {\it all} majorisation
	conditions (Eqs. (\ref{eqn60}) to (\ref{eqn64}))  hold good.
	Clearly, we have $J_c \geq \Phi/6$.

	\begin{figure}[h!]
		\centering
		{\includegraphics[width=0.7\linewidth]{MJ4.jpeg}}
		\caption{{Majorisation inequalities for the coupled
				$(1/2,1)$ system,
				(Eqs. (\ref{eqn60}) to (\ref{eqn64})).
				$P_{6}-P_{6}^{\prime} \geq 0$ (red, dot curve),
				$\sum_{k=5}^{6}(P_{k}-P_{k}^{\prime})\geq 0$ (green, star dashed curve),
				$\sum_{k=4}^{6}P_{k}-P_{k}^{\prime} \geq 0$ (cyan, dot-dashed curve),
				$\sum_{k=3}^{6}P_{k}-P_{k}^{\prime}\geq 0$ (purple, dashed curve),
				$\sum_{k=2}^{6}P_{k}-P_{k}^{\prime}\geq 0 $ (blue),
				are together verified in a limited range of $J$ values.
				The last inequality (blue curve) yields the
				strictest range of $J$ values, beyond which the
				majorisation relation ($P \prec P'$)  does not hold.
				The parameters are chosen as $B_{1}=5,B_{2}=3,T_{1}=6,T_{2}=3$},
			such that $B_1 > B_2$ and $B_2/T_2 > B_1/T_1$.
			Analytically, we can derive a sufficient
			condition for the majorisation relation to  hold, as
			$0\leq J \leq J_c$ where in this case $J_c=0.189$     (see Eq. (\ref{eqn69}).
		}
		\label{fig_maj21}
	\end{figure}

	\begin{figure}[h!]
		\centering
		{\includegraphics[width=0.7\linewidth]{WKMJ.jpeg}}
		\caption{Work output (Eq. (\ref{eqn58})) and the  majorisation condition
			(Eq. (\ref{eqn64})) for the coupled $(1/2,1)$ system, with the parameters set as
			$B_{1}=5,B_{2}=3,T_{1}=6,T_{2}=3$.
			The majorisation relation (blue curve, as in Fig.
			\ref{fig_maj21}) is violated beyond a certain
			value of the coupling strength $J$. However, the
			work output can be nonzero even beyond that range.}
		\label{fig_wmj}
	\end{figure}
	In summary,
	we have shown that if the majorisation
	relation ($P \prec P'$) is satisfied, then the model works like an engine (PWC).
	However, it is important to note that
	the PWC {\it can} hold
	even if the majorisation relation is not satisfied,
	as depicted in Fig. \ref{fig_wmj}.
	For the coupled $(1/2,1)$
	model, it implies that net work may be obtained
	(${\cal X} > 0$) for $J > J_c$.
	Given that $B_1 > B_2$ along with the
	conditions (\ref{bt12}) and (\ref{eqn69}), all majorisation inequalities hold  and thus
	yield sufficient criteria for PWC in the case  of $(1/2,1)$ system.
	On the other hand, for the single spin, we obtained
	sufficient {\it and} necessary conditions from majorisation.
	Finally, based on induction,
	we can infer sufficient criteria for the general case
	of $(1/2,s)$ system, which are
	the conditions: $B_1 > B_2$, $B_2/T_2 \geq B_1/T_1$
	and  $0 \leq J \leq {\cal J}_c$,
	within which PWC holds for the $(1/2,s)$ system.
	The details are mentioned in Appendix B.

	\section{Analysis of local work}

	Now, each spin constituting the bipartite coupled system
	is governed by a local Hamiltonian which also depends
	on the parameter $B$, and so undergoes a cyclic evolution.
	Thus, it is of intact to examine the local performance
	of each spin in the quantum Otto cycle. The local
	state of each spin is obtained from its
	reduced density matrix. Thus,
	upon summing over the degrees of freedom of spin-1,
	the reduced density matrix for spin-1/2, in Stage-1, is
	defined as:

		$\rho^{(1/2)}=\sum_{m=0^\prime}^{2^\prime}\bra{I\otimes
			m}\rho\ket{I\otimes m}$,

	which may be written in diagonal form,
	as: $\rho^{(1/2)} = \{ q_1, q_2 \}$,
	where $q_2$ is the occupation probability
	of the excited state, given by
	\begin{align}

		q_{2}&=\frac{2P_{2}}{3}+\frac{P_{3}}{3}+\frac{P_{4}}{3}+\frac{2P_{5}}{3}+P_{6} \label{eqn82},
	\end{align}
	and $q_1 = 1-q_2$.

	The work performed by spin-1/2 is evaluated to be:
	$W_{1/2}^{}=2(B_{1}-B_{2})(q_{2}-q_{2}^{\prime})$.
	For convenience,  we express this work in the form:
	\begin{equation} \label{eqn91}
		W_{1/2}=\frac{2}{3} (B_{1}-B_{2})({\cal X} + {\cal Z}),
	\end{equation}
	where ${\cal X}$ has been defined in Eq. (\ref{calx}), and
	\begin{equation} \label{eqn106}
		{\cal Z}=(P_{2}-P_{2}^{\prime})-(P_{4}-P_{4}^\prime).
	\end{equation}
	Similarly, the reduced density matrix for spin-1,
	$\rho^{(1)}=\sum_{m=0}^{1}\bra{m\otimes I}\rho\ket{m\otimes I}$,
	is written as: $\rho^{(1)} \equiv\{r_1,r_2, r_3  \}$, where the occupation probabilities:
	\begin{align}
		r_{1}&=P_{1}+\frac{2P_{2}}{3}+\frac{P_{3}}{3}, \label{eqn84}\\
		r_{2}&=\frac{P_{2}}{3}+\frac{2P_{3}}{3}+\frac{P_{4}}{3}+\frac{2P_{5}}{3}, \label{eqn85}\\
		r_{3}&= \frac{2P_{4}}{3}+\frac{P_{5}}{3}+P_{6}, \label{eqn86}
	\end{align}
	are ordered as: $r_1 > r_2 > r_3$.
	Then, the local work by spin-1 is given by:
	$W_{1}=2(B_{1}-B_{2})[2(r_{3}-r_{3}^{\prime}) + (r_{2}-r_{2}^{\prime})]$,
	which can be rewritten in the form:
	\begin{equation} \label{eqn105}
		W_{1}=\frac{2}{3} (B_{1}-B_{2})(2 {\cal X} - {\cal Z}),
	\end{equation}
	where we have used Eqs. (\ref{eqn85}) and (\ref{eqn86})
	along with the definitions of $\cal X$ and $\cal Z$.
	Now, it is easily verified that the local work contributions,
	Eqs. (\ref{eqn91}) and (\ref{eqn105}),
	add up to yield the global work (Eq. (\ref{eqn58})),
	i.e. $W_{1/2} + W_1 = |W|$. In Fig. \ref{fig_glw},
	we compare the local and global work output.
	It is observed that
	the spin-1/2 work vanishes prior to the spin-1
	work. This observation may be justified as
	follows.
	%

	\begin{figure}[h!]
		\centering
		{\includegraphics[width=0.7\linewidth]{work1.jpeg}}
		\caption{ Global (black, dot-dashed curve)
		versus the local work for spin-1/2 (red, dashed curve)
		and spin-1 (blue curve) for the
			coupled $(1/2,1)$ system with parameters set as:
			$B_{1}=5,B_{2}=3,T_{1}=4,T_{2}=2$. Global
			work is the sum of local contributions by
			the two spins. As $J$ value is increased,
			spin-1/2 ceases to yield work first.
			Spin-1 is able to extract work even at
			higher $J$ values. These features
			can be analytically inferred, as described in Section VI
			}
		\label{fig_glw}
	\end{figure}
	%
	As discussed earlier,
	for $B_1 > B_2$ and $T_1 > T_2$,
	the majorisation relation ($ P \prec P'$) suggests
	the following set of sufficient conditions:  $B_2/T_2  \geq B_1/T_1$
	and $0 \leq J \leq J_c$. From the majorisation
	inequalities, we have derived, in Section V.A,
	the condition ${\cal X} \geq 0$ or PWC
	for the global system.
	Furthermore, it can also be shown that sufficient conditions for PWC in spin-1/2 and spin-1 are respectively (see Appendix C.1),
	\begin{equation} \label{eqnJmx1}
		0 \leq J\leq\frac{\Phi}{6} +\frac{1}{12}\left(\frac{1}{T_2}-\frac{1}{T_1}\right)^{-1}\ln\left(\frac{2+e^{-2B_1/T_1}}{2+e^{-2B_2/T_2}}\right)\equiv J_c^{(1/2)}.
	\end{equation}
	\begin{equation} \label{eqnJmx2}
		0 \leq J\leq\frac{\Phi}{6} +\frac{1}{12}\left(\frac{1}{T_2}-\frac{1}{T_1}\right)^{-1}\ln\left(\frac{1+5e^{-2B_1/T_1}}{1+5e^{-2B_2/T_2}}\right)\equiv J_c^{(1)}.
	\end{equation}

	Now, let us analyze the behavior of
	local work in this range of $J$ values.
	From the analytic expressions as well as
	Fig. (\ref{fig_jmax}), we notice that
	$J_c^{(1/2)} \approx J_c^{(1)} \approx  J_c \approx  \Phi/6$ at low temperatures. It implies that both local work and the global work are positive in $0\leq J\leq \Phi/6$. For higher temperatures,
	when $J\geq \Phi/6$,
	 the critical values  of $J$ follow: $J_c^{(1/2)}\leq J_c\leq J_c^{(1)}$, which indicates that the work performed by spin-1/2 can be negative even when the  global work and the
	 local work due to spin-1 are positive.
	The same conclusion can be justified from the sign of $\cal Z$.
	Thus, it can be proved that
	${\cal Z}\leq 0$ for $J > \Phi/6$
	(see Appendix C.2).
	Applied to Eq. (\ref{eqn105}), this result implies
	that the local work by spin-1 is {\it always}  positive
	provided the global work is positive (${\cal X} >0$) for
	$J > \Phi/6$.
	Analogously,
	it can be inferred from Eq. (\ref{eqn91})
	that for values of $J > \Phi/6$, the work performed by spin-1/2
	{\it can} be negative even when the global work is positive.

	Thus, we can summarize that in the range
	$0\leq J \leq \Phi/6$,
	both spin-1/2 and spin-1 yield positive work, and so the
	global work is positive. Beyond this range ($J > \Phi/6$),
	spin-1/2 may yield negative work even when the global work
	is positive. However,
	spin-1 yields positive work only if global work is positive.
	The local work analysis can be extended to the general case
	of coupled $(1/2,s)$ which is described in Appendix C.1.
	\begin{figure}[h!]
		\centering
		{\includegraphics[width=0.7\linewidth]{JC1.jpeg}}
		\caption{ The critical $J$ values for which
		majorisation inequalities hold for the
		coupled $(1/2,1)$ system, plotted as
		function of the hot temperature $T_1$.
		$J_c$ values shown are for the global system  (black, dashed curve)
		and spin-1 (blue curve) and spin-1/2 (red, dot-dashed curve)
        with parameters set as: $B_{1}=5,B_{2}=3$, $T_2 /T_1 = 0.5$.
        At low temperatures, global as well as local systems
        acquire a value $J_c = \Phi/6=0.167$.
        }
		\label{fig_jmax}
	\end{figure}

	\section{Enhancement of Otto efficiency}

	\begin{figure}[h!]
		\centering
		{\includegraphics[width=0.7\linewidth]{EF42.jpeg}}
		\caption{Quantum Otto efficiency for
			the coupled $(1/2,1)$ system versus the coupling
			strength $J$ (red, dashed curve).
			The parameters are set at $B_{1}=5,B_{2}=3,T_{1}=6,T_{2}=3$. The horizontal line is
			uncoupled efficiency, $\eta_0 = 1-B_2/B_1$,
			and Carnot bound is equal to 0.5. The
			upper bound, $\eta_{\rm ub}^{}$, is shown as the blue,
			dot-dashed curve.}
		\label{fig_eff}
	\end{figure}
	%
	The quantum feature of exchange coupling
	as a resource to enhance the Otto efficiency
	has been studied in earlier works.
	In Ref. \cite{PhysRevE.83.031135}, an upper bound for
	Otto efficiency which is tighter than the
	Carnot limit was derived. In the following,
	we revisit the feature of efficiency
	enhancement and provide justification
	for the upper bound based on the majorisation
	approach.

	The efficiency of the QOE, $\eta={|W|}/{Q_{1}}$
	can be expressed using Eqs. (\ref{eqn56}) and (\ref{eqn58}) as:

	%
	\begin{equation} \label{eqn73}
		\eta= {\eta_{0}}
		\left({1-\frac{6J{\cal Y}}{B_{1}{\cal X}}}
		\right)^{-1},
	\end{equation}
	where $\eta_{0}=1-{B_{2}}/{B_{1}}$
	is the efficiency of uncoupled system
	which is the same
	regardless of the magnitudes of individual spin.
	We have seen that owing to
	consistency with the second law,
	the quantum Otto efficiency is bounded from above by
	the Carnot value, $1-T_2/T_1$.
	However, heuristics can be applied to this system \cite{2021}
	to obtain a tighter upper bound on the efficiency,
	which is given by:
	\begin{equation} \label{eqn79}
		\eta_{\rm ub}^{} = {\eta_{0}^{}}
		\left({1-\frac{6J}{B_{1}}} \right)^{-1} \leq 1-\frac{T_2}{T_1}.
	\end{equation}
	Thus it can be seen that in the range
	$0\leq J \leq \Phi/6$, the above bound
	is lower than the Carnot value, and
	is saturated for $J = \Phi/6$, as shown in Fig. \ref{fig_eff}.

	Now, we show how the majorisation relation
	reveals this upper bound for Otto efficiency.
	We have seen ${\cal X}$ is positive
	in the range $0 \leq J\leq J_c$.
	The expression for $\eta$ suggests
	that in the presence of coupling ($J>0$),
	the condition ${\cal Y} = (P_{2}-P_{2}^{\prime})+(P_{4}-P_{4}^{\prime})>0$ enhances the efficiency
	over the uncoupled model ($\eta > \eta_0$).
	Now, combining $\mathcal{Y} >0$  with global majorisation condition (Eq. (\ref{eqn66})), we can write:
	\begin{equation}
		\frac{P_{2}+P_{4}}{P_{1}}\geq \frac{P_{2}^{\prime}+P_{4}^{\prime}}{P_{1}^{\prime}},
	\end{equation}
	implying that
 	\begin{align*}
		\frac{\exp{\big (\frac{B_{1}+12J}{T_{1}}\big)}+\exp{\big (\frac{-B_{1}+12J}{T_{1}}\big)}}{\exp{\big (\frac{3B_{1}}{T_{1}}\big)}}&\geq\frac{\exp{\big (\frac{B_{2}+12J}{T_{1}}\big)}+\exp{\big (\frac{-B_{2}+12J}{T_{2}}\big)}}{\exp{\big (\frac{3B_{2}}{T_{2}}\big)}},
	\end{align*}
which is satisfied in the range
$0\leq J\leq J_c$, indicating
the efficiency enhancement region.
The Otto efficiency is depicted in Fig. \ref{fig_eff}.

	Now, consider the expression,
	$({\cal X} - {\cal Y}) \equiv
	3(P_{6}-P_{6}^{\prime})+2(P_{5}-P_{5}^{\prime})+(P_{4}-P_{4}^{\prime})+(P_{3}-P_
	{3}^{\prime})$. Adding the inequalities (\ref{eqn60}), (\ref{eqn61})
	and (\ref{eqn63}), and upon rearranging terms,
	we obtain ${\cal X} - {\cal Y} \geq 0$, or
	${\cal X} \geq {\cal Y}$. Then, from Eq. (\ref{eqn73}),
	this implies that in the domain of the majorisation relation
	($P \prec P'$),
	the efficiency has the  upper bound given by $\eta_{\rm ub}$.
	The proof can be extended to the general case of $(1/2,s)$
	system, which is discussed in Appendix D.

\section{Conclusions}
QOE is one of the most well studied models of a quantum heat engine. It is based on generalizations of the classical adiabatic and  isochoric processes. Further, a spin-based working medium provides a convenient platform to investigate quantum features and their advantages
for a QOE.
Thus a QOE based on spin-1/2 particle is well known to have an efficiency, $\eta_0=1-B_2/B_1$ \cite{kieu2004second}. Also, $B_2/T_2>B_1/T_1$
is a necessary condition, given that $T_1>T_2$ and $B_1>B_2$. The additional condition guarantees PWC, which can be seen from
the expression for work extracted in a QOE:
\begin{equation}
	\vert W \vert=2(B_1-B_2)\left( \frac{1}{1+e^{2B_1/T_1}}-\frac{1}{1+e^{2B_2/T_2}}\right) \geq 0.
\end{equation}
We have shown that PWC for a QOE based on an arbitrary spin can be derived from
the concept of majorisation between $P$ and $P'$, whereby the relation $P\prec P'$ provides a necessary and sufficient condition for PWC. Further, we have expressed the work output in a QOE in terms of the relative entropies $D(P||P')$ and $D(P'||P)$. It is clarified
in general that $S(P)>S(P')$ is necessary, but not a sufficient condition for PWC. Further, the total entropy generated in a quantum Otto
cycle is given by the sum total of the two relative entropies.

Then, we have considered a spin-1/2 interacting with a spin-$s$ via 1-d Heisenberg exchange interaction with isotropic coupling strength $(J>0)$. In this case, majorisation yield PWC, provided we impose $B_2/T_2>B_1/T_1$ and additionally restrict $0\leq J\leq J_c$. The critical value  $J_c$ provides a sufficient range for parameter $J$ such that the majorisation inqualities hold good. We have treated the $s=1$ case in detail and provided expressions for the case with a general $s$ value, in the Appendix.
It is important to remark that in Ref. \cite{2021}, an extreme case scenario was used to infer the permissible range of $J$ as $0\leq J\leq \Phi/6$. Since, $J_c \geq \Phi/6$, so the range inferred in this paper extends the previous range of Ref. \cite{2021}.

Using the global majorisation inequalities, we have also investigated the local thermodynamics of spins. Thereby, it is possible to infer that spin-1/2 ceases to yield work at a certain $J$ value, while spin-1 continues to output work at larger $J$ values. The global work vanishes in between these two values, giving $J_c^{(1/2)}\leq J_c\leq J_c^{(1)}$.
Again, in previous works \cite{PhysRevE.83.031135, 2021}, an enhancement of efficiency was reported for the coupled model and an upper bound for the Otto efficiency was inferred which is tighter than the Carnot value. In the present work, we have justified this bound using majorisation inequalities. All these results can be extended to a (1/2,$s$) system.

In conclusion, our analysis shows that majorisation relation
can usefully characterize the operational conditions for
a quantum Otto engine based on spins as the working medium.
The approach based on majorisation is able to highlight key qualitative features of  even the local spin cycles.
 It will be interesting to extend the analysis to other interacting models  \cite{zhang2008entangled, azimi2014quantum, Halpern2019}.  Thus,
 majorisation is expected to serve as a key heuristic in inferring the thermodynamic features of complex working media in a quantum heat engine.
	\section*{Acknowledgment}
	SS acknowledges financial support in the form
	of Senior Research Fellowship from the Council for Scientific and
	Industrial Research (CSIR) via Award No. 09/947(0250)/2020-EMR-I India.