text:
\section{Introduction}
\label{sect:intro}
The Infrared Array Camera (IRAC) \cite{2004ApJS..154...10F} was operational on the {\it Spitzer} Space Telescope \cite{2004ApJS..154....1W} from 2003 - 2020 with four broad mid-infrared bands with response covering 3.15 - 9.25 $\mu$m (at 3.6, 4.5, 5.8, and  8.0 $\mu$m; also denoted channels 1 -- 4, respectively).  These filter names are labels and not the actual effective wavelengths (for more detailed information on filter transmission see Ref.~\citenum{2008PASP..120.1233H}).  After depletion of the cryogens, from mid 2009 until January 2020, only [3.6] and [4.5] were available for observations.

In order to efficiently calibrate the telescope, a discrete set of observing templates was enabled for IRAC observers.  Available options to define the templates included array mode, frame time, and dithering strategy.  The array mode could be specified either as full array (meaning, the entire IRAC 256$\times$256 array of 1.2\,arcsecond pixels was read out after each exposure) or subarray, in which case 64 consecutive 32$\times$32 pixel images were taken at a higher readout rate without moving the telescope (0.01~s  single readout for subarray vs. 0.2~s for full array).  A number of preset frame times were available; those and their corresponding exposure times are listed in Table~\ref{tab:frametimes}.  Frame time is a measure of time elapsed between successive array resets.  The exposure time (effective integration times)  is the time elapsed between the first pedestal and the first signal read, not including resets or multiple reads.  High dynamic-range exposures, in which paired short and long exposures were acquired sequentially within a single full-array observing template, were also available to IRAC observers, but did not add diversity to the range of possible observations.  Lastly, observers were able to choose from a discrete set of dither patterns, make their own mapping strategy, or employ a staring-mode.  Dithering (small position changes to reduce noise) and mapping (position changes of order the size of the array to cover an area greater than the size of the array) were both possible in the full array mode.  Dithering was available in the subarray mode between 64-frame sets, but not within the sets.  Observations which do not move positions between frames are known as ``staring-mode.''  Staring-mode was most often used for high-precision time-series observations of brown dwarfs and exoplanets.

\begin{table}[ht]
\caption{Available IRAC Frame Times}
\label{tab:frametimes}
\begin{center}
\begin{tabular}{|l|l| l|}
\hline
\rule[-1ex]{0pt}{3.5ex}  Frame Time & Exposure Time & Observing Mode \\
\rule[-1ex]{0pt}{3.5ex}  (seconds) & (seconds) &  \\

\hline\hline
\rule[-1ex]{0pt}{3.5ex}  0.02 & 0.01  & subarray\\
\hline
\rule[-1ex]{0pt}{3.5ex}  0.1  & 0.08 & subarray  \\
\hline
\rule[-1ex]{0pt}{3.5ex}  0.4 & 0.36 & subarray   \\
\hline
\rule[-1ex]{0pt}{3.5ex}  2 & 1.92 & subarray  \\
\hline \hline
\rule[-1ex]{0pt}{3.5ex}  0.4 & 0.2 & full array  \\
\hline
\rule[-1ex]{0pt}{3.5ex}  2 & 1.2 & full array \\
\hline
\rule[-1ex]{0pt}{3.5ex}  6 & 4.4 & full array\\
\hline
\rule[-1ex]{0pt}{3.5ex}  12 & 10.4 & full array \\
\hline
\rule[-1ex]{0pt}{3.5ex}  30 & (23.6, 26.8) \tablefootnote{ Entries with two exposure times are for [3.6] and [4.5] respectively.} & full array \\
\hline
\rule[-1ex]{0pt}{3.5ex}  100 & (93.6, 96.8) & full array \\
\hline

\end{tabular}

\end{center}
\end{table}

The rich {\sl Spitzer}/IRAC archive contains many observations of the same target taken with different array modes, frame times, or dithering strategies.  This is often the case for serendipitous observations but also for observations designed intentionally in different modes.  Possible science cases for this include, but are not limited to, 1) initial dithered observations to find proper motions of brown dwarfs followed by staring-mode data to refine the characterization of their atmospheres, 2) dithered debris disk searches looking for IR excess followed by a staring-mode study of variations in specific debris disks, 3) archival dithered observations followed by targeted subarray staring observations for confirmation, or 4) observations that had either a too low signal-to- noise ratio (SNR) or were saturated at one exposure time were followed up later in the mission using a different exposure time.  These science cases therefore require combining photometry from different modes to arrive at scientific conclusions.

In this work we examine how IRAC photometry depends on observing strategy.  We emphasize that these are very small measured variations (a few percent at most) and so will likely only affect science users who require high-precision photometry from multiple different observing modes.  A full set of calibration observations to test for variations in measured fluxes among observing modes were only taken during the warm IRAC mission. We consequently do not discuss data taken during the cryogenic mission (which ended in 2009), or in the [5.8] or [8.0] channels.

In Section \ref{sec:data} we describe the archival data used for this project.  Section~\ref{sec:photometry} covers our methods for reducing the data, carrying out photometry, and applying photometric corrections.  Section~\ref{sec:results} discusses the different potential correlations and photometry effects.  We make concluding remarks in Section~\ref{sec:conclusions}.

\section{The Data}
\label{sec:data}

\subsection{Observations}
\label{sec:observations}
For this work we used the calibration observations taken for {\sl Spitzer} Program ID (PID) 1336 and PID 1367.  Specifically, we observed a set of ten stars with varying exposure times, in both subarray and in full array, and in staring and dithering modes.  Almost always, staring-mode observations were taken on the same pixel ("the sweet spot"), as that pixel was the best characterized pixel on the array.  Table~\ref{tab:stars} lists the stars.  We searched the archive for additional data usable for this analysis but did not find anything suitable, combining non-variable stars with observations in all available modes and having a sufficient number of images to achieve statistical significance.

Figure~\ref{fig:flux_DN} shows a visualization of this dataset.  We plot the exposure time(s) vs.\ aperture flux for just the [4.5] channel.  Similar observations were made for [3.6]. All stars were observed with multiple exposure times.  Not all of the ten stars could be observed at all exposure times due to SNR and saturation concerns.  This plot illustrates how the range of stellar brightnesses and frame times in the
sample filled the available phase space.

We considered how close to saturation(well depth) a star is  as a means of interpreting our results.  Looking at our sample as a whole, the range of possible well depths is not well sampled;  having a median fractional well depth of 0.04 and 0.02 at [3.6] and [4.5], where a fractional well depth of 1.0 indicates saturation. The right panel of Figure~\ref{fig:flux_DN} shows well depths of the sample. The median well depths of our sample correspond to a SNR of 39 and 24 at [3.6] and [4.5], respectively. This dataset was not designed with well depth in mind, and was instead designed to find stars which would sample the available exposure time parameter space.

Although available, we reject the 0.02~s photometry because of its large scatter. We have also rejected any data where the well depth is greater than the listed saturation limit in the IRAC Instrument Handbook.\footnote{\label{handbook}\url{https://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/iracinstrumenthandbook/}}  The saturation limit results in the rejection of around 250 photometry points at [4.5].   All observations used in this work have SNRs of six or greater.

For a subset of three stars, we made a larger number of observations in all possible modes  with 2~s frame times in both channels.  This allowed us to build up statistically meaningful samples where we hold the frame time constant while varying the other observing parameters. The three stars are NPM1+66.0584, NPM1+66.0578, and KF03T2.

The 2~s frame time is the closest set of IRAC frame times available for holding exposure time constant between the full and subarray. Unfortunately, the 2~s subarray and 2~s full-array frame times do not actually correspond to the same exposure times; 1.92~s and 1.2~s, respectively.  Similarly, frame times of 0.4~s in full and subarray also have different exposure times, so that frame time was not a good alternative.   This difference in the exposure times for the same frame time comes from a different number of Fowler numbers (N), wait ticks, and readout times for the different modes as listed in the IRAC Instrument Handbook. Briefly, Fowler sampling is a way of observing by taking N non-destructive reads at the beginning of the observation and another N non-destructive reads at the end of the observation, the difference of which is the flux measured per pixel.

All of the stars chosen for this work were vetted as potential calibration stars for IRAC.  Seven of the ten stars are published as primary or secondary calibrators for IRAC known to not exhibit flux variability in the IRAC bands \cite{2005PASP..117..978R}. The three stars not included in that reference are NPM1+74.0514, NPM1+57.0835, and NPM1+66.0584.  Specifically NPM1+66.0584 is one of the three stars in the subset taken with more observations in 2~s frame times.  We cannot use this dataset to determine both if the stars are time variable, and if they vary as a function of the other parameters studied herein.  While vetted,  because these have not specifically been published as calibration stars, we experimented with removing these three stars from our sample.  All plots look similar (albeit with larger scatter due to fewer data points), and conclusions remain the same if we remove those stars from the sample.  We therefore choose to keep these in the sample for the remainder of this work.

  Each of the [3.6] and [4.5] datasets include about 80,000 total individual observations. Specifically, for each channel we have roughly 65,000 subarray and 15,000 full array observations.

\begin{figure}
\begin{center}
\begin{tabular}{c}
\includegraphics[height=9cm]{flux_welldepth_ch2.png}
\end{tabular}
\end{center}
\caption
{ \label{fig:flux_DN}
Description of the dataset used in this work for the [4.5] channel, showing the aperture flux as a function of the exposure time on the left panel and aperture flux as a function of well depth on the right panel.  The [3.6] channel has a similar observation set but for brevity is not shown. Points are color coded by star names, and symbol shapes denote full (filled circles) or subarray (crosses).  Each data point represents one observation, and data are not binned, but there are overlapping data points. }
\end{figure}

\begin{table}[ht]
\caption{The Calibration Star Sample}
\label{tab:stars}
\begin{center}
\begin{tabular}{|l l l l l l l|}
\hline
\rule[-1ex]{0pt}{3.5ex}  Name & R.A. & Dec. & Spec.  & [3.6] & [4.5] & Alternate \\
\rule[-1ex]{0pt}{3.5ex}  &  (h m s) &  (deg' '') & Type & (mJy) & (mJy) & Name\\

\hline\hline
\rule[-1ex]{0pt}{3.5ex} NPM1+61.0569 & 17 23 25.9 & +61 12 40.7 &K0.5 III & 363	& 208 & \\
\hline
\rule[-1ex]{0pt}{3.5ex} NPM1+74.0514 &	19 02 53.5& 	+74 14 43.6	&K0.5 III&	323&	194&	BD+74 804  \\
\hline
\rule[-1ex]{0pt}{3.5ex}  NPM1+61.0582 &	17 36 55.6&	+61 40 58.1	&K1 III &	274&	162 & \\
\hline
\rule[-1ex]{0pt}{3.5ex}  NPM1+66.0507 &	17 31 22.1&	+66 46 35.3&	K2 III&	250&	145& \\
\hline
\rule[-1ex]{0pt}{3.5ex}  NPM1+66.0578 &	19 25 32.2 &	+66 47 38&	K1 III &	127&	75 & \\
\hline
\rule[-1ex]{0pt}{3.5ex}  NPM1+57.0835	&18 04 12.6	&+57 42 18.6&	&	95&	54& \\
\hline
\rule[-1ex]{0pt}{3.5ex}  TYC4213-241-1&	18 03 45.53	&+66 56 03.7&	K1 III &	74	&42& KF01T3\\
\hline
\rule[-1ex]{0pt}{3.5ex}  KF03T2 &	17 57 51.4&	+66 31 03.0 &	K1.5 III&	70&	40 &\\
\hline
\rule[-1ex]{0pt}{3.5ex}  NPM1+66.0584&	19 36 07.39&	+66 21 54.2&	F0	&56	&35	&BD+66 1222 \\
\hline
\rule[-1ex]{0pt}{3.5ex}  NPM1+60.0581 &	17 24 52.3&	+60 25 50.8	&A1 V &	37&	23&	BD+60 1753 \\
\hline

\end{tabular}

\end{center}
\end{table}

\subsection{Photometry}
\label{sec:photometry}

We briefly describe our pipeline here, emphasizing where it is differs from previous work.  For an overview of the {\sl Spitzer}/IRAC absolute photometric calibration, see Refs.~\citenum{2005PASP..117..978R, 2012SPIE.8442E..1ZC}.  Specifics of the photometry pipeline used here are included in Ref.~\citenum{2021AJ....161..177K}. We used the basic calibrated data (BCD) FITS files (suffix bcd.fits) from the Spitzer Heritage Archive (SHA).  These BCD files have had a dark correction, flat-field correction, linearity correction, and conversion to flux units already applied.

One somewhat novel aspect of this work is that we apply a different dark correction to staring-mode data than we do to dithered observations.  Staring-mode data are processed in the standard pipeline in the same manner as dithered data, including using a dark image which was made by dithering.  The delay time between frames affects the bias level in the frame (including in the dark); this is known as the ``first-frame effect.''  Therefore, dithered observations will have different bias levels and patterns than staring-mode observations since staring-mode observations have shorter delay times between consecutive frames, although that pattern is constant if the delay time is constant.  This effect adds both noise and systematics to the photometry.  For this reason, dithered darks are inappropriate for staring-mode science frames when precision photometry is required.

We used PID 1345 to make our own staring-mode dark suite for all subarray frame times. We began by removing the dithered dark correction from the PID 1345 data.  Because the dark correction is not the last correction made to the BCD files, care was taken to first back out the other corrections, apply the staring-mode dark, then  re-apply the other pipeline corrections. We then created a median image for each exposure time (0.02, 0.1, 0.4, 2.0 s) and used this median frame as the staring-mode dark.  Applying this new staring-mode dark to the staring-mode data has a measurable impact on derived fluxes. We recommend that anyone doing precision photometry with staring-mode data use a staring-mode dark instead of the pipeline-provided, dithered dark.  While the IRAC pipeline will not include these starting-mode darks, code is available for users to to change which darks are used in a BCD frame on the contributed-code section of the {\it Spitzer} IRAC website. \footnote{\url{https://irsa.ipac.caltech.edu/data/SPITZER/docs/dataanalysistools/tools/contributed/irac/change_dark_calibrate/}}

To measure flux, we use our appropriately dark-corrected BCD exposures and make the following corrections in order.  We first convert images into units of electrons to enable a statistical calculation of uncertainties.  Second, we use a center-of-light method to find stellar centroids. \footnote{\url{https://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/calibrationfiles/pixelphase/box_centroider.pro}}  Third, we do aperture photometry with a three-pixel radius aperture and (3-7) pixel background annulus.  The small aperture size is chosen to reduce noise and the number of cosmic rays in the aperture.  Fourth, we make a correction for pixel-phase using pixel\_phase\_correct\_gauss.pro. \footnote{\url{https://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/warmfeatures/pixel_phase_correct_gauss.pro}}  The pixel-phase correction accounts for gain changes  as  a  function  of the position  within  a pixel, coupled with the undersampling of a point source by IRAC.  Fifth, we make a correction for array location. The array location-dependent correction takes into account the variation in system response of the instrument across the field of view, which is primarily  due  to  the  change  in  the angle  of  incidence  of light through the bandpass filter as a function of position on the array. Lastly, we discard the first frame of every subarray FITS file and of every full array AOR.  These frames are affected by the first-frame effect discussed above and are likely to have measured fluxes that differ from those of subsequent images. We do not apply an aperture correction since the same aperture is used for all photometry regardless of observing mode.

To compare photometry for all stars on the same plots, we normalize the stars to the same absolute level by dividing all photometry by the median stellar fluxes.  The distributions of fluxes per star are somewhat skewed, so a mean does not capture the peak of the distribution. Having skewed distributions causes the mean levels in the subsequent plots to differ from unity.

\section{Results and Discussion}
\label{sec:results}
This section describes how the measured IRAC photometry differs as a function of observing mode and exposure time. Sec.~\ref{sec:disentangle} examines the effects of observing mode and exposure time together. Finally, we consider the impact of staring and dithering on IRAC photometry.

We choose to present our results mainly with box plots to show the distributions per star as a function of both array mode (full or subarray) and exposure time for both channels, discussed individually below in Sections~\ref{sec:readmode}~and~\ref{sec:exposure time}.  A box plot shows the median of the distribution as the solid line in the middle of the box. The box top and bottom indicate the top and bottom quartiles and the caps at the end of the lines show the maximum and minimum values in the distributions. These box plots include the entire dataset (all stars, all exposure times, full and subarray, staring and dithering).  On average, each box contains a few thousand data points. Generally the boxes with smaller quartile ranges and max/min values have fewer data points. The plots are color-coded by star for each of the ten stars.  Stars are listed in order from the brightest (NPM1+61.0569) to the faintest (NPM1+60.0581).  Color-coding is consistent for all box plots and labeled in Figure~\ref{fig:readmod_ch1}.

\subsection{Array Mode - Full vs. Subarray }
\label{sec:readmode}

Figures~\ref{fig:readmod_ch1} shows the [3.6] and [4.5] photometry (top and bottom respectively) of the stars in Table~\ref{tab:stars} as a function of array mode, revealing that the measured fluxes appear on average to be lower for the full-array observations than the subarray observations of the same stars. We attempt to confirm this statistically by using an Anderson two-sample test, per star, to see if the distributions for full array and subarray are drawn from the same population. This statistical method considers the vertical distance between the two cumulative distributions. For all ten stars we can say that the full array and subarray data are not drawn from the same distribution at the 25\% significance level (maximum possible significance).  Thus, the difference between the fluxes measured in full array and subarray modes are statistically significant. Using all ten stars, the median difference between the full and subarray flux measurements is 1.6\% at [3.6] and 1.0\% at [4.5].

\begin{figure}
\begin{center}
\includegraphics[width = 0.75\linewidth]{Figure2.pdf}
\end{center}
\caption
{ \label{fig:readmod_ch1}
Box plot of the normalized fluxes of all the stars divided into full and subarray observations for [3.6] on the top and [4.5] on the bottom. Stars are normalized by dividing by the median flux.  The solid line in the middle of the boxes shows a median value, while the top and bottom of the box are quartiles, and the endcaps indicate the full range of the data.  Stars are color coded according to the legend and listed in order of brightness from brightest to faintest. On average, the measured fluxes in full-array observations are smaller than the measured fluxes in the subarray observations. }
\end{figure}

\subsection{Frame Time}
\label{sec:exposure time}

Figure~\ref{fig:exptime_ch1} shows the distributions of normalized calibrator star fluxes as a function of Frame time.  [4.5] includes more data on the six-second full-array frame time, so the bottom plot includes that frame time whereas the [3.6] plot does not. The median normalized flux of all the stars changes as a function of frame time such that the measured fluxes are larger at lower frame times.  This is most strongly evident in [3.6].  Figure \ref{fig:frametime_effect} shows a measure of the strength of this effect as the median difference between the 2~s and 0.4~s frame time fluxes for all the stars, per channel, per observing mode.

As in Section~\ref{sec:readmode}, we used the Anderson two-sample test to compare the full array 0.4~s observations with the full array 2~s observations.  For both channels, we find the difference between the normalized flux distributions at the two frame times to be statistically significant.  The same is true for fluxes from the subarray in 0.4~s and 2~s frame time observations.

We tested whether the trend seen in flux as a function of exposure time depends on stellar brightness such that the brighter or the fainter stars would be more or less likely to show this effect.  We do this by calculating the difference between the 2~s and the 0.4~s normalized fluxes.  A difference between the fluxes implies that different fluxes are measured using different exposure times.  Because we know there is an effect with array mode, we divide the sample into full and subarray data.  Finally, looking for a trend in the stellar brightness, we divided the sample into bright and faint stars, and re-calculated the median values.  Figure~\ref{fig:frametime_effect} shows the results, where the circular points are the median values of the difference in fluxes between the 2~s and 0.4~s frame-times for both the full array only data (blue) and the subarray only data (red). The bright and faint subsets are shown with the appropriately colored triangles. The differences between the bright and faint samples are within one sigma of each other. Here sigma is calculated as the standard deviation between stars, which means that each bin contains only a few stars (and not thousands of data points).  Our sample size is not large enough to make conclusions based on star brightness.

\begin{figure}
\begin{center}
\includegraphics[width = 0.75\linewidth]{Figure3.pdf}
\end{center}
\caption
{ \label{fig:exptime_ch1}
Box plot of the distributions of normalized fluxes of all the stars as a function of exposure time at [3.6] on the top and [4.5] on the bottom. Frame time labels also indicate full or subarray.  Color coding is the same as for Figure~\ref{fig:readmod_ch1}.  On average, fluxes measured at longer frame times are lower than those measured at shorter frame times.}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{c}
\includegraphics[height=8cm]{frametime_effect.png}
\end{tabular}
\end{center}
\caption
{ \label{fig:frametime_effect}
Exposure-time effect for both channels divided into bright and faint stars to look for an effect with star brightness.  [4.5] is shown in the top half of the diagram, while [3.6] plotted in the lower half.  Blue circles show the median difference for all the stars observed in the full array between the 2~s and 0.4~s frame time normalized fluxes, including one sigma error bars.  Blue right and left triangles show the median values for the bright and faint populations of stars.  The same quantities are plotted for subarray data in red.}
\end{figure}

\subsection{Disentangling Exposure Time and Array Mode Effects}
\label{sec:disentangle}
We examine the possibility of disentangling the effects of exposure time and array mode. Instead of dividing the dataset by star, here we consider all stellar photometry as a single dataset, and divide the dataset into four categories depending on the combination of array mode (full array vs.\ subarray) and exposure time (short vs.\ long). Each distribution has between 3000 and 50000 stars (subarray exposure times have lots more frames than full array). Exact exposure times cannot be compared in this analysis because they differ significantly for full and subarray modes (see Table~\ref{tab:frametimes}).  For that reason, the division between the short and long exposure times is set at 0.3~s, to construct statistically significant samples.

\begin{figure}
\begin{center}
\includegraphics[width = 0.75\linewidth]{Figure5.pdf}
\end{center}
\caption
{ \label{fig:corrflux_dist_ch1}
Distribution of normalized star fluxes for [3.6] on the top and [4.5] on the bottom.  The legend lists the different observing modes (full array vs subarray and long exposure times vs. short exposure times) as well as the mean well depth per distribution.  Solid dark lines show the kernel density estimation (KDE) overplotted on the histograms. We see a difference in the absolute photometry of stars taken in  available observing modes and exposure times shown here in the colored distributions.  }
\end{figure}

Figure~\ref{fig:corrflux_dist_ch1} reveals differences between the overall distributions of photometry taken in full and subarray mode while accounting for exposure time, with enough data points to accurately reflect the shapes of the distributions.  Longer exposure times are in red and light blue; shorter exposure times are in orange and dark blue.  Especially at [3.6], an effect is apparent with both exposure time and observing mode. For the stars which have full-array observations in both short and long exposure times, the median difference in the flux is  3.4\% at [3.6] (for six stars) and  0.6\% at [4.5] (for ten stars). For subarray we find the median difference in flux between short and long exposure times to be 2.9\% at [3.6] and 1.5\% at [4.5].

\subsection{Well Depth}
\label{sec:well_depth}
The legend to Figure~\ref{fig:corrflux_dist_ch1}  lists median well depths for each of the distributions.  Well depth is one physical feature which correlates with the difference between short and long and full and subarray exposures.  We do expect that longer exposure times on the same set of stars will have larger fractional well depths, so it makes sense that the full long and sub long distributions have the larger median well depths in both channels.  Also, subarray mode has the possibility of shorter frame times than full array, so we would expect that subarray would have lower well depths than full array.  Both channels show this behavior.

Figure~\ref{fig:welldepth_normflux} shows the distributions of normalized fluxes for the entire sample of stars in [3.6] and [4.5] divided into three well-depth bins.   These bins do not have the same number of photometry points in them as we have  many more low well-depth observations than intermediate well-depth observations.  Because of this change in the number of data points per bin, we do not quote values for the peak of each bin, but rather address the trends revealed by the data in hand [3.6] shows a clear trend such that observations at intermediate well depths have lower normalized fluxes.  [4.5] shows little, if any trend, consistent with trends seen in the observing mode and exposure time plots.  While there are sources at higher well depth than 50\% full well, their distribution overlaps those of the 20-50\% bin. We conclude that only fractional well-depths less than a few tens of percent are affected by this well-depth effect.  Consideration of well-depth does not help to disentangle the effects we have seen with observing mode and exposure time, but it is a clue that potentially the source of some of what we are seeing is non-linearity in the low well-depth regime.

\begin{figure}
\begin{center}
\includegraphics[width = 0.75\linewidth]{well_depth_normflux.png}

\end{center}
\caption
{ \label{fig:welldepth_normflux}
Distributions of fluxes in three fractional well depth bins for [3.6] on the left and [4.5] on the right. Colors for the distributions get darker as the fractional well depth increases.}
\end{figure}

\subsection{Dithering vs. Staring-mode}
\label{sec:stare_dither}
 We do not detect a difference in photometry between the dithered and staring-mode observations.  For three stars we have observations in all the modes for the 2~s exposure times (sub staring, sub dithering, full staring, full dithering). Recall that subarray observations taken with 2~s frame-time have 1.92~s exposure times and 2~s frame time full-array observations have 1.2~s effective exposure time, so the 2~s exposure time subarray and full array observations are not equivalent, but they are the closest to equivalent that exist in the archive.

Figure~\ref{fig:cleveland_dot} shows a Cleveland dot plot of the median and standard deviation of the normalized fluxes for each mode, for each of the three stars in this sample. The legend also lists the number of data points per mode.  The full and subarray measurements differ as before( see Section~\ref{sec:readmode}). The staring-mode photometry is consistent with the dithered mode photometry within one sigma in both channels.

\begin{figure}
\begin{center}
\begin{tabular}{c}
\includegraphics[height=8cm]{sub_full_stare_dither_noff_staredark.png}
\end{tabular}
\end{center}
\caption
{ \label{fig:cleveland_dot}
Cleveland dot plot comparing normalized fluxes from different observing modes and dithering vs. staring for three stars listed on the y-axis.  [3.6] is on the left, [4.5] is on the right.  The legend is for both plots, and lists the number of data points for each mode.  }
\end{figure}

\subsection{Multiple Regression}
Finding a relation among the observing parameters (array mode, exposure time, staring/dithering) may help both to a) understand the effects found in this paper and b) correct for them. We therefore use the statistical technique of multiple regression to search for any relations.   Specifically, we used Multiple Linear Regression, i.e., multiple independent variables, to predict the value of the dependent variable.  To correct IRAC photometry for these effects, we have tried ordinary least squares (OLS) as a multiple linear-regression technique.  We used array mode, exposure time, and stare/dither as independent variables, and flux as the dependent variable.  We used two different modules in Python for this work, statsmodels\cite{seabold2010statsmodels} and sklearn\cite{scikit-learn}.  We are unable to find successful models.  The $R^2$ goodness of fit is 0.037, when "good" models should have values close to 1.0.  The failure of OLS in this situation could imply that the relationship between the independent and dependent variables is nonlinear.

\section{Conclusions}
\label{sec:conclusions}
We document exposure time, array mode (full vs. subarray), and fractional well depth flux dependence in IRAC photometry. The full array has median fluxes higher than the subarray  by 1.6\% at [3.6] and 1.0\% at [4.5].  Dividing the sample further, the long exposures have a lower median flux than shorter exposures.  The difference at full array is 3.4\% at [3.6] and 0.6\% at [4.5]. For subarray we find the median difference in flux is 2.9\% at [3.6] and 1.5\% at [4.5]. These noted effects are only relevant for a small fraction of IRAC high precision users with data in multiple modes who should include these values in their uncertainty calculations. We posit two potential causes of the noted low-level differences in photometry.

Overall, the normalized fluxes decrease as well depth increases up to a few tens of percent full well.  Well depth is correlated with the differences in the flux distributions between full and subarray and long and short exposures.  While these correlations with array mode and exposure time are expected, it potentially indicates the presence of non-linearities at the low well depths sampled in this work.

Linearity corrections were made for warm IRAC to correct for a known effect where an increase in incoming photons does not correspond to an increase in counts (or data numbers DN).  This occurs because filling the well decreases the potential, which in turn results in a less responsive system.  The linearity correction for the IRAC InSb arrays inherently assumes no linearity correction at low well depths. Therefore, time-intensive observations were not made during the mission to include extremely low well-depth observations. Instead, calibration observations focused on the moderate ($>$ 20\%) to high well depths where the linearity deviated most significantly from a straight line correlation between photons and DN (see IRAC data handbook for a description of the derivation of the warm linearity correction. \footnote{\url{https://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/iracinstrumenthandbook/}})  One explanation for why this low count linearity effect could exist at [3.6], but not [4.5] is that the applied biases are different between [3.6] and [4.5] implying that the electric potentials are different, which could explain the stronger effect at [3.6] than [4.5].

Understanding the root cause of non-linearities is beyond the scope of this paper.  Many complicating issues are hiding under that designation including the 3D structure of a pixel, how the linearity interacts with Fowler sampling (specifically how to translate corrections derived at one set of Fowler sampling parameters to those in another, which would get worse as the integration times become comparable to the time spent reading the detector), the speed at which the electric fields are changing in these very short exposures compared to the timescales of the exposures or the readouts and persistent image trap filling.

Linearity can in some detectors depend on the flux of the source in the sense that it matters not only how many photons come in to the detector, but also the rate at which they fill the well.  We do not have enough observations to know if the effect we are seeing is dependent on the brightness of the stars.

A second possible explanation for the difference between full and subarray photometry is that the array is resetting faster between consecutive frames for the subarray.  A faster reset applies a stronger reverse bias on the array more frequently, which could affect the distribution of photoelectron traps, and therefore could have a low-level effect on photometry. We know that other similar effects ( e.g., "first-frame effect") are more significant at [3.6] than [4.5].\footnote{\url{https://irsa.ipac.caltech.edu/data/SPITZER/docs/files/spitzer/som12.2.pdf}}

No flux difference is apparent between the staring and dithering mode observations after correcting the staring-mode data by a staring-mode dark.  Anyone doing absolute photometry with staring-mode data should be using a staring-mode dark.  The full and subarray staring-mode datapoints considered here are all taken on the same array pixel (the sweet spot of the subarray), which is not true for the dithered positions.  Thus, residual pixel-phase uncertainties cannot be the cause of the measured flux differences between the full and subarray staring-modes, otherwise we would see this effect in this work when comparing staring and dithering modes.

We have no evidence for low level persistent images being the cause of the differences in photometry measured here.  We know that the persistent images are stronger at [3.6].  However, for [3.6] the same fraction of observations were taken at the sweet-spot pixel in both full and subarray images.  The sweet spot is the pixel at the center of the array where most observations are conducted in the subarray because it has been the best characterized for the pixel phase effect.  This means that the sweet spot pixel is more likely to have frequent low-level persistent images than other pixels on the array.  If a larger fraction of observations in the subarray had been taken at the sweet spot than the full array, we might have expected persistent images to be the culprit.  On the contrary, we see no evidence for this.

Our analysis implies that these differences in fluxes are systematic; they do not average out with more observations.  All full-array photometry will be different than all subarray photometry, no matter how many observations are taken in any given observing mode.

\subsection* {Acknowledgments}
We thank the anonymous referees for their time and care in providing very useful comments on this manuscript.
This work is based in part on observations made with the {\it Spitzer} Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This research has made use of NASA's Astrophysics Data System, the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration, and the SIMBAD database, operated at CDS, Strasbourg, France.  The acknowledgements were compiled using the Astronomy Acknowledgement Generator.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
\label{Intro}

The search for space-time dependence of fundamental constants plays a fundamental role in the continuous efforts to put in firmer empirical grounds our current physical theories and, at the same time, explore the possibilities of exotic physics that might become manifest trough small deviations. The experimental research can be grouped into astronomical and local methods. The latter ones include geophysical methods such as the natural nuclear reactor that operated about $1.8\times 10^9$ years ago in Oklo, Gabon \citep{DD96,Petrov06,Gould06}, the analysis of natural long-lived $\beta$ decays in geological minerals and meteorites \citep{Olive04b} and laboratory measurements such as comparisons of rates between clocks with different atomic numbers \citep{PTM95,Sortais01,Marion03,Bize03,Fischer04,Peik04}. The astronomical methods are based mainly on the analysis of high-redshift quasar absorption systems. Most of the reported data are, as expected, consistent with null variation of fundamental constants. Nevertheless, there are reports of intriguing results. For instance \citet{Webb99} and \citet{Murphy03b} have reported observations made with the Keck telescope which suggest a smaller value of the fine structure constant ($\alpha$) at high redshift as compared with its local value. However, an independent analysis performed with VLT/UVES data gave null results \citep{Srianand04}. Furthermore, a recent analysis using VLT/UVES data suggests also a variation in $\alpha$ but in the opposite sense, that is, $\alpha$ appears to be larger in the past \citep{Murphy10}. The discrepancy between Keck/HIRES and VLT/UVES is yet to be resolved. In particular, the two studies rely on data from different telescopes observing different hemispheres and it was pointed out that the Keck/Hires and VLT/UVES observations can be made consistent in the case where the fine structure constant is spatially varying \citep{Murphy10}.

Focusing on a different quantity, observations of molecular hydrogen in quasar absorption systems can be used to set constraints on the electron-to-proton mass ratio $\mu \equiv m_e/m_p$ at high redshift \citep{King08, Thompson09,Malec10}, while the present value of $\mu$ can be constrained using comparisons of different transitions in atomic clocks \cite{Bize03,Fischer04,Peik04}. Furthermore, the observed temperature isotropy of microwave background radiation can be used to set bounds on the spatial variation of $\mu$ at extragalactic scales \citep{Barrow05}. Surprisingly, a recent analysis of ammonia spectra in the Milky Way suggests a spatial variation of $\mu$ \citep{MLK09,Levshakov10b,Levshakov10}. The study, comparing the spectral lines of the ammonia inversion transition and rotational transitions of other molecules with different sensitivities to the parameter $\mu$, finds a statistically significant velocity offset that when interpreted in terms of a variation in $\mu$ gives $\Delta \mu/\mu = (2.2 \pm 0.7) \times 10^{-8}$. This will be the focus of the present paper. If we assume that the latter is not the result of some fluke and systematic experimental error, and thus take the result quite seriously, we are naturally led to the following question: What would be the simplest modification of our present physical theories that might account for such phenomena? One of the simplest possibilities one can think of is that the effective value of the coupling constants changes with space-time location. In this sense we note that, within the context of theories that are at the fundamental level background independent, the study of possible space-time dependence of fundamental constants is often considered as equivalent to the search for the existence of dynamical fields which couple to the gauge fields and/or to ordinary matter in ways that mimic the ordinary coupling constants.

There have been several proposals along those lines with various different motivations. Some of them arise from proposals for basic theories that arise in the search for unification of the four fundamental laws of physics such as string-derived field theories \citep{Wu86,Maeda88,Barr88,DP94,DPV2002a,DPV2002b}, related brane-world theories \citep{Youm2001a,Youm2001b,branes03a,branes03b} and Kaluza-Klein theories \citep{Kaluza,Klein,Weinberg83,GT85,Weson97}. There are also phenomenological models where a scalar field $\phi$ couples to the Maxwell tensor $F_{\mu\nu}$ and are characterized by Lagrangian density terms such as $- B_F(\phi) F_{\mu\nu}F^{\mu\nu}/4$ \citep{Bekenstein82,BSM02,OP02,Wetterich03,Dent09} and/or to the matter fields $\Psi_i$ ($i$ labeling the field flavor) as $ B_i(\phi) \bar \Psi_i {\Psi}_i $ \citep{OP02,KW04a,Brax04,BM05,OP08,Wetterich03,Dent09}. If these terms were added to the standard Lagrangian densities, it is quite clear that, in the first case, it would result in something like an effective fine structure constant $1/ e^2_{effective} = 1/ e^2 + B_F(\phi)$ while in the second case the effective masses of the elementary particles would be given by $m_i^{eff} = m_i+ B_i(\phi)$. Then, if the field took space-time dependent values, a feature that often requires the new field to be quite light so its value is not too rigidly tied to the minima of any self-interaction potential, then the effective fine structure constant and/or effective masses might look space-time dependent. This part of the story is quite clear, however, one can not focus on just this aspect of the theory when considering it. In fact, it is often the case that the most important bounds on the theory do not arise from the direct search for this dependence but from the effects of the direct exchange of quanta of this putative field would have on the behavior of ordinary matter. The fact that the scalar field must be light, as we have just described, indicates that this quanta exchange would not be drastically suppressed by a large mass in its propagator \citep{OscillatingBounds}. This generically leads to modification of the free fall and very often to signals that would mimic violations of the weak equivalence principle (WEP).

In fact, the connection between the theories involving spacetime dependence of coupling constants and the WEP was recognized already in the 1960s by Dicke \citep{dewitt1964relativity} (see page 163). More recently, \citet{Damour09} noted the extreme difficulty in mimicking the behavior of the mass of an object with any coupling that is not the gravitational one. This means that if the modified theory is required to be covariant, then the only possibility to ensure an exact compliance with the WEP is essentially to restrict the new fields to couple to matter in the same way as gravity, and this is impossible if we want the effective coupling constants to be scalar functions of some new fields. Consequently, one finds that the two issues are generically considered simultaneously in attempts to deal with possible variations of fundamental constants \citep{DP94,DPV2002a,DPV2002b,Bekenstein82,Barr88,Weson97,Scoccola07,BM05,OP02,Wetterich03,Dent09}. Thus, it is clear that we must face this connection in any attempt to deal with the problem at hand, and that in so doing we must consider the most modern and stringent bounds that are currently available on the possible violations of the WEP.

The basic idea of this manuscript is that rather than considering new fields which play the role of modifying the effective value of the fundamental constants, we can introduce non-standard aspects of well known fields in order to play that role. The long range fields in nature are the electromagnetic and gravitational ones. The use of the former in the desired context does not seem as a promising possibility because, for one it is very well understood and tested over very wide class of regimes, even at the quantum level (\textit{i.e.}, QED), and its enormous strength implies that any small modification would have very noticeable effects. The latter, on the other hand, is still far from being well understood (particularly its quantum aspects), and secondly, it seems conceivable that an exotic type of its coupling to matter might exist without having been detected so far. Considerations along these lines have led to proposals where the curvature of space-time might affect the propagation of matter fields in rather unusual ways \citep{QGP}, which might be viewed as violating of the strict letter of the equivalence principle without destroying overall general covariance\footnote{The proposals are motivated by searches for possible granular structure of space-time which might conflict with the ultra-locality that is implicitly assumed in the latter principle.}. In this manuscript we will explore this issue and show that, despite this early optimistic assessment, the value of the model's parameters needed to explain the observed variation of $\mu$ in the Milky Way can be ruled out by the bounds on those parameters emerging from experimental tests of the WEP.

The paper is organized as follows. In section \ref{datos} we discuss the astronomical data that suggest a variation of $\mu$ in the Milky Way. In Section \ref{models} we describe a general theoretical model with an ``exotic'' coupling between matter and gravitation; in section \ref{ginv} we present a subclass of models where gauge invariance is preserved. Section \ref{fits} is devoted to determine the value of a combination of the free parameters for the models presented in section \ref{models} as implied by the astronomical data discussed in section \ref{datos}. In section \ref{fits2} we obtain bounds on another combination of the free parameters of the models using the latest tests of the WEP. We end with a brief discussion and some conclusive remarks in section \ref{discusion}.

\section{Data discussion}\label{datos}

Astronomical spectroscopy can probe physical constants which describe atomic and molecular discrete spectra. It is the case of the electron-to-proton mass ratio $\mu$. Recently, \citet{Levshakov10b,Levshakov10} reported new bounds on $\mu$ obtained through the ammonia method. Previous bounds with the same method were obtained by \citet{MLK09}. The method consists in comparing the observed frequency of the ${\rm NH_3}$ inversion transition with a suitable rotational frequency of another molecule arising co-spatially with ammonia. In particular, the authors used precise molecular lines observed in Milky Way cold dark clouds to compare the apparent radial velocity for the ${\rm NH_3}$ inversion transition, $V_{inv}$, with the apparent radial velocity, $V_{rot}$, for rotational transitions in ${\rm HC_3N}$ and ${\rm N_2H^+}$ arising from the same molecular cloud. The method provides a relation between the shift radial velocity and the relative variation of $\mu$:
\begin{equation}
\frac{\Delta \mu}{\mu}=0.289\frac{V_{rot}-V_{inv}}{c}\equiv 0.289\frac{\Delta V}{c}.
\end{equation}

In the first report \citep{Levshakov10}, three radio telescopes were used to obtain the data: 32-m Medicina, 100-m Effelsberg and 45-m Nobeyama. The authors found several problems in treating the data: i) not all the molecular profiles can be described adequately with a single component Gaussian model and ii) molecular cores are not ideal spheres and being observed at higher angular resolutions exhibit frequently complex substructures. The line profiles may be asymmetric due to non-thermal bulk motions. Therefore the authors selected 23 pair molecular lines from the initially 55 molecular pairs observed. The weighed mean and errors as well as the robust M-estimate of the mean reported for the 23 data, for 100-m Effelsberg and 45-m Nobeyama data sets are shown in table \ref{velocidad}. In order to check if the weighed mean is a representative value of the mean value of each data set of table \ref{velocidad} we have calculated $\chi^2 = \sum_i p_i (x_i-x_W)^2/\sum_j p_j$ ($x_W$ is the weighed mean, $p_i = 1/\sigma_i^2$ and $\sigma_i$ refers to the $1 \sigma$ error reported in \citet{Levshakov10}) and compared it with the expected value of $\chi^2$ for a Gaussian distribution. In all cases the obtained value of $\chi^2$ is large compared with the expected value of $\chi^2$ for a normal (Gaussian) distribution. Therefore, the distribution of the data does not seem to be normal (or Gaussian). The authors of \citet{Levshakov10} have also calculated the robust M-estimate of the mean. From table \ref{velocidad} it follows that there is a significant difference between the Nobeyama weighed mean and robust M-estimator, while for the complete data set and Effelsberg data set both estimators are consistent. On the other hand, although there is a good agreement between the robust M-estimate of the three data sets, the Nobeyama data set has a larger systematic error due to the lower accuracy at the rest frequencies of the $\rm{N_2^+ H}$ transition. Furthermore, from the $23$ sources of the complete data set, $2$ are in common for the $3$ telescopes, and $3$ are in common for Effelsberg and Nobeyama. \citet{Levshakov10} treated them as independent \textit{observations} because the data was obtained with different instruments, \textit{so such data points have different systematic errors}; however, these are not independent as measurements of the physical quantity.

In a recent paper, the authors map four molecular cores selected from their previous sample in order to estimate systematic effects in the velocity offset and to check the reproducibility of the velocity offset on the year-to-year time base \citep{Levshakov10b}. Observations were performed with the 100-m Effelsberg telescope. In two cores the velocity offset can be explained by the observed kinematic structure. In the other two cores, they obtain a statistically significant positive velocity offset which is shown in table \ref{velocidad}. The differences between densities in clouds are not significant, thus the Effelsberg robust M-estimate reported in \citet{Levshakov10} is statistically more significant than the value reported in \citet{Levshakov10b}. We will consider the value of the Effelsberg robust M-estimate reported in \citet{Levshakov10} as the final value of the velocity offset between rotational and inversion transitions. However, it should be noted that taking the other values reported in the paper or taking the value reported in \citet{Levshakov10b} will not produce a significant change in the conclusions of the theoretical models we are testing in this paper.
\begin{table}[!ht]

\begin{center}
\renewcommand{\arraystretch}{1.5
\caption{Results for the relative radial velocity for all the data, only data obtained with Effelsberg radio telescope and only data obtained with Nobeyama radio telescope \citep{Levshakov10}. Entry name Effelsberg (2) refers to results reported in \citet{Levshakov10b}. $\left< \Delta W \right>_W$ refers to the weighed mean and corresponding error, $\left< \Delta V\right>_M$ refers to the robust M-estimate of the mean. All velocities are in units of \rm{m/s}.}
\label{velocidad}
\begin{tabular}{lcc}
\hline
\hline Data & $\left< \Delta V \right>_W$&$\left< \Delta V \right>_M$ \\ \hline
all&$20.7 \pm 3 $& $21.5 \pm 2.8$ \\
Effelsberg& $23.0 \pm 3.1$ &$23.2 \pm 3.8$\\
Nobeyama &$14.0 \pm 7.2$ & $22.9 \pm 4.2$ \\
Effelsberg (2)& & $26.9 \pm 4.2$ \\
\hline
\end{tabular}
\end{center}
\end{table}

\section{The Weyl Model}
\label{models}

The basic idea of the Weyl models involves considering the effect described in the previous section as due to a non-minimal and rather exotic coupling of the matter fields with gravity through the Weyl tensor. At first sight this may seem as an unnatural proposal as gravity is usually neglected in these regimes and, furthermore, one generally does not feel that there might be a fundamental reason to couple gravity with matter fields in exotic ways. However, in contrast to other models \citep{KW04a,Brax04,MS07,OP08}, in this scheme there is no need to invoke new unobserved dynamical fields, or non-dynamical fields that break Poincar\'e invariance, or other such problems, which in our opinion are ruled out by the analysis of its consequences on virtual particles \citep{Collins04, Collins06}. Moreover, as described in \citet{QGP}, Bonder and Sudarsky \cite{BS08,BS09,BS10}, the generic view we adopt in considering these sort of models is that gravity as ``the curvature of a manifold'' is only an effective description of more fundamental degrees of freedom (from an unknown quantum theory), and thus, the unnaturalness of the coupling terms is tied to the need to use a metric description rather than the still unknown language of the quantum gravity theory.

The task is to find a way in which gravity and matter could interact in a phenomenological level causing the described change in the matter's observed mass. We study fermionic matter fields $\Psi_i$ where $i$ labels the field's flavor. The usual mass term in this case is $m_i \bar{\Psi}_i \Psi_i$, therefore, in order to explain the observations we need to replace $m_i$ by some scalar depending on the gravitational environment. This should be implemented by
\begin{equation}
m_i \rightarrow m_i \left[1+\xi_i f\left(\frac{R}{\Lambda^2}\right)\right],
\end{equation}
where $\xi_i$ are small phenomenological parameters which may be different for each flavor, $f$ is a function of the curvature tensor and $\Lambda$ is an energy scale. The first scalar function that comes to mind is the Ricci scalar $R$. However, the Ricci tensor at a space-time point $x$, and thus the Ricci scalar, is completely determined by the matter at the same point $x$, which implies that coupling $\Psi_i(x)$ with $R(x)$ is a self-coupling that, for phenomenological purposes, is not interesting. Thus, we should build $f$ with what is left when removing from the Riemann tensor the part determined by the Ricci tensor: The Weyl tensor, $W_{abcd}$. It is trivial to note that it is not possible to construct a Lorentz scalar out of one power of $W_{abcd}$, thus, the simplest scalar one can write is $f=\frac{W_{abcd} W^{abcd}}{\Lambda^4}\equiv\frac{W^2}{\Lambda^4}$. Note that $\Lambda^4$ has (mass) dimensions $4$ ($\hbar$ and $c$ are taken to be $1$) and thus $\xi_i$ are dimensionless parameters.

\subsection{Gauge invariance and the Weyl proposal}
\label{ginv}

We might be concerned that if the coupling to different fermion flavors are completely arbitrary, the new proposal conflicts with gauge invariance destroying the renormalizability of the theory. However, the fact that this modification is assumed to be tied to gravitational sector itself, known to be non-renormalizable\footnote{Nowadays this is simply taken as indication that a theory must be
considered as an effective theory.}, places these concerns in the proper perspective. In addition, a simple restriction can remove this issue altogether: The quantity that we need to vary with location is the ratio of electron and proton masses $\mu$, and the two quantities have a rather different origin; one arising mainly a from the strong interaction of the quark and gluon constituents, while the other is due to the Yukawa term and the vacuum expectation value of the Higgs field. We can construct a rather simple theory with the required features by focusing on this difference. Namely, we assume that the exotic coupling to Weyl tensor enters into the theory only in the fermionic Yukawa terms as
\begin{equation}
{\cal L}_{Yukawa}=\sum_{a,b} \Gamma_{ab} (W^2) \bar{\psi}_{a}^{R} \Phi \Psi_b^{L} + h.c.
\end{equation}
where the sum is over the standard SU(2) left handed fermion field doublets $\Psi_b^{L}$ and right handed fermion filed singlets $\psi^R_b$ of the electroweak theory, $\Phi$ is the standard doublet Higgs doublet field, and $\Gamma_{ab}(W^2) $ are the Yukawa coefficients which are now considered to depend on the local value of the magnitude of Weyl tensor. Considering the lowest order of a Taylor expansion we can write $\Gamma_{ab} (W^2) = (1+{\xi}_{(ab)} W^2 /\Lambda^4 )\Gamma_{ab} (0)$, where the parameters $\xi_{ab}$ characterize the leading order corrections and, in general, depend on the flavor indices $a$ and $b$. As we only focus on the first generation and the quark sector is, in any event, quite different phenomenologically from the lepton sector due to the connection with the strong (or color) force, we ignore that possible dependence and simply write $\xi$. It is thus clear that as the Higgs field acquires a vacuum expectation value, the resulting fermionic mass matrices are multiplied by the factor $1+ \xi W^2/\Lambda^4$ and thus the lepton and quark masses appear in the effective Lagrangian as multiplied by that factor.

It is then evident that, as long as gravitation is considered as a fixed background and one does not attempt to introduce radiative corrections involving gravitons, the theory preserves the same renormalizability properties as the standard model. Of course, when attempting to include radiative corrections involving gravitons, the theory will encounter the usual problems facing the quantization of gravitation. On the other hand, the mass of any lepton appears in the low energy theory as
\begin{equation}
m_l\rightarrow m_l \left(1+ \frac{\xi}{\Lambda^4}W^2 \right),
\end{equation}
while in the case of hadrons only the valence quark mass, which represents a rather small part of the total hadron mass, suffers such modification, thus we can expect that the change in the total hadron mass is much smaller, \textit{i.e.},
\begin{equation}
m_h\rightarrow m_h \left(1+ a_h \frac{\xi}{\Lambda^4}W^2 \right),
\label{hadron}
\end{equation}
where the parameter $ a_h$ is of order $10^{-3}$ (the valence quark proportion of the hadronic masses). In other words, we can take the electron's $\xi_e$ parameter to be identical to the fundamental parameter $\xi$, while the parameter corresponding to the proton $\xi_p$ is given by $ a_p \xi$ where $ a_p$ is of order $10^{-3}$. The model described in this subsection is less general than the case where the proton and electron coupling to the Weyl tensor are not related. In the following sections we consider the theoretical predictions of both models and compare it with astronomical and experimental data. We call model I to the model where the couplings of leptons and quarks are not related and model II to the case where the exotic coupling enters only in the Yukawa terms.

\section{Estimates on the Weyl model parameters from electron-to-proton mass ratio}\label{fits}

In a model where the interaction between gravity and matter at a phenomenological level is such as described in section \ref{models}, the low energy limit of the interaction Lagrangian density can be written as
\begin{equation}\label{lag}
{\cal L}_{int}= \frac{\xi_e}{M_P^4} m_e W^2 \bar \psi_e \psi_e
+\frac{\xi_p}{M_P^4} m_p W^2 \bar \psi_p \psi_p +\frac{\xi_n}{M_P^4} m_n W^2 \bar \psi_n \psi_n,
\end{equation}
where subscripts $e,p,n$ respectively refer to electrons, protons and neutrons and we set $\Lambda$ equal to the Planck scale $M_P = 1.22 \times 10^{19}$ GeV. Thus, the effective masses of the particles can be expressed as $m^{eff}_i = m_i\left(1 + \frac{\xi_i}{M_P^4} W^2\right)$. Note that we have restored in the notation $\xi_i$ the possible dependence of the parameter on the particle type. The point being, as discussed in the previous section, that even if $\xi$ had a single value throughout the first generation of quarks and leptons, the composite nature of the hadrons, as well as the gluon contribution to the mass of these particles, would make the effective value quite different from that corresponding to electrons. The difference between neutrons and protons can similarly be expected to arise as a result of the electromagnetic differences between them.

According to the model we are proposing, the observable quantity is
\begin{equation}
\mu^{eff} \equiv \frac{m^{eff}_e}{m^{eff}_p}=\frac{m_e}{m_p}\left(\frac{1+\frac{\xi_e}{M_P^4}
W^2}{1+\frac{\xi_p}{M_P^4} W^2}\right) \simeq \frac{m_e}{m_p}\left(1+\alpha \frac{W^2}{M_P^4}\right),
\end{equation}
where in the last step we use the fact that $\xi_i/M_P^4$ are small and we define $\alpha\equiv \xi_e -\xi_p $ for model I and ${\alpha\equiv (1-a_h)\xi}$ for model II. Thus, the observed electron-to-proton mass ratio value in cold molecular clouds respect to the same value at Earth is
\begin{equation}\label{Delta mu eff}
\frac{\mu^{eff}_{cl}-\mu^{eff}_{\oplus}}{\mu^{eff}_{\oplus}} \simeq \frac{\alpha}{M_P^4} (W^2_{cl} -W^2_{\oplus}),
\end{equation}
where the subscripts $cl$ and $\oplus$ stand respectively for the interstellar clouds and the Earth.

In order to compute the Weyl tensor both, in the cloud and at the laboratory, we consider a sphere of constant density $\rho$ and radius $R$, surrounded by vacuum. Outside the sphere, namely at a distance $r\geq R$ from its center, we get $W^2=\frac{48 G^2 M^2}{r^6}$ where $G$ is Newton's constant and $M \equiv 4 \pi \rho R^3/3$ is the mass of the sphere. Due to the dependence of $W^2$ on $r$ we realize that we have to be careful while considering the contributions to $W^2$. For example, contributions from massive bodies near to the laboratory such as a wall may be greater than the contribution of the entire Earth. Therefore, we also consider  the contribution of the Earth and from walls of $4$ m height, $4$ m width and $0.5$ m depth made of cement and iron located $0.1$ m from the experiment.  The results of the calculations show that the contribution to $W^2$  of an iron wall   located very close to the experiment ($W^2_{I}= 1.8 \times 10^{-44} {\rm m^{-4}}$) is   greater than the effect of the Earth ($W^2_{\oplus} = 1.4 \times 10^{-44} {\rm m^{-4}}$)  or a cement wall ($W^2_{C} = 1.4 \times 10^{-44} {\rm m^{-4}}$) .
Given that we have not the exact details of the laboratory where the rest wavelengths are measured and since we are dealing with a positive detection of $\Delta\mu^{eff}$, we consider the iron wall to estimate a lower bound on $\alpha$. On the other hand, the contribution of the entire Earth can not be neglected and therefore we consider this contribution to obtain an upper bound on $\alpha$. Even though the wall is not a sphere, in the regime we are working on it is possible to use the approximation of linearized gravity where the superposition principle is valid. Thus, the contribution of the wall can be regarded as the sum of the contribution of a great number of spheres. We estimate that differences with the exact calculation may be at most of order $10$.

We model the interstellar clouds as spheres with constant density $\rho$ and radios $R$. In order to calculate $W^2$ we use the inner Schwarzschild solution, which leads to $W^{2}_{cl}=\frac{4^{6}\pi^{4}r^{4}\rho^4 G^{4}(3-8\pi
R^{2}\rho G)^{2}}{3(3-8\pi r^{2}\rho G)^2(3-2\pi \rho G(3r^{2}+R^{2}))^{2}}$. Taking for the clouds a mean density $\rho_{cl}= 1.5\times 10^{27}$ m$^{-4}$, a mean radius $R_{cl} = 0.052$ Pc \citep{MLK09} and at $r=R/2$ we obtain $W^{2}_{cl}\simeq 1.5 \times 10^{-178}\ \rm{m}^{-4}$.(We also calculated the value of $W^{2}_{cl}$ for $r=0$ and $r=R$ and the results do not differ significantly from the case $r=R/2$). Taking the value of Effelsberg robust mean discussed in section \ref{datos} we get
\begin{equation}\label{ineq alpha1}
1.22 \times 10^{36}\rm{m}^4 \le \frac{|\alpha|}{M_P^4} \le 1.57 \times 10^{36}\rm{m}^4
\end{equation}
On the other hand, if we consider the model described in section \ref{ginv}, where hadron masses suffer much smaller modifications than lepton masses, in order for the model to explain the observations described in section \ref{datos} the following condition must be fulfilled:
\begin{equation}\label{ineq alpha2}
1.22 \times 10^{36}\rm{m}^4 \le \frac{|\xi|}{M_P^4} \le 1.57 \times 10^{36}\rm{m}^4
\end{equation}
The next section is dedicated to study if the relations (\ref{ineq alpha1}) and (\ref{ineq alpha2}) are compatible tests of the WEP.

\section{Bounds from E\"otvos type experiments}\label{fits2}

The gravitational potential of an object composed of $N$ atoms with atomic number $Z$ and baryon number $B$ can be written by taking into account that the effective masses are modified in this model according to equation (\ref{lag}):
\begin{eqnarray}
V &=& N Z \frac{\xi_e}{M_P^4} m_e W^2 + N Z \frac{\xi_p}{M_P^4} m_p W^2 + N
(B-Z) \frac{\xi_n}{M_P^4} m_n W^2 \nonumber \\
&=& N \frac{\alpha'} {M_P^4} W^2 m_p,
\end{eqnarray}
where
\begin{equation}
\alpha' =Z \left(\xi_e \frac{m_e}{m_p} + \xi_p - \xi_n \frac{m_n}{m_p}
\right)+B\xi_n \frac{m_n}{m_p}.
\end{equation}
The force acting on a freely falling body of mass $M_b$ can be obtained from $\vec F= -\vec \nabla V$, thus, the respective acceleration is $\vec{a} = -\frac{N \alpha' m_p \vec{\nabla} W^2}{M_P^4 M_b}$. The differential acceleration of two bodies with different composition but the same number of atoms $N$ is $\Delta \vec{a}= -\frac{\vec{\nabla} W^2\ m_p (\delta_1\ \alpha_1' + \delta_2 \alpha_2')}{M_P^4}$, where we assume that the mass of the body can be expressed as $N m$ with $m$ the atomic mass of each body and we define
\begin{equation}
\alpha_1'\equiv \xi_e \frac{m_e}{m_p} + \xi_p - \xi_n \frac{m_n}{m_p}, \qquad
\alpha_2'\equiv \xi_n \frac{m_n}{m_p},
\end{equation}
and also $\delta_1\equiv (Z_1/m_1)-(Z_2/m_2)$ and $\delta_2\equiv (B_1/m_1)-(B_2/m_2)$, the subscripts indicating which object we are considering. On the other hand, if hadron masses suffer smaller modifications as suggested in section \ref{ginv}, the prediction of the violation of WEP can be written in terms of a single parameter $\xi$: $\alpha_1'\equiv 10^{-3} \xi$, $\alpha_2'\equiv 10^{-3} \xi$.

Limits on violations of the WEP come from E\"otvos-Roll-Krotkov-Dicke and Braginsky-Pannov measurements of the differential acceleration of test bodies. The most stringent limits are obtained from measurements of differential acceleration towards the Sun. However, since the force resulting from Weyl models is a short range force, the relevant bounds to test such models are provided by measurements towards the Earth. In this kind of experiments a continuously rotating torsion balance instrument is used to measure the acceleration difference of test bodies with different composition. In table \ref{tabla1}, we summarize the current bounds considered in this paper as well as the composition of the test bodies in each case. All bounds were obtained from an experiment made at the University of Washington Nuclear Physics Laboratory. The authors of the experiment mention two sources for the relevant signals and in our case these would be the sources for $\vec{\nabla} W^2$, a hillside of $26$ m located closed to the laboratory (we estimate $1$ m) and a layer of cement blocks added to the wall of the laboratory \citep{Adelberger90}. Here we model these sources simply as spherical masses at a given distance. Barring some fortuitous cancellation among the various contributions in the laboratory due to their detailed location and geometry, these simplification cannot amount to more than a change by a geometrical factor of order one, which we ignore. Moreover, the expression for $\left|\vec{\nabla} W^2\right|$ reads $\left|\vec{\nabla} W^2\right|=\frac{288}{r^7} \left(\frac{4}{3}\pi G \rho R^3\right)^2$, and therefore the contribution from the hillside and cement layer can be estimated to be, respectively, $\left|\vec{\nabla} W^2_{\rm{H}}\right|= 7.9 \times 10^{-45} \rm{m}^{-5}$, and $\left|\vec{\nabla} W^2_{\rm{L}}\right|=2.9 \times 10^{-45} \rm{m}^{-5}$. Since both contribution are of the same order, we continue our analysis considering only the contribution of the hillside. In order to estimate the value of $\alpha'_1$ and $\alpha'_2$, we perform a least square minimization using the calculated expression for the acceleration and the data from table \ref{tabla1}. We obtain for model I:
\begin{eqnarray}
\frac{\alpha_1'}{M_P^4}&=&(3.5 \pm 7.8) \times 10^{13}\ \rm{m}^4,\\
\frac{\alpha_2'}{M_P^4}&=&(-0.7 \pm 2) \times 10^{15}\ \rm{m}^4.
\end{eqnarray}
Moreover, using the data from table \ref{tabla1} we get for the model II:
\begin{equation}
\frac{\xi}{M_P^4} = (1.7 \pm 6) \times 10^{16} \rm{m}^4.
\end{equation}
The contribution of other massive sources such as other hills close to the laboratory would give lower values of $\left|\vec{\nabla} W^2\right|$ resulting in higher values for $\xi$. However, we cannot ignore the hill of $26$ m and the layer of cement mentioned by \citet{Adelberger90} and therefore the limit on $\xi$ is not a lower bound in the sense that it has been discussed in section \ref{fits} for the astronomical data.

\begin{table*}[!ht]
\renewcommand{\arraystretch}{1.3}
\caption{Bounds on the equivalence principle considered in this paper.$(B/{\rm m})_i$ and $M({\rm g})$ refer to the baryon number and total mass of the test bodies }
\label{tabla1}
\begin{center}
\begin{tabular}{ccccccc}
\hline
\hline
$\Delta a \,\, (\rm{m}^{-1})$&$Z_1$&$Z_2$&$(B/{\rm m})_1$&$(B/{\rm m})_2$&$M
({\rm g})$ & Reference\\\hline
 $(2.5 \pm 7.23) \times 10^{-31}$&$ 4$&$13$&$0.998648$&$1.000684$&$10$
&\cite{Su94}\\
 $(5.3 \pm 6.45) \times 10^{-31}$&$4$&$29$&$0.998648$&$1.001117$&$10$
&\cite{Su94}\\
 $(0.67 \pm 3.45) \times 10^{-32}$&$4$&$22$&$0.99868$&$1.001077$&$4.8$ &
\cite{Schlamminger08} \\
 $(2.3 \pm 2.3) \times 10^{-30}$&$4$&$13$&$0.998648$&$1.000684$&$10$
&\cite{Adelberger90}\\
 $(0.9 \pm 1.89) \times 10^{-30}$&$4$&$29$&$0.998648$&$1.001117$&$10 $ &
\cite{Adelberger90}\\
 $(1.11 \pm 3.11) \times 10^{-32}$&$29$&$82$&$1.0011166$&$1.0001694$&$10$&
\cite{Smith00}\\\hline
\end{tabular}
\end{center}
\end{table*}

\section{Discussion and Conclusions}
\label{discusion}

We have considered a model where a non-minimal coupling of Weyl tensor to matter would result in an effective mass for fermionic fields which would be space-time dependent. The model might be considered as a possible explanation for the recent reported observations of a space-time variation of the electron-to-proton mass ratio in \cite{MLK09,Levshakov10b,Levshakov10}. We have developed the model in some detail and extracted the range of values for a combination of the parameters which would be necessary to account for the ``exotic'' observation. We have also considered the constrains on the model that arise from consideration of precision tests of the WEP and used some of the most modern relevant data to contrast the two results.

In section \ref{fits} we analyzed the variation of the effective electron-to-proton mass ratio between the environments corresponding to the Earth's surface and a molecular cloud and showed that consistency with the reported data requires $\alpha = \xi_e - \xi_p \sim 10^{36}\  \rm{m}^4\  M_P^4$. On the other hand, the result of the statistical analysis performed in section \ref{fits2} with bounds on the WEP constrain the value of $\alpha_1'=\xi_e m_e/m_p + \xi_p - \xi_n m_n/m_p$ to be of order $10^{13}\  \rm{m}^4\  M_P^4$. Here one might be inclined to note that the two situations are sensitive to slightly different combination of the fundamental parameters, and is thus conceivable that $\alpha$ might be as large as required to account for the astronomical observations while $\alpha_1'$ is as small as needed to conform with the laboratory bounds. We view such possibility as very unlikely as it would imply that the particular choice of materials compared in the the laboratory tests were coincidentally those for which the signal resulting from the generic couplings happened to cancel out almost exactly (at the level of one part in $10^{17}$). Moreover, as similar tests with slightly lower precision do exist for other materials, taking this line of reasoning would only lead to a reduction by at most a couple of orders of magnitude in the constraint, something which would still be sufficient to rule the model out.

If we consider the model described in subsection \ref{ginv}, consistency with the astronomical data requires $|\xi| \sim 10^{36}\  \rm{m}^4\  M_P^4$ whereas the result of the statistical analysis performed in section \ref{fits2} with bounds on the WEP constrain the value of $\xi $ to be of order $10^{16}\  \rm{m}^4\  M_P^4$. We thus have found that barring some miraculous cancellation or some unnatural fine tuning of the experimental conditions and/or of the model, the two sets are incompatible, and that a model where the variation of the electron, proton and neutron effective mass is driven by the scalar magnitude of the Weyl tensor can not account for the experimental constraints and the observational data and should be also ruled out.

If the observations of \cite{MLK09,Levshakov10b,Levshakov10} were to be further confirmed and the evidence for a change in the value of the electron-to-proton mass ratio became incontrovertible, one would need some different sort of explanation, however, due to the connection between space-time dependency of parameters and the couplings of ordinary matter with dynamical fields, which appears inherent of background independent theories, it seems very unlikely that one might find a model where the bounds imposed by tests of the WEP would not be of great relevance and impact. Nonetheless, we should stress our belief that this kind of models should be further explored, not only as potential explanatory grounds for atypical observations, but also as leading to robust constraints on the possible nontrivial couplings of matter and gravitation.

\section*{{\bf Acknowledgments}}
We gladly acknowledge very useful discussions with E. Fischbach and P. Molaro. Y.B. and D.S. acknowledge support from PAPPIT project IN 119808 and CONACYT project No 101712. D.S. was supported in part by sabbatical fellowships from CONACYT and DGAPA-UNAM and the hospitality of the IAFE. S. J. L. and F.A. T. P. are supported by PICT 2007-02184 from Agencia Nacional de Promoci\'{o}n Cient\'{i}fica y Tecnol\'{o}gica, Argentina and by PIP N 11220090100152 from Consejo Nacional de Investigaciones Cient\'{\i}ficas y T\'{e}cnicas, Argentina.