\documentclass[fleqn,10pt]{wlscirep}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
%\usepackage{hyperref}% add hypertext capabilities
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math
%\linenumbers\relax % Commence numbering lines
\usepackage[utf8]{inputenc}
\usepackage{graphicx}
\usepackage{color}
%\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{subfigure}
\usepackage{epsfig}
\usepackage{float}
\usepackage{lipsum}
\usepackage{cases}
\usepackage{appendix}
%\usepackage{TikZ}
\usepackage{cancel}
\usepackage{soul}
\usepackage{cite}
%-------------------------------------------------------------------------------
\newcommand{\sig}{\vec{\sigma}}
\newcommand{\kapp}{\vec{\kappa}}
\newcommand{\pdf}{\Psi(\{\textbf{r}_i\},t)}
\newcommand{\om}{\Omega}
\newcommand{\omdag}{{\Omega}^{\dagger}}
\newcommand{\omeq}{\Omega_{\rm eq}}
\newcommand{\omeqdag}{{\Omega_{\rm eq}}^{\dagger}}
\newcommand{\delom}{\delta\Omega}
\newcommand{\delomdag}{{\delta\Omega}^{\dagger}}
\newcommand{\gammadot}{\dot{\gamma}}


\newcommand{\cc}{\boldsymbol{c}}
\newcommand{\bb}{\boldsymbol{b}}
\newcommand{\jj}{\vec{j}}

\newcommand{\D}{\boldsymbol{D}}
\newcommand{\sR}{\vec{\mathcal{R}}}

%\newcommand{\comment}[1]{{\color{red}{\bf #1}}}
\newcommand{\G}{\boldsymbol{G}}
%\newcommand{\Gama}{\mathbb{D}}
\newcommand{\Gama}{{\mathpalette\makeGama\relax}}
\newcommand{\makeGama}[2]{%
  \raisebox{\depth}{\scalebox{1}[-1]{$\mathsurround=0pt#1\mathbb{L}$}}%
}
\newcommand{\nd}{\boldsymbol{S}}
\newcommand{\Q}{\boldsymbol{Q}}
\newcommand{\F}{\boldsymbol{F}}
\newcommand{\f}{\boldsymbol{f}}
\newcommand{\M}{\mathrm{\mathbb{M}}}
\newcommand{\sigg}{\vec{\sigma}}
\newcommand{\siggh}{\hat{\boldsymbol{\sigma}}}
\newcommand{\B}{\boldsymbol{B}}

\newcommand{\T}{\boldsymbol{\mathrm{T}}}
\newcommand{\A}{\boldsymbol{A}}
\newcommand{\1}{\mathbb{I}}
\newcommand{\LL}{\boldsymbol{\mathrm{L}}}
\newcommand{\Lambdaa}{\boldsymbol{\Lambda}}

\newcommand{\CC}{\boldsymbol{C}}

\newcommand{\diag}{\boldsymbol{diag}}
\newcommand{\kk}{\boldsymbol{k}}
\newcommand{\qq}{\boldsymbol{z}}
%\newcommand{\pp}{\boldsymbol{p}}
\newcommand{\nn}{\boldsymbol{n}}
%\newcommand{\pp}{\boldsymbol{p}}
%\newcommand{\nn}{\boldsymbol{n}}
%\newcommand{\J}{\boldsymbol{J}}
\newcommand{\J}{\boldsymbol{\mathrm{J}}}
%\newcommand{\jj}{\boldsymbol{j}}
\newcommand{\E}{\textbf{E}}
\newcommand{\XX}{\boldsymbol{\mathcal{X}}}
\newcommand{\AAA}{\boldsymbol{\mathcal{A}}}
\newcommand{\BB}{\boldsymbol{\mathcal{B}}}
\newcommand{\Cov}{\boldsymbol{\mathrm{S}}}
\newcommand{\Z}{\mathrm{Z}}
\newcommand{\dif}{\mathop{}\!\mathrm{d}}
\newcommand{\Det}{\mathop{}\!\mathrm{Det}}
\newcommand{\K}{\mathop{}\!\mathrm{K}}
\newcommand{\V}{\mathop{}\!\mathrm{V}}
\newcommand{\tr}{\mathop{}\!\mathrm{tr}}
%\newcommand{\diag}{\mathop{}\!\boldsymbol{\mathrm{diag}}}
\newcommand{\sign}{\mathop{}\!\boldsymbol{\mathrm{sign}}}
\newcommand{\uu}{\boldsymbol{u}}
\newcommand{\phii}{\boldsymbol{\Phi}}

\newcommand{\gam}{\dot{\gamma}}
\newcommand{\strain}{\boldsymbol{\varepsilon}}
\newcommand{\joe}{\color{red} }
%\newcommand{\p}{\boldsymbol{p}}
\newcommand{\p}{\boldsymbol{p}}
\newcommand{\pii}{\boldsymbol{\pi}}
\newcommand{\vv}{\boldsymbol{v}} 
\newcommand{\ww}{\boldsymbol{w}}
%\newcommand{\rr}{\boldsymbol{r}}
\newcommand{\rr}{\boldsymbol{r}}
\newcommand{\tet}{\boldsymbol{\theta}}

\newcommand{\0}{\boldsymbol{0}}
\newcommand{\ov}{\boldsymbol{s}}
\newcommand{\xx}{\boldsymbol{x}}
\newcommand{\aalpha}{\boldsymbol{\alpha}}

\newcommand{\zz}{\boldsymbol{\mathrm{z}}}
\newcommand{\yy}{\boldsymbol{y}}
\newcommand{\nab}{\boldsymbol{\nabla}}
\newcommand{\xxi}{\boldsymbol{\xi}}

\newcommand{\xchi}{\boldsymbol{\chi}}
\newcommand{\eeta}{\boldsymbol{\chi}}

%\newcommand{\xxi}{\boldsymbol{\xi}}

\newcommand{\R}{\boldsymbol{\mathrm{R}}}
\newcommand{\U}{\boldsymbol{\mathrm{U}}}
%\newcommand{\G}{\boldsymbol{\mathrm{G}}}
\newcommand{\Phik}{\Phi_{\boldsymbol{k}}}
\newcommand{\Gammak}{\Gamma_{\boldsymbol{k}}}
\newcommand{\mk}{m_{\boldsymbol{k}}}

\newcommand{\irow}[1]{% inline row vector
  \begin{matrix}(#1)\end{matrix}%
}
\newcommand{\rhoq}{\rho_{\boldsymbol{q}}^{\phantom{\ast}}}
\newcommand{\rhok}{\rho_{\boldsymbol{k}}^{\phantom{\ast}}}
\newcommand{\rhop}{\rho_{\boldsymbol{p}}^{\phantom{\ast}}}
\newcommand{\rhol}{\rho_{\boldsymbol{l}}^{\phantom{\ast}}}
\newcommand{\rhosq}{\rho^{\ast}_{\boldsymbol{q}}}
\newcommand{\rhosk}{\rho^{\ast}_{\boldsymbol{k}}}
\newcommand{\rhosp}{\rho^{\ast}_{\boldsymbol{p}}}
\newcommand{\rhosl}{\rho^{\ast}_{\boldsymbol{l}}}

\newcommand{\rhoqq}{\rho_{\boldsymbol{q}'}^{\phantom{\ast}}}
\newcommand{\rhokk}{\rho_{\boldsymbol{k}'}^{\phantom{\ast}}}
\newcommand{\rhopp}{\rho_{\boldsymbol{p}'}^{\phantom{\ast}}}
\newcommand{\rholl}{\rho_{\boldsymbol{l}'}^{\phantom{\ast}}}
\newcommand{\rhosqq}{\rho^{\ast}_{\boldsymbol{q}'}}
\newcommand{\rhoskk}{\rho^{\ast}_{\boldsymbol{k}'}}
\newcommand{\rhospp}{\rho^{\ast}_{\boldsymbol{p}'}}
\newcommand{\rhosll}{\rho^{\ast}_{\boldsymbol{l}'}}
\newcommand{\dbar}{\lower0.08ex\hbox{$\mathchar'26$}\mkern-12mu \dif}
\newcommand{\cmnt}[1]{{\color{red}{\bf #1}}} %I am not sure or I do not know
\newcommand{\amnt}[1]{{\color{blue}{#1}}} %amendment
\newcommand{\gl}[1]{{\color{green}{\bf #1}}} %needs your greenlight

\newcommand{\hcmnt}[1]{{\color{cyan}{\bf #1}}}


%-------------------------------------------------------------------------------

\title{Supplemental Materials: Tunable Brownian magneto heat pump}

\author[1]{Iman Abdoli}
\author[2]{Ren\'e Wittmann}
\author[3]{Joseph Michael Brader}
\author[1,4]{Jens-Uwe Sommer}
\author[2]{Hartmut L\"owen}
\author[1,4,*]{Abhinav Sharma}
\affil[1]{Leibniz-Institut  f\"ur Polymerforschung Dresden, Institut Theorie der Polymere, Dresden, 01069, Germany}
\affil[2]{Institut f\"ur Theoretische Physik II, Weiche Materie, Heinrich-Heine-Universit\"at D\"usseldorf, D\"usseldorf, 40225, Germany}
\affil[3]{Department de Physique, Universit\'e de Fribourg, CH-1700 Fribourg, Switzerland}
\affil[4]{Technische Universit\"at Dresden, Institut f\"ur Theoretische Physik, Dresden, 01069,  Germany}
\affil[*]{sharma@ipfdd.de}

%\keywords{Keyword1, Keyword2, Keyword3}

\begin{abstract}
%Brownian gyrator as a minimal micro heat engine is made by anisotropically trapping a single diffusing particle with thermal fluctuations of unequal strength along each position degree of freedom. The engine is based on a misalignment between the  principal axes of the anisotropic potential and the temperature axes, which are aligned with the coordinate axes, giving rise to a systematic torque generation due to the dissipation from the different thermostats.    Here we show that in contrast to previous studies a torque can be generated by applying a uniform magnetic field to the system even when the potential aligns with the temperature axes. In addition, we show that a nontrivial flux-free steady state can be created due to the Lorentz force when the misalignment exists. By applying a linear external nonconservative force we calculate the average mechanical power of the work done by the force. 
%derive the work out of the engine and the dissipated (extracted) heat into (from) the cold (hot) heat bath. %However, the efficiency of the engine needs more information. 
\end{abstract}

\begin{document}

\flushbottom
\maketitle
In the Supplemental Materials,  we first investigate the overdamped approximation of a Brownian magneto-gyrator: We take the low-mass approach to derive the proper overdamped Langevin equation describing the dynamics of such a microengine. We calculate the stationary-state quantities of the microengine,  such as the average exerted torque by the microengine on the confining potential and mechanical power, and study the tunability of the engine via the applied magnetic field. Next, we obtain the exact expressions for the efficiency of a magneto heat engine and the COP of a magneto-refrigerator in the underdamped regime, as discussed in the main text.

%take the low-mass approach to derive the overdamped Langevin equation. We then use this equation to study the model of a Brownian magneto heat pump as a tunable microengine (i.e., in the overdamped regime). Next, we investigate the stationary-state behavior of the microengine. We obtain the stationary-state quantities of the system such as the average exerted torque by the microengine on the confining potential and mechanical power. Finally, we calculate the exact expressions for the efficiency of a mesoscopic magneto heat engine and the COP of a magneto-refrigerator.  

%\section*{Brownian Magneto-Gyrator as a Tunable Microengine}

\section*{Low-mass Limit of Langevin Equation}
In the presence of a magnetic field, the overdamped Langevin equation can not be derived by setting the inertial term to zero. 
Anomalous fluxes has been observed in overdamped Brownian dynamics by setting the inertial term to zero in the velocity Langevin equation due to Lorentz force~\cite{vuijk2019anomalous}.
Here we take the low-mass approach~\cite{chun2018emergence} to derive the overdamped Langevin equation describing the dynamics of a charged Brownian particle with mass $m$ and charge $q$ subjected to a magnetic field $B$ in the $\hat{z}$ direction. The velocity Langevin equation for the position $\rr=(x, y)^\top$ and the velocity $\vv=(v_x, v_y)$ of the particle under the effect of the external force $\f(\rr)=(f_x(\rr), f_y(\rr))^\top$, can be written as
\begin{equation}
\label{underdampedlangevin}
m\dot{\vv} = \f(\rr(t))-\G\vv(t)+\xxi(t),
\end{equation}
where  $\xxi(t)=(\xi_x(t), \xi_y(t))^\top$ is the Gaussian white noise with zero mean and Dirac delta correlation $\langle\xxi(t)\xxi^\top(t')\rangle=2\gamma\T\delta(t-t')$ with $\gamma$ being the friction coefficient and $\T=\diag(T_c, T_h)$. The matrix $\G$ is defined as
\begin{equation}
\label{matrixG}
\G = \gamma\left( \begin{array}{cc}
1 & -\kappa \\
\kappa & 1 \\
\end{array}\right),
\end{equation}
with $\kappa=qB/\gamma$ being the diffusive Hall parameter which quantifies the strength of the Lorentz force relative to the frictional force. 
By omitting the transient term $e^{-\G t/m}\vv(0)$, which is negligible in the low-mass limit, the formal solution to Eq.~\eqref{underdampedlangevin} can be written as
\begin{equation}
\label{formalsolution}
\vv(t)=\frac{1}{m}\int_0^t\dif t' e^{-\frac{\G}{m}(t-t')}[\f(\rr(t'))+\xxi(t')].
\end{equation}
which yields the following stochastic integrodifferntial equation
\begin{equation}
\label{integrodifferential}
\dot{\rr}(t)=\frac{1}{m}\int_0^t\dif t' e^{-\frac{\G}{m}(t-t')}\f(\rr(t'))+\eeta_m(t),
\end{equation}
where $\eeta_m(t)$ is Gaussian non-white noise given as
\begin{equation}
\label{nonwhitenoise}
\eeta_m(t)=\frac{1}{m}\int_0^t\dif t' e^{-\frac{\G}{m}(t-t') }\xxi(t').
\end{equation}

The main purpose of this section is to derive the statistical properties of this noise and thereby the overdamped Langevin equation. From Eq.~\eqref{nonwhitenoise} it is clear that the first moment of the noise is zero, i.e., $\langle\eeta_m(t)\rangle=0$. The interesting property of the noise is encoded in its time correlation $\langle\eeta_m(t)\eeta_m^\top(s)\rangle=\boldsymbol{\Sigma}_m(t, s)$ in the low-mass limit. Here $\boldsymbol{\Sigma}_m(t, s)$ is the correlation matrix and the subscript indicates the mass dependence of the quantity. The correlation matrix can be written as 
\begin{equation}
\label{correlationmatrix1}
\boldsymbol{\Sigma}_m(t, s)=\frac{1}{m^2}\int_0^s\dif u\int_0^t\dif t' e^{-\frac{\G}{m}(t-t') }\langle\xxi(t')\xxi^\top(u)\rangle e^{-\frac{\G^\top}{m}(s-u)},
\end{equation}
where using the properties of the Gaussian white noise $\xxi(t)$, in the limit that $t\rightarrow\infty$ and $s\rightarrow\infty$ and for a finite $t-s$, it can be rewritten as
\begin{equation}
\label{correlationmatrix2}
\boldsymbol{\Sigma}_m(t,s)=\frac{2\gamma}{m^2}e^{-\frac{\G}{m} (t-s)}\lim_{s\rightarrow\infty}\int_0^s\dif u \left[e^{-\frac{\G}{m}(s-u)}\T e^{-\frac{\G^\top}{m}(s-u)}\right],
\end{equation}
if $t\geq s$ and 
\begin{equation}
\label{correlationmatrix3}
\boldsymbol{\Sigma}_m(t,s)=\frac{2\gamma}{m^2}\left[\lim_{t\rightarrow\infty}\int_0^t\dif t' e^{-\frac{\G}{m}(t-t')}\T e^{\frac{-\G^\top}{m}(t-t')}\right]e^{-\frac{\G^\top}{m}(s-t)},
\end{equation}
%\section*{Fokker-Planck equation} 
if $t<s$. The integrals in Eq.~\eqref{correlationmatrix2} and Eq.~\eqref{correlationmatrix3} can be easily calculated, which yield
\begin{equation}
  \boldsymbol{\Sigma}_m(t,s)=\begin{cases}
    \frac{1}{m}e^{-\frac{\G}{m} (t-s)} \CC, & \text{$ t \geq s$},\\
    \frac{1}{m}\CC e^{-\frac{\G^\top}{m}(s-t)}, & \text{$ t < s$}.
  \end{cases}
  \label{magneticfield}
\end{equation}
where the matrix $\CC$ is defined as
\begin{equation}
\label{matrixC}
\CC = \frac{1}{2(1+\kappa^2)}\left( \begin{array}{cc}
2T_c+\kappa^2(T_c+T_h) & \kappa(T_h-T_c) \\
\kappa(T_h-T_c) & 2T_h+\kappa^2(T_c+T_h) \\
\end{array}\right).
\end{equation}
Note that the cross-correlation components of the correlation matrix in Eq.~\eqref{matrixC} have their origin in the broken time-reversal symmetry of the magnetic field~\cite{chun2018emergence,vuijk2019anomalous}. 

The correlation matrix $\boldsymbol{\Sigma}_m$ depends on $s'\equiv t-s$ and satisfies $\boldsymbol{\Sigma}_m(-s')=\boldsymbol{\Sigma}_m^\top(s')$. The elements of the correlation functions are given by
 \begin{equation}
\label{matrixSigma}
 \boldsymbol{\Sigma}_m(s')= \frac{1}{m}\CC^{1-a} e^{-\frac{\gamma}{m}|s'|}\left( \begin{array}{cc}
\cos\left(\frac{qB}{m}s'\right) & \sin\left(\frac{qB}{m}s'\right) \\
-\sin\left(\frac{qB}{m}s'\right) & \cos\left(\frac{qB}{m}s'\right) \\
\end{array}\right)\CC^{a}.
\end{equation}
where $a=0$ if $s^\prime<0$ and $a=1$ otherwise. As $m$ decreases, the correlation functions become singular with diverging oscillation frequency $qB/m$, diverging amplitude proportional to $1/m$, and vanishing decay time $m/\gamma$. It can be shown that for any function $h(s')$ with no singularity around $s'=0$, the low-mass limit yields

\begin{eqnarray}
     \label{lowmass1}
	\lim_{m\rightarrow 0}\int_0^\infty \dif s' \,\, h(s')\boldsymbol{\Sigma}_m(s')=h(0)\G^{-1}\CC,\\
	\label{lowmass2}
	\lim_{m\rightarrow 0}\int_{-\infty}^0 \dif s' \,\, h(s')\boldsymbol{\Sigma}_m(s')=h(0)\CC(\G^{-1})^\top.
\end{eqnarray} 

Therefore, the correlation matrix in the low-mass limit, i.e., $m\rightarrow 0$ can be written as
\begin{equation}
\label{correlationmatrixfinal}
\boldsymbol{\Sigma}(s')\equiv\lim_{m\rightarrow 0}\boldsymbol{\Sigma}_m(s')= \G^{-1}\CC\delta_{+}(s')+ \CC(\G^{-1})^\top\delta_{-}(s'),
\end{equation}
where the notations $\delta_{\pm}(s')$ indicate the different Dirac delta functions which are zero for $s'\neq 0$ while $\int_0^\infty \dif s'\delta_+(s')=\int_{-\infty}^0\dif  s'\delta_-(s')=1$ and $\int_0^\infty \dif s'\delta_-(s')=\int_{-\infty}^0\dif s'\delta_+(s')=0$.

To complete the derivation of the overdamped Langevin equation we now evaluate the first term on the right hand side of Eq.~\eqref{integrodifferential}. Using the change of variable from $t'$ to $s'=t-t'$ it can be written as
\begin{equation}
\label{firstterm}
\frac{1}{m}\int_0^t\dif t' e^{-\frac{\G}{m}(t-t')}\f(\rr(t'))=\left[\int_0^t\dif s'\,\,\,\boldsymbol{\Sigma}_m(s')\f(\rr(t-s'))\right]\CC^{-1},
\end{equation}
where using Eq.~\eqref{lowmass1}, in the low-mass limit converges to $\G^{-1}\f(\rr(t))$. Thus the Langevin equation in the low-mass limit corresponding to the overdamped dynamics can be written as
\begin{equation}
\label{overdampedLE}
\dot{\rr}(t)=\G^{-1}\f(\rr(t))+\eeta(t),
\end{equation}
where $\eeta(t)$ is Gaussian non-white noise with zero mean and the time correlation as in Eq.~\eqref{correlationmatrixfinal}.

\begin{figure}[t]
\centering
\includegraphics[width=0.5\linewidth]{figure01a}
\includegraphics[width=0.475\linewidth, height=7cm]{figure01b}
\caption{  A single charged Brownian particle, steered by an external magnetic field $B$ and trapped in a rotationally asymmetric potential, performs gyration when subjected to different thermal noises from cold, $T_c$, and hot, $T_h$, heat baths coupled to its $x$ and $y$ degrees of freedom, respectively. Figure (a) schematically depicts the diffusion of the particle in a magnetic field in the $z$ direction and under the influence of the harmonic potential $V(x,y)=k[(x^2+\alpha y^2)/2+uxy]$ with the parameters $k$ and $\alpha$. Here $u$ is the potential coupling parameter which correlates the spatial degrees of freedom. The generic scenario is considered when the principal axes of the potential, shown by dashed lines, are misaligned with the temperature axes, namely if $u\neq 0$. Note that this condition (i.e., $u\neq 0$) is not necessary for the system to operate as an engine, yet is considered for further tunability of the microengine. (b) Schematic illustration of the average torque generated by the microengine for different values of the diffusive Hall parameter $\kappa=qB/\gamma$ and the potential coupling parameter $u$. The solid, oblique line represents zero torque in the system corresponding to stalled engine. The reddish region corresponds to the clockwise torque while the gray region depicts the counterclockwise torque. As can be seen, a torque can be generated even if the potential coupling parameter is zero upon applying an external magnetic field. The generated torque can be tuned by varying the magnitude and direction of the applied magnetic field: it changes the direction by reversing the direction of the applied magnetic field if $u=0$ or by additionally changing the sign of the potential coupling parameter when it is not zero. Note that latter can be also done by further tunability of the applied magnetic field for a fixed $u$.}

\label{fig1:schema}
\end{figure}
\section*{Brownian Magneto-Gyrator as a Tunable Microengine}
A Brownian magneto-gyrator is made of a single Brownian particle with charge $q$ steered by a constant magnetic field $B$ and trapped by the conservative force $\F_c=-\nabla V(\rr)$. We consider a magnetic field in the $\hat{z}$ direction such that $\B = B\hat{z}$ and hence the particle's motion along this direction is not affected by the applied magnetic field. Consequently, we effectively have a two-dimensional system with the particle's position $\rr=(x, y)^\top$ where $\top$ indicates the transpose.   
The thermal fluctuations of unequal strength, proportional to the cold and hot heat bath temperatures $T_c$ and $T_h$, are supplied along the two Cartesian coordinates $\hat{x}$ and $\hat{y}$, respectively. Here $V(\rr)=\frac{1}{2}\rr^\top\cdot\U\cdot \rr$ is the potential where the matrix $\U$ is given as
\begin{equation}
\label{matrixU}
\U = k\left( \begin{array}{cc}
1 & u \\
u & \alpha \\
\end{array}\right), 
\end{equation}
where $k$ is the stiffness of the potential, $\alpha$ is a dimensionless parameter quantifying the difference in the stiffness in the $x$ and $y$ directions, and $u$ is the potential coupling parameter which correlates the spatial degrees of freedom. Note that the stability condition implies that $|u|< \sqrt{\alpha}$. The eigenvectors of the matrix $\U$ correspond to the principal axes of the potential, which are shown by dashed lines in Fig.~\ref{fig1:schema} (a). This figure shows a schematic of the system where the alignment of the principal axes of the potential with the temperature axes is specified by the potential coupling parameter. In the case of $u=0$ the principal axes of the potential align with the temperature axes. The overdamped dynamics of the system can be described by the following Langevin equation 
\begin{equation}
\label{langevinequaion:overdamped}
\dot{\rr}(t) = -\A\rr(t)+\xchi(t), 
\end{equation}
where $\A=\G^{-1}\U$ and $\xchi(t)$ is Gaussian nonwhite noise~\cite{chun2018emergence} with zero mean and time correlation as in Eq.~\eqref{correlationmatrixfinal}. % with $\G$ being
%\begin{equation}
%\label{matrixG}
%\G = \gamma\left( \begin{array}{cc}
%1 & -\kappa \\
%\kappa & 1 \\
%\end{array}\right),
%\end{equation}
%where $\gamma$ is the friction coefficient and  $\kappa=qB/\gamma$ is the diffusive Hall parameter which quantifies the strength of the Lorentz force relative to the frictional force. 
%Here $\xchi(t)$ is Gaussian nonwhite noise~\cite{chun2018emergence} with zero mean and time correlation as in Eq.~\eqref{correlationmatrixfinal}.
%\begin{equation}
%\label{timecorrelation}
%\langle\xchi(t)\xchi^{\top}(t')\rangle= \G^{-1}\CC\delta_{+}(t-t')+ \CC(\G^{-1})^\top\delta_{-}(t-t'),
%\end{equation}
%where the notations $\delta_{\pm}(s')$ with $s'=t-t'$ indicate the modified Dirac delta functions which are zero for $s'\neq 0$ while $\int_0^\infty \dif s'\delta_+(s')=\int_{-\infty}^0\dif s'\delta_-(s')=1$ and $\int_0^\infty \dif s'\delta_-(s')=\int_{-\infty}^0\dif s'\delta_+(s')=0$. The matrix $\CC$ is given as %\cmnt{all matrices we denote by non-bold symbols and all vectors by bold symbols? Then we need to change G.}
%\begin{equation}
%\label{matrixC}
%\CC = \frac{1}{2(1+\kappa^2)}\left( \begin{array}{cc}
%2T_c+\kappa^2(T_c+T_h) & \kappa(T_h-T_c) \\
%\kappa(T_h-T_c) & 2T_h+\kappa^2(T_c+T_h) \\
%\end{array}\right).
%\end{equation}

%Note that the cross-correlation components of the correlation matrix in Eq.~\eqref{timecorrelation} have their origin in the broken time-reversal symmetry of the magnetic field~\cite{chun2018emergence,vuijk2019anomalous}. 


\subsection*{Stationary-State Properties}
\begin{figure}[t]
\centering
\includegraphics[width=1\linewidth]{figure02}
\caption{The left figure shows a panel of the stationary-state probability density of the particle's position and fluxes in the system for different values of the diffusive Hall parameter $\kappa$ and the potential coupling parameter $u$. The results are shown from Eq.~\eqref{PDF} in the stationary state for the probability density and substituting that equation into Eq.~\eqref{flux} for the fluxes with $T_h=4T_c=4.0$ and $\alpha=4.0$. The figures correspond to $u=-0.5, 0.0, 0.5$ from  left to right for $\kappa=-0.6$ (top) and  $\kappa=0.6$ (bottom). The Brownian magneto-gyrator operates as a microengine even in the absence of the parameter $u$ due to the exerted torque by the particle on the potential. In this case, the direction of the gyration can be reversed by reversing the direction of the magnetic field resulting in a rotation in the probability density, shown in the middle figures. If $u\neq 0$ this can be done by additionally changing the sign of the potential coupling parameter or  tuning further the applied magnetic field. The direction of the fluxes are shown by arrows and the magnitude is color coded.
The right panel shows that the microengine can be stalled, which means that $\J=\0$, either
 if $\alpha=1$ and $u=0$ or $\kappa=2u/(\alpha-1)$ and $u\neq 0$ since the particle does not exert any torque. These cases are shown at the top and bottom, respectively. Note that the particle does not perform any gyrating motion if its motion is not correlated, namely when $u$ and $\kappa$ are zero. }
\label{fig:PDF}
\end{figure}

%Now we take the Fokker-Planck equation approach in order to provide a probabilistic description of the system. 
In order to investigate the stationary-state properties of the system, we use the generalized Fokker-Planck equation corresponding to the overdamped Langevin equation~\eqref{langevinequaion:overdamped}, which can be written as

\begin{equation}
\label{FPE}
	\frac{\partial \rho(\rr, t)}{\partial t} = -\nabla\cdot\J(\rr, t),
\end{equation}
where $\rho(\rr, t)$ is the probability density of finding the particle at position $\rr$ at time $t$ and $\J(\rr, t)$ is the probability flux, given as
\begin{equation}
\label{flux}
\J(\rr, t)= -\D\nabla \rho(\rr, t)- (\A\cdot\rr)\rho(\rr, t),
\end{equation}
where the matrix $\D$ can be derived from Eq.~\eqref{langevinequaion:overdamped}~\cite{chun2018emergence} or alternatively by using a first-priciple approach ~\cite{abdoli2020correlations}, which is given as
\begin{equation}
\D
= \frac{1}{\gamma} 
\left( \begin{array}{cc}
\frac{T_c + \kappa^2 T_h}{\left(1+\kappa^2\right)^2} & \frac{\kappa(T_h - T_c)}{\left(1+\kappa^2\right)^2}+\frac{\kappa(T_c+T_h)}{2(1+\kappa^2)} \\
\frac{\kappa(T_h - T_c)}{\left(1+\kappa^2\right)^2}- \frac{\kappa(T_c+T_h)}{2(1+\kappa^2)} & \frac{T_h + \kappa^2 T_c}{\left(1+\kappa^2\right)^2} \\
\end{array}\right).
\label{diffusionmatrix}
\end{equation}
Surprisingly, despite the overdamped motion of the particles, the dynamics preclude a purely diffusive description as can be seen in the unusual structure of $\D$ above. In contrast to typical diffusion tensors, $\D$ has antisymmetric elements. These give rise to additional rotational fluxes in the system, which are perpendicular to the typical diffusive fluxes~\cite{vuijk2019anomalous, vuijk2019effect, abdoli2020nondiffusive, abdoli2020stationary, abdoli2021stochastic}. One of the known features of Lorentz forces is the reduction of the diffusion by a factor of $1/(1+\kappa^2)$, as can be seen from Eq.~\eqref{diffusionmatrix} for $T_c=T_h=T$.  Note that under Lorentz force, the diffusion is spatially correlated as can be seen from the symmetric, off-diagonal terms which are non-zero only for $T_h \neq T_c$ and $\kappa \neq 0$. Throughout this work we set the Boltzmann constant $k_B$ to unity.

%The broken time-reversal symmetry gives rise to unusual dynamics in colloidal systems~\cite{ jimenez2006brownian,filliger2007kramers, jimenez2008brownian, tothova2010hydrodynamic, lisy2013brownian, jimenez2013brownian,chun2018emergence, vuijk2019anomalous, matevosyan2021non}. Surprisingly, despite the overdamped motion of the particles, the dynamics preclude a diffusive description. Lorentz forces give rise to additional rotational (or Lorentz) fluxes in the system, which are perpendicular to the typical diffusive fluxes~\cite{vuijk2019anomalous, vuijk2019effect, abdoli2020nondiffusive, abdoli2020stationary, abdoli2021stochastic}. These fluxes fundamentally alter the dynamics of the system without affecting its equilibrium properties~\cite{van1921problemes}. 

 Equation~\eqref{FPE} is a linear Fokker-Planck equation whose Gaussian solution can be written as
%\begin{equation}
%\label{PDF}
%\rho(x, y) = \frac{1}{Z}e^{-\frac{1}{2}(\mu_1x^2+\mu_2y^2+\mu_3xy)},
%\end{equation}
\begin{equation}
\label{PDF}
\rho(\rr, t) = \frac{1}{Z}e^{-\frac{1}{2}(\rr^\top\cdot\M_t^{-1}\cdot\rr)},
\end{equation}
where $Z=2\pi\sqrt{\Det(\M_t)}$  is the normalization factor implying the probability density function with total probability of one and the covariance matrix $\M_t$ satisfies the following Lyapunov equation
\begin{equation}
\label{methods:lyapunovequation}
\frac{\dif\M_t}{\dif t}=\A\M_t+\M_t\A^\top+2\D_s
\end{equation}
where $\D_s=(\D+\D^\top)/2$ is the usual (symmetric) diffusion matrix. Equation~\eqref{methods:lyapunovequation} can be solved for the stationary state by setting $\dif\M_t/\dif t$ to zero. This gives the stationary-state covariance matrix, denoted by $\M$, as
\begin{equation}
\label{methods:covariancematrix}
\M = \left( \begin{array}{cc}
\frac{\alpha^2(1+\kappa^2)T_c+u^2(T_h-T_c)+\alpha (T_c+\kappa^2T_h)}{k(1+\alpha)(\alpha-u^2)(1+\kappa^2)} & -\frac{\kappa (\alpha-u^2)(T_h-T_c)-u(1+\kappa^2)(\alpha T_c+T_h)}{k(1+\alpha)(\alpha-u^2)(1+\kappa^2)}  \\
-\frac{\kappa(\alpha-u^2)(T_h-T_c)-u(1+\kappa^2)(\alpha T_c+T_h)}{k(1+\alpha)(\alpha-u^2)(1+\kappa^2)} & \frac{(1+\kappa^2)T_h-u^2(T_h-T_c)+\alpha(\kappa^2T_c+T_h)}{k(1+\alpha)(\alpha-u^2)(1+\kappa^2)} \\
\end{array}\right), 
\end{equation}
%where $\kappa=qB/\gamma$ is the diffusive Hall parameter which quantifies the strength of the Lorentz force relative to the frictional force. Plugging Eq.~\eqref{methods:covariancematrix} into Eq.~\eqref{methods:gaussiansolution} gives 

The stationary-state probability density can be written as $\rho(\rr) = e^{-\frac{1}{2}(\rr^\top\cdot\M^{-1}\cdot\rr)}/Z$ where $\M^{-1}$ and $Z$ are given as 
%\begin{equation}
%\label{PDFstationary}
%\rho(\rr) = \frac{1}{Z}e^{-\frac{1}{2}(\rr^\top\cdot\M^{-1}\cdot\rr)},
%\end{equation}


%\bigskip \bigskip \bigskip


%where $\M$ is the stationary-state covariance matrix and $Z$ is the normalization factor implying the probability density function with total probability of one, given as
\begin{equation}
\label{normalization}
Z=\frac{2\pi\sqrt{[u(T_h-T_c)+\kappa(\alpha T_c+T_h)]^2+(1+\alpha)^2T_cT_h}}{k(1+\alpha)\sqrt{(1+\kappa^2)(\alpha-u^2)}}, \hspace{1cm} \text{and} \hspace{1cm}
\M^{-1}=\left( {\begin{array}{*{20}c}
   \mu_1 & \mu_3  \\
   \mu_3 & \mu_2  \\    
 \end{array} } \right),
\end{equation}
where 
%\rho(x, y) = \frac{(1+\alpha)\sqrt{(1+\kappa^2)(\alpha^2-u^2)}e^{-\frac{1}{2}(\mu_1x^2+\mu_2y^2+\mu_3xy)}}{2\pi\sqrt{[u(T_h-T_c)+\kappa(\alpha T_c+T_h)]^2+(1+\alpha)^2T_xT_y}}
%[u(T_h-T_c)+\kappa(\alpha T_c+T_h)]^2+(1+\alpha)^2T_xT_y
\begin{eqnarray}
     \label{mu1}
	\mu_1 &=& k(1+\alpha)\frac{(1+\kappa^2)T_h-u^2(T_h-T_c)+\alpha(\kappa^2 T_c+T_h)]}{[u(T_h-T_c)+\kappa(\alpha T_c+T_h)]^2+(1+\alpha)^2T_cT_h},\\\label{mu2}	
	\mu_2 &=& k(1+\alpha)\frac{\alpha^2(1+\kappa^2)T_c+u^2(T_h-T_c)+\alpha(T_c+\kappa^2T_h)}{[u(T_h-T_c)+\kappa(\alpha T_c+T_h)]^2+(1+\alpha)^2T_cT_h},\\
	\label{mu3}
	\mu_3 &=& k(1+\alpha)\frac{u(\alpha T_c+T_h)(1+\kappa^2)-\kappa(\alpha-u^2)(T_h-T_c)}{[u(T_h-T_c)+\kappa(\alpha T_c+T_h)]^2+(1+\alpha)^2T_cT_h}.
\end{eqnarray} 

Note that in the case of an isotropic potential with $u=0$ and $\alpha=1$ Eqs.~\eqref{normalization} to \eqref{mu3} reduce to the results reported in Ref.~\cite{abdoli2020correlations} .  The probability fluxes in the stationary state of the system can be determined by substituting the stationary-state probability density into Eq.~\eqref{flux}.

In the left panel of Fig.~\ref{fig:PDF} we show the stationary-state probability density of the position of the particle and the fluxes in the system for different values of the diffusive Hall parameter $\kappa$ and the potential coupling parameter $u$.  The fluxes are shown on top of the probability density. The Brownian magneto-gyrator does not essentially require a positional coupling due to the potential, as shown in the middle figures for $u=0$. However, the generic scenario in the presence of the potential coupling parameter, corresponding to $u\neq 0$, is also investigated that is shown in the left and right figures. The direction of the gyration of the microengine can be reversed by tuning the strength and direction of the applied magnetic field. In doing so, one rotates the probability density as well.  

%The particle performs gyrating motion around the potential energy minimum as long as the system is rotational asymmetric about the origin. 

\subsection*{Exerted Torque}

%In order to have more insight into the magneto-gyrator and its tunability 

To obtain a simple scalar quantifier for the strength of the gyrating current field, which further emphasizes the tunability of the magneto-gyrator, we investigate the average torque on the potential. The average exerted torque by the particle on the potential $V$, denoted by $\langle\tau\rangle$, can be written as  



\begin{equation}
\label{torque:definition}
\langle\tau\rangle = \int \rho(\rr) (\rr\times\F_c) \dif\rr,
\end{equation}
which is exactly equal to the opposite torque which the particle exerts via the friction forces on the thermal environment.
% or any physical objects containing the baths.
The substitution of Eq.~\eqref{PDF} for the stationary-state into Eq.~\eqref{torque:definition} gives the average exerted torque on the potential by the particle as %(see Sec.~\ref{methods})
\begin{equation}
\label{torque}
\langle\tau\rangle = \frac{(T_h-T_c)[2u+(1-\alpha)\kappa]}{(1+\alpha)(1+\kappa^2)}.
\end{equation}
\begin{figure}[t]
\centering
\includegraphics[width=0.95  \linewidth]{figure03}
%\includegraphics[width=0.5\linewidth]{g1936}
\caption{ Average exerted torque by the microengine as a function of (a) diffusive Hall parameter $\kappa$ for different values of $\alpha$ and (b) the parameter $\alpha$ for different values of $\kappa$ with $T_h=4T_c=4.0$. The main figures show the results for $u=0.5$ while the insets depict those for $u=0$. As shown in (a), the average torque from Eq.~\eqref{torque} shows a maximum at an optimal $\kappa$, which can be reversed by tuning the strength and direction of the applied magnetic field. (b) In the generic scenario when $u\neq 0$, with increasing $\alpha$ the average torque increases if $\kappa<-0.5$ and decreases if $\kappa>-0.5$ till the exerted torque changes the direction resulting in the reversal of the gyration. It is independent of $\alpha$ if $\kappa=-0.5$. For $u=0$ with increasing $\alpha$, $\langle\tau\rangle$ increases if $\kappa<0$, decreases if $\kappa>0$ till the torque changes the sign, and is zero when $\kappa=0$. In addition, there is no torque if $\alpha=1$ and $\kappa\neq 0$, which correspond to the stalled microengine. }
\label{fig:torque}
\end{figure}
 

%For the former we have studied the system in Ref.~\cite{abdoli2020correlations} without having any clear reason for the absence of the fluxes. Here we find that this is due to the cessation  of the average torque. The same can be seen for the latter from Eq.~\eqref{torque}. The results for these three cases are shown in Fig.~\ref{fig:PDF}, respectively. 


Note that Eq.~\eqref{torque} implies that even if $u=0$ there still exists a systematic torque giving rise to a nonequilibrium stationary state carrying fluxes, corresponding to an operating magneto-gyrator. This is due to the applied magnetic field ($\kappa\neq 0$) and the asymmetry of the potential about the origin ($\alpha\neq 1$). In this case, the torque exerted by the microengine on the potential can be reversed by reversing the direction of the applied magnetic field. This results in the reversal of the gyration, which is shown in the middle figures of the left panel of Fig.~\ref{fig:PDF}. In the presence of the potential coupling parameter, this can be done by additionally changing the sign of the potential coupling parameter or tuning further the applied magnetic field, shown in the left and right figures of Fig.~\ref{fig:PDF}. 
The average exerted torque is zero and consequently the microengine is stalled in two cases: (i) in the absence of the potential coupling parameter if $\alpha=1$  and (ii) if the potential is anisotropic with a nonzero $u$, but $\kappa=2u/(\alpha-1)$. The absence of fluxes in these cases are due to the cessation of the average torque which correspond to stalled microengines. The results for the  two cases are shown in the right panel of Fig.~\ref{fig:PDF}.


In Fig.~\ref{fig:torque} (a), we show that there exists a maximum torque exerted by the microengine at an optimal magnetic field. The average torque can be reversed by reversing the direction of the magnetic field when $u=0$, which is shown in the inset. In the presence of the potential coupling parameter, shown in the main figure, this can be done by tuning further the applied magnetic field. 
Figure~\ref{fig:torque} (b) shows how the exerted torque by the microengine varies by tuning the stiffness of the potential for different values of the diffusive Hall parameter.
The main figure represents the results for the Brownian magneto-gyrator without the potential coupling parameter while those for the generic scenario, i.e., in the presence of $u$, are shown in the inset. From Eq.~\eqref{torque} it is clear that the average torque is independent of the parameter $\alpha$ if $u+\kappa=0$, which is shown in the main panel. The exerted torque decreases (increases) with increasing $\alpha$ if $\kappa>-u$ ($\kappa<-u$) till it changes the direction, corresponding to the change of the direction of the gyration. A similar explanation holds when $u=0$. However, there is no torque if $\kappa=0$ or $\alpha=1$ and consequently the microengine is stalled.  
   
%\newpage

\subsection*{Mechanical Power}

Now we consider the microengine in the presence of a load in order to determine the delivered mechanical power by the engine. To do this, we apply a linear external nonconservative force ($\nabla\times\F_{nc}\neq\0$) of the form $\F_{nc} = \epsilon (-y, x)^\top$ with a parameter $\epsilon$, yielding a torque in the $z$ direction, whose sign is chosen such that the resulting torque in the $z$ direction is opposed to $\langle\tau\rangle$ in Eq.~\eqref{torque}. The goal is to calculate the average power of the work done by this force. 
The average extracted mechanical power $P=-\langle\F_{nc}^\top.\dot{\rr}\rangle$ in the stationary state can be rewritten as $P=-\epsilon\langle xv_y-yv_x\rangle$. Therefore, one needs to calculate  the stationary-state position-velocity correlation matrix, that is $\lim_{t\rightarrow\infty}\langle\rr(t)\dot{\rr}^\top(t)\rangle$. Under this condition, the motion of a particle with mass $m$ and velocity $\vv=(v_x, v_y)^\top$ can be described by the generalised Langevin equation written as $m\dot{\vv} = -\G\cdot\vv -\nabla V(\rr)+\F_{nc}+\xxi(t)$ where $\xxi(t) =(\xi_x, \xi_y)^\top$ is Gaussian white noise with zero mean and time correlation $\langle\xxi(t)\xxi^{\top}(t')\rangle= 2\gamma\diag(T_c, T_h)\delta(t-t')$ where $\diag(T_c, T_h)$ is a diagonal matrix. In this case, the conservative and nonconservative forces can be combined and the overdamped Langevin equation can be written in the form of  Eq.~\eqref{langevinequaion:overdamped} with a modified matrix $\U_l$ under loading, defined as


%To do this, one uses the overdamped Langevin equation~\eqref{langevinequaion:overdamped} for a modified matrix $\U_l$ under loading as
 \begin{equation}
\label{matrixU:modified}
\U_l = k\left( \begin{array}{cc}
1 & u+\epsilon^\prime \\
u-\epsilon^\prime & \alpha \\
\end{array}\right), 
\end{equation}
where $\epsilon^\prime=\epsilon/k$ is a dimensionless parameter. Using the equation of motion in Eq.~\eqref{langevinequaion:overdamped} it can be shown that 
\begin{equation}
\label{correlationmatrix}
\langle\rr\dot{\rr}^\top\rangle = -\langle\rr\rr^\top\rangle\A_l^\top + \langle\rr\xchi^\top\rangle,
\end{equation}
where $\A_l=\G^{-1}\U_l$. Note that by applying a prescribed external vortex
flow field such as a rotating bucket to an underdamped Brownian particle one can induce similar terms to the nonconservative force in the overdamped limit~\cite{liebchen2019optimal}. Moreover, it is known that in experiments using optical tweezers, optical scattering forces actually generate a nonconservative component to the overall force  exerted by the trap and do not simply furnish a simple potential~\cite{mangeat2019role}. Here we take a linear nonconservative force for the convenience of analytical calculations.
The position correlation and the position-noise correlation can be written as
\begin{eqnarray}
     \label{methods:positioncorrelation}
	\langle\rr\rr^\top\rangle &=& \lim_{t\rightarrow\infty}\int_0^t\dif t'\int_0^t\dif t''e^{-\A_l(t-t')}\langle\xchi(t')\xchi^\top(t'')\rangle e^{-\A_l^\top(t-t'')},\\\label{methods:positionnoisecorrelation}	
	\langle\rr\xchi^\top\rangle &=& \lim_{t\rightarrow\infty}\int_0^t\dif t' e^{-\A_l(t-t')}\langle\xchi(t')\xchi^\top(t)\rangle.
\end{eqnarray} 
The substitution of the properties of the Gaussian nonwhite noise from Eq.~\eqref{timecorrelation} into Eq.~\eqref{methods:positioncorrelation} and Eq.~\eqref{methods:positionnoisecorrelation} gives
\begin{eqnarray}
     \label{methods:2positioncorrelation}
	\langle\rr\rr^\top\rangle &=& \int_0^\infty\dif s e^{-\A_l s}[\G^{-1}\CC+\CC(\G^{-1})^\top]e^{-\A_l^\top s},\\\label{methods:2positionnoisecorrelation}	
	\langle\rr\xchi^\top\rangle &=& \CC(\G^{-1})^\top ,
\end{eqnarray} 
where Eq.~\eqref{methods:2positioncorrelation} is the solution to the following Lyapunov equation
\begin{equation}
\label{methods:lyapunov2}
\A_l\langle\rr\rr^\top\rangle+\langle\rr\rr^\top\rangle\A_l^\top=-\Xi,
\end{equation}
where $\Xi=\G^{-1}\CC+\CC(\G^{-1})^\top$. The stationary-state position-velocity matrix can be calculated and written as
 \begin{equation}
\label{correlationmatrix}
\langle\rr\dot{\rr}^\top\rangle = \left( \begin{array}{cc}
0 & \frac{[u-\frac{\kappa}{2}(\alpha-1)](T_h-T_c)+\epsilon^\prime(1+\kappa^2)(T_c+T_h)}{\gamma(1+\kappa^2)(1+\alpha-2\epsilon^\prime\kappa)} \\
-\frac{[u-\frac{\kappa}{2}(\alpha-1)](T_h-T_c)+\epsilon^\prime(1+\kappa^2)(T_c+T_h)}{\gamma(1+\kappa^2)(1+\alpha-2\epsilon^\prime\kappa)} & 0 \\
\end{array}\right),
\end{equation}
and the average mechanical power reads
\begin{equation}
\label{methods:mechanicalpaower}
P=2k\epsilon^\prime\frac{[u-\frac{\kappa}{2}(\alpha-1)](T_h-T_c)+\epsilon^\prime(1+\kappa^2)(T_c+T_h)}{\gamma(1+\kappa^2)(1+\alpha-2\epsilon^\prime\kappa)},
\end{equation}
which can be rewritten as 
\begin{equation}
\label{mechanicalpower}
P = \frac{h}{1-\lambda \frac{\epsilon}{\epsilon_s}}\frac{\epsilon}{\epsilon_{s}}\left(1-\frac{\epsilon}{\epsilon_s}\right),
\end{equation}
where $h=2\epsilon_s^2(T_c+T_h)/(k\gamma(1+\alpha))$, $\lambda=2\kappa\epsilon_s/(k(1+\alpha))$, and $\epsilon_s$ is the stall parameter which quantifies the maximum strength of the nonconservative force that one can apply to the engine without driving it in the inverse direction. Hence, the system operates as a heat engine, delivering mechanical work for $0 < \frac{\epsilon}{\epsilon_s} < 1$. The stall parameter is given as
\begin{equation}
\label{stallparameter}
\epsilon_s=-k\frac{\left (u-\frac{\kappa}{2}(\alpha-1)\right)\eta_c}{(2-\eta_c)(1+\kappa^2)},
\end{equation}
where $\eta_c=1-T_c/T_h$ is the Carnot efficiency. Having demonstrated in details the tunability of a Brownian magneto-gyrator in the presence of a rotated anisotropic potential, we focus from this point onwards on pure magneto-gyration, setting $u=0$.
\begin{figure}[t]
\centering
%\includegraphics[width=0.45\linewidth]{torque_main}
\includegraphics[width=0.8\linewidth]{figure04}
%\includegraphics[width=0.5\linewidth]{g1936}
\caption{The power of a Brownian magneto-gyrator with $T_h=4T_c=4.0$ and $u=0.0$. The main panel shows the scaled power $P/P_{max}$ as a function of $\epsilon/\epsilon_s$. The mechanical power vanishes when $\epsilon$ is zero or equal to the stall value. In addition, it has a maximum at $\epsilon/\epsilon_s=(1-\sqrt{1-\lambda})/\lambda$ where $\lambda=2\kappa\epsilon_s/(k(1+\alpha))$ which can be shifted by tuning the magnetic field. For the upper inset we use Eq.~\eqref{stallparameter} to plot the stall parameter with respect to the diffusive Hall parameter with $\alpha=4.0$. The stall parameter approaches zero with increasing magnitude of the applied magnetic field. In the lower inset, we show how the maximum power varies with respect to the diffusive Hall parameter for different values of $\alpha$.  $P_{max}$ is plotted in units of $P_0=k(T_c+T_h)/\gamma$. There is an optimal value of $\kappa$ for which the maximum power takes its highest value. The lines show the theoretical predictions and the symbols are simulation results.
}

\label{fig:work}
\end{figure}
%is defined as 
%\begin{equation}
%\label{maximumpower}
%P_{max}=\epsilon_s^2\frac{T_c+T_h}{2\gamma(2\epsilon\kappa+\alpha+1)}.
%\end{equation}

The main panel of Fig.~\ref{fig:work} shows the scaled power $P/P_{max}$ delivered by the magneto-gyrator as a function of scaled loading $\epsilon/\epsilon_s$. As expected, the mechanical power vanishes when the engine is unloaded, i.e., $\epsilon=0$ or when the engine is stalled at $\epsilon=\epsilon_s$. The engine operates at maximum power at an intermediate loading $\epsilon/\epsilon_s=(1-\sqrt{1-\lambda})/\lambda$ which can be shifted by tuning the applied magnetic field. As can be seen in the inset to Fig.~\ref{fig:work}, there exists an optimal value of $\kappa$ for which the microengine delivers the maximum mechanical power. Qualitatively, this can be understood as follows. In the limit of vanishing magnetic field, there is no spatial correlation and hence no gyration. In the opposite case of very large magnetic field, the correlations vanish due to the large reduction in the diffusion coefficient of the particle. 






\section*{Efficiency of a Mesoscopic Brownian Heat Engine}
%The exchange of energy between a Brownian particle and its surrounding environment becomes stochastic at the microscale and yet work, heat, and efficiency can be defined through the framework of stochastic thermodynamics~\cite{seifert2012stochastic}. In this section, we aim to determine such quantities. 
In this section, we calculate the stationary-state covariance matrix whose elements are needed to measure the efficiency of the heat engine and the COP of a magneto-gyrator. 
The Langevin equation~\eqref{underdampedlangevin} for the particle under a load of the form of the nonconservative force $\F_{nc}=\epsilon(-y, x)^\top$ with the parameter $\epsilon$ and in the harmonic potential $\V(x, y)=k(x^2+\alpha y^2)/2 $ can be rewritten as 
\begin{equation}
\dot{\qq}(t) =  -\F\qq(t) + \tilde{\xxi}(t),
\label{langevinequation:q}
\end{equation}
where $\qq(t)=(x(t), y(t), v_x(t), v_y(t))^\top$ and $\tilde{\xxi}(t)=(0, 0, m^{-1}\xi_x(t), m^{-1}\xi_y(t))^\top$ is Gaussian white noise with zero mean and time correlation $\langle\tilde{\xxi}(t)\tilde{\xxi}^{\top}(t')\rangle= (2\gamma/m^2)\T\delta(t-t')$ where $\gamma$ is the constant friction coefficient. Here $k$ is the stiffness of the potential, $\alpha$ is a dimensionless parameter, and $\T=\diag(0, 0, T_c, T_h)$ is a diagonal  matrix. Note that $\T$ is a two-dimensional matrix in Eq.~\eqref{underdampedlangevin}. The matrix $\F$ is defined as
\begin{equation}
\label{methods:matrixF}
\F = \frac{1}{m}\left( \begin{array}{cc}
\0 & -m\1 \\
\U_l & \G \\
\end{array}\right), 
\end{equation}
where $\1$ is the identity matrix and
\begin{equation}
\label{methods:matrixG}
\G = \gamma\left( \begin{array}{cc}
1 & -\kappa \\
\kappa & 1 \\
\end{array}\right), \hspace{1cm}
\U_l = k\left( \begin{array}{cc}
1 & \epsilon^\prime \\
-\epsilon^\prime & \alpha \\
\end{array}\right). 
\end{equation}
where $\epsilon^\prime=\epsilon/k$ is a dimensionless parameter. In order to calculate the efficiency of the microengine the calculation of the average rate of the heat out of the hot bath, $ \langle\,\dbar Q_c/dt\rangle$ and the average mechanical power, $P$ is needed, which can be determined by the steady-state covariance matrix $\Cov=\lim_{t\rightarrow\infty}\Cov(t)$, where $\Cov(t)=\langle\qq(t)\qq^\top(t)\rangle$. The change in the covariance matrix in the time interval $\dif t$ is given by
   \begin{equation}
\dif\Cov(t)=-[\F\Cov(t)+\Cov(t)\F^\top]\dif t +\int_t^{t+\dif t}\dif t^\prime\int_t^{t+\dif t}\dif t^{\prime\prime}\langle\tilde{\xxi}(t')\tilde{\xxi}^\top(t^{\prime\prime})\rangle,
\label{covariancematrixt}
\end{equation}
which using the property of the noise $\tilde{\xxi}$, the time evolution of the covariance matrix can be written as
   \begin{equation}
\frac{\dif\Cov(t)}{\dif t}=-\F\Cov(t)-\Cov(t)\F^\top+\frac{2\gamma}{m^2}\T.
\label{covariancematrixt}
\end{equation}

The steady-state covariance matrix can be calculated by setting $\dif\Cov(t)/\dif t$ to zero, which consists of the steady-state position-position, position- and velocity-velocity correlations. The solution to Eq.~\eqref{covariancematrixt} gives the steady-state covarinace matrix with the following elements

\begin{eqnarray}
\label{variances}
	\label{S11}
	\langle x^2\rangle &=&\frac{[(\alpha-\alpha^2+\epsilon^{\prime 2})(1-\alpha-2\epsilon^\prime\kappa)-2\alpha\epsilon^{\prime 2}]M+2[\epsilon^{\prime 2}+(1+\kappa^2)\alpha^2+\alpha]}{S_0k(\epsilon^{\prime 2}+\alpha)}T_c \nonumber \\ &+& \frac{\epsilon^{\prime 2}(1+\alpha+2\epsilon^\prime\kappa)M+2(1+\kappa^2)\epsilon^{\prime 2}+2\alpha\kappa^2}{S_0k(\epsilon^{\prime 2}+\alpha)}T_h, \\
	\label{S22}
	\langle y^2\rangle &=&\frac{\epsilon^{\prime 2}(1+\alpha+2\epsilon^\prime\kappa)M+(1+\kappa^2)\epsilon^{\prime 2}+\alpha\kappa^2}{S_0k(\epsilon^{\prime 2}+\alpha)}T_c \nonumber \\ &+&\frac{[(1-\alpha^2)-\epsilon^\prime(1-\alpha)(\epsilon^\prime-2\epsilon^\prime\kappa)-2\epsilon^{\prime 2}(1+\epsilon^\prime\kappa)]M+2(1+\alpha+\epsilon^{\prime 2}+\kappa^2)}{S_0k(\epsilon^{\prime 2}+\alpha)}T_h, \\
	\label{S33}
	 \langle v_x^2\rangle &=&k\frac{[(\alpha-1)^2-2\epsilon^\prime(\epsilon^\prime+(1-\alpha)\kappa)]M+(2+\kappa^2)(1+\alpha-2\epsilon^\prime\kappa)}{\gamma^2 M S_0}T_c+k\frac{\kappa^2(1+\alpha-2\epsilon^\prime\kappa)+2\epsilon^{\prime 2}M}{\gamma^2 M S_0}T_h,
 \\
	\label{S44}
	 \langle v_y^2\rangle &=&k\frac{\kappa^2(1+\alpha-2\epsilon^\prime\kappa)+2\epsilon^{\prime 2}M}{\gamma^2 M S_0}T_c +k\frac{[(\alpha-1)^2-2\epsilon^\prime(\epsilon^\prime-(1-\alpha)\kappa)]M+(2+\kappa^2)(1+\alpha-2\epsilon^\prime\kappa)}{\gamma^2 M S_0}T_h, \\ 
\end{eqnarray} 

where $M=km/\gamma^2$ is a dimensionless parameter and the non-diagonal elements are given as 
\begin{eqnarray}
\label{variances}
	 \label{S12}
	 \langle xy\rangle &=&\frac{[\alpha(\alpha-1+2\epsilon^\prime\kappa)-2\epsilon^{\prime 2}]\epsilon^\prime M+ 2[\alpha(\epsilon^\prime-\kappa+\epsilon^\prime\kappa^2)-\epsilon^{\prime 2}]}{S_0k(\epsilon^{\prime 2}+\alpha)}T_c \nonumber \\&+& \frac{(2\epsilon^{\prime 2}+\alpha-1-2\epsilon^\prime\kappa)\epsilon^\prime M+2[\kappa(\alpha+\epsilon^{\prime 2})-\epsilon^\prime(1+\kappa^2)]}{S_0k(\epsilon^{\prime 2}+\alpha)}T_h, \\
	\label{S21}
	\langle yx\rangle &=& \langle xy\rangle ,\\
	\label{S13S31}
	\langle xv_x\rangle &=& \langle v_xx\rangle= 0,\\
	\label{S24S42}
	\langle yv_y\rangle &=&\langle v_yy\rangle= 0,\\
    \label{S14}
	\langle xv_y\rangle &=&-\frac{\kappa(1-\alpha)(T_h-T_c)+2\epsilon^\prime(1+\kappa^2)(T_c+T_h)}{ S_0},\\
	\label{S32}
	\langle v_yx\rangle &=&\langle xv_y\rangle,\\
 	\label{S14S41}
	\langle yv_x\rangle &=&-\langle xv_y\rangle,\\
	\label{S14S41}
	\langle v_x y\rangle &=&\langle yv_x\rangle,\\
	\label{S34}
	 \langle v_xv_y\rangle &=&k\frac{[(\alpha-1)(T_c+T_h)-2\kappa\epsilon^{\prime 2}(T_h-T_c)]\epsilon^\prime M+\kappa(1+\alpha-2\epsilon^\prime\kappa)(T_h-T_c)}{\gamma M S_0}, \\
	\label{S34S43}
	 \langle v_yv_x\rangle &=& \langle v_xv_y\rangle,\\
\end{eqnarray}
where $S_0=[(\alpha-1)^2-4\epsilon^{\prime 2}(1+\kappa^2)]M+2(1+\kappa^2)(1+\alpha-2\epsilon^\prime\kappa)$.  Note that taking the small-mass limit of the covariance matrix gives the correlations in the overdamped regime. Since we finally need the kinetic energy of the particle which is proportional to $M(\langle v_x^2\rangle + \langle v_y^2\rangle)$, in the overdamped regime the small-mass limit of $M\langle v_x^2\rangle$ and $M\langle v_y^2\rangle$ is required which makes Eq.~\eqref{S33} and Eq.~\eqref{S44} well-defined. However, for the cross-correlation in Eq.~\eqref{S34} the explicit knowledge of the particle's mass is needed even in the overdamped limit.


\begin{figure}[t]
\centering
\includegraphics[width=1\linewidth]{figure05}
\caption{The stationary-state probability density of the particle's position and fluxes in a system with $M=1.0$, $\alpha=4.0$, $\kappa=1.0$, $u=0.0$, and $T_h=4T_c=4.0$. The results are shown from Brownian dynamics simulations to confirm the stability of the microengine for $\epsilon = 0.0$, $\epsilon<\epsilon_s$, and the stalled engine $\epsilon=\epsilon_s$ from left to right, respectively.  }
\label{fig:stability}
\end{figure}


The efficiency is given by the ratio between the average rate of heat extracted from the hot bath and the total energy transferred. The average rate of the heat out of the hot bath,  $Q_h=\gamma\kappa\langle v_xv_y\rangle-\epsilon\langle xv_y-yv_x\rangle/2$, can be calculated using Eq. \eqref{S34} which reads as $Q_h=f(\kappa, M, \epsilon)/M$ where $Q_h = \langle\,\dbar Q_h/dt\rangle$ and
\begin{equation}
\label{heatloss}
f(\kappa, M, \epsilon)=k\frac{[(\kappa(\alpha-1)+2\epsilon^\prime(1+\kappa^2))(T_c+T_h)+\kappa(1-\alpha-2\kappa\epsilon^{\prime 2})(T_h-T_c)]\epsilon^\prime M+\kappa^2(1+\alpha-2\epsilon^\prime\kappa)(T_h-T_c)}{\gamma S_0},
\end{equation}
and similarly the average mechanical power, $P=-\epsilon\langle xv_y-yv_x\rangle$, can be written as
\begin{equation}
\label{mechanicalpower}
P = 2k\epsilon^\prime\frac{\kappa(1-\alpha)(T_h-T_c)+2\epsilon^\prime(1+\kappa^2)(T_c+T_h)}{\gamma S_0}.
\end{equation}
The efficiency of the magneto-gyrator, defined as $\eta=P/Q_h$, can be calculated and written as
\begin{equation}
\label{efficiency}
\eta = \frac{2\epsilon^\prime M[\kappa(1-\alpha)\eta_c+2\epsilon^\prime(1+\kappa^2)(2-\eta_c)]}{\epsilon^\prime M[(\kappa(\alpha-1)+2\epsilon^\prime(1+\kappa^2))(2-\eta_c)+\kappa(1-\alpha-2\kappa\epsilon^{\prime 2})\eta_c] +\kappa^2(1+\alpha-2\epsilon^\prime\kappa)\eta_c},
\end{equation}
where $\eta_c=1-T_c/T_h$ is the Carnot efficiency. The efficiency as a function of the diffusive Hall parameter is plotted in Fig. $3$ (d) in the main text.


The existence condition of the steady state can be found by investigating the stability condition, which for our proposed magneto-gyrator can be done via studying the divergence of mechanical power in Eq.\eqref{mechanicalpower}. This leads to $S_0>0$, which can be written as
\begin{equation}
\label{stabilitycondition}
[(\alpha-1)^2-4\epsilon^{\prime 2}(1+\kappa^2)]M+2(1+\kappa^2)(1+\alpha-2\epsilon^\prime\kappa) > 0,
\end{equation}
which in the case of $\alpha=1.0$ reduces to $1+\epsilon^\prime\kappa-M\epsilon^{\prime 2}>0$ which has been reported in Ref.~\cite{lee2019nonequilibrium}. In Fig.~\ref{fig:stability} we show the stationary-state probability density of the particle's position and fluxes in a system with the suggested parameters for the possible experimental realisation in the main text. The simulation results are in agreement with our theoretical stability condition and confirm that the magneto-gyrator is stable for all chosen parameters in this work.   


\section*{COP of a Brownian Magneto-Refrigerator}
A Brownian magneto-gyrator can be operated as a refrigerator: by applying the nonconservative force $\F_{nc} = \epsilon (-y, x)^\top$ to the system with $\epsilon<0$. The sign of the parameter $\epsilon$ is chosen such that the resulting torque in the $z$ direction is in favour of the existed torque by the Lorentz force, namely we enhance the gyration along the direction of the unloaded magneto-refrigerator. Moreover, we operate the refrigerator in a range of $\epsilon$ where the system is stable. By doing so, we extract heat from the cold bath and pump it into the hot bath. We obtain the coefficient of performance (COP) of refrigerator by dividing the average rate of the heat extracted from the hot bath by the corresponding input work $Q_h/P$, given in Eq.~\eqref{heatloss} and Eq.~\eqref{mechanicalpower}, respectively. Therefore, the COP of the magneto-refrigerator can be written as
\begin{equation}
\label{COP}
\text{COP}=\frac{\epsilon^\prime M[(\kappa(\alpha-1)+2\epsilon^\prime(1+\kappa^2))(T_c+T_h)+\kappa(1-\alpha-2\kappa\epsilon^{\prime 2})(T_h-T_c)] +\kappa^2(1+\alpha-2\epsilon^\prime\kappa)(T_h-T_c)}{2\epsilon^\prime M[\kappa(1-\alpha)(T_h-T_c)+2\epsilon^\prime(1+\kappa^2)(T_c+T_h)]}.
\end{equation}   

The COP as a function of the parameter $\epsilon$ is plotted in Fig.~$4$ in the main text. The COP has a maximum at an optimal $\epsilon$ which is presented in Fig.~$5$ in the main text, as well.

\providecommand{\noopsort}[1]{}\providecommand{\singleletter}[1]{#1}%
\begin{thebibliography}{1}
\urlstyle{rm}
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\expandafter\ifx\csname doiprefix\endcsname\relax\def\doiprefix{DOI: }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem{vuijk2019anomalous}
\bibinfo{author}{Vuijk, H.~D.}, \bibinfo{author}{Brader, J.~M.} \&
  \bibinfo{author}{Sharma, A.}
\newblock \bibinfo{journal}{\bibinfo{title}{Anomalous fluxes in overdamped
  Brownian dynamics with Lorentz force}}.
\newblock {\emph{\JournalTitle{Journal of Statistical Mechanics: Theory and
  Experiment}}} \textbf{\bibinfo{volume}{2019}}, \bibinfo{pages}{063203}
  (\bibinfo{year}{2019}).

\bibitem{chun2018emergence}
\bibinfo{author}{Chun, H.-M.}, \bibinfo{author}{Durang, X.} \&
  \bibinfo{author}{Noh, J.~D.}
\newblock \bibinfo{journal}{\bibinfo{title}{Emergence of nonwhite noise in
  Langevin dynamics with magnetic Lorentz force}}.
\newblock {\emph{\JournalTitle{Physical Review E}}}
  \textbf{\bibinfo{volume}{97}}, \bibinfo{pages}{032117}
  (\bibinfo{year}{2018}).

\bibitem{abdoli2020correlations}
\bibinfo{author}{Abdoli, I.} \emph{et~al.}
\newblock \bibinfo{journal}{\bibinfo{title}{Correlations in multithermostat
  Brownian systems with Lorentz force}}.
\newblock {\emph{\JournalTitle{New Journal of Physics}}}
  \textbf{\bibinfo{volume}{22}}, \bibinfo{pages}{093057}
  (\bibinfo{year}{2020}).

\bibitem{abdoli2021stochastic}
\bibinfo{author}{Abdoli, I.} \& \bibinfo{author}{Sharma, A.}
\newblock \bibinfo{journal}{\bibinfo{title}{Stochastic resetting of active
  Brownian particles with Lorentz force}}.
\newblock {\emph{\JournalTitle{Soft Matter}}} \textbf{\bibinfo{volume}{17}},
  \bibinfo{pages}{1307--1316} (\bibinfo{year}{2021}).
  
\bibitem{vuijk2019effect}
\bibinfo{author}{Vuijk, H.~D.}, \bibinfo{author}{Brader, J.~M.} \&
  \bibinfo{author}{Sharma, A.}
\newblock \bibinfo{journal}{\bibinfo{title}{Effect of anisotropic diffusion on
  spinodal decomposition}}.
\newblock {\emph{\JournalTitle{Soft Matter}}} \textbf{\bibinfo{volume}{15}},
  \bibinfo{pages}{1319--1326} (\bibinfo{year}{2019}).
  
\bibitem{abdoli2020nondiffusive}
\bibinfo{author}{Abdoli, I.}, \bibinfo{author}{Vuijk, H.~D.},
  \bibinfo{author}{Sommer, J.-U.}, \bibinfo{author}{Brader, J.~M.} \&
  \bibinfo{author}{Sharma, A.}
\newblock \bibinfo{journal}{\bibinfo{title}{Nondiffusive fluxes in a Brownian
  system with Lorentz force}}.
\newblock {\emph{\JournalTitle{Physical Review E}}}
  \textbf{\bibinfo{volume}{101}}, \bibinfo{pages}{012120}
  (\bibinfo{year}{2020}).
  
\bibitem{abdoli2020stationary}
\bibinfo{author}{Abdoli, I.} \emph{et~al.}
\newblock \bibinfo{journal}{\bibinfo{title}{Stationary state in Brownian
  systems with Lorentz force}}.
\newblock {\emph{\JournalTitle{Physical Review Research}}}
  \textbf{\bibinfo{volume}{2}}, \bibinfo{pages}{023381} (\bibinfo{year}{2020}).

\bibitem{liebchen2019optimal}
\bibinfo{author}{Liebchen, B.} \& \bibinfo{author}{L{\"o}wen, H.}
\newblock \bibinfo{journal}{\bibinfo{title}{Optimal navigation strategies for
  active particles}}.
\newblock {\emph{\JournalTitle{EPL (Europhysics Letters)}}}
  \textbf{\bibinfo{volume}{127}}, \bibinfo{pages}{34003}
  (\bibinfo{year}{2019}).
  
\bibitem{mangeat2019role}
\bibinfo{author}{Mangeat, M.}, \bibinfo{author}{Amarouchene, Y.},
  \bibinfo{author}{Louyer, Y.}, \bibinfo{author}{Gu{\'e}rin, T.} \&
  \bibinfo{author}{Dean, D.~S.}
\newblock \bibinfo{journal}{\bibinfo{title}{Role of nonconservative scattering
  forces and damping on Brownian particles in optical traps}}.
\newblock {\emph{\JournalTitle{Physical Review E}}}
  \textbf{\bibinfo{volume}{99}}, \bibinfo{pages}{052107}
  (\bibinfo{year}{2019}).

\bibitem{seifert2012stochastic}
\bibinfo{author}{Seifert, U.}
\newblock \bibinfo{journal}{\bibinfo{title}{Stochastic thermodynamics,
  fluctuation theorems and molecular machines}}.
\newblock {\emph{\JournalTitle{Reports on Progress in Physics}}}
  \textbf{\bibinfo{volume}{75}}, \bibinfo{pages}{126001}
  (\bibinfo{year}{2012}).

\bibitem{lee2019nonequilibrium}
\bibinfo{author}{Lee, S.} \& \bibinfo{author}{Kwon, C.}
\newblock \bibinfo{journal}{\bibinfo{title}{Nonequilibrium driven by an
  external torque in the presence of a magnetic field}}.
\newblock {\emph{\JournalTitle{Physical Review E}}}
  \textbf{\bibinfo{volume}{99}}, \bibinfo{pages}{052142}
  (\bibinfo{year}{2019}).
  

  
\end{thebibliography}



\end{document}