\section{Introduction}
Climate predictions generated by general circulation models should be associated with an estimate for the uncertainties in order to provide meaningful and practical information.
A common method for assessing the uncertainties in climate predictions is by means of ensembles. The ensemble may consist of different values of the parameters representing their uncertainties, different climate models, different initial conditions, and different initialization methods, as well as any other group of predictions representing the uncertainties in our knowledge of the climate system or the characteristics of the stochastic processes involved.
Often, the ensemble average is used to generate the predictions, and the ensemble spread (the spread is usually quantified either by the standard deviation (STD) for large ensembles or by the entire range spanned by the ensemble predictions for small ensembles) is used to estimate the uncertainties \citep{Deque2007,Smith2013,Strobach_2015b,Woldemeskel2016,Hawkins2016}. In order to consider the ensemble spread as a meaningful estimate for the uncertainties, one must assume a certain probability density function (PDF) of the predicted variable and must consider the predictions of the ensemble members as independent and identically distributed (i.i.d.) variables. These assumptions are not expected to be valid because the different predictions are not truly independent and the real PDF of the variable often deviates from the assumed Gaussian distribution. Moreover, there is no guarantee that the ensemble predictions represent a sample drawn from the actual PDF of the predicted variable. In the context of predictions, the uncertainties are often estimated by the error of the model in predicting past conditions \citep{Deque2007}. The past error is then used to estimate the uncertainties of future predictions. This estimation removes the need for assumptions regarding the distribution of the ensemble predictions, but it still assumes that the errors in the past can serve as a good estimator for future predictions.

Several studies in the past suggested various methods to empirically improve the estimation of the uncertainties.
These methods can be divided into two main groups--regression methods and dressing kernel methods.
The regression methods ``weight'' the ensemble members in order to minimize a predefined loss function that accounts for the ensemble error and spread. These methods rely on the assumption that the ensemble predictions are i.i.d.
Common regression methods are the \textit{logistic regression} \citep{Wilks_2006,Wilks_2007} and the \textit{non-homogeneous Gaussian regression} \citep{Gneiting_2005}.
The dressing kernel methods assume a specific parametric PDF for each of the ensemble members (the characteristics of the PDF are based on the ensemble errors in the past), and the uncertainty is aggregated from the PDFs of all the ensemble members. The dressing kernel methods differ in their assumptions regarding the shape and the width of the ensemble member PDFs. For example, Roulston and Smith \citep{Roulston_2003} assumed Gaussian PDFs and determined their variance from the error spread of the best ensemble member. Wang and Bishop \citep{Wang_2005} determined the variance from the difference between the variance of the ensemble mean errors and the ensemble variance, and Raftery et al. \citep{Raftery_2005} used Bayesian inference to optimize the variance of the dressing kernel. A more thorough  discussion on the various methods can be found in \citep{Wilks_2011}.

Here, we suggest a method for improving the uncertainty estimation by bridging the spread- and error-based estimations. Our method does not rely on any assumption regarding the distribution of the ensemble member predictions or their i.i.d. characteristics. We use the past observed relation between the ensemble mean error and the ensemble spread in order to estimate future prediction uncertainties. Our method relies on the assumption that the relation between the spread and the error does not change significantly between the period used to determine the relation and the period for which the predictions are made. The method was tested on both equally weighted and learning weighted multi-model ensembles of decadal climate predictions \citep{Smith2007,meehl_decadal_2009,meehl_decadal_2014}. It is important to note that decadal climate predictions are based on initialized (using interpolated observed conditions) simulations of climate models \citep{Smith2007,meehl_decadal_2009,meehl_decadal_2014}. The relatively short duration of the simulation (and prediction) period justifies our assumption that the characteristics determined by the results of the initial simulation period will remain valid during the remaining simulation period. We found that the method results in more reliable predictions than those generated by methods that rely on more assumptions.

\section{Methods}
\subsection{Ensemble of climate model predictions}
The ensemble considered here includes simulations of eight climate models from the CMIP5 decadal experiments (see Table \ref{Models_table}).
We arbitrarily chose one realization of each model. All the models were interpolated to the same spatial resolution of the reanalysis data.
In order to ensure a long enough prediction time series, we focused on the 30-year simulations for the period 1981 to 2011.
The monthly means of the surface temperature were considered as the variable of interest throughout the paper.
The results of a similar analysis for the surface zonal wind appear in the supplementary information.
The first 10 years of the simulations were used to apply the different methods considered here and to tune their parameters while the last 20 years of the predictions were used to validate and measure the performance of the methods.
\begin{table*}
\caption{\label{Models_table} Model Availability}
\centering
\begin{tabular}{ p{3cm}  p{3cm}  p{6cm} p{2.5cm}}
  \hline
  Institute ID  & Model Name & Modeling Center (or Group) & Grid (lat X lon)\\
   \hline
  BCC  &   BCC-CSM1.1  &   Beijing Climate Center, China Meteorological Administration & 64 X 128\\
    CCCma & CanCM4 & Canadian Centre for Climate Modelling and Analysis & 64 X 128\\
    CNRM-CERFACS & CNRM-CM5 & Centre National de Recherches Meteorologiques / Centre Europeen de Recherche et Formation Avancees en Calcul Scientifique & 128 X 256\\
    LASG-IAP  & FGOALS-s2  & LASG, Institute of Atmospheric Physics, Chinese Academy of Sciences & 108 X 128\\
    IPSL* & IPSL-CM5A-LR  & Institute Pierre-Simon Laplace & 96 X 96\\
    MIROC & MIROC5 ~~~~~~~~~~~~~~~ MIROC4h &  Atmosphere and Ocean Research Institute (The University of Tokyo), National Institute for Environmental Studies, and Japan Agency for Marine-Earth Science and Technology & 128 X 256

     320 X 640\\
      MRI & MRI-CGCM3  & Meteorological Research Institute & 160 X 320\\
 \hline
 \end{tabular}
\begin{flushleft}
{\footnotesize * not available for surface zonal wind}
\end{flushleft}
\end{table*}
The NCEP/NCAR reanalysis \citep{kalnay_ncep/ncar_1996} data (with a spatial resolution of $2.5^{\circ}X2.5^{\circ}$) were considered here as observations for the purpose of tuning and validating the methods.
The ERA-interim reanalysis \citep{1321008426928} data were also used, with two different spatial resolutions ($0.75^{\circ}X0.75^{\circ}$ and $2.5^{\circ}X2.5^{\circ}$). We found that the performances of the different methods were not sensitive to the reanalysis data set that we used or to its spatial resolution (obviously, the latter statement is limited to the spatial resolutions that we tested and may not be generally true; see the supplementary information for the results using the ERA-interim reanalysis).
\subsection{Weighted ensemble}
Predictions that are based on an ensemble of models are, in general, a weighted average of the ensemble members. The prediction for time $t$ was defined as:
\begin{equation}
p_t \equiv \sum_{E=1}^N w_{E,t} f_{E,t}.
\end{equation}
Here, $f_{E,t}$ is the prediction of model $E$ for time $t$. $E\in[1..N]$ and $t\in[1..n]$ where $N$ is the number of models in the ensemble and $n$ is the number of time points for which predictions are made. For simplicity, we restricted our attention to the case in which all the models provide forecasts for the same period with equal time steps.
If there is no a priori knowledge, the weight of all members may be equal. Otherwise, the weight of the models may be based either on previous knowledge or on the past performances of the models.
In order to demonstrate the generality of the methods presented here, we used both an equally weighted ensemble and a weighted ensemble for which the weights were generated by a learning algorithm.
The learning algorithm that we used is the Exponentiated Gradient Average (EGA)\citep{kivinen_exponentiated_1997,cesa2006prediction} that was shown to outperform the equally weighted ensemble in decadal climate predictions \citep{Strobach_2015,Strobach_2016}.
The learning algorithm used the first 10 years of the simulations to assign weights to the different models, and those weights were then used to generate the predictions for the following 20 years of the simulations. In this case, the weights during the last 20 years were independent of time. It is important to note that the weights were assigned independently for each grid cell.
\subsection{Uncertainty estimation}
One of the main advantages of using an ensemble of models is the fact that one can not only obtain better predictions but also an estimate for the uncertainty range.
The uncertainty range is usually defined as the range within which there is probability $c$ to find the variable.
We defined the variance of the ensemble prediction at time $t$ as:
\begin{equation}\label{eq:var}
\sigma_t^2 \equiv \sum_{E=1}^N w_E (f_{E,t}-p_t)^2.
\end{equation}
In the equation above and in what follows, we assumed that the weights are independent of time. Because the predictions of different models are not necessarily i.i.d., we still defined the variance of the ensemble according to Eq. \ref{eq:var}, without any prefactor due to the weights of the ensemble members.

If we assume that the model predictions at each time step are i.i.d. random variables, which is obviously not the case for a multi-model ensemble of climate predictions, then for many models, their average is normally distributed according to the central limit theorem.
For a standard normal distribution (with a mean equal to zero and an STD equal to one), the interval that includes the variable probability $c$ can be derived from the \textit{probit} function:
\begin{equation}
\delta_G =\sqrt{2} \mathit{erf}^{-1}(c).
\end{equation}
The error function was defined as:
\begin{equation}
\mathit{erf}(s) \equiv \frac{1}{\sqrt{\pi}} \int_{-s}^{s} e^{-x^2} dx.
\end{equation}
Following the above assumptions, we defined the confidence interval for probability $c$ and for time $t$ as:
\begin{equation}
Pr \lbrace \left(p_t - \delta_G \sigma_t\right) \leq y_t \leq \left(p_t + \delta_G \sigma_t\right) \rbrace = c.
\end{equation}
$y_t$ denotes the value of the variable at time $t$. In what follows, we refer to this definition of the confidence interval as the Gaussian estimation.

In most ensembles of climate models, the basic assumptions required for the central limit theorem to hold are not valid and are rarely verified. The number of models in the ensemble is often limited, and the model predictions are not independent and are not necessarily identically distributed. Therefore, the above method for estimating the uncertainty range is not expected to perform well. In order to overcome the abovementioned problems, we suggest here two alternative methods for estimating the uncertainties.

The first method is based on the same assumption of the normal distribution of the ensemble average. Using the ratio between the root mean squared error (RMSE) and the STD of $p_t$ during the learning period (i.e., the first 10 years), a multiplication constant was found to derive a better estimate of the uncertainty range. For an unbiased estimator of $p_t$, the RMSE should be equal to the STD. Therefore, the uncertainty range was derived from a corrected STD, which was defined as $\delta_R\equiv\delta\gamma$ and the correction factor $\gamma$ was defined as:
\begin{equation}
\gamma \equiv \sqrt{\frac{\sum_{t=1}^n (p_t-y_t)^2}{ \sum_{t=1}^n {\sigma_t}^2}} .
\end{equation}
Using this estimation for the uncertainty range, we obtained
\begin{equation}
Pr \lbrace \left(p_t - \delta_R \sigma_t\right) \leq y_t \leq \left(p_t + \delta_R \sigma_t\right) \rbrace = c.
\end{equation}
In what follows, we refer to this uncertainty estimation method as the RMSE-corrected method.
Note that this method is based on the assumption that the ratio between the squared error and the variance does not change much in time and, therefore, can be replaced by its temporal average.

The second method that we introduce here further relaxes the assumption that the prediction is normally distributed by considering an asymmetric uncertainty range around the average.
In order to define the asymmetric range, we calculated two coefficients: $\gamma_u$ and $\gamma_d$. $\gamma_u$ is a coefficient multiplying the STD of the ensemble in order to set the upper limit of the uncertainty range such that there is a probability of $(1-c)/2$ to find the variable above this upper limit. Similarly, $\gamma_d$ sets the lower limit of the range such that there is a probability $(1-c)/2$ to find the variable below this lower limit. Both coefficients were determined by choosing the minimal range such that during the learning period, the probability of finding the variable between the lower and upper limits is $c$. Mathematically, the coefficients were defined as:
\begin{equation}
\gamma_u=\inf \lbrace \gamma_u \in \Re : \sum_{t=1}^n \Theta((p_t+\gamma_u \sigma_t)-y_t) \geq \frac{1+c}{2} \rbrace,
\end{equation}
and
\begin{equation}
\gamma_d=\inf \lbrace \gamma_b \in \Re : \sum_{t=1}^n \Theta(y_t-(p_t-\gamma_d \sigma_t)) \geq \frac{1+c}{2} \rbrace.
\end{equation}
Using these definitions, the uncertainty range is given by:
\begin{equation}
Pr \lbrace p_t - \gamma_d \sigma_t \leq y_t \leq p_t + \gamma_u \sigma_t \rbrace \approx c
\end{equation}
The non-equality sign here stems from the fact that the probability can change in steps of $1/n$ ($n$ is the number of points in the time series used for calculating $\gamma_u$ and $\gamma_d$). Practically, it implies that for a fine resolution in setting the confidence level $c$, the length of the learning time series, $n$, should be large. In what follows, we refer to this estimation of the uncertainty range as the asymmetric range method. Note that this method defines the range such that upper and lower tails of the distribution have an equal weight of $(1-c)/2$.

It is important to note that all three methods described above base the estimation of the uncertainty range on the STD of the ensemble. Moreover, they all assume that the ratio between the uncertainty range and the STD does not vary much in time. The Gaussian method relies on the assumption of normally distributed model predictions at each time step, the RMSE-corrected method relaxes this assumption and attempts to define a range that is centered on the ensemble prediction, minimizing the prediction bias, and the asymmetric range method attempts to determine a minimal range such that the probability of finding the variable within this range is $c$, regardless of the details of the model predictions distribution.
Figures 1-5 of the supplementary information present the predictions and their associated uncertainty ranges as estimated by the three methods described above.

The performances of the different methods are measured by the relation between the desired confidence level, $c$, and the actual fraction of the validation time series that fell within the estimated uncertainty range. This fraction is defined as:
\begin{equation}
 F_q=\frac{1}{n_v}\sum_{t=1}^{n_v}\Theta\left(y_t-(p_t-\delta_q\sigma_t)\right)\Theta\left((p_t+\delta_q\sigma_t)-y_t\right),
\end{equation}
where $q$ takes the values $G$ and $R$ for the Gaussian and RMSE-corrected methods, respectively. $n_v$ is the number of time points in the validation series, and the index $t$ measures the time points from the beginning of the validation time series (the period following the learning).
For the asymmetric method, the fraction is defined as:
\begin{equation}
 F_A=\frac{1}{n_v}\sum_{t=1}^{n_v}\Theta\left(y_t-(p_t-\delta_d\sigma_t)\right)\Theta\left((p_t+\delta_u\sigma_t)-y_t\right).
\end{equation}

\section{Results}
We first investigated the ratio between the uncertainty ranges of the surface temperature estimated by the RMSE-corrected and the asymmetric range methods ($ \Delta_R\equiv 2  \delta_R  \sigma_t$ and $\Delta_A\equiv\left(\gamma_u+\gamma_d\right) \sigma_t$, respectively) and those estimated by the Gaussian method ($\Delta_G\equiv2 \delta_G \sigma_t$) for a confidence level of $c=0.9$.
\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=39pc]{Fig1.pdf}
\caption{\label{fig:tas_ap=1_gamma-crop} The log of the ratio between the uncertainty ranges of the surface temperature with a confidence level of $0.9$ as estimated by the RMSE-corrected, asymmetric and Gaussian methods. The two left panels show the ratio for the predictions of the equally weighted ensemble, and the two right panels show it for the ensemble weighted by the EGA forecaster. The two upper panels show $\mathit{ln}\left(\Delta_R/\Delta_G\right)$, and the two lower panels show $\mathit{ln}\left(\Delta_A/\Delta_G\right)$. Above each panel, the global mean and the spatial STD are denoted. }
\end{center}
\end{figure*}
Figure \ref{fig:tas_ap=1_gamma-crop} shows the spatial distribution of the natural logarithm of the ratio. The log scale was chosen to emphasize values above and below one (zero in the log scale). Positive values (the ratio is larger than one) correspond to uncertainty ranges (estimated by the RMSE-corrected and the asymmetric range methods) that are larger than the range estimated by the Gaussian method. Similarly, negative values (the ratio is smaller than one) correspond to smaller uncertainty ranges than the range estimated by the Gaussian method. The two left panels show the ratio for the equally weighted ensemble, and the two right panels show it for an ensemble weighted according to the EGA learning algorithm.
The two upper panels show $\mathit{ln}\left(\Delta_R/\Delta_G\right)$, and the two lower panels show $\mathit{ln}\left(\Delta_A/\Delta_G\right)$.
Above each panel, the global average (on the left) of the ratio and its STD (on the right) are denoted.

The upper panels of Fig. \ref{fig:tas_ap=1_gamma-crop} show that, in general, our ensemble forecast is overconfident, i.e., the spread of the ensemble predictions is smaller than the typical error, and it becomes even more overconfident when the weighting is done according to the EGA learning algorithm. The spatial variability of the ratio is not very large.

The performances of the different methods of estimating the uncertainty range can be tested by comparing the desired confidence, $c$, to the fraction of the observations within the predicted range, $F_q$.
$F_q<c$ implies overconfidence (i.e., too narrow a range) and vice versa.
\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=39pc]{Fig2.pdf}
\caption{\label{fig:tas_ap=1_rel_diagram_global9-crop} The spatial distribution of the difference between the fraction of the observations that were outside the predicted range of the $c=0.9$ confidence level and the predicted $0.1$ fraction ($c-F_q$), for the surface temperature. The three left panels show the difference for the equally weighted ensemble, and the three right panels show it for the EGA forecaster.
The two upper panels show ($c-F_G$), the two center panels show ($c-F_R$), and the two lower panels show ($c-F_A$). The spatial average and the STD of the difference are denoted above each panel.
}
\end{center}
\end{figure*}
Figure \ref{fig:tas_ap=1_rel_diagram_global9-crop} shows the difference between the fraction of the observations that were outside the predicted range of the $0.9$ confidence level and the predicted $0.1$ fraction ($c-F_q$).
The three left panels show the difference for the equally weighted ensemble, and the three right panels show it for the EGA forecaster.
The two upper panels show $c-F_G$, the two center panels show $c-F_R$, and the two lower panels show $c-F_A$.
Above each panel, we denote the spatial average (on the left) and the STD (on the right) of the difference.
The figure shows that in most regions, $\Delta_G$ is too small ($c>F_G$), which justifies the larger values of $\Delta_R$ and $\Delta_A$. The underestimation of $\Delta$ by $\Delta_G$ is more severe for the EGA than for the equally weighted ensemble.

Figure \ref{fig:tas_ap=1_rel_diagram_global9-crop} demonstrates that the RMSE-corrected and asymmetric methods outperform the Gaussian method. Moreover, it shows that the asymmetric method outperforms the RMSE-corrected method. One can also notice that the methods work better for the equally weighted ensemble than for the EGA forecaster. This better performance is due to the fact that the ensemble STD, $\sigma_t$, of the EGA is defined by weights that minimize the forecast error and not by weights that maximize the forecast reliability.
In the supplementary information, we provide similar results for the surface zonal wind (Fig. 7) and for the temporal evolution of $F_q$ (Figs. 8-9). The results accord with the results presented here for the surface temperature and with the basic assumptions of the methods for estimating the uncertainty range.
\begin{figure*}[!ht]
\begin{center}
\includegraphics[width=27pc]{Fig3.pdf}
\caption{\label{fig:tas_ap=1_rel_diagram-crop}  Reliability diagrams for the different estimation methods of the surface temperature uncertainty range.
The observed frequency represents the spatial average over all the grid cells. The left panels show the reliability for the equally weighted ensemble, and the right panels show it for the EGA forecaster.
The upper panels show the reliability of the Gaussian method, the center panels show the reliability of the RMSE-corrected method, and the lower panels show the reliability of the asymmetric method.}
\end{center}
\end{figure*}
These methods can be used for different confidence levels. Tests of the method performances that span the whole range  of confidence levels, $c=[0,1]$, provide more information about the relation between the estimated and observed distributions of the variable.
Figure \ref{fig:tas_ap=1_rel_diagram-crop} shows the observed frequency versus the expected one (reliability diagrams) for different confidence levels.
The observed frequencies represent the spatial average over the whole globe.
The left (right) panels show the reliability diagrams for the equally weighted ensemble (EGA forecaster).
The upper panels show the reliability for the Gaussian method, the middle panels show the reliability for the RMSE-corrected method, and the lower panels show the reliability for the asymmetric method.

Figure \ref{fig:tas_ap=1_rel_diagram-crop} shows that the RMSE-corrected and asymmetric methods perform much better than the Gaussian method. The improvement is more pronounced for the EGA forecaster. One can also observe that the asymmetric method is more reliable for higher confidence levels.
Similar results were obtained for the surface zonal wind. However, we found that in most cases, the uncertainty range for the surface zonal wind was underestimated by all the methods (see Fig. 8 of the supplementary information).

\section{Summary and Discussion}

A new method for estimating the uncertainty of ensemble-based climate predictions was suggested.
This method is based on learning the relations between the prediction errors and their spread in the past and using these relations in order to estimate the uncertainties of future predictions for which only the spread is known.
The method does not rely on any assumptions regarding the distribution of the ensemble member predictions or that they are i.i.d.
It also has the advantage of estimating separately the upper and lower uncertainties.
The inherent assumption of the method is that the relation between the spread of the ensemble predictions and the errors does not significantly change during the period spanned by the predictions. Moreover, it assumes that this relation has relatively small fluctuations, which allows it to be considered as independent of time.

The excellent performance of the method in estimating the uncertainties during a validation period, based on the relations found during an earlier learning period, suggests that these assumptions are valid for a multi-model ensemble of decadal climate predictions. Moreover, it was shown that the method performs well for both the equally weighted ensemble and the ensemble weighted according to a learning algorithm applied during the training period. The performance of the asymmetric method was better than those of the RMSE-corrected method and the method based on the assumption of a Gaussian distribution of the ensemble member predictions.
It was also found that, in general, estimations of the uncertainties that are based on the ensemble spread without correction result in overconfident predictions (therefore, the correction factors are larger than one). The overconfidence is even stronger for the weighted ensemble. It is possible that weighting schemes that are based on the reliability of the predictions and not just on their error may result in a spread that is closer to the errors; however, a deeper discussion of this issue is beyond the scope of this paper.
The method suggested here is not limited to estimating one measure of the uncertainties (for example, the variance); applying the method with different confidence intervals was shown to provide an excellent estimation of the whole PDF of the variables.
The asymmetric range method is not limited to decadal climate predictions and is likely to be useful for climate predictions of shorter time scales or for weather predictions.
In particular, the method may be useful for estimating the uncertainties for data assimilation in numerical weather predictions \citep{Houtekamer1996,Kalnay2003,Warner2010}.

\clearpage

\acknowledgments
The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant number [293825].

We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for the CMIP, and we thank the climate modeling groups (listed in Table \ref{Models_table} of this paper) for producing and making available their model outputs. For the CMIP, the U.S. Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and leads the development of software infrastructure in partnership with the Global Organization for Earth System Science Portals.

\bibliographystyle{unsrt}


\section{Introduction}
\label{sec:intro}
BL Lac objects constitute the most elusive class of radio-loud active galactic nuclei \citep[AGNs;][]{urry95}. Hosted in elliptical galaxies \citep[see e.g.,][]{falomo14}, their emission extends from radio frequencies up to the highest TeV energies and is characterized by rapid and large-amplitude flux variability, flat radio spectra, even below $\sim$1GHz \citep[see e.g.,][]{massaro13a,massaro13b}, peculiar infrared (IR) colors \citep{paper1,paper2} and significant and variable optical polarization (i.e., $P\sim10\%$) \citep[see e.g.,][and references therein]{andruchow05,smith07}.

BL Lacs show featureless optical spectra, sometimes having weak emission lines of rest-frame equivalent width ($EW_r$) lower than 5\AA\ or absorption lines due to their host galaxies, both superimposed to a blue continuum \citep[see e.g.,][]{laurent99}. BL Lac objects are one of the two subclasses of the radio-loud AGNs known as blazars \citep[see e.g.,][for a recent review]{massaro09}. The other subclass is constituted by the flat spectrum radio quasars (FSRQs), characterized by properties similar to those of BL Lacs but having a quasar-like optical spectra \citep[i.e., with broad emission lines; see e.g.,][]{veron00} and radio spectral indices\footnote{Spectral indices, $\alpha$, are defined by flux density, S$_{\nu}\propto\nu^{-\alpha}$} lower than 0.5 between $\sim$300 MHz and a few GHz \citep[see][and references therein]{massaro14a}.

Over the entire electromagnetic spectrum both BL Lacs and FSRQs feature a typical double-bumped spectral energy distribution (SEDs) with the first component peaking between the IR and the X-rays and the second one in the gamma-rays \citep[see e.g.,][]{padovani95,paggi09,abdo10a}. Mid-IR colors clearly helps to determine the energy peak location of the low energy component \citep{paper1}. BL Lac SEDs and their broad-band features were described in terms of emission arising from particles accelerated in relativistic jets directly pointing towards the observer \citep{blandford78}.

BL Lacs are traditionally discovered using radio, X-ray surveys and/or their combination \citep[see e.g.,][]{stickel91,sambruna96,giommi05}, since in the optical energy range the lack of spectral features combined with their flux variability makes a color-color selection extremely challenging \citep{jannuzi93,londish02,collinge05,plotkin08}.  However, recent observations to search for optical polarization in BL Lacs selected using deep optical survey as the Sloan Digital Sky Survey \citep[see e.g., SDSS DR12][and references therein]{alam15} and the Two-Degree Field Quasar Redshift Survey \citep[2QZ][and references therein]{smith05} validated this strategy of investigating proper motion and broad-band color properties to minimize contamination by featureless Galactic objects while searching for quasi-featureless spectra \citep{smith07}.

At the current moment, the most efficient survey for discovering new BL Lacs is the {\it Fermi}\ gamma-ray one \citep{abdo10b,ackermann11a,nolan12}, in particular when combining gamma-ray observations with mid-IR colors analyses \citep{massaro12a,massaro12b} and follow up optical spectroscopic observations \citep{paggi13,shaw13,massaro15,refined}. This occurs because BL Lacs are the largest known population emitting at MeV-GeV energies, constituting more than 30\% of the associated sources in the latest release of the {\it Fermi}\ catalog \citep[The {\it Fermi}\ - Large Area Telescope Third Source Catalog, 3FGL][]{acero15}. Moreover, recent optical spectroscopic campaigns aimed to confirm the nature of blazar candidates of uncertain type listed in the 3FGL showed that a significant fraction of those observed to date are also BL Lacs \citep[see e.g.][]{crespo16a,crespo16b,crespo16c}.

However, nowadays, all the known {\it Fermi}\ BL Lacs have radio counterparts and radio surveys are therefore crucial to associate $\gamma$-ray objects with their counterparts. Thus, one of the major open issues for this rare AGN class is proving the existence of radio-weak BL Lac objects (RWBLs). We indicate them as radio-weak sources instead of using the term {\it ``radio-quiet''}, generally adopted for the AGN unification scenario. This designation is more useful because radio follow up observations could reveal some radio emission. However, the current lack of counterparts in the major radio surveys available to date place them in the tail of the distribution of the IR-to-radio flux ratios for the whole BL Lac population.

At the current moment only a handful of sources have been claimed to be RWBLs on the basis of their broad band properties as featureless optical spectra combined with X-ray emission and and their lying within the positional uncertainty regions of unidentified/unassociated gamma-ray sources \citep[UGSs;][]{paggi13,ricci15,landoni16}.

Their existence will have a crucial impact (i) on the estimates of their source counts and luminosity functions \citep{ajello12,ajello14}, a task already complicated by the challenge of measuring their redshifts \citep{falomo14}, (ii) on the gamma-ray association methods all currently based on radio surveys \citep[see e.g,][and references therein]{ackermann15}, if proved that RWBLs could be associated to {\it Fermi}\ sources, as well as (iii) on the unification scenario of AGNs \citep{urry95}, to name a few implications. An accurate accounting is crucial to compute correctly their contribution to the extragalactic gamma-ray background for which, together with FSRQs \citep{ajello12} and radio galaxies \citep{inoue11,dimauro14,lobes} they provide the most significant ones. Additional implications are related to jet physics. Radio emission in extragalactic radio-loud AGNs permits us to investigate the ``large-scale'' (i.e. from pc to Mpc size) evolution of relativistic jets and their interactions with the environments \citep[see e.g.,][and references therein]{liuzzo13,nulsen13,fabian12}. The existence of RWBLs could be crucial to address different phases of jet formation in cases when relativistic jets are not able to form extended radio structures \citep{blandford79,begelman84}, as it seems from low radio frequency (i.e., below $\sim$1 GHz) spectral studies \citep[see][and references therein]{massaro13b,ugs6lowfrq,mwabl}.

There are two major problems to solve before claiming the existence of RWBLs, even when featureless optical spectra are available. They can be confused with or misinterpreted as hot white dwarfs \citep[WDs; see e.g.][]{kepler15} and/or weak emission line quasars \citep[WELQs;][and references therein]{ds9,plotkin10} when optical spectroscopic data are available. { Optical polarimetric observations could shed light on their nature. However, unfortunately, optical polarization in BL Lacs can be variable and it is not always observed thus not representing the most definitive stamp to confirm their nature \citep[see e.g.,][]{londish04}.}

In the present paper we explore the use mid-IR colors as a diagnostic tool to disentangle between BL Lacs, WDs and WELQs. This new tool could potentially open a new window to search for RWBLs. This investigation is motivated by the recent discovery that {\it Fermi}\ BL Lacs show peculiar mid-IR colors \citep[see also,][and references therein]{gir}, result made possible by observations performed with the {\it WISE}\ all-sky survey \citep{wright10}, covering a region well isolated from that of other Galactic and extragalactic objects. The main explanation underlying this distinction is the BL Lac non-thermal IR emission \citep[see e.g.,][]{paper2,ugs1}.

This manuscript is organized as follows. In \S~\ref{sec:samples} we present all the samples used to perform our mid-IR color-color analysis. In \S~\ref{sec:analysis} we compare the mid-IR colors of BL Lacs with those of WELQs and WDs, while in \S~\ref{sec:rwbls} we search for RWBLs. \S~\ref{sec:variability} is dedicated to the analysis of the IR variability for the best candidates selected on the basis of the mid-IR colors ad the optical spectra. Then in \S~\ref{sec:swift} we also investigate the {\it SWIFT}\ observations available for these candidates while in \S~\ref{sec:seds} we discuss on the SEDs. Summary and conclusions are given in \S~\ref{sec:summary}

Unless stated otherwise, for our numerical results, we use cgs units. {\it WISE}\ magnitudes used here are in the Vega system and are not corrected for the Galactic extinction. As shown in our previous analyses, such corrections mostly affect the magnitude at 3.4$\mu$, for sources lying at low Galactic latitudes and it ranges only between 2\% and 5\% \citep[see e.g.][]{wibrals}, thus being negligible. In addition this allows us to compare plots in the present analysis with those of previous investigations \citep{ugs2}.

\section{Sample selection}
\label{sec:samples}
To perform our analysis we built three main samples. The first two composed of WDs and WELQs, respectively; both used in comparison with the third one of {\it Fermi}\ blazars. We searched for mid-IR counterparts of all samples in the latest release of the {\it WISE}\ catalog\footnote{http://wise2.ipac.caltech.edu/docs/release/allwise/}. All our final samples of WDs, WELQs and gamma-ray blazars only include sources having a {\it WISE}\ counterpart within 3\arcsec.3 that is detected at 3.4$\mu$m, 4.6$\mu$m and 12$\mu$m so allowing us to build at least one of the mid-IR color-color plots used in our previous blazar analyses \citep{wibrals,gir}. The choice of 3\arcsec.3 association radius is the same adopted in our previous investigations \citep{ugs1} and it was computed on the basis of Montecarlo procedure corresponding to a probability of having spurious associations of the order of $\sim$1\% \citep[see also][for additional details]{wibrals}.

\subsection{BL Lac objects and Flat Spectrum Radio Quasars}
The first sample used in our analysis lists all the blazars, both BL Lacs and FSRQs, included in the {\it Roma-BZCAT}\footnote{http://www.asdc.asi.it/bzcat/} \citep{massaro15c} and having a counterpart in the latest release of the {\it Fermi}\ catalog \citep{acero15}. This sample has been recently used to perform the analysis on the {\it IR--gamma-ray connection} \citep{gir}. According to the nomenclature proposed in the {\it Roma-BZCAT}\ hereinafter we will refer to BL Lac objects as BZBs and to the FSRQs as BZQs. We choose this sample because it has the largest fraction of BZBs with a mid-IR counterpart in {\it WISE}\ corresponding to 603 out of 610 (i.e., $\sim$99\%) within our association radius; this fraction is similar to that occurring for the BZQs with 419 out of 426 (i.e., $\sim$98\%).

\subsection{White Dwarfs}
Hot WDs could be misclassified as BZBs because they could show almost featureless optical spectra and X-ray radiation \citep[see e.g.,][]{heise85,bilikova10}. Thus in the present analysis we want to show that the mid-IR colors of the {\it WISE}\ counterparts of WDs is different from that of BZBs. We expect that for WDs mid-IR emission arises from the presence of a companion brown dwarf or the presence of a dusty disk\footnote{http://www.stsci.edu/$\sim$debes/wired.html}, both expected to have thermal IR colors \citep{debes11} while those of BZBs are dominated by the non-thermal synchrotron radiation of particles accelerated in their jets.

To achieve this goal we considered Sloan Digital Sky Survey DR7 White Dwarf Catalog (SDSSDR7WD) \citep{kleiman13}. This new catalog of spectroscopically confirmed WDs from the Sloan Digital Sky Survey (SDSS) Data Release 7 (DR7) listing more than 19000 stars. Crossmatching the SDSSDR7WD catalog with the latest release of the {\it WISE}\ catalog we found 4125 have a unique mid-IR potential counterpart within 3\arcsec.3. However, only for 255 out of 4125 sources their mid-IR counterparts were detected by {\it WISE}\ in all three of its short-wavelength bands, and so we only do a comparison with this subset.

We assume that mid-IR correspondences found for the SDSSDR7WD are the real counterparts of the WDs. However in the worst scenario if none of the WDs listed in the SDSSDR7WD has mid-IR emission this will automatically confirm that their mid-IR spectral behavior is different from that of {\it Fermi}\ BZBs known to date. We also computed the chance probability of having spurious associations within 3\arcsec.3 when crossmatching the All{\it WISE}\ source catalog with the SDSSDR7WD; the result is less than 4\% \citep[see e.g.,][for more details on the method adopted for the pchance calculation]{ugs1}.

To verify the reliability of our choice for the WD sample we also used { the latest release of the McCook \& Sion White Dwarf Catalog\footnote{http://www.astronomy.villanova.edu/WDCatalog/index.html} available to date} \citep[see][for details]{mccook99}. This catalog has been also used to search for dusty disks around WDs \citep{hoard13} using {\it WISE}\ observations. However in this sample the number of sources detected at 3.4$\mu$m, 4.6$\mu$m and 12$\mu$m for which is possible to built the mid-IR colors is only 91. The chance probability of spurious associations between the McCook and Sion WD catalog and the All{\it WISE}\ survey is lower than 1\%.

\subsection{Weak emission line quasars}
In the last two decades an intriguing and also elusive subclass of radio-quiet quasars was discovered: the weak emission line quasars (WELQs), featuring exceptionally weak or completely missing broad emission lines in the ultraviolet (UV) rest-frame energy range \citep[see e.g,][]{collinge05,shemmer06,ds9}.

The original criterion to define WELQs with respect to other members of the same class (i.e., normal quasars) is based on the $EW_r$ distribution of the Ly$\alpha$ $\lambda$1216+N V $\lambda$1240 blend for a sample of quasars in the Sloan Digital Sky Survey \citep[SDSS;][]{} at redshift $z>$3, being lower than 15$\AA$ as occurs for the $>3\sigma$ weak tail of the $EW_r[Ly\alpha+NV]$ distribution. Using this definition, a first sample of $\sim$70 WELQs was defined \citep{ds9}. Since at low $z$ it is not always possible to apply the original criterion on the $EW_r[Ly\alpha+NV]$ because these emission lines lie outside the optical energy range covered by the SDSS spectrograph an equivalent criterion was defined selecting WELQs on the basis of the $EW$ of Mg II $\lambda$2800, C III] $\lambda$1909, and/ or C IV $\lambda$1549 \citep[see e.g.,][]{m14}.

When WELQs were discovered the possibility of being classified BZBs was explored and discarded \citep[see e.g.][for a recent review]{plotkin15}, however a comparison of their mid-IR properties with those of confirmed BZBs has not yet been completely performed. Here we selected a sample of WELQs to show that their mid-IR colors are generally similar to those of normal quasars but not to BZBs. Additional analyses based on the SEDs and/or broadband spectral indices (i.e., radio-to-optical $\alpha_{ro}$ and optical-to-Xray $\alpha_{ox}$) were also used to distinguish between WELQs and BL Lacs \citep{plotkin10,lane11,wu12} when multifrequency observations are available. Thus in the following sections we also investigate the broad band SEDs of RWBL candidates selected on the basis of the combination of the mid-IR colors and the optical spectra.

We considered two catalogs of WELQs, spectroscopically selected from the SDSS. The first catalog lists 73 WELQs in the redshift range between 3.03 and 5.09 \citep[][hereinafter DS9 subsample]{ds9} while the second lists 46 sources with 0.614$<z<$3.351 \citep[][hereinafter M14 subsample]{m14}. We excluded 6 sources out of these 119 WELQs, namely: SDSS J074451.36+292006.0 and SDSS J092145.37+233548.1 from the M14 subsample and SDSS J141318.86+450522.9 from DS9 subsample because they are all listed in the {\it Roma-BZCAT}\ and SDSS J104831.29+211552.2, SDSS J115326.70+361726.3, SDSS J144204.03+132916.0 again from the M14 subsample since they show a flat radio spectrum that is widely interpreted as jet signature so suggesting that their quasi-featureless optical continuum marks non-thermal synchrotron emission.

Only 90 out of the remaining 113 sources arising from the combination of the DS9 and the M14 subsamples have a counterpart in the All{\it WISE}\ catalog. These 90 sources thus constitute our WELQ reference sample to carry out the comparison with BZB mid-IR colors. We estimated the chance probability of spurious associations between the SDSSDR7WD catalog and the {\it WISE}\ survey as being lower than 1\%.

\section{BZBs, WDs and WELQs at mid-IR wavelengths}
\label{sec:analysis}
In Fig.~\ref{fig:strip} we show the mid-IR colors built with the in the 3.4$\mu$m, 4.6$\mu$m and 12$\mu$m magnitudes of the BZBs and the BZQs in our {\it Fermi}\ blazar sample in comparison with generic infrared sources selected at high Galactic latitudes ($|b|>$50\degr), to highlight their peculiar spectral behavior. Since we know that BZQs generally lie in the region occupied by ``normal'' quasars, at the top end of the so called ``{\it WISE}\ Gamma-ray Strip'' \citep{paper1}, to highlight the mid-IR color region where BZBs reside, we will only consider RWBL candidates those {\it WISE}\ sources lying in the dashed polygon shown in Fig.~\ref{fig:strip} where the contamination of the BZB region by BZQ is lower than 5\%, lacking a radio counterpart. It is worth noting that the selected region also includes TeV BL Lacs.
\begin{figure}[!h]
\includegraphics[height=8.cm,width=8.5cm,angle=0]{./strip.pdf}
\caption{The mid-IR color-color plot built with the {\it WISE}\ magnitudes at 3.4$\mu$m, 4.6$\mu$m and 12$\mu$m. BZBs and BZQs are displayed as cyan circles and orange squares respectively, thus corresponding to the ``{\it WISE}\ Gamma-ray Strip'' \citep{paper1}. Generic mid-IR sources (grey circles) of about 3000 sources selected at high Galactic latitudes are also show. The separation between the mid-IR colors of the {\it Fermi}\ blazars and that of other sources randomly chosen appears evident. The presence of BZBs in the region of the ``{\it WISE}\ Gamma-ray Strip'' occupied by the BZQs is only limited to a handful of sources \citep[see e.g.,][for more details]{ugs1,wibrals}. The dashed polygon marks the region used to select potential RWBLs.
}
\label{fig:strip}
\end{figure}

Fig.~\ref{fig:wds} shows the comparison between {\it WISE}\ sources associated with WDs and our blazar sample. It is evident that only 4 out of the 255 sources in the SDSSDR7WD catalog have mid-IR colors consistent with the main BZB region of the ``{\it WISE}\ Gamma-ray Strip'', while the remaining $\sim$99\% of the WDs lie well separated from this region. Three out of four sources are extragalactic, while WISE J122859.87+104032.8 has a clear WD optical spectrum, making it the only WD contaminant of the BZB region in the mid-IR color-color diagram. Moreover WISE J122859.87+104032.8 has a negative value of the $u-r$ a situation that never occurs for the known BZBs \citep{massaro12c,massaro14b}. WISE J172633.50+530300.3 also shows a featureless optical spectrum with a peculiar continuum shape; this could be due to the presence of a nearby source clearly visible in the SDSS image but not in the {\it WISE}\ one. This implies a relatively small contamination of the mid-IR color selection.

Both WISE J001736.91+145101.9 and WISE J121348.83+642520.0 are BL Lacs with a radio counterpart in addition to WISE J121348.83+642520.0 listed in the latest release of the {\it Roma-BZCAT}\ catalog.

{ To complete our comparison between WDs and BZBs we also show the mid-IR color-color plot with the 91 sources selected out of the McCook \& Sion WD catalog in Fig.~\ref{fig:wds}. There are also in this case only two objects that could be confused as BZBs on the basis of their mid-IR colors, namely WISE J085506.13+063904.1 and WISE J115642.21+130603.9, both lying in the footprint of the SDSS. Thus in both cases we checked the optical images using the SDSS Navigate tool\footnote{http://skyserver.sdss.org/dr12/en/tools/chart/navi.aspx}. We found that about 6\arcsec\.6 from the location of WISE J085506.13+063904.1 there is a clear WD, namely SDSS J085506.62+063904.7, for which an optical spectrum is available being most likely the star listed in the McCook \& Sion WD catalog and lacking a {\it WISE}\ counterpart. This mis-association is due to the positional uncertainty of the sources in the McCook \& Sion WD catalog. Thus we conclude that WISE J085506.13+063904.1 is a simple contaminant unrelated to the WD selection. Then for the latter, WISE J115642.21+130603.9 the lack of its optical spectrum combined with the uncertainties of the SDSS photometric observations does not allow us to verify its nature. However, this object is classified as {\it galaxy} thus unlikely to be a WD that can contaminate our RWBL selection.}

The comparison between mid-IR colors of BZBs and WDs prove the usefulness of the mid-IR color-diagram to distinguish the two classes while searching for RWBLs, once optical spectra are available.
\begin{figure*}
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./wds.pdf}
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./wds2.pdf}
\caption{
The comparison between the mid-IR colors of the BZBs and the BZQs as shown in Fig.~\ref{fig:strip} and those of the WDs (black and blue circles) from the SDSSDR7WD (left panel) and the McCook \& Sion WD (right panel) catalogs, respectively. In both comparisons the mid-IR colors of blazars appear well distinct from those of WDs, with the only few exceptions discussed in \S~\ref{sec:analysis} { and marked by the red arrows}. Generic mid-IR sources are also shown in the background in both panels. Uncertainties on the mid-IR colors of BZBs and BZQs are not reported here as in Fig.~\ref{fig:strip} for sake of simplicity since they are shown here for comparison only.
}
\label{fig:wds}
\end{figure*}

Fig.~\ref{fig:wlqs} shows the mid-IR comparison between the WELQs and the BZBs. Here too is clear the neat separation between these two classes. It is also worth noting that sources in the DS9 subsample appear to lie in a different region of the color parameter space with respect to the M14 subsample. The reason underlying this effect is due to their different redshift distributions. We note that there is also a marginal difference between the mid-IR colors of the WELQs in the M14 sample and that of the BZQs but it is less useful to search for RWBLs hidden among this sample.
\begin{figure}
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./wlqs.pdf}
\caption{
The comparison between the mid-IR colors of the BZBs and the BZQs as shown in Fig.~\ref{fig:strip} and those of the WELQs. The WELQs belonging to the M14 subsample are shown in magenta triangles (up) while those of the DS9 subsample in black (down). In both comparisons the mid-IR colors of blazars appear well distinct from those of WELQs, with the only few exceptions discussed in \S~\ref{sec:analysis} { and marked here by the black arrows}. Generic mid-IR sources are also shown in the background in both panels. Uncertainties on the mid-IR colors of BZBs and BZQs are not reported here as in Fig.~\ref{fig:strip} for sake of simplicity since they are shown here for comparison only.
}
\label{fig:wlqs}
\end{figure}

In this case we highlight the fact that two WELQs lie below the threshold marking the BZB region of the color-color plot and additional three objects are close to it. These two sources are WISE J001514.87-103043.4 and WISE J150427.69+543902.2; both showing MgII$\lambda$2798 emission line of rest-frame equivalent width $EW_r$ consistent with the 5\AA\ threshold, usually adopted to classify BL Lacs, within 3$\sigma$, thus being ``border line cases''. The remaining three sources, belonging to the M14 subsample (see Fig.~\ref{fig:rwbls}, where they are shown in magenta), all have a radio counterpart at 1.4 GHz in the Faint Images of the Radio Sky at Twenty centimeter \citep[FIRST;][]{helfand15} catalog, but the lack of additional radio data did not allow us to verify the flatness of their radio spectra, although they are definitively not RWBLs. The location of all five sources with respect to that of the BZBs in the mid-IR color-color plot is shown in Fig.~\ref{fig:rwbls}.

Finally, we conclude that both WDs and WELQs have distinct mid-IR colors with respect to the BL Lac population and [3.4]-[4.6]-[12]$\mu$m color-color diagram can be used as diagnostic tool to confirm the source nature having featureless optical spectra. Sources that appear to be contaminants of the mid-IR selection procedure are ruled out when {\it WISE}\ data are combined with optical spectroscopic observations.

\section{Radio Weak BL Lacs}
\label{sec:rwbls}
In the literature there are five sources tentative classified as RWBLs.

\begin{enumerate}
\item 2QZ J215454.3-305654 \citep{londish02} for which several multifrequency follow up observations have been performed to date was the first. However its nature is still uncertain since it could be a RWBL or a ``lineless'' quasar \citep{londish04}.

\item Paggi et al. (2014), carried out optical spectroscopic observations of blazar candidates, potential UGS counterparts and discovered WISE J064459.38+603131. This mid-IR sources lies within the positional uncertainty region of 2FGL J0644.6+6034.7 (3FGL J0644.6+6035). This object shows a single optical emission line of EW consistent with the 5\AA\ threshold, arbitrarily set in the literature for the BL Lac class definition, similar to the two cases discussed at the end of \S~\ref{sec:analysis}.

\item In 2015 WISE J173052.85-035247.2 was discovered during the optical spectroscopic follow up observations of UGSs \citep{ricci15}. This source could be the potential counterpart of 2FGL J1730.6-0353 (a.k.a.  3FGL J1730.6-0357).

\item More recently Landoni et al. (2016) claimed that also SDSS J004054.65-0915268 could be a tentative RWBL.

\item Finally, Marchesini et al. (2016) discovered the last known candidate of RWBL, WISE J141046.00+740511.2. This source was also found while carrying out spectroscopic observations of UGS listed in the First Fermi-LAT Catalog of Sources above 10 GeV \citep[1FHL;][]{ackermann13} selected on the basis of the {\it SWIFT}\ X-ray observations \citep{landi15}. It has been indicated as potential counterpart of 1FHL J1410.4+7408 (a.k.a. 3FGL J1410.9+7406).
\end{enumerate}
\begin{figure*}[]
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./extend.pdf}
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./rwbls_kno.pdf}
\caption{
Left panel) The mid-IR color-color plot where the positions of the blazars recognized in the WD and the WELQ samples are highlighted with respect to the {\it Fermi}\ blazars. WISE J122859.87+104032.8 and WISE J172633.50+530300.3 are marked with a red arrow since additional optical information allow us to discard them as potential RWBLs. The former has the unusual negative value of the $u-r$ never seen for BL Lacs while the mid-IR colors of the latter are contaminated by the presence of a nearby source visible in the optical image (see \ref{sec:analysis} for more details). Right panel) The positions of the candidates RWBLs in the mid-IR color-color plot are shown in comparison to that of {\it Fermi}\ blazars belonging to the ``{\it WISE}\ Gamma-ray Strip'' and those of generic mid-IR sources (same color code of previous figures). Source details for the objects indicated in both panels are given in \S~\ref{sec:analysis}. Uncertainties on the mid-IR colors of BZBs and BZQs are not reported here as in Fig.~\ref{fig:strip} for sake of simplicity since they are shown here for comparison only.
}
\label{fig:rwbls}
\end{figure*}

In Fig.~\ref{fig:rwbls} we show the positions of the five RWBLs overlay to the ``{\it WISE}\ Gamma-ray Strip'' subregion of BZBs.
SDSS J004054.65-0915268 (a.k.a WISE J004054.63-091526.8) has an upper limit at 12$\mu$m and being a high $z$ source it seems to be more consistent with the location of the DS9 subsample of WELQs rather than the BZB population. A similar situation occurs for 2QZ J215454.3-305654 (a.k.a. J215454.32-305653.2) that only a follow up observations at 12$\mu$m could potentially reveal if the source is consistent with the mid-IR colors of the BZBs. WISE J064459.38+603131.7, being selected on the basis of its mid-IR colors is definitively consistent with the BZB population on the 3.4$\mu$m - 4.6$\mu$m - 12$\mu$m diagram, while WISE J173052.85-035247.2 appears to have such colors more similar to the WELQs listed in the M14 subsample.

The remaining case of WISE J141046.00+740511.2 is the most intriguing one. This object is potentially associated to a 1FHL object and its potential {\it WISE}\ counterpart lies in the region of the ``{\it WISE}\ Gamma-ray Strip'' where all known TeV BL Lacs are located \citep{massaro13d,arsioli15}.
The lack of any radio counterpart in both the NRAO VLA Sky Survey \citep[NVSS;][]{condon98} and the FIRST catalog \citep{helfand15} that cover the field of WISE J141046.00+740511.2 makes this source the most promising candidate to be the first, genuine, RWBL. Nevertheless is it worth emphasizing that the distance between this mid-IR source and the gamma-ray position of 3FGL J1410.9+7406 is $\sim$90\arcsec well below the average of the associated {\it Fermi}\ blazars listed in the {\it Roma-BZCAT}\ corresponding to $\sim$130\arcsec. At the angular separation of $\sim$90\arcsec, the association probabilities of counterparts listed in the latest releases of the {\it Fermi}\ catalogs is typically grater than 99.5\% \citep{ackermann15}.

Applying this mid-IR color-based analysis, in combination with optical spectra, we conclude that WISE J064459.38+603131 and WISE J141046.00+740511.2, lying within the positional uncertainty regions of two distinct UGSs listed in the 3FGL, are the first two RWBLs detected so far. Both of them can also be associable with 1FHL sources thus strengthening our interpretation, since the high energy sky above 10 GeV is mostly dominated by BZBs. In the following we will also show that additional multifrequency observations can strengthen our results.

\section{mid-IR variability}
\label{sec:variability}
BL Lac objects are known to be among the most variable extragalactic sources at all frequencies, showing variations either in spectral shape or in flux/luminosities. This peculiar behavior could in principle mine the mid-IR color selection proposed here and in our previous works. However this situation does not seem to occur.

BL Lacs appear to be variable in the mid-IR band \citep[see e.g.,][]{paper2} although their color variations still reside within the boundaries of the ``{\it WISE}\ Gamma-ray Strip'' \citep[see][for more details]{buson16}. Here we show an example of two BZBs: 5BZBJ1542+6129 and 5BZBJ1959+6508, the latest also detected at TeV energies \citep{aharonian03}, for which more than 50 {\it WISE}\ observations are available, when both sources are detected at least in the first three mid-IR filters at 3.4$\mu$m, 4.6$\mu$m and 12$\mu$m. In Fig.~\ref{fig:moving1} it is clear how both BZBs during their single-epoch {\it WISE}\ observations move in the mid-IR color-color plot in the same region occupied by the BL Lac objects. In addition their colors are also consistent with the boundaries chosen for our RWBL selection procedure.
\begin{figure}[!ht]
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./movingonstrip1.pdf}
\caption{The mid-IR color-color plot built with the {\it WISE}\ magnitudes at 3.4$\mu$m, 4.6$\mu$m and 12$\mu$m. BZBs are displayed as cyan circles and correspond to the ``{\it WISE}\ Gamma-ray Strip'' \citep{paper1}. We show two cluster of points relative to the {\it WISE}\ single epoch observations of the BZBs: 5BZBJ1542+6129 (black circles) and 5BZBJ1959+6508 (red circles). It appears evident that even mid-IR color variation states of these two BL Lac objects are always consistent with the boundaries of the ``{\it WISE}\ Gamma-ray Strip''. The two crosses point to the location of 5BZBJ1542+6129 (magenta) and 5BZBJ1959+6508 (orange) on the mid-IR color diagram when integrated values for the entire {\it WISE}\ all-sky survey are used.
}
\label{fig:moving1}
\end{figure}

According to the {\it WISE}\ catalog one of our best RWBL candidate: the WISE J141046.00+740511.2 is flagged as variable in the mid-IR. This is in agreement with the fact that such emission it is unlikely due to dust emission, as occurs in normal quasars or in white dwarfs \citep{debes11}. Fig.~\ref{fig:lightcurves} shows the {\it WISE}\ magnitude at 3.4$\mu$m, 4.6$\mu$m and 12$\mu$m as function of time for two different epochs; mid-IR variability, not only between the two observing periods, but also in each single set of observations is found. We also present the lightcurve for the second RWBL candidate WISE J064459.38+603131 that clearly shows infrared variability being detected more times at 12$\mu$m during the {\it WISE}\ single-epoch observations with respect to WISE J141046.00+740511.2.
\begin{figure*}
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./lightcurveJ0644.pdf}
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./lightcurveJ1410.pdf}
\caption{The light curve of WISE J064459.38+603131 (left panel) and WISE J141046.00+740511.2 (right panel), candidate RWBLs, tentatively associated to 3FGL and 1FHL {\it Fermi}\ sources. These are {\it WISE}\ single-epoch observations for two different time periods. Mid-IR flux density is expressed using the {\it WISE}\ magnitude (not corrected for absorption), as function of time in Modified Julian Date. Each epoch (top and lower panel) correspond to approximately 2 days. Infrared variability is clearly notable for both sources. Triangles are used to mark {\it WISE}\ observations that corresponds to upper limits at a given wavelength.
}
\label{fig:lightcurves}
\end{figure*}
In Fig.~\ref{fig:colorcurve} we also show their mid-IR color variability as function of time. In the same figure we also present the single epoch observations when both sources are detected in the first three {\it WISE}\ filters overlaid to the ``{\it WISE}\ Gamma-ray Strip''. { As for the other two BZBs mentioned above these variations are clearly consistent with the region chosen to select RWBLs for WISE J141046.00+740511.2 while for WISE J064459.38+603131 there is partial agreement}. It is worth mentioning that the location of integrated mid-IR colors for both WISE J064459.38+603131 and WISE J141046.00+740511.2 could appear different from that where they exhibit variability in the mid-IR color diagram. The underlying reason is that there are many {\it WISE}\ single-epoch observations when sources (see Fig.~\ref{fig:lightcurves}) are not detected that contribute to the color change for the integrated values. It is also worth noting that for the two RWBL candidates the uncertainties on the {\it WISE}\ magnitudes are greater than for the two bright BZBs selected above as example.
\begin{figure*}
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./colorcurve.pdf}
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./movingonstrip2.pdf}
\caption{Left panel) Mid-IR colors reported as function of time in Modified Julian Date for both RWBL selected in our analysis: WISE J064459.38+603131 (top) and WISE J141046.00+740511.2 (bottom). For WISE J064459.38+603131 we only considered {\it WISE}\ observations when the source is detected in all the three filters while for WISE J141046.00+740511.2 we computed the [4.6]-[12]$\mu$m color-lightcurve also using the upper limits (red triangles down). Right panel) Same of Figure~\ref{fig:moving1} but for the two RWBL selected: WISE J064459.38+603131 (purple squares) and WISE J141046.00+740511.2 (magenta diamond). The variations found in the single-epoch {\it WISE}\ observations when both sources are detected { in the first three mid-IR filters is clearly consistent with the boundaries of the portion of the ``{\it WISE} Gamma-ray strip'' for WISE J141046.00+740511.2 and partially in agreement for WISE J064459.38+603131.}
}
\label{fig:colorcurve}
\end{figure*}

On the basis of this investigation of mid-IR variability and spectral behavior we conclude that the evidence supports the interpretation of these sources as RWBLs.

\section{{\it SWIFT}\ observations}
\label{sec:swift}
Several analyses were carried out to study the nature of WELQs in the X-rays \citep[see e.g.,][]{shemmer09,wu12}. These observations were also used to build their broad band spectral indices and/or their SEDs \citep[see e.g.,][]{plotkin10} to verify their potential similarities with BL Lac objects.
Since the two sources identified here as RWBLs were discovered during a campaign aimed to search for UGS counterparts and all the {\it Fermi}\ UGSs are in a {\it SWIFT}\ observational program \citep{stroh13}\footnote{http://www.swift.psu.edu/unassociated/}, we reduced the {\it SWIFT}\ data available to date to build the RWBL SEDs.

WISE J064459.38+603131 was observed by {\it SWIFT}\ both with the X-Ray Telescope (XRT), for $\sim$ 3 ksec, and with the Ultraviolet/Optical Telescope (UVOT) instruments on board but for the latter the only ultraviolet filter used was M2. On the other hand, for WISE J141046.00+740511.2, there are optical data in the U band and all the UV filters available, together with about 11 ksec of X-ray data.

Here adopted the standard XRT and UVOT data reduction procedure, extensively used in our previous analyses \citep[see e.g.][and references therein]{specevbl,BATcan,3cunid} while only the very basic details are described below. XRT data were processed with the XRTDAS software package (v.3.0.0) developed at the ASI Science Q9 Data Center (ASDC) and distributed within the HEASOFT package (v.6.18) by the NASA High Energy Astrophysics Archive Research Center (HEASARC). All the XRT observations were carried out in the photon counting (PC) readout mode. Event files were calibrated and cleaned applying standard filtering criteria with the {\it XRTPIPELINE} task and using the latest calibration files available in the Swift CALDB. Events in the energy range 0.3--10 keV with grades 0-12 were selected in the analysis according to the standard guidelines. Exposure maps were also created using {\it XRTPIPELINE}. The {\it SWIFT}\ observations available for both sources have variable exposures,  we merged cleaned event files obtained with the above procedure using {\it XSELECT} task and considering only {\it SWIFT}\ observations with telescope aim point falling in a circular region of 12 arcmin radius centered in the median of the individual aim points, to get a uniform exposure. The corresponding merged exposure maps were then generated by summing the exposure maps of the individual observations with {\it XIMAGE} package. Source detection of in the XRT images was performed using the detection algorithm {\it DETECT}.

We collected 13 counts for WISE J064459.38+603131 in the merged XRT event file while there are 36 counts for WISE J141046.00+740511.2, making both sources well detected above the 5$\sigma$ threshold. However the poor number of counts available, corresponding to a count rate of 0.004 cts/s and 0.003 cts/s, respectively, did not permit us to carry out a detailed spectral analysis. Thus assuming an absorbed power-law model with a photon index of 2.2 and with the Galactic hydrogen column density \citep{kalberla05} and estimated the unabsorbed X-ray flux of in the 0.5-10 keV for WISE J064459.38+603131 and of for WISE J141046.00+740511.2 using WebPIMMS\footnote{https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3pimms/w3pimms.pl}.

Details on the UVOT reduction procedure used here are described in Massaro et al. (2008b) and Maselli et al. (2008). We performed the photometric analysis using a standard task {\it UVOTSOURCE} available in the HEASOFT package (v.6.18). Counts were extracted from a 5 arcseconds radius circular region in the V, B, and U filters as well as for the UV filters (UVW1,UVM2, and UVW2). The count rate was corrected for coincidence loss and the background subtraction was performed by estimating its level in an offset annular region centered on the source position with inner radius of 7.5 arcseconds and outer radius of 15 arcseconds. Since estimate of magnitude errors is complex due to possible instrumental systematics and calibration, here we only showed a typical uncertainty of 5\% for all filters.

The values of the observed AB magnitude for WISE J064459.38+603131 in the UV band at 1.36$\times$10$^{15}$Hz is 19.69 mag (M2 UVOT filter) while for WISE J141046.00+740511.2 the AB magnitudes measured in the UVOT filters are: 19.23 mag (U), 19.85 mag (W1), 20.24 mag (M2) and 20.31 mag (W2), respectively.

In Fig.~\ref{fig:swift} we show the XRT event file in the 0.3-10 keV energy range and one of the UVOT images for both the RWBL selected. We overlaid to the {\it SWIFT}\ images the positional uncertainty region, at 95\% level of confidence, from the 2FGL and the 3FGL catalogs for the nearby UGSs. WISE J064459.38+603131 lies within the error ellipse of the 2FGL, as occurred when it was selected \citep{massaro12a} but not within the one of the 3FGL, while WISE J141046.00+740511.2 that is associable to the 3FGL source, it was not found a potential counterpart of the 2FGL UGS.
\begin{figure}
\includegraphics[height=6.6cm,width=8.8cm,angle=0]{./wisej064459swift.pdf}
\includegraphics[height=6.6cm,width=8.8cm,angle=0]{./wisej141046swift.pdf}
\caption{Top panels) The XRT (left) and the UVOT (M2 filter) images of the field where WISE J064459.38+603131 lies. The yellow dashed ellipse marks the positional uncertainty region of the closest 2FGL unassociated source while the white one the error ellipse from the 3FGL catalog, both at 95\% level of confidence; { green circles mark the location of the X-ray counterpart.} Bottom panel) Same of top panels for WISE J141046.00+740511.2. The UVOT image in this case is in the optical energy range (i.e., U filter) instead of the UV.
}
\label{fig:swift}
\end{figure}

{\it SWIFT}\ observations, both in the UV and in the X-ray energy ranges are used in the following to build the SEDs for both the RWBLs selected.

\section{Spectral Energy Distributions}
\label{sec:seds}
For the two RWBL selected we also built the SEDs to check if they are similar or not to those of WELQs \citep{lane11}. Analysis of the SEDs is equivalent to comparing broad band spectral indices that could give an additional indication of the differences between BL Lacs and WELQs \citep[see e.g][]{wu12}.

In the case of WISE J064459.38+603131, the source has an upper limit in the NVSS catalog of 0.13 mJy/beam while WISE J141046.00+740511.2 1.69 mJy/beam. The former is detected in all the {\it WISE}\ bands and in the 2MASS. From the {\it SWIFT}\ observations we also got an estimate of the UV flux density in the M2 filter and the X-ray flux. The latter source still has a {\it WISE}\ counterpart detected at all mid-IR frequencies but it does not have a 2MASS correspondence. {\it SWIFT}\ observations provided us flux density estimates in the optical (U) and in all the UV filters (W1, M2, W2) together with the X-ray one.

We corrected the all the infrared, optical and ultraviolet flux densities for the Galactic absorption following the Draine (2003) and the reddening corrections of Schlafly \& Finkbeiner (2011), with the only exception for the {\it WISE}\ magnitudes at 12$\mu$m and 22$\mu$m.

In Fig.~\ref{fig:sed} we show the SEDs for both WISE J064459.38+603131 and WISE J141046.00+740511.2. For the former SED the infrared emission appear to be compatible with a log-parabolic shape \citep{massaro04,massaro06} with the only exception of the UVOT flux density in the M2 filter. This could be due to the presence of the host galaxy and/or to variability since UV and IR observations were not taken simultaneously and both of them are integrated over different epochs. The latter source: WISE J141046.00+740511.2 has also an SED well described by a log-parabolic shape but in this case the {\it SWIFT}\ optical and UV data seems to mark the presence of an host galaxy.
In both panels of Fig.~\ref{fig:sed} we also report the unabsorbed X-ray flux computed from the source counts of the XRT event file. Since this X-ray flux is an estimate we neglect second order corrections to convert the X-ray flux into a $\nu F_{\nu}$ SED value.
\begin{figure*}
\includegraphics[height=6.6cm,width=8.8cm,angle=0]{./wisej064459sed.pdf}
\includegraphics[height=6.6cm,width=8.8cm,angle=0]{./wisej141046sed.pdf}
\caption{The SEDs of both WISE J064459.38+603131 (left panel) and WISE J141046.00+740511.2 (right panel). Data points are shown as black circles and labels indicate the telescope used to obtain them. The black line marks the broad-band log-parabolic function generally adopted to described non-thermal SEDs of BL Lacs. For WISE J141046.00+740511.2 the dashed red line points towards the host galaxy interpretation for the optical-UV data. Unabsorbed X-ray fluxes, in the 0.5-10 keV energy range, estimated from the XRT source counts are also shown as red triangles in both panels.
}
\label{fig:sed}`
\end{figure*}

The agreement between the multifrequency data and the broad band log-parabolic SED support the interpretation of WISE J064459.38+603131 and WISE J141046.00+740511.2 as RWBLs.

\section{Summary and Conclusions}
\label{sec:summary}
\begin{itemize}
\item The mid-IR colors of BL Lacs appear to be clearly different from those of WDs and WELQs. This strongly indicate that the {\it WISE}\ color-color diagram built with the mid-IR magnitudes at 3.4$\mu$m, 4.6$\mu$m and 12$\mu$m together with optical spectroscopic observations are a powerful diagnostic tool to distinguish these three source classes. We only found few contaminants within the WD and the WELQ catalogs used in our analysis. However, for all these cases, the combination of mid-IR and optical spectroscopic data allowed us to exclude them as BL Lacs.
\item Applying this mid-IR color-based analysis we conclude that WISE J064459.38+603131 and WISE J141046.00+740511.2, lying within the positional uncertainty regions of two distinct UGSs, are the first two RWBLs detected so far. Both of them are also associable with 1FHL sources thus strengthening our interpretation, since the high energy sky above 10 GeV is mostly dominated by BZBs. Both sources appear to be variable in the mid-IR as expect for BZBs and they also show a broad-band SEDs, built thanks to the {\it WISE}\ and the {\it SWIFT}\ observations, consistent with a log-parabolic shape.
\item In addition to them 2QZ J215454.3-305654 could also potentially be a RWBL, but lacking of the mid-IR detection at 12$\mu$m it is not possible to confirm its nature using the mid-IR color-color plot. In addition the two contaminants: WISE J001514.87-103043.4 and WISE J150427.69+543902.2, found in the WELQ sample used here, deserve a deeper investigation to clarify their nature.
\end{itemize}

We highlight that the strategy used to search for potential counterparts of UGSs via X-ray follow up observations has been crucial in both cases to discover WISE J064459.38+603131 and WISE J141046.00+740511.2, thus emphasizing that combining X-ray and mid-IR observations could be crucial for future searches of low energy counterparts of gamma-ray sources. However, we also stress that such discovery has been done thanks to optical spectroscopic follow up campaigns that finally revealed the source nature. Since both these two sources are not only associable with 3FGL objects but also with {\it Fermi}\ sources detected above 10 GeV their existence challenges current gamma-ray association methods mostly based on radio surveys that do not allow to catch these potential counterparts.  The existence of RWBLs will also strong impacts on jet physics since they could allow us to study sources whose emission is dominated by relativistic jets that do not tend to form extended radio structures as their parent population of radio galaxies.

Finally, we remark that deeper radio follow up observations than those available to date could potentially reveal their low energy emission at GHz frequencies. However the ratio between the infrared and the radio emission of so called RWBLs, being below the flux limit of current radio surveys is unexpectedly unusual. In Fig.~\ref{fig:ratios} we show the distribution of such ratio to highlight how the two selected RWBLs lie in the tail of the BL Lac population being below the 0.5 value.
\begin{figure}[!ht]
\includegraphics[height=8.cm,width=8.8cm,angle=0]{./ratios.pdf}
\caption{The distributions of log of the ratio between the radio flux densities at 1.4 GHz and the mid-IR one at 3.4$\mu$m shown for BZBs (blue) and BZQs (red) in our sample. The two selected RWBLs lie at least on the low tail of the blue distribution drawn for the BZBs thus being extreme sources even if follow up observations will reveal their radio emission.
}
\label{fig:ratios}
\end{figure}

\acknowledgements
We thank the anonymous referee for useful comments that led to improvements in the paper.
F.M. gratefully acknowledges the financial support of the Programma
Giovani Ricercatori -- Rita Levi Montalcini -- Rientro dei Cervelli (2012) awarded by the Italian Ministry of Education, Universities and Research (MIUR).
H.A.S. acknowledges partial support from NASA Grants NNX15AE56G and NNX14AJ61G.
This research has made use of data obtained from the high-energy Astrophysics Science Archive
Research Center (HEASARC) provided by NASA's Goddard Space Flight Center.
The NASA/IPAC Extragalactic Database
(NED) operated by the Jet Propulsion Laboratory, California
Institute of Technology, under contract with the National Aeronautics and Space Administration.
Part of this work is based on the NVSS (NRAO VLA Sky Survey):
The National Radio Astronomy Observatory is operated by Associated Universities,
Inc., under contract with the National Science Foundation and on the VLA Low-frequency Sky Survey (VLSS).
This publication makes use of data products from the Wide-field Infrared Survey Explorer,
which is a joint project of the University of California, Los Angeles, and
the Jet Propulsion Laboratory/California Institute of Technology,
funded by the National Aeronautics and Space Administration.
TOPCAT\footnote{\underline{http://www.star.bris.ac.uk/$\sim$mbt/topcat/}}
\citep{taylor05} for the preparation and manipulation of the tabular data and the images.