\section{Introduction}
There are many settings in random matrix theory for which the eigenvalues (assumed real) can be scaled in relation to the
matrix size in such a way that the limiting support is compact. This is referred to as a global scaling. As some concrete examples,
let $X$ be an $N \times N$ standard complex Gaussian matrix, and construct from this the Hermitian matrices
$H_1 = {1 \over 2} (X + X^\dagger)$ and $H_2 = X^\dagger X$. The set of matrices $H_1$ ($H_2$) are said to form the
Gaussian unitary ensemble (special case of the Laguerre unitary ensemble), and have joint eigenvalue probability
density function (PDF) proportional to
\begin{equation}\label{1.1}
\prod_{l=1}^N w(x_l) \prod_{1 \le j < k \le N} (x_k - x_j)^2, \qquad w(x) = \left \{ \begin{array}{cc} e^{-x^2}, & {\rm matrices} \: H_1 \\
 e^{-x} \chi_{x > 0}, & {\rm matrices} \: H_2; \end{array} \right.
 \end{equation}
see e.g.~\cite{Fo10,PS11}. Here $\chi_A = 1$ for $A$ true, $\chi_A = 0$ otherwise.

Scaling the eigenvalues $x_j \mapsto \sqrt{2N} x_j $ (matrices $H_1$) and $x_j \mapsto 4 N x_j $ (matrices $H_2$), it is a standard
result that as $N \to \infty$ the spectrum is supported on the intervals $(-1,1)$ and $(0,1)$ respectively. Among the endpoints of
the intervals of support, the point $x = 0$ for the global scaling of the matrices $H_2$ is special. Thus the region $x < 0$
to the other side of this endpoint has strictly zero eigenvalue density for all values of $N$, because $H_2$ is positive
definite. For this reason the endpoint $x = 0$ in this example is called a hard edge. The hard edge notion extends beyond the class of
matrix ensembles permitting a global scaling to include heavy tailed distributions --- an example of the latter is given in
Section \ref{S3.3} below. The essential point then is that the limiting eigenvalue density is nonzero for $x>0$, and strictly
zero for $x<0$.

In this paper our interest is in the approach to a limiting hard edge state for various ensembles of positive definite matrices.
A hard edge state refers to the statistical distribution formed when the eigenvalues are scaled to have nearest neighbour
spacing of order unity as $N \to \infty$. For the matrices $H_2$, or more generally the ensemble of matrices with
weight function
\begin{equation}\label{Lw}
w(x) = x^a e^{-x} \chi_{x > 0}
 \end{equation}
 (Laguerre weight, realised for $a = n - N \in \mathbb Z_{\ge 0}$ as the eigenvalue PDF of
 matrices $X^\dagger X$ with $X$ an $n \times N$ complex standard Gaussian matrix) with parameter $a > - 1$,  this takes place for the scaling of the eigenvalues
$x_j \mapsto x_j / 4N$, and gives rise to the hard edge state specified by the $k$-point correlations (see \cite[\S 7.2]{Fo10})
\begin{equation}\label{1.2}
\rho_{(k)}^{\rm hard}(x_1,\dots,x_k) = \det [ K^{\rm hard}(x_j, x_l;a) ]_{j,l=1}^k,
 \end{equation}
 where, with $J_a(u)$ denoting the Bessel function,
\begin{equation}\label{1.3}
K^{\rm hard}(x,y;a)
 = {1 \over 4} \int_0^1 J_a(\sqrt{xt})  J_a(\sqrt{yt}) \, dt.
\end{equation}

For finite $N$ the $k$-point correlation function is defined in terms of the joint eigenvalue PDF, $P_N$ say,
according to
\begin{equation}\label{1.4}
\rho_{(k)}(x_1,\dots,x_k) = {N! \over (N - k)!} \int_{-\infty}^\infty dx_{k+1} \cdots  \int_{-\infty}^\infty dx_N \, P_N(x_1,\dots, x_N).
 \end{equation}
 For eigenvalue PDFs of the form (\ref{1.1}), the correlation function (\ref{1.4}) admits the determinant evaluation (see e.g.~\cite[\S 5.1]{Fo10})
\begin{equation}\label{1.5a}
\rho_{(k)}(x_1,\dots,x_k) =  \det [ K_N(x_j, x_l ) ]_{j,l=1}^k,
 \end{equation}
 where
 \begin{align}\label{1.5b}
  K_N(x,y)  = \Big ( w(x) w(y) \Big )^{1/2} \sum_{n=0}^{N-1} {1 \over h_n} p_n(x) p_n(y) \nonumber \\
 \end{align}
 In (\ref{1.5b}) $\{p_n(x)\}$ refers to the set of orthogonal polynomials with respect to the weight function $w(x)$ ---
 $p_n$ of degree $n$ and chosen to be monic for convenience --- with norm $h_n$,
 \begin{equation}\label{1.5c}
 \int_{-\infty}^\infty w(x) p_m(x) p_n(x) \, dx = h_n \delta_{m,n}.
  \end{equation}
  In the case of the Laguerre weight, the polynomials $p_n(x)$ are proportional to the Laguerre
  polynomials $L_n^{(a)}(x)$.

  Recently, attention has been given to  the rate of convergence to the hard edge limiting kernel
  (\ref{1.3}). One line of motivation came from a question posed by
 Edelman, Guionnet and P\'ech\'e \cite{EGP16}. These authors, taking a viewpoint in numerical analysis, took up the problem of studying
 finite $N$ effects in the hard edge scaling of the distribution of the smallest singular value of a (complex) standard Gaussian matrix.
 With $E^{\rm LUE}(0;(0,s))$ denoting the probability that there are no eigenvalues in the interval $(0,s)$ of the LUE, it was conjectured
 in \cite{EGP16} that
 \begin{equation}\label{pk5}
E^{\rm LUE}(0;(0,s/(4N))) =  E^{\rm hard}(0;(0,s)) + {a \over 2N} s {d \over ds} E^{\rm hard}(0;(0,s)) + O \Big ( {1 \over N^2} \Big ),
 \end{equation}
 where
 $$
 E^{\rm hard}(0;(0,s)) = \lim_{N \to \infty}  E^{\rm LUE}(0;(0,s/(4N))),
 $$
 and thus  \cite{Bo16,PS16},
  \begin{equation}\label{pk7}
 E^{\rm LUE}\bigg ( 0; \Big ( 0, {s \over 4 N + 2a} \Big )  \bigg ) =  E^{\rm hard}(0;(0,s)) + O \Big ( {1 \over N^2} \Big ),
 \end{equation}
which moreover is the optimal rate of convergence.

 Subsequently Bornemann \cite{Bo16}
 provided a proof of (\ref{pk5}) which involved extending the limit formula (\ref{1.2}) to the large $N$ expansion
 \begin{align}\label{pk6}
 {1 \over 4 N} K_N^{(L)} \Big ( {X \over 4 N} , {Y \over 4 N} \Big ) & =  K^{\rm hard}(X,Y)  + {1 \over N} {a \over 8} J_a(\sqrt{X})  J_a(\sqrt{Y})  + O \Big ( {1 \over N^2} \Big ) \nonumber \\
 &  =  K^{\rm hard}(X,Y)  + {1 \over N}  {a \over 2} \bigg ( x {\partial \over \partial x} +   y {\partial \over \partial y}  + 1 \bigg )  K^{\rm hard}(X,Y) + O \Big ( {1 \over N^2} \Big ),
 \end{align}
 valid uniformly for $X,Y \in [0,s]$. In fact knowledge of (\ref{pk6}) is sufficient to establish (\ref{pk5}).
 We remark too that analogous to (\ref{pk7}), it follows from  (\ref{pk6}) that
\begin{equation}\label{pk8}
 {1 \over 4 N + 2a } K_N^{(L)} \Big ( {X \over 4 N + 2a} , {Y \over 4 N + 2a} \Big )  =  K^{\rm hard}(X,Y)  + O \Big ( {1 \over N^2} \Big ),
 \end{equation}
 and this implies (\ref{pk7}).

 Our aim in this work is to extend hard edge scaling results of the type (\ref{pk6}) to examples of a recently isolated structured class of random
 matrices known as P\'olya ensembles \cite{KK16}. The definition of these ensembles, which include the Laguerre unitary ensemble, the Jacobi unitary ensemble,
 products of these ensembles, and their Muttalib-Borodin generalisations, will be given in Section \ref{S2.1}. The benefit of the structures provided
 by the P\'olya ensemble class is seen by our revision of the key formulas in Section \ref{S2.2}, where we also extend the theory by exhibiting differential
 recurrences satisfied by the associated biothogonal pair, and a differential identity satisfied by the correlation kernel. In Section \ref{S2.3} we make note of some
 asymptotic formulas relating to ratios of gamma functions which will be used in our subsequent large $N$ hard edge analysis. The latter is undertaken is Section \ref{S3},
 starting with products of Laguerre ensembles, then the Laguerre Muttalib-Borodin ensemble, and finally products of Laguerre ensembles and their inverses,
 with the latter including as a special case the Jacobi unitary ensemble.

 The Jacobi unitary ensemble is specified by the eigenvalue PDF (\ref{1.1}) with weight
 \begin{equation}\label{A.2}
 x^a (1 - x)^b \chi_{0 < x < 1}.
 \end{equation}
 Our results of Section \ref{S3.3} imply that
  \begin{align}\label{pk6a}
 {1 \over 4 N^2} K_N^{(J)} \Big ( {X \over 4 N^2} , {Y \over 4 N^2} \Big ) & =  K^{\rm hard}(X,Y)  +  { a+b \over 2 N} J_a(\sqrt{X})  J_a(\sqrt{Y})  + O \Big ( {1 \over N^2} \Big ) \nonumber \\
 &  =  K^{\rm hard}(X,Y)  +   { a +b \over N}  \bigg ( x {\partial \over \partial x} +   y {\partial \over \partial y}  + 1 \bigg )  K^{\rm hard}(X,Y) + O \Big ( {1 \over N^2} \Big ),
 \end{align}
and thus
 \begin{equation}\label{A.2a}
  {1 \over 4 \tilde{N}^2} K_N^{(J)} \Big ( {X \over 4 \tilde{N}^2} , {Y \over 4 \tilde{N}^2} \Big )  \bigg |_{\tilde{N} = N + (a+b)/2} =
   K^{\rm hard}(X,Y)  +   O \Big ( {1 \over N^2} \Big ).
   \end{equation}
   This gives an explanation for recent results in \cite{MMM19} relating to the large $N$ form of the distribution
   of the smallest eigenvalue in the Jacobi unitary ensemble. In Appendix A large $N$ expansions of the latter
   quantity are extended to all Jacobi $\beta$-ensembles with $\beta$ even.

  \section{Preliminaries}
  \subsection{P\'olya ensembles --- definitions}\label{S2.1}
  The Vandermonde determinant identity tells us that
 \begin{equation}\label{2.1}
 \det [ x_k^{j-1} ]_{j,k=1}^N = \det [p_{j-1}(x_k) ]_{j,k=1}^{N} = \prod_{1 \le j < k \le N} (x_k - x_j),
   \end{equation}
   where $\{p_{l}(x)\}_{l=0}^{N-1}$ are arbitrary monic orthogonal polynomials, $p_l$ of degree $l$.
   A generalisation of (\ref{1.1}) is therefore an eigenvalue PDF proportional to
  \begin{equation}\label{2.2}
   \det [p_{j-1}(x_k) ]_{j,k=1}^{N} \det [ w_{j-1}(x_k) ]_{j,k=1}^N
 \end{equation}
 for some polynomials $\{p_{l}(x)\}_{l=0}^{N-1}$ and functions
 $\{ w_j(x) \}_{j=0}^{N-1}$ --- note though that in general there is no guarantee
   (\ref{2.2}) will be positive. In \cite{KZ14} eigenvalue PDFs (\ref{2.2}) were given the name
 polynomial ensembles.

 In \cite{KK16,KK19} a further specialisation of (\ref{2.2}),
  \begin{equation}\label{2.3}
\det [p_{j-1}(x_k) ]_{j,k=1}^{N}   \det \bigg [ \Big ( - x_k {\partial \over \partial x_k} \Big )^{j-1} w(x_k) \bigg ]_{j,k=1}^N,
 \end{equation}
 was proposed. Assuming all the eigenvalues are positive, it was  shown that this class of eigenvalue
 PDF is closed under multiplicative convolution. At first  PDFs of the form (\ref{2.3}) were referred to as
 polynomial ensembles of derivative type, but subsequently with the requirement that they be non-negative,
 it was pointed out in \cite{FKK17} that it is more apt to use the term P\'olya ensemble.  The invariance of a
 determinant under the elementary row operation of adding one multiple of a row to another shows
   \begin{align}\label{2.4}
 \det \bigg [ \Big ( - x_k {\partial \over \partial x_k} \Big )^{j-1} w(x_k) \bigg ]_{j,k=1}^N
 & =     \det \bigg [  \prod_{l=1}^{j-1} \Big ( - x_k {\partial \over \partial x_k} - l \Big ) w(x_k)  \bigg ]_{j,k=1}^N \nonumber \\
 & =  \det \bigg [ {\partial^{j-1} \over \partial x_k^{j-1}} \Big ( (-x_k)^{j-1} w(x_k) \Big ) \bigg ]_{j,k=1}^N.
 \end{align}
 In relation to the second line, note that it is in fact an equality that
   \begin{equation}\label{2.5}
    \prod_{l=1}^{j-1} \Big ( - x {\partial \over \partial x} - l \Big ) w(x) = {d^{j-1} \over d x^{j-1}} \Big ( (-x)^{j-1} w(x) \Big ).
 \end{equation}

 The differential operator on the RHS of (\ref{2.5}) reveals that the Laguerre unitary ensemble fits the framework of
P\'olya ensembles. Thus choosing $w(x)$ to be given by (\ref{Lw}), the Rodrigues formula for the Laguerre polynomials
tells us that
  \begin{equation}\label{2.6}
  {d^{j-1} \over d x^{j-1}} \Big ( (-x)^{j-1} w(x) \Big ) = (-1)^{j-1} (j-1)! w(x) L_{j-1}^{(a)}(x),
  \end{equation}
  and so, up to proportionality, (\ref{2.3}) reduces to
   \begin{equation}\label{2.7}
   \prod_{l=1}^N x_l^a e^{-x_l} \det [ p_{j-1}(x_k) ]_{j,k=1}^N    \det [ L_{j-1}^{(a)}(x_k) ]_{j,k=1}^N.
   \end{equation}
  In view of (\ref{2.1}), this corresponds to the eigenvalue PDF for the Laguerre unitary ensemble. The advantage in  working
  within the P\'olya ensemble framework is that it reveals a mechanism to obtain the  asymptotic expansion of the
  correlation kernel  (\ref{1.5b}) at the hard edge, which applies at once to a much wider class of random matrix ensembles.
  The reason for this are certain general structural formulas applicable to all P\'olya ensembles. These will be revised next.

  \subsection{P\'olya ensembles --- biorthogonal system and correlation kernel}\label{S2.2}
  It is standard in random matrix theory that the ensembles (\ref{2.2}) are determinantal, meaning that the $k$-point correlation
  functions have the form (\ref{1.5a}). Moreover, if the polynomials $\{p_l(x) \}_{l=0}^N$ and the functions $\{q_j(x)\}_{j=0}^N$ --- the latter
  chosen from ${\rm span} \, \{w_j(x) \}_{j=0}^{N}$ --- have the
  biorthogonal property
     \begin{equation}\label{2.8}
 \int_{-\infty}^\infty p_m(x) q_n(x) \, dx = \delta_{m,n},
 \end{equation}
 then the correlation kernel has the simple form
    \begin{equation}\label{2.9}
    K_N(x,y) = \sum_{j=0}^{N-1} p_{j}(x) q_{j}(y);
 \end{equation}
 see e.g.~\cite[\S 5.8]{Fo10}.
 While in general computation of  the LU (lower/ upper triangular) decomposition of a certain inverse matrix
  used to construct the biorthogonal
 functions (see e.g.~\cite[Proof of Prop.~5.8.1]{Fo10}), this cannot be expected to result in
a tractable    formula for (\ref{2.9}), permitting large $N$ analysis, without further structures.
It is at this stage that the utility of P\'olya ensembles shows itself: special functional forms for the
 biorthogonal system hold true, and moreover there is a summed up form of the kernel
 as an integral analogous to (\ref{1.3}), which together facilitate a large $N$ analysis.

 The formulas, which are due to Kieburg and K\"osters \cite{KK16}, involve the Mellin transform
 of the weight $w$ in (\ref{2.3}),
   \begin{equation}\label{Me}
   \mathcal M[w](s) := \int_0^\infty y^{s-1} w(y) \, dy.
  \end{equation}
One has that the polynomials $\{p_l(x)\}_{l=0}^N$ in the biorthogonal pair $\{p_j, q_k\}$ are specified by
   \begin{equation}\label{Me1}
 p_n(x) = (-1)^n  n!   \mathcal M[w](n+1) \sum_{j=0}^n{  (-x)^{j}  \over j! (n-j)!   \mathcal M[w](j+1) } ,
 \end{equation}
and that the functions $\{q_l(x)\}_{l=0}^N$  --- chosen from the span of the functions specifying the columns in (\ref{2.3}) ---
are specified by the Rodrigues type formula
  \begin{equation}\label{Me2}
  q_n(x) = {1 \over n! \mathcal M[w](n+1)} {d^n \over d x^n} \Big((-x)^n w(x)\Big).
 \end{equation}
 Moreover, the correlation kernel can be written in a form generalising the final expression in (\ref{1.3}),
  \begin{equation}\label{Me3}
  K_N(x,y) = - N { \mathcal M[w](N+1) \over  \mathcal M[w](N)}  \int_0^1 p_{N-1}(xt) q_N(yt) \, dt.
   \end{equation}

  In \cite{KK16} the integral form (\ref{Me3}) of the correlation kernel was derived by first converting
  (\ref{Me1}) and (\ref{Me2}) to integral forms, which allow for the summation to be carried out
  in closed form. The identification with the RHS of (\ref{Me3}) then follows after some manipulation.
  In a special case this strategy was first given in \cite{KZ14}.
  An alternative method of derivation is also possible, as we will now show, which involves first
  identifying differential recurrences satisfied by each of the $p_n(x)$ and $q_n(x)$.
  (We remark that other examples of differential recurrences can be found in a number of recent studies in random matrix
  theory \cite{Ku19,FK19,FK20,Fo20a,Fo20b}.)

  \begin{proposition}\label{PR1}
  Let $p_n(x)$ and $q_n(x)$ be specified by (\ref{Me1}) and (\ref{Me2}).
  These functions satisfy the differential recurrences
  \begin{align}
  x {d \over dx} p_n(x) & =  n p_n(x) + n {\mathcal M [w](n+1) \over \mathcal M [w] (n) } p_{n-1}(x) \label{PQ1} \\
  x {d \over dx} q_n(x) & =  - {(n+1) \mathcal M[w](n+2) \over  \mathcal M[w](n+1)} q_{n+1}(x) + (n+1) q_n(x).\label{PQ2}
  \end{align}
  A corollary of these recurrences is the differential identity
   \begin{equation}\label{Me4}
  \Big ( x{\partial \over \partial x} + y{\partial \over \partial y}  + 1 \Big )   K_N(x,y) = -  N { \mathcal M[w](N+1)  \over  \mathcal M[w](N)}  p_{N-1}(x) q_N(y),
    \end{equation}
    which implies (\ref{Me3}).
  \end{proposition}

  \begin{proof}
  From the formula (\ref{Me1}),
  $$
  x {d \over dx}  p_n(x) = (-1)^n  n!   \mathcal M[w](n+1) \sum_{j=0}^n (-1)^{j} {j \over j! (n-j)!   \mathcal M[w](j+1) } x^j.
 $$
 Rewrite the $j$ in the denominator of this expression as $n - (n - j)$, and use this to decompose the sum into two.
 Upon some simple manipulation, the identity  (\ref{PQ1}) results.

 According to (\ref{2.5}), the  formula (\ref{Me2}) can be rewritten
 $$
  q_n(x) = {1 \over n! \mathcal M[w](n+1)}     \prod_{l=1}^{n} \Big ( - x {\partial \over \partial x} - l \Big )  w(x).
$$
 Acting on both sides with $- x {d \over dx} - (n+1)$ shows
 $$
 \Big (  - x {d \over dx} - (n+1) \Big ) q_n(x) = {(n+1) \mathcal M[w](n+2) \over  \mathcal M[w](n+1)} q_{n+1}(x).
 $$
This gives (\ref{PQ2}).

With the differential recurrences (\ref{PQ1}) and (\ref{PQ2}) established, we can use them in the expression (\ref{2.9}) to give
\begin{multline}
\Big ( x{\partial \over \partial x} + y{\partial \over \partial y} \Big ) K_N(x,y)  \\
= \sum_{n=0}^{N-1} \Big ( n p_n(x) + n {\mathcal M [w](n+1) \over \mathcal M [w] (n) } p_{n-1}(x) \Big )
\Big (  - {(n+1) \mathcal M[w](n+2) \over  \mathcal M[w](n+1)} q_{n+1}(y) + (n+1) q_n(y) \Big ).
\end{multline}
Simple manipulation reduces this to (\ref{Me4}).

In (\ref{Me4}) scale $x$ and $y$ by writing as $xt$ and $yt$ respectively. The LHS of (\ref{Me4}) can then
be written
\begin{equation}\label{2.17a}
{d \over dt} t K_N(tx,ty) = -  N { \mathcal M[w](N+1)  \over  \mathcal M[w](N)}  p_{N-1}(tx) q_N(ty).
 \end{equation}
Integrating both sides from $0$ to $1$, on the LHS noting $\lim_{t \to 0^+} t K_N(tx,ty) = 0$
as follows from (\ref{2.9}), reclaims  (\ref{Me3}).

\end{proof}

\begin{remark}\label{R2}
We show in Appendix B how (\ref{2.17a}), combined with a recurrence formula of fixed depth of
$t p_{N-1}(t)$ known to hold for a number of the specific P\'olya ensembles considered
in Section \ref{S3}, provides a combinatorial based method to compute the leading
large $N$ form of the moments of the spectral density.
\end{remark}

  \subsection{Asymptotics of ratios of gamma function}\label{S2.3}

  The gamma function $\Gamma(z)$ is one of the most commonly occurring of special functions \cite{andrews99}, analytic in the
complex plane except for poles at 0 and the negative integers.
Since $\Gamma(z+1) = z \Gamma(z)$ and $\Gamma(1) = 1$, for $n$ a non-negative integer
\begin{equation}\label{0.1}
 \Gamma(n+1) = n!,
 \end{equation}
 and so gives meaning to the factorial for
general complex $n$. Historically \cite{Pe24} Stirling's formula for the gamma function is the large $n$ approximation to the
factorial $n! \approx \sqrt{2 \pi } n^{n + 1/2} e^{-n}$, later extended to the asymptotic series \cite{WW65}
\begin{equation}
n! =  \sqrt{2 \pi n}  \Big ( {n \over e} \Big )^n \bigg ( 1 + {1 \over 12 n} + {1 \over 288 n^2} + O \Big ( {1 \over n^3} \Big ) \bigg ).
\end{equation}
Using (\ref{0.1}) and
truncating this asymptotic series at $O(1/n)$ leads to the large $|z|$ asymptotic expansion \cite{tricomi51}
\begin{align}\label{eq1}
\frac{\Gamma(z+a)}{\Gamma(z+b)}=z^{a-b}\left(
1+\frac{1}{2z}(a-b)(a+b-1)+O(z^{-2})
\right),\quad |z|\to\infty
\end{align}
valid for $|{\rm arg} \, z|  < \pi$ and $a,b$ fixed.  Furthermore,
specify $(u)_\alpha := \Gamma(u+\alpha)/\Gamma(u)$, which for $\alpha$ a positive integer corresponds to the
product $(u)_\alpha  = (u) (u+1) \cdots (u+ \alpha - 1)$. From this definition, and under the assumption that
$\alpha$ is a positive integer, we see
\begin{align}\label{eq2}
(-N+k)_{\alpha}=   (-1)^\alpha {\Gamma(N - k + 1) \over \Gamma(N - k + 1 - \alpha) } = (-N)^{\alpha}\left(
1-\frac{\alpha(2k+\alpha-1)}{2N}+O(N^{-2})
\right),\quad N\to\infty,
\end{align}
where the large $N$ form follows from (\ref{eq1}).
Our analysis of the rate of convergence for hard edge scalings will have use for both
 (\ref{eq1}) and (\ref{eq2}).

 \section{Hard edge scaling to $O(1/N)$ for some P\'olya ensembles}\label{S3}
 \subsection{Products of Laguerre ensembles}\label{S3.1}
The realisation of the Laguerre unitary ensemble with $a = n - N$ noted below (\ref{Lw}) can equivalently be expressed
as  being realised by the squared singular values of an $n \times N$ standard complex Gaussian matrix. A natural
generalisation, first considered in \cite{AKW13,AIK13}, is to consider the squared singular values of the product
of say $M$ rectangular standard complex Gaussian matrices (assumed to be of compatible sizes). Since each ensemble
in the product is individually a P\'olya ensemble, the closure property
of P\'olya ensembles under multiplicative convolution from \cite{KK16} tells us that the product ensemble can be formed by simply replacing $w(x)$
in (\ref{2.3}) by
\begin{equation}\label{W1}
w^{(M)}(x) := \int_0^\infty dx_1 \cdots dx_M \, \delta  \Big ( x - \prod_{j=1}^M x_j \Big ) \prod_{l=1}^M w_l(x_l), \quad w_j(x) = {1 \over \Gamma(a_j + 1)}  x^{a_j} e^{-x} .
\end{equation}
For the Mellin transform we have the factorised gamma function evaluation
\begin{equation}\label{W2}
\mathcal M[w^{(M)}](s) = \prod_{j=1}^M {\Gamma(a_j + s) \over  \Gamma(a_j + 1)}.
\end{equation}
The formula for the inverse Mellin transform then gives
\begin{align}\label{W3}
w^{(M)}(x)  & =  \Big ( \prod_{j=1}^M {1 \over  \Gamma(a_j + 1)} \Big )   {1 \over 2 \pi i }  \int_{c - i \infty}^{c + i \infty}  \prod_{j=1}^M \Gamma(a_j - s)  \, x^{s} \, ds \nonumber \\
& =   \prod_{j=1}^M {1 \over  \Gamma(a_j + 1)}  \,
\MeijerG{M}{0}{0}{M}{-}{a_1,\ldots,a_M}{x}.
\end{align}
 Here $c$ is any positive real number, and $G^{0,M}_{M,0}$ denotes a particular Meijer G-function; see \cite{MSH09}.

 Substituting (\ref{W2}) in (\ref{Me1}) and (\ref{W3}) in (\ref{Me2}) shows \cite{AIK13}
 \begin{align}\label{W4}
 p_n(x) & =   (-1)^n n! \prod_{j=1}^M \Gamma(a_j + n + 1) \sum_{j=0}^n { (-x)^j \over j! (n - j)! \prod_{l=1}^M (a_l +  1)_j}  \nonumber\\
& =
 (-1)^n \prod_{j=1}^M {\Gamma(a_j + n + 1) \over \Gamma(a_j + 1)} \,
 {}_1 F_M \bigg ( \begin{array}{cc} -n \\ a_1 + 1,\dots, a_M + 1 \end{array} \Big | x \bigg ),
 \end{align}
 with $ {}_1 F_M$ the notation for the particular hypergeometric series, and
\begin{align}\label{W5}
 q_n(x) & =   {(-1)^n \over n!}  \prod_{j=1}^M {1 \over \Gamma(a_j + n + 1)}  {1 \over 2 \pi i }
   \int_{c - i \infty}^{c + i \infty} {\Gamma(n+s+1) \over \Gamma(s+1) }  \prod_{j=1}^M \Gamma(a_j - s)  \, x^{s} \, ds \nonumber \\
& =  {(-1)^n \over n!}  \prod_{j=1}^M {1 \over  \Gamma(a_j +n+ 1)} \,
\MeijerG{M}{1}{1}{M+1}{-n}{a_1,\ldots,a_M,0}{x}.
 \end{align}

 According to (\ref{Me4}) and (\ref{Me3}), $K_N(x,y)$ is fully determined by $p_{N-1}(x)$ and $q_N(y)$. Since our aim is to
 expand $K_N(x,y)$ for large $N$ with hard edge scaled variables, it suffices then to compute the hard edge expansion
 of these particular biorthogonal functions.

 \begin{proposition}\label{P3.1}
 Denote
\begin{equation}\label{W6}
{}_0 F_M  \bigg ( \begin{array}{cc} -  \\ a_1 + 1,\dots, a_M + 1 \end{array} \Big | -x \bigg ) = \sum_{j=0}^\infty {(-x)^j \over j! \prod_{s=1}^M (a_s + 1)_j},
\end{equation}
as conforms with standard notation in the theory of hypergeometric functions. We have
\begin{multline}\label{W7}
 {}_1 F_M \bigg ( \begin{array}{cc} -N + 1 \\ a_1 + 1,\dots, a_M + 1 \end{array} \Big | {x \over N} \bigg ) \\
 = \bigg ( 1 - {1 \over 2 N} \Big ( x {d \over dx} +
 \Big (  x {d \over dx} \Big )^2 \Big ) \bigg ) \, {}_0 F_M  \bigg ( \begin{array}{cc} -  \\ a_1 + 1,\dots, a_M + 1 \end{array} \Big | -x \bigg ) + O \Big ( {1 \over N^2} \Big ).
 \end{multline}
 Also
 \begin{multline}\label{W8}
 {1 \over N!} \MeijerG{M}{1}{1}{M+1}{-N}{a_1,\ldots,a_M,0}{{x \over N}} \\
 = \bigg ( 1 + {1 \over 2 N} \Big ( x {d \over dx} +
 \Big (  x {d \over dx} \Big )^2 \Big ) \bigg )  \MeijerG{M}{0}{1}{M+1}{-}{a_1,\ldots,a_M,0}{{x }} + O \Big ( {1 \over N^2} \Big ).
 \end{multline}
 In both (\ref{W7}) and (\ref{W8}) the bound on the remainder holds uniformly for $x \in [0,s]$, for any fixed $s \in \mathbb R_+$.
 \end{proposition}

 \begin{proof}
 In the summation (\ref{W4}) defining  the LHS of (\ref{W7}) the only $N$ dependence is the factor
 $$
 {(-N + 1)_j \over N^j} = (-1)^j \bigg ( 1 -\frac{j(j+1)}{2N} +  O \Big ( {1 \over N^2}  \Big )\bigg ),
 $$
 where the expansion follows from (\ref{eq2}). This result, valid for fixed $j$, can nonetheless be substituted in
 the summation since the factor in the summand $(-N+1)_j/j! N^j$ is a rapidly decaying function of $j$. Doing this shows
 \begin{multline*}
  \sum_{j=0}^\infty {(-x)^j \over j! \prod_{s=1}^M (a_s + 1)_j}  \bigg ( 1 - \frac{j(j+1)}{2N}  +  O \Big ( {1 \over N^2}  \Big ) \bigg ) \\ =
\bigg (   1 - {1 \over 2 N} \Big ( x {d \over dx} + \Big ( x {d \over dx} \Big )^2 \Big ) \bigg ) \, {}_0 F_M  \bigg ( \begin{array}{cc} -  \\ a_1 + 1,\dots, a_M + 1 \end{array} \Big | -x \bigg )  +  O \Big ( {1 \over N^2}  \Big ),
\end{multline*}
with the bound on the RHS uniform for $x \in [0,s]$.

In relation to (\ref{W8}), after multiplying through the prefactor $1/ N!$ inside the integrand of the integral (\ref{W5}) defining the LHS,  we see
the only dependence on $N$ is the factor
$$
{\Gamma(N + s +1) \over N^s \Gamma(N+1)} = 1 + {s(s+1) \over 2N}  + O \Big ( {1 \over N^2}  \Big ),
$$
 where the expansion follows from (\ref{eq2}). The result (\ref{W8}) now follows by noting
  \begin{multline*}
 {1 \over 2 \pi i }   \int_{c - i \infty}^{c + i \infty} {1 \over \Gamma(s+1) }  \prod_{j=1}^M \Gamma(a_j - s) \bigg (  1 + {s(s+1) \over 2N}  +  O \Big ( {1 \over N^2}  \Big )  \bigg )  \, x^{s} \, ds
   \\ = \bigg ( 1 + {1 \over 2 N} \Big ( x {d \over dx} +
 \Big (  x {d \over dx} \Big )^2 \Big ) \bigg )  \MeijerG{M}{0}{1}{M+1}{-}{a_1,\ldots,a_M,0}{x } +  O \Big ( {1 \over N^2}  \Big ),
 \end{multline*}
 and arguing in relation to the error term as above.
\end{proof}

Substituting the results of Proposition \ref{P3.1} in (\ref{W4}) with $n=N-1$ and in (\ref{W5}) with $n=N$, then substituting in (\ref{Me3}) shows
\begin{multline}\label{FG}
{1 \over N} K_N(x/N,y/N) \\
= \int_0^1 \bigg ( 1 - {1 \over 2 N} \Big ( x {d \over dx} +
 \Big (  x {d \over dx} \Big )^2 \Big ) \bigg ) F(xt)
 \bigg ( 1 + {1 \over 2 N} \Big ( y {d \over dy} +
 \Big (  y {d \over dy} \Big )^2 \Big ) \bigg ) G(yt) \, dt +   O \Big ( {1 \over N^2}  \Big ),
 \end{multline}
 where $F$ denotes the function $ {}_0 F_M$ in (\ref{W7}) and $G$ denotes the function $G_{1,M+1}^{M,0}$ in (\ref{W8}).
 Note that the error bound from asymptotic forms in Proposition \ref{P3.1} persist because the error bounds therein
 are uniform with respect to $x,y$ when these variables are restricted to a compact set; see \cite{Bo16} on this
 point in relation to (\ref{pk5}).

 Independent of the details of these functions, the structure (\ref{FG}) permits simplification.

 \begin{proposition}\label{P3.2}
 The expression (\ref{FG}) has the simpler form
 \begin{equation}\label{FG1}
{1 \over N} K_N(x/N,y/N)
=  \int_0^1  F(xt)  G(yt) \, dt  - {1 \over 2N} \Big ( x {\partial \over \partial x} -  y {\partial \over \partial y} \Big ) F(x) G(y) +   O \Big ( {1 \over N^2}  \Big ).
 \end{equation}
 \end{proposition}

 \begin{proof}
 At order $1/N$ the RHS of (\ref{FG}) reads
 $$
 - {1 \over 2 N}   \int_0^1 G(yt)  \Big (  x {d \over dx} + \Big ( x {d \over dx} \Big )^2 \Big ) F(xt) \, dt
 + {1 \over 2 N}   \int_0^1 F(xt)  \Big (  y {d \over dy} + \Big ( y {d \over dy} \Big )^2 \Big ) G(yt) \, dt.
 $$
 In this expression, both the derivatives with respect to $x$, and the derivatives with respect to $y$ can
 be replaced by derivatives with respect to $t$. Performing one integration by parts for each of
 the terms involving the second derivative, (\ref{FG1}) results.
 \end{proof}

 Recalling (\ref{1.5a}), we see from (\ref{FG1}) that in general for  products of Laguerre unitary ensembles, the pointwise rate of
convergence to the hard edge limiting $k$-point correlation is $O(1/N)$. On the other hand, as noted in the text around
(\ref{pk6}), earlier works
\cite{EGP16,Bo16,PS16,HHN16,FT19} have demonstrated that for the Laguerre unitary ensemble itself (the case $M = 1$),
with the hard edge scaling variables as used in (\ref{FG1}), and with the Laguerre parameter $a=0$, the convergence
rate is actually $O(1/N^2)$. Moreover, these same references found that the  $O(1/N^2)$ rate holds for general
Laguerre parameter $a>-1$ if each $N$ on the LHS of (\ref{FG1}) is replaced by $N+a/2$.

From the viewpoint of (\ref{FG1}), the special feature of the case $M=1$ is that then $F$ and $G$ are related
by
 \begin{equation}\label{FGa}
G(x) = x^a F(x),
 \end{equation}
as follows from the final paragraph of Section \ref{S2.1}. The term $O(1/N)$ in (\ref{FG1})
can therefore be written to involve only $F$,
 \begin{equation}\label{FG2}
 - {1 \over 2N} y^a   \Big ( -a F(x) F(y) +  \Big ( x {\partial \over \partial x} -  y {\partial \over \partial y} \Big ) F(x) F(y)  \Big )  \Big |_{M=1} .
 \end{equation}
 Substituting in (\ref{FG1}), then substituting the result in (\ref{1.5a}), we factor $x_l$ from each column to
 effectively remove $y^a$ from (\ref{FG2}). The term involving partial derivatives in the latter is then antisymmetric,
 and so does not contribute to an expansion of the determinant at order $1/N$, telling us that
  \begin{align}\label{FG3}
 &{1 \over N^k}  \rho_{(k)}\Big ( {x_1 \over N} ,\dots, {x_k \over N} \Big ) \Big |_{M=1}
\nonumber  \\
 & \quad = \prod_{l=1}^k x_l^a \det \Big [ \Big ( \int_0^1 t^a F(x_jt) F(x_l t) \, dt + {a \over 2N}F(x_j) F(x_l) \Big )  \Big |_{M=1}  \Big ]_{j,l=1}^k
  + O \Big ( {1 \over N^2}  \Big ) \nonumber \\
 & \quad =  \det \Big [ \Big ( \int_0^1 \tilde{F}(x_jt) \tilde{F}(x_l t) \, dt + {a \over 2N}\tilde{F}(x_j) \tilde{F}(x_l) \Big )  \Big |_{M=1}  \Big ]_{j,l=1}^k
  + O \Big ( {1 \over N^2}  \Big ),
   \end{align}
   where $\tilde{F}(x) = x^{a/2} F(x)$, and the second equality follows from the first by multiplying each row $j$ by $x_j^{a/2}$ and each column
   $k$ by $x_k^{a/2}$. In this latter form the kernel is symmetric. Comparison with (\ref{1.2}) and (\ref{1.3}) then
   shows
   $$
   \tilde{F}(x)  \Big |_{M=1} = J_a( \sqrt{4x} ), \qquad \int_0^1 \tilde{F}(xt) \tilde{F}(yt) \, dt \Big |_{M=1} = 4 K^{\rm hard}(4x,4y)
   $$
   (the reason for the factors of 4 comes from the choice of hard edge scaling $x \mapsto x/4N$ in  (\ref{1.2}), (\ref{1.3})
   rather than $x \mapsto x/N$ as in (\ref{FG3})). This is in agreement with the references cited above relating
   to the hard edge expansion of the Laguerre unitary ensemble correlation kernel up to and including the
   $O(1/N)$ term, and so has the property
   that upon replacing $N$ by $N+a/2$ on the LHS, the convergence has the optimal rate of $O(1/N^2)$.

   \subsection{Laguerre Muttalib-Borodin model}\label{S3.2}
   The Laguerre Muttalib-Borodin model \cite{Mu95,Bo98,FW15,Zh15}, defined as the eigenvalue PDF
   proportional to
   \begin{equation}\label{MB1}
   \prod_{l=1}^N x_l^a e^{-x_l}  \prod_{1 \le j < k \le N} (x_j - x_k ) (x_j^\theta - x_k^\theta),
    \end{equation}
   with each $x_l$ positive is, with $\theta = M$ and upon the change of variables $x_l \mapsto x_l^{1/\theta}$,
   known to be closely related to the product of $M$ matrices from the LUE. Specifically, there is a choice
   of the Laguerre parameters $a_l$ for which the joint PDF of the latter reduces to this transformation of
   (\ref{MB1}) \cite{KS14}. In particular, it follows that in the case $\theta = M$ at least, (\ref{MB1}) corresponds to a
   P\'olya ensemble.
   In fact it is known from \cite{KK16} that (\ref{MB1}) is an example of a P\'olya ensemble for general $\theta > 0$.
  We can thus make use of the theory of Section \ref{S2.2} to study the hard edge expansion of the correlation
  kernel.

  The normalised weight function corresponding to (\ref{MB1}) after the stated  change of variables is
  \begin{equation}\label{MB2}
  w^{({\rm MB}, L)}(x) = {1 \over \theta \Gamma(a+1)} x^{-1 + (a+1)/\theta} e^{- x^{1/\theta}},
  \end{equation}
  which has Mellin transform
    \begin{equation}\label{MB3}
    {\mathcal M} [ w^{({\rm MB}, L)} ](s) = {\Gamma(\theta(s-1) + a + 1) \over \Gamma (a + 1) }.
 \end{equation}
 Hence the polynomials $p_n(x)$ in (\ref{Me1}) read
   \begin{equation}\label{MB4}
   p_n^{({\rm MB}, L)}(x) = (-1)^n \Gamma(\theta n + a + 1) \sum_{j=1}^n { (-n)_j x^j \over j! \Gamma(\theta j + a + 1)},
 \end{equation}
 first identified in the work of Konhauser \cite{Ko67}.

 Taking the inverse Mellin transform of (\ref{MB3}) gives the integral form of the weight,
 $$
  w^{({\rm MB}, L)}(x) = {1 \over  \Gamma(a+1)}  {1 \over  2 \pi i} \int_{c - i \theta}^{c + i \theta} \Gamma(-\theta (s + 1) + a + 1) x^s \, ds ,
  $$
  valid for $c > 0$. Using this in (\ref{Me2}) shows
    \begin{equation}\label{MB5}
    q_n^{({\rm MB}, L)}(x)  = {(-1)^n \over n! \Gamma(\theta n + a + 1)}
     {1 \over  2 \pi i }   \int_{c - i \theta}^{c + i \theta} { \Gamma(s+n+1) \over \Gamma(s) } \Gamma(-\theta (s + 1) + a + 1) x^s \, ds.
  \end{equation}

 The dependence on $n$ in the summand of (\ref{MB4}) and integrand of (\ref{MB5})  is precisely the same as in
 (\ref{W4}) and (\ref{W5}) respectively. Applying the working of Proposition \ref{P3.1} then gives hard
 edge asymptotics that is structurally identical to $p_n(x)$ and $q_n(x)$ for products of Laguerre ensembles.
 From this we conclude a formula structurally identical to (\ref{FG1}) for the hard edge asymptotics of the kernel.

 \begin{proposition}
 Define
 $$
 \tilde{p}_n^{({\rm MB}, L)}(x) = {(-1)^n \over \Gamma(\theta n + a + 1)} p_n^{({\rm MB}, L)}(x), \qquad
  \tilde{q}_n^{({\rm MB}, L)}(x) = {(-1)^n  \Gamma(\theta n + a + 1)} q_n^{({\rm MB}, L)}(x).
  $$
  Also define
  $$
  F^{({\rm MB},L)}(x) = \sum_{j=0}^\infty {x^j \over j! \Gamma(\theta j + a + 1)}, \qquad
 G^{({\rm MB},L)}(x) =   {1 \over 2 \pi i}    \int_{c - i \theta}^{c + i \theta} { \Gamma(-\theta (s + 1) + a + 1) \over \Gamma (s) } x^s \, ds.
 $$
 We have
 \begin{align*}
 \tilde{p}_{N-1}^{({\rm MB}, L)}(x/N) & = \bigg ( 1 - {1 \over 2 N} \Big ( x {d \over dx} + \Big (  x {d \over dx} \Big )^2 \Big )   \bigg ) F^{({\rm MB},L)}(x)  + O \Big ( {1 \over N^2} \Big ) \\
  \tilde{q}_N^{({\rm MB}, L)}(x/N) & = \bigg ( 1 + {1 \over 2 N} \Big ( x {d \over dx} + \Big (  x {d \over dx} \Big )^2 \Big )  \bigg ) G^{({\rm MB},L)}(x)  + O \Big ( {1 \over N^2} \Big ) ,
  \end{align*}
  and furthermore
  \begin{multline*}
  {1 \over N} K_N^{({\rm MB}, L)}(x/N,y/N) \\
=  \int_0^1  F^{({\rm MB}, L)}(xt)  G^{({\rm MB}, L)}(yt) \, dt  - {1 \over 2N} \Big ( x {\partial \over \partial x} -  y {\partial \over \partial y} \Big ) F^{({\rm MB}, L)}(x) G^{({\rm MB}, L)}(y) +   O \Big ( {1 \over N^2}  \Big ).
\end{multline*}
  \end{proposition}

  As in the discussion following Proposition \ref{P3.1}, this tells us that the rate of convergence to the hard edge scaled limit of the
  $k$-point correlation is $O(1/N)$, with the case $\theta = 1$ (corresponding to the LUE) an exception, where by appropriate choice
  of scaling variables, the rate is $O(1/N^2)$.

  \subsection{Products of Laguerre ensembles and inverse Laguerre ensembles}\label{S3.3}
  In the guise of the square singular values for the product of complex Gaussian matrices, times the inverse of a further
  product of complex Gaussian matrices, the study of the eigenvalues of a product of Laguerre ensembles and  inverses
  was initiated in \cite{Fo14}. This was put in the context of P\'olya ensembles in \cite{KS14}. Moreover, in the case that there
  are equal numbers of matrices and inverse matrices, such product ensembles can be related to a single weight function,
  as we will now demonstrate.  The essential point
  is that the eigenvalues of $X_{b_1}^{-1} X_{a_1}$, where $X_{a_1}$, $X_{b_1}$  has eigenvalues from the Laguerre unitary ensemble
  has eigenvalue PDF proportional to (see e.g.~\cite[Exercises 3.6 q.3]{Fo10})
  \begin{equation}\label{Kq1}
  \prod_{l=1}^N {x_l^{a_1} \over (1+ x_l)^{b_1 + a_1 + 2N}} \prod_{1 \le j < k \le N} ( x_k - x _j)^2
  \end{equation}
  and that this in turn is an example of a  P\'olya ensemble (\ref{2.3}) with
   \begin{equation}\label{Kq2}
   w^{(\rm I)}(x) = {x^{a_1} \over  (1 + x)^{b_1 + a_1+N+1}} \chi_{x > 0}
    \end{equation}
    (here the superscript (I) indicates `inverse').
  Structurally, a key distinguishing feature relative to the weight (\ref{Lw}) is that (\ref{Kq2}) depends on $N$.
  After normalising (\ref{Kq2}), proceeding as in the derivation of (\ref{W1}) shows the weight function
  for the P\'olya ensemble of the corresponding product ensemble is
   \begin{equation}\label{Kq3}
   \mathcal M [   w^{({\rm I}, M)}](s)  =  \prod_{l=1}^M {  \Gamma(a_l+s) \Gamma(b_l + N  + 1- s) \over \Gamma(a_l+1) \Gamma(b_l + N) }.
      \end{equation}
   Use of (\ref{Kq3}) in  (\ref{Me1}) shows
   \begin{multline}\label{Kq4}
{(-1)^n  \over      \prod_{l=1}^M    \Gamma(a_l+n+1) \Gamma(b_l + N  - n ) }
 p_n^{({\rm I}, M)}(x)
 = \sum_{j=0}^n  {(-n)_j \over j!}   {x^j \over  \prod_{l=1}^M    \Gamma(a_l+j+1) \Gamma(b_l + N  - j )  } .
 \end{multline}
 Further, using  (\ref{Kq3}) to write $w^{(\rm I)}(x) $ as an inverse Mellin transform shows from (\ref{Me2}) that
    \begin{multline}\label{Kq5}
  {(-1)^n    \over  \prod_{l=1}^M    \Gamma(a_l+n+1) \Gamma(b_l + N  - n ) }q_n^{({\rm I}, M)}(x)  \\ = {1 \over 2 \pi i}  {1 \over n!} \int_{c - i \theta}^{c + i \theta}
   { \Gamma(s+n) \over \Gamma(s) }  \Big (  \prod_{l=1}^M
   \Gamma(a_l-s) \Gamma(b_l + N  +1 + s)  \Big ) x^s \, ds.
 \end{multline}

  Proceeding as in the derivation of Proposition \ref{P3.1}, and making use in particular of the asymptotic formula
  (\ref{eq1}) for the ratio of two gamma functions, the large $N$ forms of (\ref{Kq4}) and (\ref{Kq5}) as relevant to (\ref{Me3})
 can be deduced. This allows for the analogue of (\ref{FG}) to be deduced, which then proceeding as in the derivation of
 Proposition \ref{P3.2} gives the analogue of (\ref{FG1}).

 \begin{proposition}\label{P3.4}
 Denote the LHS of (\ref{Kq4}) with $n = N -1$, and multiplied by $\prod_{l=1}^M \Gamma( N + b_l)$, by
 $\tilde{p}_{N-1}^{({\rm I}, M) }(x)$,  and let $F$ be specified as below (\ref{FG}). Also, denote the LHS of
(\ref{Kq5}) with $n = N $, and divided by $\prod_{l=1}^M \Gamma( N + b_l)$, by  $\tilde{q}_{N}^{({\rm I}, M) }(x)$,
and let $G$ be as specified below (\ref{FG}). We have
\begin{equation}
p_{N-1}^{({\rm I}, M)}\Big ( {x \over N^{M+1}} \Big )
= \bigg ( 1   - {1 \over 2 N} \bigg (
\Big ( 1 + M - 2 \sum_{l=1}^M b_l  \Big ) x {d \over dx} +
(1 + M) \Big ( x {d \over d x} \Big )^2  + O\Big ( {1 \over N^2} \Big ) \bigg ) \bigg ) F(x),
\end{equation}
\begin{multline}
 {1 \over N^M}  q_{N}^{({\rm I}, M)}\Big ( {x \over N^{M+1}} \Big ) \\
= \bigg ( 1 +  {1 \over N} \sum_{l=1}^M  b_l  +{1 \over 2N}  \bigg (
\Big ( 1 + M + 2 \sum_{l=1}^M b_l  \Big ) x {d \over dx} +
(1 + M) \Big ( x {d \over d x} \Big )^2  + O\Big ( {1 \over N^2} \Big )  \bigg ) \bigg ) G(x)
\end{multline}
and
\begin{multline}\label{3.26}
{1 \over  N^{M+1}}  K_N \Big ( {x \over N^{M+1}}, {y \over N^{M+1}} \Big ) = \int_0^1 F(xt) G(yt) \, dt \\
 - {1 \over 2 N} (1 + M) \Big ( G(y) x {d \over dx} F(x) - F(x) y {d \over d y} G(y) \Big )
+ {1 \over  N} \Big ( \sum_{l=1}^M b_l  \Big ) F(x) G(y) + O \Big ( {1 \over N^2} \Big ).
\end{multline}
\end{proposition}

The expansion (\ref{3.26}) shows that in general the leading correction to the
hard edge scaled limit of the $k$-point correlation
in the case of $M$ products of random matrices formed from the multiplication of a
Laguerre unitary ensemble and inverse Laguerre unitary ensemble is $O(1/N)$.
However, as for products studied in Section \ref{S3.1}, the case $M=1$ is special,
as then the relation (\ref{FGa}) between $F$ and $G$ holds. The $O(1/N)$ term in
(\ref{3.26}) the simplifies to read
\begin{equation}\label{3.27}
{1 \over N} y^a \bigg ( (a_1 + b_1)  F(x) F(y)  - \Big ( x {\partial \over \partial x} - y {\partial \over \partial y} \Big ) F(x) F(y) \bigg ) \bigg |_{M=1}
\end{equation}
Proceeding now as in the derivation of (\ref{FG3}), and with the same meaning of $\tilde{F}$ used
therein, we thus have
 \begin{align}\label{FG3a}
 &{1 \over N^{2k}}  \rho_{(k)}\Big ( {x_1 \over N^2} ,\dots, {x_k \over N^2} \Big ) \Big |_{M=1}
\nonumber  \\
 & \quad = \prod_{l=1}^k x_l^{a_1} \det \Big [ \Big ( \int_0^1 t^{a_1} F(x_jt) F(x_l t) \, dt + {a_1 + b_1  \over N}F(x_j) F(x_l) \Big )  \Big |_{M=1}  \Big ]_{j,l=1}^k
  + O \Big ( {1 \over N^2}  \Big ) \nonumber \\
 & \quad =  \det \Big [ \Big ( \int_0^1 \tilde{F}(x_jt) \tilde{F}(x_l t) \, dt + {a_1 + b_1  \over N}\tilde{F}(x_j) \tilde{F}(x_l) \Big )  \Big |_{M=1}  \Big ]_{j,l=1}^k
  + O \Big ( {1 \over N^2}  \Big ).
   \end{align}
  As in the discussion below (\ref{FG3}), it follows that if on the LHS $N$ is replaced by $N + (a_1 + b_1)/2$, the convergence to
  the hard edge limit has the optimal rate of $O(1/N^2)$.

  \begin{remark}
  1.~Changing variables $x_l = y_l/(1 - y_l)$, $0 < y_l < 1$ in (\ref{Kq1}) gives the functional form
  \begin{equation}\label{Ja}
  \prod_{l=1}^N y_l^{a_1} ( 1 - y_l)^{b_1} \prod_{1 \le j < k \le N} (y_k - y_j)^2,
  \end{equation}
  which up to proportionality is the eigenvalue PDF for the Jacobi unitary ensemble. In the recent work
  \cite{MMM19} the corrections to the hard edge scaled limit of the distribution of the smallest
  eigenvalue have been analysed, with results obtained consistent with (\ref{FGa}). In Appendix
  A we present a large $N$ analysis of this distribution for the Jacobi $\beta$-ensemble (the Jacobi unitary ensemble
  is the case $\beta = 2$) for general even $\beta$. \\
  2.~The case $b_1 = 0$ of the Jacobi unitary ensemble is closely related to the Cauchy two-matrix
  model \cite{BGS14}. The latter is determinantal, but since the PDF consists of two-components,
  the determinant has a block structure. Nonetheless, each block can be expressed in terms of just
  a single correlation kernel. The hard edge scaling of the latter has been undertaken in \cite{BGS14},
  with a result analogous to (\ref{FG3a}) with $b_1 = 0$ obtained. Closely related to the Cauchy two-matrix
  matrix model is the Bures ensemble, as first observed in \cite{BGS09}, and further developed
  in \cite{FK16}, with a Muttalib-Borodin type extension given in \cite{FL19}. Since the elements of the
  correlation kernel for the Bures ensemble (which is a Pfaffian point process)  are given in terms of
  the correlation kernel for the Cauchy two-matrix matrix model, it follows that by tuning the scaling
  variables at the hard edge, an optimal convergence rate of $O(1/N^2)$ can be achieved.
   \\
    3.~A Muttalib-Borodin type generalisation of (\ref{Kq1}) is known \cite[Jacobi prime case]{FI18}. Working
  analogous to that of Section \ref{S3.2} could be undertaken, although we refrain from doing that here.
  It would similarly be possible to obtain the analogue of Proposition \ref{P3.4} for the singular values
  of products of truncations of unitary ensembles \cite{KKS15}, which we know from \cite{KK16} can
  be cast in a P\'olya ensemble framework as products of Jacobi unitary ensembles.
  \end{remark}

   \subsection*{Acknowledgements}
	This research is part of the program of study supported
	by the Australian Research Council Centre of Excellence ACEMS. We thank Mario Kieburg for
	feedback on a draft of this work.


\section{Introduction}
Microscopic heat engines and refrigerators have been in the focus for the last two decades due to their possible applications across various fields. Micro-engines can find potential uses in the field of nanomachines and nanobots, which in turn would find applications in many areas, especially in medicine \cite{vyas2014}.
An efficient microscopic engine needs to be developed in order to power such small machines. Several models have been proposed in the literature and have been realized experimentally, that would act as tiny heat engines \cite{sei08_epl,saha2019,lahiri2020,roldan2016,bechinger2012,lutz2016}. In a recent experiment \cite{sood2016},  using a trapped colloidal particle in a bacterial bath, the authors have demonstrated that the engine can attain much higher efficiencies than the one that uses normal thermal baths, even in the quasi-static limit.
A nice review on microscopic heat engines has been provided in \cite{roldan2017}. The extension to quantum dynamics has been done \cite{bender2000,quan2007}.

 It was shown in \cite{saha2019,lahiri2020} that an engine using an active Orntein-Uhlenbeck particle (AOUP) trapped in a harmonic potential as its working system turns out to be more efficient than its passive counterpart. The activity was introduced by means of an exponentially correlated noise that induced finite persistence in the motion of the particle.

 Now, if such engines are driven in a time-reversed manner, then they can lead to microscopic refrigerators \cite{rana2016,pal2016}. Such refrigerators have also been modelled by using a sawtooth potential in presence of an external drive and have been studied in \cite{ai2006,lin2009,chen2010}. The basic principle of ordered motion in this case was studied in \cite{leo2010}.
 Refrigeration in presence of velocity-dependent feedback control has been studied in \cite{liang2000,kim2004,kim2007}. In \cite{liang2000}, the authors demonstrated how the accuracy of an AFM cantelever can be enhanced by reducing its thermal noise. Another potential application is in intramolecular cooling which, if achieved, can be beneficial in selecting certain reaction channels over others or for increasing the efficiency of catalysis \cite{arxiv2008}. Study of self-contained microscopic heat engines and refrigerators have been done for quantum systems \cite{linden2010,brunner2014}.

 In this work, we use an active run-and-tumble (RT) particle as the working system of the engine and the refrigerator, unlike the active Ornstein-Uhlenbeck particle used in \cite{lahiri2020}. The RT particle is placed in a harmonic trap, whose stiffness parameter is varied cyclically with time, so as to mimic the steps of a macroscopic Stirling engine \cite{bechinger2012}.
RT particles provide an alternative model for activity, which is considered to be a more realistic model for describing motion of flagellar cells like the E. Coli in a stimulus-free environment \cite{berg1972}. They swim forward while their flagella spin counter-clockwise. However, once the direction of spin becomes clockwise, the flagella unbundle and the cell undergoes a rotation about its centre of mass. In \cite{condat2005,buceta2017,buceta2018}, the authors have provided  Langevin equations that capture the statistics of such motion. An effective Smoluchowski equation for RT motion in presence of taxis effects was developed in \cite{schnitzer1993}.
The RT motion in one-dimension has been modelled by a Langevin equation involving a telegraphic noise in \cite{dhar2019,chaki2019}.
The dynamics has been contrasted with that of the active Brownian particles (ABPs) in \cite{cates2013}. For a recent review on the physics of such microswimmers, see \cite{elgeti2015}. Recently, it has been shown that RT particles can be prepared artificially by using synthetic Janus particles subjected to an electric field \cite{sano2017}.

We explore the behaviour of the RT particle in order to study its effects on the performance of thermal engines and refrigerators. An enhanced performance as compared to the passive engine/refrigerator should be of considerable practical significance.
At first we show that the results of our simulations are in agreement with the analytical results of \cite{soto2020,condat2005}. Then we use simulations to study the system under the more general conditions where the stiffness parameter of the trapping potential changes periodically with time.
In addition to the usual RT motion, the particle is subjected to thermal noise due the presence of the two heat baths at different temperatures that are required for the working of the engine/refrigerator. We closely follow the approach of \cite{condat2005} in modelling our dynamics for an active particle in a thermal environment. The stochasticity of the system entails that the particle behaves in an engine mode or refrigerator mode only in a suitably chosen range of parameters \cite{saha2019,rana2016}.  We observe that although the engine mode is generally less efficient than that of the passive engine, the same system driven under a refrigerator protocol can yield enhanced performance. This has qualitative similarity with \cite{wijland2017_entropy}, where non-Gaussian noise was used to model activity. This implies that the nature of activity has a strong influence of the behaviour of the qualitative behaviours of the thermal machines operating out of equilibrium.

In section 2, we describe the detailed model for the active engine, both in one and two dimensions. In section  3, we provide comparison of our numerical simulations of an RT particle with known analytical results, and then study numerical results for the efficiency of the active engines and their comparison with the corresponding passive ones. Section 4 describes the model for the refrigerator in one and two dimensions. Section 5 discusses the results for the active particle subjected to the refrigerator protocol in one and two dimensions. Finally, the conclusions are provided in section 6.

\section{The Engine Model} \label{sec:eng_model}

An RT particle proceeds through a series of ``runs'' followed by ``tumbles'', whose motion was first described in \cite{berg1972}. During the runs, the particle does not change direction and follows the path dictated by deterministic dynamics. However, at the end of a single run, the particle suddenly changes its direction randomly, and then sets into motion for another run along this newly chosen direction. The particle is thus said to tumble and choose a new direction in-between two run-times. Since the tumble times are typically an order of magnitude smaller than the typical run-times \cite{buceta2017}, we will approximate it to zero in our analysis. A typical trajectory of such a particle in absence of thermal noise is shown in figure \ref{fig:traj2d}, in the space of normalized position coordinates (see the discussion below Eq. \eqref{theory}). The symbols appear after every time step of integration. The persistence in motion over each runtime can be clearly observed.

\begin{figure}[!ht]
\centerline{\psfig{file=fig1.eps,width=3.0in}}
\vspace*{8pt}
\caption{A typical trajectory of a 2d RT particle. The symbols represent individual time steps.}
\label{fig:traj2d}

\end{figure}

The run times and tumble times have been shown to roughly follow  exponential distributions \cite{berg1972,berg1974}.
We therefore sample our run-times from the distribution
\begin{align}
 P(\tau_r) &= \frac{1}{\tau_p} e^{-t/\tau_p},
 \label{expdist}
\end{align}
where $\tau_p\equiv \langle\tau_r\rangle$ is the average run-time, which we refer to as the persistence time.
An exact mathematical description is lacking, although the asymptotic behaviour (in the limit of large time) can be shown to reach the diffusive regime \cite{soto2020}. Tumbling characteristics can broadly be of two types:
\begin{enumerate}
 \item The particle undergoes complete reorientation (it can  choose all possible deflection angles with equal probability).

 \item The particle undergoes partial reorientation, where the probabilities of larger deflections are smaller. In particular, we choose the cosine distribution (see Eq. \eqref{cosine_distribution} below).
\end{enumerate}
Our dynamics consists of the RT motion of the particle, as well as the thermal motion generated due the finite temperature of the medium in which the particle is present. We now consider the one dimensional model in detail, and compare the results with the two dimensional one.

\subsection{1d model for RT engine}
\label{sec:1d model}

The particle experiences a potential $V(x)=k(t)x^2/2$ created by an optical trap.
The harmonic trap undergoes four steps corresponding to a Stirling cycle \cite{bechinger2012,lahiri2020}, given below. Figure \ref{fig:stiffness}(a) shows the variation in the shape of the potential during the cycle, while functional form of $k(t)$ is shown in figure \ref{fig:stiffness}(b).
\begin{enumerate}
 \item {\bf Step 1 (isothermal expansion)}: the temperature is fixed at the value $T_h$ (corresponding to temperature of the hot thermal bath), and $k(t)$ varies linearly with time from the initial stiffness value $k_0$ to the final value $k_0/2$:
 \begin{align}
  k(t) &= k_{\mathrm{exp}}(t) =  k_0\left(1-\frac{t}{\tau}\right), \hspace{0.5cm} 0<t\le \tau/2
 \end{align}
 The randomness in the dynamics is caused by the thermal noise of strength $D_h = 2\gamma k_B T_h$, where $\gamma$ is the friction coefficient of the medium, and $k_B$ is the Boltzmann constant.

 \item {\bf Step 2 (isochoric temperature decrease):} The value of stiffness parameter is held constant, and the temperature is suddenly changed (in comparison to the sampling frequency, as per the set-up used in \cite{bechinger2012}) to the lower value $T_c$.

 \item {\bf Step 3 (isothermal compression):} The temperature is held constant at $T_c$, while the stiffness parameter changes from $k_0/2$ to $k_0$ as per the following linear function:
 \begin{align}
  k(t) &= k_{\mathrm{com}}(t) = k_0\left(\frac{t}{\tau}\right), \hspace{1cm} \tau/2<t\le \tau.
 \end{align}
 In this step, in addition to activity, a thermal noise of strength $D_c=2\gamma k_B T_c$ is present due to the finite temperature of the medium.

 \item {\bf Step 4 (isochoric temperature increase):} In the final step, the temperature is suddenly increased to $T_h$, keeping the stiffness parameter fixed, so as to allow the engine to complete a cycle.\\

\end{enumerate}

\begin{figure}[!ht]
\centering
\subfigure{\psfig{file=fig2a,width=2.4in}}
\hfill
\subfigure{\psfig{file=fig2b,width=2.4in}}
\vspace*{8pt}
\caption{(a) Stirling cycle under engine protocol. (b)Variation of stiffness parameter with time in the engine protocol.}
\label{fig:stiffness}

\end{figure}

In the case of total reorientation, the particle moves either along the $+x$ or $-x$ direction for a time duration $\tau_r$, and then suddenly changes direction ($+x\to -x$ or vice versa) with probability 1/2. The persistence time is chosen to be typically much smaller than the cycle time, so that each isothermal step consists of several run-and-tumble events. When the particle undergoes partial reorientation, it switches direction more reluctantly, with a probability $p<1/2$. With the switch in direction, the active force also changes sign.

During the expansion step, in order to make comparisons with our results in \cite{lahiri2020}, we keep the particle inactive. This can be achieved, for instance, by photoactivation \cite{buttinoni2012,uspal2019}.
The equation of motion is given by (see \cite{condat2005})
\begin{align}
 m\dot v &= -\gamma v - k_{\rm exp}(t)x + (\sqrt{D_h})\xi(t), \hspace{1cm} 0\le t<\tau/2
 \label{1d_RT_exp}
\end{align}
where $\xi(t)$ is a normally distributed white noise: $\langle\xi(t)\rangle=0$, $\langle\xi(t)\xi(t')\rangle=\delta(t-t').$ In a few cases, however, we will compare the results of this model with the one in which the activity is present throughout the circle.

During the compression step, the equation of motion during the $i^{th}$ run-time is given by
\begin{align}
 m\dot v_i &= -\gamma v_i - k_{\rm com}(t)x + F + (\sqrt{D_c})\xi(t),\nonumber\\
 &\hspace{1.5cm} \tau/2 \le t_i \le t < t_i+\tau_{r,i}\le\tau.
 \label{1d_RT_com}
\end{align}
$\tau_{r,i}$ represents the $i^{th}$ run-time, while the active force $F$ propels the particle in the chosen direction during a run \cite{condat2005}. Here, $t_i$ is the initial time for the $i^{th}$ run, which ends at time $t_i+\tau_{r,i}$. If activity is present in the expansion half as well, then a similar equation of motion will take the place of Eq. \eqref{1d_RT_exp}.

\subsection{2d model}
\label{sec:2d_model}

In this case, the particle is trapped in a 2d harmonic potential, $V(x,y)=k(t)(x^2+y^2)/2$. We can describe the position vector of the particle by ${\bf r}=x\hat x+y\hat y$ during a run-time ($\hat x$ and $\hat y$ are the unit vectors in the $x$ and $y$ directions, respectively), which changes to ${\bf r}'=x'\hat x+y'\hat y$, and the angle $\theta$ between ${\bf r}$ and ${\bf r}'$ is the deflection attained due to the sudden tumble at the end of a run-time.

For the simple case of total reorientation, we use a uniform random number $\eta$ at each tumble event in the range $\eta\in[-\pi,\pi]$, and allow the value of the angle (with the $x-$axis, for convenience) $\theta$ to change to $\theta+\eta$.
For the case of partial reorientation, we choose the distribution of the deflection to be of the form
\begin{align}
P(\eta)=\frac{1}{4}\cos\left(\frac{\eta}{2}\right),
\label{cosine_distribution}
\end{align}
where the range is $\eta\in[-\pi,\pi]$. The method of obtaining the distribution has been given in Appendix A
The stiffness parameter $k(t)$ undergoes the same set of steps as mentioned in section \ref{sec:1d model}.

Let the total velocity be given by ${\bf v}=v_x\hat x + v_y\hat y$. Then the equations of motion are given by
\begin{align}
 m\dot v_{x} &= -\gamma v_{x} - k_{\rm exp}(t)x + (\sqrt{D_h})\xi(t), \hspace{0.5cm} 0<t\le\tau/2;\nonumber\\
 m\dot v_{x,i} &= -\gamma v_{x,i} - k_{\rm com}(t)x + F_x + (\sqrt{D_c})\xi(t), \nonumber\\
 &\hspace{3.5cm} \tau/2<t_i< t \le t_i+\tau_{r,i}\le\tau,
 \label{2d_RTmodel}
\end{align}
where $\tau_{r,i}$ is the $i^{th}$ run-time. ${\bf F}$ is a self-propelling force that is directed along the axis of the particle, having components $F_x=F\cos\theta$ and $F_y=F\sin\theta$. Similar equations can be written for motion in the $y$-direction.
\section{Results and Discussions for stochastic engine using RT particle}
\label{sec:active_engine}

\subsection{Comparison with theory}

We first test the outcomes of our simulation against theory. It was shown in \cite{soto2020} that the mean squared displacement of a two-dimensional RT particle (in the absence of any thermal bath or any potential) leads to the relations
\begin{align}
 \langle r^2\rangle_{\mathrm{tot}} &= 2v_0^2\tau_p^2\left(\frac{t}{\tau_p} + e^{-t/\tau_p}-1\right);\nonumber\\
 \langle r^2\rangle_{\mathrm{par}} &= \frac{2v_0^2\tau_p^2}{1-\sigma_1}\left[\frac{t}{\tau_p} + \frac{e^{-(1-\sigma_1 )t/\tau_p}-1}{1-\sigma_1}\right].
 \label{theory}
\end{align}
Here,  $\sigma_1 = \langle\cos(\eta)\rangle$ is the Fourier transform of the tumbling kernel (see \cite{soto2020}). The parameter $v_0$ is the constant magnitude of the total velocity. On the other hand, in presence of thermal noise, the change in steady-state variance of an RT particle (without any potential) can be obtained using the approach of \cite{condat2005}. If the friction coefficient of the thermal bath is $\gamma$ and its thermal noise strength is $D$, while the self-propelling force $F$ is set to zero for analytical convenience, then we arrive at the expression (see Appendix B for details)
\begin{align}
 \langle r^2(t)\rangle &= \left[\frac{2D\tau_p}{\gamma(m+\gamma \tau_p)}\right]t + \left[\frac{2Dm\tau_p^2}{\gamma(m+\gamma\tau_p)^2}\right]\nonumber\\
 &\times\exp\left\{-\frac{(m+\gamma\tau_p)t}{m\tau_p}\right\} - \frac{2Dm\tau_p^2}{\gamma(m+\gamma\tau_p)^2}.
 \label{var_condat}
\end{align}

We have shown these results in figure \ref{fig:analytics}, where the simulation results have been obtained by using the Heun's method of integration (see \cite{man2000}) and by averaging over $\sim 10^5$ realizations. In the figures, all times are normalized by the cycle time $\tau$. All positions are normalized by the factor $\sqrt{D_c\tau}$ ($D_h$ and $D_c$ are the diffusion constants for the hot and the cold reservoirs, respectively). This leads to a characteristic mass scale of $\gamma\tau$ and characteristic force scale of $\gamma\sqrt{D_c/\tau}$. In figure \ref{fig:analytics}(a), equations \eqref{theory} have been compared with the results of simulations for total and partial reorientations. The symbols represent results of simulations, while the solid lines are the analytical results.

Figure \ref{fig:analytics}(b) compares Eq. \eqref{var_condat} with simulation outcomes (single bath with a noise strength of $D$). In all the cases, the results of simulation are in good agreement with theoretical results. Unfortunately, these analytical treatments become intractable in presence of external time-dependent potentials. The nature of velocity correlations have been discussed in Appendix C.

\begin{figure}[!ht]

\centering
\subfigure{\psfig{file=fig3a.eps,width=2.4in}}
\subfigure{\psfig{file=fig3b.eps,width=2.4in}}
\vspace*{8pt}
\caption{(a) Plots of the analytical solution in comparison to the numerical results for total reorientation (black solid line and black open circles, respectively) and for partial reorientation (red solid line and red open triangles respectively) of a 2d RT particle. The parameters are: $m=0.2, ~\tau_p=0.1$.(b) Plots for the change in variance of the particle as a function of time in the steady state when thermal noise is present. The parameters are $m=0.2,~\tau_p=0.1,~D=1$.}
\label{fig:analytics}
\end{figure}

\subsection{Efficiency of the 1d active RT engine}

	We now proceed with the calculation of efficiency of an engine using a 1d active particle by simulating the dynamics of the engine (see Eqs. \eqref{1d_RT_exp} and \eqref{1d_RT_com}).
	 The work done on the particle is calculated by using the definition as prescribed in Stochastic Thermodynamics \cite{sek98,sekimoto}:
\begin{align}
 W(t) &= \frac{1}{2}\int_0^t \dot k(t')x^2(t')dt'.
\end{align}
The average work extracted can be calculated from the variance, defined as $\sigma(t)=\langle x^2(t)\rangle$, by means of the relation
\begin{align}
 \langle W(t)\rangle &= \frac{1}{2}\int_0^t \dot k(t')\sigma(t')dt'.
\end{align}
The average heat absorbed in the expansion step is given by the first law:
\begin{align}
 \langle Q\rangle_h &= \langle W\rangle_h - \langle\Delta E\rangle_h,
\end{align}
where the subscript $h$ implies that only the expansion step (when system is in contact with hot reservoir) is considered. The sign convention is such that the work done \emph{on} the system and the heat dissipated \emph{by} the system are positive, which is why  the extracted work and absorbed heat must be given by $-\av W$ and $-\av{Q}_h$, respectively. The average change in internal energy  during the expansion step is given by
\begin{align}
 \langle\Delta E\rangle_h &= \frac{1}{2}k(\tau/2)\sigma(\tau/2) - \frac{1}{2}k(0)\sigma(0).
 \label{DeltaEh}
\end{align}
The efficiency is obtained using the standard relation
\begin{align}
 \eta &= \frac{\av W}{\av{Q}_h}= \frac{\av W}{\av{W}_h-\av{\Delta E}_h}.
\end{align}

\begin{figure}[!ht]
\centering
\subfigure{\psfig{file=fig4a.eps,width=2.4in}}
\subfigure{\psfig{file=fig4b.eps,width=2.4in}}
\vspace*{8pt}
\caption{(a) Plot of efficiency with mass in 1d for the passive engine (black solid line), RT engine, with activity only in the compression cycle, with total reorientation (magenta dashed line) and with partial reorientation (red dotted line). The parameters are $D_{c}=1$, $D_{h}=2$, $k_0=50$, $\tau_p=0.2$, $F=0.1$.  (b)Similar plot when the activity is present in the full engine cycle.}
\label{fig:RT1_eff_m}

\end{figure}
We now compare the efficiencies of the active and passive engines, when the engine has settled into a time-periodic steady state, i.e. when the variance becomes periodic in time. We show this in figure \ref{fig:RT1_eff_m}(a), where we  observe the general trend that the efficiency of the passive engine is slightly higher than the active engine with partial reorientation, the latter being in turn slightly higher than the active engine with total reorientation. However, the differences are quite small in our case.
In figure \ref{fig:RT1_eff_m}(b), we have plotted similar graphs for the case where activity is present throughout the full cycle, and the difference between the passive and active engines is more pronounced in this case. Consequently, we conclude that unlike in the case of AOUP (see \cite{lahiri2020}, and Appendix D), the passive particle turns out to produce a more efficient engine than an RTP. This is likely due to the degradation of persistence owing to the tumble events for smaller values of $F$, and the decrease in compressibility in the compression step for higher values of $F$. The general trend can be found to be a decreasing efficiency with increase in the mass of the particle. This is expected, because more work is to be pumped into the system during the compression cycle when the particle is more massive, thereby reducing the efficiency. We have separately checked the variation of efficiency with $\tau_p$ and with $F$, but in all cases the passive engine outperforms its active counterpart.

\subsection{Efficiency of a 2d RT engine}

We obtain similar trends even when the engine is two-dimensional and follows the equations of motion as outlined in section \ref{sec:2d_model}, namely the active engines are not able to outperform the passive ones.
Thus, given that the active engine does not provide any major advantage over its passive counterpart, we next focus our attention on the refrigerator protocol.

\section{The Refrigerator Model} \label{sec:ref_model}

In this section, we study the behaviour of the RT particle undergoing partial reorientations in one and two dimensions, by reversing the engine protocol, as shown in figure \ref{fig:stiffness_ref}(b).

\begin{figure}[!ht]
\centering
\subfigure{\psfig{file=fig5a.eps,width=2.4in}}
\hfill
\subfigure{\psfig{file=fig5b.eps,width=2.4in}}
\vspace*{8pt}
\caption{(a) Stirling cycle under refrigerator protocol. (b)Variation of stiffness parameter in the  refrigerator protocol.}
\label{fig:stiffness_ref}

\end{figure}

The time-dependences are given by
\begin{align}
 k_{\mathrm{com}}(t) &= k_0\left(\frac{1}{2}+\frac{t}{\tau}\right), \hspace{1cm} 0< t\le \tau/2; \nonumber\\
 k_{\mathrm{exp}}(t) &= k_0\left(\frac{3}{2}-\frac{t}{\tau}\right) \hspace{1cm} \tau/2< t\le \tau.
\end{align}
The sequence of steps are given below. The equations of motion can be written down in the same way as in the engine protocol (see sec. \ref{sec:eng_model}).
\begin{enumerate}
 \item {\bf Isothermal compression} at temperature $T_h$, changing stiffness constant from $k_0/2$ to $k_0$. In this step, activity is switched on.
 \item {\bf Isochoric decrease of temperature} from $T_h$ to $T_c$. The stiffness parameter remains fixed at $k_0$.
 \item {\bf Isothermal expansion} at temperature $T_c$,  changing the stiffness constant from $k_0$ to $k_0/2$. The particle is inactive in this step.
 \item {\bf Isochoric increase of temperature} from $T_c$ to $T_h$. The stiffness parameter remains fixed at $k_0/2$.
\end{enumerate}

The \emph{coefficient of performance} (COP) of a refrigerator is defined as the ratio of the heat removed from the cold reservoir to the work done in the full cycle:
\begin{align}
 \mathrm{COP} &= -\frac{\av{Q}_c}{\av{W}} = -\frac{\av{W}_c-\av{\Delta E}_c}{\av{W}}.
\end{align}
The subscript $c$ refers to the fact that the quantities have been computed when the particle is in contact with the cold reservoir.

For a reversible engine, we find that
\begin{align}
 \mathrm{COP} &= \frac{T_c}{T_h-T_c} = \frac{T_c}{\eta_{\mathrm{rev}} T_h},
 \label{COP_rev}
\end{align}
where $\eta_{\mathrm{rev}} \equiv (1-T_c/T_h)$ is the Carnot efficiency in the engine mode. Thus, in the nonequilibrium regime, it would be reasonable to expect that a more efficient engine when run in reverse would give a refrigerator with smaller COP, and vice versa. Keeping this in mind, we study the thermodynamics of the active refrigerator and check whether it can have a higher value of COP than a passive refrigerator in a suitable range of parameters.

\section{Results and Discussions for Stochastic Refrigerator using RT particle} \label{sec:ref_discussions}
\subsection{COP of a 1d RT refrigerator}

\begin{figure}[!ht]
\centering
\subfigure{\psfig{file=fig6a.eps,width=2.4in}}
\subfigure{\psfig{file=fig6b.eps,width=2.4in}}
\vspace*{8pt}
\caption{(a) Plot of the coefficient of performance of the 1d refrigerator with total (magneta dashed line) and partial (red dotted line) reorientation as a function of the mass of the particle. The black solid line is the COP for a passive refrigerator. The activity is present only in the compression step. The parameters are: $F=0.1$ (for active particle), $D_h=1.25,~D_c=1,~\tau_p=0.2$. (b)Similar plot when the activity is present in the full refrigerator cycle.}
\label{fig:RT1_cop_m}
\end{figure}

In figure \ref{fig:RT1_cop_m}(a), we plot the variation of the COP with the mass for 1d active (total as well as partial reorientation) and passive particles. For a suitable range of parameters (see figure caption), the active refrigerator is observed to have a higher COP than the passive one. We will focus only on the range where the system acts as a refrigerator ($\av{Q}_c<0,~\av{W}>0$), which happens when the mass of the particle is less than approximately 0.09 in the diagram, where the other parameters are as mentioned in the caption. For sake of comparison, we also show the behaviour of COP with mass in the case of a fully active cycle (i.e. the RT motion is present in both halves of the cycle), which show a similar trend as in (a). Since the partial reorientation is observed to produce little difference from total reorientation, we will henceforth deal only with a particle that undergoes a partial reorientation.

\begin{figure}[!ht]
\centering
\subfigure{\psfig{file=fig7a.eps,width=2.4in}}
\subfigure{\psfig{file=fig7b.eps,width=2.4in}}
\vspace*{8pt}
\caption{Phase plot of COP as a function of $D_h$ and $m$ for (a) 1d passive refrigerator(b)1d active refrigerator. Parameters are $D_c=1,~k_0=50,\tau_p=0.2$.}
\label{fig:ref_phaseplot}

\end{figure}

In figure \ref{fig:ref_phaseplot}, we show the phase plot indicating the variation of the COP of the passive  and active refrigerators (figures \ref{fig:ref_phaseplot}(a) and (b), respectively) as a function of the thermal noise $D_h$ of the hot bath and of the mass $m$ of the particle. One can observe the reduced value of COP when the thermal drive is smaller, in conformity with our expectations as discussed above (if engine is inefficient, COP tends to be higher). The highest values of COP is obtained when both $m$ and $D_h$ are small. Although we find similar qualitative trends for both passive and RT particles, the COP is observed to be larger in the latter case.
Overall, we find that although the engine using RT particle as its working substance finds it hard to beat the performance of a passive engine, the situation is quite different when one deals with the corresponding refrigerators by time-reversing the protocols. The RT refrigerator in a suitably chosen range of parameters can yield much better results than a passive refrigerator.

\begin{figure}[!ht]
\centering
\subfigure{\psfig{file=fig8a.eps,width=2.4in}}
\subfigure{\psfig{file=fig8b.eps,width=2.4in}}
\vspace*{8pt}
\caption{(a) Phaseplot of 1d active refrigerator, where the COP has been plotted as a function of the active force $F$ and the persistence time $\tau_p$. The parameters used are : $k_0=50,~D_c=1,~D_h=1.25,~m=0.01$.  (b) Variation of COP with $F$ for mass=0.01.}
\label{fig:comp_F_taup}

\end{figure}

Figure \ref{fig:comp_F_taup}(a) shows the variation in COP with the mean run-time $\tau_p$ and the active force $F$. We find that the COP generally decreases with increase in run-time. However, as a function of $F$, the behaviour turns out to be non-trivial at smaller values of $\tau_p$, where a non-monotonic variation of the COP can be observed. This region has been explored in figure \ref{fig:comp_F_taup}(b), where the maximum of COP is found roughly at around $F=0.5$, where the other parameters are as mentioned in the figure caption.

Thus, it is possible to tune the external parameters to make the refrigerator run at a high COP. The general rule would be to allow the system to function at relatively smaller values of mass and temperature difference, and control the value of the active force $F$, which turns out to be more important in determining the efficacy of performance of the refrigerator, in contrast to the persistence time $\tau_p$ that seems to produce a much smaller influence.

\subsection{COP of a 2d RT refrigerator}

We now extend our study to that of 2d RT refrigerators (see section \ref{sec:2d_model}). Figure  \ref{fig:m_cop_comp_2d} shows the variations in COP with particle mass for a 2d passive refrigerator as well as for the 2d RT refrigerator (partial and total reorientation). It is apparent that the behaviour resembles that of the 1d refrigerator (see figure \ref{fig:RT1_cop_m}). The total reorientation case seems to offer an advantage over the partial reorientation case, although the differences in COP are not significant. Similar qualitative behaviours are obtained when the activity is present in the full cycle, which we have checked separately.

\begin{figure}[!ht]
\centerline{\psfig{file=fig9.eps,width=2.4in}}
\vspace*{8pt}
\caption{Plots showing the variation in the COP with mass for (a) a passive (black solid line), (b) an RT refrigerator (partial reorientation, red dotted line) and (c) an RT refrigerator (total reorientation, magenta dashed lines) in two dimensions. The parameters used are: $F=0.1$ (for active particle), $D_c=1,D_h=1.25,~k_0=50, \tau_p=0.2$.}
\label{fig:m_cop_comp_2d}
\end{figure}

Let us now compare the functional dependence of the COP on mass and the temperature of the hot bath, for the passive and the active 2d refrigerators (see figure \ref{fig:RT2_m_Dh}). It can be observed that like in the 1d case, we generally have higher values of COP when the refrigerator is active. However, the parameter range within which the RT system acts as a refrigerator ($\av W>0$, $\av{Q_c}<0$) is smaller than in the case of the passive one.

\begin{figure}[!ht]
\centering
\subfigure{\psfig{file=fig10a.eps,width=2.4in}}
\subfigure{\psfig{file=fig10b.eps,width=2.4in}}
\vspace*{8pt}
\caption{(a)Phase plot showing variation of COP of 2d passive refrigerator with mass and the noise strength of hot bath.(b)Similar phase plot for 2d RT active refrigerator. The parameters are: $D_c=1,~k_0=50,\tau_p=0.2$.}
\label{fig:RT2_m_Dh}

\end{figure}

We finally discuss about the effect of changes in the active force $F$ and persistence time $\tau_p$ on the performance of the 2d refrigerator. The corresponding phase plot is shown in figure \ref{fig:comp_F_taup_2d}. We find that the non-monotonicity that was present for small values of persistence time for the 1d refrigerator (see figure \ref{fig:comp_F_taup}) does not show up in the 2d case. Instead, to obtain better performance, one should increase the persistence time .

\begin{figure}[!ht]
\centerline{\psfig{file=fig11.eps,width=2.4in}}
\vspace*{8pt}
\caption{Phase plot showing the dependence of COP of 2d RT refrigerator on the parameters $F$ and $\tau_p$. The parameters are: $D_c=1,~k_0=50,~m=0.01$.}
\label{fig:comp_F_taup_2d}

\end{figure}

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

Microscopic thermal machines have been studied extensively in recent years. Modelling of such machines using a trapped particle in a confining potential has been studied. Recent works show that if a heat engine at such small scales uses an active Ornstein-Uhlenbeck particle as its working system, it can lead to a higher efficiency than that using a passive Brownian particle. We study a microscopic Stirling engine using an active particle, but the activity is introduced by run-and-tumble motion of the particle. In this case,  we find that the efficiency is less than that of a passive engine.
However, if we apply the refrigerator protocol, then in a suitably chosen parameter regime, the RT refrigerator can give rise to a higher coefficient of performance as compared to its passive counterpart. Both efficiency (for the engine) and COP (for the refrigerator) are higher when the particle mass becomes smaller. An RT motion that involves partial reorientations at the tumble events gives only a minor improvement over that involving total reorientations. For the one-dimensional refrigerator, there is a clear non-monotonic behaviour of the COP when plotted as a function of the active force. It also becomes evident that depending on the nature of the activity (like active Ornstein-Uhlenbeck particle, active Brownian particle, run-and-tumble particle, etc.), the qualitative properties of the nonequilibrium thermal machines are required to be analyzed on a case-by-case basis.