\section{Instructions for RA-M}

	\section{Introduction}
	\IEEEPARstart{U}{ntil} now, object rigidity is one of the common assumptions in robotic grasping and manipulation. Strictly speaking, all objects deform upon force interaction. Rigidity is a valid assumption when object deformation can be neglected in the task. Nevertheless, many objects that need to be manipulated by robots present non-negligible deformation: from micro surgical operation to challenging industrial assemblies.

	Robots need to be capable of manipulating deformable objects to operate in human environments. This capability would benefit many application fields, while posing fundamental research challenges. In this article, we consider a generalized concept of manipulation where grasping is also part of the task. We will refer to the problem as deformable object manipulation (DOM).

	The tasks involved in DOM cover a broad spectrum (see Fig. \ref{fig:real-world-applications}). These include: dressing assistance in elderly care, cable harnessing in industrial automation, harvesting and processing fruit and vegetables in agriculture, surgical operations in medical services, to name a few.
	\begin{figure}[t]
		\centering
		\includegraphics[width=0.48\textwidth]{real-world-applications.pdf}
		\caption{Applications involving manipulation of deformable objects. Clockwise from top left: dressing assistance~\cite{assist_dressing}, cable harnessing~\cite{cable}, fruit harvesting~\cite{fruit_harvesting}, suturing~\cite{suture}}
		\label{fig:real-world-applications}
	\end{figure}

	On the technical side, addressing deformation introduces the following technical challenges:
	\begin{itemize}
		\item the complication of sensing deformation,
		\item the high number of degrees of freedom of soft bodies,
		\item the complexity of non-linearity in modeling deformation.
	\end{itemize}

	We believe that overcoming these challenges is not only beneficial to DOM, but can further push towards developing autonomous robots which can operate in unstructured environments.
	\begin{figure*}[t!]
		\centering
		\includegraphics[width=0.6\textwidth]{system-level-overview-3.pdf}
		\caption{A typical robotic framework for handling deformable objects. In this particular example, the framework addresses wire harness~\cite{zhu2019robotic}.}
		\label{fig:system-level-overview}
	\end{figure*}
	In recent years, there have been a few surveys on robotic manipulation of deformable objects.
	Some surveys focus on specific areas of DOM.
	The survey from Jimenez \cite{jimenez2012survey} focuses on model-based manipulation planning. More recently, Herguedas et al. \cite{herguedas2019survey} review works using multi-robot systems for DOM while the work of \cite{nadon2018multi} considers multi-modal sensing. The authors of \cite{arriola2020modeling} present the state-of-the-art on deformable object modeling for manipulation.
	There are also two comprehensive surveys in the area.
	The survey in \cite{sanchez2018robotic} reviews and classifies the state-of-the-art according to the object's physical properties. Lately, \cite{yin2021modeling} reports most recent advances in modeling, learning, perception, and control in DOM.

	In contrast with the mentioned surveys, which either focus on reporting the progress of the field or a specific area, this article aims at identifying \emph{scientific challenges} introduced by object deformations and at projecting crucial future research directions. As DOM is an emerging field of research where there is still much to be done, in this paper, previous works and open problems are given equal weights. In addition, we dedicate one section to discussing practical challenges in various applications of DOM.
	We believe the paper is a first of its kind, in the field of DOM.

	A robotic framework designed to handle deformable objects usually consists of five key components: \emph{gripper and robot design}, \emph{sensing}, \emph{modeling}, \emph{planning} and \emph{control} (See Fig. \ref{fig:system-level-overview}). To position the current research and identify future trends, we conducted a survey on the future perspective of deformable object manipulation\footnote{Link to the survey: \url{https://forms.gle/XCv2CV79yvRP5Gsd7}}. We shared the survey with people working in related field, at various career stages. They were asked to rate the importance and research maturity of each of the five identified key components, from $1$ to $4$, with $1$ being not important/low maturity and $4$ being very important/high maturity. We received 31 answers that are summarized on Fig. \ref{fig:maturity_importance}.

	\begin{figure}[h!]
		\centering
		\subfloat[Highest qualifications of the respondents.]{\includegraphics[width=0.3\textwidth]{qualification.png}\label{fig:qualifications}}
		\hspace{2mm}
		\subfloat[Means and variances of Importance and research maturity ratings of each key component.]{\includegraphics[width=0.5\textwidth]{maturity_important.png}}
		\caption{Summary of the outcomes of the survey on DOM. We received in total $31$ answers. The respondents cover different level of qualifications ranging from master students to full professors.}
		\label{fig:maturity_importance}
	\end{figure}

	We consider the promising direction of research as the ones that have the highest significance and the lowest research maturity. Based on the survey, sensing is the most promising one among all subareas. This is probably due to current booming trend in Deep Learning has offered many new methods for processing the sensory data. In addition, sensing is the prerequisite for subsequent steps such as modeling, planning and control.

	Accordingly, the following sections of this paper each present one of these five research directions. In each section, we review recent works in the field and then comment the outlook and challenges ahead. Then, Sect. \ref{sec:application} tries to provide a link from research to practical applications in the context of DOM.  Finally, we summarize key messages in Sect. \ref{sec:summary}.

	\section{Gripper and robot design}\label{sec:hardware}
	\subsection{Current capability}\label{sec:hardware_current}

	Does manipulation of deformable objects demand specific grippers as compared to manipulation of rigid objects? Generally, yes (see Fig. \ref{fig:hardware}). Unlike rigid objects (which are mostly handled by standard grippers), deformable objects are handled with custom (and often designed ad-hoc) grippers, e.g. a 3D printed gripper that enables cable sliding \cite{zhu2019robotic}, a flat clip for holding towels \cite{hu20193}, a cylindrical tool for pushing and tapping plastic materials \cite{cherubini2020model}, a soft hand for manipulating organs \cite{liu2020musha}.
	Such diversity in grippers is a result of the large variety of deformable objects, which require different actions during manipulation. To avoid designing task-specific grippers for DOM, human-like dexterity and compliance is desired. Recent works in this directions consider compliant design \cite{della2015dexterity, abondance2020dexterous} and show good potential for DOM tasks.

	As for the robot itself, it is rigid in most works. In some cases, as in the surgical application showcased in~\cite{alambeigi2018autonomous} (Fig. \ref{fig:hardware}, bottom right), both robot and object are deformable to ensure safety of manipulation.

	\begin{figure}[!thpb]
		\centering
		\includegraphics[width=0.48\textwidth]{hardware.pdf}
		\caption{Various robot grippers for DOM. Clockwise from top left: a tool for pushing and tapping on plastic materials \cite{cherubini2020model}, flat clips for holding a towel \cite{hu20193}, a gripper allowing a cable to slide \cite{zhu2019robotic}, a soft continuum manipulator interacting with a deformable material  \cite{alambeigi2018autonomous} and a soft hand for manipulating organs \cite{liu2020musha}.}
		\label{fig:hardware}
	\end{figure}
	\subsection{Challenges and outlook}\label{sec:hardware_outlook}

	Improving dexterity is core to robot manipulation. The improvement can come from different research domains, such as accurate in-hand sensing or robust control, two aspects which we will detail in Sections \ref{sec:sensing} and \ref{sec:control}, respectively. In this section, our focus is on gripper/robot hardware aspects.

	One way of achieving such dexterity is to reproduce by design the most dexterous gripper -- the human hand. An open question is whether anthropomorphic design is in itself the optimal solution in all cases, especially in the context of DOM.

	While having one dexterous gripper which can handle a variety of DOM tasks is appealing, it should be noted that additional constraints need to be considered in the design process, for hygiene/safety in tasks such as food handling or surgery. For instance, for surgical applications we are limited by biocompatibility of the materials and actuators and by the reduced available space in minimally invasive surgery. In theses cases, designing task-specific grippers is more appropriate. Non-anthropomorphic soft grippers are another emerging area of research \cite{Hao2021}. These grippers are promising, to overcome the challenges associated with traditional fingered grippers in grasping rigid objects; yet, to date, their application to DOM receives little attention.

	Otherwise, one may use a standard gripper, and provide the robot with suitable tools to be grasped and used according to the type of task at stake. This demands breakthroughs on the algorithmic side, to make the robot capable of reasoning on the proper tools for different tasks. Training the robot to have task-specific reasoning will enhance autonomy and make robots realize more complex tasks.

	\begin{figure}[b]
		\centering
		\subfloat[Picking raspberries~ \cite{strawberry-picking}.]{\includegraphics[width=0.215\textwidth]{berry-picking.png}\label{fig:strawberry}} \hspace{2mm}
		\subfloat[A custom 3D printed soft robotic gripper, grasping mushroom coral~\cite{Vogt2018-po}]{\includegraphics[width=0.215\textwidth]{mushroomCoral.png}\label{fig:coral_reef}}
		\caption{Two examples of interaction with fragile objects, which could benefit from the use of soft robots.}
		\label{fig:soft-robots}
	\end{figure}

	Another area worth investigating is that of soft robots/grippers, since these have great potential for manipulating fragile materials, such as organs or food, or for collecting biological samples or fruits (see Fig. \ref{fig:soft-robots}). While traditional rigid robots need to exhibit compliant behavior when interacting with these objects, the inherent compliance of soft robots makes the task safe. This unconventional paradigm of using soft robots to manipulate soft objects will bring new challenges in modeling and control as both the robot and the object are under-actuated and difficult to model. One pioneering work in this direction is \cite{Ficuciello2018}, which adapts the finite element modeling (FEM) based inverse soft robot model with contact handling (proposed in \cite{Coevoet2017}) for deformable objects manipulation using soft robots.

	An interesting research question to consider is whether methods can be transferred from one field to the other. To be more specific, can methods for controlling/modeling soft robots be applied to manipulating deformable objects and vice versa? If so, as a community, it may be valuable to obtain a unified approach for working with both soft robots and deformable objects.

	\section{Sensing}\label{sec:sensing}

	\subsection{Current capability}\label{sec:sensing_current}

	In this section, we consider \textit{visual}, \textit{tactile} and \textit{force} sensing for DOM. Existing research relies on these three modes to estimate the state of deformable objects. In most cases, vision provides global information about shapes on a large scale, while force and tactile provide local information on both shape and contact. At the end of this section, we also discuss the research in contrast with this common practice, where global deformation properties are recovered using tactile sensing. It should also be noted that force information is particularly important in industrial settings, e.g., for assembly \cite{luo2018deep}, \cite{hayami2019multi}.

	Vision is used in tasks such as rope manipulation \cite{nair2017combining,yan2020self} or cloth unfolding \cite{sun2015accurate,li2015regrasping}, where the object exhibits large global deformation. In these works, configurations of deformable objects were obtained from raw image readings. Although vision offers a global perspective of the object configuration, visual data can be noisy in unstructured environments, it is then important to manage occlusions~\cite{hu20193,cheng2019occlusion}.
	Most of above-mentioned works are based on 2D vision, 3D perception of deformable objects is more challenging. Existing works employ FEM~\cite{petit2017tracking} or a combination of Growing Neural Gas and Particle Graph Networks \cite{valencia2019toward} for better tracking the deformation. In a more recent study~\cite{rouhafzay2021transfer}, it has been shown that a deep convolutional neural network for processing vision data can be used with small variations to process tactile data for deformable objects recognition.

	Objects made of soft materials, such as human tissues and fruits, have a special force-displacement correlation upon contact. As a result, tactile sensing can be used to estimate the stiffness. In~\cite{yuan2016estimating}, the GelSight~\cite{yuan2017gelsight}, a vision-based high-resolution tactile sensor, measures the 3D geometry of the contact surface, and the normal/shear forces.

	Note that the division of \textit{vision for global deformation} and \textit{tactile sensing for local deformation} is not absolute. The authors of~\cite{guler2014s} present to use vision to estimate the local deformation of objects during grasping, and classify objects accordingly. In~\cite{yuan2018active}, high-resolution tactile sensing is used to estimate the physical properties of clothing materials through squeezing, assuming the robot can learn from the data about indicating global properties of clothing according to a local sampling point. In ~\cite{she2019cable} an example of servoing along a cable based on high-resolution tactile sensing is presented. Although vision is not used, the precise measurement of the local cable shape provides enough information to guide the robot motion on a small scale.

	\subsection{Challenges and outlook}\label{sec:sensing_chal_n_out}
	Here, the main challenges are: selecting appropriate sensors for the DOM task and using the measurements to obtain meaningful object representations.

	Considering the high number of degrees of freedom (DoF) of the deformable bodies, fusion of different sensing modalities (vision, force and tactile) may be a promising direction to pursue in future research.

	Another research question to be answered is: what yields a good representation of the object configuration? We (acknowledgedly) do not have a complete answer to this; rather, we will elaborate on considerations when designing the representation.

	The representation need to be robust to noise and useful for reconstructing the objects' configuration -- even when data are partially unavailable. In vision, the most common noise is occlusion. How to generate a meaningful representation of these objects under self occlusion is still an open problem in research.
	For rigid objects, one can carefully design the environment to avoid it. For deformable objects that exhibit large global deformation such as clothes, bed sheets etc, self occlusion is inevitable during manipulation. A promising direction to deal with occlusion and noise is using active/interactive perception. With vision data from different perspectives, we might be able to reconstruct the object's configuration accurately even under occlusion and noise.

	Apart from the above mentioned challenges, choosing a good representation also involves leveraging two aspects:
	\begin{enumerate}
		\item the dimensionality of the representation,
		\item the accuracy of the representation.
	\end{enumerate}
	Usually the trade-off depends on the task, relies on human intuition and involves a trial and error process.
	In end-to-end reinforcement learning settings, sensory data can be mapped directly to robot actions without explicit feature representations \cite{levine2016end}. Human demonstrations can be used for making end-to-end learning more efficient. One example is reported in \cite{matas2018sim}. The authors use an improved version of Deep Deterministic Policy Gradients, trained with $20$ demonstrations, to make robots manipulate cloth. However, since such settings often require a manually designed cost/reward function for learning, human demonstrations in this context can also be used for recovering the reward, with inverse reinforcement learning.

	\section{Modeling}\label{sec:current-model}
	\subsection{Current capability}
	For robots to perform deformation tasks using sensory data, we need a model that captures the relationship between sensor information and robot motion.
	A linear model characterized by Young's modulus can be employed for describing elastic deformation. The two other classes of deformation are: plastic, and elasto-plastic deformations. This classification serves well. Yet, since the model should be used for control, in this section, we prefer to distinguish between local and global models --  a taxonomy which has clearer implications for control. We introduce the corresponding research  and  -- at the end of the section -- we discuss the limitations of these models and present works that address them.

	Most local models approximate the perception/action relationship via a Jacobian Matrix (called \textit{Interaction Matrix} in visual servoing). Such a model is linear and can be computed in real-time with a small amount of data. Yet, since it is a local model, it should be continuously updated during task execution. Model updating methods include: Broyden rule \cite{alambeigi2018autonomous}, receding horizon adaption \cite{zhu2021vision}, local gradient descent \cite{navarro2018fourier}, QP-based optimization \cite{lagneau2020automatic}, and Multi-armed Bandit-based methods \cite{McConachie2016mab}. Another advantage of the Jacobian model is that one can design a simple controller by inverting it. However, since this controller is local, it should operate via a series of intermediate target shapes \cite{zhu2021vision}, \cite{lagneau2020automatic}.

	On the other hand, global models can be approximated with Finite Element Methods \cite{yoshida2015} and also (deep) neural networks. In contrast to simple linear models, (D)NN-based approaches benefit from stronger representation power, in terms of accuracy and robustness~\cite{valencia2020combining}. Moreover, they can incorporate physics models and reason about object interaction \cite{battaglia2016interaction}. These models can approximate highly nonlinear systems and have a larger validity range, solving (to some extent) the locality issue of the linear models. Nevertheless, these complex nonlinear representations demand large amounts of data (which might not be available in some cases).

	Yet, whether we use analytical or learned models, their predictive power will be limited. They are either specialized to some class of tasks or learned from a set of training data. Especially for the learned models, we can never hope to collect enough data to produce an accurate model in the entire state space (which is high dimensional). Thus \cite{McConachie2020} and \cite{Power2021} have developed methods to reason about the validity of a (learned) model for a given state and action, and have used these methods to reason about model uncertainty in planning and control. However, when the model is not precise, a re-planning/recovery might be desirable. The authors of \cite{mitrano2021learning} introduces two neural networks for learning and re-planning the motion when the model is unreliable.
	\subsection{Challenges and outlook}\label{sec:model-outlook}
	The complexity of modeling is manifested in the lack of simulators. While most existing robotic simulators are capable of producing rigid body kinematics and dynamic behaviours, only a fraction of them can handle deformation. One recent work, Softgym \cite{lin2020softgym} was proposed for bench-marking DOM based on Nvidia Flex. In the soft robotics community, SOFA \cite{coevoet2017software} and Chainqueen \cite{hu2019chainqueen} are example simulators. In Sect. \ref{sec:hardware_outlook}, we considered the interaction between soft robots and deformable objects. Thus, a unified simulator that is able to handle soft robots and objects, and model their interaction might be desirable.

	When choosing a model for control, one challenge of data-driven deformation modeling is to balance region of validity with number of data required for training. One possible direction is to combine a simple model with a complex nonlinear model to form a hierarchical model. An example of such structures is exploited in \cite{li2020learning} for robust in-hand manipulation.

	For DOM tasks, we can have a linear model at lower level, and a (D)NNs learning the full model at higher level.
	The lower level model can be learned in few iterations to enable instant interaction between robot and object. The higher level (D)NN can collect data and improve the model to enhance global convergence.

	\section{Planning}
	\subsection{Current capability}
	Planning aims at finding a sequence of valid (robot/object) configurations and contributes to solving the problem of limited validity of local models, as discussed in Sect. \ref{sec:current-model}.

	Planners can operate in the objects' configuration space, and sometimes rely heavily on physic-based simulation. While the obtained plan can be visually plausible, it may be unrealizable for a specific object. Recently, McConachie et. al. presented a framework which combines global planning without physics simulation, with local control~\cite{mcconachie2020manipulating}. For an elastic object, considering its energy is another way to do planning; in this direction, Ramirez-Alpizar et al. \cite{ramirez2014} proposed a dual-arm manipulation planner optimizing the elastic energy, for elastic ring-shaped objects manipulation.
	For DOM tasks involving multiple robots, planning is important for coordination. Alonso-Mora et al. employed a distributed receding horizon planner for transporting tasks that require multiple robots \cite{alonso2015local}. More recently, \cite{9410363} learns a latent representation for semantic soft object manipulation that enables (quasi) shape planning with deformable objects.

	With Learning from Demonstration (LfD), the robot can be trained to manipulate deformable objects by an expert (usually a human). LfD encodes the robot trajectory and interaction force from human demonstrations \cite{lee2015learning_1}, thus avoiding explicitly planning the motion. More recently, Wu et al. have proposed a reinforcement learning scheme for DOM, which does not require initial demonstrations \cite{wu2019learning}.
	\subsection{Challenges and outlook}
	A rigid object configuration can be described in space with 6 DoF, whereas a deformable object configuration has a much higher number of DoF. To address this from the sensing algorithm side, one can find a compact representation from sensory data, as discussed in Sect. \ref{sec:sensing_chal_n_out}. An alternative, which receives much less attention, is the use of environmental contacts to constrain some DoF of deformable objects. Examples include the use of contact points in cable harness or that of flat surfaces when folding clothes. We argue that instead of planning to avoid contacts as most planners do, for deformable objects, we need to plan to make contact, since this constrains the configuration, and therefore simplifies the task.

	Planning to grasp the correct point is often crucial in DOM tasks. For instance, grasping at convex vertices of the clothes guarantees stability and facilitates the task~\cite{van2010gravity}. Re-grasp planning is highly relevant when considering tasks which require multiple robotic arms. Additional challenges come from perception, since as soon as the robot releases one or more grasp(s), the object is likely to change its configuration. We rely on sensing to track configuration changes and then plan accordingly.

	Another important future work in planning is reasoning about a deformable object at a semantic level. What does it mean for a cloth to be \textit{folded}? What does it mean for an object to be \textit{wrapped} in paper? We cannot manually specify all the configurations of the deformable object to use as goals in these kinds of tasks. Instead, we need a way to learn the meaning of semantic concepts, such as \textit{folded} or \textit{wrapped}, so that we can determine if a given configuration of the object is a valid goal.

	\section{Control}\label{sec:control}

	\subsection{Current capability}

	Control aims at designing inputs for the robot to realize the planned motion. The type of controllers is decided usually by the task. For instance, the authors employed a data-driven model predictive control \cite{9035011} for cutting considering its predictive nature and the lower demand for manual tuning. For safe interaction in minimally invasive surgery, the authors of \cite{su2019improved} used a fuzzy compensator with impedance control. For controlling large deformation, Aranda et al., proposed a Shape-from-Template algorithm concerning its low dimensional representation (using the template) and robustness against occlusion \cite{aranda2020monocular}.

	A number of works focus on shape control.

	While global models directly map sensor data to robot motion, local models must be inverted to design the robot motion controller (see Sec.~\ref{sec:current-model}).
	Several applications of the control scheme for robotic manipulation of deformable objects can be found in 3C manufacturing \cite{li2018pcb,li2018usb}, where vision-based controllers were proposed to drive the robot to automatically grasp/contact the deformable object, then carry out the task of active deformation or separation/sorting.
	Other works consider the concept of diminishing rigidity to do deformation control \cite{nadon2020grasp}, \cite{ruan2018accounting}.

	\subsection{Challenges and outlook}

	Feedback control has been commonly used in most DOM works, by referring to the state of the object, to achieve the task. Note that such state is retrieved from the output of its deformation model and measured with sensors, and that output and state do not necessarily have the same representation and dimension.
	Furthermore, we can distinguish between model-based and model-free control. Due to the complexity of modeling the deformation, when using the model to derive control policies, the controller has to take into account that the model will be inaccurate or even wrong.

	Model-free approaches do not require information about the deformation parameters or the structure of the deformation model.
	Examples include LfD or (Deep) reinforcement learning, where the challenges are: efficient use of data, and policy generalization. To address these issues, we can combine the offline and online learning methods. In the offline phase, the supervised network can be trained to estimate the model, by collecting pairs of a series of predefined inputs (e.g., the velocity of the robot end-effector) and the deformation of the object. The estimated model in the offline phase can be further updated online during the control task with adaption techniques (e.g., the adaptive NNs), to compensate the errors due to insufficient training in the offline phase or the changes of the deformation model. Hence, both complement each other.

	When multiple features on the deformable object are controlled in parallel, the system becomes under-actuated, with less control inputs than error outputs. Then, the robot controller should be able to deal with the conflicts between multiple features or decouple the control of multiple features in a sequential manner, to guarantee controllability.

	In addition, due to the deformation during control, the contact between robot end-effector and deformable object may not always be maintained. Most existing systems require a certain level of human assistance to initiate the contact or to re-establish it, if it is lost during the task. To improve autonomy, the robot controller should automatically grasp or touch the object first, whenever physical contact is lost, laying the foundation of the subsequent manipulation task. Such a capability would allow the robot to effectively deal with the unforeseen changes due to deformation.

	\section{Practical applications}\label{sec:application}
	In previous sections, we centered our discussions from a scientific point of view, here, we instead discuss challenges in various applications where DOM can be translated to solutions.

	\textbf{Automatic laundry:} A typical domestic application of DOM is laundry folding. A Tokyo-based company unveiled its prototype laundry-folding robot in 2015 (Fig. \ref{fig:applications}a). However, the company was announced bankrupt in 2019 due to lack of funding for development and difficulties in improving the robot to reach a satisfactory level \cite{nagata_2019}. Although cloth folding has been tackled in a few previous research \cite{verleysen2020video,miller2012geometric,tsurumine2019deep,borras2020grasping}, it remains largely a laboratory product (limited to structured environments, certain types of the clothes, etc). Commercializing the technology seems requiring a substantial efforts.

	\textbf{Assistive dressing:} Robotic dressing assistance has the potential to become an important technology due to the pressing needs for ageing society support.
	Research can roughly be categorized into simulation-based learning \cite{clegg2018learning,clegg2020learning} and imitation learning \cite{joshi2019framework} approaches. Examples are dressing support for shoes \cite{canal2019adapting}, shirts \cite{Li-RSS-21,tamei2011reinforcement,zhang2019probabilistic} and pants.
	However, several technical and societal challenges have to be addressed before robot-assisted dressing will become a broadly used DOM technology: physical safety for the human, modeling and prediction of the human-robot interaction, robustness for large variations of geometric and dynamic properties of textiles, low-cost high-reliable robot hardware, human acceptance of such technologies.

	\textbf{Surgical robotics:} Soft tissue manipulation is mainly performed with tele-operation solely using visual feedback. Autonomous manipulation, however, still has a long way to go and demands developing various DOM hardware and software (Fig. \ref{fig:applications}d). The biggest concern for an autonomous solution is the safety of operation. A soft robot with intrinsic compliance will probably enhance the safety.

	\textbf{Food production \& Retail:} Handling deformable objects is a major challenge in the whole chain from production to sales. In an agricultural setting, automated harvesting of fruits and vegetables requires interactions with deformable objects that are at the same time easy to damage, which immediately decreases their value and shelf live. Frequently, these  products also undergo an intermediate processing step (e.g., filleting and packaging meat). More generally, deformable products (e.g., everything packaged in flexible bags, (Fig. \ref{fig:applications}c)) need to be handled in warehouses, in order picking, and in restocking. Solutions for specific applications and products have been developed, but more complex objects and operations still are frequently handled by human workers.

	\textbf{Marine robotics:} Underwater grasping has been led by oil and gas industry for decades, resulting in heavy machines with strong grippers for inspection and maintenance tasks (Fig. \ref{fig:applications}e). Gradually the demands turned to more detailed tasks in marine biology, sedimentology and archaeology (Fig. \ref{fig:applications}f). Another DOM application can be found in tethered robot umbilical modeling and control. Negative buoyancy cable can be modeled in real time as a simple catenary shape and tracked to control a tethered ROV~\cite{LARANJEIRA2020107018}.

	\begin{figure}
		\centering
		\includegraphics[width=0.48\textwidth]{applications.pdf}
		\caption{Various applications of DOM -- (a): laundry-folding robot from Seven Dreamers Laboratories Inc. \cite{laundry}, (b): A mock-up for robotics dressing assistance, (c): a robot picking a flexible bag of goods on the shelf, courtesy AIRLab Delft \cite{AIRLab}, (d): autonomous surgical manipulation by the dVRK system \cite{alambeigi2018robust}, (e): ROV Victor 6000 sampling black smokers (IFREMER/GENAVIR) courtesy D. Desbruyères, (f):  Ultra soft underwater gripper for jellyfish \cite{Sinatra2019-sr}}.
		\label{fig:applications}
	\end{figure}

	\section{Summary and key messages}\label{sec:summary}

	The revolution of robots from automating repetitive tasks to humanizing robot behaviours is taking place with better hardware,  robust sensing capabilities, accurate modeling, increasingly versatile planning and advanced control. Manipulation of deformable objects breaks fundamental assumptions in robotics such as rigidity, known dynamics models and low dimensional state space. It therefore requires breakthroughs in all the areas mentioned above, and serves as a great test-bench for novel ideas in both robotic hardware and software. A summary of challenges and ideas discussed are presented in Fig. \ref{fig:DOM_challenges}.
	\begin{figure}
		\centering
		\includegraphics[width = 0.49\textwidth]{DOM_challenges.png}
		\caption{A summary of open research problems and ideas/methods to pursuit discussed in this paper. Research problems in each subarea are written with bold black texts, whereas ideas/methods to resolve them are marked with bold white texts in red eclipses.}
		\label{fig:DOM_challenges}
	\end{figure}

	In terms of hardware, recently, the community has been shifting more and more from rigid to soft robots. Robotic manipulation is also gradually shifting from rigid to deformable objects. One open question is if some of the algorithms in one field are transferable to the other? We believe the interaction between a soft robot and a deformable object will bring more challenges to the robotic community.

	Sensing plays a vital part in robotics manipulation of deformable objects. Depending on the nature and complexity of the task, one or multiple fused sensing modes may be needed. Machine learning will facilitate the development of robust algorithms to process data from different sensors, to generate meaningful representations of deformation.

	\emph{All models are wrong, some are useful}. We do not believe there exists the \say{best} model for deformation. While more and more models tend to be data-driven, we would like to draw the readers' attention to the importance of physical models for studying interactions.

	For planning, current research lacks a high level semantic reasoning of the DOM task. Furthermore, while often the purpose of planning is to avoid contact and collision, we argue that for DOM, it can be very useful to plan for contact.

	Under-actuation is a key challenge of DOM, due to the deformable bodies' high DoF. Another practical issue introduced with deformation is contact loss during manipulation; future controllers should be able to detect contact loss and to react accordingly.


\section{Introduction}
Recently, superconducting hybrid systems have drawn attention due to their possible device applications as sensitive detectors\cite{hei,gua,pso}, low-temperature sensitive thermometers\cite{zgi,wang,luu}, heat valves\cite{yan,soth,dutta}, solid-state quantum machines\cite{mar,vis,kar}, solid-state micro-refrigerators\cite{mml,ngu,bos} and thermoelectric generators\cite{gma,huss,bra}.
Quantum thermodynamics implies the study of thermodynamic processes from the principle of quantum mechanics\cite{rko}, while refrigeration means transfer of heat from low to high temperature region\cite{kos} aided by work done on system. In this context  thermodynamic\cite{mjm,est,for,afo} properties of a Josephson Stirling engine have been discussed in Ref.~\onlinecite{bsc}, wherein a quantum spin Hall insulator based Josephson junction is shown to act as a quantum heat engine. Josephson Stirling engines aren't the only game in town, diffusive SNS junctions have been shown to operate as a Josephson-Otto or as Josephson-Stirling engines in Ref.~\onlinecite{carr}.

In this paper, by doping a spin flipper in a 1D Josephson junction(JJ) loop which is in turn attached to two thermal reservoirs at in-equivalent temperatures via thermal valves, we show that this device can be employed both as a quantum heat engine as well as a refrigerator and can also work as a Joule pump or cold pump even.

The rest of the paper is arranged as follows. In section II the model is introduced via Hamiltonian, wave functions and boundary conditions so as to calculate different thermodynamic quantities. In section III the thermodynamic processes involved in the Josephson-Stirling cycle are discussed, then work done and heat exchanged during each of these processes is calculated. Next, in section IV the results are analyzed and different operating modes for Josephson-Stirling cycle are discussed. Our paper ends with an experimental realization of the proposed device and the take home messages are aggregated in section V.
\section{Theory}
The model device, depicted in Fig.~1, is formed from a 1D superconducting loop\cite{butt} interrupted by a spin flipper. An external magnetic flux $\Phi$ controls the superconducting phase difference across the spin-flipper. The JJ loop is attached, via two thermal valves $v_{L}$ and $v_{R}$, to reservoirs at either end which in turn are at temperatures $T_{L}$ and $T_{R}$. The two reservoirs can exchange heat $Q_{L}$ and $Q_{R}$ with the JJ loop. The scattering problem is solved using BTK approach\cite{BTK} for superconductor-spin flipper-superconductor junction as shown in dashed line box of Fig.~1. The two reservoirs control temperature of JJ loop via thermal valves $v_{L}$ and $v_{R}$. When valve $v_{R}$ is opened and $v_{L}$ is closed, JJ loop is in thermal contact with right reservoir at temperature $T_{R}$. Similarly, when valve $v_{L}$ is opened and $v_{R}$ is closed, JJ loop is in thermal contact with left reservoir at temperature $T_{L}$. On the other hand, phase difference across JJ loop is controlled via magnetic flux $\Phi$ enclosed by the loop. Thus, by controlling both temperature and phase difference, the JJ device can be driven from one state to another. We discuss this in more detail for a Stirling cycle in section III.
\begin{figure}[ht]
\centering{\includegraphics[width=0.5\textwidth]{spin-heat-engine.pdf}}
\caption{\small \sl 1D Josephson junction (JJ) loop (circumference $L_{S}$, in orange) doped with a spin flipper and attached to two thermal reservoirs at temperatures $T_{L}$ and $T_{R}$ via two thermal valves $v_{L}$ and $v_{R}$. A magnetic flux $\Phi$ controls phase difference $\varphi$ across spin flipper in JJ loop.}
\end{figure}
\subsection{Hamiltonian}
In our proposed set-up, a spin-flipper is embedded in the JJ loop of circumference $L_{S}$, see Fig.~1. We used BTK approach\cite{BTK} to solve the scattering problem. The spin-flipper is a point like magnetic impurity, see Ref.~\onlinecite{AJP}. Hamiltonian for spin-flipper, from Refs.~\onlinecite{AJP,Liu,Maru,FC,ysr} is,
\begin{equation}
H_{\mbox{spin-flipper}}=-J_{0}\vec{s}.\vec{S}.
\label{spinflipper}
\end{equation}
Bogoliubov-de Gennes (BdG) Hamiltonian for the JJ loop is given as
\begin{equation}
H_{BdG}(x)=
\begin{pmatrix}
H\hat{I} & i\Delta(x) \hat{\sigma}_{y}  \\
-i\Delta^{*}(x) \hat{\sigma}_{y}   &  -H\hat{I}
\end{pmatrix},
\label{ham}
\end{equation}
where $H=p^2/2m^{\star}-J_{0}\delta(x)\vec s.\vec S-E_{F}$, with $p^2/2m^{\star}$ represents kinetic energy of electron-like or hole-like quasiparticle of mass $m^{\star}$, and $E_{F}$ is Fermi energy. Superconducting gap $\Delta(x)$ is of the form $\Delta(x)=\Delta[e^{i\varphi_{L}}\Theta(-x)+e^{i\varphi_{R}}\Theta(x)]$. For simplicity, we consider superconducting gap to be independent of temperature, i.e., $\Delta(T)\approx \Delta(T=0)$ with $\frac{\partial\Delta}{\partial T}\approx0$, which is valid if $T\ll T_{c}$, $T_{c}$ being superconducting critical temperature. $\varphi_{L}$ and $\varphi_{R}$ are superconducting phases for left and right superconductors respectively as shown in dashed line box of Fig.~1. In the manuscript, we use dimensionless parameter $J=\frac{m^{\star}J_{0}}{\hbar^2k_{F}}$ as a measure of strength of exchange interaction\cite{AJP} between quasi-particles and spin-flipper.
\subsection{Wavefunctions and boundary conditions}
Diagonalizing BdG Hamiltonian (Eq.~\eqref{ham}) one gets the wavefunctions in superconducting regions of our system for electron/hole like quasiparticle incidence.
For electronlike quasiparticle with spin up incident from left superconductor, wave function for left superconductor is\cite{LINDER}-
\begin{eqnarray}
\psi_{S_{L}}(x)=&&\begin{pmatrix}u\\
                              0\\
                              0\\
                              v
                              \end{pmatrix}e^{iq_{+}x}\phi_{m'}^{S}+r_{ee}^{\uparrow\uparrow}\begin{pmatrix}
                              u\\
                              0\\
                              0\\
                              v
                             \end{pmatrix}e^{-iq_{+}x}\phi_{m'}^{S}\nonumber\\&&+r_{ee}^{\uparrow\downarrow}\begin{pmatrix}
                             0\\
                             u\\
                            -v\\
                             0
                             \end{pmatrix}e^{-iq_{+}x}\phi_{m'+1}^{S}+r_{eh}^{\uparrow\uparrow}\begin{pmatrix}
                             0\\
                            -v\\
                             u\\
                             0
                             \end{pmatrix}e^{iq_{-}x}\phi_{m'+1}^{S}\nonumber\\&&+r_{eh}^{\uparrow\downarrow}\begin{pmatrix}
                             v\\
                             0\\
                             0\\
                             u
                             \end{pmatrix}e^{iq_{-}x}\phi_{m'}^{S}, \mbox{ for $x<0$ }.
                             \label{lsi}
                             \end{eqnarray}
$r_{ee}^{\uparrow\uparrow},r_{ee}^{\uparrow\downarrow},r_{eh}^{\uparrow\uparrow},r_{eh}^{\uparrow\downarrow}$ are normal reflection amplitude without any flip, normal reflection amplitude with flip, Andreev reflection amplitude with flip and Andreev reflection amplitude without any flip respectively.\\
The corresponding wave function in right superconductor,
\begin{eqnarray}
\psi_{S_{R}}(x)=&&t_{ee}^{\uparrow\uparrow}\begin{pmatrix}
                              ue^{i\varphi}\\
                              0\\
                              0\\
                              v
                             \end{pmatrix}e^{iq_{+}x}\phi_{m'}^{S}+t_{ee}^{\uparrow\downarrow}\begin{pmatrix}
                             0\\
                             ue^{i\varphi}\\
                             -v\\
                             0
                             \end{pmatrix}e^{iq_{+}x}\phi_{m'+1}^{S}\nonumber\\&&+t_{eh}^{\uparrow\uparrow}\begin{pmatrix}
                             0\\
                             -ve^{i\varphi}\\
                             u\\
                             0
                             \end{pmatrix}e^{-iq_{-}x}\phi_{m'+1}^{S}+t_{eh}^{\uparrow\downarrow}\begin{pmatrix}
                             ve^{i\varphi}\\
                             0\\
                             0\\
                             u
                             \end{pmatrix}e^{-iq_{-}x}\nonumber\\&&\phi_{m'}^{S},\mbox{ for $x>0$ },
                             \label{rsi}
                             \end{eqnarray}
$t_{ee}^{\uparrow\uparrow},t_{ee}^{\uparrow\downarrow},t_{eh}^{\uparrow\uparrow},t_{eh}^{\uparrow\downarrow}$ being transmission amplitudes, corresponding to reflection processes described above and $\varphi=\varphi_{R}-\varphi_{L}$ is phase difference between right and left superconductors. $\phi_{m'}^{S}$ represents eigenspinor of spin flipper, with the $S^{z}$ operator of spin-flipper acting as, $S^{z}\phi_{m'}^{S} = m'\phi_{m'}^{S}$, $m'$ being spin magnetic moment of spin-flipper. BCS coherence factors are $u=\sqrt{\frac{1}{2}\Bigg(1+\frac{\sqrt{E^2-\Delta^{2}}}{E}\Bigg)}$, $v=\sqrt{\frac{1}{2}\Bigg(1-\frac{\sqrt{E^2-\Delta^{2}}}{E}\Bigg)}$. Wavevectors for electron-like quasiparticles ($q_{+}$) and hole-like quasiparticles ($q_{-}$) are
$q_{\pm}=\sqrt{\frac{2m^{\star}}{\hbar^2}(E_{F}\pm \sqrt{E^2-\Delta^2})}$. Andreev approximation\cite{Kri} gives $q_{+}=q_{-}=k_{F}$, with $k_{F}$ being Fermi wavevector, and $E_{F}>>\Delta$. Imposing boundary conditions on (\ref{lsi}, \ref{rsi}) at $x=0$, gives\\
\begin{equation}
\psi_{S_{L}}(x)=\psi_{S_{R}}(x),\,\, \frac{d\psi_{S_{R}}}{dx}-\frac{d\psi_{S_{L}}}{dx}=-\frac{2m^{\star}J_{0}\vec s.\vec S}{\hbar^2} \psi_{S_{L}},
\label{bc}
\end{equation}
where $\vec s.\vec S=s^{z}S^{z}+\frac{1}{2}(s^{-}S^{+}+s^{+}S^{-})$ represents exchange operator in Eq.~\eqref{spinflipper}, with
\begin{eqnarray}
&&\vec s.\vec S \psi_{S_{L}}(x=0)=\vec s.\vec S\begin{pmatrix}u\\
                              0\\
                              0\\
                              v
                              \end{pmatrix}\phi_{m'}^{S}+r_{ee}^{\uparrow\uparrow}\vec s.\vec S\begin{pmatrix}
                              u\\
                              0\\
                              0\\
                              v
                             \end{pmatrix}\phi_{m'}^{S}\nonumber\\&&+r_{ee}^{\uparrow\downarrow}\vec s.\vec S\begin{pmatrix}
                             0\\
                             u\\
                            -v\\
                             0
                             \end{pmatrix}\phi_{m'+1}^{S}+r_{eh}^{\uparrow\uparrow}\vec s.\vec S\begin{pmatrix}
                             0\\
                            -v\\
                             u\\
                             0
                             \end{pmatrix}\phi_{m'+1}^{S}+r_{eh}^{\uparrow\downarrow}\vec s.\vec S\begin{pmatrix}
                             v\\
                             0\\
                             0\\
                             u
                             \end{pmatrix}\phi_{m'}^{S}.\nonumber\\
\end{eqnarray}
$s^{\pm} = s_{x}\pm is_{y}$ and $S^{\pm} = S_{x}\pm iS_{y}$ are spin raising and spin lowering operators for electronlike quasiparticle/holelike quasiparticle and spin-flipper respectively. In Appendix A we give a detailed calculation as to how spin flipper interacts with electron-like and hole-like quasiparticles.
\subsection{Andreev bound states}
To calculate Andreev bound states (see Ref.~\cite{annu} for details of the method) we neglect contribution from incoming quasiparticle and insert wavefunctions (\ref{lsi}, \ref{rsi}) into boundary conditions (\ref{bc}). We get a homogeneous system of $8$ linear equations for the scattering amplitudes,
\begin{equation}
Py=0
\end{equation}
where $y$ is a $8\times1$ column matrix and given as $y=[r_{ee}^{\uparrow\uparrow},r_{ee}^{\uparrow\downarrow},r_{eh}^{\uparrow\uparrow},r_{eh}^{\uparrow\downarrow},t_{ee}^{\uparrow\uparrow},t_{ee}^{\uparrow\downarrow},t_{eh}^{\uparrow\uparrow},t_{eh}^{\uparrow\downarrow}]$, $P$ being a $8\times8$ matrix. For nontrivial solution of this system, the determinant of $P=0$ and we get Andreev bound states as a function of phase difference $\varphi$ between two superconductors, i.e., Andreev bound state energy spectrum $E_{j}$, $j=\{1,...,4\}$\cite{Been}. We find that $E_{j}(\varphi)=E_{\sigma}^{\pm}(\varphi)=\pm E_{\sigma}(\varphi), (\sigma=\uparrow,\downarrow)$ and
\begin{equation}
E_{\sigma}^{\pm}(\varphi)=\pm\Delta\sqrt{1+\frac{A(\varphi)}{C}+\rho_{\sigma}\frac{\sqrt{B(\varphi)}}{C}},\label{Eq:eb}
\end{equation}
\begin{eqnarray}
\mbox{ where, } A(\varphi)=&&J^2(2+p^4J^2+2p^2(-2+J^2m'(1+m'))\nonumber\\
                           &&+m'(1+m')(4+J^2m'(1+m')))\nonumber\\
                           &&+2(8+J^2(1-2p^2+2m'(1+m')))\cos(\varphi),\nonumber
\end{eqnarray}
\begin{eqnarray}
B(\varphi)=&&2J^2(64p^4J^2+3(J+2Jm')^2+4p^2\nonumber\\
                           &&(16+J^2(5+4m'(1+m')))+4J^2(-4p^2+16p^4-\nonumber\\
                           &&(1+2m')^2)\cos(\varphi)+((J+2Jm')^2-4p^2(16+\nonumber\\
                           &&(J+2Jm')^2))\cos(2\varphi)),\nonumber\\
C=(1&&6+J^4(p^2+m'+m'^2)^2+J^2(4+8p^2+8m'(1+m'))),\nonumber
\end{eqnarray}
and $\rho_{\uparrow(\downarrow)}=+1(-1)$.\\ In absence of spin flip scattering ($p=0$), Eq.~\ref{Eq:eb} reduces to,
\begin{equation}
E(\varphi)=\pm\Delta\sqrt{\frac{4\cos^2(\varphi/2)+J^2m'^2}{4+J^2m'^2}}.
\end{equation}
\subsection{Superconductor-Spin flipper-Superconductor junction as a thermodynamic system}
Before analyzing the thermodynamic behavior of setup (Fig.~1), we introduce the different thermodynamic quantities. From Andreev bound state energies we can determine Free energy of JJ\cite{golu} (inset of Fig.~1) as
\begin{eqnarray}
F(\varphi,T)&&=-\frac{1}{\beta}\ln\Big[\prod_{i}(1+e^{-\frac{E_{j}(\varphi)}{k_{B}T}})\Big],\nonumber\\
            &&= -\frac{2}{\beta}\sum_{\sigma}\ln\Big[2\cosh\Big(\frac{E_{\sigma}(\varphi)}{2k_{B}T}\Big)\Big],
\label{ff}
\end{eqnarray}
where $k_{B}$ is Boltzmann constant. From Free energy we can calculate total Josephson current\cite{deG} as
\begin{equation}
I(\varphi,T)=\frac{2e}{\hbar}\frac{\partial F(\varphi,T)}{\partial \varphi},
\label{ff1}
\end{equation}
where $e$ is charge of electron. Entropy of our device can then be calculated from Free energy as
\begin{equation}
S(\varphi,T)=-\frac{\partial F(\varphi,T) }{\partial T}.
\end{equation}
From entropy $S$ one determines heat capacity of JJ as
\begin{equation}
C(\varphi,T)=T\frac{\partial S(\varphi,T)}{\partial T}.
\end{equation}
\section{Josephson-Stirling Cycle}
Josephson-Stirling cycle\cite{carr,bsc} represented in Fig.~2, involves two isothermal and two isophasic processes. States $\textbf1$,$\textbf{2}$ involve right reservoir, while states $\textbf3$,$\textbf{4}$ involve left reservoir. Below we summarize these different processes-
\begin{itemize}
 \item \textbf{Isothermal process ($\textbf{1}\rightarrow\textbf{2}$)}: Thermal valve $v_{R}$ is open while $v_{L}$ is closed, thus system is in thermal contact with right reservoir at temperature $T_{R}$. The device or system  goes from state $\textbf{1}\equiv(\varphi=0,T_{R})$ to state $\textbf{2}\equiv(\varphi=\varphi_{f},T_{R})$.
 \item \textbf{Isophasic process ($\textbf{2}\rightarrow\textbf{3}$)}: Thermal valve $v_{L}$ is open while $v_{R}$ is closed and system is driven from state $\textbf{2}\equiv(\varphi=\varphi_{f},T_{R})$ to state $\textbf{3}\equiv(\varphi=\varphi_{f},T_{L})$.
 \item \textbf{Isothermal process ($\textbf{3}\rightarrow\textbf{4}$)}: In this stage valve $v_{R}$ is closed and $v_{L}$ is open. System is transferred from state $\textbf{3}\equiv(\varphi=\varphi_{f},T_{L})$ to state $\textbf{4}\equiv(\varphi=0,T_{L})$.
 \item \textbf{Isophasic process ($\textbf{4}\rightarrow\textbf{1}$)}: Final stage of cycle involves closing valve $v_{L}$ and opening valve $v_{R}$, with system being driven from state $\textbf{4}\equiv(\varphi=0,T_{L})$ to state $\textbf{1}\equiv(\varphi=0,T_{R})$.
\end{itemize}
\begin{figure}[ht]
\centering{\includegraphics[width=.99\linewidth]{Thermodynamic-cycle.pdf}}
\caption{\small \sl Josephson-Stirling cycle in TS (T is temperature and S is entropy) plane. Enclosed area in $T-S$ plane corresponds to  total heat exchanged $Q$, which is equal to total work done $W$ during the cycle.}
\end{figure}

When the device or system is driven from state $i$ to state $f$ ($i\rightarrow f$) during a quasi-static process, work done and heat released during the process are given as $W_{if}=-\frac{\hbar}{2e}\int_{\varphi_{i}}^{\varphi_{f}} I(\varphi,T)d\varphi$ and $Q_{if}=\int_{S(\varphi_{i},T)}^{S(\varphi_{f},T)} T dS$ respectively, where $I(\varphi,T)$ is Josephson current and $S$ is entropy of setup depicted in Fig.~1. In expressions for work done and heat released below, we have taken sign convention such that $W_{if}$ is positive when work is released to the universe while $Q_{if}$ is positive when heat is absorbed from the universe.
Work done and heat released for an isothermal process, where phase difference $\varphi$ changes from $\varphi_{i}\rightarrow\varphi_{f}$ at constant temperature $T$ is given as $W_{if}=-[F(\varphi_{f},T)-F(\varphi_{i},T)]$ and $Q_{if}=T[S(\varphi_{f},T)-S(\varphi_{i},T)]$ respectively. For an isophasic process, temperature changes from $T_{i}\rightarrow T_{f}$ at constant phase difference $\varphi$, $W_{if}=0$ and $Q_{if}=\int_{T_{i}}^{T_{f}}C(\varphi,T) dT$.
\subsection{Work done and heat exchanged in Josephson-Stirling Cycle}
We can now explicitly calculate total work done and heat exchanged during each stage of the Josephson-Stirling cycle, shown in Fig.~2. We also distinguish between four distinct modes of operation of the Josephson-Stirling cycle. In quantum heat engine mode: $W>0$ and $Q_{R}-Q_{L}>0$, in refrigerator mode: $W<0$ and $Q_{R}-Q_{L}>0$, while in Joule pump mode: $W<0$ and $Q_{R}+Q_{L}<0$, further in cold pump mode: $W<0$ and $Q_{R}-Q_{L}<0$. The total work done per cycle is $W=W_{12}+W_{34}$. Since $2\rightarrow3$ and $4\rightarrow1$ are isophasic processes, therefore $W_{23}=W_{41}=0$. Thus, work done $W$, can be calculated as,
\begin{widetext}
\footnotesize{
\begin{eqnarray}
W=&&\Delta\Bigg\{\Bigg(Y_{R}\ln\Bigg[4\cosh\Bigg(\frac{X_{R} E_{\uparrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]+Y_{R}\ln\Bigg[\cosh\Bigg(\frac{X_{R} E_{\downarrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]-Y_{L}\ln\Bigg[4\cosh\Bigg(\frac{X_{L} E_{\uparrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]-Y_{L}\ln\Bigg[\cosh\Bigg(\frac{X_{L} E_{\downarrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]\Bigg)\nonumber\\
&&-\Bigg(Y_{R}\ln\Big[\cosh\Big(X_{R}\Big)\Big]+Y_{R}\ln\Big[4\cosh\Big(X_{R} C^{\prime}\Big)\Big]-Y_{L}\ln\Big[\cosh\Big(X_{L}\Big)\Big]-Y_{L}\ln\Big[4\cosh\Big(X_{L}C^{\prime}\Big)\Big]\Bigg)\Bigg\},\label{work}\nonumber\\\\
\mbox{ where, }&&X_{L}=\frac{\Delta}{2k_{B}T_{L}},\,\, X_{R}=\frac{\Delta}{2k_{B}T_{R}},\,\, Y_{L}=\frac{1}{X_{L}},\,\, Y_{R}=\frac{1}{X_{R}},\,\, C_{1}=\sqrt{16+J^4 (p^2+m'+m'^2)^2+J^2 (4-8p^2+8m'(1+m'))},\nonumber\\
C_{2}=&&\sqrt{16+J^4(p^2+m'+m'^2)^2+J^2(4+8p^2+8m'(1+m'))},\, \mbox{ and } C^{\prime}=\frac{C_{1}}{C_{2}}.\nonumber
\end{eqnarray}}
\normalsize
In absence of spin flip scattering ($p=0$), work done, \textquoteleft$W$\textquoteright{} in Eq.~\eqref{work}, simplifies to,
\begin{equation}
W=\Delta\Bigg\{\Bigg(Y_{L}\ln\Big[\cosh\Big(X_{L}\Big)\Big]-Y_{R}\ln\Big[\cosh\Big(X_{R}\Big)\Big]\Bigg)-\Bigg(Y_{L}\ln\Bigg[\cosh\Bigg(\frac{X_{L}E(\varphi_{f})}{\Delta}\Bigg)\Bigg]-Y_{R}\ln\Bigg[\cosh\Bigg(\frac{X_{R}E(\varphi_{f})}{\Delta}\Bigg)\Bigg]\Bigg)\Bigg\},\label{worknoflip}
\end{equation}
where, $E(\varphi)=\Delta\sqrt{\frac{4\cos^2(\varphi/2)+J^2m'^2}{4+J^2m'^2}}$ is Andreev bound state energy for no flip case.

The heat exchanged with right reservoir at temperature $T_{R}$, $Q_{R}=Q_{12}+Q_{41}$ is
\small{
\begin{eqnarray}
Q_{R}=&&\Delta\Bigg(\ln\Bigg[\cosh\Bigg(\frac{X_{R} E_{\uparrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]+\ln\Bigg[4\cosh\Bigg(\frac{X_{R} E_{\downarrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]\Bigg)Y_{R}-E_{\uparrow}(\varphi_{f})\tanh\Bigg(\frac{X_{R}E_{\uparrow}(\varphi_{f})}{\Delta}\Bigg)-E_{\downarrow}(\varphi_{f})\tanh\Bigg(\frac{X_{R}E_{\downarrow}(\varphi_{f})}{\Delta}\Bigg)\nonumber\\
     &&-\Delta\Bigg(\ln\Big[\cosh\Big(X_{L} C^{\prime}\Big)\Big]+\ln\Big[4\cosh\Big(X_{L}\Big)\Big]\Bigg)Y_{L}
     +\Delta C^{\prime}\tanh\Big(X_{L} C^{\prime}\Big)+\Delta\tanh\Big(X_{L}\Big)\nonumber\\
     &&+\Delta\Bigg\{\Bigg(\ln\Big[\cosh\Big(X_{L}\Big)\Big]+\ln\Big[4\cosh\Big(X_{L}C^{\prime}\Big)\Big]\Bigg)Y_{L}-\Bigg(\ln\Big[\cosh\Big(X_{R}\Big)\Big]+\ln\Big[4\cosh\Big(X_{R} C^{\prime}\Big)\Big]\Bigg)Y_{R}\Bigg\},\nonumber\\
\label{hrr}
\end{eqnarray}}
\normalsize and for no flip process ($p=0$), Eq.~\eqref{hrr} reduces to,
\begin{equation}
Q_{R}=\Delta\Bigg(-\ln\Big[\cosh\Big(X_{R}\Big)\Big]+\ln\Bigg[\cosh\Bigg(\frac{X_{R}E(\varphi_{f})}{\Delta}\Bigg)\Bigg]\Bigg)Y_{R}+\Delta\tanh\Big(X_{L}\Big)-E(\varphi)\tanh\Bigg(\frac{X_{R}E(\varphi_{f})}{\Delta}\Bigg).
\end{equation}
Finally, heat exchanged with left reservoir at temperature $T_{L}$, $Q_{L}=Q_{23}+Q_{34}$ is
\footnotesize{
\begin{eqnarray}
Q_{L}=&&-\Delta\Bigg(\ln\Bigg[\cosh\Bigg(\frac{X_{R} E_{\uparrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]+\ln\Bigg[4\cosh\Bigg(\frac{X_{R} E_{\downarrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]\Bigg)Y_{R}+E_{\uparrow}(\varphi_{f})\tanh\Bigg(\frac{X_{R}E_{\uparrow}(\varphi_{f})}{\Delta}\Bigg)+E_{\downarrow}(\varphi_{f})\tanh\Bigg(\frac{X_{R}E_{\downarrow}(\varphi_{f})}{\Delta}\Bigg)\nonumber\\
     &&+\Delta\Bigg(\ln\Big[\cosh\Big(X_{L} C^{\prime}\Big)\Big]+\ln\Big[4\cosh\Big(X_{L}\Big)\Big]\Bigg)Y_{L}-\Delta C^{\prime}\tanh\Big(X_{L} C^{\prime}\Big)-\Delta\tanh\Big(X_{L}\Big)\nonumber\\
     &&-\Delta\Bigg\{\Bigg(\ln\Bigg[4\cosh\Bigg(\frac{X_{L} E_{\uparrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]+\ln\Bigg[\cosh\Bigg(\frac{X_{L} E_{\downarrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]\Bigg)Y_{L}-\Bigg(\ln\Bigg[4\cosh\Bigg(\frac{X_{R} E_{\uparrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]+\ln\Bigg[\cosh\Bigg(\frac{X_{R} E_{\downarrow}(\varphi_{f})}{\Delta}\Bigg)\Bigg]\Bigg)Y_{R}\Bigg\},\nonumber\\
\label{hlr}
\end{eqnarray}}
\normalsize and for no flip process ($p=0$), Eq.~\eqref{hlr}, reduces to,
\begin{equation}
Q_{L}=\Delta\Bigg(\ln\Big[\cosh\Big(X_{L}\Big)\Big]-\ln\Bigg[\cosh\Bigg(\frac{X_{L}E(\varphi_{f})}{\Delta}\Bigg)\Bigg]\Bigg)Y_{L}-\Delta\tanh\Big(X_{L}\Big)+E(\varphi)\tanh\Bigg(\frac{X_{R}E(\varphi_{f})}{\Delta}\Bigg).
\end{equation}
\end{widetext}
Upper limit for $T_{L}$ and $T_{R}$ is constrained by superconducting transition temperature $T_{c}$. For \textquoteleft Pb\textquoteright{} superconductor, $T_{c}$ is $7$K, therefore upper limit for $T_{R}$ and $T_{L}$ is taken as $6.96$K ($<T_{c}$). In Figs.~3-8, we consider lead ($Pb$) superconductor, thus upper limit for $T_{R}$, $T_{L}$ is fixed at $6.96$K. From conservation of energy $W=Q$, with $Q=Q_{R}+Q_{L}$ being total heat exchanged during the cycle. When Josephson-Stirling cycle acts as a heat engine, work done $W$ and efficiency of cycle is $\eta=W/Q_{R}$. In refrigerator mode, the work done $W$ and coefficient of performance\cite{carr} of Josephson-Stirling cycle is given as $\mbox{COP}=Q_{R}/|W|$. Finally, in Joule pump mode the $COP$ is\cite{biz} $(|Q_{L}|+|Q_{R}|)/|W|$ with $Q_R$ being heat released to right reservoir and $Q_L$, being the heat released to left reservoir as Joule pump coverts the work to heat and transfers it to both reservoirs. A cold pump, on the other hand, has $COP=|Q_R|/|W|$\cite{biz}, as in cold pump $Q_R$ is negative implying heat energy is released to right reservoir which is cold reservoir because $T_L  > T_R$. A cold pump is used to heat the colder reservoir by transferring heat from hotter reservoir in turn cooling the hotter reservoir. Please also note that the expressions for COP in case of Joule pump  differs slightly from that in Ref.~\cite{biz} as we have taken the sign convention of heat flow or work done as in Ref.~\cite{bsc}.
\section{Results and Analysis}
\subsection{Josephson-Stirling cycle as quantum heat engine}
\begin{figure}[ht]
\centering{\includegraphics[width=.99\linewidth]{fig3-jos.pdf}}
\caption{\small \sl Total work done $W$ and efficiency $\eta$ of a Josephson-Stirling engine as a function of $\varphi_{f}$ (for (a,c)) and as a function of the right reservoir temperature $T_{R}$ (for (b,d)). Parameters are: Flip: $S=-m'=1/2$, $p=1$; No Flip: $S=m'=1/2$, $p=0$. Rest of parameters: $J=0.1$, $\varphi_{i}=0$, $\varphi_{f}=\pi$ (for (b,d)), $k_{B}T_{R}=0.6\Delta$ (for (a,c)), $k_{B} T_{L}=0.01\Delta$, $k_{B}T=0.1\Delta$. The black dotted line represents Carnot limit, $\eta_{C}=1-\frac{T_{L}}{T_{R}}$.}
\end{figure}
In Fig.~3 we plot work done $W$ (Eq.~\eqref{work}) and efficiency $\eta$ of Josephson quantum heat engine and compare no flip with spin flip processes. In Figs.~3(a) and 3(c), $W$ and $\eta$ are plotted as function of maximal phase change $\varphi_{f}$ during the cycle. The reason we plot work done, efficiency and COP as function of phase is because phase lends itself to external control via magnetic flux enclosed in Josephson junction loop, see Fig.~1. We consider here right reservoir to be hot and left reservoir to be cold, i.e., $T_{R}>T_{L}$. In the inset of Fig.~3(a), we show Free energy as function of $\varphi_{f}$. In Fig.~3(a) we see that $W$ is maximum at $\varphi_{f}=\pi$, irrespective of spin-flip scattering, however magnitude of $W_{max}$ for no flip process is much smaller than spin-flip process. In presence of spin-flip scattering magnitude of $W$ at $\varphi_{f}=\pi$, i.e., $W=1.15\Delta$ is maximum work done by our system. From inset of Fig.~3(a), we notice that Free energy is minimum at $\varphi_{f}=0$, i.e., when system shows $0$-junction behavior. From Fig.~3(c), we notice that at $\varphi_{f}=\pi$, efficiency $\eta$ for no flip case is a little bit larger than spin-flip case. Thus, spin flip process is better for work done, but no flip process is better for efficiency. In Figs.~3(b) and 3(d) we show dependence of work done and efficiency on $T_{R}$ for fixed $k_{B}T_{L}=0.01\Delta$ and compare the cycle's performance for no flip and spin flip scattering. In addition we plot Carnot efficiency in Fig.~3(d). Three important take home messages from Fig.~3 are (i) work done is maximum at $\varphi_{f}=\pi$ (see Fig.~3(a)), (ii) spin-flip scattering enhances work done (see Fig.~3(a)) and (iii) at $k_{B}T_{R}=0.6\Delta$, work done is maximum (see Fig.~3(b)).

From Fig.~3(d), we notice that for small values of $T_{R}$, efficiency of a Josephson-Stirling engine is equal to Carnot efficiency $\eta_{C}$ regardless of spin-flip scattering, while at large values of $T_{R}$ in absence of spin-flip scattering, engine is more efficient and efficiency is close to Carnot limit. The reason we take low values of exchange coupling $J$ ($=0.1$) is because for these values of $J$ we get maximum work output \textquoteleft$W$\textquoteright{} and efficiency \textquoteleft$\eta$\textquoteright.
\begin{figure}[ht]
\centering{\includegraphics[width=.99\linewidth]{fig4-jos.pdf}}
\caption{\small \sl  Total work done $W$ and efficiency $\eta$ of a Josephson-Stirling engine as function of $\varphi_{f}$ (for (a,c)) and as function of $T_{R}$ (for (b,d)) in presence of spin-flip scattering. Other parameters are: $J=0.1$,  $\varphi_{i}=0$, $\varphi_{f}=\pi$ (for $p=1$), $\varphi_{f}=0.92\pi$ (for $p=3$), $\varphi_{f}=0.85\pi$ (for $p=5$), $k_{B}T_{L}=0.01\Delta$, $k_{B}T_{R}=0.6\Delta$ (for (a,c)), $k_{B}T=0.1\Delta$.}
\end{figure}
Next, in Fig.~4 we plot total work done $W$ and efficiency $\eta$ of a Josephson-Stirling engine for different levels of spin-flip scattering. In Figs.~4(a) and 4(c) we plot $W$ and $\eta$ as function of $\varphi_{f}$. We see that at $\varphi_{f}=\pi$ and for low values of spin-flip scattering, both work done and efficiency attain their maximum. In inset of Fig.~4(a), we notice that minimum of Free energy occurs at $\varphi_{f}=0$, i.e., system shows $0$-junction behavior when it operates as a Josephson quantum heat engine. In Figs.~4(b) and 4(d) we plot $W$ and $\eta$ as function of $T_{R}$ for different levels of spin flip scattering at fixed $k_{B}T_{L}=0.01\Delta$. We choose $\varphi_{f}$ in such a way that $W$ becomes maximum. We see that for low values of spin-flip scattering, both work done and efficiency are larger.
\subsubsection{Condition for optimality}
Further we find in Fig.~4 that, at $k_{B}T_{R}=0.6\Delta$, both $W$ and $\eta$ are maximum. Thus, for spin-flip probability, $p=1$, $J=0.1$, $\varphi_{f}=\pi$, $k_{B}T_{L}=0.01\Delta$ and $k_{B}T_{R}=0.6\Delta$ both work done and efficiency of our Josephson quantum heat engine are larger and therefore this regime is an optimal operating range of our Josephson quantum heat engine. This condition of optimality for our case is contrasted with that in Ref.~\cite{bsc}. In Ref.~\cite{bsc}, parity-conserving heat engine operates in an optimal operating range at $\varphi_{f}=2\pi$, while non-parity-conserving heat engine works in an optimal operating range at $\varphi_{f}=\pi$.

\subsection{Josephson-Stirling cycle as quantum refrigerator}
Till now we have discussed work done and efficiency when Josephson-Stirling cycle acts as a quantum heat engine. In Fig.~5 we show action of Stirling cycle as a Josephson quantum refrigerator and calculate work done and coefficient of performance (COP). In Figs.~5(a) and 5(b), total work absorbed by refrigerator and COP are plotted as function of $\varphi_{f}$, for $T_{R}<T_{L}$ ($k_{B}T_{L}=0.4\Delta$ and $k_{B}T_{R}=0.39\Delta$). We see that for low values of spin flip scattering both $|W|$ and COP are larger and the product of COP and $|W|$ is maximum at $\varphi_{f}=\pi$. Thus,  Josephson quantum refrigerator works in an optimal operating range for parameters of Fig.~5. For optimal parameters, $|W|=0.017\Delta$ and $\mbox{COP}=38.48$. Although, our system exhibits maximum $|W|$ of $0.65\Delta$ at $k_{B}(T_{L}-T_{R})=0.3\Delta$ and a maximum COP of $58.27$ at $k_{B}(T_{L}-T_{R})=0.01\Delta$ as shown in Figs.~5(c) and 5(d) respectively.
\begin{figure}[ht]
\centering{\includegraphics[width=.99\linewidth]{fig5-jos.pdf}}
\caption{\small \sl (a,c) Total work $W$ absorbed, and (b,d) COP of quantum refrigerator as function of $\varphi_{f}$ in presence of spin-flip scattering. Other parameters are: $J=0.1$, $k_{B}T_{L}=0.4\Delta$ (for (a,b)), $k_{B}T_{R}=0.39\Delta$ (for (a,b)), $S=-m'=1/2$ (for (c,d)), $\varphi_{i}=0$.}
\end{figure}
\subsection{Phase diagram of Josephson-Stirling cycle}
In Figs.~3,4,5 we depict the performance of this device as a quantum heat engine or quantum refrigerator depending on relative temperatures of reservoirs. In this subsection we will discuss phase diagram of a Josephson-Stirling cycle as function of phase difference $\varphi_{f}$ during the cycle. In Fig.~6, we plot $W$, $Q_{R}$, $Q_{L}$, $\eta$ and COP as functions of $\varphi_{f}$ and $T_{R}$ for low values of spin flip scattering ($S=1/2$, $m'=-1/2$, $p=1$) and high values of exchange coupling ($J=2$) of spin flipper. We fix $k_{B}T_{L}=0.01\Delta$ and consider $T_{R}>T_{L}$. From Fig.~6(a) we notice that with change of  $\varphi_{f}$ there is a sign change in $W$. $W$ changes sign from $W>0$ to $W<0$ between $\varphi_{f}=0$ and $\varphi_{f}=\pi$ and again changes sign from $W<0$ to $W>0$ between $\varphi_{f}=\pi$ and $\varphi_{f}=2\pi$. The regime where $W$ is positive ($W>0$) in Fig.~6(a), $Q_{R}$ is positive (Fig.~6(b)) but $Q_{L}$ is negative (Fig.~6(c)), in these regimes Josephson-Stirling cycle acts as an engine. Similarly in regime, wherein both $W$ and $Q_{L}$ are negative, cycle acts as a refrigerator if $Q_{R}$ is positive while cycle acts as a Joule pump if $Q_{R}$ is negative. A Joule pump transforms work into heat released to the reservoirs. Thus, there is a transition from engine mode to refrigerator or Joule pump mode with change of $\varphi_{f}$. Since, $\varphi_{f}$ is controlled by magnetic flux $\Phi$, therefore by changing the enclosed flux $\Phi$ in the Josephson junction loop we can tune our system from engine mode to refrigerator/Joule pump mode which is an attractive feature of our device.
\begin{figure}[ht]
\centering{\includegraphics[width=.99\linewidth]{fig6-jos.pdf}}
\caption{\small \sl (a) Total work $W$ (b) $Q_{R}$ and (c) $Q_{L}$, (d) $\eta$ or COP as as function of $\varphi_{f}$ and $T_{R}$ in presence of spin flip scattering. Other parameters are: $S=1/2$, $m'=-1/2$, $J=2$, $k_{B}T_{L}=0.01\Delta$, $\varphi_{i}=0$.}
\end{figure}

\subsection{Why is spin-flip scattering necessary?}
In absence of spin-flip scattering ($p=0$), we find work done (see Eq.~\eqref{worknoflip}) $W>0$ since,
\begin{eqnarray}
(Y_{L}\ln[\cosh(X_{L})]-Y_{R}\ln[\cosh(X_{R})]) &>& (Y_{L}\ln[\cosh(X_{L}E(\varphi_{f})/\Delta)]\nonumber\\
-Y_{R}\ln[\cosh(X_{R}E(\varphi_{f})/\Delta)]),
\end{eqnarray}
for $T_{R}>T_{L}$ and regardless of $\varphi_{f}$. Thus, $W$ is always positive ($W>0$) in absence of spin-flip scattering irrespective of $\varphi_{f}$ and the Josephson-Stirling cycle operates solely as a quantum heat engine. However, in presence of spin-flip scattering ($p\neq0$) and with $T_{R}>T_{L}$ , we find $W$ (see Eq.~\eqref{work}) is not always positive and depending on other parameter values, $W$ can be negative ($W<0$) in certain range of $\varphi_{f}$. In the parameter regime wherein $W$ is negative and $Q_{L}$ is negative with $Q_{R}$ being positive, the Josephson-Stirling cycle operates as a refrigerator. Similarly, in parameter regime wherein $W$, $Q_{R}$ and $Q_{L}$ are all negative, Josephson-Stirling cycle acts as a Joule pump. Thus, in presence of spin-flip scattering when $W$ is negative, the Josephson-Stirling cycle can act as a refrigerator or Joule pump even when $T_{R}>T_{L}$. This is unique to this proposal since in other proposals\cite{bsc,carr}, $W$ is always positive and Josephson-Stirling cycle only acts as a heat engine when $T_{R}>T_{L}$.
\begin{figure}[ht]
\centering{\includegraphics[width=.99\linewidth]{fig7-jos.pdf}}
\caption{\small \sl (a) Total work done $W$ (b) $Q_{R}$ and (c) $Q_{L}$  (d) efficiency $\eta$ as function of maximal phase change $\varphi_{f}$ during Stirling cycle and right reservoir temperature $T_{R}$ in absence of spin-flip scattering. Other parameters are: $S=m'=1/2$, $J=2$, $k_{B}T_{L}=0.01\Delta$, $\varphi_{i}=0$.}
\label{spinflip-fig}
\end{figure}
In Fig.~\ref{spinflip-fig}, $W$, $Q_{R}$, $Q_{L}$ and $\eta$ are plotted as function of $\varphi_{f}$ and $T_{R}$ in absence of spin-flip scattering ($S=m'=1/2$, $p=0$). We see that $W$ does not change sign with $\varphi_{f}$ and thus there is no transition from engine mode to refrigerator or Joule pump mode. Therefore, we can conclude that it is spin flip scattering, which is responsible for tuning our system from engine mode to refrigerator/Joule pump mode via magnetic flux $\Phi$.
\subsection{Josephson-Stirling cycle as Joule pump and cold pump}
Finally, our device doesn't only exhibit the engine, refrigerator and Joule pump phases, it also exhibits a cold pump phase again in presence of spin-flip scattering only. The thin sliver in Fig.~\ref{cold-pump}(a) at bottom right corner shows the cold pump phase. In Fig.~\ref{cold-pump}, we plot $W$, $Q_{R}$, $Q_{L}$, efficiency of  engine and COP of refrigerator as functions of $T_{R}$ and $T_{L}$ for low values of spin flip scattering ($S=1/2$, $m'=-1/2$, $p=1$) and low values of exchange coupling ($J=0.1$) of spin flipper. We choose $\varphi_{f}$ such that $|W|$ will be maximized. We see that for $T_{R}>T_{L}$, the Josephson-Stirling cycle acts as an engine. At $k_{B}T_{L}=0.01\Delta$ and $k_{B}T_{R}=0.6\Delta$, the work done and efficiency in the engine mode attain their maximum values which are $W_{max}=1.15\Delta$ and $\eta_{max}=0.923$. For $T_{R}<T_{L}$, we notice that the machine acts as a refrigerator or as Joule pump or cold pump depending on the relative values of $Q_{R}$ and $Q_{L}$. The work done in the refrigerator mode is maximum at $k_{B}T_{L}=0.44\Delta$ and $k_{B}T_{R}=0.14\Delta$ and the maximum value is $|W|_{max}=0.65\Delta$, while the COP of the refrigerator is maximum at $k_{B}T_{L}=0.6\Delta$ and $k_{B}T_{R}=0.59\Delta$ and the maximum value is $COP_{max}=58.27$. Further, in the Joule pump mode the work done is maximum at $k_{B}T_{L}=0.58\Delta$ and $k_{B}T_{R}=0.01\Delta$ and the maximum value is $|W|_{max}=1.13\Delta$, while the COP of the Joule pump is maximum at $k_{B}T_{L}=0.16\Delta$ and $k_{B}T_{R}=0.01\Delta$ and the maximum value is $COP_{max}=1.2$. If $Q_{R}<0$, $W<0$ and $Q_{L}>0$, the cycle acts as a cold pump which transfers heat from the hot reservoir to the cold reservoir. The work done in the cold pump mode is maximum at $k_{B}T_{L}=0.6\Delta$ and $k_{B}T_{R}=0.01\Delta$ and the maximum value is $|W|_{max}=1.15\Delta$, while the COP of the cold pump is maximum at $k_{B}T_{L}=0.6\Delta$ and $k_{B}T_{R}=0.01\Delta$ and the maximum value is $COP_{max}=0.55$. Thus, by tuning the temperature of the reservoirs our system can be tuned from engine mode to refrigerator or Joule pump or cold pump mode.
\begin{figure}[ht]
\centering{\includegraphics[width=.99\linewidth]{fig8-jos.pdf}}
\caption{\small \sl (a) Total work done $W$, (b) heat exchanged $Q_{R}$ and (c) $Q_{L}$ and  (d) efficiency $\eta$ or COP as function of temperatures $T_{R}$ and $T_{L}$ in presence of spin-flip scattering. Parameters are: $S=1/2$, $m'=-1/2$, $J=0.1$, $\varphi_{i}=0$, $\varphi_{f}=\pi$.}
\label{cold-pump}
\end{figure}
\section{Experimental realization and Conclusions}
Our proposed Josephson quantum spin thermodynamic device acts as a heat engine or refrigerator or Joule pump or as a cold pump via spin-flip scattering and is experimentally realizable. Doping a magnetic adatom or spin-flipper in an one dimensional superconducting loop shouldn't be difficult, especially with a $s$-wave superconductor like Lead or Aluminum it should be perfectly possible. Further, work done per cycle during each isothermal process can be experimentally determined via the current-phase relation. For a Josephson Stirling cycle work done during each isothermal process is given as $W_{if}=-\frac{\hbar}{2e}\int_{\varphi_{i}}^{\varphi_{f}}I(\varphi,T)d\varphi$, where $\varphi_{i}$ and $\varphi_{f}$ are the initial and final phases respectively. Thus, to measure the work done experimentally one needs to control phase difference $\varphi$ across the junction, regulated by flux $\Phi$, as well as to know about the Josephson current $I$ flowing through the loop. This can be done using a scanning superconducting quantum interference device (SQUID) microscope to perform the measurements of current-phase relation of the junction\cite{bsc,iso,sha}. By applying a current through the field coil of SQUID sensor, the magnetic flux through Josephson junction loop can be tuned. This magnetic flux, controls the phase difference across JJ and induces a supercurrent in Josephson junction loop. The supercurrent leads to a signal which is measured by pickup loop of SQUID sensor\cite{iso,iko}.

We have compared our proposal with other Josephson quantum heat engines and refrigerators in Table I. While, Refs.~\onlinecite{bsc}, \onlinecite{carr} are very important works which laid down the principle of Josephson quantum heat engine and refrigerator, efficiency and coefficient of performance (COP) of these can none-the-less be still enhanced, and the tunability of the device improved as we show in this work.   In both Refs.~\onlinecite{bsc} and \onlinecite{carr}, the Josephson quantum heat engine and refrigerator can be tuned from engine mode to refrigerator mode via changing the temperatures of left and right reservoirs only, while in this proposal it is seen that heat engine to refrigeration mode transition or for that matter to Joule or cold pump can be effected by either the magnetic flux enclosed or by tuning the temperature of reservoirs. This makes spin flipper doped Josephson junction loop much more versatile as regards possible applications. As a quantum heat engine, the proposed device is much more efficient than previously proposed Josephson quantum heat engines\cite{bsc,carr}. Further, when operating as a quantum refrigerator, COP of the proposed device is higher than that seen in Refs.~\cite{bsc,carr}. In Table I, we see that although work done by Josephson heat engine proposed in Ref.~\onlinecite{carr} is larger than ours, but efficiency of the Josephson quantum heat engine discussed in this work is much larger than those of Josephson heat engines of Refs.~\onlinecite{bsc,carr}. Further, as a Josephson quantum refrigerator, COP of the proposed device is huge compared to other Josephson quantum refrigerator proposals.
\begin{table}[ht]
\caption{Comparative analysis of Josephson junction based quantum heat engines and refrigerators}
\begin{tabular}{|p{1.7cm}|p{1.7cm}|p{1.67cm}|p{1.38cm}|p{1.68cm}|}
\hline
& \multicolumn{2}{|c|}{Heat engine mode} & \multicolumn{2}{|c|}{Refrigerator mode} \\
\hline
& $W_{max}$ & $\eta_{max}$ & $|W_{max}|$ & $\mbox{COP}_{max}$ \\
\hline
The JJ device (Fig.~1) & $1.15\Delta$ ($\Delta=1meV$) (Fig.~3(a)) & $0.958$ (Fig.~3(c)) & $0.65\Delta$ (Fig.~5(c))  & $58.27$ (Fig.~5(d)) \\
\hline
Topological Josephson heat engine (Ref.~\onlinecite{bsc}) & $2\Delta$ ($\Delta=150\mu eV$) (Fig.~4(c) of Ref.~\onlinecite{bsc}) & $0.8$ (Figs.~4(b,d) of Ref.~\onlinecite{bsc}) & \hspace{0.5cm}$-$ & $20$  (Figs.~4(b,d) of Ref.~\onlinecite{bsc}) \\
\hline
Josephson heat engine (Ref.~\onlinecite{carr}) & $0.38\Delta$ ($\Delta=180eV$) (Fig.~14(a) of Ref.~\onlinecite{carr})  & $0.5$ (Fig.~15(a) of Ref.~\onlinecite{carr}) & \hspace{0.5cm}$-$ & $10$ (Fig.~15(a) of Ref.~\onlinecite{carr}) \\
\hline
\end{tabular}
\end{table}
To conclude, in this work we have shown that a 1D Josephson junction loop doped with a spin flipper, and attached to two thermal reservoirs at inequivalent temperatures via thermal valves can act as a Josephson quantum heat engine or  quantum refrigerator or as a Joule pump or even as a cold pump. The proposed device can be tuned from engine mode to other operating modes like refrigerator or Joule pump/Cold pump not only via changing the temperatures of the reservoirs, but also via the enclosed flux in Josephson junction loop which is the most lucrative aspect of work, since this fact alone implies a much greater possibility of our proposal being experimentally realized.
\acknowledgments
This work was supported by the grants: 1. Josephson junctions with strained Dirac materials and their application in quantum information processing, SERB Grant No. CRG/20l9/006258, and 2. Nash equilibrium versus Pareto optimality in N-Player games, SERB MATRICS Grant No. MTR/2018/000070.