text:
\section{Introduction}

Recently \citet{Brownetal2009} presented submillimetre images of three transition discs LkH$\alpha$\,330, SR\,21N, HD\,135344B. Using a 2D radiative transfer model \citep{DullemondDominik2004}, they showed that the observed inner discs are cleared from dust as the optical through millimetre-wave spectral energy distribution (SED) can be fitted with a disc having an inner cavity. The non-axisymmetric intensity distributions (see fig. 1 of \citealt{Brownetal2009}), however, require further assumptions such as the presence of an embedded planet at large distances from the central star ($\ga50$\,AU) that sufficiently perturbs the disc dust distribution (see, e.g., \citet{Wyattetal1999}). A massive enough planet opens a gap, and the disc's material will be piled up just outside of the giant planet's orbit \citep{deValBorroetal2007}. Although the giant planet induced non-axisymmetric density distribution hypothesis is challenging, we propose in this Paper that the observed asymmetries are due to large observable vortices formed in these discs during their early evolutionary stages.

However, the current most accepted scenarios for giant planet formation, the core-accretion \citep{BodenheimerPollack1986,Pollacketal1996} and the gravitational instability \citep{Boss2001,Boss2006a} face problems forming and keeping giant planets at such a large distances. Within the core-accretion scenario, the optimal places for planet formation is the vicinity of the snow line where the condensation of volatiles into ices increases the surface density of solids by a large factor, increasing the isolation mass above the critical core mass. The location of the snow line determined by the opacity of the dust grains and the mass accretion rate is recently found to be $\la 16\,\rm AU$ \citep{Minetal2011}. For significantly larger distances, the growth rate of the planetary core is strongly limited. On the other hand, forming giant planets on wide orbits via gravitational instability does not face many problems (assuming rapid cooling of the gas), except that it is not yet clear how to keep them there. \citet{Dodson-Robinsonetal2009} investigated three possible mechanisms -- core-accretion (with or without migration), scattering from the inner disc, or gravitational instability -- which might result in massive planets at wide ($\ga35$\,AU) orbits discovered by high-contrast direct imaging observations of Fomalhaut \citep{Kalasetal2008mn}, HR\,8799 \citep{Maroisetal2008} or CD-352722 \citep{Wahhajetal2011mn} for instance. Dodson-Robinson et al. found that the only possible scenario to form stable wide-orbit planets is the gravitational instability, but planets formed this way were recently shown (by \citealp{Baruteauetal2011}) to migrate quite rapidly inwards -- puts in question its feasibility. Note that \citep{Cridaetal2009} have shown that runaway outward migration by planets in mean motion resonance carving a common gap in the protoplanetary disc is an other possible mechanism to place giant planets on wide orbits, if the inner one is significantly more massive than the outer one.

In recent years the formation of vortices in protoplanetary discs have been investigated thoroughly, since they may play an important role in the formation of planets in the core-accretion scenario. Observation of such vortices would have a great importance providing further informations and support for planetary formation theories (see, e.g., \citet{KlahrHenning1997,BargeSommeria1995}).

The classical approach of accretion discs \citep{ShakuraSunyaev1973} postulates a simple scaling relation between the gas viscosity and the level of angular momentum transport parametrized by $\alpha$. An accretion disc coupled dynamically to a weak magnetic field is subject to the so-called magnetorotational instability (MRI) resulting in outward angular momentum transport, which can be described by an effective $\alpha$-viscosity \citep{BalbusHawley1991}. The disc becomes MRI-active, if the ionisation state of the gas is sufficient. Investigating the ionisation processes in MRI-active discs, \citet{Gammie1996} proposed the existence of a radially extended region in the disc midplane, where the accretion is strongly reduced. Near the disc inner edge, the midplane gas is collisionally ionised, owing to the high temperatures \citep{PneumanMitchell1965,Umebayashi1983}. Further out, the temperature is not high enough to collisionally ionise the gas, while the gas is yet dense enough to shield the interstellar cosmic-rays or diluted stellar X-rays. The magnetic field of the ionised disc corona might even deflect the cosmic rays. As a result, only a thin surface layer of the disc is ionised at this region, while at even larger distances, the cosmic rays can penetrate the tenuous disc which might ionise the disc midplane again \citep{Glassgoldetal1997,IgeaGlassgold1999}. As a consequence, the turbulent viscosity is low in the poorly ionised region of the disc midplane, resulting in the formation of the so called {\it dead zone}. Within this region the viscous accretion is taking place mostly in a narrow upper layer of the disc, while it is  stalled in the disc's midplane.

\citet{VarniereTagger2006} and \citet{Terquem2008} demonstrated that density and pressure enhancements (`bumps') appear at the boundaries of a \emph{dead zone} because of the sharp jump in the viscosity. These pressure bumps are found to be unstable to the Rossby wave instability (RWI) \citep{Lovelaceetal1999,Lietal2000}, leading to the formation of a large-scale anticyclonic vortices.

\citet{KretkeLin2007} have already shown that pressure-gradient inversion formed at the inner edge of the dead zone or at the snow line causes the accumulation of dust. The existence of large-scale vortices can also accelerate the formation of planetesimals and planetary embryos in the core-accretion paradigm. \citet{BargeSommeria1995} and \citet{KlahrHenning1997} have already shown that particles tend to get trapped in anticyclonic vortices. At the radial location of the vortex the type-I migration is halted due to radial pressure-gradient inversion, thus not only dust is retained here but larger bodies too. Particles of about centimetre to metre in size subsequently drift into the pressure maxima, i.e., to the vortices \citep{Lyraetal2009,Katoetal2010,Dzyurkevichetal2010} and form gravitationally bound clumps of solids, which can coalesce forming embryos between the masses of the Moon and Mars. The swarm of these embryos evolves further by mutual N-body collisions forming massive cores ($\sim 10 M_{\oplus}$) of giant planets in $\leq 5\times10^5\,\rm yr$ \citep{Sandoretal2011}. Thus, large-scale anticyclonic vortices formed in discs are extremely important for the planetary formation helping to overcome the meter-size barrier problem \citep{BlumWurm2008}, and the time-scale problem of oligarchic growth. From this perspective these vortices can be regarded as `planetary factories'.

The micron sized dust is dynamically coupled to the gas as it is indirectly supported against the settling by gas pressure via the gas drag, thus their spatial distribution is expected to be the same.
Assuming that the three transition discs presented in \citet{Brownetal2009} contain significant amount of gas -- supported by H\,$\alpha$ \citep{Fernandezetal1995}, [OI] \citep{vanderPlasetal2008} and fundamental band CO emissions \citep{Pontoppidanetal2008,Salyketal2009} --, these vortices can be imaged in submillimetre wavelengths as overdense regions in the dust. The inner holes in the submillimetre images of the three discs may indicate that dust is either consumed here by a planet forming process or by coagulation to larger sized grains, which are drifted rapidly into the star. Note that coagulation without fragmentation would make the vortex also invisible as the dust opacity is decreasing with grain size. If the dust trapped inside giant vortices \citep{InabaBarge2006} is accountable for the observed features, fragmentation should be happening inside them, in order to replenish the population of small dust grains.

Inspired by these ideas, we explore, whether a large-scale, long-lived (on $10^6\,\rm yr$ time-scale) non-axisymmetric gas-dust perturbation formed at the outer edge of the disc dead zone could explain the non-axisymmetric submillimetre images. To calculate the density perturbation, we run long-term numerical hydrodynamical simulations on a GPU aided computer. In \S\,2, we present our hydrodynamical disc model, and a study on the robustness of the formation of large-scale vortices. Submillimetre dust continuum images calculated by a 3D radiative transfer using our hydrodynamical results are presented in \S\,3. The Paper closes (\S\,4) with the discussion and concluding remarks, along with a short discussion on the reasonability of our findings.

\section{Vortex formation}

We revisit the development of RWI, \citep{Lovelaceetal1999}, by 2D, grid-based, global hydrodynamical disc simulations, assuming the $\alpha$-type viscosity prescription \citep{ShakuraSunyaev1973}. To model the dead zone, the value of $\alpha$, is smoothly reduced such that $\alpha(R)=\alpha\delta_\alpha(R)$ and
\begin{equation}
	\label{eq:deltaalpha}
	\delta_\alpha(R)=1-0.5\left(1-\alpha_\mathrm{mod}\right)\left[-\tanh\left(\frac{R-R_\mathrm{dze}}{\Delta R_\mathrm{dze}}\right)+1\right],
\end{equation}
where $\alpha_\mathrm{mod}$ is the depth of turbulent viscosity reduction, $R_\mathrm{dze}$ and $\Delta R_\mathrm{dze}$ are the radial distance and the width of the transition region of turbulent viscosity at the dead zone outer edge, respectively. Regarding the present study the effect of the inner edge of the dead zone is not relevant, therefore it is not investigated here. Examples of $\alpha(R)$ used in our models are shown in Fig.\,\ref{fig:alpha_damping}.

\begin{figure}
	\centering
	\includegraphics[width=\columnwidth]{fig1.pdf}
	\caption{Example of smooth viscosity reduction according to Eq. (\ref{eq:deltaalpha}) used in our models. $\delta_\alpha(R)$ for models $\alpha=0.01$, $\alpha_\mathrm{mod}=0.01$, $R_\mathrm{dze}=50\,\rm AU$, and $\Delta R_\mathrm{dze}=1,3,5$\,AU are shown with \emph{black solid}, \emph{blue dashed}, and \emph{red dotted} curves, respectively.}
	\label{fig:alpha_damping}
\end{figure}

To quantify $R_\mathrm{dze}$, we adopt the results of \citet{MatsumuraPudritz2005}. According to their calculations, the outer edge of the dead zone lies between 12\,AU  and 36\,AU, depending on their various model parameters. To induce the vortex formation at distances suggested by submillimetre images of \citep{Brownetal2009}, we set the outer edge of the dead zone to slightly larger distance ($R_\mathrm{dze}=50$\,AU). Note that the exact location of the outer edge of the dead zone is very uncertain, because the physics of the ionisation (depending on the ionisation rate of the stellar X-rays, interstellar cosmic-rays, proton radiation from the high temperature disc corona, decay of radionuclides such as $\rm ^{40}K$ or $\rm ^{26}Al$) and the dependence of the magnetorotational instability on this ionisation degree is not yet well understood (see, e.g., \citealp{TurnerDrake2009}).

\subsection{Numerical simulations}

To follow the evolution of the RWI triggered vortices on several $10^6\,\rm yr$ time-scale, we used {\small GFARGO}\footnote{http://fargo.in2p3.fr/spip.php?rubrique18}, the GPU\footnote{Graphical Processor Unit constructed with several hundred procesors (240 for GTX-280, which we used for the simulations).} version of the {\small FARGO} code \citep{Masset2000}. We slightly modified the CUDA\footnote{http://www.nvidia.com/object/what\_is\_cuda\_new.html} kernel of {\small GFARGO} code in order to handle the viscosity reduction according to Eq.\,(\ref{eq:deltaalpha}).

The dimensionless continuity and Navier-Stokes equations
, describing the surface mass density evolution in an accretion disc, are numerically solved in two-dimensional cylindrical coordinate system. We use the locally isothermal approximation, in which the vertically integrated gas pressure is $p(R)=c_\mathrm{s}(R)^2\Sigma(R)$, where $c_\mathrm{s}(R)$ is the local sound speed, and $\Sigma(R)$ is the vertically integrated surface mass density of gas. The initial surface mass density profile follows the power-law distribution, $\Sigma(R)\sim R^{-0.5}$. We use radially constant disc aspect ratio $h=H(R)/R$, where $H(R)$ is the pressure scale height at the radial distance $R$, for which case $c_\mathrm{s}(R)=\Omega H = \sqrt{G M_*/R^3}(h R)\propto R^{-1/2}$.

\begin{figure}
	\centering
	\includegraphics[width=\columnwidth]{fig2.pdf}
	\caption{The distribution of $Qh$ (where $Q$ is the Toomre parameter and $h$ is the disc aspect ratio) calculated for the base model, when a large-scale vortex is fully-developed. Since $Qh$ is well above unity, even inside the vortex eye, seems reasonable to neglect the self-gravity.}
	\label{fig:Q2d}
\end{figure}

The disc self-gravity can be neglected as long as the Toomre parameter \citep{Toomre1964} is high, $Q=h M_*/(\pi R^2\Sigma)\gg1$, which holds in our unperturbed discs. Note, however, that $Q$ is relevant only for small scale, i.e. high $m$ modes, for global modes (e.g., m=1) the parameter that might really measure the relevance of self-gravity is $Qh$ \citep{Papaloizou2002}. Since $Qh$  might be smaller than unity, particularly inside the vortex owing to the density enhancement, we calculate $Qh$ for the base model ($h=0.05$ and $M_\mathrm{disc}=0.01\,M_\odot$, and neglecting self-gravity) by the time when the large-scale vortex have fully developed, see Fig.\,\ref{fig:Q2d}. As one can see, $Qh$ never reaches unity ($Qh|_\mathrm{min}=3.68$ inside the vortex eye), therefore we neglected the disc self-gravity in our models. Note that the disc mass $M_\mathrm{disc}=0.01\,M_\odot$ that we used in our simulations is in the range ($0.006\leq M_\mathrm{disc}\leq 0.026$) as has been derived by SED modelling for the three transition discs by \citet{Andrewsetal2011}. Nevertheless, $Qh$ might be smaller than unity for a disc which is three times more massive than the one we used, thus for that case one should take into account the self-gravity.

The disc inner and outer boundaries are fixed at 4\,AU\ and 100\,AU, respectively. The computational domain contains 256 radial logarithmically and 512 azimuthal uniformly distributed grid cells. The numerical convergence with the above resolution is found to be satisfactory, see Sect.\,\ref{sect:robust}. At the inner and outer boundaries the so called `damping boundary condition' are applied \citep{deValBorroetal2006mn}, where a solid wall is assumed with wave killing zones next to the boundaries.

\subsection{Development of a large-scale vortex}
\label{sec:vortex-devel}
In what follows we give a short description of the formation and evolution of a large-scale vortex in our base disc model ($M_\mathrm{disc}=0.01\,M_\odot$, $h=0.05$, $\alpha=0.01$, $\alpha_\mathrm{mod}=0.001$, $R_\mathrm{dze}=\mathrm{50\,AU}$, and $\Delta R_\mathrm{dze}=1\,\mathrm{AU}$). The major evolutionary stages of the surface mass density distribution is shown in Fig.\ref{fig:vortex_evol}, while the azimuthally averaged surface mass density profiles are shown in Fig.\,\ref{fig:densprofile_evol}.

To excite RWI, a small mass orbiting body (negligible to the local disc mass) was put in the disc at 50\,AU. The orbiting perturber excites low-$m$ mode density waves emerging from its Lindblad resonances, and these in turn trigger the RWI. Triggering RWI by this method is appropriate, because RWI is a low-$m$ mode instability. We observed, though, that for models where the RWI is strong -- e.g., for sharp dead zone edge models ($\Delta R_\mathrm{dze}=1\,\rm AU$), where the radial gradients of pressure and density are steep -- the numerical noise in the radial velocity component can be large enough to trigger the vortex formation.
This can be explained by the fact that the magnitude of radial velocity component is very small in Keplerian discs, thus a small numerical noise (which is in the order of $10^{-6}$ in the radial velocity component for our simulations done by {\small GFARGO}) is enough to excite RWI. For this reason, we investigated whether the mode of triggering RWI (by perturbing the radial velocity component or the surface mass density with superposition of sines or white noise; igniting density waves by a small mass planet orbiting at $\rm50\,AU$) has any effect on the vortex evolution. Studying carefully the results of our simulations, we have not found any significant differences, see details in Appendix\,B.

\begin{figure}
	\centering
	\includegraphics[width=\columnwidth]{fig3.pdf}
	\caption{Major evolutionary stages of the surface mass density evolution in our base model ($M_{disc}=0.01\,M_\odot$, $H/R=0.05$, $\alpha=0.01$, $\Delta R_\mathrm{dze}=1\,\mathrm{AU}$, $\alpha_\mathrm{mod}=0.001$, ). \emph{Panel A:} excitation of the RWI at $\sim4\times10^3\,\rm yr$ ($P_{50,1}=12$). \emph{Panel B:} only one large-scale vortex remains after $\sim1.5\times10^4\,\rm yr$ ($P_{50,1}=40$). \emph{Panel C:} vortex evolution results in the formation of horseshoe-shaped (azimuthal wave number $m=1$ mode) surface mass density distribution at $\sim10^5\,\rm yr$ ($P_{50,1}=300$). \emph{Panel D:} vortex wings are interlocked at $\sim3\times10^5\,\rm yr$ ($P_{50,1}=900$). \emph{Panel E:} ring-like structure forms with small vortices on the top of that at $\sim5\times10^5\,\rm yr$ ($P_{50,1}=1400$). \emph{Panel F:} reappearance of the horseshoe-shaped vortex at $\sim6\times10^5\,\rm yr$ ($P_{50,1}=1700$).}
	\label{fig:vortex_evol}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=\columnwidth]{fig4.pdf}
	\caption{Evolution of the azimuthally averaged surface mass density profile in the base model. The \emph{solid black curve} represents the initial surface mass density profile. The density profiles shown with \emph{coloured dashed curves} are calculated from the surface mass density distribution presented in Fig.\ref{fig:vortex_evol}.}
	\label{fig:densprofile_evol}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=\columnwidth]{fig5.pdf}
	\caption{Azimuthal spectral power of the surface mass density distribution against time measured in $P_{50,1}$ for model presented in Fig\,\ref{fig:vortex_evol}. The Fourier amplitudes are normalised by the total power measured in the $m=0-6$ modes.}
	\label{fig:power_01}
\end{figure}

According to our simulation, the gas accumulates soon ($\leq10^2\,\rm yr$) at slightly smaller distance than the outer edge of the dead zone (see, e.g., Fig.\,\ref{fig:densprofile_evol} solid line curve). The RWI is excited at $\sim4\times10^3\,\rm yr$ ($12$ orbits at $50\,\rm AU$ assuming a Solar mass central star, i.e. $P_{50,1}=12$) forming several anticyclonic vortices with azimuthal wavenumber $m=6$ (Fig.\,\ref{fig:vortex_evol}, Panel A). These vortices coalesce, while the density bump continuously strengthens (see, e.g., Fig.\,\ref{fig:densprofile_evol} dashed line curves). As a result, only one large-scale vortex ($m=1$ mode vortex) remains by $\sim1.5\times10^4\,\rm yr$ ($P_{50,1}=40$) (Fig.\,\ref{fig:vortex_evol}, Panel B). Both the azimuthal and radial extensions of vortex increase, producing a horseshoe-shaped density perturbation by $\sim10^5\,\rm yr$ ($P_{50,1}=300$) (Fig.\ref{fig:vortex_evol} panel C). The vortex wings finally locks by $\sim3\times10^5\,\rm yr$ ($P_{50,1}=900$) (Fig.\,\ref{fig:vortex_evol}, Panel D). Investigating the following evolution of large-scale vortex in the base model, we found that the horseshoe-shaped density perturbation disappears, and a ring like density perturbation develops at $5\times10^5\,\rm yr$ ($P_{50,1}=1400$) (Fig.\,\ref{fig:vortex_evol}, Panel E). On the top of the ring-like density perturbation, small-scale vortices form, which are also subject to the vortex coalescing process (see in the next paragraph). After $\sim10^5\,\rm yr$ ($P_{50,1}=300$) of evolution, the $m=1$ mode (horseshoe-shaped) vortex reappears (Fig.\,\ref{fig:vortex_evol}, Panel F).

To illustrate this vortex mode oscillation we calculated the azimuthal spectral power of the surface mass density distribution against time measured in the orbital period $P_{50,1}$ shown in Fig\,\ref{fig:power_01}. First, a periodic Fourier transform of the density in the ring (period of $2\pi$ along $\phi$-direction at the radius where the vortices are) was calculated. Then the Fourier amplitudes of the $m=0-6$ modes are calculated normalised by the total power in these seven modes. It is clearly seen that  the spectral power is mainly in the $m=1$ mode (black curve) throughout the whole simulation. However,  by $5\times10^5\,\rm yr$ ($P_{50,1}=1400$) of disc evolution, the spectral power in $m=1$ mode temporarily drops. When the $m=1$ mode vortex decays, the spectral power of $m=4,5$ modes abruptly rises, i.e. RWI is re-excited. These newly formed vortices are also subject to coalescing and form again an $m=1$ mode vortex  by $\sim10^5\,\rm yr$ ($P_{50,1}=300$). Thus, the typical vortex mode is $m=1$ through $\sim 1.8\times10^6\,\rm yr$ ($P_{50,1}=5200$) of vortex evolution, which is interrupted by short-lived higher mode ($1<m<5$) vortex configurations. Note that the spectral power in $m=0$ characterising the average surface mass density drops off immediately as RWI is excited.

As shown in Fig.\,\ref{fig:densprofile_evol}, the density maximum is pushed inward from its original location, i.e., the $m=1$ vortex is slightly shrinked due to the extra pressure caused by the density build-up. As a result, the fully-developed $m=1$ mode vortex forms at slightly smaller distance ($R\simeq40\,\mathrm{AU}$) to the central star than the distance where the viscosity reduction ($R_\mathrm{dze}=50\,\mathrm{AU}$) begins. It is also appreciable that the density enhancement inside the vortex eye increases with time. Finally, emphasise that the surface mass density of the disc is significantly lowered next to the dead zone edge. This density rarefaction continuously strengthens with time, forming a deep hole in the disc by $\sim1.5\times10^6\,\rm yr$ ($P_{50,1}=4250$).

\subsection{Robustness of large-scale vortex development}
\label{sect:robust}

To examine the robustness (sensitivity to disc properties) of the development of the large-scale vortex, we performed several hydrodynamic simulations. By harnessing the power of the GPU, each simulation required only a day of computer time. For this parameter study the boundary conditions, the disc mass, and the disc aspect ratio were the same as in the previously shown base model. In what follows we summarise briefly the results of these simulations.

\subsubsection{Numerical resolution}

To check whether our numerical solution is satisfactory with the original grid ($N_R = 256, N_\phi=512$), we recalculated the base model with higher resolution, namely $N_R=512, N_\phi=1024$. We found that the results are practically identical, even by frame-to-frame comparison. Hence, the spatial resolution of the base model (256$\times$512) is sufficient for the following simulations.

\begin{figure}
	\centering
	\includegraphics[width=\columnwidth]{fig6.pdf}
	\caption{Fully-developed large-scale vortices at $3.5\times10^5\,\rm yr$ ($P_{50,1}=1000$) in models using $\Delta R_\mathrm{dze}=1,2$, and $3\,\mathrm{AU}$ dead zone edge width shown in panels A, B, and C, respectively. The axisymmetric $m=0$ configuration for model $\Delta R_\mathrm{dze}=5\,\mathrm{AU}$ shown in panel D is just the result of the viscosity reduction at the dead zone edge stalling the mass inflow, with subsequent pileup of gas.}
	\label{fig:vortex_dR1235}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=\columnwidth]{fig7a.pdf}
	\includegraphics[width=\columnwidth]{fig7b.pdf}
	\caption{Azimuthally averaged surface mass density profiles at $3.5\times10^5\,\rm yr$ ($P_\mathrm{50,1}=1000$) in models using $\Delta R=1,2,3,\mathrm{and}\,5\,\mathrm{AU}$ dead zone edge width and $\alpha_\mathrm{mod}=0.001$ (\emph{upper panel}). Models with lower viscosity reduction $\alpha_\mathrm{mod}=0.01$ (\emph{lower panel}) show no significant departures to previous simulations.}
	\label{fig:densprofile_dRa}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=\columnwidth]{fig8.pdf}
	\caption{Azimuthal spectral power of the surface mass density distribution against time measured in $P_{50,1}$ for model $\alpha=0.01$, $\Delta R_\mathrm{dze}=3\rm \,AU$, $\delta\alpha=0.01$. The Fourirer amplitudes are normalised by the total power measured in the $m=0-6$ modes.}
	\label{fig:power_02}
\end{figure}

\subsubsection{Numerical precision}
Scientific calculations with conventional GPU cards (e.g., NVIDIA GTX-280 that we used) lacks the usual double precision at the moment.\footnote{Note that there exist scientific grade graphics cards in the NVIDIA Tesla series which are capable of double precision calculations.} To check the validity of our results provided by the single precision GPU we compare them to those of given by the CPU based {\small FARGO} code in the fastest evolving model ($\alpha=0.1$, $\Delta R_\mathrm{dze}=1\rm \,AU$, $\alpha_\mathrm{mod}=0.001$). The two runs give qualitatively the same result: the evolution of the large-scale vortex (growth rate, vortex coalescing, wing interlock, and oscillation) is found to be practically the same (see the details in Appendix A).

\subsubsection{Dead zone edge width}

\citet{Lietal2000} have found that the excitation of the RWI occurs later in time if the transition region of the viscosity reduction (i.e., the width of the dead zone edge) is wider. We confirmed this by setting the dead zone edge width $\Delta R_\mathrm{dze}$ between $1\,\mathrm{AU}-5\,\mathrm{AU}$ for models using $\alpha=0.01$ and $\alpha_\mathrm{mod}=0.01$. The excitation of the RWI occurs about $4$, $8.5$, and $20\times10^3\,\rm yr$ ($P_{50,1}=12,24$ and $57$) for models with $1$, $2$, and $3\,\rm AU$ dead zone edge width, respectively. However, in model with $\Delta R_\mathrm{dze}=5\,\rm AU$, the RWI is failed to be excited and only axisymmetric ($m=0$) density perturbation develops as a result of the viscosity reduction at the dead zone edge stalling the mass inflow, with subsequent pileup of gas. Thus, the RWI excitation occurs only in $\Delta R_\mathrm{dze}\leq2H$ models in a good agreement with \citet{Lyraetal2009}.

Figures\,\ref{fig:vortex_dR1235} and \ref{fig:densprofile_dRa} show the surface mass density contours and  the azimuthally averaged surface mass density profiles developed at $3.5\times10^5\,\rm yr$ ($P_\mathrm{50,1}=1000$) for various dead zone edge width models. Although the vortex structures are different (the larger the dead zone edge width, the smaller the azimuthal extension of the $m=1$ mode vortex, see Fig.\,\ref{fig:vortex_dR1235}, Panel A, B, C), the azimuthally averaged surface mass density profiles are very similar (Fig.\,\ref{fig:densprofile_dRa}, dashed line curves).

The vortex mode oscillation behaviour is completely absent in models with $\Delta R_\mathrm{dze}\geq2\,\mathrm{AU}$ during the whole simulation covering $\sim5\times10^6\,\rm yr$ ($P_{50,1}=14000$). Figure\,\ref{fig:power_02} show the azimuthal spectral power of the surface mass density distribution against time calculated as described before. For this case, the largest spectral power is in $m=1$ mode and remains nearly constant throughout the whole simulation after the initial vortex coalescing process.

\subsubsection{Depth of viscosity reduction}

We investigated the effect of the magnitude of $\alpha_\mathrm{mod}$ being in the range of $0.001-0.01$ for models with $\alpha=0.01$. We found that the time required to ignite the RWI, and to form the fully-developed large-scale vortex are independent of $\alpha_\mathrm{mod}$. Moreover, the magnitude of density contrast developed in the large-scale vortex is also independent of $\alpha_\mathrm{mod}$. As a consequence, the azimuthally averaged surface mass density profiles (Fig.\,\ref{fig:densprofile_dRa} calculated at $\sim3.5\times10^5\,\rm yr$ ($P_{50,1}=1000$) in weak viscosity reduction ($\alpha_\mathrm{mod}=0.01$, \emph{dashed line curves})) models are very similar to that of in models with stronger viscosity reduction ($\alpha_\mathrm{mod}=0.001$, \emph{solid line curves})). The life-time of the $m=1$ vortex, i.e., the vortex mode oscillation behaviour, is also found to be independent of $\alpha_\mathrm{mod}$: for sharp dead zone edge models the vortex mode oscillation do occur, while for thick dead zone edge models the $m=1$ mode vortex is long-lived.

\subsubsection{Magnitude of global viscosity}

We found that the vortex evolution is sensitive to the magnitude of the global viscosity $\alpha$ being in the range of $0.001-0.1$. The time required to excite RWI and to form the fully-developed large-scale $m=1$ mode vortex decreases with increasing $\alpha$. Both the radial extension of the large-scale vortex and the density contrast in the vortex increase with $\alpha$. These can be interpreted by the accelerated vortex evolution in models with high global viscosity.

For high viscosity models ($\alpha=0.1$), the inner disc becomes ever denser during the later ($t\gg5\times10^5\,\rm yr$) epochs, and the density contrast drops between the large-scale vortex and the inner disc. In consequence, the $m=1$ mode vortex survives only up to $\sim10^6\,\rm yr$ in these models. For low viscosity discs ($\alpha=0.001$), the vortex evolution slows down considerably, thus the excitation of RWI occurs later and the life-time of $m=1$ mode vortex can reach $\sim5\times10^6\,\rm yr$.

We found that for low viscosity discs the vortex mode oscillation is completely missing, even for sharp dead zone edge ($\Delta R_\mathrm{dze}=1\,\rm AU$) models. Contrary to this, for high viscosity discs, the vortex mode oscillation does occur even for thick dead zones ($\Delta R_\mathrm{dze}<5\,\rm AU$) models. Emphasise that the $m=1$ mode (horseshoe-shaped) vortex might survive up to the life-time of the disc for $\alpha=0.001$, independently of the dead zone edge width being smaller than the twice of the local disc height ($\Delta R_\mathrm{dze}\leq5\,\rm AU$).

\section{Submillimetre images}

In order to show the resemblance of the observations of \citet{Brownetal2009} and signatures of large-scale Rossby-vortices in protoplanetary discs, we calculated submillimetre images using the 3D radiative transfer code RADMC-3D\footnote{http://www.ita.uni-heidelberg.de/$\sim$dullemond/software/radmc-3d}. First, the temperature structure of the disc was determined in a Monte Carlo simulation then the images at 340\,GHz (880\,$\mu$m ) were calculated via raytracing. Submillimetre Array (SMA) observations were simulated in the Common Astronomy Software Applications package (CASA).\footnote{http://casa.nrao.edu/index.shtml}

Snapshots of surface mass density distributions shown in Fig.\,\ref{fig:model_frame} from the hydrodynamic simulations were converted into volume density assuming a vertically Gaussian density distribution with the same constant aspect ratio ($h=0.05$) as used in the hydrodynamic simulations. The geometrical dimensions of the sources as well as their distances were adapted from \citet{Andrewsetal2011}. The hydrodynamical frames had to be rotated and spatially scaled to approximately match the disc shapes and sizes on the measured submillimetre images. The used rotation angles were $330^\circ$, $320^\circ$, and $170^\circ$ clockwise, while the linear scaling factors were 2, 1.2, and 1.7 for LkH$\alpha$\,330, SR\,21N, and HD\,135344B, respectively. Since there are evidences for a fairly large gap in all discs \citep{Brownetal2009}, we reduced the dust density within a certain radius according to \citet{Andrewsetal2011}. Inside the dead zone, dust grains may grow and drift rapidly inwards leaving a reduced density region, i.e. a gap, behind. In anticyclonic vortices, however, dust grains get trapped and cannot drift inwards, which would qualitatively explain why the outer radius of the gap is close to the location of the RWI. Note that there is also evidence for the presence of an inner disc in these sources \citep{Fernandezetal1995,vanderPlasetal2008,Pontoppidanetal2008,Salyketal2009}. Since these inner discs have an outer radius probably smaller than the inner radius of our hydrodynamic simulations (4\,AU) we neglected them.

In the radiative transfer simulations we used the hydrodynamic grid in the radial and azimuthal directions while  N$_{\theta}$=100 grid points were used in the $\theta$ direction. The polar grid was a uniformly distributed, but we split the $-\pi/2\leq\theta\leq +\pi/2$ domain into two regions in order to resolve the disc vertically, while keeping the simulation computationally feasible. We placed 80 grid points in the $|\theta|\leq 0.35$ domain and 20 grid points were placed in the $|\theta|\geq 0.35$ region. In the thermal Monte Carlo simulations $2\times10^7$ photons were used to determine the temperature structure.

\begin{figure}
	\centering
	\includegraphics[width=\columnwidth]{fig9a.pdf}
	\includegraphics[width=\columnwidth]{fig9b.pdf}
	\caption{Hydrodynamic frames that were used to calculate the synthetic submillimetre images presented in Fig.\,\ref{fig:mm_images}. For modelling LkH$\alpha$\, 330 and SR\,21 frame shown in \emph{upper panel}, while for HD\,135344B frame shown in \emph{lower panel} were used. The frames were calculated in two different models: parameters were $R_\mathrm{dze}=50\,\mathrm{AU}$, $\Delta R_\mathrm{dze}=1\,\mathrm{AU}$, $\alpha=0.01$, $\alpha_\mathrm{mod}=0.001$ for \emph{upper panel}, and $R_\mathrm{dze}=50\,\mathrm{AU}$, $\Delta R_\mathrm{dze}=2\,\mathrm{AU}$, $\alpha=0.1$, $\alpha_\mathrm{mod}=0.01$) for \emph{lower panel}. To match the observed submillimetre images of \citet{Brownetal2009} frames had to rotate by $330^\circ$, $320^\circ$, and $170^\circ$ clockwise and spatially scale by 2, 1.2 and 1.7 for LkH$\alpha$\,330, SR\,21N, and HD\,135344B, respectively.}
	\label{fig:model_frame}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=9.2cm]{fig10a.pdf}
	\includegraphics[width=9.2cm]{fig10b.pdf}
	\includegraphics[width=9.2cm]{fig10c.pdf}
	\caption{Submillimetre images of LkH$\alpha$\,330 (\emph{top}), SR\,21N (\emph{middle}) and HD\,135344B (\emph{bottom}) calculated by the 3D radiative transfer code, from the frames of the hydrodynamical calculations presented in Fig.\,\ref{fig:model_frame}. The synthesised beam is shown in the lower left corners in blue. The position angle of the discs, i.e. the vortex eye position was rotated corresponding to the observation of \citet{Brownetal2009}.}
	\label{fig:mm_images}
\end{figure}

The dust in our model consisted of ``astronomical'' silicate grains \citep{WeingartnerDraine2001}, whose absorption and scattering cross sections were calculated using Mie theory from the optical constants. Dust grains had an MRN size distribution \citep{MRN1977}, $n(a)\propto a^{-3.5}$, sampled with 10 logarithmic spaced bins with minimum and maximum grain size of 0.1\,$\mu$m and 1000\,$\mu$m, respectively. Dust particles with different sizes have different temperatures, thus they are not thermally coupled to each other. The dust-to-gas ratio was assumed to be 0.01 and constant within the disc outwards of the outer radius of the gap. Inside the gap the dust-to-gas ratio should be lower than 0.01 if the gap is formed due to the growth and inward drift of dust grains. Finally, the simulation of the observation (Fourier transformation, u-v sampling, cleaning) was done with CASA. The resulting images are shown in Fig.\,\ref{fig:mm_images}. We also calculated SED from our synthetic disc model for all three sources and compared them to the observations (see Fig.\,\ref{fig:SEDs}).

Our synthetic model images (Fig.\,\ref{fig:mm_images}) resemble the most important features of the submillimetre continuum map of the three transition discs presented in \citet{Brownetal2009}. The disc emission is clearly non-axisymmetric, horseshoe-shaped, and it shows single arc-like maximum at position angle around 225$^\circ$ for LkH$\alpha$\,330, 120$^\circ$ for HD\,135344B, while two maxima are visible at $120^\circ$ and $270^\circ$ position angles for SR\,21N. The SEDs (Fig.\,\ref{fig:SEDs}) calculated from our disc models match the observations shortwards of about 1\,$\mu$m and longwards of about 10\,$\mu$m. Between 1 and 10\,$\mu$m, however, our models under predict the observed fluxes, which is a natural consequence of neglecting the inner disc emission. We concluded that our disc models with large-scale Rossby vortices fit the observations of the mid-infrared to submillimetre emitting regions of all three discs.

Given the fact that our model is vertically optically thin at submillimetre wavelengths the observed asymmetry in the continuum emission can arise from temperature, density perturbations or the combination of both. Since our disc is optically thin at 880\,$\mu$m the increase in the column density will increase the observed flux. The higher density regions, however, can be shielded from the heating radiation (higher optical depth) and therefore the temperature can decrease. The temperature perturbations in our model are less than a factor of two, which means the same factor for the flux perturbation caused by temperature effects, since we are in the Rayleigh-Jeans regime. The density perturbations on the other hand are much larger in amplitude (factor of 6--8, see, e.g., the red dashed curve in Fig.\,\ref{fig:vortex_evol}) and are therefore most likely responsible for the observed asymmetries in the image. A detailed study of these effects in non-axisymmetric discs will be the topic of a forthcoming paper.

\begin{figure*}
	\centering
	\includegraphics[width=5.8cm]{fig11a.pdf}
	\includegraphics[width=5.8cm]{fig11b.pdf}
	\includegraphics[width=5.8cm]{fig11c.pdf}
	\caption{The observed (black dots) and the theoretical SEDs calculated in the ray-tracing models of LkH$\alpha$\,330 (\emph{left}), SR\,21N (\emph{middle}), and HD\,135344B (\emph{right}). The failure to fit the NIR flux is due to the fact that, for computational feasibility reasons, the inner disc emission is neglected.}
	\label{fig:SEDs}
\end{figure*}

\section{Discussion \& Conclusion}

In this Paper we propose an interesting physical phenomenon that can be responsible for the large non-axisymmetric density perturbations revealed by submillimetre images of three transition discs of \citet{Brownetal2009}. Based on the results of our long-term, two-dimensional, grid-based, global hydrodynamic simulations (using compressible-gas model, neglecting the thermodynamical effects, and the disc self-gravity), we suggest that the observed asymmetries are caused by a large-scale vortex developed where the turbulent viscosity has a steep gradient, e.g. at the outer edge of the disc dead zone. Vortices are triggered by the Rossby wave instability (RWI), which later are subject to a coalescing process and form a large-scale anticyclonic vortex. Note that \citet{WolfKlahr2002} have previously shown that such large-scale anticyclonic vortices formed in the inner disc ($<10\,\rm AU$) can causes detectable intensity asymmetry in millimetre wavelengths ($345\,\rm GHz$ and $900\,\rm GHz$) which are proposed to be observed by the fully operable Atacama Large Millimeter/submillimeter Array (ALMA) with the longest available baseline.

After the initial excitation of RWI, a large-scale anticyclonic vortex develops by the coalescing of small-scale vortices at $R\simeq50\,\mathrm{AU}$, near the outer edge of disc dead zone, (Fig.\,\ref{fig:vortex_evol} and Fig.\,\ref{fig:densprofile_evol}), which grows both in size and density contrast. The time required to form a $m=1$ mode large-scale vortex is found to be $\sim1.5-4\times 10^4\,\rm yr$ ($P_{50,1}=40-115$), depending on the dead zone edge width $\Delta R_\mathrm{dze}$ (the sharper the dead zone edge, the faster the vortex evolution) and the magnitude of global viscosity (the larger the viscosity, the faster the vortex evolution). The $m=1$ mode vortex temporarily disappears and reforms through RWI re-excitation and subsequent vortex coalescing processes (Fig.\,\ref{fig:power_01}) with a cadence of $\sim10^5\,\rm yr$ ($P_{50,1}=300$) for sharp dead zone edge models ($\Delta R_\mathrm{dze}<2\,\mathrm{AU}$) and large global viscosity ($\alpha \geq 0.01$). However, the $m=1$ mode vortex remains stabilised for $\sim5\times10^6\,\rm yr$ ($P_{50,1}=14000$) for thick dead zone edge models ($2\,\mathrm{AU}\leq\Delta R_\mathrm{dze}<5\,\mathrm{AU}$) or small global viscosity ($\alpha\leq0.01$). For even thicker dead zone edge models ($\Delta R_\mathrm{dze}\geq 5\,\mathrm{AU}$), $m=1$ mode vortex does not form, only $m=0$ density perturbation develops. In short, we conclude that the horseshoe-shaped vortex ($m=1$ mode) formed due to the viscosity jump near the disc dead zone edge being at $\sim 50\,\rm AU$ can survive up to the disc life-time ($t\simeq5\times10^6\,\rm yr$), if the turbulent viscosity is small ($\alpha=0.001$), and the width of the dead zone edge is smaller than twice of the local disc height ($\Delta R_\mathrm{dze}<2H$).

To demonstrate the observability of the large vortices, we calculated synthetic submillimetre images with a 3D radiative transfer code using the surface mass density distribution provided by our hydrodynamic simulations. The distribution of the dust, emitting in the submillimetre wavelengths, is assumed to follow that of the gas, owing to the dust collecting effect of the pressure bump. The submillimetre intensity distributions given by our ray-tracing models -- fitting the observed SEDs in the stellar and submillimetre emitting domain (Fig.\,\ref{fig:SEDs}) -- show significant azimuthal asymmetries (Fig.\,\ref{fig:mm_images}). Comparing our submillimetre predictions to the observations of three transition discs LkH$\alpha$333, SR\,21N, and HD\,135344B \citep{Brownetal2009}, we found astonishing similarities. Note that there exist further examples of discs showing non-axisymmetric submillimetre images, such as LkCa\,15 \citep{Pietuetal2006}, GM\,Aur \citep{Hughesetal2009}, SR\,24S \citep{Andrewsetal2010}, and J1604-2130 \citep{WilliamsCieza2011} for which the same phenomenon might be responsible. Since the formation process of planetesimals and planetary embryos could be accelerated by the dust collecting and embryo trapping effect of the anticyclonic vortices caused by the vortical motions and the pressure bump \citep{BargeSommeria1995, KlahrHenning1997}, the non-axisymmetric submillimetre images presented by \citet{Brownetal2009} might be interpreted as snapshots of active planet-forming regions.

Note that there are other possible routes to generate vortices in protoplanetary discs, which also might result in $m=1$ vortex formation. For the presence of gap-carving giant planets the gas is expelled from the planet orbit and settles at the gap edges resulting in density maximum which is RWI-unstable and develops into vortices as shown by \citet{deValBorroetal2007}. Another way to generate vortices is the baroclinic instability, proposed by \citep{KlahrBodenheimer2003}, and recently confirmed by several works (e.g., \citealp{Petersenetal2007a,Petersenetal2007b,LesurPapaloizou2010,LyraKlahr2011}).

As it was mentioned before, not only enhanced dust coagulation, fragmentation should also happen inside the vortices if the observed non-axisymmetric features of submillimetre images are indeed due to dust trapped inside them. Otherwise, the submilimetre intensity would be weak due to the decreasing dust opacity of large grains. Thus, this would constitute indirect evidence of elliptic instability \citep{LesurPapaloizou2009} inside vortices leading to 3D turbulence in the vortex core, which \citet{LyraKlahr2011} quantified in compressible simulations to be at the level of 10 percent of the sound speed.

Here we have to mention some of the caveats of our calculations. Two-dimensional and three-dimensional turbulence behave quite differently regarding the energy cascade, which is form small to large $m$ (direct cascader) in 3D but can be form large to small $m$ (inverse cascade) in 2D. 3D global disc simulations of \cite{Meheutetal2010} revealed that although the RWI results in significant vertical streaming, the growth rate and evolution of vortices in 3D is very similar to that of in 2D. This may indicate that the vorticity generation by the RWI is strong enough to counter the 3D vortex stretching (the mechanism responsible for the direct cascade), so that the vortices survive and cascade inversely via viscous merging.

Our simulations cover several million years of the disc evolution, yet, the disc photoevaporation acting on the same time-scale was neglected for simplicity. Since photoevaporation process disperse the gas content of the disc, this might cause inward migration of the dead zone edge. Thus, photoevaporation process might play an important role in vortex evolution, which will be addressed in a forthcoming paper.

\citet{LinPapaloizou2011} have investigated the role of self-gravity on vortex instabilities generated at the outer edge of the gap opened by a giant planet in locally isothermal disc models. They found that the unstable modes shifts to higher azimuthal wavenumbers assuming self-gravity: the wavenumber of the most unstable mode increases as Q decreases. For weak self-gravity ($M_\mathrm{disc}\leq0.024\,M_*$), a single vortex forms through coalescing, which process is delayed as $Q$ decreases. However, for higher disc mass ($M_\mathrm{disc}=0.031\,M_*$) the final configuration is a vortex pair, $m=2$ mode. Since we used $M_\mathrm{disc}=0.01\,M_\odot$, our models are in the weak self-gravity regime, where the coalescing process leads to the formation of $m=1$ vortex.

Self-gravitating disc simulations require the inclusion of thermodynamics (the solution of energy equation) too \citep{Pickettetal1998,Pickettetal2000,Gammie2001}, otherwise the gas concentration would be over predicted. \citet{Lyraetal2009} found that assuming adiabatic equation of state for the gas, the pressure bump (that gives rise to RWI) sharpens considerably compared to their locally isothermal models, due to the high temperatures associated with compression. As a result, RWI might be triggered easier in adiabatic than in isothermal model. \citet{Petersenetal2007b} investigated the baroclinic feedback, in which a vortex enhances azimuthal temperature gradients to reinforce the vortex itself, in global disc simulations with anelastic-gas model. They concluded that this baroclinic feedback mechanism producing stronger local temperature gradients results in long-lived vortices if the radial temperature gradient is sufficiently large and either the radiative cooling time is small, or the thermal dissipation is high. On the other hand, vortices in the course of vortex coalescing may grow massive enough and reach the Jeans length, which provides a natural barrier to growth as shown by \citet{MamatsashviliRice2009}. In their shearing sheet simulations, assuming a simple cooling law with a constant cooling time and including self-gravity, the vortex starts to shock and radiate away the excess vorticity around the Jeans length. As a result, the vortices can not grow above the Jeans length and short-lived, if there is no continuous vorticity production like for RWI, for which case the shear is continuously being converted into vorticity \citep{Lyraetal2009}. In our simulations, however, the vorticity is not generated by the local advective heat transport enhancing the azimuthal temperature gradients (see, e.g., \citealp{Petersenetal2007b}), but the steep pressure bump and azimuthal density perturbations formed near the dead zone edge. In short, investigating the vortex evolution and longevity in a self-gravitating global disc model, considering non-isothermal and compressible-gas model, would be important, and will be addressed in a forthcoming paper.

Finally, worth noting that multiple viscosity jumps might exist in the disc, e.g., at the inner edge of the dead zone or at the snow line \citep{KretkeLin2007}, which could also result in the formation of large-scale vortices. In a non-magnetised disc its self-gravity causes outward angular momentum transport via quasi-steady self-gravitating turbulence, which results in the generation of effective viscosity \citep{Pringle1981}. Since only the outer disc ($R\ga$100\,AU) is expected to become gravitationally unstable in a standard protoplanetary discs, the self-gravity caused effective viscosity has a gradient there. In such systems the large-scale vortex formed beyond $\sim$100\,AU might result in vortex aided planet formation process at extremely high distances, such as observed around Fomalhaut \citep{Kalasetal2008mn}. As a conclusion, planets might be assembled in diverse fashion by these ``planetary factories'' located at different distances to their host star. In protoplanetary discs where the RWI is excited at smaller radii than we investigated, the development of large-scale vortices might lead to the formation of a planetary system similar to the Solar System.

\section*{Acknowledgments}

This research has been supported in part by DAAD-PPP mobility grant P-M\"OB/841/ and `Lend\"ulet' Young Researcher Program of the HAS. We are grateful to F. Masset who helped us making some necessary modification in his {\small GFARGO} code. We thank M.R. Hogerheijde and Ch. Brinch for the help and discussion on CASA data processing. We also thank to the referee Wladimir Lyra for thoughtful comments that helped to significantly improve the quality of the paper.

\bibliographystyle{mn2e}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

Finding new fundamental physics beyond the standard model and general relativity has
long been a topic of great interest. One way to identify definitive evidence of new
physics is to look for exotic phenomena that
simply cannot occur in the standard physical
theories. One such exotic possibility that has provoked significant interest is
having time-dependent fundamental
constants~\cite{ref-fischer,ref-peik,ref-tzanavaris,ref-gould,ref-cingoz,ref-blatt,ref-murphy,ref-tobar}.

The possibility of time-varying fundamental constants opens up many new lines of
inquiry. The most commonly studied theory with varying constants is one
with a time-dependent fine structure constant $\alpha$. A quantum
approach to this subject can be found in~\cite{ref-ferrero1}, while other approaches
treat the problem rather phenomenalistically (e.g.~\cite{ref-flambaum1}).
Another possibility that has been fairly well studied is a change in the quantum chromodynamics scale
$\Lambda_{QCD}$ and therefore the electron-proton mass ratio.

If some of the fundamental constants in nature are allowed to change over time,
the effects of these changes might be detected. Searches for such temporal variations
entail measurements of the same experimental observable at different points in time,
to see if there is any change. However, since the changes involved
should be tiny, laboratory experiments searching for them demand very high
precision. Less precision is required if large time scales are
involved; a very slow time variation can be observed if the times used
to detect it are long enough. Thus, cosmological observations, which give us a view of
how physics operated millions or billions of years in the past, seem to be good
candidates to test for time variations in fundamental constants.

The purpose of this paper is to study a possible change in the Fermi
constant $G_{F}$. This quantity is defined to be $G_{F}=
\frac{\sqrt{2}}{8}\frac{g^{2}}{M_{W}^{2}}$, with $g$
the standard model's $SU(2)_{L}$ gauge coupling constant and $M_{W}$ the mass of the
$W$ boson. Since the vector boson mass is $M_{W}=\frac{1}{2}gv$, $G_{F}$ is really a measure of
$v$, the Higgs field vacuum expectation value.
The masses of fundamental fermions depend on $v$ and the values
of the fermion-Higgs Yukawa couplings $\lambda_f$, with $m_f=\lambda_f v$.
So a changing value of $G_F$ would also affect the fermion masses unless the
couplings $\lambda_f$ also changed to compensate. Since there is no evidence of any time variation
in $m_e$ in atomic spectra, we shall assume that such cancellation between the time
variation of $v$ and $\lambda_f$ occurs, even though that is not a particularly natural scenario.

$G_{F}$ also governs the strength of
charged current weak interactions at low energies.
It can be measured with relatively high precision by studying the
lifetimes of muons, and its accepted value is
$G_{F}=(1.6637\pm0.0001)\times 10^{-5}$
GeV$^{-2}$~\cite{ref-chitwood,ref-barczyk}.
However, laboratory measurements
of decay rates are not precise enough to offer a particularly promising way of
measuring $\dot{G}_{F}$. On the other hand, the light coming from supernovae
explosions can provide a useful way to study weak interactions that occurred at much
earlier times.

In Section~\ref{sec-SN}, we shall discuss the physics of supernovae and how their
light curves are related to the Fermi coupling constant $G_{F}$.
Section~\ref{sec-method} covers the details of the calculational methodology used to
relate the supernova light curve data to $\dot{G}_{F}$. In Section~\ref{sec-data},
we apply this methodology to data from the Harvard-Smithsonian Center
for Astrophysics (CfA) Supernova
Archive~\cite{ref-riess,ref-jha} and state our results. Section~\ref{sec-concl}
presents conclusions and outlook.

\section{Supernova Structure}
\label{sec-SN}

Supernovae are catastrophic variable stars that, having previously produced energy
through steady-state fusion processes, explode and release huge additional amounts
of energy. The amount of energy released is so vast that even cosmologically
distant supernovae can be easily detected in the sky. Normal stars avoid
gravitational collapse because nuclear fusion processes are
constantly generating energy in their cores, where the
high temperature and pressure enable light nuclei to
overcome their electrostatic repulsion and fuse. A large red giant star may burn a
significant portion of its mass into isotopes of iron, beyond which further fusion is
no longer exothermic. However, smaller stars, without thick external layers to
compress the core to the
high temperatures and pressures needed to produce iron, may cease burning when the
core is mostly carbon and oxygen, without having extracted all of this fuel's
available energy. In either case, the core of the star settles down to become
a white dwarf, supported against gravitational collapse by the degeneracy
pressure of its electrons.

If the mass of a white dwarf star is larger than approximately $1.4M_\odot$
(the Chandrasekhar mass $M_{Ch}$), the electron degeneracy pressure can support the star
against gravitational collapse. However, if the mass of the star grows larger than
this limit, the electron pressure can no longer prevent a collapse. Such
a collapse triggers the release of a tremendous amount of energy and the star
becomes a supernova. Iron white dwarfs at the cores of massive stars collapse to
form even more compact objects---neutron stars or black holes. However, of interest
here are the type Ia supernovae produced from carbon-oxygen white
dwarfs. These stars may accrete mass from (or possibly merge with) their partners
in binary systems, until the white dwarf's mass becomes sufficient to trigger a
supernova. When these stars exceed $M_{Ch}$, they begin to collapse.
However, the collapse phase is short-lived;
as the star contracts, there is an increase in temperature and density in
the core of the star, which allows the ignition of carbon fusion. Within a short
period of time, a huge amount of fusion energy is released, blowing the star apart
and increasing its luminosity to enormous values. (The total luminosity may be that
of millions of ordinary suns).

In this paper, we shall be analyzing the light curves
of type Ia supernovae. The near equality of the masses of the white dwarf progenitors
gives them all similar energy outputs and
makes them extremely important standard candles in cosmology.
Type Ia supernovae all seem to follow a characteristic light curve. (Some corrections to
the light curves are required in order for the population of type Ia supernovae to act as
genuinely standard candles, but these are not important for our present purposes.)
The relative uniformity of this curve structure for different stars also
makes them a very useful tool for studying certain physical properties.

The light curve is a plot of the apparent magnitude $m=-2.512\log\frac{f}{f_{0}}$
of the source versus time; $f$ is the radiation flux from the source, and $f_{0}$ is
a reference value of the flux defining 0 magnitude. Using different kinds
of filters, which cover different parts of the spectrum, several different magnitudes
may be measured for a single source. There is a separate light curve for each band,
but the curves for different colors show similar behavior, which
is again consistent across different sources.

At early times, the brightness of a supernova increases rapidly until it reaches a
maximum, which is followed by two regions in which $m$ decreases linearly
(corresponding to exponential decays in the flux).
Since type Ia supernovae are ultimately powered by uncontrolled nuclear fusion,
they produce large quantities of the heavy elements,
especially those elements with binding energies close to the maximum
value. Specifically, each supernova produces a great deal of the
double magic nucleus $^{56}$Ni. This nucleus decays by electron capture
(EC) to $^{56}$Co with a lifetime of 6.077 d, and this decay provides much
of the energy released by the supernova during its earlier stages. The first
linear region along the light curve corresponds to the period during which
$^{56}$Ni decay is the dominant source of heating. The daughter nucleus
$^{56}$Co also decays by EC to $^{56}$Fe, with a lifetime of 77.27 d. At
later times (in the second, more slowly decaying linear region), it is this decay
that primarily fuels the supernova's luminosity. This region, whose extent is
determined by the decay rates of $^{56}$Ni and $^{56}$Co, can be found between
approximately 35 and 100 days after the maximum magnitude is reached.
During this time, some months after a supernova reaches its peak, the light curve
approximately tracks the quantity of $^{56}$Co still present. The falloff in
the luminosity can therefore be used to measure the lifetime of this
nuclide at the time the supernova occurred. This gives direct information
about the value of several standard model parameters at the time of the explosion's
occurrence, including $G_F$. This was first suggested in ~\cite{scher-sper}, when much less data
was available.

An EC decay rate depends on fundamental constants in the following
fashion. Being a tree-level weak decay, the rate is naturally proportional to
$G_{F}^{2}$.
However, it also depends on the probability that an $s$-state electron will be found
close enough to the nucleus to be taken up by a proton. The nuclear volume
is proportional to $\Lambda_{QCD}^{-3}$.
The square of the electron wave function in the vicinity of the
nucleus is proportional to $a_{0}^{-3}$, where $a_{0}=\frac{1}{m_{e}e^{2}}$ is
the Bohr radius.
The inverse $\Gamma$ of the mean life of a nucleus that is unstable against EC
is therefore proportional to the combination of fundamental quantities
$G_{F}^{2}\alpha^{3}\left(\frac{m_{e}}{\Lambda_{QCD}}\right)^{3}$. (The rate also
depends on the $Q$-value of the decay, which is further dependent on
$\alpha$ and $\frac{m_{e}}{\Lambda_{QCD}}$, in more complicated fashion.)
However, changes in $\alpha$ and $\frac{m_{e}}{\Lambda_{QCD}}
\propto\frac{m_{e}}{m_{p}}$ are much easier to measure than would
be changes in $G_{F}$. $\alpha$ and $\frac{m_{e}}{m_{p}}$ affect atomic
energy levels directly, so their time derivatives can be constrained by comparing
spectra taken at different times. Cosmological and laboratory measurements have
constrained these dimensionless quantities to vary only by parts in $10^{16}$ or
less per year.

A supernova's light curve and total luminosity are also affected by the
value of the gravitational constant $G_{N}$. The amount of $^{56}$Ni
produced in the explosion is proportional to the total mass of the
collapsing star, $M_{Ch}\propto G_{N}^{-3/2}$~\cite{ref-gaztanaga, ref-riazuelo}.
Observations of supernovae have been used to place bounds on the
variation of $G_{N}$ at the
$\left|\frac{\dot{G}_{N}}{G_{N}}\right|\lesssim10^{-11}$ yr$^{-1}$ level.
We choose here to look at a model that neglects variations in $G_{N}$, although an
analysis that considered simultaneously varying $G_{N}$ and $G_{F}$ would
also be interesting. We shall also neglect any possible time dependences of
$\alpha$ and $\frac{m_{e}}{\Lambda_{QCD}}$, in light of the strong bounds
on such time dependences.
This means there is a simple relationship between $\dot{\Gamma}$ and $\dot{G}_{F}$,
\begin{equation}
\frac{\dot{G}_{F}}{G_{F}}=\frac{1}{2}\frac{\dot{\Gamma}}{\Gamma}.
\end{equation}
Bounds on changes in the $^{56}$Co decay rate may thus
be interpreted as bounds on a change in the Fermi constant.

\section{Method}
\label{sec-method}

In the light curve of a type Ia supernova, the slope of the second linear region
is related to the decay rate of $^{56}$Co.
The main goal of this paper is to compare this slope across different
supernovae of different ages. The slope for each one is related to the Fermi
constant at a particular time. We will make this comparison using data from the
CfA Supernova Archive.

Ideally, one could simply measure the slopes of the supernova light curves in their
$^{56}$Co-dominated regions and then determine whether the slopes are correlated with
the ages of the events. However, several corrections must be made before this can
be done.

The light curves do not give a pure measurement of the $^{56}$Co lifetime;
what is observed is the luminosity of the supernova remnant. The
EC decay is the object's main energy source, but the decay energy is
moderated through the complicated interactions of the expanding gas, so
the luminosity does not track perfectly with the $^{56}$Co decay profile.
The details of the remnant dynamics differ a bit from supernova to supernova.
It is not possible to correct for this lack of complete uniformity in the
structures, and this will necessarily be a source of random error.

Some of the corrections may be handled more simply, however.
This is the case for the correction arising from the well-known Doppler shift
associated with
the expansion of the Universe. Since the inverse of the mean life,
$\Gamma$, measures how fast $^{56}$Co is decaying (and therefore,
the frequency at which EC processes are occurring), the observed
value of $\Gamma$ will be subject to the same cosmological redshift as any other
frequency. In terms of the red shift $z$, the measured value $\Gamma$ on Earth and
the proper decay rate $\Gamma_{0}$ at the exploding star  are related by
\begin{equation}\label{a1}
\Gamma_{0}=(1+z)\Gamma.
\end{equation}
The physical effect of the redshift is to slow down the decay of the light curve.

We shall denote by $\dot{m}$ the slope of the second linear region of
each particular light curve, as observed on the Earth. This slope and the $^{56}$Co decay rate
should be directly proportional, so we have
\begin{equation}\label{a2}
\Gamma_{0}=(1+z)\Gamma\equiv A(1+z)\dot m= A\dot{m}_{0},
\end{equation}
where $\dot{m}_{0}$ is the proper decay rate of the supernova light
curve, and $A$ is a proportionality constant related to the complex physics of the explosion.
The quantity $\dot{m}_{0}$ is proportional to the square of the Fermi coupling
constant at the time $t_{0}$ that the explosion occurred,
$\dot{m}_{0}\propto\Gamma_{0}\propto G_{F}^{2}(t_0)$. If $G_{F}$ has changed
over time, we would expect to see a change in the decay rates $\dot{m}_{0}$
observed for different supernovae. Comparing
the light curves for supernovae of different ages gives a good way to measure
a change in the Fermi coupling constant.

The light from a particular supernova is emitted and collected at the times
$t_{0}$ and $t$, respectively. Denoting  the scale factor of the Universe by $a(t)$
and the Hubble constant by $H$, we have the standard relations
\begin{eqnarray}\label{a4}
z & = & \frac{a(t)}{a(t_{0})}-1 \\
\label{eq-H}
H & = & \frac{\dot{a}}{a}.
\end{eqnarray}
The factor $H$ converts derivatives with respect to $z$ into ones with respect to
time.
Based solely on the relations (\ref{a4}) and (\ref{eq-H}), one may be tempted to
conclude that
$\frac{\dot\Gamma_0}{\Gamma_0}=-\frac{H}{\langle\dot m_0\rangle}\frac{d\dot m_0}{dz}$
(where $\langle\dot m_0\rangle$ is the average value of $\dot m_0$ as computed
from the different sources). However,
a new complication arises because of the differences in the behaviors of the
$\dot m_0$
observed in different wavelength bands, and this relation for
$\frac{\dot{\Gamma}_{0}}{\Gamma_{0}}$
is subject to further corrections.
Looking at the light from different color regions, we see
slightly different slopes associated with the different bands. At first glance this
might not seem to be a problem,
since the value of the proportionality constant $A$ defined in eq. (\ref{a2}) could just be
redefined for each band. The value $\frac{\dot\Gamma_0}{\Gamma_0}$ does not
depend on $A$, so it would not affect the final result.

Nevertheless, the fact that $A$ does not remain constant
across the various bands introduces a new correction. To understand
why a further correction is necessary, let us consider light that is
detected on Earth using a red filter. This light has been redshifted
by the expansion of the Universe. If the slopes for the different bands were the
same this would not be a problem, but the red light that is detected
actually had a shorter wavelength when it was emitted during the explosion.
All the red light that is detected on Earth was some amount bluer at the time of
its emission; and we should ideally be comparing light from different sources that
had the same wavelength at the time of emission---not at the time of observation.
If the light curve steepens with increasing wavevector $k$, the redshift would
transform this into a steepening (within a given wavelength band) as a function of
$z$, because the older light originally had a higher frequency. What is required is
similar to a light curve $K$ correction~\cite{ref-hubble}. However, rather than using tabulated
$K$ corrections, we shall, for this low redshift sample of supernovae, calculate
our correction directly from the data set. This is justified because the
technique minimizes our reliance on the results of unrelated fits and because we
shall find that the correction it produces is relatively small.

To account for the redshift's effect, we
suppose that the dependence of $\dot m_0$ on $k$ is weak,
and make a linear approximation
\begin{equation}\label{a4a}
(1+z)\dot m_X=(\dot m_X\!)_0+\frac{d\dot m_0}{dk}(k_0-k_X)=(\dot m_X)_0+\frac{d\dot m_0}{dk}(zk_X).
\end{equation}
$\dot m_X$ is the observed slope in the band $X$ with mean wavenumber $k_X$.
Once corrected by $(1+z)$, $\dot{m}_{X}$ gives the proper decay rate $\dot m_0$,
but for a band
centered at a higher wavenumber $k_0=(1+z)k_X$. The $z$ dependence of $k_0$ generates
a spurious dependence on $z$ in the observational data, which we
must subtract away.

Unlike the simple $(1+z)$ factor relating $\dot m$ to $\dot m_0$ (which is an
exact correction), there is
no particular reason for this effect to be linear in $z$. The wavelength dependence
of the proportionality constant $A$ is a consequence of the complicated dynamics of
the supernova remnants; its detailed form is not known.
If $\frac{d\dot m_0}{dk}$
is itself a function of $z$, there will be $\mathcal O(z^2)$ corrections to (\ref{a4a}).
At the same order, it would also be necessary to correct the average $\langle(\dot m_X)_0\rangle$
to account for the fact that we are actually averaging decay rates in different color regions.
We are neglecting these possibilities (which would correspond to corrections that are
higher order in $z$ in our final formula for $\frac{\dot{\Gamma}_{0}}{\Gamma_{0}}$),
because the redshifts in the sample we shall study
range only from $3.9\times 10^{-3}$ to $2.4\times 10^{-2}$. However, these corrections can be
computed, and working with sources at higher $z$ would require us to have a better
knowledge of how $\dot m_0$
depends on $k$. Yet there is a limit to the precision that can be derived from more detailed
models, since there are data available for only a handful of wavelength
bands (five in this case), and data points at multiple frequencies would be required
to map out the detailed dependence of $A$ on $k$.

In principle, we could try to correct the values of $\dot m_0$ for each supernova
individually, by studying the dependence of its light curves on $k$.
We could then use the corrected values obtained separated for each source
to compute $\frac{1}{\langle\dot m_0\rangle}\frac{d\dot m_0}{dz}$.
However, this would require additional fitting for each supernova,
which may introduce additional errors and which cannot be practically implemented for many of
the sources because of the limited data available. We
shall therefore make only a single calculation of
$\frac{d\langle\dot m_0\rangle}{dk}$, using the combined data from all the sources.

Using eq. (\ref{a4a}) we find our final formula for $\frac{\dot\Gamma_0}{\Gamma_0}$,
\begin{equation}\label{a4c}
\frac{\dot\Gamma_0}{\Gamma_0}=-\frac{H}{\langle\dot m_0\rangle}\left[\frac{d\dot m_0}{dz}
-\frac{2\pi}{\lambda}\left(\frac{d\langle\dot m_0\rangle}{dk}\right)\right],
\end{equation}
including the linear correction to eliminate the spurious $z$ dependence caused by a
nonzero $\frac{d\dot m_0}{dk}$. Evaluating this requires multiple fits. The order of the
fitting procedures
is as follows. First, $\dot m_0$ is determined for each supernova in each band separately.
These are then fitted as a linear function of the redshift $z$,
giving both an average $\langle\dot m_0\rangle$ and a slope $\frac{d\dot m_0}{dz}$ for each band.
From the fit of the five average values $\langle\dot m_0\rangle$
versus $k$, the parameter $\frac{d\langle\dot m_0\rangle}{dk}$ can be extracted.
Then $\frac{\dot{\Gamma}_{0}}{\Gamma_{0}}$ can be evaluated according to
eq.~(\ref{a4a}); this evaluation can be done separately in each of the five bands,
although the single fitted value for the quantity
$\frac{d\langle\dot{m}_{0}\rangle}{dk}$ is common to all of them.

The overall negative sign in eq.~(\ref{a4c}) merits some explanation. Because the light curves are
decaying in the region of interest, $\dot{m}_{0}$ is always negative. If $\dot{m}_{0}$ is
more negative at greater values of $z$ ($\frac{d\dot{m}_{0}}{dz}<0$ and
$\frac{1}{\langle \dot{m}_{0}\rangle}\frac{d\dot{m}_{0}}{dz}>0)$, the decay rate was
faster at earlier times. This means $\Gamma_{0}$ and thus $G_{F}$ are decreasing as
time passes, leading to the negative sign in eq.~(\ref{a4c}).

\section{Data Analysis}
\label{sec-data}

\begin{figure}
\begin{center}
\includegraphics[scale=1.2]{gr1998aqCol}
\caption{Light curves for SN 1998aq. The magnitudes for the different
bands are $m_U$, $m_B-2$, $m_V-4$, $m_R-6$,
and $m_I-8$. The black lines represent linear fits between
50 and 130 days after the peak.}
\label{f1}
\end{center}
\end{figure}

Using the method we have discussed, we can now place a bound on the time
variation rate of $G_{F}$ using data from the CfA Supernova Archive.
For purposes of illustration, we shall initially focus on one particular
supernova, SN 1998aq~\cite{ref-riess2}, for which $z=3.95\times10^{-3}$.
Its light curves are shown in fig.~\ref{f1}, plotting the supernova's magnitude as
a function of time $t$ in days.
The linear behavior arising from the $^{56}$Co decay
can be seen in a region covering approximately 50--130 days after the maximum magnitude
was recorded.
Although the data has small errors related to the measurements, these
have been omitted for clarity of presentation. To find the slopes in the different bands, we used a weighted linear fit in the region of interest, minimizing
the $\chi^2$ per degree of freedom. The results of these fits are also shown in
fig.~\ref{f1} and the fitted slopes are listed in table~\ref{t1}. In order to
transform these results into the reference frame of supernova, we just use the
relation $\dot{m}_{0}=(1+z)\dot{m}$.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|}
\hline
\rule{0pt}{2.4ex}\!\!
Band &  $\dot{m}$ (10$^{-4}$\,d$^{-1}$)\\
\hline
\rule{0pt}{2.5ex}\!\!
  U  &  $-214\pm7$\\
  B  &  $-138\pm2$ \\
  V  &  $-261\pm4$\\
  R  &  $-321\pm4$\\
  I  &  $-415\pm6$\\
\hline
\end{tabular}
\end{center}
\caption{Fitted values of the slope $\dot{m}$ for SN1998aq, in units of $10^{-4}$
d$^{-1}$.}
\label{t1}
\end{table}

We have applied the same procedure to multiple sources. Not all sources in the
CfA Supernova Archive had sufficient data to be fitted reliably, and not all bands
were represented for all sources.
In table \ref{t2}, we summarize the results found for the different
supernovae we analyzed.
The slopes $\dot{m}$ are in the observer reference frame on Earth (no $1+z$
factor included). There was one further complication associated with fitting some
of the light curves. For some sources, there is limited data at early times, and
it was not possible to resolve the date at which the light curve achieved its
absolute maximum. In many cases, we could identify the proper linear region by
inspection, but if there was any doubt whether a given data point should be
included, the fitted region was chosen so as to produce the best linear fit.
\begin{table}[h]
\begin{center}
\begin{scriptsize}
\begin{tabular}{|l|r|l|l|l|l|l|}
\hline
\rule{0pt}{2.3ex}
Source   &        $ 10^5\,z $         & $\dot m_U$ ($10^{-4}\,$d$^{-1}$) & $\dot m_B$ ($10^{-4}\,$d$^{-1}$)
& $\dot{m}_V$ ($10^{-4}\,$d$^{-1}$) & $\dot m_R$ ($10^{-4}\,$d$^{-1}$) & $\dot m_I$ ($10^{-4}\,$d$^{-1}$) \\
\hline
\rule{0pt}{2.3ex}\!\!
1995bd& $160 $ & \hspace{0.95cm}--         & \hspace{0.3cm}$-171\pm7 $ & \hspace{0.3cm}$-271\pm3 $ & \hspace{0.95cm}--         & \hspace{0.95cm}--          \\
1998aq& $395 $ & \hspace{0.3cm}$-214\pm8 $ & \hspace{0.3cm}$-138\pm2 $ & \hspace{0.3cm}$-261\pm4 $ & \hspace{0.3cm}$-321\pm4 $ & \hspace{0.3cm}$-415\pm6  $ \\
1994ae& $427 $ & \hspace{0.95cm}--         & \hspace{0.3cm}$-145\pm2 $ & \hspace{0.3cm}$-224\pm4 $ & \hspace{0.3cm}$-286\pm7 $ & \hspace{0.3cm}$-359\pm8  $ \\
1995al& $514 $ & \hspace{0.95cm}--         & \hspace{0.3cm}$-152\pm8 $ & \hspace{0.3cm}$-246\pm4 $ & \hspace{0.3cm}$-311\pm3 $ & \hspace{0.3cm}$-366\pm1  $ \\
1995D & $654 $ & \hspace{0.95cm}--         & \hspace{0.3cm}$-131\pm5 $ & \hspace{0.3cm}$-248\pm6 $ & \hspace{0.3cm}$-333\pm35$ & \hspace{0.3cm}$-384\pm24 $ \\
1999gh& $767 $ & \hspace{0.3cm}$-182\pm20$ & \hspace{0.3cm}$-175\pm5 $ & \hspace{0.3cm}$-251\pm5 $ & \hspace{0.3cm}$-322\pm12$ & \hspace{0.3cm}$-299\pm33 $ \\
2000cx& $793 $ & \hspace{0.3cm}$-226\pm4 $ & \hspace{0.3cm}$-187\pm2 $ & \hspace{0.3cm}$-257\pm8 $ & \hspace{0.3cm}$-312\pm8 $ & \hspace{0.3cm}$-347\pm19 $ \\
1999ac& $949 $ & \hspace{0.3cm}$-240\pm6 $ & \hspace{0.3cm}$-175\pm14$ & \hspace{0.3cm}$-231\pm11$ & \hspace{0.3cm}$-313\pm13$ & \hspace{0.3cm}$-317\pm40 $ \\
1998bp& $1042$ & \hspace{0.95cm}--         & \hspace{0.3cm}$-187\pm13$ & \hspace{0.95cm}--         & \hspace{0.95cm}--         & \hspace{0.3cm}$-289\pm42 $ \\
1998es& $1056$ & \hspace{0.95cm}--         & \hspace{0.3cm}$-187\pm16$ & \hspace{0.3cm}$-314\pm15$ & \hspace{0.3cm}$-376\pm7 $ & \hspace{0.95cm}--          \\
1995E & $1156$ & \hspace{0.95cm}--         & \hspace{0.95cm}--         & \hspace{0.3cm}$-253\pm20$ & \hspace{0.3cm}$-321\pm17$ & \hspace{0.3cm}$-441\pm23 $ \\
1999dq& $1432$ & \hspace{0.3cm}$-254\pm1 $ & \hspace{0.3cm}$-163\pm15$ & \hspace{0.3cm}$-250\pm17$ & \hspace{0.3cm}$-332\pm18$ & \hspace{0.3cm}$-387\pm33 $ \\
1999aa& $1443$ & \hspace{0.3cm}$-210\pm2 $ & \hspace{0.3cm}$-130\pm8 $ & \hspace{0.3cm}$-248\pm7 $ & \hspace{0.3cm}$-316\pm6 $ & \hspace{0.3cm}$-414\pm7  $ \\
1993ae& $1904$ & \hspace{0.95cm}--         & \hspace{0.95cm}--         & \hspace{0.3cm}$-231\pm41$ & \hspace{0.3cm}$-316\pm24$ & \hspace{0.3cm}$-428\pm123$ \\
1994M & $2296$ & \hspace{0.95cm}--         & \hspace{0.3cm}$-149\pm7 $ & \hspace{0.3cm}$-294\pm18$ & \hspace{0.3cm}$-399\pm2 $ & \hspace{0.3cm}$-454\pm1  $ \\
1995ak& $2300$ & \hspace{0.95cm}--         & \hspace{0.3cm}$-243\pm20$ & \hspace{0.3cm}$-310\pm12$ & \hspace{0.3cm}$-418\pm25$ & \hspace{0.3cm}$-489\pm57 $ \\
\hline
\end{tabular}
\end{scriptsize}
\end{center}
\caption{Fitted values of $\dot{m}$, in units of $10^{-4}$ d$^{-1}$. Sources and
bands for which the data were insufficient are marked ``--''.}
\label{t2}
\end{table}

The data available for some supernovae and some wavelength bands were not of
sufficient quality to be useful. We omitted those sources for which the standard
error in the fit of $\dot{m}$ was larger than 0.015 d$^{-1}$ and those with fewer
than three data points in the region of interest. The former was indicative of
substantial deviations from linear behavior (that is, nonexponential decay in
the luminosity), and the later simply represented too few points to be useful.

The $z$ dependence of $\dot{m}_{0}$ for each wavelength band is shown in
fig.~\ref{f2}. It is evident that $\frac{d\dot{m}_{0}}{dz}$ is negative in each of
the five
distinct wavelength regimes. The values of the slopes, as well as some other
relevant quantities, are collected in table~\ref{t5}.

\begin{figure}
\begin{minipage}[b]{0.5\linewidth}
\centering
\includegraphics[scale=0.6]{g-Ultcol}
\end{minipage}
\hspace{0.5cm}
\begin{minipage}[b]{0.5\linewidth}
\centering
\includegraphics[scale=0.6]{g-Blucol}
\end{minipage}
\begin{minipage}[b]{0.5\linewidth}
\centering
\includegraphics[scale=0.6]{g-Viscol}
\end{minipage}
\hspace{0.5cm}
\begin{minipage}[b]{0.5\linewidth}
\centering
\includegraphics[scale=0.6]{g-Redcol}
\end{minipage}
\begin{minipage}[b]{0.5\linewidth}
\centering
\includegraphics[scale=0.6]{g-Infcol}
\end{minipage}
\caption{Plots of $\dot{m}_{0}$ versus $z$ for the U, B, V, R, and I bands
(denoted by colors as in fig.~\ref{f1}). In each plot, the dashed horizontal
line marks the mean $\langle\dot{m}_{0}\rangle$; the substantial difference in
this mean for different wavelength regions indicates a nonzero
$\frac{d\langle\dot{m}_{0}\rangle}{dk}$.}
\label{f2}
\end{figure}

\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|}
\hline
\rule{0pt}{2.8ex}
Band & $\lambda$ (nm)& $\frac{d\dot m_0}{dz}$ ($10^{-2}\,$d$^{-1}$) & $\langle\dot m_0\rangle$ ($10^{-3}\,$d$^{-1}$)\\
\hline
\rule{0pt}{2.7ex}\!\!
 U   & $ 350.0 $ & $ -23\pm27$  & $ -22\pm3 $ \\%& $9.30\pm 4.87$ \\
 B   & $ 438.0 $ & $ -19\pm12$  & $ -17\pm3 $ \\%& $9.30\pm 4.87$ \\
 V   & $ 546.5 $ & $ -28\pm10$  & $ -26\pm3 $ \\%& $9.30\pm 4.87$ \\
 R   & $ 647.0 $ & $ -46\pm11$  & $ -34\pm4 $ \\%& $9.30\pm 4.87$ \\
 I   & $ 786.5 $ & $ -52\pm19$  & $ -39\pm7 $ \\%& $9.30\pm 4.87$ \\
\hline
\end{tabular}
\end{center}
\caption{Fitted means and $z$-dependences of $\dot{m}_{0}$ for the five
wavelength bands. $\lambda$ is the central wavelength for each band.}
\label{t5}
\end{table}

The last fit required was for the wavelength dependence of $\langle\dot m_0\rangle$.
The result was $2\pi\frac{d\langle\dot m_0\rangle}{dk}=(9.3\pm4.8)$ nm/d.
The correction associated with this factor will make only a small contribution to the
final value of $\frac{\dot\Gamma_0}{\Gamma_0}$. This correction, along with the other
factors necessary to produce the final computed value of
$\frac{\dot{\Gamma}_{0}}{\Gamma_{0}}$, are given in table~\ref{t1a}. The fourth
column of this table collects all the corrections for each band, giving the full
$z$-dependence of the $^{56}$Co decay rate $\Gamma_{0}$; only a factor
of the Hubble constant $H$ is needed to convert this into a time dependence.

\begin{comment}
{\bf One important question we have neglected is how the main correction applied
to the light curves to make them standard candles---a stretch factor $s$
---would affect our fittings. First, including $s$
should make it easier to determine $d_0$.
Second, the modification would change the plots shown in figure \ref{f2} and
therefore the inferred value of $\frac{\dot\Gamma_0}{\Gamma_0}$.
When $s$ is included, the slopes for individual sources seem to lie closer
to the average. This is particularly true of the outlier sources
sn1994M and 1995ak at high redshifts, since these sources have been
assigned small values of $s$ (see figure \ref{f3}). Including the effects of $s$ could improve
the apparent values on $\frac{\dot G_F}{G_F}$
by around one order of magnitude, but we feel that doing so at this
stage is unjustified. The physics behind $s$ is not very well understood,
and it is still unclear if $s$ should also affect the latter time of the
light curves---where the $^{56}$Co decay prevails.}
\end{comment}

\begin{figure}
\centering
\includegraphics[scale=0.65]{xgra1col}
\caption{V-band light curves of SN 1995ak and SN 1998aq, showing
linear fits in the $^{56}$Co decay region.
The maximum intensity occurs at $t=0$, where $m_\textrm V$ is normalized to
$-12.5$ and $-12$, respectively.}
\label{f3}
\end{figure}

\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
\rule{0pt}{3ex}
Band & $-\frac{1}{\langle\dot m_0\rangle}\frac{d\dot m_0}{dz}$ &
$\frac{k}{\langle\dot m_0\rangle}\frac{d\langle\dot m_0\rangle}{dk}$
& $\frac{1}{H}\frac{\dot\Gamma_0}{\Gamma_0}$ \\
\hline
\rule{0pt}{2.5ex}\!\!
 U   & $ -10.5\pm12.1 $ & $ -1.2\pm0.6 $ & $ -11.7\pm12.1 $ \\
 B   & $ -11.5\pm7.7  $ & $ -1.3\pm0.7 $ & $ -12.8\pm7.8 $ \\
 V   & $ -10.6\pm4.0  $ & $ -0.7\pm0.3 $ & $ -11.2\pm4.0 $ \\
 R   & $ -13.6\pm3.8  $ & $ -0.4\pm0.2 $ & $ -14.1\pm3.8 $ \\
 I   & $ -13.4\pm5.5  $ & $ -0.3\pm0.2 $ & $ -13.7\pm5.5 $  \\
 \hline
\end{tabular}
\end{center}
\caption{Contributions to $\frac{\dot\Gamma_0}{\Gamma_0}$ using supernovae in the domain
$3.9\times 10^{-3}<z<2.4\times 10^{-2}$.}
\label{t1a}
\end{table}

Using the values from table~\ref{t1a}, we may compute the weighted average of
$\frac{\dot{\Gamma}_{0}}{\Gamma_{0}}$ for all the frequency bands.
The result is $\frac{1}{H}\frac{\dot\Gamma_0}{\Gamma_0}=-12.7\pm7.4$, so
\begin{equation}\label{a11}
\frac{\dot{G}_{F}}{G_{F}}=(-4.8\pm2.8)\times 10^{-10}\,\textrm{yr}^{-1}.
\end{equation}
where we have used $H=(74.2\pm3.6)$ km/s/Mpc~\cite{ref-riess3}, the uncertainty in
which
introduces a small additional error.

\section{Conclusions}
\label{sec-concl}

Our main result, eq.~(\ref{a11}), shows that the Fermi constant $G_{F}$ may be
decreasing with time; however, the evidence for a change
is at less than the $2\sigma$ level. The error is too large to conclude that
there is any systematic change in $\Gamma_{0}$ or $G_{F}$. Stating our results
as a $2\sigma$ bound, we have
\begin{equation}
\left|\frac{\dot{G}_{F}}{G_{F}}\right|<10^{-9}\,\textrm{yr}^{-1}.
\end{equation}

The main source of
error was the fitting errors associated with the
determination of $\frac{d\dot m_0}{dz}$.
Nevertheless, the results for the variation of the decay rate of the $^{56}$Co
isotope seem to be roughly consistent across the five different wavelength bands.

Most of the statistical pull toward a
negative value of $\frac{\dot G_F}{G_F}$ comes from two sources, SN1995ak and SN1994M,
at comparatively high redshifts. The V-band light curves for SN 1995ak and the more typical
SN 1998aq are shown in figure \ref{f3}. The two sources at the highest $z$ values show markedly more
rapid light curve decays than do the lower-$z$ supernovae. When sources such as these are used for cosmic distance scale
measurements, the more rapid decay is removed by applying a stretch factor $s$
to the time axis~\cite{ref-stretch}. However, the use of such an empirically determined factor appears to be
inappropriate here; it could rescale away any actual signal from variations in $G_F$ that was
seen from the data. If $G_F$ is not varying, some other systematic difference between these
sources and the others is probably responsible, but the nature of this difference is unknown.

The second term in eq.~(\ref{a4c}), which comes from the $k$ dependence of $\dot m_0$,
is considerably smaller than the first term, which accounts for its direct dependence on
$z$. Thus, the poorly understood physics that causes the light curve to depend on
wavelength contributes only a small adjustment to the inferred value of
$\frac{\dot{G}_{F}}{G_{F}}$. This was found to be true in all the bands.
Since the leading order correction due to the $k$-dependence is small, as is $z$ for the
sources we have considered, we do not expect there would be any important changes to
$\frac{\dot\Gamma_0}{\Gamma_0}$ if higher order in $z$ corrections were taken into account.
However, if sources with greater redshifts were to be analyzed, the higher order
corrections would need to be considered and their contributions could be important.

More sources should be analyzed for this method to produce its best possible result
for the time variation of the Fermi constant.
Although the CfA Supernova Archive contained plenty of information for small redshift
supernovae ($z<2.4\times 10^{-2}$), more large redshift
supernovae should also be analyzed in order to obtain better results. At present,
the $10^{-9}$ yr$^{-1}$ bound on the fractional time variation can be seen as
ruling out roughly a 10\% change in $G_{F}$ over about 100 million years of time.
Although the late
time behavior of light curves is more difficult to observe for more distant
supernovae, using older sources could significantly improve these results.
It would also be desirable to have more data with well-defined locations for their
maxima; light curves with more
data points in the linear region of interest would also be particularly valuable.

We found no significant evidence for a change in $G_{F}$ over the cosmological
interval we studied, $3.9\times 10^{-3}< z< 2.4\times10^{-2}$. However, if we
did see a systematic shift in the decay rate of the supernova light
curves, it could probably be more reasonably explained by something other
than a change in the strength of the weak interaction. We have not
measured $G_{F}$ (or the half-life of $^{56}$Co) directly. Changes in supernova
composition over the lifetime of the sample could potentially have
affected the decay rate of the light curve, changing the way the $^{56}$Co
decay products heat the expanding remnant. Statistically significant evidence for
systematic changes in the supernova light curves over the lifetime of the universe would
be suggestive, but it would by no means be compelling evidence for changing fundamental
constants.

This analysis has determined the rate of change of the dimensional quantity
$G_{F}$. In general, only the time dependence of dimensionless numbers can
be defined unambiguously. A quantity with dimensions must be compared with
some similarly dimensioned reference unit, which may itself depend on time
in a correlated fashion; what reference is chosen will affect how the
measured quantity appears to depend on time. However, there is
essentially no problem in this instance. We can calculate how the
dimensional parameter $\Gamma$ (or $G_{F}$) depends on time with no
meaningful ambiguity---just as one can observe the growth in a child's weight over
time without confronting any ambiguity---since any reference
with which the dimensional observable is implicitly compared is
independent of $G_{F}$. All possibly relevant time standards with which
$\Gamma_{0}$ could be compared have the same dependence on $G_{F}$---none.

In order to extend these bounds to data sets with older, more distant supernovae, we
would need to account for several additional complications. For very distant sources
(and especially for late periods along the light curve), the quality of the magnitude
data diminishes. In addition to limiting the amount of usable data, this could
introduce unexpected biases into the data. (For example, light curves with more rapid
$^{56}$Co decay rates might be underrepresented in the set of usable data, because their
luminosity falls below the threshold for detection more quickly, providing too few
data points to map out the relevant region of their light curves accurately.)
A more systematic $K$ correction would probably also be required at high redshifts,
and it would be desirable to make use of more detailed models of type Ia supernovae
structure, so as to get at the $^{56}$Co lifetime more directly.

These are not the first published bounds on $\dot{G}_{F}$. Although the Fermi coupling
constant cannot straightforwardly be observed through its effects on atomic spectra, it
does effect nuclear energy levels in a small way. Aside from astrophysical observations
of ancient events and precision laboratory experiments, the most successful method
for placing bounds on the time variations of fundamental constants has been to look at
the isotopes produced by the natural reactors that operated near Oklo, Gabon about
$2\times 10^{9}$ years ago. A bound on the fraction change in the dimensionless quantity
$\beta=G_F m_p^2$ at the $10^{-11}$ yr$^{-1}$ level was based on
the production of samarium isotopes at Oklo~\cite{ref-damour}. However, bounds based on
the Oklo data require a very detailed understanding of the operation of the reactor, as
well as complicated nuclear structure calculations, whose reliability as regards the
effects of $G_{F}$ is somewhat questionable. Therefore, our technique, which can be
extended to older supernovae, may offer the best possibility for placing robust bounds
on the time variation of the Fermi constant.

\section*{Acknowledgments}
The authors wish to thank Santiago Gonzalez for his helpful comments.
This research has made use of the CfA Supernova Archive, which is funded in part by the
National Science Foundation through Grant No. AST 0606772.