\section{Introduction}
Duality \cite{krammers_wannier,wannier,onsager,syozi} is an important symmetry of
Ising models,
as it implies a relationship between the partition functions of two Ising models with Hamiltonians
 $\mathscr{H}_{1}\left(J\right)$ and $\mathscr{H}_{2}\left(K\right)$
\begin{equation}
 \frac{Z_{1} \left(\mathscr{H}_{1}\left(J\right)\right)}{f_{1}\left(J\right)}
   =
 \frac{Z_{2} \left(\mathscr{H}_{2}\left(K\right)\right)}{f_{2}\left(K\right)},
\end{equation}
where $f_{1}$ and $f_{2}$ are some nonsingular functions of interaction constants $J,K$
appearing in Hamiltonians $\mathscr{H}_{1}$ and $\mathscr{H}_{2}$, respectively.
Similar functional relations between the partition functions of two different  Ising  models
can be also obtained with the help of the  decoration, star-triangle or star-square transformations \cite{syozi}.
Duality combined with those transformations enables one to derive the expressions
for the critical temperature for a wide range of Ising models.
It should be stressed, however, that in the case of duality the total
number of spin interactions must be the same for both Hamiltonians,
$\mathscr{H}_{1}$ and $\mathscr{H}_{2}$,  which allows one to introduce the concept of dual lattice.
For example, the triangular lattice turns out dual to the honeycomb lattice
and the rectangular lattice is dual to itself.
Using the concept of dual lattice, the duality relation can be interpreted
as a relationship between high-temperature (low-temperature) expansion
for the model with Hamiltonian $\mathscr{H}_{1}$  and the low-temperature (high-temperature) expansion
for the model described  by $\mathscr{H}_{2}$.

The concept of duality has been extended by Wegner \cite{wegner} for a wider class of Ising models $M_{dn}$
on d-dimensional lattices characterized by a number $n = 1, 2, \ldots , d$,
where $n = 1$ corresponds to the Ising model with a two-spin interaction.
In this paper we construct and investigate analytically the dual model to the Ising model on a rectangular lattice
with the nearest- and next-nearest-neighbor interactions (NNNI model) \cite{zandvlietOriginal,bittner,lee,kim}.
Although this simple extension of the Onsager model has been widely investigated for decades now,
no analytical expression for the partition function of the model has been found, nor its dual model has been presented.
To this end we use the usual algebraic approach and express the partition function as the trace of a power of the transfer matrix.
In \Sref{sec:1} we employ the classic 2D Ising model on a rectangular lattice to introduce the formalism and some basic formulas necessary to investigate the duality of a more complex model in \Sref{sec:union-jack}. In particular, we show how the self-dual symmetry of this model can be related to the similarity transformation ($U$) of the transfer matrix. We show that the difference between two forms of the transfer matrix
can be related to a change of the direction in which the matrix was constructed and present the explicit form of $U$.
In \Sref{sec:union-jack} we apply the similarity transformation to the transfer matrix of the NNNI model.
We show that there exists  a simple  Ising model for which the transformed matrix is  the transfer matrix.
The Hamiltonian of this model contains a two-spin nearest-neighbor interaction between nodes of the Brick-Wall lattice
(which is topologically equivalent to a honeycomb lattice) and additional four-spin interactions (henceforth this model will be called the (2 + 4) BWI model).
The number of the four-spin interactions is equal to the number of bricks in the lattice.
Since  the two Hamiltonians have the same number of all interactions, the  model can be regarded
as the dual model to the NNNI model.
As a result we obtain the duality relation between the partition functions of these  two models.
Moreover,  we show that  the dual model  becomes  self-dual in the presence of an external magnetic field.
In ~\Sref{sec:tc} the critical temperature of the (2 + 4) BWI model is investigated.
We show that an increase in the critical temperature of the dual model caused by the four-spin interactions
is limited to a finite range. We derive the equations for both the lower and upper limits of this temperature range and find that they depend on the magnitude of the two-spin interactions only.

\section{Selfduality of the Ising model on the rectangular lattice \label{sec:1}}
Consider the classical Ising model with the nearest-neighbor interaction constants $J_{1}$ and $J_{2}$ (for simplicity, we assume that $J = J_{1} = J_{2}$)
on a rectangular lattice with  $N$ columns and $2M$ rows.
We assume the cyclic boundary conditions in the direction of the transfer matrix action
and in the perpendicular direction we adjust  the boundary conditions to the transformations
performed on the transfer matrix to get the result in a closed form
(we adopt this convention throughout this paper).
The partition function for this model reads \cite{schultz, onsager, thompson}
\begin{eqnarray}
\label{Z_1}
Z_{N,2M}(\kappa) &=& \operatorname{Tr}V^{N} \nonumber \\
               &=& \left(2 \sinh 2\kappa\right)^{MN} \operatorname{Tr}\left(e^{\kappa A}e^{\tilde{\kappa}C}\right)^{N},
\end{eqnarray}
where $\kappa = \beta J = J/k_{\mathrm{B}} T$,
$\tilde{\kappa} = -\frac{1}{2} \ln \tanh \kappa$,
and $A$ and $C$ are $2^{2M}\times 2^{2M}$ matrices defined by
\begin{eqnarray}
  A &=& \sigma_{1}^{x} + \sum\limits_{j = 1}^{2M - 1}{\sigma_{j}^{x} \sigma_{j+1}^{x}},\\
  C &=& \sum\limits_{j = 1}^{2M}{\sigma_{j}^{z}},
\end{eqnarray}
with
\begin{equation}
  \sigma_{j}^{\alpha} =
   \underbrace{ \mathds{1} \otimes \mathds{1} \otimes \ldots \otimes \mathds{1} }_{j-1}
      \otimes \sigma^{\alpha} \otimes
   \underbrace{\mathds{1} \otimes \ldots \otimes \mathds{1}}_{2M-j},
\end{equation}
where the symbol $\otimes$ denotes the tensor product, $\mathds{1}$ is the $2 \times 2$ unit matrix,
and $\sigma^{\alpha}$ ($\alpha = x, y, z$) represent the Pauli matrices
\begin{equation}\sigma^{x} =
\left(
\begin{array}{c c}
		0 & 1 \\
		1 & 0 \\
		\end{array}
\right),
\sigma^{y} =
\left(
\begin{array}{c c}
		0 & -i \\
		i & 0 \\
		\end{array}
\right),
\sigma^{z} =
\left(
\begin{array}{c c}
		1 & 0 \\
		0 & -1 \\
		\end{array}
\right),
\end{equation}
and $j$ runs through the $2M$ nodes of a lattice column.

It turns out that $A$ and  $C$ a related to each other by a similarity transformation $U$,
\begin{equation}
\label{eq:similarity}
   UAU^{-1} = C,  \qquad
   UCU^{-1} = A.
\end{equation}
Relations of this type were found to be  essential in proving the duality
or self-duality of many Ising models, including the self-duality of
the 1D Ising model with an external transverse magnetic field \cite{fradkin_susskind,kogut,turban}.

Applying \Eref{eq:similarity} under the trace operator in \Eref{Z_1}, one obtains
\begin{eqnarray}
\label{Eq:up_form_z}
 Z_{N,2M}(\kappa) &=&
   (2 \sinh 2\kappa)^{MN} \operatorname{Tr} \left( U \left(e^{\kappa A} e ^{\tilde{\kappa} C}\right)^{N} U^{-1} \right) \nonumber \\
  &=&  (2 \sinh 2\kappa)^{MN} \operatorname{Tr}  \left(e^{\kappa U A U^{-1}} e^{\tilde{\kappa} U C U^{-1}}\right)^{N} \nonumber \\
  &=& (2 \sinh 2\kappa)^{MN} \operatorname{Tr}  \left(e^{\kappa C} e^{\tilde{\kappa} A}\right)^{N}.
\end{eqnarray}
Using the cyclicity of the trace operator, \Eref{Eq:up_form_z} can be rewritten in a well-known form
\begin{eqnarray}
 \label{eq:prof-2D}
  Z_{N,2M}(\kappa) &=&  (2 \sinh 2\kappa)^{MN} \operatorname{Tr}  \left(e^{\tilde{\kappa} A} e^{\kappa C}\right)^{N} \nonumber \\
  & = &  \left(\frac{2 \sinh 2 \kappa}{2 \sinh 2 \tilde{\kappa}}\right)^{MN} (2 \sinh 2 \tilde{\kappa})^{MN} \operatorname{Tr} \left(e^{\tilde{\kappa} A} e^{\kappa C}\right)^{N}  \nonumber\\
  & = &  \left(\frac{2 \sinh 2 \kappa}{2 \sinh 2 \tilde{\kappa}}\right)^{MN} Z_{N,2M}(\tilde{\kappa})  \nonumber\\
  & = &  (\sinh 2 \kappa)^{2MN} Z_{N,2M}(\tilde{\kappa}),
\end{eqnarray}
which proves the self-duality of the discussed model.

\Eref{eq:prof-2D} was usually derived by constructing the transfer matrix
along different directions of the lattice, without deriving the explicit form of $U$,
which however will be necessary in our study further below.
To obtain it, one can assume that it can be expressed as a product of two $2^{2M} \times 2^{2M}$ matrices
\begin{equation}
 \label{eq:SR}
   U = SR,
\end{equation}
where $R$ is a permutation matrix,
whereas $S$ exchanges Pauli matrices  $\sigma^{x}_{i}$ with $\sigma^{z}_{i}$ and changes the signs
of the interaction coefficients appropriately.
A particularly simple form of $R$ reads
\begin{equation}
 \label{eq:R_1}
   R = \prod\limits_{j = 1}^{2M - 1}{R_{j,j + 1}},
\end{equation}
where
\begin{equation}
  R_{j, j + 1} =
  \underbrace{\mathds{1} \otimes \mathds{1} \otimes \ldots  \otimes \mathds{1}}_{j-1} \otimes r \otimes
  \underbrace{\mathds{1} \otimes \ldots \otimes \mathds{1}}_{2M-j-1}
\end{equation}
and
\begin{equation}
	r = \left(
		\begin{array}{c c c c}
		0 & 1 & 0 & 0 \\
		1 & 0 & 0 & 0 \\
		0 & 0 & 1 & 0 \\
		0 & 0 & 0 & 1
		\end{array}
	\right ).
\end{equation}
It follows that $r^{-1} = r$, $[R_{j,j+1},R_{j+1,j+2}] \neq 0$, and
\begin{equation}
  \label{eq:R-1_1}
    R^{-1} = \prod\limits_{l = 1}^{2M - 1}{R_{2M - l,2M - (l - 1)}}.
\end{equation}
The following relations
\begin{eqnarray}
R_{l, l + 1} \sigma_{l}^{x} \sigma_{l+1}^{x} R_{l, l + 1}^{-1} &=& \sigma_{l + 1}^{x} \\ \nonumber
R_{l, l + 1} \sigma_{l}^{x} R_{l, l + 1}^{-1} &=& \sigma_{l}^{x} \\ \nonumber
R_{l, l + 1} \sigma_{l+1}^{x} R_{l, l + 1}^{-1} &=& \sigma_{l}^{x} \sigma_{l+1}^{x} \\ \nonumber
R_{l, l + 1} \sigma_{l}^{z} \sigma_{l+1}^{z} R_{l, l + 1}^{-1} &=& - \sigma_{l}^{z} \\ \nonumber
R_{l, l + 1} \sigma_{l}^{z} R_{l, l + 1}^{-1} &=& -\sigma_{l}^{z} \sigma_{l+1}^{z} \\ \nonumber
R_{l, l + 1} \sigma_{l+1}^{z} R_{l, l + 1}^{-1} &=& \sigma_{l+1}^{z}
\end{eqnarray}
can be used to show how $A$ and $C$ transform under  $R$,
\begin{eqnarray}
\label{eq:R_2}
 RAR^{-1} &=&
    \sum\limits_{j = 2}^{2M} \sigma_{j}^{x}, \nonumber\\
 RCR^{-1} &=&
    \sum\limits_{j = 1}^{2M - 1} \sigma_{l}^{z}\sigma_{j + 1}^{z} + \sigma_{2M}^{z}.
\end{eqnarray}

The cyclic boundary conditions in  $A$ are mathematically troublesome.
However, since the thermodynamics of the system does not depend on the boundary conditions in the thermodynamic limit,
they can be neglected \cite{fradkin_susskind,kogut}.
As for $S$, it can be  expressed as
\begin{equation}
   S = \mathscr{P}\prod\limits_{j = 1}^{2M} s_{j} ,
\end{equation}
where
\begin{equation}
\label{eq:sj}
  s_{j} = \underbrace{\mathds{1} \otimes \mathds{1} \otimes \ldots \mathds{1}}_{j-1}
     \otimes s \otimes
  \underbrace{\mathds{1} \otimes \ldots \otimes \mathds{1}}_{2M-j},
\end{equation}
with
\begin{equation}
\label{eq:s}
s = \frac{1}{\sqrt{2}}\left(
\begin{array}{c c}
		1 & 1 \\
		1 & -1 \\
		\end{array}
\right),
\end{equation}
whereas $\mathscr{P}$ is a transformation exchanging matrices  $\sigma^{\alpha}_{i}$ with $\sigma^{\alpha}_{2M - i + 1}$,
\begin{equation}
\mathscr{P} = \prod_{i = 1}^{M}{P_{i,2M - i + 1}},
\end{equation}
where
\begin{equation}
P_{i,j} = \frac{1}{2}\left(1 + \sigma_{i}^{x}\sigma_{j}^{x} + \sigma_{i}^{y}\sigma_{j}^{y} + \sigma_{i}^{z}\sigma_{j}^{z}\right)
\end{equation}
and
\begin{equation}
P_{i,j} \sigma_{i}^{\alpha} P_{i,j}^{-1} = \sigma_{j}^{\alpha}, \mbox{ } \alpha = x,y,z \mbox{ and } i,j = 1,2,\ldots, 2M.
\end{equation}
On applying the full transformation $U = SR$  to $A$ and $C$, one arrives at \Eref{eq:similarity}.

\section{The dual model to the NNNI Model \label{sec:union-jack}}
The transfer matrix for the Ising model with the nearest  and next-nearearest neighbor
interactions was found in 1956 by Temperley \cite{temperley_book},
who  constructed it along the diagonal of the basic rectangular lattice
defined by the nearest-neighbor interactions.
Although this matrix is not of the lowest possible rank,
its form is particularly simple. We will use this matrix for the rhomboidal
lattice presented in \Fref{fig:lattice_1_2_3}.
\begin{figure}
\centering
  \includegraphics[width=0.9\columnwidth]{figure1.eps}
  \caption{The Ising model with the nearest and next-nearest neighbor interactions
  on a rhomboidal lattice (the NNNI model).
  $\bar{J_{1}}$, $\bar{J_{2}}$, and  $\bar{J_{3}}$ are the interaction constants for the diagonal,
  horizontal and vertical directions, respectively.  }
  \label{fig:lattice_1_2_3}
\end{figure}
While in most papers devoted to the NNNI model only two interaction constant are used \cite{bittner, lee},
$\bar{J_{1}}$, $\bar{J_{2}} = \bar{J_{3}}$, here we consider the general case
with three arbitrary constants controlling  the ferromagnetic interactions
$\bar{J_{1}}, \bar{J_{2}}, \bar{J_{3}}$ ($\bar{J_{i}} < 0$).

The transfer matrix and the partition function for this model on the rhomboidal  lattice
consisting of $2M \times N$
nodes satisfy
\begin{equation}
\frac{Z_{N,2M}}{\left(2 \sinh 2 L_{2}\right)^{NM}} =
 \operatorname{Tr} \left(e^{L_{1}A} e^{L_{3} B_{2\mathds{Z}-1}}
    e^{\tilde{L_{2}} C_{2\mathds{Z}}} e^{L_{1}A} e^{L_{3} B_{2\mathds{Z}}} e^{\tilde{L_{2}} C_{2\mathds{Z}-1}}\right)^{N},
\end{equation}
where $L_{i} = \beta \bar{J_{i}}$ $(i=1,2,3)$,
$\tilde{L}_{2} = -\frac{1}{2} \ln \tanh L_{2}$,
 $A$, $B_{2\mathds{Z}}$, $B_{2\mathds{Z}-11}$, $C_{2\mathds{Z}}$, and $C_{2\mathds{Z}-1}$
are matrices of size $2^{2M} \times 2^{2M}$, defined as follows,
\begin{eqnarray}
 A                 &=& \sigma_{1}^{x} + \sum\limits_{k = 1}^{2M - 1} \sigma_{k}^{x} \sigma_{k+1}^{x},  \\
 B_{2\mathds{Z}}   &=& \sigma_{2}^{x} + \sum\limits_{k = 1}^{M - 1} \sigma_{2k}^{x} \sigma_{2k+2}^{x}, \\
 B_{2\mathds{Z}-1} &=& \sum\limits_{k = 1}^{M - 1} \sigma_{2k - 1}^{x} \sigma_{2k+1}^{x} + \sigma_{1}^{x} \sigma_{2M - 1}^{x} \sigma_{2M}^{x}, \\
 C_{2\mathds{Z}}   &=& \sum\limits_{k = 1}^{M} \sigma_{2k}^{z}, \\
 C_{2\mathds{Z}-1} &=& \sum\limits_{k = 1}^{M} \sigma_{2k - 1}^{z},
\end{eqnarray}
with $2\mathds{Z}$ denoting the set of even numbers.
The correspondence between these matrices and the interaction constants are shown in  \Fref{fig:lattice_1_2_3_d}.
\begin{figure}
  \centering
  \includegraphics[width=1.0\columnwidth]{figure2.eps}
  \caption{Correspondence between the two-spin interactions and the matrix operators.}
  \label{fig:lattice_1_2_3_d}
\end{figure}
Similar transfer matrices for the NNNI model were already used in \cite{oitmaa},
however, without giving their explicit forms.

On applying the similarity transformation $U$
to the matrices of the NNNI model we arrive at
\begin{eqnarray}
UAU^{-1}                 &=& \sum\limits_{k = 1}^{2M} \sigma_{k}^{z} = C,  \\
UB_{2\mathds{Z}-1}U^{-1} &=& \sum\limits_{k = 2}^{M} \sigma_{2k - 2}^{z} \sigma_{2k-1}^{z} + \sigma_{1}^{z} \sigma_{2M}^{z} = D_{2\mathds{Z}}, \\
UC_{2\mathds{Z}}U^{-1}   &=& \sigma_{1}^{x} + \sum\limits_{k = 1}^{M - 1} \sigma_{2k}^{x} \sigma_{2k+1}^{x} = A_{2\mathds{Z}}, \\
UB_{2\mathds{Z}}U^{-1}   &=& \sum\limits_{k = 1}^{M} \sigma_{2k-1}^{z} \sigma_{2k}^{z} = D_{2\mathds{Z}-1}, \\
UC_{2\mathds{Z}-1}U^{-1} &=& \sum\limits_{k = 1}^{M} \sigma_{2k - 1}^{x} \sigma_{2k}^{x} = A_{2\mathds{Z}-1},
\end{eqnarray}
which also defines $A_{2\mathds{Z}}$,  $A_{2\mathds{Z}-1}$, $D_{2\mathds{Z}}$, and $D_{2\mathds{Z}-1}$.
These relations lead to
\begin{equation}
  \label{eq:duality-UJ}
  \frac{Z_{N,2M}}{\left(2 \sinh 2 L_{2}\right)^{NM}} = \operatorname{Tr} \tilde{V}^{N},
\end{equation}
where
\begin{equation}
	\label{eq:tildeV}
  \tilde{V} =
    e^{L_{1}C}e^{L_{3}D_{2\mathds{Z}}}e^{\tilde{L_{2}}A_{2\mathds{Z}}}
    e^{L_{1}C}e^{L_{3}D_{2\mathds{Z}-1}}e^{\tilde{L_{2}}A_{2\mathds{Z}-1}}.
\end{equation}

Equation (\ref{eq:duality-UJ}) ensures that if there exists an Ising model whose transfer matrix
is equal to $\tilde{V}$, this model will be dual to the NNNI model.
This is a nontrivial requirement, as while the  transfer matrix of any Ising model can be expressed as
a product of exponential matrices,
the product of  exponential matrices is hardly ever the transfer matrix of an Ising model.
In the case considered here, however, such a model exists and is described by the Hamiltonian
\begin{equation}
\mathscr{H} = - \sum\limits_{<i,j>} J_{ij}\sigma_{i}\sigma_{j} - J_{4} \sum\limits_{<i,j,k,l>}
\sigma_{i}\sigma_{j}\sigma_{k}\sigma_{l}
\label{eq:H_J_J4}
\end{equation}
defined on the Brick-Wall lattice \cite{temperley_chapter}, which is topologically equivalent
to the honeycomb lattice.
As explained in \Fref{fig:square_j4_j1_j2_b},
\begin{figure}
	\centering
		\includegraphics[width=0.75\columnwidth]{figure3.eps}
	\caption{The (2 + 4) BWI Model.}
	\label{fig:square_j4_j1_j2_b}
\end{figure}
$J_4$ is the four-spin interaction constant for the four spins occupying the right-hand side of a brick,
and the values of the nearest-neighbor interaction constants $J_{ij}$ are non-zero only if
$i$ and $j$ are on the same side of a brick, in which case they are equal $J_1$ and $J_2$
for interactions along the vertical and horizontal direction, respectively.
All interaction constants are ferromagnetic in nature ($J_1, J_2, J_4 < 0$).

The dual lattice is composed of $2M$ horizontal chains, each containing $2N$ nodes.
The partition function of the model satisfies
\begin{eqnarray}
\label{eq:Z_K1_K2_K4}
& &\frac{Z_{2N,2M} \left(K_{2},K_{4},K_{1}\right)}{2^{2MN}\left[\sinh^{2}
    2 K_{2} \left(e^{4K_{4}}\cosh^{2} 2 K_{2} - 1\right)\right]^{MN/2}} = \nonumber\\
& & = \operatorname{Tr} \left(e^{\tilde{K_{2}}C}e^{\tilde{K_{4}}D_{2\mathds{Z}}}e^{K_{1}A_{2\mathds{Z}}}e^{\tilde{K_{2}}C}e^{\tilde{K_{4}}D_{2\mathds{Z}-1}}e^{K_{1}A_{2\mathds{Z}-1}}\right)^{N},
\end{eqnarray}
where $K_{i} = \beta J_{i}$, $(i = 1,2,4)$ and
\begin{eqnarray}
 \label{eq:tilde_K2}
  \tilde{K_{2}} &=& - \frac{1}{4} \ln \frac{\cosh 2K_{2} - e^{-2K_{4}}}{\cosh 2 K_{2} + e^{-2K_{4}}}, \\
 \label{eq:tilde_K4}
  \tilde{K_{4}} &=& - \frac{1}{4} \ln \frac{\cosh^{2} 2K_{2} - 1}{\cosh^{2} 2 K_{2} - e^{-4K_{4}}}.
\end{eqnarray}

The transfer matrix  in \Eref{eq:Z_K1_K2_K4} acts on the vertical groups of nodes and advances
the partition function in the horizontal direction \cite{baxter}.
Its explicit form can be obtained from Onsager's classical transfer matrix for the rectangular lattice
\cite{onsager}, see \Eref{Z_1}. First, we add the four-spin interactions that modify
only two horizontal two-spin bonds lying one above the other.
Next, half of the vertical two-spin interactions connecting the horizontal
chains has to be removed.
While adding the four-spin interactions to the transfer matrix, we  can use a relationship
\begin{equation}
\Delta^{1/4} e^{\tilde{K_{2}}\left(\sigma_{k}^{z} + \sigma_{k+1}^{z}\right)}e^{\tilde{K_{4}}\sigma_{k}^{z}\sigma_{k+1}^{z}}
= \cosh K_{4} c_{k} c_{k+1} + \sinh K_{4} \sigma_{k}^{x} \sigma_{k+1}^{x} c_{k}c_{k+1} \sigma_{k}^{x} \sigma_{k+1}^{x},
\end{equation}
where
$\Delta = \frac{1}{16} \sinh^{2} 2K_{2} \left(e^{4K_{4}} \cosh^{2} 2K_{2} - 1\right)$ and
\begin{equation}
c_{k} = \underbrace{\mathds{1} \otimes \mathds{1} \otimes \ldots \otimes \mathds{1}}_{\mbox{k - 1}}
  \otimes c \otimes
        \underbrace{\mathds{1} \otimes \ldots \otimes \mathds{1}}_{\mbox{2M - k}},
\end{equation}
\begin{equation}
c =
\left(
\begin{array}{c c}
		\cosh K_{2} & 0 \\
		0 & \sinh K_{2} \\
		\end{array}
\right),
\end{equation}
and $k=1,\ldots,\mbox{2M}$.
Finally, we find the duality relationship between the NNNI model, \Fref{fig:lattice_1_2_3},
and the (2 + 4) BWI model, \Fref{fig:square_j4_j1_j2_b},
\begin{equation}
 \label{eq:duality-UJ-BM}
\frac{Z_{N,2M}\left(L_{1},L_{3},L_{2}\right)}{\left(2 \sinh 2 L_{2}\right)^{MN}} =
\frac{Z_{2N,2M} \left(K_{2}, K_{4}, K_{1}\right)}{2^{2MN} \left[\sinh^{2} 2K_{2}
 \left(e^{4K_{4}}\cosh^{2} 2 K_{2} - 1\right)\right]^{MN/2}},
\end{equation}
where the relations between the interaction constants for both model are given by
\begin{equation}
 \label{LLL-KKK}
  L_{1}         = \tilde{K_{2}}, \qquad
  L_{3}         = \tilde{K_{4}}, \qquad
  \tilde{L_{2}} = K_{1}.
\end{equation}

\section{The critical temperature for the (2 + 4) BWI model \label{sec:tc}}
While the Ising model on the Brick-Wall lattice can be solved exactly \cite{wannier_brickwall}, addition of the four-spin interaction
makes it intractable analytically.
However, it turns out that  it is still possible to draw some important conclusions about the way the  critical temperature
of the Ising model defined by Hamiltonian (\ref{eq:H_J_J4})
depends on the strength of the four-spin interactions.

In the Ising model with nearest  and next-nearest  neighbor interactions (NNNI model, \Fref{fig:lattice_1_2_3})
there are two important limiting cases: $\bar{J_{3}} = 0$ and $\bar{J_{1}} = 0$.
In the former case
the NNNI model reduces
to the Ising model with nearest-neighbor interactions on the  triangular lattice with the interaction constants
$(\bar{J_{1}}, \bar{J_{2}})$.
Since in this case $L_{3} = \beta \bar{J_{3}} = 0$,
matrices $e^{L_{3}D_{2\mathds{Z}}}$ and $e^{L_{3}D_{2\mathds{Z} - 1}}$
appearing in \Eref{eq:tildeV} become the  identity matrices and
the critical point of this model can be determined exactly.
It is worth noticing  that the duality transformation converts the transfer matrix of this model
into the transfer matrix of the Ising model with two nearest-neighbor interactions
on a honey comb lattice equivalent to the Brick-Wall lattice.
This explains why the dual model to NNNI model is defined on the  Brick-Wall lattice.
In the second case ($\bar{J_{1}} = 0$),
the NNNI model splits into two identical independent self-dual Ising models
with nearest-neighbor interactions on the rectangular lattice.
Similarly to the previous case, matrices $e^{L_{1}A}$ in \Eref{eq:tildeV}
become the identity matrices and the critical point for this model is exactly determinable.

\Eref{LLL-KKK},
which describes relationships between  the interactions of the NNNI  and (2+4) BWI models, reduces to
\begin{equation}
\label{eq:L3K4}
0 = L_{3} = \tilde{K_{4}} = - \frac{1}{4} \ln \frac{\cosh^{2} 2K_{2} - 1}{\cosh^{2} 2 K_{2} - e^{-4K_{4}}}
\end{equation}
for  $\bar{J_{3}} = 0$  and to
\begin{equation}
\label{eq:L1K2}
0 = L_{1} = \tilde{K_{2}} = - \frac{1}{4} \ln \frac{\cosh 2 K_{2} - e^{-2K_{4}}}{\cosh 2 K_{2} + e^{-2K_{4}}}
\end{equation}
for  $\bar{J_{1}} = 0$.
\Eref{eq:L3K4} is equivalent to
$K_{4} = 0$, whereas \Eref{eq:L1K2} can be satisfied only in the limit of $K_{4} \rightarrow \infty$.
The critical temperatures satisfy \cite{syozi}
\begin{eqnarray}
\label{eq:e_cosh_tanh}
e^{-4 \frac{J_{4}}{kT_{C}}} &=& \cosh^{2} \left(\frac{2 J_{2}}{kT_{C}}\right) \frac{\tanh^{2}
  \left(\frac{J_{1}}{kT_{C}}\right)}{\tanh^{2} \left(\frac{J_{1}}{kT_{C}}\right) + 1},          \\
\label{eq:e_Cosh_e_sinh}
e^{-4 \frac{J_{4}}{kT_{C}}} &=&
  \cosh^{2} \left(\frac{2 J_{2}}{kT_{C}}\right) - e^{4 \frac{J_{1}}{kT_{C}}} \sinh^{2} \left(\frac{2 J_{2}}{kT_{C}}\right),
\end{eqnarray}
for $\bar{J_{3}} = 0$  and  $\bar{J_{1}} = 0$, respectively.

Figure~\ref{fig:J4_T} presents a sketch of the
relation between $J_{4}$ and the critical temperature $kT_\mathrm{C}$
in a particular case of $J_{1} = J_{2} = 1$.
\begin{figure}
	\centering
		\includegraphics[width=0.9\columnwidth]{figure4.eps}
	\caption{Relation between the critical temperature ($kT_\mathrm{C}$) and the four-spin interaction constant ($J_4$)
             in the dual model with $J_1=J_2=1$.
             The solid and dotted lines are the solution of Eqs. (\protect\ref{eq:e_cosh_tanh}) and
             (\protect\ref{eq:e_Cosh_e_sinh}), respectively,
             and the crosses indicate a possible location of the phase transition.
             $kT_\mathrm{CBW}$ is the critical temperature for the Ising model on a Brick-Wall lattice
             and $kT_\textrm{Cmax}$ is the maximum critical temperature for the ferromagnetic case of the (2 + 4) BWI model.
             }
	\label{fig:J4_T}
\end{figure}
The solid line represents the asymptotic solution for the case $K_{4} \rightarrow 0$ obtained from \Eref{eq:e_cosh_tanh},
and the dotted line depicts the asymptotic solution of \Eref{eq:e_Cosh_e_sinh} obtained for $K_{4} \rightarrow \infty$.
The solution of \Eref{eq:e_Cosh_e_sinh}
has only the left-hand vertical asymptote at the temperature $kT_\mathrm{Cmax}$  given as the solution to
\begin{equation}
\label{eq:Cosh_e_Sinh}
  \cosh \left(\frac{2}{kT_\mathrm{Cmax}}\right)
    =  \exp\left(\frac{2}{kT_\mathrm{Cmax}}\right) \sinh \left(\frac{2}{kT_\mathrm{Cmax}}\right),
\end{equation}
which yields $kT_\mathrm{Cmax} \approx 3.28204$.
The fact that \Eref{eq:e_Cosh_e_sinh} has only the left-hand side asymptote at the critical point  defined by
$K_{4} \rightarrow \infty$
is in contradiction to the Griffith's inequalities \cite{griffiths}, which imply  that
if an additional ferromagnetic interaction is introduced into a ferromagnetic Ising model,
the critical temperature  $T_\textrm{c}$ will increase with the magnitude of this additional interaction.
Therefore, this solution has to be rejected as  nonphysical.
Our preliminary computer simulations (which will be discussed elsewhere)
suggest that the relation between  $J_{4}$ and $T_{C}$ is located along the line
marked with crosses in \Fref{fig:J4_T},
with the right vertical asymptote at the point given by \Eref{eq:Cosh_e_Sinh}.

Our considerations show that the additional four-spin interaction
in the Ising model on a Brick Wall lattice does not yield an infinite increase of the critical temperature.
Instead, the critical temperature is always limited by a condition involving two-spin interactions only.
In order to confirm this qualitative result obtained for a particular choice of $J_1$ and $J_2$,
we plan to perform detailed computer simulations of the model.

\section{Self-duality of the (2 + 4) BWI model in the presence of an external magnetic field}
As in other pairs of Ising models, the inclusion of an external magnetic field makes
the duality relationship between the NNNI  and the (2 + 4) BWI models cease to exist:
duality relations (\ref{eq:duality-UJ-BM}) -- (\ref{LLL-KKK}) were obtained only in the absence of an external magnetic field.
However, the Ising model described by the Hamiltonian (\ref{eq:H_J_J4}) has an additional internal symmetry.
If one adds to the Hamiltonian an interaction with an external magnetic field $-H \sum\limits_{i} \sigma_{i}$,
then
 the partition function will take the form
\begin{equation}
\frac{Z_{2N,2M}\left(K_{2},K_{4},h,K_{1}\right)}{2^{2MN} \left[\sinh^{2} 2K_{2}
    \left(e^{4K_{4}} \cosh^{2} 2K_{2} - 1\right)\right]^{MN/2}} =
       \operatorname{Tr} \left(V \right)^{N},
\end{equation}
where
\begin{equation}
V = e^{\tilde{K_{2}}C}e^{\tilde{K_{4}}D_{2\mathds{Z}}}e^{hE}e^{K_{1}A_{2\mathds{Z}}}
e^{\tilde{K_{2}}C}e^{\tilde{K_{4}}D_{2\mathds{Z}-1}}e^{hE}e^{K_{1}A_{2\mathds{Z}-1}}
\end{equation}
is the transfer matrix, $h = \beta H$, and
$E$ is a matrix defined as
\begin{equation}
E = \sum\limits_{k = 1}^{2M}\sigma_{k}^{x}.
\end{equation}

Let
\begin{equation}
  \mathscr{S} = \prod\limits_{i = 1}^{2M} s_{i},
\end{equation}
where $s_{i}$ are given by \Eref{eq:sj}.
This operator exchanges $\sigma_{k}^{x}$ with $\sigma_{k}^{z}$. Applying it as a similarity transformation for $V$
leads to
\begin{equation}
  V' = \mathscr{S} V \mathscr{S}^{-1},
\end{equation}
where
\begin{equation}
  V' = e^{\tilde{K_{2}}E}e^{\tilde{K_{4}}A_{2\mathds{Z}}}
          e^{hC}e^{K_{1}D_{2\mathds{Z}}}e^{\tilde{K_{2}}E}e^{\tilde{K_{4}}A_{2\mathds{Z}-1}}e^{hC}e^{K_{1}D_{2\mathds{Z}-1}}.
\end{equation}
Thus, $V'$ can be interpreted as the transfer matrix of the (2 + 4) BWI model in which the interactions have been changed
as follows,
\begin{eqnarray}
h &\rightarrow& \tilde{K_{2}},  \nonumber\\
K_{1} &\rightarrow& \tilde{K_{4}}.
\end{eqnarray}
The second difference is that while in the initial Ising model the four-spin interactions
act on the right-hand side of the bricks,
the four-spin interactions appear on their left-hand sides of the transformed model.
The self-duality relation for this  model is described by the equation
\begin{equation}
Z_{2N,2M} \left( K_{2}, K_{4}, h, K_{1} \right) =
\alpha Z_{2N,2M} \left( \tilde{K_{2}}, \tilde{K_{4}}, \tilde{h}, \tilde{K_{1}}\right),
\end{equation}
where
\begin{equation}
  \alpha = \left[\frac{\sinh^{2} 2 K_{2} \left( e^{4 K_{4}} \cosh^{2} 2 K_{2} -1 \right)}{\sinh^{2} 2 \tilde{h}
    \left(e^{4 \tilde{K_{1}}} \cosh^{2} 2 \tilde{h} - 1 \right)}\right]^{\frac{MN}{2}}
\end{equation}
and
\begin{eqnarray}
\label{eq:RSD}
 \tilde{K_{2}} &=& -\frac{1}{4} \ln \frac{\cosh 2K_{2} - e^{-2K_{4}}}{\cosh 2K_{2} + e^{-2K_{4}}}, \nonumber \\
 \tilde{K_{4}} &=& -\frac{1}{4} \ln \frac{\cosh^{2} 2K_{2} - 1}{\cosh^{2} 2K_{2} - e^{-4K_{4}}}, \nonumber\\
 \tilde{K_{1}} &=& -\frac{1}{4} \ln \frac{\cosh^{2} 2h - 1}{\cosh^{2} 2h - e^{-4K_{1}}}, \nonumber \\
 \tilde{h} &=& -\frac{1}{4} \ln \frac{\cosh 2h - e^{-2K_{1}}}{\cosh 2h + e^{-2K_{1}}}.
\end{eqnarray}

Unfortunately the model under consideration is a ferromagnetic one,
so the phase transition in an external magnetic field does not exist.
Self-duality enables one  only  to find a number of interesting relations between
correlations function for the (2 + 4) BWI model.

\section{Summary \label{sec:summary}}

The main goal of this paper was to find and investigate the dual model for the
2D Ising model with  the isotropic nearest and anisotropic next-nearest neighbor interactions  on a rectangular lattice.
The dual model  turned out to be a 2D Ising model with two-spin nearest-neighbor anisotropic
interactions and additional four-spin interactions on the Brick-Wall lattice.
The appearance of the four-spin interactions is associated with the presence of non-planar
next-nearest neighbor interactions  in the original model.
The way we constructed the dual model is quite general, does not require
the introduction of the dual lattice concept and can
be applied to other Ising models.

Investigation of the impact of the four-spin interactions on the critical temperature
of the dual model revealed that while the four-spin interactions can increase
the critical temperature of the model, this increase is restricted to a finite temperature range.
Moreover, the limits of this range are bounded by the two-spin interaction constants of the dual model.
This rather unexpected result has to be analysed further, e.g.\ by computer simulation, for example,
 using Landau's approach \cite{landau}.
Computer simulations are also necessary to investigate the exact location of the phase transition as well as its
critical properties, including the universality class.
They will also be useful in relating the results for the dual model with the properties of the original Ising model.

\ack
We would like to thank J. J\c{e}drzejewski  for his valuable comments,
especially on the important consequences of  the Griffiths inequality.
\section*{References}


\section{Introduction}

The conventional picture of the interaction between galaxies and their environments is based on the idea that galaxies enter dense environments with a reservoir of gas.  The removal of this reservoir, either from the stellar disk (ram pressure stripping), or the galaxy's halo (strangulation), leads to a decline in the star formation rate. This view is out of date with our understanding of the dynamic baryon cycle in galaxies. Cosmological simulations, as well as semi-analytic models, show that galaxies grow as a result of gas infall from surrounding filaments. At least at high redshift, rapid star formation is fed by a supply of fresh, infalling gas which dominates over the consumption of any reservoir \citep[eg.,][]{papovich_2011_short}. Taking this into account leads to a significant update to our picture of how galaxies interact with their environment. If this scenario is correct, simply removing the supply of fresh material will have a dramatic effect on the properties of satellite galaxies, without the need for gas stripping processes. In this Letter, we use a compilation of observed satellite star formation timescales at a range of redshifts, to examine the mechanisms which quench satellites and the insights this gives to the baryon cycle of all galaxies.

In recent years, many authors have combined direct observations of satellite quenched fractions with cosmological infall rates of galaxies to estimate the time it takes a satellite to quench; the so-called ``quenching time'' \citep[eg.,][]{mcgee_11_short, delucia_infall, wetzel_timescale}. This is interpreted as the length of time after ``accretion'', where accretion is loosely associated with the time a galaxy is within the virialized region of some larger structure, until the galaxy is completely quenched. At low redshift, this timescale is long (several Gyrs) \citep{wetzel_timescale}, leading to the suggestion that gentle stripping of the gas reservoir is the mechanism responsible. However, such a picture naturally leads to most satellite galaxies having relatively paltry star formation rates as their gas is dwindled away.

In direct contrast, observational studies of satellites have been hampered by the lack of such `smoking guns' of galaxy transformation. Indeed, many recent studies have shown that those satellites which are forming stars are doing so at a rate indistinguishable from central galaxies of similar mass \citep[eg.][]{peng10, mcgee_11_short, wetzel, muzzin_gclass}. To parametrize this behaviour, \citet{wetzel_timescale} introduced the notion of a `delay time' -- that is, the length of time after infall for which a satellite galaxy appears to track the star formation rate of an analogous central galaxy\footnote{The existence of this delay time is still a matter of some debate as several authors have found that star forming galaxies have lower specific star formation rates in satellites than in centrals \citep[eg.][]{patel_2009, vulcani_2010, Haines_2013}. However, some of this discrepancy results from the definition of star forming galaxies, where authors who have more inclusive definitions find the discrepancy. Nonetheless, the existence of a tail of low star formation rate galaxies does not affect the overall observation that most satellites are forming stars at normal rates for a significant time.}.  In context of a rapid baryon cycle, the length of this delay time holds clues for our understanding of galaxy formation as it is expected that the fresh gas supply is also cut off near infall. Dark matter simulations show that the accretion of mass into a subhalo stops when it reaches 1.8 R$_{\mathrm{vir}}$ (as defined in \citep{Bryan_norman_98}) of a larger halo \citep{behroozi_accretion}. This limits the reservoir of gas which can be used for future star formation of the satellite, even in the absence of stripping.

The final relevant timescale of satellite quenching is a `fading time'. This is equal to the quenching time minus the delay time, and is interpreted physically as the time it takes a satellite to be quenched once its star formation rate deviates from that of an analogous, central galaxy. This phase is short ($<$ 1 Gyr) and is reflected in the relative lack of `green' or intermediate star forming galaxies.

In this letter, we compile measurements of the quenching and delay times of satellites at a range of redshifts to determine if `orbit-based' or `outflow-based' models dominate the quenching. Further, we can directly put limits on the outflow rates allowed by such delay times. We finish by exploring the difficulty of achieving long delay times in standard models of galaxy formation. In this letter, we adopt a \LCDM cosmology with the parameters: $\Omega_{\rm m} = 0.27$, $\Omega_{\Lambda}=0.73$ and $h=H_0/(100 \kmsmpc)=0.70$. All stellar masses assume a \citet{chabrier_imf} initial mass function.

\section{Observed satellite quenching timescales}\label{sec-obsquench}

The best current measurement comes from the analysis of \citet{wetzel_timescale} from a low redshift (z$\sim$0.05) SDSS sample of galaxy groups, where the activity level of the galaxies is determined from a H$\alpha$-based spectroscopic measurement. Wetzel et al. have carefully tracked the infall history of satellites at this epoch, accounting for the evolution in infalling field star formation rates, to find that a satellite galaxy with a stellar mass of 10$^{10.5}$ \Msun has a total quenching time of 4.4 $\pm$ 0.4 Gyrs. They find that 3.8 $\pm$ 0.4 Gyrs of this time is due to the `delay time', while the remaining $\sim$ 0.6 Gyrs of the quenching time is the `fading' stage. \citet{wetzel_timescale} define the moment a galaxy becomes a `satellite' as the time the subhalo crosses 200 times the matter background density (R$_{200,b}$), not the critical density (R$_{200,c}$). This R$_{200,b}$ definition is similar to the R$_{\mathrm{vir}}$ used in \citep{behroozi_accretion}. These densities evolve differently with redshift, so we consistently assume the Wetzel et al. definition and make corrections when the observational results were presented differently. At z=0.0 this matter density radius is 1.92 times larger than the critical density radius, although by z=1.0 it is only 1.15 times larger.

Using data from the GEEC2 sample of spectroscopic galaxy groups at z $\sim$ 0.9, \citet{mok_timescale} presented an estimate of the quenching and delay timescales by examining the relative fractions of `star forming', `intermediate' and `quiescent' satellites as defined by their locations in a colour-colour plane designed to separate the contributions from age and dust. Using an infall model and parametrization of the star formation histories of central galaxies, these authors concluded that these z $\sim$ 0.9 satellites are best fit by a total quenching time of 1.05 $\pm$ 0.25 Gyrs, made up of a 0.8 $pm$ 0.25 Gyr delay time and the remaining  $\sim$ 0.25 Gyr due to the fading time. This timescale applies to 10$^{10.5}$ \Msun\ galaxies, but is not seen to be a strong function of stellar mass for this sample.

From a recent comparison between the phase space distribution of post-starburst galaxies in z $\sim$ 1 galaxy clusters from the GCLASS sample and the phase space of subhalos accreted into cosmological-zoom simulations of galaxy clusters, \citet{muzzin_timescale} find that the total timescale for quenching is $\sim$ 1 $\pm$ 0.25 Gyr. The characteristic phase space features of these post-starburst galaxies can be best represented if quenching occurs quickly (0.1 - 0.5 Gyrs) following the first passage of 0.25 - 0.5 R$_{200,c}$. This is consistent with an analysis of the spectral features of the post-starbursts, selected to have low D$_n$(4000) and a lack of [OII] emission, which suggests that they have a short fading time ($\sim$ 0.1 - 0.5 Gyr). This implies that the delay time is the travel time from R$_{200,b}$ to this position. Directly from this simulation, this delay time is then 0.45$\pm$0.15 Gyrs depending on which inner distance the quenching occurs at. While the Muzzin et al. timescale is not presented in bins of stellar mass, they do state that most of the relevant galaxies are in the range of 10$^{10.5}$ \Mstel, and thus should be comparable to our previous timescales.

We also add a quenching time derived from the \galex-based star formation rates of galaxy groups at z$\sim$0.4 from \citet{mcgee_11_short}. We can define the excess quenching between groups and the field in a manner similar to \citet{vandenbosch_2008} as $\frac{f_{q,grp} - f_{q,fld}}{1 - f_{q,fld}}$, where $f_{q,grp}$ and $f_{q,fld}$ are the quenched fraction of the groups and the field, respectively. In 10$^{10.5}$\Msun galaxies, this excess quenching is 0.32 $\pm$ 0.1. By comparison with the accretion histories as compiled in \citet{mcgee_accretion_short}, we find that at z=0.4, 30$\%$ of galaxies became satellites more than 4.0$\pm$ 0.6 Gyrs ago. Although these authors do not explicitly measure the fading time, for our purposes we assume it is also 0.5 Gyrs. We caution that even minor differences in the criteria to demarcate between star forming and passive can lead to large differences within different surveys. In the future, a more comprehensive and homogeneous analysis of the existing data will be useful.

\section{Analysis}

\subsection{Orbit-based or outflow-based quenching?}

The time it takes a galaxy to proceed from the virial radius of a galaxy cluster to a specific location (e.g. perigee) in its orbit depends on redshift. Perhaps the simplest hypothesis is that the efficiency of gas stripping depends only on the orbit. That is, the same satellite will be stripped of the same fraction of gas at perigee, in the same orbit, regardless of whether this occurs at $z=0$ or $z=1$. In this case, the quenching timescale of galaxies from `orbit-based' quenching mechanisms will evolve as the ratio of the inverse densities, $\propto$ (1 + z)$^{-3/2}$.\footnote{The crossing time of a cluster is the ratio of the radius, R, and the velocity, V. For a cluster in virial equilibrium the cluster mass, M, scales as RV$^2$ and for a spherical cluster as R$^{3}\rho$. So, the dynamical time, R/V, scales as $\rho^{-1/2}$, while the density, $\rho$, scales as (1+z)$^3$.}

In contrast, an outflow-based quenching mechanism should have a timescale which scales as the star formation rate. Once the cosmological accretion of gas is halted from penetrating into the sphere of influence of a satellite galaxy, the timescale on which the galaxy quenches is equal to the total gas available at this point divided by the `gas consumption' rate. While this gas consumption rate is usually taken to be equal to the star formation rate of the galaxy, this ignores the outflows of gas associated with this star formation.  It is partly for this reason that satellite galaxies in semi-analytic models are quenched more quickly than the star formation rate alone would suggest. In this scenario, satellites prematurely halt their growth because they continue vigorous star formation, and drive associated outflows, despite a significant restriction on their available resources.  We therefore borrow a term from the field of ecology and refer to this quenching mechanism as ``overconsumption''.

For our purposes in this section, the exact outflow rate is not important if it is invariant in redshift for a given star formation rate (SFR) and stellar mass (M$_{stel}$). In that case, any outflow-based quenching would require the quenching time to scale with the star formation rate, or the specific star formation rate (sSFR = SFR/M$_{stel}$). Several authors have presented fitting formula to the observed star forming `sequence' of galaxies over cosmic time. In particular, \citet{peng10}, found that
\beq
\mathrm{sSFR} = 2.5 \left(\frac{M_{stel}}{10^{10} M_\odot}\right)^\beta \left(\frac{t}{3.5 \mathrm{Gyr}}\right)^{-2.2}
\eeq
where $t$ is the age of the universe and $\beta$ describes the tilt of the specific star formation rate relation, which we take to be 0. Similarly, \citet{whitaker_2012} parametrized the star formation rate as
\beq
\mathrm{log}(\mathrm{SFR}) = \alpha(z)(\mathrm{log} M_{stel} - 10.5) + \beta(z),
\eeq
where $\alpha(z) = 0.70 - 0.13z$ and $\beta(z) = 0.38 +1.14z - 0.19z^2$. We include these two relations as a measure of the uncertainty in such evolutionary measurements of the sSFR, and note that such relations rely heavily on the assumed functional form at very low and very high redshift where the fit is not well constrained by the data \citep[eg.,][]{Behroozi_2013}.

In Figure \ref{fig_scaling}, we show the evolution of quenching time expected from orbit-based and outflow-based quenching (using either the Peng et al. or Whitaker et al. star formation rates with \Mstel\ = 10$^{10.5}$\Mdot) when normalized to a quenching time of 4.4 Gyr at z=0.05. These are shown along with the observed quenching times described in \textsection \ref{sec-obsquench}. While the lower redshift point is roughly consistent with the dynamical time evolution, the z $\sim$ 1 quenching times seem to be lower than expected. This could be due either to an onset of outflow-based quenching or an enhanced stripping efficiency at high redshift. We note that \citet{tinker_wetzel_2010} and \citet{tinker_2013_short} have concluded that their halo model/clustering measurements are consistent with dynamical evolution, although they do not compare with SFR-derived timescales. In the next section, we will determine the outflow rates required to explain these quenching times.

\begin{figure}
\hspace{-0.7cm}
\leavevmode \epsfysize=6.3cm \epsfbox{figs/orbits_outflows_quenchtimes.ps}
\caption{The evolution of satellite quenching time with redshift. The lines represent the expectation from an `orbit-based' or `outflow-based' quenching mechanism. The data points represent the observed quenching times of \citet{wetzel_timescale}, \citet{mcgee_11_short}, \citet{mok_timescale} and \citet{muzzin_timescale} as described in \textsection \ref{sec-obsquench}.}
\label{fig_scaling}
\end{figure}

\subsection{Limits on the outflow rate}

As discussed in the introduction, one of the most curious observed aspects of satellite galaxies is that those which are forming stars appear to have similar star formation rates to central galaxies of similar mass -- requiring a `delay time'. Given the halting of cosmological accretion, the largely unchanging star formation rate during the delay time must result from an unchanged fuelling rate of the cold gas disk from the halo's gas reservoir, and will result in an unchanged outflow rate driven by this star formation. Thus, even if the delay time is established by the onset of orbit-based quenching, there must be a sufficient reservoir at satellite infall to maintain this fuelling throughout. Given the delay time, we can determine the maximum allowed outflow rate from basic baryon accounting measures.

The inflow and outflow of gas from a galaxy halo can be parametrized as
\beq
\dotMin + \dotMres = \dotMstar (1 - R) + \dotMout
\label{equilibrium_equation}
\eeq
where $\dot{M}_x$ is the rate of gas accretion due to cosmological infall (x = $in$), the rate of change of the gas reservoir of the system (x = $res$), the star formation rate (x = $stel$) and the rate of outflows from feedback (x = $out$). $R$ is fraction of gas returned to the ISM through supernova, stellar winds, etc., and we assume that this is done instantaneously. This basic equation has been the basis of galaxy formation models for some time \citep[eg.,][]{Colemodel}, and has recently been recast as the `equilibrium condition' by several authors \citep[eg.,][]{Bouche_2010, Dave_2012, Lilly_2013}, although the precise definition varies depending on the definition of the system (eg., as either the galaxy alone, or the galaxy and its extended halo) and whether a `reservoir' exists. Notice that a non-zero delay time immediately suggests that a sizeable reservoir must exist upon satellite infall.

Assuming that cosmological accretion is halted on infall, then $\dotMin$ = 0. As noted previously, there are simulations which suggest this happens to the dark matter prior to R$_{\mathrm{vir}}$ and, in such cases, we would overestimate the maximum outflow rate. The maximum time that the star formation rate can be maintained at a constant rate is

\beq
T_{\mathrm{delay}} = \frac{M_{res}}{\dotMres}.
\eeq
Taking account of the possible locations of the halo's baryon content at the point of infall, M$_{res}$ can be given from a conservation equation,
\beq
M_{res} = M_{baryon} - M_{cold} - M_{stel} - M_{ej} - M_{strip},
\label{conservation_equation}
\eeq
where $M_{baryon}$ is the total mass of baryons that {\it belong} to the system, and can be taken as the universal fraction of baryons multiplied by the satellite halo mass ($f_{baryon}*M_{halo}$). $M_{ej}$ is the mass of baryons that have been ejected from this system prior to infall. In our discussion, we will assume that baryons ejected are lost from the system, with no re-accretion. By this definition the integral of $\dotMout$ from its formation time to infall is equal to $M_{ej}$. This is largely just a definition of our terms, and therefore defines the outflow rate to be the rate at which gas is expelled from the system forever. $M_{strip}$ is the mass of gas which could be stripped during the lifetime of the quenching phase, and can be held to zero to find the maximum possible outflow rate. $M_{cold}$ is the cold gas currently locked up in the galaxy. It is likely that this mass must be approximately constant throughout the delay time in order to produce an approximately constant star formation rate according to the Kennicutt-Schmidt law \citep{Schmidt_1959, kennicutt}. Perhaps it is the eventual depletion of this gas that results in the `fading time'.

In the case of outflows produced from star formation, either directly through radiation pressure or from subsequent supernova explosions, the outflow rate is related to the star formation rate and is often parametrized as a mass-loading factor, $\eta$ ($\equiv$ $\dotMout$/$\dotMstar$). Using the equilibrium equation (equation \ref{equilibrium_equation}), along with the conservation equation (equation \ref{conservation_equation}) and the definition of the mass loading factor, the delay time is

\beq
T_{\mathrm{delay}} = \frac{M_{baryon} - M_{cold} - M_{stel} - M_{ej} - M_{strip}}{M_{stel} (1 - R+\eta) sSFR}
\eeq
where we have defined the specific star formation rate (sSFR = $\dotMstar$/M$_{stel}$). For the purposes of parametrizing our knowledge of the baryonic components it is convenient to define each component with respect to the halo mass at infall, so that $f_x$ is $M_x$/$M_{halo}$. $R$, the returned gas fraction, is usually on the order of $\sim$ 0.25-0.4, depending on the choice of initial mass function and the time since the stellar population was born. We assume a value relevant for a Chabrier IMF with a 2 Gyr population (R = 0.4).

We can now estimate each of these components. In our assumed cosmology, $f_{baryon}$ = 0.17. There have been many attempts to quantify $f_{stel}$, largely through halo abundance matching techniques. It has been shown that this fraction is constant with redshift, assuming a universal initial mass function, although it does strongly depend on halo mass \citep[eg.,][]{moster_hod, guo_2011_efficiency, Behroozi_2013, mcgee_imf}. For our canonical 10$^{10.5}$\Msun\ stellar mass galaxy, using the functional form presented by \citet{moster_hod}, we find that this stellar fraction is 0.03. Although we adopt $f_{stel}$ =0.03, it is worth noting that \citep{Behroozi_2013} find similar stellar fractions in the range 0.024 to 0.028 to z$\sim$ 3 at this stellar mass. \citet{carilli_walter_2013} have compiled existing measurements of the H$_{2}$ content of strongly star forming galaxies at a range of redshifts and found that $f_{cold}$ = 0.1 * (1+z)$^2$ $f_{stel}$. As the H$_{2}$ is expected to correlate well with the star formation rate, this is an appropriate value for the cold gas content, rather than the sum of molecular and atomic hydrogen gas. In our parametrization atomic hydrogen will be part of the gas reservoir. If we assume that the ejected gas is a result of past outflows from star formation rates, then $f_{ej}$ = $\eta$ $f_{stel}$ (1+$R$), where the $R$ is needed to recover the total past star formation rate. This model of past gas ejection is clearly an approximation, as earlier star formation occurred as the halo was forming and thus $\eta$ is unlikely to be constant throughout. Indeed, given the lower potential wells, and the need to widely eject metals into the intergalactic medium to match observation, this $\eta$ may be an underestimate of the amount of ejected material. As well, this assumes that the star formation--driven outflows dominate the ejected gas budget.

Using these approximations, the delay time is

\beq
 T_{\mathrm{delay}} = \frac{f_{baryon} - f_{stel}(0.1*(1+z)^2 + \eta (1+R) + 1) - f_{strip}}{f_{stel} (1-R+\eta) sSFR}
\label{equ-delay}
\eeq

We note that Equation \ref{equ-delay} has only one free parameter, $\eta$, and so can be completely characterized by a measure of the satellite delay time. In Figure \ref{fig_delay}, we show the maximum delay times (assuming $f_{strip}$ = 0) as a function of redshift for $\eta$ = 1, 2, 2.5, 3. We have also plotted the delay times from the observations, and a dynamical scaling which fits the delay time at z=0.05. Notice that the scaling with any particular $\eta$ is much steeper than the dynamical scaling because the evolution of the observed sSFR is very rapid. The turnover in delay time at $\eta$ = 3 and z $>$ 1.0, or $\eta$ = 2.5 and z $>$ 2.0, results from the need to maintain ever larger fractions of cold gas in the disk, while past outflows caused by the high $\eta$ make this impossible. The incompatible slope of the star formation rate and that implied by the dynamical time is a long-standing problem for models, and may suggest an underlying change in the efficiency of star formation, either through varying wind re-accretion or small scale cloud conditions.

Interestingly, in Figure \ref{fig_delay}, we see that the z $\sim$ 1 observations are consistent with a maximum mass loading of 2.5. Observations of strongly star forming galaxies suggest that the mass-loading is on the order of unity \citep{Heckman_2000, martin_outflows_2005}, although the measurements are difficult to interpret. While the ultimate fate of this outflowing gas is unclear, some simulations suggest that up to a third of the gas is re-accreted \citep{oppenheimer_dave_2008}, although this does not necessarily apply in satellites. In semi-analytic models, which are designed to reproduce the luminosity/mass function, the mass loading factors are typically much greater than one, and are a strong function of halo mass \citep[eg.,][]{guo_2011, Bower2012}. Although much of this gas may also be re-accreted, it is again not clear if this is the case in satellites \citep{font}. However, notice that with even moderate mass loading of 1, the expected delay time is shorter than expected from a scaling of the dynamical time in z $>$ 1.5 satellite galaxies. This is a potentially new quenching mechanism for satellite galaxies, and calls for a systematic observational study of such galaxies.

\begin{figure}
\hspace{-0.7cm}
\leavevmode \epsfysize=6.0cm \epsfbox{figs/orbits_outflows_delaytimes.ps}
\caption{The maximum possible delay times inferred as a function of redshift for a given mass loading factor ($\eta$ = 1, 2, 2.5, 3) and assuming a \citet{whitaker_2012} star formation history of star forming galaxies. The dynamical time scale is also shown. }
\label{fig_delay}
\end{figure}

\subsection{Implications for the cooling model}

In the previous section, we have implicitly assumed that the only necessary condition to sustain star formation at a constant rate is sufficient gas in the reservoir. However, such a constant SFR also requires that the cooling rate of gas from the reservoir is unchanged throughout the delay time. In a common model of cooling gas, the bulk of the gas reservoir is in the form of a hot halo, with a temperature that is a result of shock heating to the virial temperature upon accretion. This gas cools, largely through bremsstrahlung, and the rate depends on the square of the gas density. In simple models, as the gas near the central, dense region cools, the gas will respond to the flow quickly, so that the gas density profile shape will be driven by the dark matter potential. This seems like the best case scenario, as it is expected that the hot gas may build up outside the cooling radius \citep{fabian_1977}.) Thus, while the profile shape is constant, the decline in the total gas content will cause the central density to decline. The cooling rate is then directly linked to the amount of the gas in the reservoir. As the hot gas reservoir begins to drain, the cooling rate will go down, the subsequent amount of cold gas will decline and thus, star formation will also decline.

One solution to maintain uniform cooling rates is to change the shape of the gas profile as it is drained. It is possible that the immersion of a satellite galaxy into the hot, dense intracluster medium will compress the hot halo. However, given the need to match the compression rate with the cooling rate regardless of the exact density of the ICM, this suggests a fine tuning problem. Unfortunately, it is very difficult to observe X-ray halos of galaxies, and thus to determine this directly. One alternative to such a fine tuning is that the accreted gas is never shock heated to the virial temperature, and simply adds to the atomic gas reservoir. The thermal history of gas in galaxies has long been studied \citep{white_rees}, and is still debated \citep[eg.,][]{keres_2005, Dekel_birnboim}. Interestingly, though, this may hint that the cooling rate is governed by the interplay between atomic and molecular gas rather then the hot gas cooling. As we have mentioned, our definition of the gas `reservoir' includes the HI content of gas. In principle, if this is the primary gas reservoir, we would expect to see a decrease in the HI content for satellite galaxies which are strongly star forming in relation to that seen in central galaxies. In a future paper, we will directly measure the HI content of star forming satellites for indications that this is being drained, and thus place stronger constraints on the cooling mode.

\section{Conclusions}
We have compiled measurements of the quenching and delay times of satellite galaxies as a function of redshift from the literature. We have used these measurements, in combination with a simple model for the baryon cycle in galaxies, to arrive at the following conclusions:

\begin{itemize}
\item While our compilation of measurements may be subject to systematic effects, the evolution of quenching timescales is consistent with the scaling expected if the quenching was driven not by `orbit-based' stripping events like ram pressure stripping or strangulation but rather by secular outflows in light of a halt of cosmological accretion (overconsumption).  Alternatively, it may be that `orbit-based' stripping is significantly more effective at high redshift than a dynamical scaling would suggest.

\item The delay time of satellite galaxies places an upper limit on the rate of outflow driven gas, or the mass-loading, which from current data is $\eta$ = 2.5 in \Mstel\ = 10$^{10.5}$\Msun\ galaxies. Further measurements at higher redshift and lower masses will put even tighter constraints on the outflow rate.

\item The delay time requires an uninterrupted rate of cooling gas from the gas reservoir. This presents a challenge to typical models of cooling gas for which the cooling rate is strongly dependent on the density of the reservoir, and thus declines quickly as the reservoir drains.

\item We predict that at z$>$1.5 satellite galaxies are quenched through secular outflow processes instead of external stripping events. Future observations of such galaxies will put constraints on the outflow rates of all star forming galaxies and probe this unique quenching mechanism.

\end{itemize}
\section*{Acknowledgments}
SLM acknowledges support from an NWO grant to George Miley. RGB is supported by STFC grant ST/I001573/1. MLB acknowledges support from an NSERC Discovery grant, and from NOVA and NWO visitor grants that supported his sabbatical visit to Leiden Observatory, where this work was undertaken.