\section{Introduction}
It is found that the power spectra of the {\rm H~{\sc i}~} column density fluctuations in our Galaxy \citep{1983A&A...122..282C, 1993MNRAS.262..327G} follow power law  at length scales ranging a few parsec to few hundred parsecs.  Different probes to the statistics of the turbulent velocity fluctuations of our Galaxy; statistics of centroid of velocities \citep{2009RMxAC..36...45E}, velocity coordinate spectrum or VCS \citep{2009RMxAC..36...54P}, velocity channel analysis or VCA \citep{2000ApJ...537..720L}, spectral correlation function \citep{2001ApJ...555L..33P} show  that at  similar scales, the power spectrum of the velocity fluctuations also follow power laws. Theoretically, these structures are understood in terms of compressible fluid turbulence \citep{2004ARA&A..42..211E} driven by supernovae shocks, self gravity etc.  Recently, {\rm H~{\sc i}~} column density power spectrum from external dwarf \citep{2006MNRAS.372L..33B, 2008MNRAS.384L..34D} and spiral galaxies \citep{2008MNRAS.384L..34D, 2009MNRAS.397L..60D, 2013NewA...19...89D} are also found to follow power law at  length scales ranging from a few parsec to $\sim 10$ kpc.
Generating mechanism of these large scale structures are yet to be understood \citep{2014MNRAS.441..525W}. Estimating the velocity fluctuation power spectrum  would give more clue here. For this purpose, VCA is the most efficient technique and has been used by several authors successfully for our galaxy and nearby dwarf galaxies \citep{2001ApJ...551L..53S,2015ApJ...810...33C}.  However, for the external spiral galaxies, this technique has limited use \citep{2015MNRAS.452..803D}.

In this letter we propose to use moments of the observed visibilities to perform a model dependent estimate of the statistics of the line of sight component of turbulent velocities from external galaxies. We outline an estimator for the power spectrum of the turbulent velocity fluctuations and  access it's efficacy for a noise free simulated data.
Section (2) defines visibility moments and  the estimator. Simulation parameters are given  in section (3) and in section (4) we explain the procedure to use the estimator on the simulated data. We discuss the results  and conclude in the last section (5).

\section{Power spectrum of visibility moments}
Specific intensity of radiation with rest frequency $\nu_{0}$ from {\rm H~{\sc i}~} gas at $\vec{r}= (x, y, z) = (\vec{R}, z)$ at a distance $D$ and  having temperature $T$ can be written for small optical depth limit as  \citep{2011piim.book.....D}
\begin{equation}
I (\vec{\theta}, v)\ =\ I_{0}\,  \int d\,  z\, n_{HI} (\vec{r})\, \phi (v).
\end{equation}
Here $v = c ( \nu_{0} - \nu)/nu_{0}$,  $\nu$ is the frequency of observation,  $I_{0} = \frac{3 h \nu_{0} A_{21}}{16 \pi}$ and $z$ is the line of sight direction. $\vec{R}$ is in the plane of the sky, angular separation of a direction in the field of view with respect to the field centre is $\vec{\theta} = \frac{\vec{R}}{D}$. The line shape function $\phi(v)$ can be written as
\begin{equation}
\phi(v) = \phi_{0}\, \exp \left [ - \frac{ \left ( v - v_{z} (\vec{r}) \right ) ^{2} } {2 \sigma^{2}(\vec{r}) } \right ],
\end{equation}
were $\sigma(\vec{r}) = \sqrt{\frac{ k_{b} T}{m_{HI}}}$ \footnote{$k_{b}$ : Boltzmann constant, $m_{HI}$ : mass of hydrogen atom} is the thermal velocity dispersion and $v_{z}(\vec{r})$ is the line of sight velocity component of the {\rm H~{\sc i}~} gas.

We use radio interferometers to observe {\rm H~{\sc i}~} emissions from the  spiral galaxies.
Every antenna pair of an interferometric array measures the quantity termed as  visibility $V({\vec{ U}}, v)$ at a baseline ${\vec{ U}}$ given by the projected antenna separation at the plane of the sky in units of observing wavelength.  Visibilities are direct Fourier transform of the specific intensity distribution, i.e,
\begin{equation}
V({\vec{ U}}, v) \ =\ \int d {\vec \theta}\  e^{2 \pi i {\vec \theta} . {\vec{ U}}} I (\vec{\theta}, v).
\end{equation}
Depending on the antenna positions of the particular telescope array, integration time of observation and declination of the source, visibilities can only be measured at a finite set of  points in the ${\vec{ U}}$ plane. While reconstructing the image by performing an inverse Fourier transform of the observed visibilities, above sampling results in a complicated telescope beam pattern, usually known as dirty beam \citep{1999ASPC..180.....T}. Hence, image obtained by an inverse Fourier transform of the visibilities is a convolution of the dirty beam with the sky specific intensity. Various  deconvolution techniques are used  to estimate the specific intensity from the image.  However, all  these algorithms effectively interpolate the visibilities in the unsampled part of the  ${\vec{ U}}$ plane. This gives rise to  an unknown but  correlated noise   in the reconstructed image \citep{2011arXiv1102.4419D}. Hence, any estimate of power spectrum of  column density or  line of sight velocity,  evaluated using the reconstructed specific intensity is expected to  have an unknown bias.

We define the zeroth moment of the visibility function as follows
 \begin{equation}
{\cal V}_{0}(\vec{U}) = \int \ dv\  V ({\vec{ U}}, v).
\end{equation}
The  velocity  integral above has to be performed over the entire spectral range of {\rm H~{\sc i}~} emission. Clearly the moment ${\cal V}_{0}(\vec{U}) $ is proportional to the the Fourier transform of columns density $N_{HI}({\vec \theta}) = \int dz\, n_{HI}(\vec{R}, z)$. In fact, the quantity ${\cal P}_{0}(\vec{U}) = \langle \mid {\cal V}_{0}(\vec{U}) {\cal V}_{0}(\vec{U})^{*} \mid \rangle $ has been used in literature \citep{1993MNRAS.262..327G, 2013MNRAS.436L..49D} to estimate the power spectrum of the column density fluctuations of our Galaxy and external galaxies. We define the first moment of visibility function as
\begin{eqnarray}
 {\cal V}_{1}(\vec{U}) &=& \int \ dv\  v\ V ({\vec{ U}}, v) \\ \nonumber
 &=& \int d {\vec \theta}\ e^{2 \pi i {\vec \theta} . {\vec{ U}}}\  I_{0}\,  \int d\,  z\, n_{HI} (\vec{r})\, v_{z}({\vec r}).
\end{eqnarray}
Here $v_{z} ({\vec r})$, the  line of sight velocity, has two components:  $v^{\Omega}_{z}({\vec \theta})$,  systematic rotation of the galaxy,  and  $v^{T}_{z}({\vec \theta})$,  random motion because of turbulence. Hence, $v_{z}({\vec r}) = v_{z}^{\Omega}({\vec r}) + v_{z}^{T}({\vec r})$. In case of the spiral galaxies with relatively small inclination angle, we may assume that the line of site component of the rotational velocity is independent of the distance along the line of sight, then
\begin{equation}
{\cal V}_{1}(\vec{U}) \ =\ {\cal V}_{0}(\vec{U}) \otimes  \tilde{v_{z}^{\Omega}} ({\vec{ U}}) \ +\ {\tilde{\chi}}(\vec{U}),
\end{equation}
where $\tilde{v_{z}^{\Omega}} ({\vec{ U}})$ is the Fourier transform of $v_{z}^{\Omega}({\vec r})$ and ${\tilde{\chi}}(\vec{U})$ is the Fourier transform of
\begin{equation}
{\chi}(\vt) = I_{0}\ \int dz\  n_{HI}({\vec r}) v_{z}^{T}({\vec r}).
\end{equation}

\subsection{An estimator of the turbulent velocity fluctuations}
It is safe to assume that the systematic rotational velocity of the galaxy and the random turbulent velocity are not correlated. We define
\begin{equation}
{\cal P}_{\chi}(\vec{U}) =\langle \mid {\tilde{\chi}}(\vec{U}) \mid^{2} \rangle =\langle  \mid {\cal V}_{1}(\vec{U}) \mid^{2} -  {\cal P}_{0}(\vec{U}) \otimes \mid \tilde{v_{z}^{\Omega}} ({\vec{ U}})  \mid ^{2} \ \rangle.
\end{equation}
The angular brackets $< ... >$ denote ensemble averaging, which has to be taken over many statistical  realisations of the observed sky. As we are provided with only one realisation in reality, this can not be performed. Turbulence generated fluctuations in density and velocity, particularly for magnetohydrodynamic (MHD) turbulence, can have preffered orientations \citep{2005ApJ...631..320E}. However, it is not clear, if at the length scales of $1$ kpc to $10$ kpc or higher one would expect to see anisotropic fluctuations in neutral hydrogen resulting because of MHD turbulence. Here, we assume that we are probing homogenious and isotropic fluctuations only. The method presented here can be generalized for the anitropic case, however, observations with higher signal to noise than presently available may be required.  With this assumptions, the ensemble average in eqn.~(8) can be  replaced by averages over annular bins in  ${\vec{ U}}$. Note that the quantities $ \mid {\cal V}_{1}(\vec{U}) \mid^{2}$  and ${\cal P}_{0}(\vec{U})$ can be estimated directly from the observed visibilities. The line of sight component of the rotational velocity $v^{\Omega}_{z}({\vec \theta})$, inclination or position angle etc depends only on the galacto-centric radius.  These quantities,  and hence also $\mid \tilde{v_{z}^{\Omega}} ({\vec{ U}})  \mid $, can be estimated from the reconstructed image without any bias. It is clear from eqn.~(7) that ${\cal P}_{\chi}(\vec{U})$ contain in it informations about  {\rm H~{\sc i}~} density as well as {\rm H~{\sc i}~} velocity fluctuation power spectra, whereas the {\rm H~{\sc i}~} density fluctuation power spectra can be estimated independently using ${\cal P}_{0}(\vec{U})$.

Next we assume that the turbulent velocity and the density fluctuations are not correlated and refer to the following two cases.

\begin{itemize}
\item If we are interested in estimating the power spectra at length scales comparable to the scale of the galaxy's disk, i.e, scales much larger than the scale height of the galaxy in consideration, we may approximate the galaxy as a thin disk. In this case, we can write
\begin{equation}
{\cal P}_{\chi}(\vec{U}) \ ={\cal P}_{0}(\vec{U}) \otimes {\cal P}_{v}^{T}(\vec{U}),
\end{equation}
where ${\cal P}_{v}^{T}(\vec{U}) = \langle \mid  v_{z}^{T} ({\vec{ U}}) \mid^{2} \rangle$ is the power spectrum of the turbulence velocity. It is to be noted that the quantities ${\cal P}_{0}(\vec{U})$ and ${\cal P}_{v}^{T}(\vec{U})$ are two dimensional power spectra of the column density and the velocity fluctuations respectively. We have already discussed that ${\cal P}_{0}(\vec{U})$ can be measured from the zeroth moment of the observed visibilities. Since we expect all three functions in the above expression to be smooth and continuous, we can estimate  ${\cal P}_{v}^{T}(\vec{U})$  either by deconvolving ${\cal P}_{0}(\vec{U})$ from ${\cal P}_{\chi}(\vec{U})$ or by regression analysis using a parametric model for ${\cal P}_{v}^{T}(\vec{U})$. Since we have an independent and exact measure of ${\cal P}_{0}(\vec{U})$ from the zeroth moment of visibilities, we can have  a complete estimate of ${\cal P}_{v}^{T}(\vec{U})$ in this case.

\item On the other hand, if we are to estimate the power spectra at scales comparable to the scale height of the galaxy's disk, we may write,
\begin{equation}
{\cal P}_{\chi}(\vec{U}) \ = P^{3D}_{n_{HI}} ({\vec{ U}}, 0) \otimes P^{3D}_{v^{T}_{z}}({\vec{ U}}, 0),
\end{equation}
where $P^{3D}_{n_{HI}} ({\vec{ U}}, 0)$ and $P^{3D}_{v^{T}_{z}}({\vec{ U}}, 0)$ are the three dimensional power spectra of $n_{HI}({\vec r})$ and $v^{T}_{z}({\vec r})$ evaluated at the zero line of sight wavenumber. It can be shown that $P^{3D}_{n_{HI}} ({\vec{ U}}, 0)$ is  proportional to  ${\cal P}_{0}(\vec{U})$ for a thick disk. Hence, just like the  case for the thin disk, we may use ${\cal P}_{0}(\vec{U})$ to find the turbulent velocity fluctuation power spectrum from above relation either by deconvolution or by regression analysis. However, our ability to estimate the amplitude of $P^{3D}_{v^{T}_{z}}$ would depend on the details of the morphology of the galaxy in the line of sight direction.
\end{itemize}

\section{Simulating {\rm H~{\sc i}~} observations from a spiral galaxy}

\begin{figure}
\begin{center}
\epsfig{file=fig1.eps,width=2.8in}
\end{center}
\caption{Radial profile (dash line), vertical profile (dot-dash line) of the {\rm H~{\sc i}~} number density  adopted for simulation are plotted against distance from the centre of the galaxy in radial and vertical directions. The solid black line represents the rotation curve used to generate $v^{\Omega}_{z}({\vec \theta})$. Note distances are normalised to the size of the simulation box.}
\label{fig:profile}
\end{figure}
A major assumption we took in deriving the expression for ${\cal P}_{\chi}(\vec{U})$  in the previous section is that $v^{\Omega}_{z}$ is independent of the line of sight direction. In case of galaxies with lower inclination angle, this assumption is justified. However, for  galaxies with higher inclination angles, line of sight pass through different galacto-centric radius and this assumption may not be justified.
We simulate   {\rm H~{\sc i}~} observation from an external  spiral  galaxy and use it to  study the effect of galaxy's inclination on the turbulent velocity power spectrum estimator  described in the previous section. Simulation of the {\rm H~{\sc i}~} signal  is based on eqn.~(1) and (2) and discussed in \citet{2015MNRAS.452..803D} in detail, we note the salient features here.  We consider an exponential radial profile and a gaussian vertical profile to model the {\rm H~{\sc i}~} distribution in the galaxy, i.e,
\begin{equation}
\begin{aligned}
n_{HI} (\vec{r_{g}})  \ =\ n_{0}\  & \exp\left (- \mid \vec{R}_{g} \mid / R_{0}  \right )\  \exp \left ( -\frac{1}{2} \left( z/z_{0} \right )^2 \right ) \\
&\left [ 1 + f_{HI}\ \delta n_{HI} (\vec{r_{g}}) \right ],
\end{aligned}
\end{equation}
where $\vec{r_{g}} = (\vec{R{g}}, z_{g})$ is a vector centred at the galactic centre with $\vec{R_{g}}$ at the plane of the galaxy. We assume a constant inclination angle $i$ between the line of sight direction and $z_{g}$ for simplicity. This gives a rotation mapping between the coordinates $\vec{r}$ and $\vec{r}_{g}$. Turbulent {\rm H~{\sc i}~} number density fluctuation $\delta_{HI}(\vec{r}_{g})$ is expected to have a power law spectrum and we quantify it with a slope $\alpha$. We assume that the tangential rotation velocity $v_{t}(R_{g})$ follow Brandt profile \citep{1960ApJ...131..293B} of the form
\begin{equation}
v_{t}(R_{g}) \ =\ \frac{ \frac{R_{g} } {R_{v}}  v_{0} }{ \left [  \frac{1}{3} + \frac{2}{3} \frac{R_{g}}{R_{v}} \right ]^{3/2} },
\end{equation}
and use this to generate the line of sight component of the systematic velocity of the {\rm H~{\sc i}~} gas across the galactic disk
\begin{equation}
v_{z}(\vec{r}) \ =\ \frac{v_{t}(R_{g}) sin(i) x }{R_{g}} +  f_{v}v^{T}_{z}(\vec{r}).
\end{equation}
We shall discuss the other parameters in above expressions shortly.

We consider $512^{3}$ cubic grids in our simulation volume with each cube representing an individual {\rm H~{\sc i}~} cloud. To accommodate the simulated galaxy's disk inside the simulation volume and adopt a generic rotation curve, we considered both $R_{0} = R_{v} = 100$ grid units. Assuming {\rm H~{\sc i}~} extent of a typical  spiral galaxy is $\sim 50$ kpc these choice gives us  a spatial resolution of $\sim 200$ pc.  Simulating a thin disk requires relatively higher computational resources. On the other hand inclination angle affects both the thick and thin disks in similar way, we opted for simulating a thick disk here and  adopt a value of $z_{0} = 100$ grid units.  We show these radial and vertical profiles normalised to their central value in Figure~(1) along with the tangential velocity profile normalised to $v_{0}$.

\begin{figure}
\begin{center}
\epsfig{file=fig2.eps,width=2.8in}
\end{center}
\caption{${\cal P}_{1}^{T}(\vec{U})$ estimated from the simulation is shown with dot-dashed line. This uses the parameters given in the bottom left corner of the image ($\beta = -2.0$). A model spectrum ${\cal P}_{\chi_{mod}}^{T}(\vec{U})$ with $\gamma = -2.5$ is also shown (dashed line) and it does not match with ${\cal P}_{\chi}^{T}(\vec{U})$. We also show the corresponding column density power spectrum  (solid black line). All three spectra are normalised at $U=1$.}
\label{fig:PS}
\end{figure}

\begin{figure*}
\begin{center}
\epsfig{file=fig3a.eps,width=2.4in}
\epsfig{file=fig3b.eps,width=2.4in}\\
\epsfig{file=fig3c.eps,width=2.4in}
\epsfig{file=fig3d.eps,width=2.4in}
\end{center}
\caption{Best fit values of $\beta$, i.e, $\hat{\beta}$  found using the estimator to the simulated data is plotted against the input values of $\beta$ with $3\sigma$ errors for four different sets of simulations. The solid line represents the ideal case of $\hat{\beta} = \beta$. Clearly at high inclination angles the estimated values of $\beta$ has relatively higher uncertainty.}
\label{fig:NHI}
\end{figure*}

 \citet{2013MNRAS.436L..49D} have reported that the amplitude of the {\rm H~{\sc i}~} column density fluctuation of the spiral galaxies in the  THINGS sample\footnote{THINGS: The {\rm H~{\sc i}~} Nearby Galaxy Survey \citep{2008AJ....136.2563W}} are approximately $1/10^{th}$ of the mean column density. Following this, we choose $f_{HI} = 0.1$. Line of sight velocity dispersion of {\rm H~{\sc i}~} for the galaxies in the THINGS sample varies in the range  $\sim 5 $ to $20$ km sec$^{-1}$ \citep{2009AJ....137.4424T}, while the flattening value of the rotation curve are found to  be in the range $\sim 100$ to $200$ km sec$^{-1}$ \citep{2008AJ....136.2648D}. Since this velocity dispersion is mostly because of turbulence, we adopt $f_{v} = v_{0}/10$ in the simulation,  the actual value of $v_{0}$ is unimportant here. \citet{2011ApJ...737L..33K} found that in our Galaxy, line of sights with column density $> 2 \times 10^{20}$ atoms cm$^{2}$ have temperature of $240$ K or lower, whereas the mean column density of the spiral galaxies in the THINGS sample is $> 3 \times 10^{20}$ atoms cm$^{2}$ \citep{2008AJ....136.2563W,2013MNRAS.436L..49D}. This suggest existence of large amount of cold gas in the spiral galaxies. In view of this, thermal velocity dispersion $\sigma$ is considered to be $v_{0}/100$. For a $ v_{0} = 200 $ km sec$^{-1}$, this corresponds to a gas  temperature $500$ K, i.e from cold neutral medium.

 {\rm H~{\sc i}~} column density power spectra of 18 spiral galaxies in the THINGS sample  are found to follow power laws \citep{2013NewA...19...89D}  for length scales more than the scale height of the disk and the slope of the power law for most of the galaxies were found to be in the range $-1.5$ to $-1.8$.   To see the effect of $\alpha$ on the efficacy of our estimator we perform the simulation for two different  $\alpha = (-2.5, -1.5)$.  We consider a nearly face on disk with $i = 10^{\circ}$ and a more edge on disk with $i = 45 ^{\circ}$.

With these choices,  we run four sets of simulations as {\bf  (a)} $( i = 10^{\circ}, \alpha = -1.5)$, {\bf (b)} $( i = 10^{\circ}, \alpha = -2.5)$, {\bf (c)} $( i = 45^{\circ}, \alpha = -1.5)$ and {\bf (d)} $( i = 45^{\circ}, \alpha = -2.5)$. Turbulent velocity component $v^{T}(\vec{r})$ is also taken to follow power law power spectrum \citep{2009ApJ...692..364F},  with a slope $\beta$.  \citet{2000ApJ...541..701G}  pointed out that it is likely that the ISM density fluctuations passively follow the velocity fluctuations and hence they would have similar power spectra. Observations  show that the column density fluctuation power spectra has a slope of $\sim -2.5$ \citep{1993MNRAS.262..327G} in our Galaxy and  $\sim -1.5$  \citep{2013MNRAS.436L..49D} in case of the external galaxies. We simulate the {\rm H~{\sc i}~} emission for $\beta$ varying in the range $-0.5$ to $-3.5$ in steps of $0.1$ and generate {\rm H~{\sc i}~} specific intensity data cubes following eqn.~(1) with a velocity resolution of $\sigma/2$. This amounts in calculating the line shape function  for at least six points within it's support, necessary and sufficient to avoid the shot noise bias discussed in \citet{2009ApJ...693.1074C}. We convert these specific intensity maps to corresponding visibilities and use them as the observed visibilities $V({\vec{ U}}, v)$ for further analysis. In this letter  we do not consider any measurement noise. Further, visibilities are  estimated at grids in $\vec{U}$ space and hence the baseline coverage is assumed to be complete.

\section{Estimating turbulent velocity fluctuation from the simulated data}

We proceed to estimate the power spectrum of turbulent velocity fluctuations through the following steps.
\begin{itemize}
\item
We estimate the visibility moments ${\cal V}_{0}(\vec{U})$ and ${\cal V}_{1}(\vec{U})$ from the visibilities $V({\vec{ U}}, v)$.
\item
We reconstruct the image data cube from visibilities and estimate the rotation curve using a standard tilted ring fit to the data. Since, in the simulation we have kept the inclination angle constant, we also keep it to be constant with radius while performing the tilted ring fit. For a real data, we would need to consider the uncertainties in the estimate of the rotation curve (and also inclination and position angles). At the present case these uncertainties are not considered. Using the rotation curve estimate we construct $\tilde{v}^{\Omega}_{z}({\vec{ U}})$.
\item
We bin the data in annular bins and estimate the azimuthally averaged ${\cal P}_{0}(U)$ from the visibility moment ${\cal V}_{0}(\vec{U})$ and ${\cal P}_{\chi}(\vec{U})$ from ${\cal V}_{1}(\vec{U})$ and  $\tilde{v}^{\Omega}_{z}({\vec{ U}})$ in the same bins.
\item
Assuming the power spectrum of the turbulent velocity fluctuations to follow power law, we consider ${\cal P}_{v_{mod}}^{T}(\vec{U}) \sim U^{\gamma}$ to model its power spectrum. We convolve ${\cal P}_{v_{mod}}^{T}(\vec{U})$ with ${\cal P}_{0}(\vec{U})$ to make the azimuthally averaged model ${\cal P}_{\chi_{mod}}^{T}(\vec{U})$. It is important to realise here that as discussed in Section~(3) our model of ${\cal P}_{\chi_{mod}}^{T}(\vec{U})$ derived this way does not have the proper amplitude. We shall restrict ourself in estimating the slope of the turbulent velocity power spectra here. The solid line in Figure~(2) shows the column density power spectra for the case (b). The dotted line corresponds to ${\cal P}_{\chi}^{T}(\vec{U})$ for $\beta = -2.0$ while the dot-dash line shows the model ${\cal P}_{\chi_{mod}}^{T}(\vec{U})$ with $\gamma=-2.5$. As expected, the model and the estimated value are in mismatch here.

\item
We vary $\gamma$ and compare ${\cal P}_{\chi_{mod}}^{T}(\vec{U})$ with   ${\cal P}_{\chi}^{T}(\vec{U})$ by a least square estimator with likelihood function of the form $\exp(-{\cal L})$, where ${\cal L}$ is defined as
\begin{equation}
{\cal L}(\gamma) = \frac{L (\gamma)}{L_{min}},\ \ L(\gamma) =\ \sum_{i = 1}^{N} \left [  \mathcal{P}_{\chi}^{T} \left (\frac{U}{U_{0}} \right )[i]  - \mathcal{P}_{\chi_{ mod}}^{T} \left ( \frac{U[i]}{U_{0}} \right ) \right ]^{2}.
\end{equation}
Here the index `i' refer to a particular annular  bin and $N$ is the total number of bins.  Note that both the power spectra are normalised at a given baseline $U_{0}$ and hence this estimator is not sensitive to the amplitude of the velocity fluctuation spectra.  We find the best fit value of $\gamma$, i.e, $\hat{\beta}$ and also error in estimating this value.

\end{itemize}

Figure~(3) summarises the  efficacy of the model dependent maximum likelihood estimation of the slope of the turbulent velocity fluctuation power spectrum from the four sets of simulations. Let us concentrate on the panel (a) first. Here $\hat{\beta}$ is plotted against the different values of input $\beta$ for $\alpha = -1.5$ and $i = 10^{\circ}$ with $3\sigma$ error bars. The solid line corresponds to an exact match between the input and estimated values of $\beta$.

Clearly, for a lower inclination angle, i.e $i = 10^{\circ}$, the estimator does a better job in accessing the real value of $\beta$. For the higher inclination angle,  i.e $i = 45^{\circ}$, for all the values of $\beta$ the best estimate has systematically higher uncertaininity. This is expected as for high inclination angle our assumption of independency of $v_{z}$ on line of sight coordinate is poorly valid. We also see for $\beta$ values $< -3.0$ the measurement uncertainty is considerably high.

\section{Discussions and Conclusions}
In this letter we define visibility moments and  an estimator of turbulent velocity fluctuations based on the visibility moments.Using  numerical simulations we find that  for lower inclination angles and  turbulent velocity fluctuation power law slope  $>-3.0$, the visibility moment estimator can  be used.  Though we have used a model dependent approach here, in principle it is possible to deconvolve the column density power spectra from ${\cal P}_{\chi}^{T}(\vec{U})$ and measure the  velocity power spectrum of the line of sight velocity fluctuations directly. Moreover, applying  this method to  the   thin disk case, we can recover the turbulent velocity power spectrum completely.  This would be particularly useful to access the energy involved in the large scale fluctuations seen in external galaxies.

We have considered several simplifying assumptions. Firstly, we  ignored the uncertainty in reconstructing the  rotation velocity of the galaxy from the data, observational noise and incomplete baseline coverage. Each of these are expected to provide additional uncertainty for a real data and need to be looked upon carefully. Particularly, recovering the turbulent velocity power spectra depends highly on the accuracy of the rotational velocity estimate. For galaxies with warp disks, estimating rotational velocity is tricky and this method may not be useful there. Random fluctuations that may arise because of different systematic effects in the galactic rotation, however, are not expected to have scale invariant nature.  An unwanted mixer of these fluctuations with the fluctuations generated because of the underlying turbulence would destroy the power law nature of the spectrum.  Hence, estimate of the velocity spectra itself can be used as a self diagnostic test.

The approach described here to estimate the line of sight turbulent velocity fluctautions from the external galaxies, apparently follow the same philosophy from the method of statistics of centroids of velocities \citep{2005ApJ...631..320E}, where  the unnormalised first moment of the specific intensity function is used to estimate  turbulent velocity statistics. Main difference in our approach is that we use the visibility function directly. Hence, it  does not have any unwanted noise bias like that arise because of  image reconstruction from incompletely sampled visibilities.

Scaling relations of compressible and highly magnetized turbulence are studied in \citet{2007ApJ...666L..69K}. \citet{2005ApJ...631..320E} has also discussed incompatibility of the velocity centroid method for high Mach number turbulence. Though the $10$ kpc scales coherent structures seen at external galaxies are most likely be the result of ISM turbulence, it is neither clear if the neutral medium ({\rm H~{\sc i}~}) at these scales would have an MHD coupling or they are highly supersonic. In either case, the applicability of the method described here needs to be investigated. Another basic assumption in developing this method has been that the {\rm H~{\sc i}~} is optically thin. The interpretation of the intensity spectra would change in case the optical depth of the {\rm H~{\sc i}~} is sufficiently higher \citep{2004ApJ...616..943L}.

Nevertheless, here we discussed an  effective way of measuring the statistics of the turbulent velocity fluctuations from external galaxies at large scales.
We plan to use this estimator to the spiral galaxies with  lower inclination angle from the THINGS sample in near future.

\section*{Acknowledgement}
PD would like to acknowledge useful discussion with Nirupam Roy. This work is supported by the DST INSPIRE Faculty Fellowship award  [IFA-13 PH 54] and performed at Indian Institute of Science Education and Research, Bhopal.
PD is thankful to Swagat Saurav Mishra for reading the earlier version of the draft and providing valuable comments.

\bibliographystyle{mn2e}


\section{Introduction}
\label{sec:intro}

Thanks to the successful running of the Large Hadron Collider~(LHC), particle physics has entered a new era, which is featured by the discovery of a new boson with a mass close to $125$~GeV~\cite{ATLASh,CMSh}. Its measured properties are so far in good agreement with those of the Standard Model~(SM) Higgs~\cite{ATLASh2,CMSh2,CDFD0h2}, suggesting that the electro-weak symmetry breaking~(EWSB) is probably realized via the Higgs mechanism implemented through a single scalar doublet. It should be noted, however, that the EWSB is not necessarily induced by just a single scalar. Interestingly, many new physics~(NP) scenarios are equipped with an extended scalar sector. The search for additional scalars is one of the important programs of the LHC experiments.

One of the extensions of SM scalar sector is the so-called two-Higgs-doublet model~(2HDM)~\cite{Lee2HDM}, in which a second  scalar doublet is added to the SM field content. To avoid the experimental constraints on flavor-changing neutral-current~(FCNC) transitions, which are forbidden at tree level in the SM due to the GIM mechanism~\cite{GIM}, two different hypotheses, natural flavor conservation~(NFC)~\cite{NFC} and minimal flavor violation~(MFV)~\cite{MFV}, have been proposed~\footnote{The NFC and MFV hypotheses are not the only alternatives to avoid constraints from FCNCs; models with controlled FCNCs have also been addressed in the literature and shown to be compatible with the data~\cite{FCNC_BGL}.}. In the NFC hypothesis, depending on the $Z_2$ charge assignments on the scalar doublets and fermions, there exist four types of 2HDM~(type-I, II, X and Y)~\cite{Branco}. In the MFV hypothesis, to control the flavor-violating interactions, all the scalar Yukawa couplings are assumed to be composed of the SM ones $Y^U$ and $Y^D$. In the ``Higgs basis"~\cite{Davidson:2005cw}, in which only one doublet gets a nonzero vacuum expectation value~(VEV) and behaves as the SM one, the allowed $SU(3)_C\otimes SU(2)_L\otimes U(1)_Y$ representation of the second scalar doublet that couples to quarks via Yukawa interactions is fixed to be either $(1,2)_{1/2}$ or $(8,2)_{1/2}$~\cite{Wise}, which implies that the second scalar doublet can be either color-singlet or color-octet. For convenience, they are referred to as type-III and type-C models~\cite{Slavich}, respectively. Examples of the former include the aligned 2HDM~(A2HDM)~\cite{A2HDM} and the four types of 2HDM reviewed in Ref.~\cite{Branco}. The scalar spectrum of the latter contains, besides a CP-even and color-singlet Higgs boson~(the usual SM one), three color-octet particles, one CP-even, one CP-odd and one electrically charged~\cite{Wise}.

Although the scalar-mediated flavor-violating interactions are protected by the MFV hypothesis, the type-III and type-C models still present very interesting phenomena in some low-energy processes, especially due to the presence of a charged Higgs boson~\cite{Wise,Slavich,Yuan}. In this paper, we shall study the $B_s^0-\bar{B}_s^0$ mixing within these two models and pursue possible differences between their effects. Since the charged Higgs contributes to the process at the same order as does the $W$ boson in the SM, the NP effects might be significant.

It is known that the $B_s^0-\bar{B}_s^0$ mixing, which is governed by a Schr\"odinger equation and induced by the $b\to s$ transition, plays an important role in accurate tests of the SM and indirect searches for NP. In terms of the off-diagonal elements of the mass and decay matrices, $M_{12}^s$ and $\Gamma_{12}^s$, the mass and width differences between the two mass eigenstates $|B_{H}\rangle$ and $|B_{L}\rangle$ are defined, respectively, by
\begin{align}\label{Dm-Dga-def}
\Delta M_s \equiv M_{H}-M_{L}=2 |M_{12}^s|\,, \quad
\Delta \Gamma_s \equiv \Gamma_{L}-\Gamma_{H}=2 |\Gamma_{12}^s|\cos\phi_s\,,
\end{align}
where $\phi_s \equiv \arg(-M_{12}^s/\Gamma_{12}^s)$ is the CP-violating phase. Such two observables have been measured precisely, with the averaged values given, respectively, by~\cite{HFAG}
\begin{align}
 \Delta M_s^{exp.}=17.757\pm0.021~{\rm ps^{-1}}\,,\quad
 \Delta \Gamma_s^{exp.}=0.081\pm 0.006~{\rm ps^{-1}}\,,
\end{align}
which are in good agreement with the recent SM predictions, $\Delta M_s^{SM}=(17.3\pm2.6)\,{\rm ps^{-1}}$ and $\Delta \Gamma_s^{SM}=(0.087\pm0.021)\,{\rm ps^{-1}}$~\cite{Lenz3}.

There are another two interesting observables related to $B_s^0-\bar{B}_s^0$ mixing, the flavor-specific CP asymmetry $a_{sl}^s$ and the CP-violating phase $\phi_s^{c\bar{c}s}$~\footnote{The phase $\phi_s^{c\bar{c}s}$ appears in tree-dominated $b\to c\bar{c}s$ $B_s$ decays, such as $B_s\to J/\psi\phi$, and is generally different from $\phi_s$ unless the terms proportional to $V_{cb}V_{cs}^{\ast}V_{ub}V_{us}^{\ast}$ and $(V_{ub}V_{us}^{\ast})^2$ in $\Gamma^s_{12}$ are neglected~\cite{Lenz1}.}, which are defined, respectively, by
\begin{equation}
a_{sl}^{s}={\rm Im}\frac{\Gamma^s_{12}}{M_{12}^s}=\frac{\Delta M_s}{\Delta \Gamma_s}\tan\phi_s\,, \qquad
\phi_s^{c\bar{c}s}={\rm arg}(M_{12}^s)\,.
\end{equation}
For $\phi_s^{c\bar{c}s}$, the SM prediction $\sim -0.036$~\cite{Lenz3} agrees with the experimental data $-0.015\pm0.035$~\cite{HFAG} within $1\sigma$ error bar. For $a_{sl}^{s}$, on the other hand, the measurement $(-0.75\pm0.41)\%$~\cite{HFAG} is significantly different from the SM estimation $\sim {\cal O}(10^{-5})$~\cite{Lenz3}, even through they are in agreement with each other at $1.5\sigma$ level due to the large experimental error bars. Under the constraints of the above four observables, the allowed NP spaces could possibly be strictly reduced. So, in this paper, we shall evaluate the effects of charged Higgs in type-III and type-C 2HDMs on $B_s^0-\bar{B}_s^0$ mixing, and pursue possible differences between these two models.

Our paper is organized as follows. In Sec.~2, after a brief review of the 2HDMs under the MFV hypothesis, we perform a complete one-loop computation of the electro-weak corrections to the amplitudes of $B_s^0-\bar{B}_s^0$ mixing within the two models. In Sec.~3, the numerical results and discussions are presented in detail. Finally, our conclusions are made in Sec.~4. Explicit expressions for the loop functions appearing in $B_s^0-\bar{B}_s^0$ mixing are collected in the appendix.

\section{Theoretical Framework}
\label{sec:theor}

\subsection{Brief review of the 2HDMs under the MFV hypothesis}
\label{subsec:2HDM}

Firstly, for convenience and consistence, we shall give a brief review of the 2HDMs under the MFV hypothesis. In the ``Higgs basis'', the Yukawa interactions of the two Higgs fields $\Phi_1$ and $\Phi_2$ with quarks are given by~\cite{Wise,Slavich}
\begin{align}\label{L:Yukawa}
-\mathcal{L}_Y=\bar q_L^0 \tilde \Phi_1 Y^U u_R^0
    +\bar q_L^0 \Phi_1 Y^D d_R^0
    +\bar q_L^0 \tilde\Phi_2^{(a)}T_R^{(a)}\bar Y^U u_R^0
    +\bar q_L^0       \Phi_2^{(a)}T_R^{(a)}\bar Y^D d_R^0
    +\rm{h.c.}\,,
\end{align}
where $q_L^0$, $u_R^0$ and $d_R^0$ are the quark fields in the interaction basis; $T_R^{(a)}$ is the $SU(3)_C$ generator and determines the color nature~(color-singlet or color-octet) of the second scalar doublet; $Y^{U,D}$ and $\bar Y^{U,D}$ denote the Yukawa couplings and are generally complex $3 \times 3$ matrices in the quark flavor space. According to the MFV hypothesis, the transformation properties of the Yukawa coupling matrices $Y^{U,D}$ and $\bar Y^{U,D}$ under the quark flavor symmetry group $SU(3)_{Q_L}\otimes SU(3)_{U_R} \otimes SU(3)_{U_D}$ are required to be the same. This can be achieved by requiring $\bar Y^{U,D}$ to be composed of pairs of the matrices $Y^{U,D}$~\cite{Slavich}
\begin{align} \label{eq:Y}
  \bar Y^U&=A_u^\ast (1+\epsilon_u^\ast Y^UY^{U\dagger}+\dotsc)Y^U\,, \nonumber\\[0.2cm]
  \bar Y^D&=A_d(1+\epsilon_d Y^UY^{U\dagger}+\dotsc)Y^D\,,
\end{align}
where the ellipses denote trivial terms involving higher power of $Y^U Y^{U\dagger}$ and powers of $Y^D Y^{D\dagger}$.

After applying the SM unitary transformations to rotate the fermionic fields from the interaction to the mass-eigenstate basis, one can finally obtain the Yukawa interactions of charged Higgs bosons with quarks in the mass-eigenstate basis~\cite{Slavich,Yuan}
\begin{align}\label{eq:LCH}
  \mathcal L_{H^+}=\frac{g}{\sqrt{2}m_W}\,\sum_{i,j=1}^3\,\bar u_iT_R^{(a)}(A_u^i m_{u_i}P_L-A_d^im_{d_j}P_R)V_{ij}d_jH_{(a)}^+ +\rm{h.c.}\,,
\end{align}
where $g$ is the $SU(2)_L$ coupling constant, $i\,,j$ the fermionic generation indices, and $m_{u,d}$ the quark masses; $V$ denotes the Cabibbo--Kobayashi--Maskawa~(CKM) matrix~\cite{Cabibbo:1963yz,Kobayashi:1973fv}, and $P_{R,L}=\frac{1\pm \gamma_5}{2}$ are the right- and left-handed chirality projectors. The couplings $A_{u,d}^i$ are generally family-dependent and read
\begin{align}
  A_{u,d}^i=A_{u,d}\left(1+\epsilon_{u,d}\frac{m_t^2} {v^2}\delta_{i3}\right),
\end{align}
where $v=\langle \Phi_1^0\rangle=174~{\rm GeV}$. Since only the couplings of charged Higgs bosons to the top quark are involved for $B_s^0-\bar{B}_s^0$ mixing, we shall drop the family index in $A_{u,d}^i$ from now on.

Following the notation used in Ref.~\cite{Slavich}, we shall denote the model with the second scalar doublet being color-singlet and the one with the second scalar doublet color-octet as the type-III and the type-C model, respectively, both of which satisfy the principle of MFV. Their explicit contributions to the $B_s^0-\bar{B}_s^0$ mixing
will be presented in the next subsection.

\subsection{$\boldsymbol{B_s^0-\bar{B}_s^0}$ mixing within the SM and the 2HDMs with MFV}
\label{subsec:wcs}

\begin{figure}[t]
\begin{center}
\subfigure[]{\includegraphics[width=3.6cm]{Topo1.pdf}} \quad
\subfigure[]{\includegraphics[width=3.6cm]{Topo2.pdf}}
\quad
\subfigure[]{\includegraphics[width=3.6cm]{Topo3.pdf}} \quad
\subfigure[]{\includegraphics[width=3.6cm]{Topo4.pdf}}
\caption{\label{topology} \small Box diagrams relevant to the $B_s^0-\bar{B}_s^0$ mixing in the unitary gauge, both within the SM~(the first one) and in the 2HDMs with MFV~(the last three ones). We have also taken into account the crossed diagrams, which are related to the original ones by interchanging the external lines.}
\end{center}
\end{figure}

Within the 2HDMs with MFV, the $B_s^0-\bar{B}_s^0$ mixing occurs through the box diagrams shown in Figs.~\ref{topology}(b)-(d), which are obtained from the SM one~(Fig.~\ref{topology}(a)) with the $W^{\pm}$ propagator(s) replaced by the charged-Higgs $H^{\pm}$ one(s). After calculating these one-loop box diagrams and applying the standard procedure of matching~\cite{Buras:1990fn,Urban:1997gw}, one can obtain the 2HDM corrections to the $B_s^0-\bar{B}_s^0$ mixing. Together with the SM contribution, the resulting effective weak Hamiltonian responsible for $B_s^0-\bar{B}_s^0$ mixing can be written as
\begin{equation}\label{eq:heff}
  \mathcal{H}_{eff}^{full}=\frac{G_F^2}{16\pi^2}\,m_W^2\,\left(V_{tb}V_{ts}^*\right)^2\,\left[\,C^{VLL}(\mu) O^{VLL}+C^{SRR}(\mu)O^{SRR}+C^{TRR}(\mu)O^{TRR}\,\right]+\text{h.c.}\,,
\end{equation}
where $G_F$ is the Fermi coupling constant, and $C_{i}(\mu)$ the scale-dependent Wilson coefficients of the four-quark operators $O_{i}$, which are defined, respectively, as
\begin{align}\label{oper}
  &O^{VLL}=\bar{s}^\alpha\gamma_\mu(1-\gamma_5)b^\alpha\,\bar{s}^\beta\gamma^\mu(1-\gamma_5)b^\beta\,,\nonumber\\[0.2cm]
  &O^{SRR}=\bar{s}^\alpha(1+\gamma_5)b^\alpha\,\bar{s}^\beta(1+\gamma_5)b^\beta\,,\nonumber\\[0.2cm]
  &O^{TRR}=\bar{s}^\alpha\sigma_{\mu\nu}(1+\gamma_5)b^\alpha\,\bar{s}^\beta\sigma^{\mu\nu}(1+\gamma_5)b^\beta\,,
\end{align}
with $\alpha,\beta$ being the color indices and $\sigma_{\mu\nu}=\frac{1}{2}[\gamma_{\mu},\gamma_{\nu}]$. In addition, the hadronic matrix elements of the four-quark operators can be parameterized as~\cite{BurasDF2}
\begin{align}
  &\langle B_s^0 | O^{VLL} |\bar{B}_s^0 \rangle=\frac{8}{3}m_{B_s}^2f_{B_s}^2B_1(\mu_b)\,,\nonumber\\[0.2cm]
  &\langle B_s^0 |  O^{SRR} |\bar{B}_s^0 \rangle=-\frac{5}{3}\Big(\frac{m_{B_s}}{\bar{m}_b(\mu_b)+\bar{m}_s(\mu_b)}\Big)^2m_{B_s}^2f_{B_s}^2B_2(\mu_b)\,,\nonumber\\[0.2cm]
  &\langle B_s^0 |  O^{TRR} |\bar{B}_s^0 \rangle=\frac{4}{3}\Big(\frac{m_{B_s}}{\bar{m}_b(\mu_b)+\bar{m}_s(\mu_b)}\Big)^2m_{B_s}^2f_{B_s}^2\left[-5B_2(\mu_b)+2B_3(\mu_b)\right]\,,
\end{align}
where $f_{B_s}$ is the $B_s$-meson decay constant, and $B_i(\mu_b)$ are the non-perturbative bag parameters calculated on the lattice at a characteristic scale $\mu_b\sim\mathcal{O}(m_b)$~\cite{Carrasco:2013zta}.

The Wilson coefficients $C_i(\mu)$ in Eq.~(\ref{eq:heff}) consist of both the SM and 2HDM contributions and the results at the initial scale $\mu_W\sim\mathcal{O}(m_W,m_t,m_{H^{\pm}})$ can be written as
\begin{align}
  &C^{VLL}(\mu_W)=C^{VLL}_{SM}(\mu_W)+C^{VLL}_{2HDM}(\mu_W)\,,\nonumber\\[0.2cm]
  &C^{SRR}(\mu_W)=C^{SRR}_{SM}(\mu_W)+C^{SRR}_{2HDM}(\mu_W)\,,\nonumber\\[0.2cm]
  &C^{TRR}(\mu_W)=C^{TRR}_{SM}(\mu_W)+C^{TRR}_{2HDM}(\mu_W)\,.
\end{align}
Within the SM, the explicit expressions of the Wilson coefficients are computed from the box diagram shown in Fig.~\ref{topology}(a), accompanied by perturbative QCD corrections up to the desired order, details of which could be found, for example, in Refs.~\cite{Buras:1990fn,Urban:1997gw,Buchalla}. Including the next-to-leading order (NLO) QCD corrections, the SM contribution is given by~\cite{Buras:1990fn,Urban:1997gw,Buchalla}
\begin{equation}
  C^{VLL}_{SM}(\mu_W)=S_0(x_t)+\frac{\alpha_s(\mu_W)}{4\pi}\,\left[S_1(x_t)+F(\mu_W)S_0(x_t)+B_tS_0(x_t)\right]\,,
\end{equation}
where $x_t=\frac{\bar{m}_t^2(\mu_W)}{m^2_W}$, and the leading order (LO) coefficient $S_0(x_t)$ is the known Inami-Lim function~\cite{Inami:1980fz}. Generally, the SM also contributes to $C^{SRR}$ and $C^{TRR}$ and, in the absence of QCD corrections, we get
\begin{align}
  C^{SRR}_{SM}(\mu_W)\,=\,&\frac{x_b}{6}\bigg[\frac{x_t^2(5-22x_t+5x_t^2)}{3(1-x_t)^4}+\frac{x_t^2(1-3x_t-3x_t^2+x_t^3)\ln{x_t}}{(1-x_t)^5}\bigg]\,,\\[0.2cm]
  C^{TRR}_{SM}(\mu_W)\,=\,&\frac{x_b}{6}\bigg[-\frac{5-15x_t+8x_t^2-15x_t^3+5x_t^4}{3(1-x_t)^4}+\frac{(1-5x_t+9x_t^2-x_t^3)\ln{x_t}}{(1-x_t)^5}\bigg]\,,
\end{align}
with $x_b=\frac{\bar{m}_b^2(\mu_W)}{m^2_W}$. It is obvious that both $C^{SRR}_{SM}(\mu_W)$ and $C^{TRR}_{SM}(\mu_W)$ are suppressed by the factor $x_b$ and are, therefore, usually neglected in the literature~\cite{Buras:1990fn}.

The charged-Higgs contributions to the Wilson coefficients are computed from the last three box diagrams shown in Fig.~\ref{topology}, and depend on the two Yukawa coupling parameters $A_u$ and $A_d$, as well as the charged-Higgs mass $m_{H^{\pm}}$~\cite{Urban:1997gw,Buras:2001mb,Jung:2010ik,Manohar:2006ga,Cheng:2015lsa}. For the most general values of these parameters, especially when $A_d/A_u\simeq m_t/m_b$, each term in Eq.~(\ref{eq:LCH}) can give a comparable contribution and should be, therefore, taken into account simultaneously. Explicitly, for the color-singlet charged-Higgs contributions~(type-III model), we get
\begin{align} \label{eq:Ci_III}
  C^{VLL}_{III}(\mu_W)=&A_uA_u^*\,f_1(x_t,x_h)+A_u^2A_u^{*2}\,f_2(x_t,x_h)\,,\\[0.2cm]

  C^{SRR}_{III}(\mu_W)=&-x_b\Big[A_uA_u^*\,f_3(x_t,x_h)+A_dA_u^*\,f_4(x_t,x_h)+A_d^2A_u^{*2}\,f_5(x_t,x_h)\,\nonumber\\[0.1cm]
                       &\hspace{1.1cm}+A_u^2A_u^{*2}\,f_6(x_t,x_h)+A_dA_uA_u^{*2}\,f_7(x_t,x_h)\Big]\,,\\[0.2cm]

  C^{TRR}_{III}(\mu_W)=&0\,,
\end{align}
where $x_h=m^2_{H^{\pm}}/m^2_W$, and the explicit expressions for $f_i(x_t,x_h)$ are collected in the appendix. For the color-octet charged-Higgs contributions~(type-C model), on the other hand, we get
\begin{align} \label{eq:Ci_C}
  C^{VLL}_{C}(\mu_W)=&\frac{1}{3}\,A_uA_u^*\,f_1(x_t,x_h)+\frac{11}{18}\,A_u^2A_u^{*2}\,f_2(x_t,x_h)\,,\\[0.2cm]

  C^{SRR}_{C}(\mu_W)=&-x_b\bigg[-\frac{5}{12}\,A_uA_u^*\,f_3(x_t,x_h)-\frac{5}{12}\,A_dA_u^*\,f_4(x_t,x_h)
                     -\frac{19}{72}\,A_d^2A_u^{*2}\,f_5(x_t,x_h)\,\nonumber\\[0.1cm]
                     &\hspace{1.2cm}-\frac{19}{72}\,A_u^2A_u^{*2}\,f_6(x_t,x_h) -\frac{19}{72}\,A_dA_uA_u^{*2}\,f_7(x_t,x_h)\bigg]\,,\\[0.2cm]

  C^{TRR}_{C}(\mu_W)=&-x_b\bigg[\frac{1}{16}\,A_uA_u^*\,f_3(x_t,x_h)+\frac{1}{16}\,A_dA_u^*\,f_4(x_t,x_h)
                     +\frac{7}{96}\,A_d^2A_u^{*2}\,f_5(x_t,x_h)\,\nonumber\\[0.1cm]
                     &\hspace{1.2cm} +\frac{7}{96}\,A_u^2A_u^{*2}\,f_6(x_t,x_h) +\frac{7}{96}\,A_dA_uA_u^{*2}\,f_7(x_t,x_h)\bigg]\,.
\end{align}
It is noted that the Wilson coefficient $C^{TRR}_{C}(\mu_W)$ is now nonzero in the type-C model. To check the gauge independence of our results, we have performed the calculation both in the Feynman and in the unitary gauge. For $C^{VLL}$, our results agree with the ones presented in Refs.~\cite{Urban:1997gw,Buras:2001mb,Manohar:2006ga}~\footnote{There are two typos in Eq.~(26) of Ref.~\cite{Manohar:2006ga}: a global factor $2$ should be added to the term proportional to $|\eta_U|^2$ and $1/2$ to the term proportional to $|\eta_U|^4$.}. For $C^{SRR}$ and $C^{TRR}$, on the other hand, in order to get a gauge-independent result, the external momenta of the heavy quarks inside the mesons should be taken into account, and the heavy-quark masses should be kept up to the second orde; our results for these two coefficients differ from the ones presented in Refs.~\cite{Urban:1997gw,Buras:2001mb,Manohar:2006ga}.

The QCD renormalization group~(RG) evolution of these Wilson coefficients from the matching scale $\mu_W$ down to the lower scale $\mu_b$ has been calculated in Ref.~\cite{BurasDF2}. One can then obtain the corresponding Wilson coefficients at the scale $\mu_b$ through~\cite{BurasDF2}
\begin{align}
  C^{VLL}(\mu_b)&=\left[\eta(\mu_b)\right]_{VLL}C^{VLL}(\mu_W)\,, \\[0.2cm]
  \left(
      \begin{aligned}
         &C^{SRR}(\mu_b)\\[0.1cm]
         &C^{TRR}(\mu_b)
      \end{aligned}
   \right)&=
   \left(
      \begin{aligned}
         &\left[\eta_{11}(\mu_b)\right]_{SRR}\,\, &\left[\eta_{12}(\mu_b)\right]_{SRR}\\[0.1cm]
         &\left[\eta_{21}(\mu_b)\right]_{SRR}\,\,
         &\left[\eta_{22}(\mu_b)\right]_{SRR}
      \end{aligned}
   \right)
   \left(
     \begin{aligned}
         &C^{SRR}(\mu_W)\\[0.1cm]
         &C^{TRR}(\mu_W)
     \end{aligned}
   \right)\,,
\end{align}
where the explicit expressions of the evolution factors $\eta$ could be found in Ref.~\cite{BurasDF2}.

Equipped with the above information, the off-diagonal mass matrix element $M_{12}^s$ is given as
\begin{equation}
 M_{12}^s = \langle B_s^0|{\cal H}_{eff}^{full}|\bar{B}_s^0 \rangle
 ={\cal A}^{VLL}+{\cal A}^{SRR}+{\cal A}^{TRR}\,,
 \end{equation}
where ${\cal A}^{VLL}$, ${\cal A}^{SRR}$ and ${\cal A}^{TRR}$ denote the contributions induced by the three four-quark operators defined by Eq.~(\ref{oper}), respectively. Within the SM, the off-diagonal decay matrix element $\Gamma_{12}^s$ can be written as~\cite{Lenz,beneke}
\begin{eqnarray}
\Gamma_{12}^s(SM)=
    - \left[\, \lambda_t^2 \, \Gamma_{12}^{cc} \; + \;
        2 \, \lambda_t\,  \lambda_u \,
        \left( \Gamma_{12}^{cc} - \Gamma_{12}^{uc} \right) \; + \;
        \lambda_u^{2} \,
        \left( \Gamma_{12}^{cc} - 2 \Gamma_{12}^{uc} + \Gamma_{12}^{uu}\right)
       \right]\,, \label{ga12t}
\end{eqnarray}
with the CKM factors $\lambda_i=V_{ib}V_{is}^*$ for $i=u,c,t$. The explicit expressions for $\Gamma_{12}^{cc,uu,uc}$ could be found in Refs.~\cite{Lenz,beneke}. It should be noted that $\Gamma_{12}^s$ is dominated by the CKM-favored tree-level $b\to c\bar{c}s$ transition within the SM, and the NP effects are generally negligible~\cite{Lenz}. Hence $\Gamma_{12}^s=\Gamma_{12}^s(SM)$ holds as a good approximation, which will be assumed throughout this paper.

\section{Numerical results and discussions}
\label{sec:results}

We now proceed to present our numerical results and discussions. Values of the relevant input parameters used throughout this paper are summarized in Table~\ref{inputs}. Our SM predictions for the observables of $B_s^0-\bar{B}_s^0$ mixing are given in the third row of Table~\ref{pred}, in which the experimental data averaged by the HFAG~\cite{HFAG} are also listed in the second row for comparison. As mentioned already in the introduction section, there is no significant deviation between the SM predictions and the experimental data for the observables at the current level of precision, even through a bit disagreement appears for $a_{sl}^s$. Therefore, these observables are expected to put strong constraints on the parameter spaces of 2HDMs with MFV.

\begin{table}[t]
\begin{center}
 \caption{\label{inputs} \small Values of the relevant input parameters throughout this paper.}
 \small
 \vspace{0.2cm}
 \doublerulesep 0.8pt \tabcolsep 0.05in
 \begin{tabular}{l}
 \hline\hline
 $|V_{us}|=0.2253\pm0.0008$\,,\,
 $|V_{ub}|=0.00413\pm0.00049$\,,\,
 $|V_{cb}|=0.0411\pm0.0013$\,,\,
 $\gamma=(68.0^{+8.0}_{-8.5})^{\circ}$\, \cite{PDG14}
 \\ \hline
 $\bar{m}_{s}(2\,{\rm GeV})=95 \pm 5$~MeV\,,\quad
 $\bar{m}_{b}(\bar{m}_{b})=4.18 \pm 0.03$~GeV\,,\quad
 $m_t=173.21\pm 0.87$ GeV\, \cite{PDG14}\\
 $\frac{\bar{m}_s(\mu)}{\bar{m}_{u,d}(\mu)}=27.5 \pm 1.0$ \cite{PDG14}\,,\quad
 $m_{b}^{\rm pow}=4.8^{+0.0}_{-0.2}$ GeV\, \cite{Lenz}
 \\ \hline
 $f_{B_s}=228 \pm5\pm6$~MeV\,,\quad
 $f_{B_s}\sqrt{B_1}=211\pm5\pm6$~MeV\,,\quad
 $f_{B_s}\sqrt{B_2}=195\pm5\pm5$~MeV\,,\\
 $f_{B_s}\sqrt{B_3}=215\pm14\pm9$~MeV\, \cite{Carrasco:2013zta} \\
 \hline \hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[t]
\begin{center}
 \caption{\label{pred} \small Numerical results for $\Delta M_s[{\rm ps^{-1}}]$, $\Delta \Gamma_s[{\rm ps^{-1}}]$, $\phi^{c\bar{c}s}_s$ and $a_{sl}^s[\%]$ within the SM. The theoretical uncertainties are obtained by varying each input parameter listed in Table~\ref{inputs} within its respective allowed range and then adding the individual uncertainty in quadrature.}
 \small
 \vspace{0.2cm}
 \doublerulesep 0.8pt \tabcolsep 0.25in
 \begin{tabular}{lccccccccccc} \hline \hline
 &    $\Delta M_s$       &     $\phi^{c\bar{c}s}_s$    &$\Delta \Gamma_s$        & $a_{sl}^s[\%]$  \\ \hline
 Exp. & $17.757\pm0.021$ & $-0.015\pm0.035$ & $0.081\pm 0.006$ & $-0.75\pm0.41$ \\
 SM   & $17.228^{+1.731}_{-1.672}$ & $-0.043^{+0.006}_{-0.006}$ & $0.082^{+0.009}_{-0.013}$
      & $0.0026^{+0.0004}_{-0.0004}$\\
 \hline \hline
\end{tabular}
\end{center}
\end{table}

From the analytic expressions of the charged-Higgs contributions to the $B_s^0-\bar{B}_s^0$ mixing calculated in section~2.2, one can find that the model parameters relevant to our study include the two Yukawa coupling parameters $A_u$ and $A_d$, as well as the charged-Higgs mass $m_{H^{\pm}}$. For the case of complex couplings, one could equivalently choose $|A_u|$ and $A_dA_u^{*}=|A_dA_u^{*}|e^{-i\theta}$ as the independent parameters, with $\theta$ being the relative phase between the two Yukawa coupling parameters. For the parameter $|A_u|$, as detailed in Ref.~\cite{Slavich}, an upper bound can be obtained from the $Z\to b\bar{b}$ decay. The parameter $A_d$ is, however, much less constrained phenomenologically~\cite{Slavich,Yuan}. Concerning the charged-Higgs mass, we shall use the LEP lower bound $m_H^{\pm}>78.6$~GeV ($95\%$ CL)~\cite{LEPHpm}, which is obtained under the assumption that $H^{\pm}$ decays dominantly into fermions and does not refer to any specific Yukawa structure. Direct searches for charged-Higgs bosons are also performed by the Tevatron~\cite{Gutierrez:2010zz}, ATLAS~\cite{ATLASHpm} and CMS~\cite{CMSHpm} collaborations. However, most of the limits on $m_H^{\pm}$ depend strongly on the assumed Yukawa structure. In this paper, we shall generate randomly numerical points for the model parameters in the ranges
\begin{align}\label{eq:range}
|A_u|\in[0,3]\,,\qquad m_H^{\pm}\in[80,500]~\mathrm{GeV}\,,
\end{align}
whereas no severe constraints for $|A_dA_u^{*}|$~(or $|A_d|$) and $\theta$.

To compare the relative strength of the charged-Higgs contributions with respect to the SM one at the scale $\mu_b=m_b$, choosing $m_{H^{\pm}}=200$~GeV and the default values of the input parameters listed in Table~\ref{inputs}, we get
\begin{align}
\label{eq:AIII1}
\frac{{\cal A}^{VLL}}{(V_{tb}V_{ts}^*)^2}\times 10^{8}\simeq&3.73+2.00|A_u|^2 + 0.45|A_u|^4\,,\\[0.2cm]
\label{eq:AIII2}
\frac{{\cal A}^{SRR}}{(V_{tb}V_{ts}^*)^2}\times 10^{12}\simeq &-3.21|A_u|^2+0.69|A_u|^4+5.99A_dA_u^{*}-3.44|A_u|^2 A_dA_u^{*}+3.44A_d^2A_u^{*2}\,,\\[0.2cm]
\label{eq:AIII3}
\frac{{\cal A}^{TRR}}{(V_{tb}V_{ts}^*)^2}\times 10^{12}\simeq&0.03|A_u|^2 -0.01|A_u|^4 -0.05A_dA_u^{*} +0.03|A_u|^2 A_dA_u^{*}-0.03A_d^2A_u^{*2}\,,
  \end{align}
in the case of type-III model, and
\begin{align}
\label{eq:AC1}
\frac{{\cal A}^{VLL}}{(V_{tb}V_{ts}^*)^2}\times 10^{8}\simeq&3.73+0.67|A_u|^2 + 0.28|A_u|^4\,,\\[0.2cm]
\label{eq:AC2}
\frac{{\cal A}^{SRR}}{(V_{tb}V_{ts}^*)^2}\times 10^{12}\simeq &1.10|A_u|^2-0.12|A_u|^4-2.05A_dA_u^{*}+ 0.61|A_u|^2 A_dA_u^{*}-0.61A_d^2A_u^{*2}\,,\\[0.2cm]
\label{eq:AC3}
\frac{{\cal A}^{TRR}}{(V_{tb}V_{ts}^*)^2}\times 10^{12}\simeq&-0.15|A_u|^2 +0.04 |A_u|^4 + 0.28A_dA_u^{*} -0.18|A_u|^2 A_dA_u^{*}+0.18A_d^2A_u^{*2}\,,
  \end{align}
in the case of type-C model, respectively. The number $3.73$ in Eqs.~\eqref{eq:AIII1} and \eqref{eq:AC1} is the SM contribution, while the SM contributions to ${\cal A}^{SRR}$ and ${\cal A}^{TRR}$ are suppressed by the factor $x_b$, making them numerically smaller by about three orders than ${\cal A}^{VLL}_{SM}$ and hence negligible. From the above numerical results, we make the following observations:
\begin{enumerate}
\item[ (i) ] In both the type-III and type-C models, the charged-Higgs contributions to ${\cal A}^{VLL}$ (Eqs.~\eqref{eq:AIII1} and \eqref{eq:AC1}) depend only on the Yukawa coupling parameter $A_u$ via $|A_u|$, and hence are always constructive to the SM one. For a value $|A_u|\sim{\cal O}(1)$, the type-III contribution could be comparable with the SM one, while the type-C model provides a relatively smaller correction.

\item[ (ii) ] Comparing Eqs.~\eqref{eq:AIII2}-\eqref{eq:AIII3} with \eqref{eq:AIII1} (for the type-III model) and Eqs.~\eqref{eq:AC2}-\eqref{eq:AC3} with \eqref{eq:AC1} (for the type-C model), one can see that the NP contributions to ${\cal A}^{SRR}$ and ${\cal A}^{TRR}$ are much smaller than to ${\cal A}^{VLL}$, especially when $|A_d|\sim|A_u|$. This is because the Wilson coefficients $C^{SRR}(\mu_W)$ and $C^{TRR}(\mu_W)$ are always suppressed by the factor $x_b$ with respect to $C^{VLL}(\mu_W)$, both within the SM and in the 2HDMs with MFV.

\item[ (iii) ] In the case with large complex values of $A_dA_u^{*}$, however, the charged-Higgs contributions to ${\cal A}^{SRR}$ and ${\cal A}^{TRR}$ could provide a large imaginary part to the off-diagonal mass matrix element $M_{12}^s$, which may result in a significant correction to the CP-violating observables, such as $\phi_s$ and $\phi_s^{c\bar{c}s}$.

\item[ (iv) ] Different from the type-C model, the type-III contribution to ${\cal A}^{TRR}$ is induced only by the RG evolution effect, and is numerically much smaller. There are, however, cancelations between the charged-Higgs contributions to ${\cal A}^{SRR}$ and ${\cal A}^{TRR}$ in the type-C model.
\end{enumerate}
It is therefore expected that the current experimental data on $B_s^0-\bar{B}_s^0$ mixing could put some constraints on the model parameters and be used to distinguish between these two models.

To get the explicitly allowed parameter spaces, we perform the analysis with the following procedure: we scan the parameter spaces within the ranges specified by Eq.~\eqref{eq:range}, with the value of $m_{H^{\pm}}$ fixed at $100$, $250$ and $500~{\rm GeV}$, respectively. At each point in the parameter spaces, we evaluate the theoretical prediction for an observable, together with the corresponding theoretical uncertainty induced by the input parameters listed in Table~\ref{inputs}. The theoretical range for an observable at each point is obtained by varying each input parameter within its respective allowed range and then adding the individual uncertainty in quadrature. If the obtained theoretical range has overlap with the $2\sigma$ range of the experimental data, the point is regarded as allowed. In addition, we consider two different cases: real and complex couplings with respect to $A_u$ and $A_d$. Under the combined constraints from $\Delta M_s$, $\phi^{c\bar{c}s}_s$ and $a_{sl}^s$, the allowed parameter spaces of the type-III and type-C models are shown in Figs.~\ref{fig:realcase} (for the case of real couplings) and \ref{fig:comcase} (for the case of complex couplings), respectively.

\begin{figure}[t]
\begin{center}
  \subfigure[]{\includegraphics[width=7.8cm]{T3Re.pdf}} \qquad
  \subfigure[]{\includegraphics[width=7.8cm]{TCRe.pdf}}
  \caption{\small Allowed spaces of the parameters $A_u$ and $A_d$ in type-III and type-C models under the combined constraints from $\Delta M_s$, $\phi^{c\bar{c}s}_s$ and $a_{sl}^s$, in the case of real couplings. The red, blue and green pointed regions are obtained with $m_{H^{\pm}}=100$, $250$ and $500~{\rm GeV}$, respectively.}
  \label{fig:realcase}
\end{center}
\end{figure}

\begin{figure}[t]
\begin{center}
  \subfigure[]{\includegraphics[width=7.8cm]{T3Com1.pdf}} \qquad
  \subfigure[]{\includegraphics[width=7.8cm]{TCCom1.pdf}}\\
  \subfigure[]{\includegraphics[width=7.8cm]{T3Com2.pdf}} \qquad
  \subfigure[]{\includegraphics[width=7.8cm]{TCCom2.pdf}}
  \caption{\small Allowed spaces of the parameters $|A_u|$, $|A_u^*A_d|$ and $\theta$ in type-III and type-C models under the combined constraints from $\Delta M_s$, $\phi^{c\bar{c}s}_s$ and $a_{sl}^s$, in the case of complex couplings. The other captions are the same as in Fig.~\ref{fig:realcase}.}
  \label{fig:comcase}
\end{center}
\end{figure}

For the case of real couplings, it can be seen from Fig.~\ref{fig:realcase} that:
\begin{enumerate}
\item[ (i) ] In the type-III model, as shown in Fig.~\ref{fig:realcase}(a), the module of Yukawa coupling parameter $A_u$ is severely constrained by the good agreement between the SM prediction and the experimental data for $\Delta M_s$; for instance $|A_u|<1$ is obtained with $m_{H^{\pm}}=500\,{\rm GeV}$. There are, however, almost no constraints on the coupling $A_d$, because the contribution involving it is negligible with respect to the one involving only $A_u$.

\item[ (ii) ] In the type-C model, because the charged-Higgs contribution to ${\cal A}^{VLL}$ is relatively small and large cancelation effects exist between the terms involving $A_d$ and $A_u$, the allowed values of $A_d$ and $A_u$ could be large simultaneously, with either the same or the opposite signs, as shown by the four ``legs" in Fig.~\ref{fig:realcase}(b).

\item[ (iii) ] Besides the ``legs" in Fig.~\ref{fig:realcase}(b), the difference between the two models is also featured by the different shapes of the allowed parameter spaces. The current data on $B_s^0-\bar{B}_s^0$ mixing generally puts a stronger constraint on the type-III model; for instance, with the assumption $|A_d|\sim |A_u|$ and choosing $m_{H^{\pm}}=500\,{\rm GeV}$, the upper bound $|A_u|\sim 1.5$ obtained in type-C model is obviously looser than the one $|A_u|\sim 1$ in type-III model.
\end{enumerate}

For the case of complex couplings, one more model parameter $\theta$ is introduced. From Fig.~\ref{fig:comcase}, it is found that:
\begin{enumerate}
\item[ (i) ] In the type-III model, as shown in Figs.~\ref{fig:comcase}(a) and \ref{fig:comcase}(c), large values of $|A_u|$ and $|A_u^*A_d|$ are still allowed around $\theta\sim \pm90^{\circ}$, which makes it different from the case of real couplings. This is due to the fact that large cancelation effects appear among the charged-Higgs contributions when $\theta\sim \pm90^{\circ}$, which can also be seen from Eqs.~\eqref{eq:AIII2} and \eqref{eq:AIII3}. Moreover, as shown in Fig.~\ref{fig:comcase}(a), an approximately linear relationship is observed between $|A_u^*A_d|$ and $|A_u|$ when $|A_u|\gtrsim0.5$.

\item[ (ii) ] As shown in Figs.~\ref{fig:comcase}(b) and \ref{fig:comcase}(d), similar observations could also be made in the type-C model, except for the fact that the constraints on the model parameters are now much looser. In addition, the cancelation effects among the charged-Higgs contributions occur around $\theta\sim 0^{\circ}$ and $\pm180^{\circ}$, which is different from that observed in the type-III model.

\end{enumerate}

From the above discussions, we conclude that the type-III and type-C models exhibit some significantly different behaviors under the experimental constraints from $B_s^0-\bar{B}_s^0$ mixing. However, due to the large theoretical and experimental uncertainties, the differences in the small $|A_u|$ and $|A_d|$ ranges are hardly to be distinguished from each other. The future refined measurement and precise theoretical evaluation for $B_s^0-\bar{B}_s^0$ mixing might show a much clearer phenomenological picture for the type-III and type-C models.

As a final comment, it should be noted that the same analysis could also be applied to the $B_{d}^0-\bar{B}_d^0$ mixing, which is another important related low-energy process. The charged-Higgs effect on it can be obtained from that on the $B_{s}^0-\bar{B}_s^0$ mixing, with the replacement $s\to d$ throughout the theoretical formulae presented in Sec.~\ref{subsec:wcs}. However, we find that the bounds on the model parameters derive from the $B_{d}^0-\bar{B}_d^0$ mixing are quite similar to the ones from the $B_{s}^0-\bar{B}_s^0$ mixing, and no any further information on the model parameters could be obtained from the former. Therefore, the constraints from $B_{d}^0-\bar{B}_d^0$ mixing will not be shown here.

\section{Conclusion}
\label{sec:conc}

In this paper, we have calculated the one-loop electro-weak corrections to the $B_s^0-\bar{B}_s^0$ mixing within the type-III and type-C 2HDMs with MFV, in which the second scalar doublet is fixed to be color-singlet and color-octet, respectively. It is noted that, in order to get a gauge-independent result, the external momenta of the heavy quarks inside the mesons should be taken into account, and the heavy-quark masses should be kept up to the second order.

Based on the obtained short-distance Wilson coefficients of the four-quark operators appearing in the effective weak Hamiltonian, and combining the up-to-date experimental data on $B_s^0-\bar{B}_s^0$ mixing, we then performed a detailed phenomenological analysis of the charged-Higgs effects on this process. Our main conclusions are summarized as follows:
\begin{enumerate}
\item[ (i) ] While the type-C model gives a nonzero contribution to the Wilson coefficient $C^{TRR}_{C}(\mu_W)$, the type-III contribution to the amplitude ${\cal A}^{TRR}$ is induced only by RG evolution effect.

\item[ (ii) ] In the case of real couplings, the allowed spaces of the Yukawa coupling parameters $A_u$ and $A_d$ in the two models are obviously different, as shown in Fig.~\ref{fig:realcase}.

\item[ (iii) ] In the case of complex couplings, due to the cancelation effects among the charged-Higgs contributions, large values of $|A_u|$ and $|A_d|$ are still allowed around $\theta\sim \pm 90^{\circ}$ in the type-III and around $\theta\sim0^{\circ}\,,\pm180^{\circ}$ in the type-C model, which is shown in Fig.~\ref{fig:comcase}.
\end{enumerate}
The observed differences could be used to distinguish the two models. It should be noted, however, that their differences in the small $|A_d|$ and $|A_u|$ ranges are hardly to be distinguished, due to the large theoretical and experimental uncertainties. More refined theoretical and experimental efforts are therefore needed for a much clearer phenomenological picture.

\section*{Acknowledgments}

We thank Antonio Pich and Martin Jung for useful discussions and cross-checks for the box Feynman diagrams. This work is supported by the National Natural Science Foundation of China (Grant Nos. 11475055, 11435003, 11105043 and 11005032). Q.~Chang is also supported by the Foundation for the Author of National Excellent Doctoral Dissertation of P.~R.~China~(Grant No.~201317) and the Program for Science and Technology Innovation Talents in Universities of Henan Province~(No.~14HASTIT036). X.~Q.~Li is also supported by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, by the Open Project Program of SKLTP, ITP, CAS, China~(No.~Y4KF081CJ1), and by the self-determined research funds of CCNU from the colleges' basic research and operation of MOE~(CCNU15A02037).

\begin{appendix}

\section{Relevant coefficients for $B_s^0-\bar{B}_s^0$ mixing}
\label{appendix:1}

Here we present the explicit expressions for $f_i(x_t,x_h)$ appearing in Eqs.~\eqref{eq:Ci_III} and \eqref{eq:Ci_C}:
\begin{align} \label{fi_def}
f_1(x_t,x_h) & = \frac{1}{2}\left[\frac{x_t^2(x_t-4)}{(x_h-x_t)(1-x_t)} +\frac{x_t^2(3x_t^2-x_h(4-2x_t+x_t^2))\ln{x_t}}{(x_h-x_t)^2(1-x_t)^2}
-\frac{(x_h-4)x_hx_t^2\ln{x_h}}{(1-x_h)(x_h-x_t)^2}\right]\,,\\[0.2cm]
f_2(x_t,x_h) & = \frac{1}{2}\left[\frac{x_t^2(x_h+x_t)}{2(x_h-x_t)^2} +\frac{x_hx_t^3\ln{x_t}}{(x_h-x_t)^3}-\frac{x_hx_t^3\ln{x_h}}{(x_h-x_t)^3}\right]\,,\\[0.2cm]
f_3(x_t,x_h) & = \frac{1}{3} \bigg[-\frac{x_t^2(x_t^2+x_h^4)(-11+7x_t-2x_t^2)+x_hx_t^3(7+53x_t-55x_t^2+19x_t^3)}{3(1-x_h)^2(x_h-x_t)^3(1-x_t)^3}\nonumber\\[0.1cm]
&\hspace{1.1cm}-\frac{x_h^2x_t^2(-2-55x_t+15x_t^2+17x_t^3-11x_t^4)+x_h^3x_t^2(19+17x_t-19x_t^2+7x_t^3)}{3(1-x_h)^2(x_h-x_t)^3(1-x_t)^3}\nonumber\\[0.1cm]
&\hspace{1.1cm}+\frac{2x_t^2(x_h^3-3x_h^2x_t+3x_hx_t^2-3x_t^4+3x_t^5-x_t^6)\ln{x_t}}{(x_h-x_t)^4(1-x_t)^4}\nonumber\\[0.1cm]
&\hspace{1.1cm}-\frac{2x_hx_t^2(x_h^2+(-3+x_h)x_hx_t+(3+(-3+x_h)x_h)x_t^2)\ln{x_h}}{(1-x_h)^3(x_h-x_t)^4}\bigg]\,,\\[0.2cm]
f_4(x_t,x_h) & = -\frac{x_t^2((x_h^2+x_t)(-3+x_t)+x_h(1+6x_t-3x_t^2))}{2(1-x_h)(x_h-x_t)^2(1-x_t)^2}\nonumber\\[0.1cm]
&\hspace{0.46cm}-\frac{x_t^2(x_h^2-2x_hx_t-(-2+x_t)x_t^3)\ln{x_t}}{(x_h-x_t)^3(1-x_t)^3}+\frac{x_hx_t^2(x_h-2x_t+x_hx_t)\ln{x_h}}{(1-x_h)^2(x_h-x_t)^3}\,,\\[0.2cm]
f_5(x_t,x_h) & = -\frac{2x_t^2}{(x_h-x_t)^2}-\frac{x_t^2(x_h+x_t)\ln{x_t}} {(x_h-x_t)^3}+\frac{x_t^2(x_h+x_t)\ln{x_h}}{(x_h-x_t)^3}\,,\\[0.2cm]
f_6(x_t,x_h) & = \frac{1}{6} \bigg[-\frac{x_t^2(5x_h^2-22x_hx_t+5x_t^2)}{3(x_h-x_t)^4} -\frac{x_t^2(x_h^3-3x_h^2x_t-3x_hx_t^2+x_t^3)\ln{x_t}}{(x_h-x_t)^5}\nonumber\\[0.1cm]
&\hspace{1.1cm}+\frac{x_t^2(x_h^3-3x_h^2x_t-3x_hx_t^2+x_t^3)\ln{x_h}}{(x_h-x_t)^5}\bigg]\,,\\[0.2cm]
f_7(x_t,x_h) & = \frac{2x_t^2}{(x_h-x_t)^2}+\frac{x_t^2(x_h+x_t)\ln{x_t}}{(x_h-x_t)^3} -\frac{x_t^2(x_h+x_t)\ln{x_h}}{(x_h-x_t)^3}\,.
\end{align}

\end{appendix}


\section*{Introduction}

For many complex diseases such as cancer, diabetes, and cardiovascular diseases, prognosis is of essential interest. In the omics era, profiling studies have been extensively conducted, searching for genetic markers associated with prognosis. It has been suggested that, beyond the main effects of genetic (G) and environmental (E) risk factors, the gene-environment (G$\times$E) interactions also have important implications. Multiple statistical methods have been developed for G$\times$E interaction analysis. For reviews, refer to~\cite{Caspi2006gene},~\cite{Cordell2009detecting},~\cite{Thomas2010methods}, and others.

Denote $T$ as the prognosis time of interest, $X = (X_1, \ldots, X_q)$ as the $q$ environmental/clinical variables, and $Z=(Z_1, \ldots, Z_p)$ as the $p$ genetic variables. Assume $n$ independent subjects. Regression-based interaction analysis, with its broad applicability, has been extensively adopted and proceeds as follows. (a) For gene $j (=1, \ldots, p)$, consider the model $T \sim \phi( \alpha_{j,0}+\sum_{k=1}^qX_k\alpha_{j,k} + Z_j\beta_j + Z_j\sum_{k = 1}^qX_k\gamma_{j,k})$, where $\phi(\cdot)$ is the known link function (for example, the Cox or exponential model), and $\alpha_{j,0},\alpha_{j,k}, \beta_{j}, \gamma_{j,k}$ are the unknown regression coefficients. As usually $q<<n$, this is a low-dimensional model and can be fitted using standard, usually likelihood-based, techniques and software. Denote $p_{j,k}$ as the p-value for $\gamma_{j,k}$. (b) With $\{p_{j,k}: j=1, \ldots, p; k=1, \ldots, q \}$, conduct multiple comparison adjustment using the Bonferroini or FDR (false discovery rate) approach, and identify important interactions.

A limitation of the above approach is a lack of robustness properties. Usually it is assumed that all subjects satisfy the same prognosis models. In practice, most genetic studies cannot afford to conduct strict subject selection. Seemingly homogeneous subjects can have different disease subtypes~\citep{Haibe-Kains2012three} and different survival patterns~\citep{burgess2011cancer}. Cause of death can be misclassified, leading to contamination in disease-specific survival~\citep{fall2008reliability}. The survival times extracted from medical records are not always reliable (\cite{bowman2011doctors},~\cite{rampatige2013assessing}). With unrobust for example likelihood-based estimation, even a single contaminated observation can result in severely biased estimates~\citep{huber2009robust}, which can lead to false marker identification. Another limitation is that, significance level-based identification, although {\it asymptotically} valid, may generate unreliable results when the sample sizes are small to moderate, as in typical profiling studies. A recent study suggests that regularized estimation can lead to more reliable estimation and hence more accurate marker identification~\citep{shi2014gene}.

With low-dimensional biomedical data, robust statistical methods have demonstrated great power. As suggested in a recent review by~\cite{wu2014selective}, with genetic data, the development is limited and unsystematic and has been mostly on the analysis of main effects but not interactions. In the literature, relevant studies include~\cite{gui2011novel}, which identifies important interactions using the multifactor dimensionality reduction (MDR) technique. However, this method is limited to categorial data (such as SNPs) and not broadly applicable.~\cite{shi2014gene} developed a rank-based method, which is robust to model mis-specification but not data contamination. The most relevant study is perhaps~\cite{wang2015identifying}, which developed a quantile-regression based method. With that method, the quantile needs to be specified, which is not a trivial task in practical data analysis. The objective function is not differentiable, causing difficulty in estimation. In addition, as to be shown in this study, its numerical performance can be less satisfactory under many data settings. With low-dimensional data, it has been shown that no robust method dominates the others. It is thus prudent to develop alternative robust methods.

Consider prognosis data with both G and E measurements. The goal is to develop a new method for identifying important G$\times$E interactions. Significantly advancing from the existing studies, the proposed method is robust to contamination in the prognosis data. In addition, our simulation suggests that, under certain scenarios, its numerical performance is better than the quantile regression analysis, which is perhaps the most popular robust method for genetic data~\citep{wu2014selective}. A penalization approach is adopted for marker identification, which differs from the significance level-based approach and can have better numerical performance when the sample size is small to moderate.

\section*{Robust identification of interactions}

\subsection*{Data and model settings}

In the literature, there are multiple prognosis models. In this study, we consider the accelerated failure time (AFT) model, which has been adopted in~\cite{shi2014gene}, \cite{liu2013identification}, and many others. Compared to alternatives such as the Cox model, the AFT model can be preferred because of its more intuitive interpretations and lower computational cost, which are especially desirable with high-dimensional genetic data. With a slight abuse of notation, still $T$ to denote the logarithm of prognosis time. For gene $j$, the AFT model assumes that
$$T = \alpha_{j,0}+\sum_{k=1}^qX_k\alpha_{j,k} + Z_j\beta_j + Z_j\sum_{k = 1}^qX_k\gamma_{j,k}+\epsilon,$$
where $\epsilon$ is the random error.

Assume $n$ independent observations. Consider the scenario where a small subset of the random errors are contaminated, leading to contamination in the prognosis times. In a typical profiling study, $q<<n$, while $p$ can be comparable to or much larger than $n$. For subject $i$, denote $C_i$ as the logarithm of censoring time and $\mathbf{x}_i$ and  $\mathbf{z}_i$ as the observed $X$ and $Z$ values, respectively. Under right censoring, we observe $(y_i=min(T_i, C_i), \delta_i=I(T_i \leq C_i), \mathbf{x}_i, \mathbf{z}_i)$. Further denote $\mathbf{u}_{i,j}=(1,\mathbf{x}^{\prime}_i, z_{i,j},z_{i,j}x_{i,1},\ldots,z_{i,j}x_{i,q})'$,  $\bm\zeta_j=(\alpha_{j,0},\alpha_{j,1},\ldots,\alpha_{j,q},\beta_j,\gamma_{j,1},\ldots, \gamma_{j,q})'$, and $\mathbf{U}_j=(\mathbf{u}_{1,j}, \ldots, \mathbf{u}_{n,j})'$ which is a $n\times (2q+2)$ matrix. For gene $j$ and subject $i$, the AFT model can now be written as
$T_i = \mathbf{u}_{i,j}^\top\bm{\zeta}_j + \epsilon_{i,j}.$ Without loss of generality, assume that the data $\{(y_i,\delta_i,\mathbf{x}_i,\mathbf{z}_i), i=1,\ldots, n\}$ have been sorted according to $y_i$'s from the smallest to the largest.

\subsection*{Penalized robust identification}

A penalized marker identification method is defined by its loss function and penalty, which are defined separately as follows.

\noindent{\bf Loss function} Before defining the robust loss function, we take one step back and consider the scenario with no data contamination. When the distribution of $\epsilon$ is not specified, likelihood-based estimation cannot be adopted. A popular approach is the weighted least squared estimation proposed in~\cite{Stute1993consistent} and proceeds as follows. First compute the Kaplan-Meier weights as $$\omega_{1}=\frac{\delta_{1}}{n},~\omega_{i}=
\frac{\delta_{i}}{n-i+1}\prod_{j=1}^{i-1}
(\frac{n-j}{n-j+1})^{\delta_{j}},~i=2,\cdots,n.$$
For gene $j(=1,\ldots, p)$, the weighted least squared objective function is defined as
\begin{eqnarray}
\sum_{i=1}^n \omega_{i} (y_i - \mathbf{u}_{i,j}^\top\bm\zeta_j)^2,
\end{eqnarray}
This loss function has the following properties. With a simple linear regression model, the most commonly adopted loss function is the least squared one. Here to accommodate censoring, a weight function is imposed to re-weigh different observations according to their observed times and event status. When there is no censoring, $\omega_i=1/n$. With the quadratic form, the loss function is not robust to data contamination. If subject $i$ is not censored, then $\omega_i\neq 0$, and an arbitrarily large $y_i$ results in arbitrarily large estimates. Biased estimation can happen with data contamination, which may lead to false marker identification.

To accommodate data contamination, for gene $j(=1,\ldots, p)$, we propose the exponential squared loss function
\begin{eqnarray}
Q_{\theta}(\bm\zeta_j|\mathbf{y}, \mathbf{U}_j,\bm\omega)=
\sum_{i=1}^n\omega_i\exp(-(y_i - \mathbf{u}_{i,j}^\top\bm\zeta_j)^2 / \theta).
\end{eqnarray}
$\mathbf{y}$ and $\bm\omega$ denote the vectors composed of $y_i$'s and $\omega_i$'s, respectively. $\theta>0$ is a tuning parameter. The rationale of this approach is as follows. For contaminated subjects with $y_i$s deviate from $\mathbf{u}_{i,j}^\top\bm\zeta_j$s (predicted values from the model), $(y_i - \mathbf{u}_{i,j}^\top\bm\zeta_j)^2$s have large values.
The exponential function down-weighs such contaminated observations. The degree of down-weighing is adjusted by $\theta$: with $\theta$ getting smaller, contaminated observations have smaller influence. To accommodate censoring, $\omega_i$'s are imposed in a similar manner as the original Stute's approach. For low-dimensional linear regression model without censoring, the exponential squared loss has been examined in~\cite{wang2013robust}. Different from the existing studies, here we consider the more challenging high-dimensional genetic data especially interactions. In addition, the Kaplan-Meier weights are introduced to accommodate censoring. As to be shown in Appendix, such differences lead to significant differences in the statistical development.

\noindent{\bf Penalized estimation}
For gene $j(=1,\ldots, p)$, we propose the penalized objective function
\begin{eqnarray}
L_{\lambda,\theta}(\bm\zeta_j|\mathbf{y}, \mathbf{U}_j,\bm\omega) = Q_{ \theta}(\bm\zeta_j|\mathbf{y}, \mathbf{U}_j,\bm\omega)-\lambda||\bm\zeta_j||_1.
\label{eqn2}
\end{eqnarray}
$\lambda>0$ is the data-dependent tuning parameter, and $||\cdot||_1$ is the $\ell_1$ norm. Denote $\tilde{\bm\zeta_j}$ as the maximizer of $L_{\lambda,\theta}(\bm\zeta_j|\mathbf{y}, \mathbf{U}_j,\bm\omega)$. Interactions (and main effects) corresponding to the nonzero components of $\tilde{\bm\zeta_j}$ are identified as important.

The strategy of using penalization for identifying important interactions has been adopted in~\cite{shi2014gene} and other studies. Here a Lasso penalty is imposed, which is the most commonly adopted penalty. $\lambda$ controls the degree of sparsity and number of identified interactions. We impose the same $\lambda$ to all genes to ensure comparability.

Multiple penalties can take place of Lasso. In some recent studies such as~\cite{bien2013lasso} and~\cite{liu2013identification}, penalties have been developed to respect the ``main effects, interactions" hierarchy, which reinforces that the main effects corresponding to the identified interactions must be identified. In a large number of studies (\cite{zimmermann2011interaction},~\cite{Caspi2006gene}), it has been observed that genes can have important G$\times$E interactions but no main effects. In addition, if the hierarchy has to be reinforced, we can identify the important interactions first, and then add back the corresponding main effects. Computationally, Lasso is much simpler than the existing alternatives. Our limited experience suggests that with the complex robust loss function, complex penalties have a considerable probability of running into convergence problems. With the above considerations, the Lasso penalty is adopted here.

\noindent{\bf Consistency properties} A significant advantage of the proposed method is that it has the much deserved consistency properties under the ultrahigh dimensional setting. For many of the existing interaction analysis methods, there is a lack of theoretical development. The rigorously established consistency properties provide a solid statistical ground for the proposed method. Details are provided in Appendix.

\subsection*{Computation}

We propose a coordinate-wise updating procedure to compute the solution to~\eqref{eqn2}. For low-dimensional data under simpler settings, an iterative approach is suggested in~\cite{wang2013robust} to select the robust tuning parameter $\theta$. However, under the present high-dimensional settings and with the coordinate-wise updating procedure, such an approach is computationally infeasible. Alternatively, for each $(\lambda,\theta)$ pair, we compute the solution to each marginal model. This way, we can obtain a solution surface over the two-dimensional tuning parameter grid. Then a prediction error based method such as cross-validation can be used to select the robust and penalization tuning parameters. Computation is conducted for each gene separately and can be realized in a highly parallel manner to reduce computer time.

Consider gene $j(=1,\ldots, p)$. Let $r_i(\bm\zeta_j)=y_i-\mathbf{u}_{i,j}^\top\bm\zeta_j$. The first and second order derivatives of $Q(\bm\zeta_j)$ are
\begin{equation}
\begin{cases}
\dot{Q}_k(\bm\zeta_j)=\frac{\partial Q(\bm\zeta_j)}{\partial \zeta_{j,k}}=2\sum_{i=1}^n\omega_iu_{(i,j)_k}r_i(\bm\zeta_j)
\exp(-r_i^2(\bm\zeta_j)/\theta) /\theta, \\
\ddot{{Q}}_{kl}(\bm\zeta_j)=\frac{\partial^2 Q(\bm\zeta_j)}{\partial \zeta_{j,k}\partial \zeta_{j,l}}=2\sum_{i=1}^n\omega_iu_{(i,j)_k}u_{(i,j)_l}\exp(-r_i^2(\bm\zeta_j)/\theta) (2r_i^2(\bm\zeta_j)/\theta-1)/\theta.
\end{cases}
\label{eqn4}
\end{equation}
For a $\bm{\zeta}_j^m$ that is in a small neighborhood of $\bm{\zeta}_j$,  $Q(\bm\zeta_j)$ in~\eqref{eqn2} can be locally approximated by
\begin{equation}
Q(\bm\zeta_j) \approx Q(\bm\zeta^m_j) + \dot{Q}(\bm\zeta^m_j)^\top(\bm\zeta_j-\bm\zeta_j^m)+
\frac{1}{2}(\bm\zeta_j-\bm\zeta^m_j)^\top\ddot{Q}
(\bm\zeta^m_j)(\bm\zeta_j-\bm\zeta^m_j).
\label{eqn5}
\end{equation}
Replacing $Q(\bm\zeta_j)$ in~\eqref{eqn2} with this approximation and taking the first order derivative of $L(\bm\zeta_j)$ with respect to the $k$th element $\zeta_{j,k}$ give us
\begin{equation}
\zeta_{j,k}^{m+1}=\zeta_{j,k}^m -\ddot{Q}^{-1}_{kk}(\bm{\zeta}^m_j)\dot{Q}_k(\bm{\zeta}^m_j)+ \ddot{Q}_{kk}^{-1}(\bm{\zeta}^m_j)\lambda.
\end{equation}
Note that when $2r_i^2(\bm{\zeta}_j)/\theta>1$, $\ddot{Q}_{kk}(\bm{\zeta}_j)\geq0$. Then \eqref{eqn2} can be maximized at infinity, and the algorithm may fail to converge. To avoid this problem, notice that $$\ddot{Q}_{kk}(\bm{\zeta}_j)\geq -2\sum_{i=1}^n\omega_iu_{(i,j)_k}^2\exp(-r_i^2(\bm\zeta_j)/\theta)/\theta$$ for $\omega_i\geq 0$. The right hand side is always non-positive. We re-define $$\ddot{Q}_{kk}(\bm{\zeta}_j)\equiv -2\sum_{i=1}^n\omega_iu_{(i,j)_k}^2\exp(-r_i^2(\bm\zeta_j)/\theta)/\theta,$$ and use the minorization-maximization (MM) algorithm to ensure convergence. The algorithm that combines the coordinate descent and MM algorithm is summarized in Algorithm~\ref{alg1}.
\begin{algorithm}
\caption{The coordinate descent MM algorithm for marginal model $j$}
\label{alg1}
\SetKwInOut{Input}{input}\SetKwInOut{Output}{output}
\Input{$\mathbf{y}, \bm{\omega}, \mathbf{U}_j$, and tuning parameters $\lambda$ and $\theta$.}
\Output{the maximizer $\tilde{\bm{\zeta}}_j$ defined in~\eqref{eqn2}.}
\BlankLine
\textbf{initialization}: let $m = 0, \bm{\zeta}^0_j=\mathbf{0}$. Normalize $\mathbf{y}$ and $\mathbf{U}_j$ such that $\sum_{i=1}^n\omega_iy_i=0$, $\sum_{i=1}^n\omega_i u_{(i,j)_k}=0$ and $\sum_{i=1}^n\omega_iu_{(i,j)_k}^2=n$ for $k=1,\ldots,2q+2$.\;
\Repeat{$||\bm{\zeta}^{m_0}_j - \bm{\zeta}^{m}_j||_2\leq$ a predefined thresholding value}{
$m_0 = m$\;
\For{$k = 1, 2, \ldots, 2q+2$}{
\Repeat{$|{\zeta}^{m}_{j,k} - {\zeta}^{m-1}_{j,k}|\leq$ a predefined thresholding value. In our numerical study, we set the thresholding =$10^{-3}$}{
$\zeta_{j,k}^{m+1} = \zeta_{j,k}^m - \ddot{Q}^{-1}_{kk}(\bm{\zeta}^m_j)\dot{Q}_k(\bm{\zeta}^m_j)+ \ddot{Q}_{kk}^{-1}(\bm{\zeta}^m_j)\lambda$\;
$\zeta_{j,l}^{m+1} = \zeta_{j,l}^m$ for $l = 1, 2, \ldots, 2q+2$, $l\neq k$\;
$m = m + 1$\;
}
}
}
\Return{$\bm{\zeta}^m_j$} as the maximizer.
\end{algorithm}

As mentioned above, we search over a two-dimensional grid of $(\lambda,\theta)$ for the optimal. The range of $\lambda$ is determined as follows. First, its upper bound $\lambda_{\max}$ is selected such that $\tilde{\bm{\zeta}}_j=\mathbf{0}$ for $j=1,\ldots,p$. With the unrobust weighted least squared loss, $\lambda_{\max} = \max_{j=1}^p\{||\mathbf{U}_j^\top W \mathbf{y}||_{\infty}\}$ where $W$ is the diagonal matrix composed of $\omega_i$s. With the proposed robust method, the derivatives in~\eqref{eqn4} can be viewed as a weighted sum of $u_{(i,j)_k}r_i(\bm{\zeta}_j)$s. Because the weight for each subject changes with $\bm{\zeta}_j$, the previously defined $\lambda_{\max}$ may not guarantee that $\tilde{\bm{\zeta}}_j=\mathbf{0}$ for $j$. After some trials, we find that $\lambda_{\max} = 20\max_{j=1}^p\{||\mathbf{U}_j^\top W \mathbf{y}||_{\infty}\}$ is in general a safe upper bound for $\lambda$. The lower bound is chosen as $\lambda_{\min}=\lambda_{\max} / 1000.$ The range of $\theta$ depends on the degree of contamination in data, which is not known in practice. In our numerical studies, we find that the estimator is not very sensitive to the value of $\theta$. For cautions, we choose a relatively wide range for $\theta$. Specifically, after centralization, we choose $\theta \in (\min_{i=1}^n{y_i^2}/100, \max_{i=1}^ny_i^2 \times 100)$.

With the quadratic approximation and minorization, optimization in~\eqref{eqn2} is non-convex. Therefore, there is a risk of converging to local maximums. Thus when trying to solve the estimate at a new $(\lambda, \theta)$ point, instead of using $\bm{\zeta}^0_j=\mathbf{0}$ as the starting value, we use the estimate at a neighboring $(\lambda, \theta)$ point as the warm starting value. This modification speeds up the computation and also improves the convergence property.

\section*{Simulation}

In simulation, we set $n= 300$, $q=3$, and $p = 500, 1000$. There are a total of 18 nonzero effects: 3 main effects of E factors, 5 main effects of G factors, and 10 interactions. The positions of nonzero main G effects and interactions are uniformly placed. The nonzero regression coefficients are randomly generated from $uniform(0.5, 1.5)$. The E and G factors are generated from multivariate normal distributions with marginal means zero, marginal variances one, and the following variance matrix structures: Independent, AR(0.2), AR(0.8), Band(0.3), and Band(0.6). Under the AR($\rho$) correlation structure, for the $i$th and $j$th factors, $corr=\rho^{|i-j|}$. Under the Band($\rho$) correlation structure, for the $i$th and $j$th factors, $corr=\rho I(|i-j|\leq 2)$, where $I(\cdot)$ is the indicator function. Under each correlation scenario, consider seven different distributions for the random error $\epsilon$: N(0,1), 0.95N(0,1)+0.05Cauchy, 0.85N(0,1)+0.15Cauchy, 0.7N(0,1)+0.3Cauchy, 0.95N(0,1)+0.05t(3), 0.85N(0,1)+0.15t(3), and 0.7N(0,1)+0.3t(3). That is, we consider the scenarios with no contamination and three different levels of contamination. Two contamination distributions are considered with different thickness of tails. The log event times are generated from the AFT models. The censoring times are generated independently from exponential distributions. The parameters are adjusted so that the censoring rates are about 25\%.

Beyond the proposed method (referred to as ``Robust"), we also consider the following three alternatives: (a) the ``Unrobust" method, which adopts the weighted least squared loss function and applies the Lasso penalization for selecting important effects; (b) the ``Stute" method, which adopts the weighted least squared loss function, does not apply any penalization, and uses the significance level (p-value) as the criterion for quantifying the importance of effects; and (c) the ``Quantile" method, which adopts the quantile regression-based robust loss function and applies the Lasso penalization for selecting important effects. The Unrobust and Stute methods, as the proposed, are built on the weighted least squared estimation. Comparing with them can establish the advantage of robust loss and penalization-based identification, respectively. The Quantile method is also robust and has been recently developed. Comparing with it can establish the advantage of the weighted exponential squared loss. We acknowledge that multiple methods are potentially applicable to the simulated data. The three alternatives have frameworks closest to that of the proposed.

With the Robust, Unrobust, and Quantile methods, the numbers of selected interactions depend on the tuning parameter values. With the Stute method, the number depends on the p-value cutoff. To eliminate the effect of tuning parameter selection on identification, we examine a sequence of tuning parameter values, evaluate identification performance at each value, and use the ROC (receiver operating characteristics) curve to evaluate the interaction identification accuracy of different methods. A representative ROC plot is shown in Figure 1. In this plot, the proposed method has the dominatingly better accuracy.

Summary AUCs based on 100 replicates are shown in Tables 1 ($p=500$) and 2 ($p=1000$). When there is no contamination, performance of the proposed method is comparable to or slightly worse than that of the unrobust method. For example with $p=500$ and the independence correlation structure, the robust and unrobust methods have mean AUCs 0.861 and 0.889, respectively. And with the AR(0.2) correlation, they have mean AUCs 0.901 and 0.881, respectively. With contamination, the robust method outperforms the unrobust method. For example with $p=500$, the AR(0.2) correlation structure, and 0.7N(0,1)+0.3Cauchy error, the robust and unrobust methods have mean AUCs 0.886 and 0.751, respectively. Under all simulation scenarios, the Stute method, which adopts the robust loss function but significance level-based identification, has inferior identification accuracy. As has been suggested in the literature, with a moderate sample size, the unregularized estimates can be less reliable, leading to inaccurate identification. When comparing the proposed method with the quantile regression-based, we see that under the majority of the settings, the proposed method has superior performance. For example with $p=500$, AR(0.2) correlation, and 0.85N(0,1)+0.15Cauchy error, the proposed method and Quantile have mean AUCs 0.892 and 0.842, respectively. However, under a small number of settings, the Quantile method excels. For example with $p=500$, AR(0.8) correlation, and 0.7N(0,1)+0.3t(3) error, the two methods have mean AUCs 0.863 and 0.933, respectively. Such observations are also reasonable. No robust approach is expected to be able to dominate all others. The proposed method outperforms the quantile regression-based under most scenarios and provides a useful alternative.

\section*{Analysis of the lung squamous cell carcinoma data}

Lung squamous cell carcinoma is the second most common lung cancer and causes around 400,000 deaths each year worldwide~\citep{comprehensive2012TCGA}. Profiling studies have been extensively conducted searching for its prognostic factors. In this section, we analyze the TCGA (The Cancer Genome Atlas) data on the prognosis of lung squamous cell carcinoma. The TCGA data were recently collected and published by NCI and have a high quality. The dataset we analyze was downloaded in April of 2015.

The prognosis outcome of interest is overall survival. Multiple environmental and clinical variables are available for analysis. We first remove variables with a low quality of measurements (for example with a high missing rate). We then select the following variables for analysis: age, gender (female is coded as baseline), smoking pack years, and smoking status (non-smoker, reformed smoker for more than 15 years, reformed smoker for less than or equal to 15 years, current smoker; coded as 0, 1, 2 and 3). These variables have been suggested in the literature (\cite{miller2004Bronchioloalveolar} and~\cite{nakachi1991genetic}). A total of 422 samples have clinical and environmental measurements available.

For the genetic part, we analyze gene expression data. A total of 18,969 measurements are available for 502 samples. When matching the clinical/environmental data with genetic data, complete data are available for 404 samples. Among them, 129 died during following. The median followup was 30 months. 275 were censored, and the median followup was 18 months.

The tuning parameter can be chosen using data-driven methods for example cross validation. However, we note that the commonly used methods are mostly based on the notion of prediction. With a large number of marginal models, the goal is to identify markers top-ranked in a marginal sense,
not prediction. Thus, following published studies, we vary the tuning so that a prefixed number of interactions are selected. In Table~\ref{tab3}, we provide results on the 33 top-ranked interactions. Longer or shorter lists of identified interactions are available from the authors.
After the interactions are identified, we refit the marginal models without penalization on the main effects to satisfy the ``main effects, interactions" hierarchy. Beyond the proposed method, we also apply the three alternatives considered in simulation. Table~\ref{tab7} in Appendix suggests that different methods identify significantly different genes and interactions. Results under the proposed method are provided in Table~\ref{tab3}. Those under the alternative methods are provided in Appendix.

To assess the stability of our findings, we apply a cross validation-based approach. Specifically, one sample is removed from data, and the proposed method is applied. This step is repeated over all samples. For each identified interaction in Table 3, we compute its probability of being identified in the reduced datasets. It can be seen that three interactions have very small stability measures. All other interactions have stability measures close to 1. We have also examined those interactions not identified and found that their stability measures are all close to 0. This analysis suggests a certain degree of stability of the proposed method.

Literature search suggests that the identified interactions and corresponding genes may have important implications. Specifically, previous studies \citep{shi2014gene} have suggested that the major G$\times$E interactions occur between genes and smoking status. Among the identified gene-smoking interactions in the current study, the CEBPB (CCAAT/enhancer-binding proteins beta) protein is a transcription factor that works with other CCAAT/enhancer-binding protein family members in the regulation of cell cycle progress, differentiation and pro-inflammatory gene expression. CEBPB has been found to be upregulated by tobacco smoke in both human lung fibroblast~\citep{Miglino2012Cigarette} and mice emphysema~\citep{hirama2007increased}. Another gene that interacts with smoking is
EFNA1 (a.k.a. ephrin-A1). A recent {\em in vitro} study found that ephrin-A1 is overexpressed in tobacco smoke treated bronchial airway epithelial cells compared to control cells~\citep{nasreen2014tobacco}. In addition, the elevated expressions of ephrin-A1 are positively associated with tumor proliferative capacity in non-small cell lung carcinoma patients~\citep{giaginis2014ephrin}. The gene that shows both a strong main effect and interaction with duration of smoking is TERF1 (telomeric repeat binding factor 1). The TERF1 expression levels have been found to be decreased in lung cancer (\cite{lin2006expression},~\cite{hu2009analysis}), as well as other types of cancer in several studies (\cite{yamada2002down},~\cite{miyachi2002correlation}). It functions as an inhibitor of telomerase and is identified as a prognostic marker for overall survival in non-small cell lung cancer~\citep{hu2009analysis}. The molecular mechanism of why and how TERF1 decreases in the process of cancer is not clear. Our results suggest that smoking can be one of the factors.
Other than smoking duration, we also find that several genes interact with smoking intensity as well. The gene that interacts with both gender and smoking intensity is KATNB1. KATNB1 encodes protein katanin p80 subunit B1, which has been found to participate in cytokinesis by interacting with tumor suppressor gene LAPSER1. The disruption of cytokinesis process may potentially cause genetic instability and cancer~\citep{Sudo2007LAPSER1}.
In addition to smoking, we find eight genes interact with gender in Table~\ref{tab3}. Among these, STRADB and CA5BP1 draw our attention. STRADB is an important gene in lung cancer progression and metastasis through the activation of LKB1. LKB1 is essential for G1 cell cycle arrest, cell polarity and stress, cell detachment and adhesion~\citep{marcus2010LKB1}. The STRADB encoded protein also interacts with the X chromosome-linked inhibitor of apoptosis protein by enhancing its anti-apoptotic activity. In addition, gene CA5BP1 is located on X chromosome and found to be gender-associated.

\section*{Discussion}
G$\times$E interactions have important implications for the prognosis of a large number of complex diseases. In this study, we have developed a novel interaction analysis method. It is robust to contamination in the prognosis outcome, which is not rare in practice but cannot be accommodated using most of the existing methods. In addition, we adopt penalization for identifying important interactions. This strategy differs from the commonly adopted significance level-based identification and can have better marker identification accuracy. Significantly advancing from most of the existing interaction studies, the consistency property has been rigorously established under the ultrahigh dimensional setting, making the proposed method one of the very few with a strong theoretical basis. An effective computational algorithm has been developed. In simulation study, the proposed method is observed to have satisfactory performance. With contamination, robust methods are needed. Under the majority of the cases, the proposed method outperforms the quantile regression based method. Thus, it provides a useful alternative to the quantile regression-based and other existing methods. In the analysis of lung cancer data, the proposed method identifies meaningful genes and interactions, which differ from those using the alternatives, with satisfactory stability.

We choose the popular Lasso penalty because of its low computational cost and satisfactory performance. A limitation of Lasso is that its results may not respect the ``main effects, interactions" hierarchy. However, as shown in Appendix, the proposed method can consistently identify the important interactions. As demonstrated in data analysis, once the important interactions are identified, if the hierarchy is desired, the main effects can be added back. With a robust loss function, the proposed method is computationally more expensive. However, as the computing can be executed in a highly parallel manner, and with the effective coordinate descent algorithm, the computer time is very much affordable. A total of 70 scenarios are considered in simulation and show the satisfactory performance of proposed method. More comprehensive comparisons of different methods, especially robust methods, have not been conducted in the literature and are of interest to future studies. Interesting findings are made in the lung cancer prognosis data analysis. Unfortunately there is still a lack of objective way of determining which set of results is more sensible. Validation of the findings will be pursued in future studies.