\section{Introduction}

In 1949, Lehner~\cite{Lehner1, Lehner2} showed that many of the Fourier coefficients $c(n)$ of the classical $j$-function $j(z) = q^{-1} + 744 + \sum c(n) q^n$, where $q = e^{2\pi i z}$ as usual, are highly divisible by small primes dividing $n$.  For instance, he showed that $c(2^a3^b5^c7^dn) \equiv 0 \pmod{2^{3a+8}3^{2b+3}5^{c+1}7^d}$ for any positive integers $a, b, c, d$.  These divisibility results were generalized to stronger congruences modulo larger powers of these primes by Aas~\cite{Aas} and Kolberg~\cite{Kolberg1, Kolberg2}, and such congruences were shown to hold for all elements of the canonical basis $\{q^{-m} + O(q)\}$ for the space of weakly holomorphic modular functions for ${\text {\rm SL}}_2(\mathbb{Z})$ by Griffin~\cite{Griffin}.

Let $N$ be a positive integer such that the group $\Gamma_0(N)$ has genus zero; thus, $N \in \{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 13, 16, 18, 25\}$.  The space of modular functions of genus zero level $N$ with poles only at the cusp at $\infty$ has a similar canonical basis, and analogous divisibility results have been shown for the Fourier coefficients of the basis elements in the genus zero prime power levels $N \in \{2, 3, 4, 5, 7, 8, 9, 13, 16, 25\}$ in~\cite{Nick,Haddock, JT2,JT}.  In this paper, we prove divisibility results for the Fourier coefficients of canonical basis elements in the remaining genus zero levels $N \in \{6, 10, 12, 18\}$.  Additionally, we show that the coefficients of these basis elements satisfy Zagier duality, and we explicitly construct their generating functions.

Let $M_k^!(N)$ be the space of weakly holomorphic modular forms for $\Gamma_0(N)$ of even integer weight $k$, and let $M_k^\sharp(N)$ be the subspace of such forms with poles allowed only at the cusp at $\infty$.  If $n_0$ is the maximal order of vanishing at the cusp at $\infty$ for a modular form in $M_k^\sharp(N)$, then the space $M_k^\sharp(N)$ has a countably infinite canonical basis $\{f_{k, m}^{(N)}(z)\}_{m \geq -n_0}$ indexed by the order of vanishing $-m$ of $f_{k, m}^{(N)}(z)$ at the cusp at $\infty$.  Specifically, define $f_{k, m}^{(N)}(z)$ to be the unique form in $M_k^\sharp(N)$ with Fourier expansion given by \[f_{k, m}^{(N)}(z) = q^{-m} + \sum_{n > n_0} a_k^{(N)}(m, n) q^n,\] so that the gap between $q^{-m}$ and the next nonzero Fourier coefficient is as large as possible.  These bases were explicitly constructed for prime and prime power levels $N$ of genus zero in~\cite{Duke, Sharon, Haddock, JT2, JT}.

Let $S_k^\sharp(N)$ be the subspace of $M_k^\sharp(N)$ consisting of forms which vanish at all cusps except possibly at $\infty$, and let $n_1$ be the maximal order of vanishing at the cusp at $\infty$ for a modular form in this space.  As we did for $M_k^\sharp(N)$, we may construct a canonical basis $\{g_{k, m}^{(N)}(z)\}_{m \geq -n_1}$ for $S_k^\sharp(N)$, where $g_{k, m}^{(N)}(z)$ is the unique form in $S_k^\sharp(N)$ with Fourier expansion \[g_{k, m}^{(N)}(z) = q^{-m} + \sum_{n > n_1} b_k^{(N)}(m, n) q^n;\] again, the gap between $q^{-m}$ and the next nonzero Fourier coefficient is as large as possible.

The first theorem of this paper is that just as in prime power genus zero levels, the Fourier coefficients of these canonical bases satisfy Zagier duality.
\begin{theorem}\label{thm:duality}
For $N = 6, 10, 12, 18$, all even integers $k$, and all integers $m, n$, we have
\[a_k^{(N)}(m, n) = -b_{2-k}^{(N)}(n, m).\]
\end{theorem}
Additionally, we have the following formula for generating functions for genus zero level basis elements, generalizing the results in~\cite{Duke, Sharon, Haddock, JT2} to all genus zero levels.

\begin{theorem}\label{thm:genfn}
For every genus zero level $N$ we have
\[\sum_{m = -n_0}^\infty f_{k, m}^{(N)}(\tau)q^m = \frac{f_{k, -n_0}^{(N)}(\tau) g_{2-k, n_0+1}^{(N)}(z)}{f_{0, 1}^{(N)}(z) - f_{0, 1}^{(N)}(\tau)}.\]
\end{theorem}

The principal results of this paper give divisibility results for coefficients for weight zero basis elements for levels $6$, $10$, $12$ and $18$ by primes dividing the level.  The following congruences hold in level $6$.
\begin{theorem}\label{thm:level6}
Let $r,s$ be integers relatively prime to $2$. Then
\begin{align}
\label{cong1} &a_0^{(6)}(2^ar,2^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }2^{a-b+2}) & \textnormal{\emph{ if }} a>b,\\
\label{cong2} &a_0^{(6)}(2^ar,2^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }2^{2}) & \textnormal{\emph{ if }} b>a.
\end{align}
Now let $r,s$ be relatively prime to $3$. Then
\begin{align}
&a_0^{(6)}(3^ar,3^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }3^{a-b+1}) & \textnormal{\emph{ if }} a>b,\\
&a_0^{(6)}(3^ar,3^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }3^{1}) & \textnormal{\emph{ if }} b>a.
\end{align}
\end{theorem}

Using the $U_2$ and $U_3$ operators, these congruences in level $6$ lead to very similar congruences for levels $12$ and $18$.
\begin{theorem}\label{thm:level1218}
For $N=12,18$, let $r,s$ be integers relatively prime to $2$. If $N=18$, assume that $3|r$. Then
\begin{align}
&a_0^{(N)}(2^ar,2^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }2^{a-b+2}) & \textnormal{\emph{ if }} a>b,\\
&a_0^{(N)}(2^ar,2^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }2^{2}) & \textnormal{\emph{ if }} b>a.
\end{align}
Now let $r,s$ be relatively prime to $3$. If $N=12$, assume that $2|r$. Then
\begin{align}
&a_0^{(N)}(3^ar,3^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }3^{a-b+1}) & \textnormal{\emph{ if }} a>b,\\
&a_0^{(N)}(3^ar,3^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }3^{1}) & \textnormal{\emph{ if }} b>a.
\end{align}
If $N = 18$ and $3 \nmid r$, or if $N = 12$ and $2 \nmid r$, we have the weaker divisibility results
\begin{align}
&a_0^{(18)}(2^ar,2^bs) \equiv 0 \pmod{2^{a-b}} & \textnormal{\emph{ if }} a > b, \\
\label{weakcong}&a_0^{(12)}(3^ar,3^bs) \equiv 0 \pmod{3^{a-b}} & \textnormal{\emph{ if }} a > b.
\end{align}

\end{theorem}

Similar congruences hold for level $10$, the last genus zero level.
\begin{theorem}\label{thm:level10}
If $r,s$ are relatively prime to $2$, then
\begin{align}
\label{cong9} &a_0^{(10)}(2^ar,2^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }2^{a-b+1}) & \textnormal{\emph{ if }} a>b, \\
\label{cong10} &a_0^{(10)}(2^ar,2^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }2^{1}) & \textnormal{\emph{ if }} b>a.
\end{align}
If $r,s$ are relatively prime to $5$, then
\begin{align}
\label{cong11} &a_0^{(10)}(5^ar,5^bs) \equiv 0 \textnormal{ }(\textnormal{\emph{mod} }5^{a-b}) & \textnormal{\emph{ if }} a>b.
\end{align}
\end{theorem}
One obvious difference in level $10$ is that there is no $b>a$ congruence given modulo powers of $5$.  This behavior is similar to the congruences for Fourier coefficients in level $13$, as described in~\cite{JT}.

We note that because of the duality given by Theorem~\ref{thm:duality}, all divisibility results in weight $0$ immediately imply similar results for the coefficients of the weight $2$ basis elements.  In higher weights, such divisibility results are not expected due to the presence of cusp forms.

In Section~\ref{sec:basis} of this paper, we construct the canonical bases in levels $6, 10, 12$, and $18$.  In Section~\ref{sec:duality}, we prove Theorems~\ref{thm:duality} and~\ref{thm:genfn}.  Sections~\ref{sec:level6}, \ref{sec:level1218}, and \ref{sec:level10} contain the proofs of Theorems~\ref{thm:level6}, \ref{thm:level1218}, and~\ref{thm:level10}.

\section{Canonical Basis Construction}\label{sec:basis}

We now explicitly construct the canonical bases for $M_k^\sharp(N)$ and $S_k^\sharp(N)$ for $N = 6, 10, 12, 18$.  We will give explicit explanations for the construction in level $6$ and more general outlines for the other levels, as all of the constructions are similar.

The congruence subgroup $\Gamma_0(6)$ has four cusps: $\infty$, $0$, $\frac{1}{2}$ and $\frac{1}{3}$. For this level, we choose the Hauptmodul
\[ \psi^{(6)}(z) = f_{0,1}^{(6)}(z) = \frac{\eta(2z)^8\eta(3z)^4}{\eta(z)^4\eta(6z)^8}-4=q^{-1}+6q+4q^2-3q^3+O(q^4) \in M_0^\sharp(6),\]
which has a pole of order $1$ at $\infty$.
We also use the form
\[ f_{2,-2}^{(6)}(z) = \frac{\eta(z)^2\eta(6z)^{12}}{\eta(2z)^4\eta(3z)^6}=q^2-2q^3+3q^4-q^6+7q^8+O(q^9) \in M_2(6),\]
which has order of vanishing $2$ at $\infty$ and has no other zeros or poles.

For nonnegative even integers $k$, the dimension of the space $M_k(6)$ of holomorphic modular forms is $k+1$, so increasing the weight by $2$ increases the dimension by $2$.  Note that multiplying or dividing a form in $M_k^\sharp(6)$ by $f_{2, -2}^{(6)}(z)$ increases or decreases both the form's weight and order of vanishing at $\infty$ by $2$, and can only introduce poles at the cusp at $\infty$.  Thus, in weight $k$ we have the initial basis element $f_{k, -k}^{(6)}(z) = \left(f_{2, -2}^{(6)}(z)\right)^\frac{k}{2}$, which has maximal order of vanishing at $\infty$ for that space. We may recursively construct additional basis elements by multiplying by the Hauptmodul $\psi^{(6)}(z)$ and subtracting earlier basis elements to obtain the largest possible gap between the leading power of $q$ and the next nonzero coefficient.  Therefore, for all $m \geq -k$, the basis element $f_{k, m}^{(6)}(z)$ has the form
\[ f_{k,m}^{(6)}(z)= q^{-m} + \sum_{n = k+1}^\infty a_k^{(6)}(m, n) q^n = \left(f_{2,-2}^{(6)}(z)\right)^\frac{k}{2}P\left(\psi^{(6)}(z)\right), \]
where $P(x)$ is a polynomial in $x$ of degree $m+k$.  We note that because $\psi^{(6)}(z)$ and $f_{2, -2}^{(6)}(z)$ have integral Fourier coefficients, the Fourier coefficients $a_k^{(6)}(m, n)$ of these basis elements and the coefficients of the polynomial $P(x)$ are also all integers.  Additionally, we have $n_0 = k$ for the largest possible order of vanishing at $\infty$ of a form in $M_k^\sharp(6)$.

The congruence subgroup $\Gamma_0(12)$ has six cusps: $\infty$, $0$, $\frac{1}{2}$, $\frac{1}{3}$, $\frac{1}{4}$ and $\frac{1}{6}$. For this level, we choose the Hauptmodul
\[ \psi^{(12)}(z) = f_{0,1}^{(12)}(z)=\frac{\eta(4z)^4\eta(6z)^2}{\eta(2z)^2\eta(12z)^4}=q^{-1}+2q+q^3+O(q^7) \in M_0^\sharp(12), \]
which may be used to generate the basis for $M_0^\sharp(12)$.

Increasing a nonnegative even weight $k$ by $2$ increases the dimension of $M_k(12)$ from $2k$ to $2k+4$.  The form
\[
\begin{split}
f_{2,-4}^{(12)}(z)= &\frac{1}{27}\frac{\eta(z)^{10}\eta(4z)\eta(6z)^9}{\eta(2z)^7\eta(3z)^6\eta(12z)^3}+\frac{11}{72}\frac{\eta(z)^7
\eta(4z)^4\eta(6z)^9}{\eta(2z)^7\eta(3z)^5\eta(12z)^4}-\frac{1}{12}\frac{\eta(z)^4\eta(4z)^7\eta(6z)^9}{\eta(2z)^7
\eta(3z)^4\eta(12z)^5}\\
&+\frac{1}{54}\frac{\eta(z)\eta(4z)^{10}\eta(6z)^9}{\eta(2z)^7\eta(3z)^3\eta(12z)}-\frac{1}{8}\frac{\eta(z)^9
\eta(4z)^3\eta(6z)^2}{\eta(2z)^6\eta(3z)^3\eta(12z)} \\
= &q^4 + O(q^5) \in M_2(12)
\end{split}
\]
has a zero of order $4$ at $\infty$.  We compute that there are no other zeros or poles at any other cusp, and note that as eta-quotients have no poles on the upper half plane, the valence formula implies that there are no other zeros.  Thus, the first basis element of $M_k^\sharp(12)$ must be $\left(f_{2, -4}^{(12)}(z) \right)^\frac{k}{2}$.
Multiplying by $\psi^{(12)}(z)$ and row reducing as before, we construct the basis elements $f_{k, m}^{(12)}(z)$ for all $m \geq -2k$, and note that $n_0 = 2k$ in level $12$.  Again, all Fourier coefficients $a_k^{(12)}(m, n)$ are integers, because $\psi^{(12)}(z)$ and $f_{2, -4}^{(12)}(z)$ have integral Fourier coefficients by the $p$-adic Sturm bound.

The congruence subgroup $\Gamma_0(18)$ has six cusps: $\infty$, $0$, $\frac{1}{2}$, $\frac{1}{3}$, $\frac{1}{6}$ and $\frac{1}{9}$. For this level, we choose the Hauptmodul
\[ \psi^{(18)}(z) = f_{0,1}^{(18)}(z)=\frac{\eta(6z)\eta(9z)^3}{\eta(3z)\eta(18z)^3}\]
and raise or lower the weight by multiplying by the form
\[
\begin{split}
f_{2,-6}^{(18)}(z) = &\frac{25}{216}\frac{\eta(z)^8\eta(6z)^2\eta(9z)^4}{\eta(2z)^4\eta(3z)^4\eta(18z)^2}  -\frac{11}{144}\frac{\eta(z)^3\eta(6z)^8\eta(9z)^7}{\eta(2z)^3\eta(3z)^6\eta(18z)^5}-\frac{121}{972}\frac{\eta(z)^6\eta(6z)^7
\eta(9z)}{\eta(2z)^3\eta(3z)^5\eta(18z)^2}\\
 &- \frac{41}{144}\frac{\eta(z)^6\eta(6z)^2\eta(9z)^6}{\eta(2z)^3\eta(3z)^4\eta(18z)^3}+\frac{67}{144}\frac{\eta(z)^4\eta(6z)^7
 \eta(9z)^3}{\eta(2z)^2\eta(3z)^5\eta(18z)^3}\\
 & +  \frac{1}{972}\frac{\eta(2z)^9\eta(3z)^8\eta(12z)}{\eta(z)^6\eta(6z)^6\eta(9z)^2}-\frac{125}{1296}\frac{\eta(z)\eta(2z)^4
 \eta(9z)^2}{\eta(3z)\eta(6z)\eta(18z)} \\
= &q^6 + O(q^7) \in M_2(18).
\end{split}
\]
Each basis element in $M_k^\sharp(18)$ will be $\left(f_{2, -6}^{(18)}(z)\right)^\frac{k}{2}$ multiplied by a polynomial in $\psi^{(18)}(z)$, and we have $n_0 = 3k$.  Again, all Fourier coefficients $a_k^{(18)}(m, n)$ are integers.

The congruence subgroup $\Gamma_0(10)$ has four cusps: $\infty$, $0$, $\frac{1}{2}$ and $\frac{1}{5}$. For this level, we use the Hauptmodul
\[\psi^{(10)}(z) = f_{0,1}^{(10)}(z)=\frac{\eta(2z)\eta(5z)^5}{\eta(z)\eta(10z)^5}-1.\]

Increasing the weight $k$ by $4$ increases the dimension of $M_k(10)$ by $6$, with the dimension of $M_0(10)$ equal to $1$ and the dimension of $M_2(10)$ equal to $3$.  Writing $k = 4\ell + k'$ with $\ell \in \mathbb{Z}$ and $k' \in \{0, 2\}$, we use the forms
\[ f^{(10)}_{2,-2} = \frac{\eta(10z)^5}{\eta(2z)} = q^3 + O(q^4) \in M_2(10), \]
\[\begin{split}
f^{(10)}_{4,-6} = &\frac{27}{10000}\frac{\eta(2z)^{14}\eta(5z)^8}{\eta(z)^{8}\eta(10z)^{6}}-\frac{13}{4000}\frac{\eta(2z)^{11}\eta(5z)^{11}}
{\eta(z)^7\eta(10z)^7} -\frac{14}{125}\frac{\eta(z)^3\eta(2z)\eta(5z)^9}{\eta(10z)^5}\\
&+\frac{167}{625}\frac{\eta(z)^2\eta(2z)^4\eta(5z)^6}{\eta(10z)^4}+\frac{1}{625}\frac{\eta(z)^{12}\eta(5z)^4}
{\eta(2z)^6\eta(10z)^2}-\frac{7}{32}\frac{\eta(z)\eta(2z)^7\eta(5z)^3}{\eta(10z)^3}+\frac{1}{16}\frac{\eta(2z)^{10}}
{\eta(10z)^2} \\
= &q^6 + O(q^7) \in M_4(10)
\end{split}\]
to see that the first basis element in weight $k = 4\ell + k'$ will be
\[\left(f_{4, -6}^{(10)}(z)\right)^\ell \left(f_{2, -2}^{(10)}(z)\right)^\frac{k'}{2}.\]
Further basis elements are obtained, as before, by multiplying by a polynomial in $\psi^{(10)}(z)$ to get the correct Fourier expansion, with $n_0 = 6\ell + k'$ and $a_k^{(10)}(m, n) \in \mathbb{Z}$.

To construct the basis elements $g_{k, m}^{(N)}(z)$ for the spaces $S_k^\sharp(N)$, we follow the same general process of finding the basis element with the largest leading power of $q$, then multiplying by the Hauptmodul and row reducing with previous basis elements.  The first basis element $g_{k, -n_1}^{(N)}(z)$ can be constructed as
\[ g_{k, -n_1}^{(N)}(z) = f_{k, -n_0}^{(N)}(z) \prod_{\frac{a}{b}} \left(\psi^{(N)}(z) - c_{\frac{a}{b}}\right),\]
where $f_{k, -n_0}^{(N)}(z)$ is the first basis element for $M_k^\sharp(N)$, the product is over all non-$\infty$ cusps of $\Gamma_0(N)$, and $c_{\frac{a}{b}}$ is the value of $\psi^{(N)}(z)$ at the cusp $\frac{a}{b}$.  It is clear that this will be a form of the correct weight which vanishes at all non-$\infty$ cusps, and it must have maximal order of vanishing at $\infty$ in $S_k^\sharp(N)$ because $f_{k, -n_0}^{(N)}(z)$ does in $M_k^\sharp(N)$, and the existence of a Hauptmodul for $\Gamma_0(N)$ means that $f_{k, -n_0}^{(N)}(z)$ cannot vanish at any cusp other than $\infty$.  Therefore, we have $n_1 = n_0 - c(N) + 1$, where $c(N)$ is the number of cusps of $\Gamma_0(N)$.

In level $6$, the Hauptmodul has the values $0$, $1$, $9$ at the three non-infinity cusps. This may be computed using Ligozat's theorem (see~\cite{Lig}) to find that the form
\[g_{2,1}^{(6)}(z) = \frac{\eta(2z)^6\eta(3z)^8}{\eta(6z)^{10}} \in M_2^\sharp(6)\]
vanishes at all non-$\infty$ cusps, and then writing the quotient $g_{2, 1}^{(6)}(z)/f_{2, -2}^{(6)}(z)$ as a polynomial in $\psi^{(6)}(z)$. The roots of the polynomial are the values at the cusps.  Thus, the first basis element of $S_k^\sharp(6)$ will be
\[g_{k, -k+3}^{(6)}(z) = \left(f_{2, -2}^{(6)}(z)\right)^{\frac{k}{2}} \psi^{(6)}(z) (\psi^{(6)}(z)-1) (\psi^{(6)}(z) - 9).\]

In level 12, the polynomial \[(\psi^{(12)}(z))(\psi^{(12)}(z)-1)(\psi^{(12)}(z)+1)(\psi^{(12)}(z)-3)(\psi^{(12)}(z)+3)\] vanishes at all non-$\infty$ cusps, while levels 18 and 10 have the polynomials
\[ (\psi^{(18)}(z))^7-7(\psi^{(18)}(z))^4-8(\psi^{(18)}(z)),\]
\[(\psi^{(10)}(z)+2)(\psi^{(10)}(z)+1)(\psi^{(10)}(z)-3),\]
which vanish at all cusps away from $\infty$.

\section{Duality and Generating Functions}\label{sec:duality}
The proof of Zagier duality for the coefficients of these canonical bases (Theorem~\ref{thm:duality}) follows the argument in~\cite{JT2} for prime power levels, which is itself a generalization of an argument credited to Kaneko by Zagier~\cite{Zagier}.

Recall that Ramanujan's $\theta$-operator acts on a weakly holomorphic modular form $f$ by $\theta f = q \frac{d}{dq} f = \frac{1}{2\pi i}\frac{d}{dz} f$ and increases the weight by $2$, preserving modularity if $f$ is weight $0$.  With the usual slash operator $f|_k \gamma = (\det \gamma)^{k/2} (cz+d)^{-k}f\left(\frac{az+b}{cz+d}\right)$ for a matrix $\gamma = \smat{a}{b}{c}{d}$, the $\theta$-operator satisfies the intertwining relation \[\theta(f|_0 \gamma) = (\theta f)|_2 \gamma.\]
For a form $f \in M_0^\sharp(N)$, we know that $f |_0 \gamma$ is holomorphic for any matrix $\gamma \in {\text {\rm SL}}_2(\mathbb{Z})$ and not in $\Gamma_0(N)$, so the intertwining relation and the fact that derivatives annihilate constants imply that $\theta$ actually maps forms in $M_0^\sharp(N)$ to forms in $S_2^\sharp(N)$.  Specifically, by examining Fourier expansions we find that
\begin{equation}\label{lemma:theta}
\theta\left(f_{0,m}^{(N)}(z)\right) =-mg_{2,m}^{(N)}(z).
\end{equation}
In weight $2$, the maximal order of vanishing $n_1$ for a form in $S_2^\sharp(N)$ is equal to $-1$ for $N = 6, 10, 12, 18$, so the functions $\theta\left(f_{0,m}^{(N)}\right) =-mg_{2,m}^{(N)}$ span the space $S_2^\sharp(N)$, and every form in $S_2^\sharp(N)$ must have no constant term in its Fourier expansion.

Now consider the product $f_{k,m}^{(N)}g_{2-k,m}^{(N)}$, which must be in the space $S_2^\sharp(N)$.  By examining the Fourier expansions of the forms in the product, we find that its constant term is $a_k^{(N)}(m,n)+b_{2-k}^{(N)}(n,m)$, which must be zero, and Theorem~\ref{thm:duality} follows.

We next prove Theorem~\ref{thm:genfn}; the proof is similar to that in~\cite{ElG}.  Fix an even integer weight $k$ and a genus zero level $N$.  Let $n_0, N_1$ be the maximal orders of vanishing at $\infty$ for forms in $M_k^\sharp(N)$, $S_{2-k}^\sharp(N)$ respectively.  Thus, we have \[f_{k, m}^{(N)}(z) = q^{-m} + \sum_{n > n_0} a_k(m, n)q^n\] for every integer $m \geq -n_0$, and
\[g_{2-k, m}^{(N)}(z) = q^{-m} + \sum_{n > N_1} b_{2-k}(m, n) q^n\] for every integer $m \geq -N_1$.
Let \[F_k(z, \tau) = \sum_{m = -n_0}^\infty f_{k, m}^{(N)}(\tau) q^m\] be the generating function for the $f_{k, m}^{(N)}(z)$.  In what follows, we often write $f_{k,m}(z)$ for $f_{k, m}^{(N)}(z)$ for simplicity.

Multiply the Fourier expansions for $f_{0, 1}(z)$ and $f_{k,m}(z)$ and write the result, which is a form in $M_k^\sharp(N)$, in terms of basis elements, remembering that the powers of $q$ less than or equal to $n_0$ completely determine the form.  Doing this, we obtain
\begin{align*}
f_{0, 1}(z) f_{k,m}(z) = &\left(q^{-1} + \sum_{n=1}^\infty a_0(1, n)q^n\right)\left(q^{-m} + \sum_{n>n_0}a_k(m, n)q^n\right) \\
= &q^{-m-1} + \sum_{n=1}^\infty a_0(1, n)q^{n-m} + \sum_{n > n_0} a_k(m, n)q^{n-1} + \sum_{n=1}^\infty \sum_{N > n_0} a_0(1, n) a_k(m, N) q^{N+n} \\
= &f_{k, m+1}(z) + a_k(m, n_0+1)f_{k, -n_0}(z) + a_0(1, n_0+m)f_{k, -n_0}(z) \\
  &+\sum_{n=1}^{n_0+m-1} a_0(1, n) f_{k, m-n}(z) \hfill
\end{align*}
for any $m \geq -n_0$, where the last sum appears only if $m \geq -n_0+2$ and the second to last term appears only if $m \geq -n_0+1$.  This gives us a recurrence relation for $f_{k, m+1}(z)$ in terms of forms $f_{k, i}(z)$ with $i \leq m$.

Using this recurrence, we may write
\begin{equation*}
\begin{split}
F_k(z, \tau) = &f_{k, -n_0}(\tau)q^{-n_0} + \sum_{m = -n_0+1}^\infty f_{k, m}(\tau) q^m = f_{k,-n_0}(z) q^{-n_0} + \sum_{m = -n_0}^\infty f_{k, m+1}(z) q^{m+1} \\
= &f_{k, -n_0}(\tau)q^{-n_0} + \sum_{m=-n_0}^\infty (f_{0, 1}(\tau) f_{k,m}(\tau)q^{m+1} - a_k(m, n_0+1)f_{k, -n_0}(\tau)q^{m+1} \\
& - a_0(1, n_0+m)f_{k, -n_0}(\tau)q^{m+1} - \sum_{n=1}^{n_0+m-1} a_0(1, n)f_{k, m-n}(\tau)q^{m+1}).
\end{split}
\end{equation*}

Using the duality $a_k(m, n_0+1) = -b_{2-k}(n_0+1, m)$, we note that the term \[-\sum_{m=-n_0}^\infty
a_k(m,n_0+1)f_{k, -n_0}(\tau)q^{m+1}\] may, if $N_1\geq -n_0-1$, be rewritten as $q \cdot f_{k, -n_0}(\tau) (g_{2-k, n_0+1}(z) - q^{-n_0-1})$.  We note also that this inequality is true (and is, in fact, an equality) for all genus zero levels.  Reversing the order of the double summation, we rewrite other terms similarly to obtain
\begin{align*}
F_k(z, \tau) = &f_{k, -n_0}(\tau)q^{-n_0} + q f_{0, 1}(\tau)F_k(z, \tau) \\
               &+f_{k, -n_0}(\tau) q (g_{2-k, n_0+1}(z) - q^{-n_0-1}) - f_{k, -n_0}(\tau) q^{1-n_0} (f_{0, 1}(z)-q^{-1}) \\
               &-q(f_{0, 1}(z) - q^{-1})(F_k(z, \tau) - f_{k, -n_0}(\tau) q^{-n_0}),
\end{align*}
which simplifies to
\[0 = qf_{0, 1}(\tau)F_k(z, \tau) + q f_{k, -n_0}(\tau)g_{2-k, n_0+1}(z)-qf_{0, 1}(z) F_k(z, \tau).\]
Therefore, we have the generating function
\[F_k(z, \tau) = \frac{f_{k, -n_0}(\tau) g_{2-k, n_0+1}(z)}{f_{0, 1}(z) - f_{0, 1}(\tau)}.\]

\section{Level 6 Congruences}\label{sec:level6}

We will now prove Theorem~\ref{thm:level6} by first proving the congruence~(\ref{cong2}), and then showing that~(\ref{cong2}) implies~(\ref{cong1}).

We use the usual slash operator defined at the beginning of section~\ref{sec:duality}, omitting the subscript $k$ for brevity.  If $f \in M_0^\sharp(6)$ and $\gamma \in {\text {\rm SL}}_2(\mathbb{Z})$, this becomes $f|\gamma = f\left(\frac{az+b}{cz+d}\right)$.
Recall that if $p | N$, the $U_p$ operator
\[U_p f = \frac{1}{p} \sum_{i=0}^{p-1} f|\smat{1}{i}{0}{p}\]
sends forms in $M_k^!(N)$ to forms in $M_k^!(N)$, while if $p^2|N$, the $U_p$ operator sends forms in $M_k^!(N)$ to forms in $M_k^!(\frac{N}{p})$.  Its action on Fourier expansions is given by $U_p \sum a(n) q^n = \sum a(pn)q^n$.
Thus, the $U_2$ operator on forms in $M_0^\sharp(6)$ can be written as
\[2U_2f(z) = f|\smat{1}{0}{0}{2} + f|\smat{1}{1}{0}{2}.\]

We are interested in the order of vanishing of $U_2f(z)$ at the other cusps of $\Gamma_0(6)$, which are $0, \frac{1}{2}$, and $\frac{1}{3}$.  Computing, we find that
\[2U_2f|\smat{0}{-1}{1}{0} = f|\smat{0}{-1}{1}{0}|\smat{2}{0}{0}{1} + f|\smat{1}{0}{2}{1}|\smat{1}{-1}{0}{2}, \]
\[2U_2f|\smat{1}{0}{2}{1} = f|\smat{1}{0}{4}{1}|\smat{1}{0}{0}{2} + f|\smat{3}{-1}{4}{-1}|\smat{1}{1}{0}{2}, \]
\[2U_2f|\smat{1}{0}{3}{1} = f|\smat{1}{0}{6}{1}|\smat{1}{0}{0}{2} + f|\smat{2}{-1}{3}{-1}|\smat{2}{1}{0}{1}.\]
Under the action of $\Gamma_0(6)$, we have $\frac{1}{4}$ equivalent to $\frac{1}{2}$, and $\frac{3}{4}$ also equivalent to $\frac{1}{2}$.  Since replacing $z$ with $\frac{az+b}{d}$ does not change holomorphicity at a cusp, it follows that if $f \in M_0^!(6)$ is holomorphic at $0$ and $\frac{1}{2}$, then $U_2 f$ is as well.  Therefore, applying $U_2$ repeatedly to a form in $M_0^\sharp(6)$ will give a form in $M_0^!(6)$ with poles that may occur only at the cusps at $\infty$ and $\frac{1}{3}$.

We note that for a form in $M_0^!(6)$, slashing with the matrix $\smat{2}{-1}{6}{-2}$ is an Atkin-Lehner involution which takes a form in $M_0^!(6)$ to a form in $M_0^!(6)$.  Computing as before, we find that
\[f|\smat{2}{-1}{6}{-2}|\smat{0}{-1}{1}{0} = f|\smat{-1}{-2}{-2}{-6} = f|\smat{-1}{0}{-2}{-1}|\smat{1}{2}{0}{2},\]
\[f|\smat{2}{-1}{6}{-2}|\smat{1}{0}{2}{1} =  f|\smat{0}{-1}{2}{-2} = f|\smat{0}{1}{-1}{1}|\smat{-2}{1}{0}{-1},\]
\[f|\smat{2}{-1}{6}{-2}|\smat{1}{0}{3}{1} = f|\smat{-1}{-1}{0}{-2}.\]
These calculations show that if $f$ is holomorphic at $\infty$, then $f|\smat{2}{-1}{6}{-2}$ is holomorphic at $\frac{1}{3}$; similarly, if $f$ is holomorphic at $0, \frac{1}{2}$, or $\frac{1}{3}$, then $f|\smat{2}{-1}{6}{-2}$ is holomorphic at $\frac{1}{2}, 0$, or $\infty$ respectively.

We now show that if $b > a$ and $r, s$ are odd integers, then the coefficient $a_0^{(6)}(2^a r, 2^b s)$ is divisible by $4$, proving equation~(\ref{cong2}).  Suppose first that $a = 0$. Applying the $U_2$ operator to the form $f_{0, r}^{(6)}(z) = q^{-r} + \sum_{n \geq 1} a_0^{(6)}(r, n) q^n$ gives us a modular form in $M_0^!(6)$ with Fourier expansion $\sum_{n \geq 1} a_0^{(6)}(r, 2n) q^n$ which is clearly holomorphic at $\infty$; from the discussion above, it must also be holomorphic at $0$ and $\frac{1}{2}$, since $f_{0, r}^{(6)}$ is as well.  Thus, $U_2 f_{0, r}^{(6)}(z)$ has poles only at $\frac{1}{3}$.

Replacing $z$ with $\frac{2z-1}{6z-2}$ in the formula for $U_2 f(z)$ above, we compute that
\[2U_2 f_{0, r}^{(6)}\left(\frac{2z-1}{6z-2}\right) = f_{0, r}^{(6)}\left(z-\frac{1}{2}\right) + f_{0, r}^{(6)}\left(\frac{8z-3}{12z-4}\right).\] The first term gives $-q^{-r} + O(1)$ in the Fourier expansion since $r$ is odd. The second term is holomorphic at $\infty$ because $\smat{8}{-3}{12}{-4} = \smat{2}{-3}{3}{-4} \smat{4}{0}{0}{1}$ and $f_{0, r}^{(6)}|\gamma$ is holomorphic at $\infty$ for any $\gamma \in {\text {\rm SL}}_2(\mathbb{Z})$ not in $\Gamma_0(6)$.  Slashing with $\smat{0}{-1}{1}{0}, \smat{1}{0}{2}{1}, \smat{1}{0}{3}{1}$ as above, we find that this form is holomorphic at 0 and $\frac{1}{2}$, and that at the cusp at $\frac{1}{3}$ it is equal to $f_{0, r}^{(6)}(\frac{z+1}{4}) + f_{0, r}^{(6)}(\frac{z+3}{4})$.  Replacing $z$ by $\frac{z+1}{4}$ and $\frac{z+3}{4}$ in the Fourier expansion of $f_{0, r}^{(6)}(z)$, we find that the negative powers of $q$ cancel because $r$ is odd, and conclude that the form is also holomorphic at $\frac{1}{3}$.  Since its Fourier expansion begins $-q^{-r}$, it must be $-f_{0, r}^{(6)}(z)$.

Let $\psi(z)$ be the Hauptmodul $\psi(z) = \frac{\eta(z)^5 \eta(3z)}{\eta(2z)\eta(6z)^5} \in M_0^\sharp(6)$; we note that $\psi(z)$ has a pole of order $1$ at $\infty$ and vanishes at the cusp at $0$.  From the construction of the basis elements $f_{0, m}^{(6)}(z)$, we may write each basis element as an integral polynomial in powers of $\psi(z)$, giving us
\[2U_2 f_{0, r}^{(6)}\left(\frac{2z-1}{6z-2}\right) = -f_{0, r}^{(6)}(z) = \sum_{i=0}^r c_i \psi^i(z),\]
with the coefficients $c_i \in \mathbb{Z}$.

Using the equalities $\eta(z+1) = e^{2\pi i/24} \eta(z)$ and $\eta(-1/z) = \sqrt{z/i}\eta(z)$, we compute that
\[\psi|\smat{2}{-1}{6}{2} = \frac{-2^3 \eta(2z)^5\eta(6z)}{\eta(3z)^5\eta(z)}.\]  We note that $\psi_{\frac{1}{3}}(z) = \frac{\eta(2z)^5\eta(6z)}{\eta(3z)^5\eta(z)}$ is a modular form in $M_0^!(6)$ with integer coefficients; it is holomorphic at $\infty$ and $0$, vanishes at $\frac{1}{2}$, and has a pole of order 1 at $\frac{1}{3}$.

Replacing $z$ with $\frac{2z-1}{6z-2}$ in our expression for $2U_2 f_{0, r}^{(6)}$ above, we find that
\[2U_2 f_{0, r}^{(6)}(z) = \sum_{i=0}^r c_i 2^{3i} \psi_{\frac{1}{3}}(z)^i.\]  Because the $c_i$ and all Fourier coefficients of $\psi_{\frac{1}{3}}$ are integers, we conclude that for $n \geq 1$, all of the Fourier coefficients $a_0^{(6)}(r, 2n)$ of $U_2 f_{0, r}^{(6)}(z)$ must be divisible by 4.

Now suppose that $a \geq 1$ and that the coefficients $a_0^{(6)}(2^\alpha r', 2^\beta s')$ are divisible by 4 for every nonnegative integer $\alpha < a$ and every $\beta > \alpha$ and odd $r', s'$.  We will show that $a_0^{(6)}(2^a r, 2^b s)$ is divisible by $4$ for every integer $b > a$, proving the theorem by induction.

Define $G(z)$ to be the form $U_2 f_{0, 2^a r}^{(6)}(z) - f_{0, 2^{a-1} r}^{(6)}(z)$.  We note that $G(z)$ is holomorphic at $\infty$ because the $q^{2^{a-1}r}$ terms cancel, that it is holomorphic at $0$ and $\frac{1}{2}$ and meromorphic at $\frac{1}{3}$, and that its $n$th Fourier coefficient is $a_0^{(6)}(2^a r, 2n) - a_0^{(6)}(2^{a-1}r, n)$.  Replacing $z$ by $\frac{2z-1}{6z-2}$ as before, we have the form
\[2 G\left(\frac{2z-1}{6z-2}\right) = f_{0, 2^a r}^{(6)}\left(z-\frac{1}{2}\right) + f_{0, 2^a r}^{(6)}\left(\frac{8z-3}{12z-4}\right) - 2 f_{0, 2^{a-1} r}^{(6)}\left(\frac{2z-1}{6z-2}\right),\]
which is holomorphic away from $\infty$.
The second two terms are holomorphic at $\infty$, and the first term gives a $q^{-2^a r}$, so this must be $f_{0, 2^a r}(z)$, which has integer coefficients.  Thus, we have
\[2 G\left(\frac{2z-1}{6z-2}\right) = \sum_{i=0}^{2^a r} c_i \psi(z)^i,\] with the coefficients $c_i \in \mathbb{Z}$.  Replacing $z$ with $\frac{2z-1}{6z-2}$ as before, we get
\begin{equation}2G(z) = \sum_{i=0}^{2^a r}c_i 2^{3i-1} \psi_{\frac{1}{3}}(z)^i. \label{sumeqn} \end{equation}
Applying the $U_2$ operator $a$ times to this equation, we get a form $U_2^a G(z) \in M_0^!(6)$ with poles only at the cusp at $\frac{1}{3}$ and with $n$th Fourier coefficient given by
\[a_0^{(6)}(2^a r, 2^{a+1} n) - a_0^{(6)}(2^{a-1}r, 2^a n).\]  By our inductive hypothesis, the second term in this expression is divisible by 4.  Since each positive power of $q^n$ in~\eqref{sumeqn} has a coefficient divisible by $4$ as well, we conclude that $a_0^{(6)}(2^a r, 2^{a+1} n)$ is divisible by $4$, proving congruence~(\ref{cong2}).

Now consider the coefficient $a_0(2^a r,2^b s)$ with $a > b$. By duality and equation~(\ref{lemma:theta}), we have \[  a_0(2^a r,2^b s) = \frac{-2^a r}{2^b s}b_2(2^a r,2^b s) = \frac{2^a r}{2^b s}a_0(2^b r,2^a s).\]
Note that $\frac{r}{s}a_0(2^b r, 2^a s)$ must be an integer because all basis elements have integer coefficients in level 6. Then by the congruence just proved, since $a>b$, we have $2^2|a_0(2^b r, 2^a s)$, and it follows that $2^{a-b+2}|a_0(2^a r, 2^b s)$.

The proof of congruences for powers of 3 in level $6$ is very similar; we will note explicit calculations that differ.
The $U_3$ operator on forms $f \in M_0^\sharp(6)$ can be written as
\[3U_3f(z) = f|\smat{1}{0}{0}{3} + f|\smat{1}{1}{0}{3} + f|\smat{1}{2}{0}{3}.\]
To check the order of vanishing of $U_3f(z)$ at the cusps $0$, $\frac{1}{2}$ and $\frac{1}{3}$ of $\Gamma_0(6)$, we compute that
\[3U_3f|\smat{0}{-1}{1}{0} = f|\smat{0}{-1}{1}{0}|\smat{3}{0}{0}{1} + f|\smat{1}{0}{3}{1}|\smat{1}{-1}{0}{3} + f|\smat{2}{-1}{3}{-1}|\smat{1}{1}{0}{3},\]
\[3U_3f|\smat{1}{0}{3}{1} = f|\smat{1}{0}{9}{1}|\smat{1}{0}{0}{3} + f|\smat{4}{-1}{9}{-2}|\smat{1}{1}{0}{3} + f|\smat{7}{3}{9}{4}|\smat{1}{-1}{0}{3},\]
\[3U_3f|\smat{1}{0}{2}{1} = f|\smat{1}{0}{6}{1}|\smat{1}{0}{0}{3} + f|\smat{1}{0}{2}{1}|\smat{3}{1}{0}{1} + f|\smat{5}{-1}{6}{-1}|\smat{1}{1}{0}{3}.\]

Since the cusps at $\frac{2}{3}$, $\frac{1}{9}$, $\frac{4}{9}$, and $\frac{7}{9}$ are equivalent to the cusp at $\frac{1}{3}$, if $f(z)$ is holomorphic at $0$ and $\frac{1}{3}$, then $U_3f(z)$ will be also. Thus, applying $U_3$ repeatedly to a form in $M_0^\sharp(6)$ gives a form with poles only at $\frac{1}{2}$ and $\infty$.

We again use an Atkin-Lehner involution $\smat{3}{1}{6}{3}$. Computing as before, we find that
\[f|\smat{3}{1}{6}{3}|\smat{0}{-1}{1}{0} = f|\smat{1}{-3}{3}{-6} = f|\smat{1}{0}{3}{1}|\smat{1}{-3}{0}{3},\]
\[f|\smat{3}{1}{6}{3}|\smat{1}{0}{2}{1} =  f|\smat{5}{1}{12}{3} = f|\smat{5}{2}{12}{5}|\smat{1}{-1}{0}{3} = f|\smat{1}{-1}{0}{3},\]
\[f|\smat{3}{1}{6}{3}|\smat{1}{0}{3}{1} = f|\smat{6}{1}{15}{3} = f|\smat{2}{3}{5}{8}|\smat{3}{-1}{0}{1}\]
for any $f \in M_0^!(6)$.

Suppose now that $b > a$ and that $r, s$ are not divisible by $3$.  In the $a=0$ case, replacing $z$ with $\smat{3}{1}{6}{3}$ in the formula for $U_3f(z)$, we compute that
\[3U_3 f_{0, r}^{(6)}\left(\frac{3z+1}{6z+3}\right) = f_{0, r}^{(6)}\left(z+\frac{1}{3}\right) + f_{0, r}^{(6)}\left(\frac{9z+4}{18z+9}\right) + f_{0, r}^{(6)}\left(z+\frac{2}{3}\right).\]
The second term is holomorphic at $\infty$ since $\smat{9}{4}{18}{9} = \smat{1}{2}{4}{9} \smat{9}{0}{0}{1}$. The first and last terms give $-q^{-r}+O(1)$ because
$r\not\equiv 0\ (\textrm{mod}\ 3) $. Slashing with $\smat{0}{-1}{1}{0}$ and $\smat{1}{0}{3}{1}$, we find that this form is holomorphic at $0$ and $\frac{1}{3}$. Slashing with $\smat{1}{0}{2}{1}$ gives $f_{0, r}^{(6)}(\frac{z+2}{9})+f_{0, r}^{(6)}(\frac{z+8}{9}) + f_{0, r}^{(6)}(\frac{z+5}{9})$, and all negative powers of $q$ cancel out since $r\not\equiv 0 \pmod{3}$. We conclude that $3U_3 f_{0, r}^{(6)}\left(\frac{3z+1}{6z+3}\right) =-f^{(6)}_{0,r}(z)$.

Now let $\psi(z)$ be the level $6$ Hauptmodul $\psi(z) = \frac{\eta(z)^5\eta(3z)}{\eta(2z)\eta(6z)^5}$. Because powers of $\psi$ form a basis for $M_0^\sharp(6)$, we can write
\[3U_3 f_{0, r}^{(6)}\left(\frac{3z+1}{6z+3}\right) = -f_{0, r}^{(10)}(z) = \sum_{i=0}^r c_i \psi(z)^i,\]
with the coefficients $c_i \in \mathbb{Z}$. We compute that
\[\psi|\smat{3}{1}{6}{3} = \frac{3 \eta(3z)^5\eta(z)}{\eta(6z)\eta(2z)^5}.\]
Note that $\psi_{\frac{1}{2}}(z) = \frac{\eta(3z)^5\eta(z)}{\eta(6z)\eta(2z)^5}$ is a modular form in $M_0^!(6)$ with integer coefficients that is holomorphic at $\infty$ and $0$, vanishes at $\frac{1}{3}$, and has a pole of order $1$ at $\frac{1}{2}$.

We then use near-identical arguments for the base case and inductive step to show that the stated congruences hold for powers of $3$.

\section{Congruences in levels $12$ and $18$}\label{sec:level1218}
To extend the level $6$ congruences to levels $12$ and $18$, we first relate basis elements in these levels to basis elements in level $6$, using the usual $U_p$ and $V_p$ operators with $p=2, 3$.  Recall that $V_p(f(z)) = f(pz)$, and that if $f \in M_k^!(N)$, then $V_p f \in M_k^!(pN)$.

\begin{lemma}\label{lemma:U_oper}
For any nonnegative integer $m$, we have
\[ U_2\left(f_{0,m}^{(12)}\right)=f_{0,\frac{m}{2}}^{(6)} \textnormal{ \emph{if} }2|m,  \]
\[ U_3\left(f_{0,m}^{(18)}\right)=f_{0,\frac{m}{3}}^{(6)} \textnormal{ \emph{if} }3|m. \]
If $m$ is odd, then $U_2(f_{0, m}^{(12)}) = 0$, and if $3 \nmid m$, then $U_3(f_{0, m}^{(18)}) = 0$.
\end{lemma}
To see this for $p=2$, let $f_{0,2m}\in M^\sharp_0(12)$. We compute that $V_2U_2$ preserves the space $M_0^\sharp(12)$ as in section 3 of \cite{Haddock}.
The first term of $f_{0, 2m}$ will be $q^{-2m}$, as will the first term of $V_2 f_{0,m}^{(6)}$. Then in the expression
\[ V_2\left(U_2\left(f_{0,2m}^{(12)}(z)\right)-f_{0,m}^{(6)}(z)\right),\]
the first terms will cancel and the resulting form will have only nonnegative powers of $q$.  Since it is in the space $M^\sharp_0(12)$, it must be identically zero.
If $2$ does not divide $m$, then $U_2 f_{0, m}^{(12)}$ will be a form in $M_0^\sharp(6)$ which vanishes at $\infty$, and must be zero; thus, when $m$ is odd, the form $f_{0, m}^{(12)}(z)$ will have only odd powers of $q$ in its Fourier expansion.  The proof of the case for level $18$ is similar.

Using Lemma~\ref{lemma:U_oper}, we can prove the congruences for levels $12$ and $18$.
Consider $f_{0,m}^{(12)}\in M_0^\sharp(12)$,
and assume that $m = 2^ar$ for some $a\geq 0$ and odd $r$. Then by Lemma~\ref{lemma:U_oper}, we have $U_2(f_{0,m}^{(12)}) = f_{0,2^{a-1}r}^{(6)}$ if $a \geq 1$ and $U_2(f_{0,m}^{(12)})=0$ if $a = 0$. If $a \geq 1$, then for some $b \geq 0$ and some odd $s$, we have the relation
\[ a_0^{(12)}(2^ar,2^bs) = a_0^{(6)}(2^{a-1}r,2^{b-1}s).\]
Note that if $b=0$, then $a_0^{(12)}(2^ar,2^bs) = 0$.
Suppose that $a>b$. If $b=0$, then the congruence holds vacuously. If $b>0$, then $a-1 > b-1$ and the congruence in level 12 follows from the congruence in level $6$. If instead $b>a$, then if $a=0$, we have $U_2 f_{0, 2^a r}^{(12)} = 0$ and the congruence vacuously holds. Otherwise, $b-1>a-1$ and the congruence again follows from the congruence in level $6$.

For the congruences modulo powers of $3$, write $m=2^a 3^c r$, with $a \geq 1$ and $r'$ not divisible by $2$ or $3$.  Using the $U_2$ operator again, we obtain the relation
\[ a_0^{(12)}(2^a3^cr,2^b3^ds) = a_0^{(6)}(2^{a-1}3^cr,2^{b-1}3^ds')\]
for appropriate $b, d, s'$ with $s'$ not divisible by $2$ and $3$. Again, the congruences follow directly from the congruences in level $6$.  The weaker congruence~(\ref{weakcong}) modulo powers of $3$ when $2$ does not divide the order of the pole can be proved in a manner similar to the way in which~(\ref{cong2}) implies~(\ref{cong1}), using duality, equation~(\ref{lemma:theta}), and the fact that all coefficients $a_0^{(12)}(m, n)$ are integers.

The proof for the level 18 congruences is similar.

\section{Congruences in Level 10}\label{sec:level10}
The proof of congruences~(\ref{cong9}) and~(\ref{cong10}) for powers of $2$ in level $10$ is very similar to the proof for level $6$ in section~\ref{sec:level6}. We will make note of explicit calculations that differ.

Let $f \in M_0^\sharp(10)$. Again, we want to find the behavior at the cusps under the $U_2$ operator.
Level $10$ has the three non-infinity cusps $0$, $\frac{1}{2}$, and $\frac{1}{5}$. The calculations showing order of vanishing for the cusps at $0$ and $\frac{1}{2}$ are similar to the level $6$ calculations. For the cusp at $\frac{1}{5}$, we find that
\[2U_2f|\smat{1}{0}{5}{1} = f|\smat{1}{0}{10}{1}|\smat{1}{0}{0}{2} + f|\smat{3}{-2}{5}{-3}|\smat{2}{1}{0}{1}. \]
It follows that applying $U_2$ repeatedly to a form in $M_0^\sharp(10)$ gives forms with poles only at $\frac{1}{5}$ and $\infty$.

For the space $M_0^!(10)$, the matrix $\smat{4}{-1}{10}{-2}$ is an Atkin-Lehner involution. Computing as before,
\[f|\smat{4}{-1}{10}{-2}|\smat{0}{-1}{1}{0} = f|\smat{-1}{-4}{-2}{-10} = f|\smat{1}{0}{2}{1}|\smat{-1}{-4}{0}{-2},\]
\[f|\smat{4}{-1}{10}{-2}|\smat{1}{0}{2}{1} =  f|\smat{2}{-1}{6}{-2} = f|\smat{1}{0}{3}{1}|\smat{2}{-1}{0}{-1},\]
\[f|\smat{4}{-1}{10}{-2}|\smat{1}{0}{5}{1} = f|\smat{-1}{-1}{0}{-2}.\]

Now if we replace $z$ with $\smat{4}{-1}{10}{-2}$ in the formula for $U_2f(z)$, we compute that
\[2U_2 f_{0, r}^{(10)}\left(\frac{4z-1}{10z-2}\right) = f_{0, r}^{(10)}\left(\frac{4z-1}{20z-4}\right) + f_{0, r}^{(10)}\left(z-\frac{1}{2}\right).\]
Note that the first term here is holomorphic at $\infty$ since $\smat{4}{-1}{20}{-4} = \smat{1}{-1}{5}{-4} \smat{4}{0}{0}{1}$ and that the second term will give $-q^{-r}+O(1)$.  After checking that this is holomorphic at the other cusps, we conclude that $2U_2 f_{0, r}^{(10)}\left(\frac{4z-1}{10z-2}\right) =-f^{(10)}_{0,r}(z)$.

Now let $\psi(z)$ be the level $10$ Hauptmodul $\psi(z) = \frac{\eta(2z)\eta(5z)^5}{\eta(z)\eta(10z)^5}$. By the construction of basis elements, we can write
\[2U_2 f_{0, r}^{(10)}\left(\frac{4z-1}{10z-2}\right) = -f_{0, r}^{(10)}(z) = \sum_{i=0}^r c_i \psi^i(z),\]
with the coefficients $c_i \in \mathbb{Z}$. We compute that
\[\psi|\smat{4}{-1}{10}{2} = \frac{-2^2 \eta(z)\eta(10z)^5}{\eta(2z)\eta(5z)^5}.\]
Note that $\psi_{\frac{1}{5}}(z) = \frac{\eta(z)\eta(10z)^5}{\eta(2z)\eta(5z)^5} $ is a modular form in $M_0^!(10)$ with integer coefficients that is holomorphic at $\infty$, vanishes at $\frac{1}{2}$ and $0$, and has a pole of order $1$ at $\frac{1}{5}$.

After these explicit calculations, we can use near-identical arguments for the base case and inductive step to show that congruence~(\ref{cong10}) holds.  Congruence~(\ref{cong9}) follows as before using duality and equation~(\ref{lemma:theta}); the congruence~(\ref{cong11}) modulo powers of $5$ is proved in the same way.

\bibliographystyle{amsplain}


\section{Introduction}
Fossil systems are classically thought to be representative of old, undisturbed galaxy systems where almost all $L^*$ members have been cannibalized by the dominant central elliptical, as dynamical friction draws the massive $L^*$ member galaxies to the center over long time scales.  As the central elliptical cannibalizes more galaxies, the magnitude gap between the central elliptical and the next most massive member widens and mass becomes more concentrated at the center until a `fosssil system' is created (Jones et al.\ 2003).  One study of the mass concentration of fossil groups using N-body simulations show some support for this assumption toward fossil systems formation  (Khosroshahi et al.\ 2007).  Due to the apparent evolved nature of fossil systems and high concentration parameters (Wechsler et al. 2002), it is possible that studying fossil systems can help us understand of properties of the early universe as well as brightest cluster/group galaxy formation.  However, exactly how common fossil systems are in the universe is not well constrained.

Jones et al.\ (2003) originally defined a fossil system to be a galaxy group or cluster with a massive brightest group or cluster galaxy (BGG or BCG) that dominates (by more than two magnitudes in the $r$-band) the rest of the galaxies within $0.5R_{200}$, defined as half of the virial radius of the system, and shows a bolometric X-ray luminosity $L_{X, bol}\geq10^{42}h^{-1}_{50}\;$ergs s$^{-1}$.  This definition has done well in identifying  massive fossil systems but has the potential to miss poorer fossil groups along with being less robust for these poor fossil groups (Dariush et al.\ 2010), since the infall of a lone luminous galaxy would destroy the system's fossil status.  To address this, Dariush et al.\ (2010) proposed altering the classic Jones optical criteria: instead of using ($\Delta m_{12}\geq2.0$), where $\Delta m_{12}$ is the $r$-band magnitude gap between the first and second rank galaxies, Dariush et al.\ (2010) requires ($\Delta m_{14}>2.5$), where $\Delta m_{14}$ is the $r$-band magnitude gap between the first and fourth rank galaxies.  This change would allow for the infall of one or two luminous galaxies without destroying the fossil status of the poorer group.  Jones et al. (2003) found that fossil systems should comprise between 8\% and 20\% of all galaxy groups.  However, this study was only done for nearby groups and clusters of $z<0.25$.  A more recent study by Gozaliasl et al.\ (2014)  shows that the fossil group fraction for massive galaxy groups ($M_{200}\sim10^{13.5}\;$M$_{\odot}$), where $M_{200}$ is the mass within a sphere of density equal to 200 times the critical density of the universe, is  $22\pm6\%$ for $z\leq0.6$ and $13\pm7\%$ for $0.6<z<1.2$ if a $\Delta m_{12}\geq1.7$ is required.

Since most fossil systems found to date lie within $z<0.2$, fossil galaxy groups could be old, undisturbed systems, as the infall of any bright galaxies has the potential to destroy their fossil statuses.  The Millennium Simulation supports this idea, as it shows fossils being formed near $z=0.9$ and subsequently being destroyed due to the infall of a bright $L^*$ galaxy which breaks the $r$-band magnitude gap requirement before $z=0$ (von Benda-Beckmann et al. 2008).  However, as some fossils are destroyed in the simulation, others are created as bright members are consumed by the BGG (Ponman et al.\ 1994) thus establishing the required $r$-band magnitude gap.  This result from the simulation suggests that fossil systems could be more of a `fossil phase' that all groups have a likelihood of passing through as opposed to a unique class of object all their own.  There also exists controversial evidence that fossils have a higher than expected mass concentration parameter (Sun et al.\ 2004; Khosroshahi et al.\ 2004, 2006), although some observations of nearby fossils dispute this (Sun et al.\ 2009).  These pieces of evidence point to fossil systems possibly having different initial conditions than most groups.  To further complicate the matter, there are many supposedly old, evolved nearby fossil systems that do not possess cool cores as would be expected from relaxed systems (Sun et al. 2004; Khosroshahi et al.\ 2004, 2006).  These fossil systems needed an energetic event of some kind (such as an AGN turning on or a group merger) in their recent past to heat up their intragalactic mediums (IGM) or destroy any cool cores which is at odds with these fossils being relaxed.

Despite the certain existence of fossil progenitors, little observational work has gone into locating any.  A progenitor to today's fossil systems would be a system at a higher redshift with ongoing or imminent major merging, in mid-assembly of the eventual fossil system's massive BGG\footnote[1]{It is important to note that it is logically possible for a BGG/BCG to be extremely efficient at turning gas into stars, essentially preventing any other large galaxies from forming in the group/cluster.  Fossil systems formed via this channel would simply be ``born that way". }.  The amount of merging concluded by $z=0$ would be sufficent to push the final BGG $r$-band magnitude two magnitudes brighter than any other remaining galaxy member within $0.5R_{200}$.  The progenitor could also be more centrally mass concentrated than other non-fossil groups at similar redshifts if fossil systems are a unique set of galaxy groups as some suggest.  We expect fossil progenitors to be distinctly different from compact groups, since compact groups do not show velocity dispersions indicative of a deep cluster-like potential well, and the formation of fossil systems is not explained by the merger of galaxies in compact groups (La Barbera et al.\ 2010).

An overarching question in the study of fossil systems is what these groups looked like in the early universe.  Are they all old, evolved systems, the inevitable end for all clusters, or are fossils simply a phase that all groups have a probability of transitioning through?  While the former explanation is possible, the likelihood that the entire fossil population is comprised of isolated, undisturbed groups seems low based on the frequency of mergers in the early universe as seen in simulations.  Studies of the Millennium Simulation also cast doubt on this being the sole explanation.  Many fossil systems in the simulation form between $0.3<z<0.6$, and the fossil system BGGs were always larger than their non-fossil counterparts (Kanagusuku et al.\ 2016), suggesting that fossil systems could begin with different initial conditions than most galaxy systems.  By studying the Cheshire Cat fossil progenitor system, it has been demonstrated that if two groups merge, the final product has the potential to be a fossil group once the BGGs of the respective groups merge (Irwin et al. 2015).  It was estimated that the first and second rank galaxies in the Cheshire Cat gravitational lens will merge in 0.9 Gyr, and once this merger concludes, the group will become a fossil system.  Further, optical and X-ray observations revealed that the system is comprised of two separate galaxy groups undergoing a line of sight merger, opening another avenue for fossil system formation: (fossil) group mergers.  This second possible formation mechanism offers an explanation for observed non-cool core fossil systems (as X-ray cooling time scales are typically longer than galaxy merger time scales.)

It is known that fossil systems house massive BGGs at their centers (Jones et al.\ 2003), implying a higher than average mass concentration when compared to normal groups of similar richness.  This enhancement could make fossil-like systems more efficient stong gravitational lenses.  Since gravitational lensing also preferentially selects merging systems along the line of sight, it follows that targeting systems with strong gravitational arcs nearby could be a more efficient way of locating fossil systems and their progenitors.  Kanagusuku et al.\ (2015) found that in the Millennium Simulation, most of today's fossils became fossils between $0.3<z<0.6$ which happens to be the optimal distance to observe groups act as strong lenses (Trentham 1995).  This motivates us to select our sample from the CAmbridge Sloan Survey of Wide ARcs in the Sky (CASSOWARY)  catalog (Belokurov et al.\ 2009, Stark et al.\ 2013) which searches for strong gravitational arcs in the Sloan Digital Sky Survey (SDSS) and has an average lens redshift of $z\sim0.4$.  Our goals in this study are to identify more examples of fossil progenitors (such as the Cheshire Cat) in the CASSOWARY catalog using SDSS photometry, form a catalog of these progenitors, and contrast their properties against fossil and non-fossil systems.  A control set of near-identical, non-lensing galaxy groups will also be analyzed to see if the presence of a strong gravitational arc near a group biases it toward being a fossil system.  We present results from our analysis of all 58 CASSOWARY members along with average cumulative luminosity functions for each category (fossil, progenitor, and normal systems).

In Section 2, we discuss the selection criteria for our sample from the SDSS archive, group scaling relations involved in determining each group's physical parameters, how fossil status is determined, and how average luminosity functions were generated.  Section 3 presents our catalog of lensing fossil systems and potential fossil progenitor systems and compares the results against a control sample of near identical, non-lensing galaxy groups.  In Section 4, we propose a possible fossil system formation timeline using data from SDSS, incorporating fossil progenitors at varying stages of BCG/BGG formation.  Section 5 summarizes our findings.  We adopt the standard  $\Lambda$CDM cosmology with $H_0=70\;$km s$^{-1}$ Mpc$^{-1}$ and $\Omega_{M}=0.286$ throughout this work.

\section{Sloan Data Analysis}
\subsection{\textbf{Selection Criteria}}

The CASSOWARY catalog (Belokurov et al.\ 2009, Stark et al.\ 2013) identified strong gravitational arcs in the SDSS DR7 archive by searching for blue companions or arcs separated from a luminous red galaxy by $\sim$3".  To date, 58 lensing systems have been identified with many having been confirmed via spectroscopic observations of both the lensing and lensed galaxies with typical lensing galaxies lying between $0.2<z<0.7$ and lensed galaxies beyond $z\sim1.5$.  While it was not required that the lens be a galaxy group, it was found that most were.  Moreover, many of the CASSOWARY groups were not previously identified in group catalogs, as they generally have few members and their higher redshifts create difficulties for automated selection techniques.

\begin{figure}[h]
\vspace{-1.25truein}
\hspace{-0.2truein}
{\includegraphics[scale=0.4625]{CSWA_2_color_regions.pdf}}
\vspace{-1.25truein}
\caption{\footnotesize
{The four color-color group inclusion parellograms created for the Cheshire Cat (CSWA 2) for red-ridge ellipticals;  K-corrections and stellar evolutionary corrections have been included. Red triangles show confirmed member galaxies and black circles show interlopers.  Error bars were calculated quadratically from the reported SDSS DR12 magnitude errors. Confirmed member galaxies outside the regions are very faint and were only confirmed via spectroscopic redshifts.}}
\label{fig:colors}
\end{figure}

Our analysis primarily used SDSS DR12 photometric data, which consists of over 100 million cataloged sources observed in the $ugriz$ bands that span over a quarter of the sky.  Unfortunately, due to the vastness of the data set, only a few percent of the brighter galaxies have available spectra, limiting our ability to know precise distances and therefore accurate group membership.  We therefore use photometric redshifts (photoZ) estimates provided in the SDSS archive which are found using the observed colors of the sources and correlate them to a database of spectroscopic redshifts (specZ) of galaxies of similar colors.  This method allows for an estimate of the source's distance, bearing in mind that improvements could be made once a spectrum is taken.

It is important to note where the uncertanities in SDSS's photoZ measures come from and how they can be minimized.  The error in the photoZ estimate directly correlates with the errors in the source's colors, meaning the fainter sources have less accurate corresponding photoZ.  Normally, at the distances of our targets, solely relying on photoZ to determine group membership is not optimal, so we developed a technique to construct reliable photoZ cuts for the other groups without spectroscopic information.  {\it Gemini} GMOS optical spectroscopic redshifts for 48 galaxies for the Cheshire Cat (also a CASSOWARY object known as CSWA 2) were available (Irwin et al. 2015) along with supplemental redshifts found in the literature; these allowed us to confidently determine group membership.  We therefore chose to use this group to construct our photometric inclusion criteria for the other CASSOWARY catalog systems.  By querying SDSS for the photoZ of each specZ confirmed member of the Cheshire Cat fossil progenitor and generating a histogram in photoZ space, we were able to quantify the spread in the photoZ distribution of known member galaxies.  We fit the distribution to a gaussian and found its standard deviation.  Using the photoZ errors provided by SDSS DR12, we determined that a window of $2\sigma$ centered on the peak of the distribution was sufficient to include all confirmed member galaxies.  If a galaxy's photoZ error bars fall within the $2\sigma$ inclusion window it was counted as a potential member and moved on to the next set of member cuts.

In order to count red-sequence (or red-ridge) elliptical galaxies for use in group scaling relations, we incorporate a series of color cuts to exclude blue, late-type galaxies as well as any red-ridge ellipticals at the wrong redshift.  The average redshift of the CASSOWARY groups, combined with most CASSOWARY member galaxies being red-ridge ellipticals, drove us to omit the $u$-band; the SDSS $u$-band errors were large for even the brightest BGG in our sample suggesting that including the $u$-band would add little information to our analysis.  Thus, every color combination using the $griz$ bands was inspected for the confirmed Cheshire Cat member galaxies.  The tighest groupings would limit the number of interlopers being accidentally counted among the red-ridge elliptical members (Figure~\ref{fig:colors}).  Four parallelograms in color-color space were chosen and are displayed in Table~\ref{tab:cuts}.  Galaxies outside any parallelogram are excluded.  These two criteria were used to determine group membership status for the rest of the CASSOWARY catalog.  The lack of spectra for each group means that when using these cuts the results will not be perfectly pure, however statistically it is complete.

To account for the relatively high redshift of the CASSOWARY groups,  K-corrections for groups of $z<0.5$ were found using the ``K-corrections Calculator" (Chilingarian et al.\ 2010)\footnote[2]{http://kcor.sai.msu.ru/ }.  For groups with $z>0.5$, K-corrections were estimated using $D_{n}4000$ measures (Westra et al. 2010) assuming a typical elliptical galaxy $D_n4000\sim1.65$.  Stellar evolutionary corrections for the $griz$ bands were taken from Roche et al.\ (2009).  Finally, we imposed a corrected $i$-band magnitude limit of $m_i\;<\;20.84$  to find $N_{200}$, defined as the number of red-ridge ellipticals with $L_{gal}>0.4L^*$ inside the group's virial radius, taken to be $R_{200}$ (Wiesner et al.\ 2012).

In order to find each group's $R_{200}$ and other physical properties, we adopted an iterative process using group scaling relations involving $N_{200}$ in Lopes et al. (2009),
\begin{equation}
\textnormal{ln}(R_{200})=0.05+0.39 \textnormal{ln}(\frac{N_{200}}{25})h_{70}^{-1} \;\textnormal{Mpc}
\end{equation}
\begin{equation}
	\textnormal{ln}(M_{200})=0.21+0.83 \textnormal{ln}(\frac{N_{200}}{25})h_{70}^{-1} \; \textnormal{M}_{\odot}
\end{equation}
We first adopted a characteristic $R_{200}$ of $1\;$Mpc as a starting point.
This size was used in an SDSS query to extract all galaxies within the circular angular region.  These galaxies were analyzed using our color/photoZ selection criteria which gave us a preliminary value for that system's $N_{200}$.  To deal with any remaining interlopers present within the extraction region, which could slightly inflate $N_{200}$ and consequently $R_{200}$, we took eight regions of angular size equal to the group's current $R_{200}$ immediately adjacent to the target group and applied the same selection process to the galaxies within these regions.  We averaged these eight results together to obtain an expected $N_{200}^{interloper}$ for each group based on the group's location in the sky.  This value was subtracted from the group's $N_{200}$ to arrive at a more accurate richness.  The new $N_{200}$ was used to calculate a new $R_{200}$, and the process was repeated until a single value for $N_{200}$ and $R_{200}$ was converged upon.  By using this method, we were able to agree with other SDSS galaxy cluster richness catalogs (Wen et al. 2012) on $N_{200}$ to within 10\% for the 17 groups with existing $N_{200}$ estimates.  Errors in galaxy colors were dealt with by adding in quadrature, and if a galaxy's color error bars pushed it into the inclusion regions, it was counted as a potential group member; this only significantly change results in the most distant groups.
\smallskip

\begin{table*}
{\scriptsize
\begin{tabularx}{\textwidth}{l c c c c}
\hline
\hline
 & Cut 1 & Cut 2 & Cut 3 & Cut 4 \\
\hline
R1&$(g-z)>-3.63(g-r)+3.0$&$(g-z)<1.05(g-r)+1.11$&$(g-z)<-3.63(g-r)+7.45$&$(g-z)>1.05(g-r)+0.35$ \\
R2&$(r-i)>-0.25(g-i)+0.28$&$(r-i)<0.81(g-i)-0.29$&$(r-i)<-0.25(g-i)+1.23$&$(r-i)>0.81(g-i)-1.04$ \\
R3&$(i-z)>-0.94(r-i)+0.20$&$(i-z)<1.95(r-i)+0.25$&$(i-z)<-0.94(r-i)+1.05$&$(i-z)>1.95(r-i)-1.73$ \\
R4&$(r-z)>-0.41(g-i)+0.63$&$(r-z)<0.90(g-i)+0.33$&$(r-z)<-0.41(g-i)+1.73$&$(r-z)>0.90(g-i)-1.18$ \\
\end{tabularx}
}
\caption{\footnotesize{Equations for lines used to build parallelograms in color space to exclude all galaxies except red-ridge ellipticals useful for group scaling relations with each inclusion region denoted by R1, R2, R3, and R4.  K-corrections and stellar evolutionary corrections are included to bring galaxies to a $z=0$ frame.}}
\label{tab:cuts}
\end{table*}

\subsection{\textbf{Determining Fossil Status}}

With color cuts, photometric redshift cuts, galaxy interloper averages, and reliable $R_{200}$ estimates in hand, we constructed galaxy membership catalogs for all 58 CASSOWARY members.  At this point, we checked the fossil status of each group using both the Jones et al. (2003) criteria of $\Delta m_{12}\geq2.0$ and the Dariush et al. (2010) criteria of $\Delta m_{14}\geq2.5$; systems which satisfied the optical criteria as-is were labeled as fossil systems.  Groups which were not classified fossil systems moved on to the next stage of analysis to determine whether they are fossil progenitors.

Taking the brightest galaxy as the center of the system, we calculated the projected separation of all members from the BGG along with their masses (assuming a mass-to-light ratio of six in the $r$-band).  We took this information and calculated an expected time scale for the galaxy to merge with the BGG for each member using $T_{merge} \approx1.6r_{25}M_{*}^{-0.3}\;$Gyr, where $r_{25}$ is the maximum projected separation of the galaxies in units of $35.7h_{0.7}^{-1}\;$kpc and $M_{*}$ is the sum of the galaxies' masses in units of $4.3\times10^{11}h_{0.7}^{-1}M_{\odot}$ (Kitzbichler \& White 2008).   We adopt this time scale since it likely overpredicts the merger time by a factor of two relative to other methods (Kitzbichler \& White 2008) ensuring the merger is most likely completed in the specified time scale.   Based on the group's look back time, we determined which member galaxies had sufficient time to be cannibalized by the BGG and their luminosities added to the BGG.  Using this `new' BGG luminosity, the fossil status of each group was again checked; if a group became a fossil via this process it was labeled as a fossil progenitor.  Finally, all groups that still did not meet either fossil criteria were labeled for this work as normal groups, as no amount of possible merging could build a large enough BGG for these groups to transition into fossil groups by $z=0$ (Table 2).
\smallskip

\subsection{\textbf{Luminosity Functions}}

Once all 58 CASSOWARY catalog members were sorted based on their fossil status, we converted galaxy apparent $r$-band magnitudes ($m_r$) into absolute magnitudes ($M_r$) using
\begin{equation}
	M_r=m_r-25-5{\textnormal{log}}(\frac{D_L}{1\;{\textnormal{Mpc}}})-K_r-0.85z
\end{equation}
where $D_L$ is the luminosity distance to the group, $K_r$ is the $r$-band K-correction, and the last term is the stellar evolutionary correction involving the group's redshift (Roche et al. 2009).  We then generated three average luminosity functions (one for fossils, one for progenitors, and one for normal groups) using all the member galaxies for each category.
Due to the low galaxy count in the brightest bins of the average luminosity functions, errors in galaxy count were handled using Poisson statistics with $\sigma\approx1+(n+0.75)^{\frac{1}{2}}$ where $n$ is the number of member galaxies within the luminosity bin (Gehrels 1986).
\smallskip

\section{Discussion}

\smallskip
\subsection{\textbf{CASSOWARY Strong Lensing Sample}}

Of the 58 CASSOWARY members, it was found that six are most likely large, lone ellipticals (possessing an $N_{200}<5$) that happen to act as strong gravitational lenses and were not included in any luminosity functions or fossil/progenitor fractions.  Of the remaining 52 strong lensing systems, we found that $13.5\pm2.8$\%\footnote[3]{Errors are reported at 1$\sigma$ confidence.}  are Jones fossils ($\Delta m_{12}\geq2.0$) and $17.3\pm2.6$\% are Dariush fossils ($\Delta m_{14}\geq2.5$), consistent with the expected 8\% to 20\% fossil system rate for randomly selected groups within $z<0.2$ (Jones et al. 2003).  We found that $23.1\pm2.5$\% of the CASSOWARY systems are Jones fossil progenitors and $28.9\pm2.5$\% are Dariush fossil progenitors.  This higher rate of fossil progenitors in the CASSOWARY sample is not surprising considering that Kanagusuku et al. (2016) found that in the Millennium Simulation, systems which were fossils at $z=0$ finished forming their BGG (creating the required $\Delta m_{12}$/$\Delta m_{14}$  $r$-band magnitude gap) between $0.3<z<0.6$ on average.  Since the average redshift of the CASSOWARY members is $z\sim0.4$, we expect to see a collection of near-fossil systems, as we are seeing analogs to today's fossil systems in mid-cannibalization of their $L^*$ members.  An interesting thing to note is that if one assumes {\it all} the CASSOWARY progenitors and fossils become/stay fossils, one arrives at a $z=0$ Jones fossil percentage of $36.6\pm2.4$\% and a Dariush percentage of $46.1\pm2.3$\%; this implies we should see far more nearby fossil systems than we currently do.  One explanation why we do not see such an overabundance of nearby fossils is that fossils are transitory in nature, and the look back time is long enough to allow some bright galaxies to fall within half of the fossil's virial radius thereby breaking the fossil's status.  Another explanation is that in the CASSOWARY strong lensing sample, we are seeing a subset of systems that are more likely to be fossil systems; these two hypotheses will be explored further in the non-lensing control sample section.

\begin{center}

\begin{longtable*}{llllp{1.0cm}|c|l|c|p{0.8cm}}
\caption{Fossil Status of All CASSOWARY Members} \\
\hline
Name	& RA & Dec &$N_{200}$&$M_{200}\times10^{14}$M$_{\odot}$&P$_J$/P$_{D}$/F$_{J}$/F$_{D}$& $\Delta m_{12}$/$\Delta m_{14}$ & $\Delta m_{12}^{merge}$/$\Delta m_{14}^{merge}$&$t_{merge}$ (Gyr)\\
\hline
\hline
CSWA 1	& 177.1381\degree& 19.5008\degree&5&0.32& x/x/\checkmark/x & {\bf2.3}/2.4 &&\\
CSWA 2	&159.6816\degree&48.8216\degree&10&0.58& \checkmark/\checkmark/x/x & 0.2/2.4 & {\bf3.3}/{\bf3.7}&0.9\\
CSWA 3	&190.1345\degree&45.1508\degree&12&0.67& x/x/x/x & 0.2/1.4 &  \\
CSWA 4	&135.3432\degree&18.2423\degree&32&1.52& x/x/x/\checkmark & 1.4/{\bf3.1} & &\\
CSWA 5	&191.2126\degree&1.1122\degree&14&0.78& x/\checkmark/x/x & 0.2/1.6 & 1.9/{\bf2.9}&4.1\\
CSWA 6 	&181.5087\degree&51.7082\degree&18&0.94& x/x/\checkmark/\checkmark& {\bf2.2}/{\bf2.8} &&   \\
CSWA 7	&174.4169\degree&49.6099\degree&20&1.05& x/x/x/x & 1.0/1.2 &&   \\
CSWA 8	&182.3487\degree&26.6796\degree&61&2.60& x/x/x/x & 1.6/2.2 && \\
CSWA 9	&186.8281\degree&17.4311\degree&23&1.16& x/x/x/x & 0.7/1.9 &&  \\
CSWA 10	&339.6305\degree&13.3322\degree&30&1.43& x/\checkmark/x/x & 1.4/2.2 & 1.9/{\bf2.8}&3.9\\
CSWA 11	&120.0544\degree&8.2023\degree&26&1.28& x/x/\checkmark/x & {\bf2.0}/2.1&& \\
CSWA 12	&173.3049\degree&50.1445\degree&45&2.00& x/x/x/x & 0.7/1.4 &   \\
CSWA 13	&189.4008\degree&55.5619\degree&26&1.27& x/x/x/x & 0.4/1.8 &   \\
CSWA 14	&260.9007\degree&34.1995\degree&18&0.93&\checkmark/\checkmark/x/x & 0.6/1.9 & {\bf2.5}/{\bf2.8}&3.6\\
CSWA 15	&152.2491\degree&19.6215\degree&51&2.24& x/x/x/x & 0.3/1.6 &&  \\
CSWA 16	&167.7653\degree&53.1486\degree&79&3.20& x/x/x/x & 0.1/1.1 &&  \\
CSWA 17	&174.5373\degree&27.9085\degree&74&3.04& x/x/x/x & 1.3/1.4 &&  \\
CSWA 18	&173.5281\degree&25.5598\degree&25&1.25& x/x/x/x & 1.1/2.0 &&  \\
CSWA 19	&135.0110\degree&22.5680\degree&17&0.89& x/x/x/x & 0.6/0.8 &&  \\
CSWA 20	&220.4548\degree&14.6890\degree&1&0.09& x/x/x/x & && \\
CSWA 21	&5.6705\degree&14.5196\degree&28&1.36& x/x/x/x & 1.0/1.3 &&   \\
CSWA 22	&26.7334\degree&-9.4979\degree&60&2.56& x/x/x/x & 1.1/1.8 &&   \\
CSWA 23	&126.8701\degree&22.5483\degree&88&3.49& x/x/x/x & 1.1/1.5 &&  \\
CSWA 24	&227.8281\degree&47.2279\degree&13&0.71&\checkmark/\checkmark/x/x & 2.0*/2.6* & {\bf2.1}/{\bf2.7}&2.2\\
CSWA 25	&162.4298\degree&44.3432\degree&50&2.21& x/x/x/x & 1.3/1.6 &&  \\
CSWA 26	&168.2944\degree&23.9443\degree&83&3.35&\checkmark/x/x/\checkmark & 1.5/{\bf2.5} & {\bf2.9}/{\bf3.4}&2.1\\
CSWA 27	&247.4773\degree&35.4776\degree&63&2.66& x/x/x/x & 0.2/0.7 &&  \\
CSWA 28	&205.8869\degree&41.9176\degree&31&1.48&\checkmark / \checkmark/x/x &1.8/2.3& {\bf2.1}/{\bf2.8}&1.6\\
CSWA 29	&130.0877\degree&10.8702\degree&2&0.18&  x/x/x/x & &&  \\
CSWA 30	&132.8604\degree&35.9705\degree&31&1.46& x/\checkmark/x/x & 0.7/2.2 &1.1/{\bf2.6}&2.0\\
CSWA 31	&140.3573\degree&18.1715\degree&5&0.32& x/x/\checkmark/x & {\bf2.1**}/{\bf**} &&  \\
CSWA 32	&153.7740\degree&55.5051\degree&1&0.09& x/x/x/x &  &&   \\
CSWA 33	&162.3475\degree&35.7447\degree&26&1.26& x/x/x/x &1.0/2.4 &&  \\
CSWA 34	&4.2564\degree&-10.1531\degree&3&0.19& x/x/x/x & &&  \\
CSWA 35	&149.4133\degree&5.1589\degree&9&0.52& x/x/x/x &1.1/1.3 &&  \\
CSWA 36	&181.8996\degree&52.9165\degree&26&1.25& x/\checkmark/x/x &0.4/1.3 &0.6/{\bf2.6}&3.0\\
CSWA 37	&199.5480\degree&39.7075\degree&19&0.98& \checkmark/\checkmark/x/x &0.8/1.7 &{\bf3.0}/{\bf3.3}&3.5\\
CSWA 38	&186.7154\degree&21.8737\degree&167&5.97& x/x/x/x &1.1/1.4 &&  \\
CSWA 39	&231.9376\degree&6.8761\degree&98&3.83& x/x/x/x &0.3/0.4 &&  \\
CSWA 40	&148.1676\degree&34.5795\degree&77&3.13& x/x/x/x & 0.1/0.3 && \\
CSWA 41	&222.6277\degree&39.1386\degree&20&1.03& x/x/\checkmark/\checkmark& {\bf2.2}/{\bf2.7} &&   \\
CSWA 102 &14.7039\degree&-7.3660\degree&5&0.30& x/x/\checkmark/\checkmark&{\bf2.5}/- && \\
CSWA 103 &26.2679\degree&-4.9311\degree&4&0.27& x/x/x/x  & &&  \\
CSWA 104 &167.5738\degree&64.9965\degree&2&0.12& x/x/x/x & && \\
CSWA 105 &168.7683\degree&16.7606\degree&7&0.41&\checkmark/\checkmark/x/x &0.9/1.8&{\bf2.8}/-&4.7\\
CSWA 107 &176.8471\degree&33.5314\degree&41&1.85& x/x/x/x &0.9/1.5 &&  \\
CSWA 108 &179.0228\degree&19.1868\degree&13&0.70& x/x/x/x &1.1/1.4 &&  \\
CSWA 111 &345.0719\degree&22.2249\degree&25&1.23& x/x/x/x &1.1/1.6 &&  \\
CSWA 116 &25.9589\degree&16.1274\degree&7&0.85&\checkmark/x/x/x &1.9/2.2 &{\bf2.9}/{\bf3.4}&2.0\\
CSWA 117 &150.5106\degree&60.3404\degree&30&1.44& x/\checkmark/x/x &0.8/1.6 &1.9/{\bf2.7}&3.4\\
CSWA 128 &299.6473\degree&59.8495\degree&42&1.90& x/\checkmark/\checkmark/x &{\bf2.0}/2.3 &{\bf2.1}/{\bf2.7}&1.9\\
CSWA 139	&121.8814\degree&44.1800\degree&8&0.45&\checkmark/\checkmark/x/x &0.8/2.0 &{\bf2.6}/-&3.4\\
CSWA 141	&131.6977\degree&4.7680\degree&22&1.11& \checkmark/x/x/x &1.4/1.8 &{\bf2.0}/2.4&0.6\\
CSWA 142	&133.6197\degree&10.1373\degree&54&2.34&\checkmark/\checkmark/x/x &0.6/1.2 &{\bf2.}1/{\bf3.1}&2.9\\
CSWA 159	&335.5358\degree&27.7596\degree&6&0.36&\checkmark/x/x/\checkmark&1.3/{\bf2.8} &-/-&3.6\\
CSWA 163	&329.6820\degree&2.9584\degree&18&0.93& x/\checkmark/x/x &1.4/2.1 &1.8/{\bf2.7}&3.2\\
CSWA 164	&38.2078\degree&-3.3906\degree&10&0.55& x/x/x/\checkmark &1.7/{\bf3.2} &&  \\
CSWA 165	&16.3318\degree&1.7489\degree&14&0.78& x/x/x/\checkmark &1.7/{\bf2.8}  &&  \\
\hline
\caption{\footnotesize{All 58 CASSOWARY members and their general properties.  Since both Jones et al.\ (2003) and Dariush et al.\ (2010) criteria were used to determine fossil status, we include the magnitude gap in the $r$-band between the first and second rank galaxy ($\Delta m_{12}$) along with the first and fourth rank galaxy ($\Delta m_{14}$)  within $0.5R_{200}$.  The P$_{J/D}$ and F$_{J/D}$ columns were added to differentiate between Jones(J)/Dariush(D) progenitors or fossils, respectively.  Bolded entries indicate optical fossil status being reached either now or after merging is completed.  Dashes under the merged column indicate all galaxies within $0.5R_{200}$ merging into one BGG by $z=0$.  $t_{merge}$ indicates the expected merger time scale from Kitzbichler \& White's (2008) relation until fossil status in acheved.  *\;Double nucleus detected in archival {\it HST} imaging negating fossil status until merging is finished.  **$\;${\it Gemini} GMOS data from Grillo et al.\ (2013).}}
\label{fig:data}
\end{longtable*}
\end{center}

To better determine if we are truly seeing the progenitors of today's fossil systems, we contrasted the galaxy luminosity functions of each category against one another, as fossil system luminosity functions show a clear deficit in $L^*$ members when compared to comparable sized normal groups and clusters (Gozaliasl et al.\ 2014).  Fossil progenitors might be expected to be a transitional step between the two extremes, losing $L^*$ galaxies as they are consumed by the BCG's.  We created three average cumulative galaxy luminosity functions (a fossil, progenitor, and normal function) from the CASSOWARY lensing sample to ensure we were comparing strong lensing systems to other strong lensing systems{\footnote[4]{In the cases of split identification (e.g., Jones fossil and Dariush progenitor), the Jones criterion was chosen for forming the luminosity functions since it is the most widely cited.}}.  Due to the large amount of overlap between the two prevailing fossil criteria (Jones/Dariush) in this sample, the cumulative luminosity functions combined both Jones and Dariush fossils/progenitors to minimize errors. It is important to note that the poor fitting at the bright end is due to the ``BCG bump," a known artifact of galaxy mergers in the centers of clusters and groups (Hansen et al.\ 2005).  Excluding the BGGs, we found that overall, the lensing fossil and normal population fits are nearly identical.  However, when the BGGs are introduced into the data set, we found that the fossil luminosity function greatly diverged from the normal luminosity function at the bright end, as expected (Figure~\ref{fig:lensing}). The progenitor luminosity function matched the normal function at the faint end.  However, the progenitor function fell between the normal and fossil functions at the bright end suggesting that fossil progenitors are currently losing their intermediate members while gaining very bright members, thereby moving the groups closer to fossil status.  This supports the notion that fossil systems form their massive BCGs via cannibalization of intermediate mass member galaxies.  SDSS images of the inner regions of groups classified as fossil progenitors also very often show an extremely crowded environment near the BCG further supporting this mechanism of fossil formation (Figure 6).

\begin{figure*}
\vspace{-1.2truein}
\hspace{-0.20truein}
\begin{tabular}{cc}
	{\includegraphics[scale=0.44]{CSWA_luminosity_functions.pdf}}&\hspace{-0.45 truein}{\includegraphics[scale=0.44]{CSWA_in_out.pdf}} \\
\end{tabular}
\vspace{-1.0truein}
\caption{\footnotesize
{{\it Left:} Cumulative luminosity functions of CASSOWARY systems separated by classification: normal, progenitor, and fossil.  The curves represent the best fit Schechter function found using least squares fitting for each population excluding the BGGs.  The progenitor function falls between the normal and fossil functions at the bright end ($L\gtrsim 10^{11}$ L$_{\odot}$), supporting the idea that fossil progenitors are a bridge between normal and fossil systems.  {\it Right:} The same luminosity functions binned by inner and outer half virial radii.  The inner regions exaggerate differences between the three populations which may be due to frequent mergers/interactions of members.  The outer regions show little statistical difference from one another, however as a whole, the outer members possess more bright galaxies than the best-fit model would suggest.}  Error bars are at $1\sigma$.}
\label{fig:lensing}
\end{figure*}

Much work has been done investigating the global deficit in intermediate-luminosity members and the value/evolution of the faint-end slope of the fossil luminosity function, finding that the deficit in $L^*$ members is likely due to cannibalization by the BCG and the faint end slope is consistent with normal groups (Lieder et al.\ 2013, Gozaliasl et al.\ 2014, Zarattini et al.\ 2015).  However, little work has gone toward investigating fossil populations in different radial bins, where initial group conditions could still be encoded (particularly in the outer regions).   Binning the average luminosity functions into inner ($r\leq0.5R_{200}$) and outer ($0.5R_{200}<r\leq R_{200}$) regions reveals that this deficit in intermediate mass galaxies/abundance of extremely bright galaxies in fossil progenitors and fossil systems is exaggerated for $r\leq0.5R_{200}$ (Figure~\ref{fig:lensing}; right).  This is likely due to the increased galaxy density near the center effectively speeding up the galaxy interaction rate,
therefore the central regions of fossil or near-fossil systems should show the most extreme differences from non-fossils.  While the inner progenitor fit is nearly identical to the normal fit, when the BGGs are included in the histogram a clear difference can be seen, placing it firmly between the normal and fossil functions.

We quantified the statistical significance of these differences between datasets via a one-sided K-S test which gives the probability that differing data sets come from the same cumulative distribution function.  For $0<r\leq 0.5R_{200}$, even with the relative lack of data points, lensing fossil systems showed only a 0.84\% chance of being identical to normal lensing systems.  Due to a larger sample size, lensing progenitors also proved to be significantly different than normal lensing systems with only a 0.01\% chance of being identical within $r\leq 0.5R_{200}$.  Unfortunately, due to insufficient galaxy counts in the lensing fossil population, we were not able to find a significant difference between lensing fossils and progenitors.  The outer half virial radii galaxies showed no statistically significant differences between the populations, suggesting that most differences in galaxy populations for fossil systems exist near the center where the most processing has occurred.  However, all lensing systems exhibit an average $\sim2\sigma$ deviation from a Schechter function in the outer $0.5R_{200}$ for galaxies brighter than $L^*$ which is surprising.  Since these galaxies are farther away from the center, one would expect them to be much less processed and therefore be better represented by the fits.  Additionally, there are no BGGs in the outer regions to skew the data away from a Schechter function.  While by no means definitive, this suggests that lensing systems may form differently from non-lensing systems of comparable mass.

\smallskip
\subsection{\textbf{Non-lensing Control Sample}}

\begin{figure*}
\vspace{-1.0truein}
\hspace{-0.20truein}
\begin{tabular}{cc}
	{\includegraphics[scale=0.44]{control_lum_functions.pdf}}&\hspace{-0.45 truein}{\includegraphics[scale=0.44]{control_in_out.pdf}} \\
\end{tabular}
\vspace{-1.0truein}
\caption{\footnotesize
{{\it Left}: Cumulative luminosity functions of our non-lensing control sample for all member galaxies.  For non-lensing groups, differences in the bright end (while still visually apparent) are not as prominent as our lensing systems, supporting the existence of a strong lensing bias.  {\it Right}: The same non-lensing luminosity functions binned by inner and outer half virial radii.  The inner regions for non-lensing fossils and progenitors also show the excess of bright member galaxies when compared to normal groups, though again, less pronounced than in our lensing sample.  The outer regions show virtually no difference between progenitors and normal groups, with fossils only housing a few more bright galaxies.  Error bars are at $1\sigma$.}}
\label{fig:control}
\end{figure*}

All members of our sample exhibit strong gravitational arcs near the BGGs, implying a high central mass concentration for these systems.  To see how/if the presence of strong gravitational arcs biases our fossil system and progenitor findings, we assembled a one-to-one random control sample of non-lensing groups from the Augmented maxBCG cluster catalog (Rykoff et al.\ 2012), Clusters of galaxies in SDSS-III (Wen et al. 2012), and Richness of Galaxy Clusters (Oguri 2014) catalogs.  Control groups were selected to match (within 10\%) each CASSOWARY group in both redshift and galaxy richness, and when multiple matches for control groups were found among the catalogs, the closest match was chosen.  To increase the accuracy of this one-to-one comparison, we found two non-lensing matches for each CASSOWARY group\footnote[5]{Only one non-lensing match was able to be found for CSWA 31 due to its high redshift ($z=0.683$) and poor member count ($N_{200}=5$).}.  Using the same photoZ and color cuts as the lensing sample, we found fossil percentages of $2.9\pm1.6$\% (Jones) and $13.6\pm1.2$\% (Dariush) and fossil progenitor percentages of $17.5\pm1.2$\% (Jones) and $25.2\pm1.1$\% (Dariush) showing that while being a strong gravitational lens does not significantly alter the Dariush fossil fraction, it does greatly improve the chances that a fossil will be a classic Jones fossil.  The progenitor fraction is consistent between the control sample and lensing sample suggesting that in general the progenitor fraction is not greatly affected by the presence of gravitational arcs.  It is important to note that while the non-lensing Dariush fossil fraction is consistent with the prediction from Gozaliasl et al.\ (2014) for $z\leq0.6$, our non-lensing Jones fossil fraction of $2.9\pm1.6$\% is far below their estimate of $22\pm6\%$.  This could be partially due to the Gozaliasl et al.\ (2014) $\Delta m_{12}\geq1.7$ being less restrictive than the classical Jones et al.\ (2003) $\Delta m_{12}\geq2.0$.  Also, while their sample includes all groups within $z\leq0.6$, we only have two groups within $z\leq0.2$ making our sample somewhat different from theirs.

To see if the lower fossil occurrence rate in our sample could be due to inadvertantly including too many bright galaxies (since spectroscopic data is lacking for these groups) we calculated the bright galaxy ($L_{\odot}>0.4L^*$) overdensity within each CASSOWARY group compared to the surrounding regions.  In non-fossil CASSOWARY groups, an overdensity of bright galaxies (sufficient to prevent that group from being a fossil) within $0.5R_{200}$ was confirmed above $13\sigma$ confidence, indicating our fossil fractions are reliable.  The apparent fossil deficit could be accounted for due to the redshift range of our sample ($0.2<z<0.7$).  Since most fossils discovered lie near $z\sim0.1$, and Kanagusuku et al.\ (2016) found in the Millennium Simulation that most $z=0$ fossils made the transition between $0.3<z<0.6$, there could be a lower fossil fraction in our samples.

Average galaxy luminosity functions were also generated for the non-lensing sample to compare against the lensing sample to see how strong lensing might affect a group's galaxy population (Figure~\ref{fig:control}). The non-lensing fossil and normal luminosity functions exhibit similar behavior as the lensing sample.  Non-lensing fossils show a 0.01\% chance of being identical to non-lensing normal {\it and} progenitor systems, confirming that on average fossil systems' inner regions house a different population of galaxies than non-fossils.  Interestingly, non-lensing progenitors showed virtually no differences from non-lensing normal systems at any radii; this reinforces our hypothesis that the presence of a strong gravitational arc marks the most extreme examples of fossil formation at all stages.  Also, the non-lensing progenitor fit falls between the other two in each radial bin.  Reintroducing the BGGs maintains this in-between state for the non-lensing progenitors.  Applying the K-S test to the non-lensing populations showed again that significant differences only appear within $0<r\leq 0.5R_{200}$.

\subsection{\textbf{Comparing Lensing vs. Non-Lensing Samples}}

\begin{figure*}
\vspace{-1.0truein}
\hspace{-0.60truein}
\begin{tabular}{cc}
	{\includegraphics[scale=0.49]{lensing_control_compare.pdf}}&\hspace{-0.6 truein}{\includegraphics[scale=0.49]{data_comparison.pdf}} \\
\end{tabular}
\vspace{-1.0truein}
\caption{\footnotesize
{{\it Left:} Contrasting the best fit Schechter functions, excluding the BGGs, of lensing/non-lensing samples at inner-half and full virial radii bins separated by fossil status.  The CSWA strong lensing fit is shown with a solid line and the non-lensing control fit is represented by the dashed line with the derived $L^*$ for each model marked.  {\it Right:}  Galaxy member data taken from SDSS with BGGs included.  Circles represent the lensing systems, and triangles mark the non-lensing systems.  In each case, lensing groups show a deficit of intermediate-luminosity members and an excess in bright members implying that lensing systems could be an example of the most extreme fossil systems along with them being more likely to become fossils as opposed to similarly sized non-lensing systems. }}
\label{fig:compare}
\end{figure*}

A comparison of the lensed vs.\ non-lensed fossil luminosity function revealed interesting distinctions; on average, lensing fossil systems lack intermediate mass galaxies and house larger BGGs than their non-lensing counterparts (Figure~\ref{fig:compare}).  While the latter is not terribly surprising (as larger galaxies are more likely to act as good gravitational lenses), the former offers help explaining why the lensing sample has significantly more Jones fossils than the non-lensing sample; the lensing sample is the extreme case in fossil formation.  In CASSOWARY fossils (and progenitors to a lesser extent), we are seeing elevated $L^*$ cannibalization resulting in intermediate galaxy deficits and an overrepresentation of extremely bright galaxies.  Comparing the non-lensing progenitor fit to its lensing counterpart reveals even sharper differences between non-lensing fossils and lensing fossils (top of Figure~\ref{fig:compare}). On average, the lensing progenitors have fewer bright $L^*$ galaxies than the non-lensing progenitors.  This could be an indication of a strong lensing selection bias. Since the presence of strong lensing indicates a high mass concentration, lensing progenitors could have already had most of their $L^*$ members consumed by the BGGs.

Since Dariush fossils do not need such a large luminosity difference between the BGG and the next ranked galaxies, the non-lensing control sample holds many more Dariush fossils than Jones fossils.
The non-lensing control sample also indicates that while the presence of a strong gravitational arc does not strongly affect the likelihood of finding the progenitors to today's fossils around $z\sim0.4$, it does appear to greatly increase the probability of locating Jones fossil systems.  Applying the K-S test, this time to lensing vs.\ non-lensing systems, again shows significant differences only within $0.5R_{200}$.  For the inner regions, normal systems proved to be consistent between lensing and non-lensing systems.  Progenitors, on the other hand, showed a 0.01\% chance of being identical; such a striking result means that a strong lensing bias may very well exist.  Since lensing progenitors on average show different galaxy populations than non-lensing progenitors, it can be inferred that the same (if not more) can be said for lensing vs.\ non-lensing fossils.  Unfortunately, errors in bin count for both fossil populations kept us from arriving at any meaningful result; this can be remedied by increasing the sample size.

When one compares the best fit Schechter functions of lensing vs.\ non-lensing systems as a whole, an interesting distinction is found: at every radial bin, systems acting as strong gravitational lenses exhibit an intermediate-luminosity member deficit regardless of fossil status.  While this is expected near the center of most groups, (the act of forming the BGG consumes many of these galaxies thereby shifting $L^*$ toward the faint end) this deficit supports the existence of a strong lensing bias toward fossil-like systems.  To test whether or not these lensing systems are preferentially selecting systems with different initial conditions from normal groups, thereby supporting the idea that some fossils are born differently than most systems, the outer regions must be probed to see if the galaxy populations differ there as well.  However the overall lack of members in the outer half virial radii of the systems made any conclusions statistically insignificant.

\section{\textbf{Progenitor - Fossil Properties, Timelines, and the Longevity of Fossil Systems}}

Since fossil systems are traditionally believed to be relaxed systems, it stands to reason that fossil IGMs should possess well-developed cool cores.  While some do{\footnote[6]{Although some fossil systems show cool cores, many are smaller than one would expect from the cooling time.}}, there are many which instead show flat radial temperature profiles or even temperature spikes near their centers (Sun et al.\ 2004; Khosroshahi et al.\ 2004, 2006).  More recently, Irwin et al.\ (2015) discovered that the Cheshire Cat fossil progenitor was likely being shock heated near its center because of a group merger, meaning that even after the BGG merger is complete it will be seen as a non-cool core fossil system for $\sim$2 Gyr.  Miller et al.\ (2012) also noticed odd asymmetries in the IGM of many nearby fossils, suggesting that their histories were complex and possibly violent as well.  Such irregularities in the hot gas could be better understood if similar irregularites are also seen in fossil progenitors. Following this thought, one of the most promising systems in the CASSOWARY catalog is CSWA 26. The $z=0.336$ fossil progenitor is six times more massive than the Cheshire Cat and shows some evidence of being another group merger, the result of which would be the most massive fossil cluster yet discovered which does not possess a cool core.  The photometric redshifts of CSWA 26's members indicate a possible double peak about the BCG and the third rank galaxy (Figure \ref{fig:CSWA_26_hist}).  The mean redshift of each peak shows a radial velocity difference of $\Delta v>1000$ km s$^{-1}$ giving weight to this being a large group merger.

\begin{figure}[h]
\vspace{-1.25truein}
\hspace{-0.25truein}
{\includegraphics[scale=0.445]{CSWA_26_histogram.pdf}}
\vspace{-1.25truein}
\caption{\footnotesize
{Histogram of all candidate members of the fossil progenitor CSWA 26.  The double peak near $z=0.34$ is suggestive of another group merger scenario similar to the Cheshire Cat (CSWA$\;2$) but on a much larger scale.}}
\vspace{-0.1truein}
\label{fig:CSWA_26_hist}
\end{figure}

Since we have many progenitors in the CASSOWARY catalog with a wide range of merging time scales until transitioning into a fossil system (Table 2), we can form a rough timeline of today's average fossil system's formation process from its beginning, through the cannibalization phase building the large BGG, and finally concluding with a fossil system housing a large BGG and possessing a deficit in bright $L^*$ galaxy members. To better illustrate the hypothesis of formation of a fossil system through the progenitor phase, we have assembled a collage of SDSS images from the CASSOWARY catalog (Figure \ref{fig:collage}).  We order them to simulate the building of a fossil BGG via cannibalization of $L^*$ members.  Early in the progenitor phase, we expect there to be many bright galaxies present in the group and concentrated near the BGG, since dynamical friction has slowed the orbits of the largest galaxies and caused them to fall inward over the group's history.  As time until fossil status is achieved shortens, more and more $L^*$ galaxies will merge with the BGG subsequently shifting the galaxy luminosity function of the group toward a fainter population leaving only one or two bright galaxies near the BGG.  Once the last bright member merges with the BGG, a fossil system will form housing an elongated (possibly asymmetric) BGG.  As the BGG begins to relax after the final major merger it will eventually settle into a massive symmetric elliptical galaxy stereotypical of fossil systems.  A follow up study of progenitors is currently being done using {\it Chandra/HST} data to better see how the hot gas evolves alongside the stellar population as a group draws closer to the fossil threshold (Johnson et al.\, in preparation).  We expect to see a correlation between a progenitor's X-ray properties and time until fossil status is achieved as well.

\begin{figure*}
\centering
\begin{tabular}{lccc}
	{a. {\includegraphics[angle=-90, scale=0.1325]{CSWA_23_2.pdf}}}&{b. \includegraphics[angle=-90, scale=0.1325]{CSWA_10_2.pdf}}&{c. \includegraphics[angle=-90, scale=0.1325]{CSWA_26_2.pdf}}&{d. \includegraphics[angle=-90, scale=0.1325]{CSWA_30_2.pdf}}\\
	{e. \includegraphics[angle=-90, scale=0.1325]{CSWA_2_2.pdf}}&{f. \includegraphics[angle=-90, scale=0.1325]{CSWA_11_2.pdf}}&{g. \includegraphics[angle=-90, scale=0.1325]{CSWA_4_2.pdf}}&{h. \includegraphics[angle=-90, scale=0.1325]{CSWA_1_2.pdf}} \\
\end{tabular}
\caption{\footnotesize
{An SDSS collage of CASSOWARY systems at various stages in BGG formation, in time order from a to h.  Appearing in order from upper left to bottom right along with time until fossil formation is complete: normal group (a.) CSWA 23, (b.) CSWA 10 (3.9 Gyr), (c.) CSWA 26 (2.4 Gyr), (d.) CSWA 30 (2.0 Gyr), (e.) CSWA 2 (0.9 Gyr), (f.) CSWA 11 (double nuclei; $<100$ Myr), (g.) CSWA 4 (young fossil with ongoing mergers), (h.) CSWA 1 (relaxed fossil).}}
\label{fig:collage}
\centering
\end{figure*}

The high redshift of the CASSOWARY sample is beneficial for locating possible precursors to nearby fossil systems. However, the look back time that allows large BGGs to form also allows new galaxies to fall within $0.5R_{200}$, potentially destroying a system's fossil status before $z=0$.  To account for this chance as conservatively as possible, we define a `danger zone' for each system; this zone is the maximum projected distance from which a galaxy could free-fall inside $0.5R_{200}$ within the look back time.  The free-fall time is given by:
\begin{equation}
t_{ff}=\frac{\pi}{2}\frac{r^{\frac{3}{2}}}{\sqrt{2G(M+m)}}
\end{equation}
where we took $r$ to be the distance from the galaxy down to $0.5R_{200}$ for the group, as this is the threshold for a member to be considered in a system's fossil classification.
In an accelerating universe, the ultimate mass of a galaxy cluster in the far future is around two times the virial mass ($M_{200}$) at $z=0$ (Busha et al.\ 2005).  Since the density inside the virial radius is constant by definition, the virial radius will scale as $M_{200}^{1/3}$ making the final radius $\sim$1.25 times the current $R_{200}$.  The turnaround radius (the point beyond which matter near an overdensity of a certain mass will not collapse inward but expand) is defined as $2R_{200}$, implying that the final turnaround radius will be $\sim2.5R_{200}$ at $z=0$.  However, since our targets are at $z>0$, their masses have grown between when we have observed then and now.  Assuming a mass growth of a factor of four between then and now, that gives a $z=0$ virial radius of $1.6R_{200}^{obs}$ and a $z=0$ turnaround radius of $3.2R_{200}^{obs}$.  For this work, we adopt a turnaround radius of $3.0R_{200}$ and set this as the upper limit to our `danger zone.'
All bright galaxies within the `danger zone' that passed our group inclusion criteria and were bright enough to violate either the Jones ($\Delta m_{12}\geq2.0$) or Dariush ($\Delta m_{14}\geq2.5$) fossil criteria were flagged.  Groups with flagged galaxies within their `danger zone' may still be fossil progenitors, however one cannot rule out the possibility that one or more of the flagged galaxies will eventually fall into the group potential.  We found that out of the 28 strong lensing fossil/progenitor systems in the CASSOWARY catalog, only two (CSWA 26 and CSWA 159) have no nearby galaxies bright enough to endanger their eventual fossil statuses making these true fossil progenitors\footnote[7]{It is interesting to note that both CSWA 26 and CSWA 159 are classified as Dariush fossils and Jones progenitors, implying that in both cases we are witnessing the formation of extremely massive BGGs compared to their respective group richnesses.}.

In simulations, the entire lifetime of fossil systems can be chronicled by observing when and if bright galaxies fall into the group.  Observationally, it is more difficult, as we do not know the proper motions of all galaxies around the group.  Spectroscopic redshifts of galaxies in and near the group can constrain the radial velocities, however the tangential components remain unknown.  This means we cannot know for certain which bright galaxy outside $0.5R_{200}$ will fall inwards.  Therefore, the result that only $\sim$$7\%$ of $0.2<z<0.7$ fossil/progenitor systems will stay fossils until $z=0$ is an extremely conservative estimate.  Spectra of observed fossils and progenitors at higher redshifts have the potential to increase this estimate, as some bright galaxies will undoubtedly be eliminated, being identified as either foreground or background.

\section{Summary}
The progenitors to today's fossil systems have been shown to exist in the universe and can be located.  Kanagusuku et al.\ (2016) found that in the Millennium simulation most of today's fossils finished forming their BGGs between $0.3<z<0.6$ which also happens to be the optimal distance to see strong gravitational lensing.  The discovery of the Cheshire Cat fossil group progenitor (CSWA 2) in the CASSOWARY catalog of strong lensing events in SDSS prompted us to analyze the remaining 57 CASSOWARY members to see if more progenitors could be located.

In this study of CASSOWARY strong lensing systems, we find fossil rates of $13.5\pm2.8\%$ and $17.3\pm2.6\%$ for the Jones et al.\ (2003) and Dariush et al.\ (2010) criteria respectively which is consistent with the expected rate of 8\% -- 20\% of all groups being fossil systems (Jones et al.\ 2003).  This contrasts our non-lensing control fossil rates of $2.9\pm1.6\%$ and $13.6\pm1.2\%$ indicating the presence of a strong lensing bias toward classical (Jones) fossil formation.  Our CASSOWARY progenitor rates of $23.1\pm2.5\%$ and $28.9\pm2.5\%$ (Jones and Dariush respectively) are also elevated compared to our non-lensing progenitor rates of $17.5\pm1.2\%$ and $25.2\pm1.1\%$.  Average galaxy luminosity functions for each class of system (normal, progenitor, and fossil systems) confirmed that fossil progenitors fall between the normal and fossil fits, indicating the formation of the $L^*$ galaxy deficit observed in fossil systems.  For the CASSOWARY sample, the progenitor luminosity function showed a slight deficit of intermediate member galaxies supporting the hypothesis that fossil BGGs are formed via cannibalization of their $L^*$ neighbors and that we are witnessing this process in some CASSOWARY fossil progenitors.

A control sample of non-lensing groups at similar redshifts and galaxy counts was complied to compare against the CASSOWARY lensing sytems to see how the conditions leading to a strong gravitational arc near the BGG could bias the sample.  It was found that while being a strong gravitational lens did not appear to affect the odds of a group being a fossil progenitor, it does increase the odds that a group will be a Jones fossil, indicating the existence of a strong lensing bias possibly linked to the inital formation of the lensing group.  Comparing the cumulative galaxy luminosity functions of the non-lensing control sample to the CASSOWARY groups showed the non-lensing progenitor function agreeing more with the non-lensing normal groups rather than transitioning to the non-lensing fossil luminosity function.  This could also be an indication of the strong lensing bias preferentially selecting the most extreme examples of fossil formation (i.e., systems with the highest mass concentration, largest intermediate mass galaxy deficit, largest BGGs).  Additionally, we observe lensing systems possessing an average of $2\sigma$ more bright galaxies than the best fit gives for galaxies outside $0.5R_{200}$.  This is not seen in the non-lensing systems, further supporting the existence of a strong lensing bias toward classical (and possibly older) fossil-like systems.

Within our sample of lensing fossil progenitors, we identified fossil progenitors which could better explain the origin of both cool core and non-cool core fossil systems observed today.  Most fossil progenitors in this study seem to be in the process of forming a massive BGG via $L^*$ member cannibalization making these progenitors of the observed cool core fossil population seen today.  We also find two fossil progenitors that show evidence of being a group-group merger (e.g. the Cheshire Cat and CSWA 26).  When the multiple BGGs merge, a fossil system will be created, however the shocked gas will still be cooling a few Gyrs after the galaxy merging is complete (Irwin et al. 2015) making the Cheshire Cat and CSWA 26 progenitors to the observed non-cool core fossil population seen today.  We are engaged in futher work on this topic using {\it Chandra} imaging of the hot gas of eight fossil progenitors in the CASSOWARY catalog at a range of stages in their evolution toward fossil status, from 100 Myr to 5 Gyr until sufficient merging has concluded to establish the required $r$-band magnitude gap, the goal being to observe the evolution of the hot gas component of a fossil progenitor and whether or not a cool core is present as the BGG forms (Johnson et al.\ in preparation).

\acknowledgments
Acknowledgements:  We thank William Keel, Jeremy Bailin, Yuanyuan Su MacLellan, Eddie Johnson, and Spring Johnson for useful conversations and constructive feedback.  This work was supported by {\it Chandra} grant GO6-17106X and {\it HST} grant HST-GO-14362.003-A.

\begin{footnotesize}