\section{Introduction}
Nowadays, the optimization of thermal heat devices is receiving a
special attention because of its straightforward relation with the
depletion of energy resources and the concerns of sustainable
development. A number of different performance regimes based on
different figures of merit have been considered with special
emphasis in the analysis of possible universal and unified features.
Among them, the efficiency at maximum power (EMP) $\eta_{\rm maxP}$
for heat engines working along cycles or in steady-state is largely
the issue most studied independently of the thermal device nature
(macroscopic, stochastic or quantum) and/or the model
characteristics~\cite{esposito09, benenti10, norma10, tu12b, wu,
durmayaz, tu08, seifert07, seifert08, seifert11, gaveau11,
esposito10a, esposito10b, abe11}. Even theoretical results have been
faced with an experimental realization for micrometre-sized
stochastic heat engines performing a Stirling cycle~\cite{blickle11}.

For most of the Carnot-like heat engine models analyzed in the
finite-time thermodynamics (FTT) framework~\cite{wu, durmayaz} the
EMP regime allows for valuable and simple expressions of the
optimized efficiency, which under endoreversible assumptions
(\textit{i.e.}, all considered irreversibilities coming from the
couplings between the working system and the external heat
reservoirs through linear heat transfer laws) recover the
paradigmatic Curzon-Ahlborn value~\cite{curzon75} $\eta_{\rm maxP}=
1-\sqrt{\tau}\equiv \eta_{\rm CA}$ using $\tau \equiv T_{\rm c}/T_{\rm h}$,
where $T_{\rm c}$ and $T_{\rm h}$ denote
the temperature of the cold and hot heat reservoirs, respectively.
A significant conceptual advance was reported by Esposito {\em et al.}~\cite{esposito10c}
by considering a low-dissipation Carnot heat engine model.
In this model the entropy production in the hot (cold) heat
exchange process is assumed to be inversely proportional to the time duration of the process
(see eqs.~(\ref{eq.low.qh}) and (\ref{eq.low.qc}) for more details).
Then the reversible regime is approached in the
limit of infinite times and the maximum power regime allows us to
recover the Curzon-Ahlborn value when symmetric dissipation is
considered. This derivation does not require the assumption of any specific
heat transfer law nor the small temperature difference.
Considering extremely asymmetric dissipation limits, these authors
also derived that EMP is bounded from the lower and upper sides as
\begin{eqnarray}
\frac{\eta_{\rm C}}{2}\le \eta_{{\rm maxP}}
\le \frac{\eta_{\rm C}}{2-\eta_{\rm C}},\label{eq.emp.bounds}
\end{eqnarray}
where $\eta_{\rm C}\equiv 1-\tau$ is the Carnot efficiency. The Curzon-Ahlborn value
$\eta_{\rm CA}$ (\textit{i.e.}, the symmetric dissipation
limit) is located between these bounds. Additional studies on
various heat engine models~\cite{yan12,tu12a,tu12c} confirmed and
generalized above results.

Work by de Tom\'as {\em et al.}~\cite{carla12} extended the low
dissipation Carnot heat engine model to refrigerators and besides
proposed a unified figure of merit ($\chi$-criterion described below) focusing the attention in the
common characteristics of every heat energy converter (the working
systems and total cycle time) instead of any specific coupling to
external heat sources which can vary according to a particular arrangement.
For refrigerators, we denote the heat absorbed by the working system from the cold
heat reservoir by $Q_{\rm c}$, the heat evacuated from the working system
to the hot heat reservoir by $Q_{\rm h}$
and the work input to the working system by
$W\equiv Q_{\rm h}-Q_{\rm c}$.
For heat engines, we may inversely choose natural directions of heat transfers and work.
That is, we denote the heat absorbed by the working system from the hot heat reservoir by $Q_{\rm h}$,
the heat evacuated from the working system to the cold heat reservoir by $Q_{\rm c}$
and the work output by $W\equiv Q_{\rm h}-Q_{\rm c}$.
Then we can define this $\chi$-criterion
as the product of the converter efficiency $z$ times the heat
$Q_{\rm in}$ exchanged between the working system and the heat reservoir,
divided by the time duration of cycle $t_{\rm cycle}$:
\begin{eqnarray}
\chi=\frac{z\, Q_{\rm in}}{t_{\rm cycle}},\label{eq.chi.def}
\end{eqnarray}
where $Q_{\rm in}=Q_{\rm c}$, $z=\varepsilon \equiv Q_{\rm c}/W$ for refrigerators and
$Q_{ \rm in}=Q_{\rm h}$, $z=\eta \equiv W/Q_{\rm h}$ for heat engines.
It becomes power as $\chi=W/t_{\rm cycle}$ when applied to
heat engines and also allows us to obtain the optimized coefficient of
performance (COP) $\varepsilon_{{\rm max\chi}}$ under symmetric
dissipation conditions when applied to refrigerators~\cite{carla12}:
\begin{eqnarray}
\varepsilon_{{\rm
max\chi}}=\sqrt{1+\varepsilon_{\rm C}}-1\equiv\varepsilon_{\rm
CA},\label{eq.maxchi}
\end{eqnarray}
where $\varepsilon_{\rm C}\equiv \tau/(1-\tau)$ is the Carnot COP.
This result could be considered as the genuine counterpart of the
Curzon-Ahlborn efficiency for refrigerators. It was first obtained
in FTT for Carnot-like refrigerators by Yan and Chen~\cite{yan90}
taking as target function $\varepsilon \dot Q_{\rm c}$, where $\dot
Q_{\rm c}$ is the cooling power of the refrigerator, later and
independently  by Velasco \textit{et al.}~\cite{velasco97a,
velasco97b} using a maximum per-unit-time COP and by Allahverdyan
\textit{et al.}~\cite{allahv10} in the classical limit of a quantum
model with two n-level systems interacting via a pulsed external
field. Very recent results by Wang \textit{et al.}~\cite{tu12d}
generalized the previous results for refrigerators~\cite{carla12} and they obtained
the following bounds of the COP considering extremely asymmetric dissipation limits:
\begin{eqnarray}
0\le \varepsilon_{{\rm max\chi}} \le
\frac{\sqrt{9+8\varepsilon_{\rm C}}-3}{2}.\label{eq.cop.bounds}
\end{eqnarray}
These bounds for refrigerators play the same role that eq.~(\ref{eq.emp.bounds}) for heat engines.

A shortcoming of the all above results is the model dependence.
Thus additional research work~\cite{broeck05, borja06, borja07,borja08}
has been devoted to obtain model-independent results on efficiency and COP using the well
founded formalism of linear irreversible thermodynamics (LIT) for
both cyclic and steady-state models including explicit calculations
of the Onsager coefficients using molecular kinetic
theory~\cite{yuki09, yuki10, yuki08}. Beyond the linear regime a
further improvement was reported by Izumida and Okuda~\cite{yuki12}
by proposing a model of minimally nonlinear irreversible heat
engines described by extended Onsager relations with nonlinear
terms accounting for power dissipation. The validity of the theory
was checked by comparing the results for EMP and its bounds with
eq.~(\ref{eq.emp.bounds}).
Now, the main goal of the present paper is
to extend the application of this nonlinear irreversible
theory~\cite{yuki12} to refrigerators and analyze their COP under maximum $\chi$-condition.
To get this we first present the model and analyze the results when the
tight-coupling condition is met. Then we also study the nontight-coupling case numerically.
Finally we extend the low-dissipation Carnot refrigerator model~\cite{tu12d} to a leaky model
and show that it serves as an example of the minimally nonlinear irreversible refrigerator
by calculating its Onsager coefficients explicitly.
\section{Minimally nonlinear irreversible refrigerator model}
Refrigerators are generally classified into two types, that is, steady-state refrigerators and
cyclic ones. Our theory can be applied to both cases.
Since the working system stays unchanged for steady-state refrigerators and comes back to the original state after one-cycle for cyclic refrigerators,
the entropy production rate $\dot{\sigma}$ of the refrigerators agrees with
the entropy increase rate of the heat reservoirs as
\begin{eqnarray}
\dot{\sigma}=\frac{\dot{Q}_{\rm c}+\dot{W}}{T_{\rm h}}-\frac{\dot{Q}_{\rm c}}{T_{\rm c}}
=\frac{\dot{W}}{T_{\rm h}}+\dot{Q}_{\rm c} \left(\frac{1}{T_{\rm h}}-\frac{1}{T_{\rm c}}\right),\label{eq.entropy-refrigerator}
\end{eqnarray}
where the dot denotes the quantity per unit time for steady-state refrigerators and
the quantity divided by the time duration of cycle $t_{\rm cycle}$ for cyclic refrigerators.
From the decomposition of $\dot{\sigma}$ into the sum of the product of thermodynamic flux $J_i$ and its conjugate thermodynamic force $X_i$ as $\dot{\sigma}=J_1X_1+J_2X_2$,
we can define $J_1\equiv \dot{x}$, $X_1\equiv F/T_{\rm h}$ and $J_2\equiv \dot{Q}_{\rm c}$,
$X_2\equiv 1/T_{\rm h}-1/T_{\rm c}$ for steady-state refrigerators, where
$\dot{W}=F\dot{x}$ with $F$ and $x$ denoting the time-independent
external generalized force and its conjugate variable, respectively.
Likewise we can also define
\begin{eqnarray}
J_1\equiv \frac{1}{t_{\rm cycle}},\ X_1\equiv \frac{W}{T_{\rm h}},\label{eq.def_J1_X1}\\
J_2\equiv \dot{Q}_{\rm c}=\frac{Q_{\rm c}}{t_{\rm cycle}},\  X_2\equiv \frac{1}{T_{\rm h}}-\frac{1}{T_{\rm c}},\label{eq.def_J2_X2}
\end{eqnarray}
for cyclic refrigerators, where $X_1$ and $X_2$ play the
role of driver and driven forces, respectively.~\footnote{To establish above
election of the thermodynamic fluxes and forces from eq.~(\ref{eq.entropy-refrigerator}) we mainly
focus on the specific function of the refrigerator system (the
extracted cooling power $\dot Q_{\rm c}$ of the low-temperature reservoir)
as a consequence of the input of an external power $\dot W$~\cite{borja08}. An
alternative starting point is to express the entropy production rate
in eq.~(\ref{eq.entropy-refrigerator}) in terms of $\dot W$ and $\dot Q_{\rm h}$. If this is done the
thermodynamic fluxes and forces are exactly the same as those
obtained for a heat engine in~\cite{yuki12} with the sole difference
coming from the sign convention: $\dot{W}$ is given in (delivered to) for
a refrigerator (heat engine) while $\dot Q_{\rm h}$ is delivered to (given
in) for a refrigerator (heat engine). We have checked that this
alternative election of fluxes and forces does not change the main results in the present paper.}
\begin{figure}
\begin{center}
\includegraphics[scale=0.15]{fig1.eps}
\end{center}
\caption{Schematic illustration of the minimally nonlinear irreversible
refrigerator described by eqs.~(\ref{eq.J2.by.J1}) and (\ref{eq.J3.by.J1}).}\label{fig.1}
\end{figure}
Anyway, we assume that the refrigerators are described by the following extended Onsager relations~\cite{yuki12}:
\begin{eqnarray}
&&J_1=L_{11}X_1+L_{12}X_2,\label{eq.minimal.J1} \\
&&J_2=L_{21}X_1+L_{22}X_2-\gamma_{\rm c} {J_1}^2, \label{eq.minimal.J2}
\end{eqnarray}
where $L_{ij}$'s are the Onsager coefficients with reciprocity $L_{12}=L_{21}$.
The nonlinear term $-\gamma_{\rm c} J_1^2$ means the power dissipation into the cold
heat reservoir and $\gamma_{\rm c}>0$ denotes its strength. We also disregard
other possible nonlinear terms by assuming that their coefficients are too small
when compared with the power dissipation term.

In the limit of $X_i$'s $\to 0$,
the thermodynamic forces can be approximated as
$X_1\simeq W/T$ ($X_1=F/T$) and $X_2\simeq \Delta T/T^2$ for the cyclic (steady-state)
refrigerators, where $\Delta T \equiv T_{\rm h}-T_{\rm c}$ and $T\equiv (T_{\rm h}+T_{\rm c})/2$.
These expressions of the forces are consistent with~\cite{broeck05,yuki09}.
Then the extended Onsager relations in eqs.~(\ref{eq.minimal.J1}) and (\ref{eq.minimal.J2})
recover the usual linear Onsager relations such as in~\cite{broeck05,yuki09}.
Thus our minimally nonlinear irreversible model naturally includes the linear irreversible model
as a special case.
The non-negativity of the entropy production rate $\dot{\sigma}=J_1 X_1+J_2X_2$,
which is a quadratic form of $X_i$'s in the linear irreversible model, leads to the restriction of the
Onsager coefficients~\cite{yuki09,O,GM} to
\begin{eqnarray}
L_{11}\ge 0, L_{22}\ge 0, L_{11}L_{22}-L_{12}L_{21} \ge 0.\label{eq.onsager.coeffi.restrict}
\end{eqnarray}
Although our model includes the power dissipation as the nonlinear
term, we also assume eq.~(\ref{eq.onsager.coeffi.restrict}) holds.
This idea that the dissipation is included in the scope of the
linear irreversible thermodynamics is also found in~\cite{CW,AOGGLa,AOGGLb}.

The power input $P$ is given by $P\equiv \dot{W}=J_1X_1T_{\rm h}>0$.
Then the heat flux into the hot heat reservoir
$\dot{Q}_{\rm h}\equiv P+\dot{Q}_{\rm c}\equiv J_3$ is given as
\begin{eqnarray}
J_3=P+\dot{Q}_{\rm c}=J_1X_1T_{\rm h}+J_2.\label{eq.minimal.J3}
\end{eqnarray}
By solving eq.~(\ref{eq.minimal.J1}) with respect to $X_1$ and substituting it
into eqs.~(\ref{eq.minimal.J2}) and (\ref{eq.minimal.J3}),
we can rewrite $J_2$ and $J_3$ by using $J_1$ instead of $X_1$ as (see fig.~\ref{fig.1}):
\begin{eqnarray}
J_2=\frac{L_{21}}{L_{11}}J_1+L_{22}(1-q^2)X_2-\gamma_{\rm c} {J_1}^2,\label{eq.J2.by.J1}\\
J_3=\frac{L_{21}}{L_{11}}\frac{T_{\rm h}}{T_{\rm c}}J_1+L_{22}(1-q^2)X_2+\gamma_{\rm h} {J_1}^2,\label{eq.J3.by.J1}
\end{eqnarray}
where $q\equiv \frac{L_{12}}{\sqrt{L_{11}L_{22}}}$ is the usual
coupling strength parameter ($|q| \le 1 $ from
eq.~(\ref{eq.onsager.coeffi.restrict})) and $\gamma_{\rm h}$ denotes
the strength of the power dissipation into the hot heat reservoir as
\begin{eqnarray}
\gamma_{\rm h} \equiv \frac{T_{\rm h}}{L_{11}}-\gamma_{\rm c}>0.\label{eq.gh}
\end{eqnarray}

Eqs.~(\ref{eq.J2.by.J1}) and~(\ref{eq.J3.by.J1}) allow a clearer
description of the refrigerator than eqs.
(\ref{eq.minimal.J1}) and (\ref{eq.minimal.J2}) by considering $J_1$
as the control parameter instead of $X_1$ at constant $X_2$. So,
each term in eqs.~(\ref{eq.J2.by.J1}) and (\ref{eq.J3.by.J1}) has
explicit physical meanings (see also~\cite{yuki12}): the two first
terms mean the reversible heat transport between the working system
and the heat reservoirs; the second ones account for coupling
effects between the heat reservoirs; and the third ones account for the power dissipation ($\sim
J_1^2$) into the heat reservoirs, which inevitably occurs in a finite-time motion ($J_1 \ne 0$).

Using eqs.~(\ref{eq.J2.by.J1}) and (\ref{eq.J3.by.J1}) instead of eqs.~(\ref{eq.minimal.J1})
and (\ref{eq.minimal.J2}), the power input $P=J_3-J_2$ and the entropy production rate
$\dot{\sigma}=\frac{J_3}{T_{\rm h}}-\frac{J_2}{T_{\rm c}}$ are written as
\begin{eqnarray}
P=\frac{L_{21}}{L_{11}}\left(\frac{1}{\tau}-1\right) J_1+\frac{T_{\rm h}}{L_{11}}{J_1}^2,\label{eq.P}
\end{eqnarray}
and
\begin{eqnarray}
\dot{\sigma}=L_{22}(1-q^2)X_2^2+J_1^2\left(\frac{\gamma_{\rm c}}{T_{\rm c}}+\frac{\gamma_{\rm h}}{T_{\rm{h}}}\right).
\end{eqnarray}

For the refrigerator, the converter efficiency $z$ in eq.~(\ref{eq.chi.def}) is the COP
$\varepsilon$ defined as
\begin{eqnarray}
\varepsilon\equiv \frac{\dot{Q}_{\rm c}}{P}=\frac{J_2}{P}=\frac{\frac{L_{21}}{L_{11}}J_1+L_{22}(1-q^2)X_2-\gamma_{\rm c}J_1^2}{\frac{L_{21}}{L_{11}}(\frac{1}{\tau}-1)J_1+\frac{T_{\rm h}}{L_{11}}J_1^2},\label{eq.cop}
\end{eqnarray}
where we used eqs.~(\ref{eq.J2.by.J1}) and (\ref{eq.P}).
We also introduce the $\chi$-criterion as a product of the COP $\varepsilon$ and the cooling power $J_2$,
identifying the cooling power $Q_{\rm in}/t_{\rm cycle}$ in eq.~(\ref{eq.chi.def})
as $\dot{Q}_{\rm c}=J_2$:
\begin{eqnarray}
\chi=\frac{{J_2}^2}{P}=\frac{\left(\frac{L_{21}}{L_{11}}J_1+L_{22}(1-q^2)X_2-\gamma_{\rm
c}{J_1}^2\right)^2}
{\frac{L_{21}}{L_{11}}\left(\frac{1}{\tau}-1\right)J_1+\frac{T_{\rm h}}{L_{11}} {J_1}^2},\label{eq.chi}
\end{eqnarray}
where we used eqs.~(\ref{eq.J2.by.J1}) and (\ref{eq.P}).
Though eq.~(\ref{eq.chi.def}) using $Q_{\rm in}/t_{\rm cycle}$ appears to be valid only for cyclic heat
energy converters, eq.~(\ref{eq.chi}) using $J_2$ is valid even for steady-state refrigerators.
In order to get a positive cooling power $J_2>0$,
we see from eq.~(\ref{eq.J2.by.J1}) that $J_1$ should be located in a finite range as
\begin{eqnarray}
&&\frac{L_{21}-\sqrt{D}}{2L_{11}\gamma_{\rm c}}<J_1< \frac{L_{21}+\sqrt{D}}{2L_{11}\gamma_{\rm c}}\label{eq.J1_range_nontight},
\end{eqnarray}
where the discriminant $D$ is given by
\begin{eqnarray}
D\equiv L_{21}^2+4L_{11}^2L_{22}\gamma_{\rm c}(1-q^2)X_2>0.\label{eq.def.D}
\end{eqnarray}
We consider the optimization of the $\chi$-criterion with respect to $J_1$ in this range.
In the following, we first analyze the case $|q|=1$ called the tight-coupling condition.
Afterward we consider the nontight-coupling case $|q|\ne 1$.
\subsection{Tight-coupling case $|q|=1$}
In the tight-coupling case, eq.~(\ref{eq.J1_range_nontight}) becomes
\begin{eqnarray}
0< J_1 < \frac{L_{21}}{\gamma_{\rm c}L_{11}}.\label{eq.J1_range_tight}
\end{eqnarray}
We find $J_1=J_1^{\rm max\chi}$ which gives the maximum of $\chi$ by solving the equation $\partial \chi/\partial J_1=0$.
Noticing that $J_1=L_{21}/(\gamma_cL_{11})$
as the solution of $J_2=0$ is always one of the solutions of  $\partial \chi/\partial J_1=0$,
we can simplify that equation to the following quadratic one:
\begin{eqnarray}
J_1^2+\frac{3L_{21}}{2\varepsilon_{\rm C}T_{\rm h}}J_1-\frac{L_{21}^2}
{2\gamma_{\rm c}L_{11}\varepsilon_{\rm C}T_{\rm h}}=0.\label{eq.quadratic}
\end{eqnarray}
Solving this equation, we obtain the sole physically acceptable solution
which satisfies ${\partial}^2 \chi/\partial {J_1}^2<0$ and
locates in the range given by eq.~(\ref{eq.J1_range_tight}) as
\begin{eqnarray}
J_1^{{\rm max\chi}}=\frac{2{L_{21}}}{3\gamma_{\rm
c}{L_{11}}+\gamma_{\rm c}{L_{11}} \sqrt{8\varepsilon_{\rm
C}\left(1+\frac{\gamma_{\rm h}}{\gamma_{\rm c}}\right)+9}}
.\label{eq.J1.optimum}
\end{eqnarray}
Substituting eq.~(\ref{eq.J1.optimum}) into eq.~(\ref{eq.cop}),
we can obtain the general expression of the COP under maximum $\chi$-condition
for the model of minimally nonlinear irreversible refrigerators with the tight-coupling condition:
\begin{eqnarray}
\varepsilon_{{\rm max\chi}}=\frac{\sqrt{9+8\varepsilon_{\rm C}\left(1+\frac{\gamma_{\rm h}}{\gamma_{\rm c}}\right)}-3}{2\left(1+\frac{\gamma_{\rm h}}{\gamma_{\rm c}}\right)}.\label{eq.chi.optimum}
\end{eqnarray}
It is a monotonic decreasing function of the dissipation ratio $\gamma_{\rm h}/\gamma_{\rm c}$.
Considering that the assumption $\gamma_{\rm c}>0$ and $\gamma_{\rm c}<T_{\rm h}/L_{11}$
due to the non-negativity of $\gamma_{\rm h}$ as
in eq.~(\ref{eq.gh}),
the range of the dissipation ratio for the tight-coupling case is restricted to
\begin{eqnarray}
0<\frac{\gamma_{\rm h}}{\gamma_{\rm c}}<\infty.\label{eq.dissipation_ratio_range_q_1}
\end{eqnarray}
Then we can obtain the lower and the upper bounds by considering
asymmetric dissipation limits $\gamma_{\rm h}/\gamma_{\rm c} \to \infty$
and $\gamma_{\rm h}/\gamma_{\rm c} \to 0$, respectively, as
\begin{eqnarray}
\varepsilon_{\rm max\chi}^- \equiv 0 < \varepsilon_{{\rm max\chi}} < \frac{\sqrt{9+ 8\varepsilon_{\rm C}}-3}{2}\equiv \varepsilon_{\rm max\chi}^+,\label{eq.cop.bounds_2}
\end{eqnarray}
in full agreement with the results reported in~\cite{tu12d} and expressed in above eq.~(\ref{eq.cop.bounds}).
Here we stress that eq.~(\ref{eq.cop.bounds_2}) has been derived in a model-independent way and
is valid for general refrigerators, regardless of the types of refrigerator (steady-state or cyclic).
It is also interesting to see from eq.~(\ref{eq.cop.bounds_2}) that the lower bound for
$\varepsilon_{\rm max\chi}$
is $0$, which is in contrast to the case of the heat engines (see eq.~(21) in~\cite{yuki12}).
One may feel that this is a counterintuitive result at a first glance since the limit of the zero dissipation
$\gamma_{\rm c}\to 0$ in the cold heat reservoir seems to be an advantageous condition for a refrigerator.
This behavior can be understood as follows: from eq.~(\ref{eq.J1.optimum}),
$J_1^{\rm max\chi}$ diverges as $\sim O(\gamma_{\rm c}^{-1/2})$.
Then we can confirm $J_2^{\rm max\chi}\sim O(\gamma_{\rm c}^{-1/2})$ and $P_{\rm max\chi}\sim O(\gamma_{\rm c}^{-1})$ from
eqs.~(\ref{eq.J2.by.J1}) and (\ref{eq.P}), respectively,
which leads to $\varepsilon_{\rm max\chi}=J_2^{\rm max\chi}/P_{\rm max\chi}\sim O(\gamma_{\rm c}^{1/2})\to 0$.
This is an anomalous situation where the refrigerator transfers
the infinite amount of heat from the cold heat reservoir
to the hot one, consuming the infinite amount of work.
\subsection{Nontight-coupling case $|q|\ne 1$}
\begin{figure}
\begin{center}
\includegraphics[scale=1.2]{fig2.eps}
\end{center}
\caption{Numerical solutions for the nontight-coupling case $|q|\ne 1$: (a) Dependence of  $\varepsilon_{\rm max\chi}$
on the dissipation ratio $\gamma_{\rm h}/\gamma_{\rm c}$ under given temperatures ($T_{\rm h}=1$ and $T_{\rm c}=0.7$) for some $q$'s. We also plot eq.~(\ref{eq.chi.optimum}) for the tight-coupling case $|q|=1$ for comparison.
We always set $L_{11}=1$, $L_{22}=1$. Giving $q$, we can determine $L_{21}$
as $L_{21}=\sqrt{q^2L_{11}L_{22}}$.
(b) Dependence of the upper bound $\varepsilon_{\rm max\chi}^{q+}$ on  $\varepsilon_{\rm C}$
for some $q$'s. We also plot the upper bound $\varepsilon_{\rm max\chi}^+$ for the tight-coupling case
$|q|=1$ in eq.~(\ref{eq.cop.bounds_2}) for comparison.
We used the same values for $L_{ij}$'s as in (a).
We fixed $T_{\rm h}=1$ and changed $T_{\rm c}$.}\label{fig.2}
\end{figure}
Comparing with the tight-coupling case,
we found that the analytic solution of the nontight-coupling case
becomes too complicated to write down explicitly.
Thus we here solve the optimization problem numerically.
We first note that eq.~(\ref{eq.def.D}) implies the restriction
\begin{eqnarray}
\gamma_{\rm c}<-\frac{L_{12}^2}{4(1-q^2)L_{11}^2L_{22}X_2}\equiv \gamma_{\rm c}^+,\label{eq.gc_restrict}
\end{eqnarray}
for given Onsager coefficients $L_{ij}$'s and $X_2$.
Such restriction does not appear in the tight-coupling case.
Besides, we note that we also have the restriction $0<\gamma_{\rm c}<T_{\rm h}/L_{11}$
due to the assumption of the non-negativity of $\gamma_{\rm c}$ and $\gamma_{\rm h}$
as in eq.~(\ref{eq.gh}).
Thus, by combining these two inequalities,
we obtain the restriction of the dissipation ratio for the nontight-coupling case,
depending on the value of $\gamma_{\rm c}^+$ as
\begin{eqnarray}
0<\frac{\gamma_{\rm h}}{\gamma_{\rm c}}<\infty \label{eq.dissipation_ratio_range_q_ne_1_1}
\end{eqnarray}
in the case of ${T_{\rm h}}/{L_{11}}\le \gamma_{\rm c}^+ < \infty$, and
\begin{eqnarray}
\frac{\frac{T_{\rm h}}{L_{11}}-\gamma_{\rm c}^+}{\gamma_{\rm c}^+}<\frac{\gamma_{\rm h}}{\gamma_{\rm c}}<\infty
\label{eq.dissipation_ratio_range_q_ne_1_2}
\end{eqnarray}
in the case of $0<\gamma_{\rm c}^+<{T_{\rm h}}/{L_{11}}$, respectively.

We can obtain $\varepsilon_{\rm max\chi}$ numerically as follows:
we fix the values of $L_{11}$, $L_{22}$, $q$, $T_{\rm h}$ and $T_{\rm c}$.
Choosing $\gamma_{\rm c}$ so as to satisfy $0<\gamma_{\rm c}<T_{\rm h}/L_{11}$ and eq.~(\ref{eq.gc_restrict}),
which also determines $\gamma_{\rm h}$ from eq.~(\ref{eq.gh}),
we can find $\varepsilon$ at the maximum value of $\chi$ by calculating eq.~(\ref{eq.chi})
as $J_1$ is varied in the range of eq.~(\ref{eq.J1_range_nontight}).

We first consider typical dependence of $\varepsilon_{\rm max\chi}$
on the dissipation ratio $\gamma_{\rm h}/\gamma_{\rm c}$ under given temperatures for some $q$'s.
By changing $\gamma_{\rm c}$,
we plotted $\varepsilon_{\rm max\chi}$ as a function of $\gamma_{\rm h}/\gamma_{\rm c}$ in fig.~\ref{fig.2} (a).
As one of the interesting characteristics,
we can see that the dependence can become non-monotonic and $\varepsilon_{\rm max\chi}$
attains the maximum value at a finite dissipation ratio for each $q$.
This contrasts to
the tight-coupling case in eq.~(\ref{eq.cop.bounds_2}),
where its upper bound $\varepsilon_{\rm max\chi}^+$
is always attained in the asymmetric dissipation limit $\gamma_{\rm h}/\gamma_{\rm c}\to 0$.

Next we determine the upper bound $\varepsilon_{\rm max\chi}^{q+}$ and
the lower bound $\varepsilon_{\rm max\chi}^{q-}$ for the nontight-coupling case under various temperatures.
By changing $\varepsilon_{\rm C}$, we evaluated the maximum value of $\varepsilon_{\rm max\chi}$
as $\varepsilon_{\rm max\chi}^{q+}$ for some $q$'s  (see fig.~\ref{fig.2} (b)).
We found that $\varepsilon_{\rm max\chi}^{q+}$ is always lower than $\varepsilon_{\rm max\chi}^{+}$ in eq.~(\ref{eq.cop.bounds_2})
and monotonically decreases as $q$ is decreased for each $\varepsilon_{\rm C}$.
Meanwhile we also found that the lower bound $\varepsilon_{\rm max\chi}^{q-}$ agrees with
$0$ for each $q$ (data not shown), which can always be attained in the asymmetric dissipation limit $\gamma_{\rm h}/\gamma_{\rm c} \to \infty$
as in the tight-coupling case in eq.~(\ref{eq.cop.bounds_2}).
We also note that in the case where $\gamma_{\rm h}/\gamma_{\rm c}$ satisfies eq.~(\ref{eq.dissipation_ratio_range_q_ne_1_2}),
the lower bound $\varepsilon_{\rm max\chi}^{q-}=0$ can be realized also at the lower endpoint of
eq.~(\ref{eq.dissipation_ratio_range_q_ne_1_2}),
that is, in the limit of
$\gamma_{\rm h}/\gamma_{\rm c}\to (T_{\rm h}/L_{11}-\gamma_{\rm c}^+)/\gamma_{\rm c}^+$.
(See the lines of $q=0.3$, $0.5$ and $0.7$ in fig.~\ref{fig.2} (a).)
In this case,
$\varepsilon_{\rm max\chi}^{q-}$ becomes $0$ since $J_2^{\rm max\chi}$ vanishes:
in the lower endpoint of eq.~(\ref{eq.dissipation_ratio_range_q_ne_1_2}),
the discriminant $D$ in eq.~(\ref{eq.def.D}) becomes $0$. Thus the solution of the inequality
$J_2>0$ in eq.~(\ref{eq.J1_range_nontight}) approaches this multiple root in that limit, where $J_2$
in eq.~(\ref{eq.J2.by.J1}) is identically $0$.
This contrasts to the case of $\gamma_{\rm h}/\gamma_{\rm c}\to \infty$, where $J_2^{\rm max\chi}$ diverges.

\section{Example: leaky low-dissipation Carnot refrigerator}
As an example of the minimally nonlinear irreversible refrigerator,
we here introduce a leaky low-dissipation Carnot refrigerator model.
The original low-dissipation Carnot refrigerator without heat leak
studied in~\cite{carla12,tu12d} is an
extension of the quasistatic Carnot refrigerator by assuming that
heat transfer accompanying finite-time operation in each isothermal process is inversely proportional
to the duration of the process $t_i$:
\begin{eqnarray}
Q_{\rm h}=T_{\rm h} \Delta S+\frac{T_{\rm h} \Sigma_{\rm h}}{t_{\rm h}},\label{eq.low.qh}\\
Q_{\rm c}=T_{\rm c} \Delta S-\frac{T_{\rm c} \Sigma_{\rm c}}{t_{\rm c}}, \label{eq.low.qc}
\end{eqnarray}
where we denote by $\Delta S>0$ the quasistatic entropy change of the working system during
each isothermal process, and by $\Sigma_i/t_i$ the corresponding
entropy production with a constant strength $\Sigma_i>0$~\cite{carla12,tu12d}.
Then we formally modify this low-dissipation Carnot refrigerator
to a leaky model by adding a heat leak term as follows:
\begin{eqnarray}
Q_{\rm h}=T_{\rm h} \Delta S+\frac{T_{\rm h} \Sigma_{\rm h}}{t_{\rm h}}-\kappa (T_{\rm h}-T_{\rm c})(t_{\rm h}+t_{\rm c}),\label{eq.low.qh.leak}\\
Q_{\rm c}=T_{\rm c} \Delta S-\frac{T_{\rm c} \Sigma_{\rm c}}{t_{\rm c}}-\kappa (T_{\rm h}-T_{\rm c})(t_{\rm h}+t_{\rm c}), \label{eq.low.qc.leak}
\end{eqnarray}
where we assume the linear Fourier-type heat transfer from the hot heat reservoir to the cold one as the
heat leak per cycle, denoting by $\kappa >0$ the effective thermal conductivity.
We here reinterpret $Q_{\rm h}$ ($Q_{\rm c}$)
as the net heat transfer into (from) the heat reservoir per cycle,
instead of the heat transfer in each isothermal process,
considering the heat leak contribution which always continues during one-cycle.

We show that this leaky low-dissipation Carnot refrigerator is an example of the
minimally nonlinear irreversible refrigerator by rewriting eqs.~(\ref{eq.low.qh.leak})
and (\ref{eq.low.qc.leak}) into the forms of eqs.~(\ref{eq.J2.by.J1}) and (\ref{eq.J3.by.J1}):
since it works as a cyclic refrigerator and according to eqs.~(\ref{eq.def_J1_X1}) and (\ref{eq.def_J2_X2}),
we can define its thermodynamic fluxes and
forces as $J_1\equiv 1/t_{\rm cycle}=1/(t_{\rm h}+t_{\rm c})=1/((\alpha+1)t_{\rm h})$,
$J_2\equiv \dot{Q}_{\rm c}=Q_{\rm c}/t_{\rm cycle}$, $X_1\equiv W/T_{\rm h}$ and $X_2\equiv
1/T_{\rm h}-1/T_{\rm c}$, with $\alpha \equiv t_{\rm
c}/t_{\rm h}$. By using these definitions and eqs.~(\ref{eq.low.qh.leak})
and (\ref{eq.low.qc.leak}), we can calculate the Onsager coefficients
$L_{ij}$'s and the strength of the dissipation $\gamma_i$'s as
\begin{eqnarray}
\left(\begin{array}{cc}
L_{11} & L_{12} \\ L_{21} & L_{22}
\end{array}\right)
=
\left(\begin{array}{cc}
\frac{T_{\rm h}}{Y}
&\frac{T_{\rm h} T_{\rm c} \Delta S}{Y} \\
\frac{T_{\rm h} T_{\rm c}\Delta S}{Y}
& \frac{{T_{\rm c}^2{T_{\rm h}}{\Delta S}^2}}{Y}+T_{\rm h}T_{\rm c}\kappa
\end{array}\right)
,\label{eq.onsager.coeffi}
\end{eqnarray}
\begin{eqnarray}
\gamma_{\rm h}=T_{\rm h}\Sigma_{\rm h}(\alpha+1),\label{eq.minimal.gammah}\\
\gamma_{\rm c}=\frac{T_{\rm c}\Sigma_{\rm c}(\alpha+1)}{\alpha},\label{eq.minimal.gammac}
\end{eqnarray}
respectively, where the common factor $Y$ in eq.~(\ref{eq.onsager.coeffi})
is given as
\begin{eqnarray}
Y\equiv (\alpha+1)\left(T_{\rm h}\Sigma_{\rm h}+\frac{T_{\rm c}\Sigma_{\rm c}}{\alpha}\right).
\end{eqnarray}
From eq.~(\ref{eq.onsager.coeffi}) we see that
the Onsager reciprocity $L_{12}=L_{21}$ holds as expected.
In general, the Onsager coefficients show the nontight-coupling condition
$|q|=\frac{L_{12}}{\sqrt{L_{11}L_{22}}} < 1$, but recover the tight-coupling condition $|q|\to 1$
in the zero leak limit of $\kappa \to 0$.
Thus our model eqs.~(\ref{eq.J2.by.J1}) and (\ref{eq.J3.by.J1}) under the
tight-coupling condition includes the low-dissipation Carnot
refrigerator~\cite{carla12,tu12d} in eqs.~(\ref{eq.low.qh}) and (\ref{eq.low.qc}) as a special case.

Finally we note that we obtained eq.~(\ref{eq.chi.optimum}) for the tight-coupling case
maximizing $\chi$ only by one parameter $J_1$, whereas in the
original low-dissipation model without the heat leak
the optimization uses the two parameters $t_{\rm h}$
and $t_{\rm c}$~\cite{carla12,tu12d}. Thus the parameter space for the maximization is different.
In the formalism presented here, a further optimization of
$\chi=\chi(\alpha)$ by $\alpha$ is still possible
after the optimization by $J_1$ and it corresponds to the two-parameter maximization done
in~\cite{carla12,tu12d}.
But, that further optimization does not change the bounds in~eq.~(\ref{eq.cop.bounds}).
In particular, under symmetric conditions ($\Sigma_{\rm
h}=\Sigma_{\rm c}$) the further optimization by $\alpha$ gives $\alpha_{\rm
max\chi}=\sqrt{1-\tau}+1,\label{eq.alpha}$ from which we can
straightforwardly obtain the symmetric optimization COP in eq.~(3)~\cite{carla12}.
\section{Summary}
We introduced the extended Onsager relations as the
model of minimally nonlinear irreversible refrigerators, which can be applied to
both steady-state and cyclic refrigerators.
We analytically derived the general expression for the
coefficient of performance (COP) under maximum $\chi$-condition when the tight-coupling
condition is met and determined its upper and lower bounds which perfectly agree with the result derived previously in
the low-dissipation Carnot refrigerator model in~\cite{tu12d}.
We also studied the nontight-coupling case numerically and found that its upper bound is
always lower than that of the tight-coupling case, whereas all the cases have the same lower bound.
Moreover, we introduce the leaky low-dissipation Carnot refrigerator.
Calculating the Onsager coefficients and the strength of the power
dissipation explicitly, we proved that this leaky low-dissipation Carnot refrigerator is an example of the minimally nonlinear irreversible refrigerator.

\acknowledgments
Y. Izumida acknowledges the financial support from a
Grant-in-Aid for JSPS Fellows (Grant No. 22-2109).
JMM Roco and A. Calvo Hern\'andez thank financial support from
Ministerio de Educaci\'on y Ciencia of Spain under Grant  No.
FIS2010-17147 FEDER.


\section{Introduction}

The Gonihedric Ising model was defined by Savvidy and Wegner \cite{1} as a cubic lattice
transcription of an earlier Gonihedric random surface model developed by Savvidy
\cite{2}  which was intended to be a discretization of string theory. In the case of the Gonihedric Ising model, the spin cluster boundaries were used to model
a gas of surfaces and the spin couplings were chosen so that edges and intersections of plaquettes rather than their surface areas were weighted, which was also the defining feature of the Gonihedric random surface model.

Such a weighting can be arranged using a one-parameter family of generalized $3D$ Ising models with fine-tuned nearest neighbour
$\langle ij \rangle$, next to nearest neighbour $\langle \langle ij \rangle \rangle$ and plaquette
interactions $[ijkl]$.  The Hamiltonians in this family are given by
\begin{equation}
\label{e1}
H_{gonihedric} = - 4 \kappa \sum_{\langle ij\rangle }\sigma_{i} \sigma_{j}  +
\kappa \sum_{\langle \langle ij\rangle \rangle }\sigma_{i} \sigma_{j}
- ( 1-\kappa )\sum_{[ijkl]}\sigma_{i} \sigma_{j}\sigma_{k} \sigma_{l} \; ,
\end{equation}
where the parameter $\kappa$ gives the relative weight of two plaquettes meeting at right angles and four-plaquette intersections in the spin cluster boundaries.
One member of this family, the purely plaquette action with $\kappa=0$
\begin{equation}
\label{e2k}
H_{\kappa=0} =  -  \sum_{[ijkl]}\sigma_{i} \sigma_{j}\sigma_{k} \sigma_{l}
\end{equation}
has merited particular attention due to its interesting properties, both static and dynamical \cite{3a,3}. It has a
strong first order transition surrounded by a region of metastability and shows evidence of glass-like behaviour at lower temperatures in spite of having no quenched disorder.

The symmetry properties of both the ground states and low temperature phase of $H_{\kappa=0}$ are  unusual. The standard Ising model with only nearest neighbour interactions
\begin{equation}
\label{e0I}
H_{Ising} =  -  \sum_{\langle ij \rangle}\sigma_{i} \sigma_{j}
\end{equation}
displays a twofold symmetry in its ground state, whereas a $3D$
$\mathbb{Z}_2$ gauge theory (where the spins live on plaquette edges rather than vertices as in $H_{\kappa=0}$)
\begin{equation}
\label{e0G}
H_{Gauge} =  -  \sum_{[ijkl]} U_{ij} U_{jk} U_{kl} U_{li}
\end{equation}
has the local gauge symmetry $U_{ij} \rightarrow s_i U_{ij} s_j $, with $s_i = \pm 1$. The symmetries of plaquette Gonihedric model
of equ.\,(\ref{e2k}) lie somewhere between these since it is possible to flip arbitrary, possibly intersecting, planes of spins
in the ground state at zero energy cost. Both Monte-Carlo simulations and low temperature expansions \cite{4} confirm that this remains a symmetry of the low temperature phase, which is thus highly degenerate.  A consequence of this symmetry is that the standard magnetization is zero in the low temperature phase. The nature of the magnetic order has thus remained unclear.

In this paper we  discuss the recent suggestion by Hashizume and Suzuki \cite{5} that
a suitable order parameter for $H_{\kappa=0}$ might be constructed using nearest-neighbour spin-spin correlations $\langle \sigma_i \sigma_j \rangle$, which was prompted
by Suzuki's earlier study of an {\it anisotropic} Gonihedric model \cite{6}. We  conduct measurements of this candidate order parameter on
a toy ensemble of ground state configurations constructed by flipping planes of spins starting from a
ferromagnetic starting configuration to mimic the ground state structure of the Gonihedric model, before moving on to  simulations of both the anisotropic and isotropic Gonihedric models themselves.

\section{The Anisotropic Gonihedric Model (a.k.a. the Fuki-Nuke Model)}

The study of a strongly anisotropic variant of the plaquette Gonihedric
Hamiltonian by Suzuki \cite{6} pre-dates the work of Savvidy and Wegner by some twenty years.
For reasons that will become apparent, Suzuki dubbed this the ``Fuki-Nuke'' or ``no roof'' model and showed that it was (non-trivially) equivalent to a stack of $2D$ Ising models. This was later rediscovered by Jonsson and Savvidy \cite{7} using a transfer matrix approach and discussed by Castelnovo
{\it et.al.} \cite{8} using methods more akin to Suzuki's original approach.

In Suzuki's  work the principal concern was to obtain solvable $3D$ spin models, and the Fuki-Nuke
model was constructed as a specific example. The starting point was an anisotropically coupled variant of
$H_{\kappa=0}$
\begin{eqnarray}
H_{aniso} = - J_x  \sum_{[ijkl]}\sigma_{i} \sigma_{j}\sigma_{k} \sigma_{l}  - J_y \sum_{[ijkl]}\sigma_{i} \sigma_{j}\sigma_{k} \sigma_{l}
 - J_z \sum_{[ijkl]}\sigma_{i} \sigma_{j}\sigma_{k} \sigma_{l}
\end{eqnarray}
where $J_x, J_y, J_z$ were the couplings for plaquettes perpendicular to the $x,y,z$ axes respectively.
Setting $J_x=J_y=1$ and $J_z=0$ gave a model where the horizontal ``roof'' plaquettes had zero coupling, hence the sobriquet of ``Fuki-Nuke''.

To show that this is equivalent to a stack of $2D$ Ising models requires a little further work, which can be done succinctly by using the non-linear $\sigma$-$\tau$ transformation, following Suzuki \cite{6} and  Castelnovo {\it et.al.} \cite{8}. We first define bond spin variables $\tau_{j} = \sigma_i \sigma_j$ that can be thought of as living at the centre of (only) the vertical lattice bonds.
Since we have condensed the site labels on both sets of spins, we  write $\sigma_k$ and $\tau_k$ to represent any of the spins in the $kth$ horizontal $2D$ layer. With this notation the direct transformation for a spin in the $kth$  layer  is
\begin{equation}
\sigma_k = \tau_1 \tau_2 \tau_3 \cdots \tau_k
\end{equation}
and the inverse transformation is given by
\begin{equation}
\tau_1 = \sigma_1, \;  \tau_2 = \sigma_1 \sigma_2 \; , \; \ldots \; , \; \tau_N = \sigma_{N-1} \sigma_N
\end{equation}
In order to obtain a one-to-one correspondence between
$\sigma$ and $\tau$ spin configurations we need to specify the value of the $\sigma_1 \, / \, \tau_1$ spins on a given (first in this case) horizontal plane,
but in all the other planes a product of two $\sigma$ spins at the end of each
vertical bond gives the $\tau$ spin.
The resulting Hamiltonian is then
\begin{equation}
H = - \sum_i \left( \tau_{i}  \tau_{i + \hat e_x} + \tau_{i} \tau_{i + \hat e_y} \right)
\end{equation}
where $\hat e_x, \hat e_y$ are unit vectors in the horizontal directions.
We have not written the contribution of the fixed plane, which will vanish in the thermodynamic limit. The $\tau_{i}$ spins are coupled to their horizontal nearest neighbours in the ${\hat e_x}$ and ${\hat e_y}$ directions
and there are  no vertical inter-plane couplings.
We have thus shown that the Fuki-Nuke model is equivalent to a stack of
uncoupled $2D$ Ising models arranged horizontally one above each other.

The magnetic order parameter for the $i$th $2D$ Ising layer is
given by the standard expression
\begin{equation}
\label{Mone}
M_{2D, \, i} =  \left< \frac{1}{N^2} \sum_{single \; plane} \tau_{i} \right>
\end{equation}
which may be rewritten in terms of the original $\sigma_i$ spins as
\begin{equation}
M_{2D, \, i} = \left<  \frac{1}{N^2} \sum_{single \; plane}   \sigma_i \sigma_{i +  \hat e_z}   \right> \; .
\end{equation}
where $\hat e_z$ is  the unit vector in the vertical direction. The nearest neighbour $ \sigma_i \sigma_{i +  \hat e_z}$ correlator is measured on all the vertical bonds bisected by the horizontal plane of the corresponding $\tau_i$ spins.
The unusual appearance of a two-spin correlator as an order parameter is simply the result
of expressing the magnetization in terms of the original $\sigma_i$ variables rather than the bond spins $\tau_{i}$. Remembering that
$\tau_i =  \sigma_i \sigma_{i +  \hat e_z}$ and using the standard
results for the $2D$ Ising model magnetization   we have \cite{6}
\begin{equation}
\label{M2D}
M_{2D, \, i} = \left<  \frac{1}{N^2} \sum_{single \; plane} \sigma_i \, \sigma_{i +  \hat e_z}  \right>  = \pm \left( 1 - \sinh^{-4}( \beta)\right)^{1/8}
\end{equation}
which will behave as $\pm \left| \beta - \beta_c \right|^{1 \over 8}$
near the critical point $\beta_c = \ln ( 1 + \sqrt{2})$.
More generally, since
$\sigma_{k} = \tau_1 \, \tau_2 \, \ldots \, \tau_{k}$ and the different $\tau_i$ layers are decoupled,
\begin{equation}
M_{2D, \, i, \, n} = \left<\frac{1}{N^2} \sum_{single \; plane} \sigma_i \, \sigma_{i +  n \hat e_z} \right> =   ( M_{2D, \, i} )^n \; .
\end{equation}
The in-plane horizontal correlations between ``dimers'' of $\sigma_i$ spins,
where the $i,j$ spins now lie on the same plane,
\begin{equation}
\left< ( \sigma_i \, \sigma_{i +  \hat e_z} ) (  \sigma_{j} \, \sigma_{j+  \hat e_z}) \right> = \left< \tau_{i} \, \tau_{j} \right>
\end{equation}
are also, of course, Ising-like.

Since the various  $2D$ Ising layers magnetize independently  we could construct a quasi-$3D$ order parameter to be
of the form
\begin{equation}
\label{Mabs}
M_{abs}  = \frac{1}{N} \sum_{xy \, planes} \Big\langle  \Big| \frac{1}{N^2} \sum_{single \; plane}   \sigma_i \sigma_{i +  \hat e_z}  \Big| \Big\rangle
\end{equation}
or
\begin{equation}
\label{Msq}
M_{sq} =  \frac{1}{N} \sum_{xy \, planes} \Big\langle  \left( \frac{1}{N^2}  \sum_{single \; plane}   \sigma_i \sigma_{i +  \hat e_z}   \right)^2 \Big\rangle \;
\end{equation}
to avoid inter-plane cancellations.
We have explicitly retained the various normalizing factors in equs.(\ref{Mabs}, \ref{Msq}) for a cubic  lattice with $N^3$ sites.
We could also define similar quantities for non-nearest neighbour correlators
\begin{equation}
\label{MInabsxyz}
M_{abs,\, n}^{} =  \frac{1}{N} \sum_{xy \, planes} \Big\langle  \Big| \frac{1}{N^2} \sum_{single \; plane}   \sigma_i \sigma_{i + n \hat e_{z} } \Big| \Big\rangle
\end{equation}
or
\begin{equation}
\label{MInsqxyz}
M^{}_{sq, \, n} =  \frac{1}{N} \sum_{xy  \, planes} \Big\langle  \left( \frac{1}{N^2}  \sum_{single \; plane}
 \sigma_i \sigma_{i +  n \hat e_{z}}   \right)^2 \Big\rangle
\end{equation}
with $n>1$. Since each plane magnetizes at the same point the resulting expressions are given by
\begin{eqnarray}
\label{Msinh}
M_{abs} &=& \left( 1 - \sinh^{-4}( \beta)\right)^{1/8}  \nonumber \\
M_{sq} &=&  \left( 1 - \sinh^{-4}( \beta)\right)^{1/4}  \nonumber \\
M_{abs, \, n} &=& \left( 1 - \sinh^{-4}( \beta)\right)^{n/8}  \nonumber \\
M_{sq, \, n} &=&  \left( 1 - \sinh^{-4}( \beta)\right)^{n/4}  \nonumber \; .
\end{eqnarray}
As the various other variants of the magnetization are given in terms of $M_{abs} = | M_{2D, \, i} |$ this is probably the most
suitable  choice of ``three-dimensional'' order parameter for the Fuki-Nuke model.
In addition, it emphasizes that  the critical behaviour of the Fuki-Nuke model is really that of the $2D$ Ising model in each layer.

In the anisotropic model the orientation of the planes in the summation is fixed by the choice of zero coupling, which we have taken to be horizontal planes in our case.  The suggestion by Hashizume and Suzuki  in \cite{5}  is that a similar order parameter measuring coplanar, nearest neighbour spin correlations could still serve for the {\it isotropic} Gonihedric model.
They considered a mean-field approach to the model and correlation function inequalities in order to corroborate this and found supporting evidence from both, as well as extending the discussion to disordered couplings.
Here, we
investigate the order parameter in the isotropic $3D$ Gonihedric model  directly using Monte-Carlo simulations, restricting discussion to the purely ferromagnetic case for simplicity.

In equs.\,(\ref{Mabs},\ref{Msq}) we have introduced either a modulus or square for the spin correlators in order to avoid inter-plane cancellations, whereas Suzuki's approach is tantamount to considering the direct  equivalent of equ.\,(\ref{Mone})
\begin{equation}
\label{MSuz}
M =  \frac{1}{N} \sum_{planes} \Big\langle   \frac{1}{N^2} \sum_{single \; plane}   \sigma_i \sigma_{i +  \hat e_z}   \Big\rangle
\end{equation}
in the presence of an external field (later taken to zero) which picks out one of the $2D$ magnetizations on each plane.
The possibility of using an order parameter akin to the $M_{sq}$ in equ.\,(\ref{Msq}) has been  mentioned previously by Lipowski \cite{3a}, who noted in a simulation of $H_{\kappa=0}$ that it  appeared to possess the correct behaviour. He also remarked on the possibility of using non-nearest neighbour correlators as the order parameter(s), which we shall discuss further below.
In the remainder of the paper we compare measurements of the candidate order parameters on toy ground state configurations, the Fuki-Nuke model and the isotropic $3D$ Gonihedric model itself.

\section{Monte Carlo investigations: Toy Model, Fuki-Nuke and the Real Thing}

\subsection{Toy Model - Ground States}
Ground states of the  Gonihedric Hamiltonian $H_{\kappa=0}$  may be obtained by flipping some set of planes of spins in a reference ferromagnetic configuration which has all spins up or down. A simple test for any candidate order parameter is then to generate an ensemble of such configurations and measure its expectation value. A typical configuration on a finite cubic lattice
generated in this manner would look something like that in Fig.\,(1).
\begin{figure}[h]
\begin{center}
\includegraphics[height=4cm]{cube_stripe3}
\label{cube_stripe}
\caption{A ``ground-state'' configuration generated by flipping arbitrary planes of spins, shown dotted,
starting in a ferromagnetic state with all the spins positive or negative.}
\end{center}
\end{figure}
We can also measure proposed order parameters on purely random (i.e. infinite temperature) configurations. Any that fail to distinguish between the
 ground states and random configurations would not, at first sight, appear to be suitable  for use in a  simulation of the full Hamiltonian.

If the suggestion of \cite{5} is correct, we would expect that the order in the
low temperature phase of the
isotropic $3D$ Gonihedric model could be discerned using an appropriate generalization of either the  $M_{abs}$ and $M_{sq}$ in equs.\,(\ref{Mabs},\ref{Msq}). For   isotropic  configurations such as that in Fig.\,(1) the system should be agnostic to the orientation
of the planes used in calculating the correlators, so any of the three choices of orientation for reference planes in
\begin{eqnarray}
\label{Mabsxyz}
M^{x,y,z}_{abs} &=&  \frac{1}{N} \sum_{yz|xz|xy \; planes} \Big\langle  \Big| \frac{1}{N^2} \sum_{single \; plane}   \sigma_i \sigma_{i +  \hat e_{x,y,z}} \Big| \Big\rangle
\end{eqnarray}
or
\begin{eqnarray}
\label{Msqxyz}
M^{x,y,z}_{sq} &=&  \frac{1}{N} \sum_{yz|xz|xy  \, planes} \Big\langle  \left( \frac{1}{N^2}  \sum_{single \; plane}
 \sigma_i \sigma_{i +  \hat e_{x,y,z}}   \right)^2 \Big\rangle
\end{eqnarray}
where $\hat e_x, \hat e_y, \hat e_z$ are the unit vectors in the $x,y,z$ directions respectively,
should be able to serve as order parameters. In all cases the correlators are nearest neighbour correlators for the $\sigma_i$ spins on bonds perpendicularly
bisected by the reference planes, as in the Fuki-Nuke model.
Reassuringly, we find that the various  $M^{x,y,z}_{abs}$ from equ.\,(\ref{Mabsxyz}) and the $M^{x,y,z}_{sq}$ from equ.\,(\ref{Msqxyz}) all give $+1$
when measured on an ensemble of flipped states such as that in Fig.\,(1) so they do, indeed, pick up the planar order
that is characteristic of the Gonihedric models. In addition, the $M^{x,y,z}_{abs}$  and $M^{x,y,z}_{sq}$ all give zero when measured on purely random configurations, so they are capable of distinguishing this order from the complete disorder at infinite temperature.

The choice of nearest neighbour correlators in defining the order parameter does not  appear to be obligatory, though given the line of reasoning leading to it from the Fuki-Nuke model it is probably the most natural. If we measure quantities such as
\begin{equation}
\label{Mnabsxyz}
M_{abs,\, n}^{x,y,z} =  \frac{1}{N} \sum_{yz|xz|xy \, planes} \Big\langle  \Big| \frac{1}{N^2} \sum_{single \; plane}   \sigma_i \sigma_{i + n \hat e_{x,y,z}}  \Big| \Big\rangle
\end{equation}
or
\begin{equation}
\label{Mnsqxyz}
M^{x,y,z}_{sq, \, n} =  \frac{1}{N} \sum_{yz|xz|xy  \, planes} \Big\langle  \left( \frac{1}{N^2}  \sum_{single \; plane}
 \sigma_i \sigma_{i +  n \hat e_{x,y,z}}   \right)^2 \Big\rangle
\end{equation}
with $n>1$ on the ensemble of flipped ground state configurations these also give $+1$ (and zero on random configurations). That this should be so is clear from starting with configurations in which planes of one orientation only, say horizontal, are flipped. For these, plane-to-plane correlations at any separation must be $\pm 1$. Taking a modulus or square for each planar sum then makes all of these $+1$. Flipping planes of spins in the two remaining directions perpendicular to this will change the signs of lines of spins
on both the horizontal planes contributing to the correlator sums and therefore not change the values of the contributing correlators.

When the standard magnetization is zero
\begin{equation}
M = \frac{1}{N^3} \sum_i \langle \sigma_i \rangle =0
\end{equation}
the magnetic susceptibility $\chi$ is given by
\begin{equation}
\label{chi}
\chi = \frac{1}{N^3} \sum_i \sum_j \langle \sigma_i \sigma_j \rangle = \frac{1}{N^3} \sum_i \sum_n \langle \sigma_i \sigma_{i+n} \rangle
\end{equation}
which is  related to the class of order parameters discussed here, since the sums may be arranged in a similar manner to those in $M_{sq}$ and $M_{abs}$. However, the expression for $\chi$ includes a summation over distances greater than one and also retains the signs of the various different contributions in the sum.
Consequently, the susceptibility does not distinguish between configurations such as those in Fig.\,(1) and random configurations, giving $+1$ in both cases, but as we shall see below  it  {\it does} appear to display the characteristics of an order parameter
in a simulation of the full isotropic Gonihedric model.

\subsection{Fuki-Nuke Model}
Since the anisotropic Gonihedric Hamiltonian, or Fuki-Nuke model, originally inspired the definition of the candidate order parameters it too can  serve as a test case for simulations. We know  that $2D$ Ising criticality is expected in each horizontal layer, so the transition point will be seen in the continuum limit at the $2D$ Ising value of $\beta_c = \ln ( 1 + \sqrt{2} ) \simeq 0.88$.
In the Fuki-Nuke model the orientation of the the $2D$ planes used to measure the nearest neighbour spin correlations is fixed by the choice of the zero couplings. In our simulations we take $J_z=0$, which gives horizontal planes in the sums
and means that $M^z_{abs}$ should be the correct choice for the order parameter.

We simulate the Fuki-Nuke model on a $20^3$ lattice with periodic boundary conditions using $10^7$ Metropolis Monte-Carlo measurement sweeps with a hot start and an equilibration time of $10^5$ sweeps as a test case.
In Fig.\,(2) we can see that both $M_{abs}$ and $M_{sq}$ give a clear signal for the pseudo-critical point, which can be observed close to the continuum value of
$\beta \sim 0.88$ already on  the $20^3$ lattice via, for example, the peak in the specific heat measurements.  As a consequence of cancellations between differently signed $2D$ layers of Ising spins in the magnetized phase
the standard magnetization $M$ remains zero for all $\beta$ in the Fuki-Nuke model
and does not give an obvious  signal for the transition.
\begin{figure}[h]
\begin{center}
\includegraphics[height=6cm]{aniso20}
\label{aniso20}
\caption{The energy and order parameters $M_{abs}$ and $M_{sq}$ on a $20^3$  lattice for the Fuki-Nuke model.}
\end{center}
\end{figure}

\subsection{(Isotropic) Gonihedric Model}

The preceding measurements reported for the toy ground state model and the Fuki-Nuke model suggest that
any of $M^{x,y,z}_{abs}$ from equ.\,(\ref{Mabsxyz}), the  $M^{x,y,z}_{sq}$ from equ.\,(\ref{Msqxyz}) or even the $M^{x,y,z}_{abs, \, n}$ and
$M^{x,y,z}_{sq, \, n}$ of equs.\,(\ref{Mnabsxyz},\ref{Mnsqxyz}) might be viable candidates for an order parameter in the isotropic $3D$ Gonihedric model itself.
To clarify this, we  carried out simulations
of the  isotropic Gonihedric model on $10^3$ and $15^3$ lattices using
$10^7$ Metropolis Monte-Carlo measurement sweeps, following a hot start and an equilibration time of $10^5$ sweeps. The various $M^{x,y,z}_{abs}$ and $M^{x,y,z}_{sq}$  along with
$M^{x,y,z}_{abs, \, n}$ and $M^{x,y,z}_{sq, \, n}$ for $n>1$ were measured on both lattices.

In Fig.\,(3) we plot the measurements of $M_{sq}^{x,y,z}$ and $M_{abs}^{x,y,z}$
on the $10^3$ lattice and in Fig.\,(4) on the $15^3$ lattice. The measurements of the $M$'s for the different $x,y,z$ orientations of the reference planes are indistinguishable from each other and the observed jump in the various $M$'s coincides with the pseudo-critical values of $\beta$ determined from the jump in the energy observed
in simulations ($\sim 0.54$ on the $10^3$ lattice and $\sim 0.57$ on the $15^3$ lattice with ``hot'' starts).  These results confirm that the orientation of the reference planes is, indeed, irrelevant in measuring these magnetizations and that all three choices are equivalent. They also show clearly that the
correlators distinguish between the disordered and ordered phases of the isotropic model as suggested by the arguments in \cite{5}.
\begin{figure}[h]
\begin{center}
\includegraphics[height=6cm]{M_10}
\label{M_10}
\caption{$M_{abs}^{x,y,z}$ (upper) and $M_{sq}^{x,y,z}$ (lower) on a $10^3$ lattice for the isotropic Gonihedric model. Lines are drawn to guide the eye.}
\end{center}
\end{figure}
\begin{figure}[h]
\begin{center}
\includegraphics[height=6cm]{M_15}
\label{Msq_20}
\caption{$M_{abs}^{x,y,z}$ (upper) and $M_{sq}^{x,y,z}$ (lower) on a $15^3$ lattice for the isotropic Gonihedric model. Again, lines are drawn to guide the eye.}
\end{center}
\end{figure}

To guard against the possibility of discarding too much information by considering a simple global order parameter we can measure the distribution  $P(M_{2D, \,  i})$ of the coplanar $2D$ spin correlations. In Figs.\,(5,6)
we have combined the three orientations of reference planes to plot a $P(M_{2D})$  in the high temperature and low temperature phases respectively. The three individual orientations contributing to this give identical distributions. In Fig.\,(5) $P(M_{2D})$ is  gaussian (narrowing with system size and increasing temperature) and centred
at $M_{2D}=0$, whereas in the low temperature phase in Fig.\,(6) at $\beta=0.58$ there are two sharp peaks at $M_{2D} = \pm 1$. We would expect the slight asymmetry between the peaks in the plot to disappear with better statistics and increasing lattice size.
The observed behaviour strongly supports the hypothesis of Fuki-Nuke type order, where coplanar correlators ``magnetize'' independently of the other planes, in the isotropic model too.
The plots also demonstrate that a $[0,1]$ type order parameter is sufficient for the isotropic Gonihedric model, since there is no non-trivial structure in $P(M_{2D})$ such as that seen in the overlap distribution $P(q)$ of (mean-field) spin glasses.
\begin{figure}[h]
\begin{center}
\includegraphics[height=5cm]{Hist20}
\label{H_10}
\caption{$P(M_{2D})$ at $\beta=0.20$ in the high temperature phase on a $10^3$ lattice}
\end{center}
\end{figure}
\begin{figure}[h]
\begin{center}
\includegraphics[height=5cm]{Hist58}
\label{H_58}
\caption{$P(M_{2D})$ at $\beta=0.58$ in the low temperature phase on a $10^3$ lattice}
\end{center}
\end{figure}

We also find, as with the toy ground state ensemble, that $M_{abs, \, n}^{x,y,z}$ and
$M_{sq, \, n}^{x,y,z}$ for $n>1$ show very similar behaviour to $M_{abs}^{x,y,z}$ and
$M_{sq}^{x,y,z}$.
In Fig.\,(7)  $M_{abs, \, 6}^z$ on a
$15^3$ lattice is plotted, along with $M_{abs}^z$ for comparison. $M_{abs, \, 6}^z$ is typical of the various $M_{abs, \, n}$ (and $M_{sq, n}$) for other $n$.
\begin{figure}[h]
\begin{center}
\includegraphics[height=6cm]{Mabs_15_06}
\label{Mabs_15_06}
\caption{$M_{abs}^{z}$ (upper) and $M_{abs, 6}^{z}$ (lower) on a $15^3$ lattice for the isotropic Gonihedric model.}
\end{center}
\end{figure}
As might be expected the magnetization constructed from the more widely spaced correlators in $M_{abs, 6}^{z}$ is lower in the disordered phase, but it too jumps sharply to one at the transition point just as does $M_{abs}^{z}$.  Since we employ periodic boundary conditions the measured magnetizations repeat from $M_{abs, 7}$ on the $15^3$ lattice as  the distances in the correlators wrap around the lattice. With the isotropic model there is no longer a mapping to the $2D$ Ising model, unlike the Fuki-Nuke model, so we do not have an exact expression for the planar average of nearest neighbour correlators as in equ.\,(\ref{M2D}), nor do we have any assurance that different (bi-)layers of spins decouple, which allows us to express everything in terms of the nearest neighbour correlators alone in the Fuki-Nuke case. Nonetheless, it would appear that the most fundamental correlations in the isotropic case are also those between nearest neighbour spins. The observed signal for the transition in the higher correlators are a consequence of that in the nearest neighbour correlations, just as with the toy ground states and the Fuki-Nuke model.

A curious feature of the isotropic Gonihedric model simulations that should be highlighted is that
the standard magnetic susceptibility, as defined in equ.\,(\ref{chi}):
$\chi = \frac{1}{N^3} \sum_i \sum_j \langle \sigma_i \sigma_j \rangle = \frac{1}{N^3} \sum_i \sum_n \langle \sigma_i \sigma_{i+n} \rangle$,
also behaves as an effective order parameter, unlike the toy ground state ensemble. This can be seen in Fig.\,(8) where $\chi$ is  $+1$  at high temperatures and drops sharply to zero
at the observed pseudo-critical points.
\begin{figure}[h]
\begin{center}
\includegraphics[height=6cm]{Susp}
\label{Susp_10_20}
\caption{The magnetic susceptibility $\chi$ on both the $10^3$ and  $15^3$ lattices is  one at high temperatures and drops sharply to zero
at the observed pseudo-critical points.}
\end{center}
\end{figure}
We have already seen in our discussion of the toy model for ground states that $\chi$ does {\it not} distinguish between random, infinite temperature configurations and the
flipped ground-state configurations, which means that the observed behaviour is not  a direct consequence of $\chi$ detecting
the coplanar order of the low temperature phase as with $M^{x,y,z}_{abs,sq}$. It appears rather to be a dynamical  ergodicity breaking  effect. At high temperatures the flip symmetry of the Hamiltonian is not broken, so $\chi$ must be one, whereas the observed value of zero for $\chi$ throughout the  low temperature (high $\beta$) phase is due to the  system ``freezing'' into one configuration as the Metropolis spin flip acceptance drops rapidly to zero at $\beta_c$. It would be interesting to explore the behaviour of $\chi$ in a multicanonical simulation such as that
carried out in \cite{Marco}. In this case the system would still be free to explore multiple ground states, so the behaviour of $\chi$ would presumably be more akin to that in the toy ground state ensemble (and less like an order parameter).

It is also worth noting that the strong hysteresis associated with the first order transition means that  a more accurate estimate of the continuum transition
point $\beta_c (\infty)$ is best done with such a multi-canonical simulation, or similar. In \cite{Marco} the value of  $\beta_c (\infty) = 0.549 25(6)$ obtained from the finite size scaling of the specific heat maxima is a little higher than that obtained
 from taking careful account of scaling corrections arising from  the boundary conditions,  $\beta_c (\infty) = 0.547 57(63)$,
using a Metropolis simulation \cite{Fixed}. The various pseudo-critical temperatures observed on differently sized lattices in a multi-canonical simulation fall between those observed for hot and cold starts in Metropolis simulations such as those performed here on the same lattices, as one might expect.

\section{Conclusions}

As Suzuki remarked many years ago \cite{6} the anisotropic Gonihedric Ising (Fuki-Nuke) model may be
reformulated in terms of new spin variables as a stack of decoupled $2D$ Ising models which magnetize independently.
When translated back into the original spin variables the magnetization of each $2D$ Ising model is expressed in terms of a
two spin nearest neighbour correlator $\langle \sigma_i \sigma_j \rangle$.
To take account of the independent ordering of planes a  quasi-$3D$ order parameter may then be constructed by taking  the modulus, $M_{abs}$
of nearest neighbour correlators summed over  each plane.

We investigated numerically
the suitability of using this and related quantities as an order parameter for the  isotropic $3D$ Gonihedric Hamiltonian, $H_{\kappa=0}$, following the suggestion of \cite{5}. The transition in this case is now first order
and hence no longer in the universality class of of the $2D$ Ising model, but the distinctive planar flip symmetry remains.
We concluded that  the $M^{x,y,z}_{abs}$  of equ.\,(\ref{Mabsxyz}) might serve as an order parameter for the isotropic model and found that the  choice of reference direction $x,y,z$ was immaterial in this case. The absence of non-trivial structure in the distribution $P(M_{2D})$ of the coplanar magnetizations suggested that such a global parameter was sufficient to capture the order.

Somewhat unexpectedly, the standard magnetic susceptibility $\chi$  also shows the characteristics of an order parameter in the simulations of the isotropic Gonihedric model since it was measured to be
one at  high temperatures and dropped sharply to zero at the pseudo-critical point. However, $\chi$  failed to distinguish between an ensemble of flipped toy ground states and purely random configurations. We suggested that the behaviour in the full Gonihedric simulations was a dynamical consequence of freezing and that a multicanonical simulation might help to clarify this.

Regarding similar order parameters in other models, the dual Gonihedric model may be written as strongly anisotropic Ashkin-Teller model and also possesses similar flip symmetries to the purely plaquette model  \cite{99}. Understanding the low temperature order in this dual model, whose standard magnetization and polarization are zero at all $\beta$, might proceed along similar lines to those considered here.

\section{Acknowledgements}
D. A. Johnston would like to thank Martin Weigel, Wolfhard Janke and Adam Lipowski for helpful discussions.

\bigskip
\bigskip
\bigskip