\section{Introduction}
Organic-based materials are attractive for optoelectronic device applications, notably due to their low fabrication cost and their relative ease to produce and characterize.\cite{Prineted_book} Not only can the structural properties of these materials be tuned through the functionalization of molecules,\cite{Cheng_ChemRev2009} but they are also composed of elements which are Earth-abundant. In contrast to conventional inorganic electronic materials, organic compounds bring in flexibility,\cite{LOGOTHETIDIS200896} biocompatibility\cite{Rivnay-NatComm:2016} and biodegradability,\cite{Lei5107} as well as self-healing properties.\cite{Oh2016, Ocheje2017} Organic conducting and semiconducting materials hold promises for several application niches, including next-generation wearable and printed photovoltaics,\cite{Facchetti_ChemMat2011,Lipomi_AdvMat2011} fuel cells, \cite{guo2012self,winther2008high} thermoelectrics, \cite{Dubey_JPolSci2011,Wang_AdvEnMat2015,Wei_Materials2015,Bubnova_NatMat2011} and other optoelectronics applications.\cite{sun2015review, Roch2017}\\

One of the fundamental challenges for the design of organic optoelectronics lies in the intrinsic structural disorder of the materials. This disorder occurs on multiple length scales starting from the conformations of single molecules and the nearest-neighbor packing to the formation of multi-molecule domains and nanocrystals. The electronic properties of organic materials are highly sensitive to the packing of composing molecules hence dependent on the processing conditions.\cite{Hinton-JACS:2018} Fast optical probing of local electronic couplings can benefit both applied and fundamental research. On one hand, such a method brings the possibility to combine continuous testing of devices with roll-to-roll device manufacturing technology.\cite{Tabor:2018, MI-report} On the other hand, optical characterization techniques can advance our understanding of charge transport in organic structures. In particular, UV-Vis, XPS, and Raman scattering measurements of thin films of conductive polymers can provide insight on composition and electronic structure, including the nature of charge carriers.\cite{Luo-JMCA:2013, Ouyang-polymer:2004} \\

The microscopic structure of molecule and polymer packing is difficult to measure directly. In contrast, obtaining optical spectra such as IR absorption, Raman scattering, electronic absorption, and fluorescence is more straightforward and requires sufficiently less experimental effort. Therefore, indirect characterization methods play a key role in evaluating the level of disorder. Both optical and electronic transport properties are influenced by the microscopic molecular packing. For example, in the simplest qualitative picture, the close proximity of two molecules yields an overlap of electronic clouds which leads to charge transfer. This proximity also leads to a F{\"o}rster coupling between electronic excitations, which can be observed as changes in the lines in the electronic absorption spectra.\cite{Saikin-Nanophotonics:2013} Moreover, weak charge-transfer excitations are observed if the electronic coupling between molecules becomes sufficient. Because both effects are caused by the molecular interactions, in principle it is possible to find a machine learning model that correlates both of them. The conventional computational approach involves three steps: (i) building physical models that describe both properties of interests; (ii) fitting the parameters of the models to experimental data, e.g., absorption spectra; and (iii) using the fitted models to describe the other property, e.g., conductivity. Such an approach might be tedious, since the interrelations between these properties can be too complex to derive a simple or tractable physical model.\\

Herein, we report an alternative approach, where the aforementioned physical model is replaced by machine learning (ML) models. To this end, we design a multi-scale computational workflow where the first three steps -- force-field calculations, molecular dynamic simulations, and quantum-based approaches -- generate the data for the ML algorithms. It is of course possible to employ experimental data in addition if it is available. In this work, and to begin with, we employed ML algorithms to identify correlations between the two properties of interests i.e., strength of intermolecular coupling and electronic absorption spectra. Then, we used the trained ML models as a relative classifier of the coupling strength of a given spectrum with respect to a reference coupling, which is to be defined by a scientist for the application at hand. As a model system demonstrating the reliability of the classifier to identify structures with strong electronic couplings from their electronic absorption spectra, we study pairs of poly(3,4-ethylenedioxythiophene) (PEDOT) oligomers. PEDOT is one of the most technologically-developed conducting polymers. Owing to its high hole conductivity and optical transparency in a doped state\cite{Groenendaal_AdvMat2000} it is has been used as transparent contacts in photovoltaic devices, touch screens, and light-emitting diodes. \cite{Lovenich_PolSci2014} PEDOT-based materials are frequently used as a mixture with polystyrene sulfonate polymers (PEDOT:PSS). In this mixture, PEDOT oligomers transfer the charges while PSS chains play the role of a solid electrolyte. This material becomes conductive at high concentrations of dopant.\cite{PEDOT_book}\\

Although multiple experimental studies have addressed the molecular organization of PEDOT-based materials,\cite{Martin-PR:2010,Groenendaal_AdvMat2000, Takano_MacroM2012, Nardes-AM:2007, Crispin-JPS:2003, Kemerink-JPCB:2004, VanDeRuit-AFM:2013, Crispin-ChemMater:2006, Luo-JMCA:2013, Ouyang-polymer:2004} the microscopic electronic states that lead to high conductance and the interplay between these states, optical properties, and the material structure have yet to be determined. As a matter of fact, the key factor for practical applications of PEDOT-based materials lies in understanding the relations between their solid-state packing and their unique electronic properties. Consequently, the main obstacle to elucidating this relationship is the strong structural disorder that appears on multiple length scales and is highly sensitive to the thin film preparation procedure.\cite{Groenendaal_AdvMat2000}\\

Hereafter, we demonstrate that our ML models
confirm the existence of correlation between the coupling strengths and the electronic absorption spectra. We also show the robustness of our ML models with respect to potential spurious statistical correlations to capture the relevant physical correlations. Finally, we use the ML models as classifiers to determine whether a given electronic absorption spectrum of interest relates to a coupling strength above or below an \textit{a priori} selected reference coupling strength. Such an approach has proven to be reliable and robust, reaching an average error rate of only 12\% when employing a Bayesian convolutional neural network. The importance of such a classifier becomes apparent in the context of the self-driving laboratories,\cite{SciRobot18:Roch,Tabor:2018,Hase-2019,MI-report,Dimitrov-ACS:2019} where the goal of the experimentation process to identify fabrication procedures for aggregates yielding high conductivities is embodied as an optimization procedure. \\

\section{Methods}
\label{sec:methods}
This section details the computational workflow designed to generate the data and to correlate electrical and optical properties using ML algorithms. The workflow is depicted in Fig.~\ref{fig:concept}. Each of the five composing steps (i.e., initial structures, refinement, MD simulation, physical models, ML models in Fig.~\ref{fig:concept}A) are described in their corresponding subsections.\\

Hereafter, we assume that the packing of conjugated oligomers in the solid PEDOT:PSS mixture depends on the initial preparation procedure and post-processing steps. For PEDOT:PSS, such post-processing steps, typically evaluated via trial and error, are critical to achieve peak performance. Furthermore, it is hypothesized that PEDOT:PSS solid films consist of grains with a hydrophobic and highly conductive PEDOT-rich core and a hydrophilic insulating PSS-rich shell.\cite{David-2016,Takano_MacroM2012} This phase segregation of PEDOT and PSS occurs on a scale beyond current computational capabilities and, thus, is not captured by our model. Nonethless, our computational workflow allows to study the disorder within each of these grains.\\

\begin{figure}[htb!]
\centering
\includegraphics[width=0.9\textwidth]{computational_pipeline_v7.png}
\caption{The computational pipeline used, from structure generation to correlating $J$ and electronic absorption spectra. (A) General workflow highlighting the steps and summarizing the methods involved. (B) Structures of PEDOT and PSS, represented in the top and bottom panels, respectively. (C) One of the ten supercells of PEDOT:PSS blends. (D) Two distinct example pairs of PEDOT oligomers extracted from the PEDOT:PSS bulk. $J_1$ and $J_2$ are the coupling strengths for each of the pairs. (E) Associated simulated electronic absorption spectra. (F) State, and (G) orbital representation of the monomers and dimers involved in the calculation of the coupling strength, $J$.}
\label{fig:concept}
\end{figure}

\paragraph*{Initial structures.} To screen the orientation stability of the PEDOT:PSS complex within the grains, 100 starting structures were generated using the Packmol software package.\cite{Packmol2009} Each of the complexes consisted of one PEDOT chain with eight 3,4-ethylenedioxythiophene units ($n=8$, Fig.~\ref{fig:concept}B) carrying two positive charges, and two PSS chains consisting of three poylstyrene sulfonate units ($m=3$, Fig.~\ref{fig:concept}B) with one negative SO$_3^-$ and two SO$_3$H groups per chain. In the generation of the initial structures, we constrained the SO$_3^-$ group of the PSS chain to point towards the positive PEDOT chain. These initial structures were then optimized using a classical force field (FF) approach.

\paragraph*{Refinement.} The 20 energetically most stable PEDOT:PSS complexes obtained from the initial structure search were relaxed at the B97-D/6-31G(d,p) level of theory in the gas phase, using the Gaussian software package.\cite{g09} Note that the influence of the solvent was found to be negligible, and that the performance of the B97-D functional on geometries has already been assessed in previous work.\cite{Roch:jctc-2017} Single point energy calculations at the HF/6-31G(d,p)//B97-D/6-31G(d) level were performed on the 20 relaxed complexes to parametrize the charges of the molecular dynamics simulation, as customary with GAFF.\cite{GAFF2004}

\paragraph*{Molecular dynamics simulation.} Ten molecular dynamics simulations were carried out on PEDOT:PSS model systems in periodic cubic boxes of size\footnote{Note: average size of the different boxes for the duration of the NPT simulations} \unit[61.32]{\AA} $\times$ \unit[61.32]{\AA} $\times$ \unit[61.32]{\AA} for a density of \textit{ca.} \unit[1.4]{g/cm$^3$}, with the LAMMPS software package.\cite{PLIMPTON19951} The GAFF force-field was chosen to describe the systems. PEDOT oligomers with eight repeat units ($n=8$, Fig.~\ref{fig:concept}B) and a +2 charge were considered, while the PSS atactic chains consisted of 20 repeat units ($m=20$, Fig.~\ref{fig:concept}B), with four deprotonated units randomly distributed in the sequence of each chain. Note that the PSS chain length was increased from $m=3$ to $m=20$ to better represent experimental blends. Additional details can be found in the Supporting Information (see section \ref{sec:si_md}).

\paragraph*{Physical models to compute $J$, and simulate the electronic absorption spectra.}
For two interacting PEDOT monomers, denoted $A$ and $B$ from hereon, the strength of the charge transfer integral $J$ can be determined in several ways.\cite{Bredas_JACS1983, Bredas_CPL2002, Bredas_JACS2008, martin2000monodisperse, bredas:jcp-2007, Apperloo_ChemEurJ2002, Beljonne_AdvFM2001, Li-JCP:2007, deibel:prb-2011, Kubas-JCP:2014, Cave-JCP:97} Bi-polaron charge transport models have previously been discussed in the case of PEDOT systems.\cite{Bredas_AChR1985} For the sake of simplicity, we assume a single-polaron transport model. Nonetheless, the ML models used for electronic coupling prediction are agnostic to the type of transport and will learn correlations between $J$ and electronic absorption spectra independently of the type of the charge transport model. To calculate the coupling $J$, we used the framework of a tight-binding formalism\cite{MolCrys-Book,Bassani-Book} as well as a Kohn-Sham orbital based method.\cite{bredas:jcp-2007,deibel:prb-2011} In both cases, orbital energies were obtained at the B3LYP/def2-SV(P) level of theory. The choice of basis set and functional balances computational cost and accuracy. Detailed results on the performance of the def2-SV(P) results can be found in the Supporting Information (see section \ref{sec:basis_set}). \\

The nearest-neighbor PEDOT dimers were extracted across the ten simulation boxes. These pairs were selected according to a distance criterion; any pair of PEDOT molecules having at least two heavy (i.e., non hydrogen) atoms at a distance closer than \unit[4]{\AA} is selected. The cutoff distance was taken to be comparable to the sum of the van der Waals radii of these heavy atoms. This procedure lead to the selection of 1,420 PEDOT pairs, ulteriorly used to model the strength of the charge transfer, and to simulate the electronic absorption spectra.\\

Tight-binding formalism is based on the change of orbital energies when going from isolated monomers to dimer systems \cite{Bredas-ChemRev:2004}

\begin{equation}
\centering
J=\sqrt[]{\left(\Delta\epsilon_{\left|AB\right>,+5}^{H,L}\right)^2 - \frac{1}{4}\left[ \left(\epsilon_{\left|A\right>,+3}^{H,H-1} + \epsilon_{\left|B\right>,+2}^{H,H-1}\right) - \left(\epsilon_{\left|A\right>,+2}^{H,H-1} + \epsilon_{\left|B\right>,+3}^{H,H-1}\right) \right]^2},
\label{eq:j}
\end{equation}

\noindent where $\Delta\epsilon_{\left|AB\right>,+5}^{H,L}$ is the splitting between the HOMO and the HOMO-1 level of the dimer $AB$ with the charge $q=+5$, and $\epsilon_{\left|i\right>,+q}^H$ is the HOMO energy of monomer $i$, carrying charge $+q$. Note that the correction due to the offset between the HOMOs of the monomers is negligible; hence, $J$ is mostly governed by the splitting $\Delta\epsilon_{\left|AB\right>,+5}^{H,H-1}$. This formalism considers frontier orbitals assuming that only highest occupied orbitals are hybridized due to the electronic coupling between the oligomers (see Fig.~\ref{fig:concept}F-G). Both the presence of non-equilibrated charges and non-zero spin makes the model rather complicated. Nonetheless, our interpretation of  eq.~\ref{eq:j} is a lower estimate for the electronic coupling between the oligomers. Note that this model can also be used to describe bi-polaron transport.\\

The second approach to calculate charge transfer integrals uses the Kohn-Sham orbitals of isolated monomers as well as the Fock matrix and the overlap matrix of the dimer systems:\cite{bredas:jcp-2007,deibel:prb-2011}

\begin{equation}
J_{\text{AB}}=\frac{F_{\text{AB}}-\frac{1}{2}(F_{\text{AA}}+F_{\text{BB}})S_{\text{AB}}}{1-S_{\text{AB}}^2}
\label{eq:j2}
\end{equation}

\noindent The matrix elements $F_{\text{AB}}=\left<\text{A}|F_{\text{dimer}}|\text{B}\right>$ and $S_{\text{AB}}=\left<\text{A}|S_{\text{dimer}}|\text{B}\right>$ are calculated using the Fock and overlap matrices of a molecular dimer system with a charge of $+4$. The states $\left.|\text{A}\right>$ and $\left.|\text{B}\right>$ are the highest occupied molecular orbitals of the doubly positively charged monomers. The distributions of electronic couplings obtained with eq. \ref{eq:j}, and with eq. \ref{eq:j2} are depicted in Fig.~\ref{sec:si_coupling_strength_distributions}.\\

The electronic absorption spectra of the 1,420 PEDOT pairs were simulated using the computed TD-CAM-B3LYP/def2-SV(P) transitions, in the gas phase. Note that for the simulation of the spectra, the influence of solvation was found to be negligible. We employed a Lorentzian distribution to broaden the TD-DFT transitions. This translates the discrete oscillator strength, $f$, and transition energies, $\omega$, to a continuous spectra to resemble experimental outcomes. The broadening was chosen to be \unit[50]{meV}.\cite{Randall:NL-2018}\\

Some insights about correlations between electronic absorption spectra and intermolecular charge transfer can be obtained only for the case when the coupling is weak. A completely relaxed doubly charged oligomer composed of eight to ten units would have a strong electronic transition at about 0.9 -- \unit[1.0]{eV}.\cite{David-2016} This transition is predominately composed of HOMO and LUMO orbitals. For an oligomer with an odd number of charges, additional HOMO-1 $\rightarrow$ HOMO transition appears at lower frequencies, at \textit{ca.} \unit[0.5]{eV}.  Therefore, the low frequency part of the spectra of dimers without the interaction should be composed of three lines -- a low-frequency, weak transition and a strong doublet. Electronic or excitonic interaction between the molecules modulate the spectra. Specifically, weak charge transfer transitions should appear at the low frequency tail of the spectra. The intensities of these transitions should be sensitive to the coupling strength, while their frequencies should be more stable as determined by the alignment of the molecular energy levels. However, this intuitive picture fails for intermediate and strong intermolecular couplings. In the latter case, the intramolecular states hybridize, which in turn leads to a realignment of the energy levels and a redistribution of the oscillator strengths among multiple transitions. Yet another advantage of our ML approach is that it allows to capture correlations between the electronic interaction and optical spectra independently of the coupling regime used as our ML independent variable.

\paragraph*{Machine learning models to identify correlations between electronic absorption spectra and coupling strengths.} Correlations between electronic absorption spectra and coupling strengths were identified with ML models at different levels of complexity. To mimic experimental conditions, we encoded the Lorentzian broadened electronic absorption spectra based on their intensities at specific frequencies, using a \unit[1]{meV} binning on the considered frequency domain (\unit[0.5]{eV} to \unit[1.5]{eV}). A total of 170 (\unit[11.97]{\%}) of the 1,420 spectra-coupling pairs were randomly selected to construct a test set. The remaining 1,250 (\unit[88.03]{\%}) of the dataset were used for ten-fold cross-validation. The size of the dataset motivates the use of Bayesian models for a robust and transferable identification of relevant physical correlations. \\

    Specifically, we employed three different classes of Bayesian models (see Fig.~\ref{fig:ml_models}): (i) Bayesian linear regression models assume a linear dependence of coupling strengths on electronic spectra and thus present the simplest approximation; (ii) Bayesian multi-layer perceptrons (MLPs) are Bayesian generalizations of conventional deterministic MLPs with similar flexibility to model non-linear relations while retaining the robustness to overfitting of Bayesian methods; and (iii) Bayesian one-dimensional convolutional neural networks (CNNs) are special cases of MLPs which have the potential to efficiently exploit spatial correlations in the presented features due to their sparse nature. All models were set up to predict coupling strengths directly from the intensities of the associated electronic absorption spectra at different frequencies. Additionally, we constructed Bayesian MLPs which are trained on a compressed representation of the electronic absorption spectra obtained from principle component analysis (PCA). All models were trained based on an early-stopping criterion. Hyperparameters for all four models are optimized in a random grid search, and the details are provided in the Supporting Information (see section \ref{sec:si_hyperparameter_optimization}).

    \begin{figure}
        \centering
        \includegraphics[width=0.75\textwidth]{ML_models_shadow-effect_4.png}
        \caption{Schematic representation of (A) the data processing and (B) the three classes of ML models, depicting Bayesian linear regression models, Bayesian multi-layer perceptrons (MLPs), and Bayesian one-dimensional convolutional neural networks (CNN). Note that the MLPs are trained either on the raw intensities (blue array) or on the compressed representation after a principle component analysis (green array).}
        \label{fig:ml_models}
    \end{figure}

\paragraph*{Machine learning models as relative classifiers.} The aforementioned ML models trained for predicting absolute values of coupling strengths from electronic absorption spectra can be used to classify the conductivity performance associated with the electronic absorption spectra of the materials. Instead of asking for the absolute value of the coupling strength, the model provides an estimate for whether the considered coupling is above or below a reference. This reference is a hyperparameter defined by a scientist as a threshold for high and low values of $J$.\\

\section{Results and discussion}
We begin by discussing our results with the performance of the ML models to identify correlations and to predict absolute values of coupling strengths from electronic absorption spectra. Fig.~\ref{fig:test_set_predictions} illustrates the accuracies of all four models to predict coupling strengths computed from the tight-binding formalism with all four models after a full hyperparameter optimization. Then, we detail the test designed to assess the performance of the ML models to capture relevant physical correlation. We continue our discussion with the results obtained when the ML models are used as relative classifier, where associated error rates are reported in Fig.~\ref{fig:contingency_plot}. We also highlight the practicality of such an approach in discovery applications with the self-driving laboratories. Finally, we discuss the robustness of our ML models upon variations of the peak broadening, which experimentally translates to noise. Fig.~\ref{fig:predictions_for_broadenings} illustrates the performance of our ML models at different broadening.

    \begin{figure}[!ht]
        \centering
        \includegraphics[width = 0.75\textwidth]{test_set_predictions_adjusted_labels.png}
        \caption{Coupling strengths predicted by the employed machine learning models in comparison to coupling strengths obtained with the tight-binding formalism. All reported predictions are shown for the test set, and were obtained from averaging $200$ prediction samples from each of the models. Prediction uncertainties are color-coded. We report three comparative metrics to assess the prediction accuracies of each model: root-mean-square deviation (RMSD), mean-absolute deviation (MAD), and coefficients of determination ($r^2$). Based on all three metrics, Bayesian CNNs provide the most accurate predictions.}
        \label{fig:test_set_predictions}
    \end{figure}

    \subsection{Correlation between electronic spectra and coupling strengths}
    All constructed ML models are able to predict coupling strengths at positive coefficients of determination ($r^2$), which indicate that the coupling strengths correlate with the electronic absorption spectra of the PEDOT dimers and that the models are capable of identifying this correlation. Bayesian CNNs provide the most accurate predictions based on all computed comparative metrics ($r^2 = 0.313$, RMSD $= 0.107$ eV, MAD $= 0.082$ eV)\footnote{RMSD: root-mean-square deviation, MAD: mean-absolute deviation, and $r^2$: coefficients of determination} and Bayesian linear regression yields the least accurate predictions ($r^2 = 0.163$, RMSD $= 0.118$ eV, MAD $= 0.092$ eV). We further observe an improved prediction accuracy when compressing the electronic absorption spectra via PCA for the Bayesian MLP models. \\

    Despite the relatively small size of the dataset, we found that the studied ML models, most notably the Bayesian CNN, present an efficient approach to identify the relevant correlations. Estimations of the sampling efficiency of the Bayesian CNN model suggest that it can be trained to reach similar prediction accuracies with only 850 instead of 1,250 training points. No significant improvement in the prediction accuracy is observed when increasing the size of the training set from beyond 850 to 1,250 examples (see Supporting Information, section \ref{sec:si_sampling_efficiency}). This observation, in conjunction with the generalization of the models observed for the test set predictions, indicates that the Bayesian CNN exploits all identifiable correlations to their full extent. \\

	To ensure that our ML models did not capture spurious correlations that could arise from the methods and formalisms employed to compute the electronic absorption spectra and model the coupling strengths, we tested for the nature of the identified correlations by training the Bayesian CNNs to predict couplings strengths obtained with the Kohn-Sham orbital formalism from the same electronic absorption spectra. Details are provided in the Supporting Information (see section \ref{sec:si_cross_predictions}). The trained Bayesian CNNs achieve prediction accuracies of $r^2 = 0.264$ on the same test set. In addition, we constructed a hybrid dataset, where half of the couplings are randomly chosen from the tight-binding formalism and the other half from the Kohn-Sham orbital formalism. Again, the trained Bayesian CNNs achieve prediction accuracies of $r^2 = 0.280$ indicating that the presented ML models do not capture potential spurious statistical correlations but extract the relevant physical correlations. \\

    \subsection{Machine learning models as relative classifiers}
    The practicality of the presented trained ML models becomes more apparent when weakening the requirement of accurate absolute predictions of the coupling strengths to accurate relative predictions, which are of interest in discovery applications. Rather than requesting an estimate for the actual numeric value of the coupling strength, the trained model is used to determine if a given electronic absorption spectrum of interest relates to a coupling strength above or below an \textit{a priori} selected reference coupling strength. For such scenarios, the prediction accuracy of the model can be assessed by treating it as a binary classifier to determine true positive and true negative rates for different reference coupling strengths.\\

We assess the prediction accuracies of such relative classifiers by estimating the probability of the model to make a correct prediction, i.e., predicting the coupling to be above the reference when it is above or predicting the coupling to be below the reference when it is below, versus the probability of the model to make an incorrect prediction, i.e., predicting the coupling to be above the reference when it is below and vice versa (see Fig.~\ref{fig:contingency_plot}). These probabilities are estimated for different reference coupling strengths, spanning the entire range of coupling strengths computed with the tight-binding formalism. \\

        \begin{figure}[!ht]
            \centering
            \includegraphics[width = 0.75\textwidth]{summed_error_rates_horizontal_v3.png}
            \caption{Probabilities to make correct predictions for $J$ regarding the order with respect to another reference $J$ (left panel) vs. making incorrect predictions (right panel).}
            \label{fig:contingency_plot}
        \end{figure}

        Fig.~\ref{fig:contingency_plot} illustrates the error rates, i.e., wrongly predicting a coupling strength to be above or below the considered reference coupling strength, for the trained ML models along with two na\"ive models for comparison. The most na\"ive model draws random samples from a uniform distribution to predict the coupling strength for a given electronic absorption spectrum (model-1, depicted in light grey in Fig.~\ref{fig:contingency_plot}). A slightly more sophisticated na\"ive model predicts by drawing random samples from the distribution of coupling strengths (model-2, depicted in dark grey in Fig.~\ref{fig:contingency_plot}).\\

	 We find that the error rates for the two na\"ive models yield the largest error rates: \unit[31.9]{\%} and \unit[20.2]{\%}. Bayesian linear regression scored an average error rate of \unit[19.0]{\%}, which, despite its simplicity, already provides an advantage over simple models and captures some of the relevant correlation in the dataset. The lowest error rate is observed for the Bayesian CNN with a \unit[12.6]{\%} average error for coupling strengths chosen within the range of smallest and largest computed coupling strengths. Additionally, it is noted that Bayesian CNN never exceeds an error of \unit[38]{\%} for any chosen reference coupling. In fact, if the focus of discovery process is to identify fabrication procedures and post-processing steps leading to large coupling strengths (above \unit[0.6]{eV}), the Bayesian CNN yields error rates of less than \unit[10]{\%}. We suggest that the trained ML models can be applied to classify the coupling strengths with respect to a reference coupling strength with reasonable confidence.\\

    \subsection{Robustness of the machine learning models}
    Finally, we estimate the dependence of the model performances on the particular choice of the peak broadening. While we demonstrated that for this particular choice, the trained ML models are indeed capable of identifying the relevant correlations between the electronic absorption spectra and the coupling strengths, experimentally obtained electronic absorption spectra might be noisy and feature peaks at slightly varying broadenings. The robustness of the model predictions for different broadenings is tested by predicting coupling strengths from electronic absorption spectra at different broadenings, ranging from \unit[5]{meV} to \unit[1,000]{meV}. Note that while a \unit[10]{meV} broadening is too small for a room temperature measurement, a \unit[200]{meV} broadening would correspond to an unphysically fast dephasing rate. Fig.~\ref{fig:predictions_for_broadenings} summarizes the coefficients of determination for predictions of coupling strengths from electronic absorption spectra generated at different broadenings. Note that all coupling strengths are predicted by models, which were trained on electronic absorption spectra at a \unit[50]{meV} broadening.

    \begin{figure}[!ht]
        \centering
        \includegraphics[width = 0.75\textwidth]{broadenings_labelled.png}
        \caption{Prediction accuracies measured with $r^2$ coefficients for all models when predicting couplings from spectra at different broadenings. All three models depicted have been trained on a spectra at a \unit[50]{meV} broadening (dashed vertical line). The regime ranging from \unit[10]{meV} to \unit[200]{meV} (highlighted in white) corresponds to the expected experimental resolution: while a \unit[10]{meV} broadening is too small for a room temperature measurement, a \unit[200]{meV} broadening corresponds to an unphysically fast dephasing rate. Note that the horizontal axis, \textit{Broadening [meV]}, is not linear.}
        \label{fig:predictions_for_broadenings}
    \end{figure}

    We find that prediction accuracies of the Bayesian linear regression model are mostly insensitive to the particular broadening value (yellow trace). Only for very small broadenings below \unit[10]{meV} and very large broadenings above \unit[200]{meV} degradations in the predictive power can be observed. Bayesian MLPs (turquoise trace) show faster degradations in their prediction accuracy for small and large broadenings, but maintain comparative predictive powers across broadenings of \unit[25]{meV} to \unit[100]{meV}. Bayesian CNNs (blue trace) are the least robust with respect to changes in the broadening, with accurate predictions only within the \unit[45]{meV} to \unit[80]{meV} interval. Nevertheless, they demonstrate their predictive power for varying broadenings despite having been trained on peaks of one particular broadening, indicating that ML models can indeed be trained to identify relevant correlations without being overly sensitive to the broadening of the peaks.  \\

\subsection{Associating ML models with experimental studies}
A long-standing goal of experimental materials characterization of organic optoelectronic materials is a map of how electronic properties are distributed in space as a result of different instantaneous molecular configurations. Scanning probe\cite{ODea-MRS:2012, Groves-ACR:2010, Eisfeld-2018} and super-resolution optical measurements\cite{Penwell-NatMater:2017} can provide readouts of electronic properties on length scales below the diffraction limit. Simultaneously, single-particle measurements\cite{Barbara-ACR:2005, Thiessen-PNAS:2013} can provide a bottom-up understanding of how optoelectronic properties evolve from molecular precursors. For conductive polymers like PEDOT:PSS, single-particle measurements have been particularly difficult to employ due to the lack of emission in these materials because of rapid quenching. Simultaneously, single-particle measurements have tremendous spectroscopic utility due to reduced inhomogeneous broadening. Use of high quality-factor optical microresonators as the readout for ultrasensitive photothermal spectroscopy\cite{Heylman-NatPhoton:2016} has allowed the first single-particle optical measurements to be performed on PEDOT:PSS,\cite{Randall:NL-2018} even down to a single or a small number of polymer strands. This study provided an experimental bound for the line broadening used in the above simulations. More recently, optical microresonator spectroscopy has been used to show how annealing processing act on single PEDOT:PSS polymer particles.\cite{Rea-InReview} A means of directly connecting spectral measurements on single PEDOT:PSS polymer strands and particles to electronic couplings would significantly amplify the information content of these experiments.\\

\section{Conclusion}
    Our findings suggest that ML models can identify physical correlations between the measurable electronic absorption spectra and the strength of intermolecular electronic couplings, which in turn determine the charge transport. While the presented ML models provide coupling strength estimates with limited accuracy, relative estimates with respect to reference coupling strengths show promising error rates. Using the trained Bayesian CNN model to classify given electronic absorption spectra above or below an \textit{a priori} selected reference coupling strength displays an error rate of only \unit[12.6]{\%}, and as low as \unit[10]{\%} at high coupling regime. With such a promising error rate, we suggest to use the trained models as classifiers to evaluate performance of fabrications procedure and post-processing steps. Further investigations towards the construction of reliable and transferable ML models, notably the usage of ensemble methods such as adaboost,\cite{Freund:1997,Freund:1996} or mixture density networks,\cite{Bishop:1994,Richter:2003} might allow for more detailed insights into the relation between couplings and electronic absorption spectra. Another important venue for improvement of our approach is the incorporation of features, such as structural information, which would introduce the notion of similarity between complexes.\\

	We believe that the combination of the developed approach with spectroscopy techniques and its integration with the self-driving laboratories\cite{SciRobot18:Roch,Tabor:2018,Hase-2019,MI-report,Dimitrov-ACS:2019} has the potential to enhance characterization and accelerate optimization of organic materials. As experimental approaches for providing optical readouts improve in sensitivity,  spatial resolution, and access to different spectral features, growth in theoretical treatments will allow one to draw deeper connections between these measurements and the underlying molecular structure. We also envision the use spectroscopic methods to measure spectra of nanoaggregates with high-finesse toroidal optical cavities.

\begin{acknowledgement}

The computations were done on the Arran cluster supported by the Health Sciences Platform (HSP) at Tianjin University, P.R. China, on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University, USA, and on the computational resource ForHLR II funded by the Ministry of Science, Research and the Arts Baden-W{\"u}rttemberg and DFG (``Deutsche Forschungsgemeinschaft''). L. M. R. and A. A.-G. acknowledge Natural Resources Canada (EIP2-MAT-001) for their financial support. F. H. acknowledges support from the Herchel Smith Graduate Fellowship and the Jacques-Emile Dubois Student Dissertation Fellowship. P. F. has received funding from the European Unions Horizon 2020 research and innovation program under the Marie Skodowska-Curie grant agreement MolDesign No 795206. L. M. R., S. K. S., F. H., P. F., and A. A.-G. thank Dr. Anders Fr{\o}seth for generous support. R. H. G. acknowledges the National Science Foundation (NSF: DMR-1610345). S. L. acknowledges the Real Colegio Complutense in Harvard for a Research Grant, and to the Spanish Ministerio de Ciencia e Innovaci{\'o}n for a Fellowship through the Salvador de Madariaga Program.

\end{acknowledgement}

\begin{suppinfo}

The following files are available free of charge.
\begin{itemize}
  \item Supporting\textunderscore Information.pdf: Supporting information to the main manuscript
\end{itemize}

\end{suppinfo}

\newpage
\setcounter{section}{0}

\begin{center}
\Large{\textbf{SUPPORTING INFORMATION}}\\
\Large{\textbf{From absorption spectra to charge transfer in PEDOT nanoaggregates with machine learning}}
\end{center}

\section{Assessment of the basis set}
\label{sec:basis_set}
\begin{table}[htb!]
\centering
\begin{tabular}{l l | c c | c}
\hline\hline
\multicolumn{2}{l}{Pair \#1} & $\epsilon_\mathrm{HOMO}$ & $\epsilon_\mathrm{HOMO-1}$ & $\Delta \epsilon$\\
\hline
\multirow{4}{*}{$\left| \mathrm{mono}_A\right>_{+2}$} & Def2-SVP & -9.134 & -10.018 & 885 \\
& Def2-SVPD & -9.254 & -10.141 & 887 \\
& Def2-SV(P) & -9.144 & -10.030 & 885 \\
& Def2-TZVP & -9.164 & -10.051 & 887 \\
\hline
\multirow{4}{*}{$\left| \mathrm{mono}_A\right>_{+3}$} & Def2-SVP & -11.850 & -12.420 & 570 \\
& Def2-SVPD & -11.934 & -12.505 & 571 \\
& Def2-SV(P) & -11.863 & -12.434 & 571 \\
& Def2-TZVP & -11.842 & -12.411 & 569 \\
\hline
\multirow{4}{*}{$\left| \mathrm{mono}_B\right>_{+2}$} & Def2-SVP & -9.139 & -10.067 & 928 \\
& Def2-SVPD & -9.263 & -10.189 & 926 \\
& Def2-SV(P) & -9.149 & -10.077 & 929 \\
& Def2-TZVP & -9.171 & -10.096 & 925 \\
\hline
\multirow{4}{*}{$\left| \mathrm{mono}_B\right>_{+3}$} & Def2-SVP & -11.932 & -12.518 & 586 \\
& Def2-SVPD & -12.017 & -12.605 & 588 \\
& Def2-SV(P) & -11.944 & -12.531 & 586 \\
& Def2-TZVP & -11.921 & -12.510 & 589 \\
\hline
\multirow{4}{*}{$\left| \mathrm{dimer}\right>_{+5}$} & Def2-SVP & -11.721 & -12.594 & 873 \\
& Def2-SVPD & -11.860 & -12.733 & 873 \\
& Def2-SV(P) & -12.319 & -12.809 & 489 \\
& Def2-TZVP & -12.358 & -12.739 & 381 \\
\hline\hline
\end{tabular}
\caption{Energies of the HOMO and HOMO-1 molecular orbitals of pair \#1, in [eV], at the CAM-B3LYP level of theory. $\Delta \epsilon$ is the energy difference between the HOMO and the HOMO-1, in [meV].}
\label{table:pair1_BSet}
\end{table}

\begin{table}[htb!]
\centering
\begin{tabular}{l l | c c | c}
\hline\hline
\multicolumn{2}{l}{Pair \#2} & $\epsilon_\mathrm{HOMO}$ & $\epsilon_\mathrm{HOMO-1}$ & $\Delta \epsilon$\\
\hline
\multirow{4}{*}{$\left| \mathrm{mono}_A\right>_{+2}$} & Def2-SVP & -9.126 & -10.047 & 921 \\
& Def2-SVPD & -9.246 & -10.166 & 921 \\
& Def2-SV(P) & -9.137 & -10.058 & 922 \\
& Def2-TZVP & -9.156 & -10.074 & 919 \\
\hline
\multirow{4}{*}{$\left| \mathrm{mono}_A\right>_{+3}$} & Def2-SVP & -12.031 & -12.551 & 520 \\
& Def2-SVPD & -12.119 & -12.645 & 525 \\
& Def2-SV(P) & -12.044 & -12.565 & 521 \\
& Def2-TZVP & -11.653 & -12.445 & 792 \\
\hline
\multirow{4}{*}{$\left| \mathrm{mono}_B\right>_{+2}$} & Def2-SVP & -9.190 & -10.206 & 1017 \\
& Def2-SVPD & -9.301 & -10.321 & 1020 \\
& Def2-SV(P) & -9.200 & -10.218 & 1017 \\
& Def2-TZVP & -9.208 & -10.227 & 1018 \\
\hline
\multirow{4}{*}{$\left| \mathrm{mono}_B\right>_{+3}$} & Def2-SVP & -12.031 & -12.668 & 637 \\
& Def2-SVPD & -12.114 & -12.749 & 635 \\
& Def2-SV(P) & -12.044 & -12.682 & 638 \\
& Def2-TZVP & -12.021 & -12.654 & 633 \\
\hline
\multirow{4}{*}{$\left| \mathrm{dimer}\right>_{+5}$} & Def2-SVP & -11.783 & -12.744 & 961 \\
& Def2-SVPD & -- & -- & -- \\
& Def2-SV(P) & -11.793 & -12.755 & 962 \\
& Def2-TZVP & -11.842 & -12.789 & 947 \\
\hline\hline
\end{tabular}
\caption{Energies of the HOMO and HOMO-1 molecular orbitals of pair \#2, in [eV], at the CAM-B3LYP level of theory. $\Delta \epsilon$ is the energy difference between the HOMO and the HOMO-1, in [meV].}
\label{table:pair2_BSet}
\end{table}

\begin{table}[htb!]
\centering
\begin{tabular}{l | c c}
\hline\hline
J & Pair \#1 & Pair \#2 \\
\hline
Def2-SVP & 870 & 959 \\
Def2-SVPD & 870 & -- \\
Def2-SV(P) & 483 & 960 \\
Def2-TZVP & 374 & 893 \\
\hline\hline
\end{tabular}
\caption{J in [meV].}
\label{table:J_BSet}
\end{table}

    \section{Molecular dynamics}
	\label{sec:si_md}
	To ensure equilibration of the systems, the following approach has been adopted. First, random starting configurations have been generated with the packmol program, and each one has been allowed to relax with an energy minimization followed by a molecular dynamics simulation of \unit[1]{ns} within the NVT ensemble at \unit[300]{K}. After that, a NPT simulation of \unit[2]{ns} at \unit[300]{K} and \unit[1]{atm} to relax the dimensions of the box. Then, the temperature of the system has been subjected to a heating up to \unit[1000]{K} in a \unit[1]{ns} run, a NVT run at \unit[1000]{K} for another \unit[1]{ns}, and a cooling down to \unit[300]{K} for \unit[10]{ns}. After that, the production run has been carried out in the NPT ensemble at \unit[300]{K} and \unit[1]{atm} for \unit[1]{ns}. To more thoroughly explore the configurational space, a total of ten different starting configurations generated by Packmol have been simulated through this procedure.

	\section{Hyperparameter optimization of correlation models}
	\label{sec:si_hyperparameter_optimization}

        Well performing model architectures were determined from a random grid search of selected hyperparameter for each model. Tab.~\ref{tab:hyperparam_ranges} summarizes the hyperparameters which were varied for each model and their respective ranges. Multi-layer perceptrons (applied to the full spectra and PCA contracted spectra) were set up with three hidden layers but varying number of neurons per layer. Activations for all model architectures were chosen to be either a leaky version of the ReLU function ($\alpha = 0.2$) or the softplus function.

        \begin{align*}
            \text{leaky\_ReLU}(x) &= \begin{cases} x, & \text{if } x > 0 \\ \alpha x, & \text{otherwise} \end{cases}  \\
            \text{softplus}(x)    &= \log(1 + \exp(x))
        \end{align*}

        Tab.~\ref{tab:hyperparam_ranges} also reports the set of hyperparameters for which the trained models achieved the highest prediction accuracies as determined from a 10-fold cross-validation protocol. Best performing sets of hyperparameters were determined from a random grid search with a total of 512 different models constructed for each model class. Training of individual models was aborted based on an early stopping criterion.

        \begin{table}[!ht]
            \centering
            \begin{tabular}{p{6cm} | p{3cm} | p{2.5cm}p{1.5cm}p{1.5cm}}
                    Model    & Hyperparameter & optimal & lower bound & upper bound \\
                \hline
                    Bayesian Linear Regression     & Regularization & $10^{-3.875}$ & $10^{-4}$ & $10^1$    \\
                             & Learning rate  & $10^{-4.75}$  & $10^{-5}$ & $10^{-1}$ \\
                \hline
                    Convolutional Neural Network  & Regularization    & $10^0$        & $10^{-4}$ & $10^1$    \\
                             & Learning rate     & $10^{-1.875}$ & $10^{-5}$ & $10^{-1}$ \\
                             & Filters           & $24$          & & \\
                             & Filter sizes      & $20, 50$      & & \\
                             & Filter activation & softplus      & \multicolumn{2}{c}{softplus or leaky ReLU} \\
                             & Neurons / layer   & $2$           & & \\
                             & Dense activation  & leaky ReLU    & \multicolumn{2}{c}{softplus or leaky ReLU} \\
                \hline
                    Multi-Layer Perceptron (full)     & Regularization    & $10^{-2}$     & $10^{-4}$ & $10^1$    \\
                             & Learning rate     & $10^{-3.5}$   & $10^{-5}$ & $10^{-1}$ \\
                             & Neurons / layer   & $91$          & $1$       & $100$     \\
                             & Activations       & softplus      & \multicolumn{2}{c}{softplus or leaky ReLU} \\
                \hline
                    Multi-Layer Perceptron (pca)     & Regularization    & $10^{-3.5}$   & $10^{-4}$ & $10^1$    \\
                             & Learning rate     & $10^{-1.5}$   & $10^{-5}$ & $10^{-1}$ \\
                             & Neurons / layer   & $41$          & $1$       & $100$     \\
                             & Activations       & leaky ReLU    & \multicolumn{2}{c}{softplus or leaky ReLU} \\
            \end{tabular}
            \caption{Hyperparameter ranges and sets of hyperparameters for which trained models achieved the highest prediction accuracies. Best performing sets of hyperparameters were determined from a random grid search. }
            \label{tab:hyperparam_ranges}
        \end{table}

    \section{Sampling efficiency}
    \label{sec:si_sampling_efficiency}

        The predictive power of a regression model depends on the how representative the training set is for the underlying (unknown) physical correlations. Typically, larger training sets with more examples of absorption spectra and associated couplings yield higher prediction accuracies. To asses how many examples are needed for sufficient prediction accuracies we train the introduced regression models on training sets of different sizes. More specifically, we employ the Convolutional Neural Network  model (see main text for details) with the determined best performing set of hyperparameters. The training set of 1,250 examples is reduced by randomly selected examples. Following this procedure, we generate a total of $20$ training subsets of different sizes. Convolutional Neural Network  models are then fully trained on their respective training set and their predictive powers are assessed via test set predictions. Note, that the test set is identical for all predictions. Results are reported in Fig.~\ref{fig:sampling_efficiency}

        \begin{figure}[!ht]
            \centering
            \includegraphics[width = 1.0\textwidth]{sample_size_prediction_accuracy.png}
            \caption{Sampling efficiency for Bayesian convolutional neural networks trained to predict tight-binding coupling strengths from Vis/NIR spectra. Colors indicate poor predictions (purple) and accurate predictions (green). }
            \label{fig:sampling_efficiency}
        \end{figure}

        We observe a steady increase in the average prediction accuracy of the Convolutional Neural Network  model when increasing the size of the training set from $100$ examples up to about $850$ examples. From there on, no significant change in the average prediction accuracy can be observed with a further increase of the training set size.

    \section{Prediction accuracies for coupling strengths computed with the Kohn-Sham orbital formalism}
    \label{sec:si_cross_predictions}

        ML models can, at best, identify statistical correlations between the features and the targets to which they are exposed during the training process. As such, the positive correlation between predicted and computed coupling strengths obtained from the tight-binding formalism does not necessarily indicate that the ML models identified relevant physical correlations, and might just be due to statistical correlations based on the impreciseness of the methods used to compute Vis/NIR spectra and/or coupling strengths. \\

        We suggest to test if the ML models in the main text only identified statistical correlations (and not the relevant physical correlations) by constructing additional models trained to predict coupling strengths computed with the Kohn-Sham orbital formalism. Similar to the ML models reported in the main text, we run a full hyperparameter optimization for the Bayesian CNN models to predict Kohn-Sham orbital formalism coupling strengths. Test set predictions of the best performing model as identified from the hyperparameter search are reported in Fig.~\ref{fig:sampling_efficiency}. We observe, that the trained Bayesian CNN is capable of identifying correlations between the Vis/NIR spectra and the coupling strengths even for this case where coupling strengths are computed based on the Kohn-Sham orbital formalism. Moreover, we observe similar prediction accuracies despite significantly different coupling strength distributions (see section \ref{sec:si_coupling_strength_distributions}).

        \begin{figure}[!ht]
            \centering
            \includegraphics[width = 0.5\textwidth]{pascal_hyperparam_optimized.png}
            \caption{Predictions for coupling strengths computed based on the Kohn-Sham orbital formalism. Coupling strengths are predicted from a Bayesian convolutional neural network for which hyperparameters have been optimized in a random grid search. }
            \label{fig:kohn_sham_predictions}
        \end{figure}

        If the predictive power of the Bayesian CNN reported in the main text relied on spurious statistical correlations between the Vis/NIR spectra and the tight-binding coupling strengths, the model presented in this section no longer has the opportunity to exploit these spurious correlations due to the change in formalism for computing coupling strengths. However, we still observe relatively high prediction accuracies for Kohn-Sham orbital couplings despite the change in the range of the couplings and the change in their distribution. We conclude that the predictive power of the Bayesian CNN model presented in Fig.~\ref{fig:kohn_sham_predictions} must either be due to the model identifying the relevant physical correlations, or other statistical correlations arising from the Kohn-Sham orbital couplings. \\

        To rule out the possibility that the Bayesian CNN presented in Fig.~\ref{fig:kohn_sham_predictions} identified other spurious statistical correlations yielding similar prediction accuracies we propose the construction of a hybrid dataset. This hybrid dataset is constructed by randomly choosing half of the coupling strengths computed with the tight-binding formalism, and the other half with the Kohn-Sham orbital formalism. Note, that the differences in the ranges of the coupling strengths are accounted for by standardizing coupling strengths based on the characteristics of each individual dataset following

        \begin{align}
            j = \frac{J - \langle J \rangle_\text{train}}{\sqrt{ \langle \left( J - \langle J \rangle_\text{train} \right)^2 \rangle_\text{train} }},
        \end{align}

        where $\langle \cdot \rangle_\text{train}$ denotes the average over the training set.

        \begin{figure}[!ht]
            \centering
            \includegraphics[width = 0.5\textwidth]{pedot_hybrid_test_predictions.png}
            \caption{Coupling strength predictions obtained from a model trained on a hybrid dataset assembled from both formalisms (equal contributions). Note, that coupling strengths have been standardized before training the model, and the model predicted coupling strengths in standardized units (shown in the plot). }
            \label{fig:kohn_sham_predictions}
        \end{figure}

        We find that the proposed Bayesian CNN architecture yields prediction accuracies comparable to prior experiments with coupling strengths computed from one or the other formalism. We interpret this observation as an indicator showing that the model does not solely rely on spurious statistical correlations caused by inaccuracies of the physical models used to compute the coupling strengths.

    \section{Distribution of electronic coupling}
    \label{sec:si_coupling_strength_distributions}
        \begin{figure}[htb!]
         \centering
         \includegraphics[width=0.7\linewidth]{J_g09_log_color.png}
         \includegraphics[width=0.7\linewidth]{J_turbomole_log_color.png}
        \caption{Distribution of electronic couplings in PEDOT:PEDOT complexes. The total number of complexes is 1,420. (A) is computed with the tight binding formalism model, and (B) with the orbital overlap method.}
        \label{fig:couplings}
        \end{figure}


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

Blazars,  some of the most powerful Active Galactic Nuclei (AGNs), have a relativistic jet pointing toward the observer \citep[e.g.,][]{abdo01,fra} and show rapid variability and high optical and radio polarization.  Such objects are the most numerous class of extragalactic sources detected by TeV telescopes, the most sensitive of which are the Imaging Atmospheric Cherenkov Telescopes (IACTs) such as  MAGIC\footnote{\url{https://magic.mpp.mpg.de/}}, H.E.S.S.\footnote{\url{https://www.mpi-hd.mpg.de/hfm/HESS/}},  VERITAS\footnote{\url{https://veritas.sao.arizona.edu/}},  and the upcoming Cherenkov Telescope Array (CTA)\footnote{\url{https://www.cta-observatory.org/}}. Despite their high sensitivity, however, observations by current IACTs are limited by their fairly small fields of view ($\sim$ few degrees), weather conditions, the need for relatively dark night skies, and by a high background that requires fairly long observations. A source with a flux of $\sim$1$\%$
of the Crab Nebula flux  requires around 50 hours of observation time for a detection at 5$\sigma$.  IACTs typically take data for only about 1200 hours per year \citep{Naurois}.
Those constraints provide a strong incentive to identify likely targets for IACT observations.

Blazar Spectral Energy Distributions (SEDs) show two
broad peaks in a $\nu f{_\nu}$ representation. The low-energy peak is attributed to synchrotron radiation, and the high-energy one is usually thought to be produced by inverse Compton radiation (IC) \citep[e.g.][]{Sikora}.  Based on the position of the  synchrotron peak ($\nu^{S}_{peak}$) in the SED, blazars are divided into three subclasses:   low-synchrotron-peaked (LSP,  with $\nu^{S}_{peak} \leq 10^{14}\,$Hz), intermediate-synchrotron-peaked (ISP,  with $10^{14}\, {\rm Hz} < \nu^{S}_{peak} \leq 10^{15}\,$Hz) and high-synchrotron-peaked (HSP, with $\nu^{S}_{peak} > 10^{15}\,$Hz) \citep{abdo01}.
HSPs, primarily BL Lac objects, represent the most numerous class of extragalactic TeV-energy sources. The TeVCat\footnote{\url{http://tevcat.uchicago.edu/}} is an online, interactive catalog for very-high-energy
(VHE energies, $E > 100$ GeV) $\gamma$-ray astronomy \citep{horan}.
This catalog reports 223 TeV sources as of this writing.  Among the 61 objects associated with blazars,  51 of them are HSPs and only 10 are LSP/ISP flat-spectrum radio quasars (FSRQs)\footnote{The rest of the sources  in TeVCat are Galactic sources or of unidentified nature.}.

All-sky observations with the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi) \citep{FLAT} at GeV  energies offer opportunities to find such targets.  An example is the Third Catalog of Hard {\it Fermi}-LAT Sources \citep[3FHL:][]{3fhl}, which reports the locations and spectra of sources significantly detected in the 10 GeV - 2 TeV energy range during the first 7 years of the {\it Fermi} mission.
From the 3FHL
it is possible to select  TeV candidates by $\gamma$-ray flux and photon index.

An alternative approach to searching for TeV candidates is to find objects belonging to a class of sources likely to be seen at TeV energies. In the case of blazars, this can be done by identifying those objects whose synchrotron emission peaks at high frequencies. An example of this approach is the second {\it WISE} High Synchrotron Peak Catalog (2WHSP) \citep{2wHSP}, which is a list of HSP candidates based on  multi-frequency analysis of $\gamma$-ray source candidates with $|b| > 10^\circ$

Here we present a third approach to search for TeV HSP source candidates. This method includes two steps: (1) we use $\gamma$-ray variability information to search out potential HSPs among the unclassified {\it Fermi}-LAT sources; and (2) we analyze $\gamma$-ray spectra of these sources using more {\it Fermi}-LAT data than used in published catalogs.
The starting point is the third {\it Fermi}-LAT all-sky catalog  of sources detected at energies between 100 MeV and 300 GeV \citep[3FGL:][]{3FGL}.
The 3FGL catalog lists 3033 $\gamma$-ray sources, of which 1745 are AGNs, mostly BL Lacs and FSRQs, and includes $\gamma$-ray source locations, energy spectra, variability information on monthly time scales, and likely associations with objects seen at other wavelengths. In this catalog 573 sources are listed as  blazars of uncertain type (BCUs) and 1010 objects lack  a plausible counterpart at other wavelengths (Unassociated Catalog Sources, UCSs)\footnote{A preliminary version of the 4FGL catalog is available at \url{https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/}, but it does not include the variability information needed for this analysis.}.

Although BCUs and UCSs often lack optical spectra and sufficient information for a rigorous classification,  statistical methods such as the Artificial Neural Network (ANN) algorithm can potentially provide  classifications of these sources \citep[e.g.,][]{bflap,pablo, zoo, leu}. In particular,  \citet{pablo} found 559 of the UCS sources have characteristics similar to those of AGN.  These UCS$_{agns}$ are combined with the original 573 BCUs  from the 3FGL catalog to provide the targets for our search
for HSP/TeV candidates.

The paper is organized as follows: in Sect.~\ref{<2>} we present the machine learning method used in this study; in Sect.~\ref{<3>} we describe the selection of HSP candidates among the uncertain 3FGL objects; and in Sect.~\ref{<4>}  we discuss the results of a dedicated {\it{Fermi}}-LAT analysis of the sources found analyzing 104 months of data. In Sect.~\ref{<5>} we examine the detectability of the targets by the present generation of IACTs and CTA. In Sect.~\ref{<6>} we summarize the conclusions of this study.

\section{The ANN  search method}
\label{<2>}
The starting point for selecting HSP candidates is the ANN method previously applied to {\it Fermi}-LAT sources to distinguish FSRQ-like sources from those with BL Lac characteristics \citep{bflap, zoo}.  The key idea is that the $\gamma$-ray flares of BL Lacs tend to be smaller and less frequent than those of FSRQs. The input is the empirical cumulative distribution function (ECDF) of monthly $\gamma$-ray flux values for each source, taken from the 3FGL. We included in the algorithm values corresponding to increments of 10\% from the 10th to the 100th percentile of the ECDF. The  algorithm computes a likelihood value arranged to have two possibilities: class $A$ or class $B$, with a likelihood ({\it{L}}) assigned to each analyzed source so that the likelihoods to belong to one or the other of the two classes are related by $L_A = 1-L_B$. In this way, the greater the value of {\it{L}$_{A}$}, the greater the likelihood that the source is a class $A$ candidate. In this case, the likelihood applies to the source having HSP characteristics, $L_{HSP}$. This approach uses the two-layer-perceptron ANN technique \citep{gis90,bio}, which is probably the most widely used architecture for practical applications of neural networks.  Data enter the neural network through  nodes in the input layer. The
information travels across the links and is processed in the nodes through an activation function. Each node in
the output layer returns the likelihood of a source to be a specific
class.
We applied the algorithm to the three  synchrotron peak subclasses as classified in the Third Catalog of Active Galactic Nuclei detected by the {\it Fermi} LAT \citep[3LAC:][]{3lac} in order to train it to distinguish each source class.

Repeating the analysis of \citet{bflap} as a cross-check on that analysis, we considered 289 HSPs and 824 non-HSP objects classified in the third {\it{Fermi}}-LAT AGN catalog \citep[3LAC;][]{3lac}. Maintaining the same ratio as in the catalog, that is, one third  HSPs and two thirds non-HSP sources,   we randomly mixed the sample and divided it  into 3 subsamples:  training,  validation, and testing. The training sample is used to optimize the network. The validation sample is used to avoid over-fitting. The testing sample is independent both of the training and validation ones and was used to monitor
the accuracy of the network.
Although the random sampling resulted in a different training set from the previous analysis, the results were the same:  for  $L_{HSP}$ $>$ 0.8, 75\% of the sources have characteristics of HSPs, while for $L_{HSP}$ $>$ 0.89, we expect 90\% of the sources to be HSP-like.

\section{Identifying HSP candidates}
\label{<3>}

We then applied the
 ANN HSP algorithm to the 573 BCUs and the 559 UCS$_{agn}$ of our sample.
 The resulting likelihood distributions plotted in Fig.~\ref{can} show, as expected, a  peak at $L_{HSP}$ = 0.0 due to the non-HSP populations (ISP and LSP), which are much more numerous than HSPs at energies covered by the {\it Fermi} LAT.
As for the original analysis \citep{bflap}, the lack of a peak at $L_{HSP}$ = 1.0
indicates that the ANN network was not able to separate HSPs cleanly, but for the purpose of selecting candidates for additional analysis we are primarily interested in finding a high fraction of the sources with the desired characteristics.
Requiring  the $L_{HSP}$ $>$ 0.8 value, we identified 48 BCU and  32 UCS$_{agn}$ as HSP candidates.
Applying the higher threshold value $L_{HSP}$ $>$ 0.89,
we identified  11 BCUs and 5 UCS$_{agn}$ as Very High Confidence (VHC) HSP candidates.
Table 1 and Table 2 show the full lists of candidates.  The names of the VHC candidates are shown in bold.

 We compared the results of our analysis with the list of blazars presented in the 2WHSP catalog \citep{2wHSP}, since both methods attempt to identify HSP objects.  Of the 11 VHC candidates from the BCU list, 6 were also identified by 2WHSP: 3FGL J0506.9$-$5435, 3FGL J0921.0$-$2258, 3FGL J1155.4$-$3417, 3FGL J1711.6+8846, 3FGL J1714.1$-$2029, and 3FGL J1944.1$-$4523.  None of the 5 VHC candidates from the UCS list were found in the 2WHSP catalog.  Similar fractions were found for the other parts of the lists, indicating that the two methods are complementary in finding candidate HSPs.

\begin{figure*}[!]
\includegraphics[width=0.5 \textwidth]{HBLbcuLog.pdf}
\includegraphics[width=0.48 \textwidth]{HBLucsLog.pdf}
\caption{\label{can} ANN likelihood to be HSP candidates of 3FGL BCUs (left) and UCS$_{agn}$  (right).
Blue bars: sources with  $0.89 >L_{HSP} > 0.8$; Red bars: VHC candidates with $L_{HSP} > 0.89$}
\end{figure*}

\section {Fermi-LAT Spectral Analysis}
\label{<4>}
 Because the 4FGL catalog will soon be available with a larger sample of sources to study, we chose to focus our spectral analysis for this new method of selecting HSP candidates on the 16 VHC sources.
We analyzed 104 months of {\it Fermi}-LAT Pass 8 \citep{P8} data
from 2008 August 4 to 2017 April 4, selecting $\gamma$-ray events in the energy range $E=[0.1,1000]$ GeV, passing standard data quality selection criteria and zenith angle cuts for AGN \citep[e.g.][]{meyer},
in order to find the $\gamma$-ray properties of our HSP candidates.
We considered events belonging to the Pass~8 {\tt SOURCE} event class and used the corresponding  instrument response functions {\tt P8R2\_SOURCE\_V6}, since we were interested in point source detection.
We used the interstellar emission model (IEM) released with Pass 8 data \citep{Acero:2016qlg} (i.e., {\tt gll\_iem\_v06.fits}).
This is the model recommended for use with Pass 8 analyses.
We also included the standard template for the isotropic emission ({\tt iso\_P8R2\_SOURCE\_V6\_v06.txt})\footnote{For descriptions of these templates, see \url{http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html}.}.

We developed an analysis pipeline using {\tt FermiPy}, a Python package that automates analyses with the {\it Fermi} Science Tools \citep{2017arXiv170709551W}\footnote{See \url{http://fermipy.readthedocs.io/en/latest/}.}.
{\tt FermiPy} includes tools that can 1) generate simulations of the $\gamma$-ray sky, 2) detect point sources, and 3) calculate the characteristics of their SEDs.

The likelihood analysis works on a square region of interest (ROI). We used a $16^{\circ}\times16^{\circ}$ ROI centered on the sources of our sample.
We analyzed each ROI separately.
In each ROI, we binned the data with a pixel size of $0.08^{\circ}$ and 8 energy bins per decade.
Our background model included the IEM, isotropic template and sources from the preliminary 8-year list, FL8Y\footnote{\url{https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/gll_psc_8year_v5.fit}}, except for the source being analyzed.
We allowed the normalization and slope of the IEM and the normalization of the isotropic template  to vary.
We first relocalized the source of interest, and then we searched for new point sources with Test Statistic $TS >25$, defined as twice the difference between the log-likelihood for the null hypothesis (no source) and the hypothesis of a source at the location.  New sources were then added to the analysis. We then analyzed each of our HSP candidates, assuming a power-law spectrum after determining that all the VHC candidates used that spectral form in the 3FGL catalog \citep{3FGL}. For the individual energy bins, we plotted as upper limits those points with a $TS <9$.

In Table 1 and Table 2 we report the following parameters for the HSP candidates: the best fit and $1\sigma$ error of the position, the $TS$ of detection,  the photon index found for a power-law SED shape, the flux, and the value of $L_{HSP}$.
The photon index parameter, if less than 2, can be a relevant indicator for an IC peak at TeV energies and therefore quite useful in selecting IACT and CTA  candidates. The mean and rms of the photon indexes of HSPs are $1.87\pm0.20$ while for LSPs and ISPs these are $2.21\pm0.18$, $2.07\pm0.20$ respectively \citep{3lac}.

We have repeated the analysis with a log-parabola intrinsic spectrum instead of a power law. Only for three sources (3FGL~J0921.0-2258, 3FGL~J0153.4+7114, and 3FGL~J0506.9-5435) did we find that a log-parabola is preferred at the 4\,$\sigma$ level, i.e. $\mathrm{TS}_\mathrm{curv} > 16$, where $\mathrm{TS}_\mathrm{curv}$ is twice the difference of the log-likelihood values of the best fit with a log-parabola and a power law, respectively.
These results are largely consistent with the spectral parameters listed in the preliminary 4FGL catalog\footnote{\url{https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/}}.  The 4FGL sources corresponding to 3FGL J0921.0$-$2258 and 3FGL J0153.4+7114 have power-law spectra instead of the curved spectra we found, although the evidence for curvature is of low significance in our analysis. If we had used the power-law fit, the extrapolations would have obviously been higher.

\section{TeV candidates}
\label{<5>}

In this section, we compare the extrapolated fluxes of the {\it Fermi}-LAT spectra of our VHC sources to the sensitivity of present IACTs and the future CTA.
In the case that the log-parabola fit to the {\it Fermi}-LAT data is preferred with ${\rm TS}_{\rm curv} > 16$, we use the log-parabola parametrization for the extrapolation and the best-fit power law otherwise.
In order to evaluate whether the VHC HSP candidates can realistically be observed with IACTs or CTA, we must take into account the interaction of $\gamma$-rays with photons of  the extragalactic background light (EBL). The EBL spans the wavelength regime from ultraviolet to far-infrared wavelengths and mainly consists of the integrated starlight emitted over the history of the Universe and starlight  absorbed and re-emitted by dust in galaxies \citep{hauser,kashlinsky}.
During the propagation of  $\gamma$ rays through the EBL, the electron-positron pairs produced via $\gamma\gamma\to e^+e^-$
leads to an attenuation of the initial $\gamma$-ray flux \citep{nikisov,gould,dwe}.
To properly evaluate the absorption effect of the EBL it is necessary to know the redshift of the analyzed $\gamma$-ray source. Since the redshifts of the selected VHC HSP candidates are unknown, we assume redshifts  of $z = 0$ and $z = 0.5$ for the calculation. These are typical values of observed $\gamma$-ray BL Lacs\footnote{We recognize that a redshift identically 0 is unrealistic. This choice provides a ``best case'' limit.}.
We use these $z$ values while acknowledging that the blazar redshift range is a very open and long-standing debate. Some authors argue that the BL Lacs without a redshift are likely much more distant than those with a measured one \citep[e.g.,][]{pad}, so it could be possible that some objects fall beyond the 0.5 redshift.

Using the EBL model of \citet{dominguez}, we extrapolate the best-fit spectra obtained with the \emph{Fermi}-LAT analysis in Sec.~\ref{<4>} up to 10\,TeV with the assumed redshift values. The results are shown in Fig.~\ref{bcu} and Fig.~\ref{ucscta} where we compare the extrapolated fluxes with the CTA sensitivity~\citep{ctasens} for 50\,hours (5\,hours) of observations as a solid (dashed) gray line.
The sensitivities for currently operating IACT arrays H.E.S.S., MAGIC, and VERITAS are shown as a red band \citep{hessperf,magicperf}\footnote{The VERITAS sensitivity curve from \url{https://veritas.sao.arizona.edu/about-veritas-mainmenu-81/veritas-specifications-mainmenu-111} is used.}.
The CTA sensitivity curves are available for the northern and southern array and for zenith angles ($Z$) of $20^\circ$ and $40^\circ$. We choose the sensitivity curve depending on the source declination $\delta$ assuming CTA site locations at 28.7${}^\circ$ northern latitude and 24.7${}^\circ$ southern latitude. For $\delta > 58.7^\circ$,  we choose the northern array with $Z = 40^\circ$, for $58.7^\circ \leqslant \delta < 2.7^\circ$ the northern array with $Z=20^\circ$, for $2.7 \geqslant Z \geqslant -54.7^\circ$ the southern array with $Z = 20^\circ$, and the southern array with $Z=40^\circ$ for $\delta < -54.7^\circ$\footnote{The sensitivity curves of the northern and southern array  are available at \url{www.cta-observatory.org}}.

 In Table~\ref{tab:zmax},  we report an estimate of the maximum redshift values of the VHC TeV candidates so that the sources are still detectable at 5\,$\sigma$ for 50-hr observations with current IACTs as well as 5- and 50-hr CTA  observations.
 For these estimates, we assume that the spectra are detectable if the extrapolation is above the sensitivity curves.
 We also report how these values change if the best-fit spectral parameters are changed by their $1~\sigma$ statistical uncertainties.
 An actual detection will depend on the spectral shape and the exact exposure time, and therefore these numbers should be regarded as a rough estimate.
 If no value is given, the source might not be significantly detected with the assumed observation within the assumed redshift range.

 Any of the VHC sources that lie at low redshift should be detectable with CTA, with the exception of  3FGL~J0921.0$-$2258, which shows strong intrinsic curvature.
 Similarly, most sources should already be detectable by the present-generation IACTs if their redshifts are close to 0.
 For a redshift of $z = 0.5$, most sources appear to be good candidates for detection with 50 hours of CTA observations, primarily in the 100-GeV energy range.  Only 3FGL~J1714.1$-$2029 seems possible as a candidate for current IACTs if its distance is at the upper end of the range considered, due to the EBL attenuation.

\begin{figure}[p]
\centering
\includegraphics[width=0.9\textwidth]{spectra_VHC_BCU_LP-crop.pdf}
\caption{\label{bcu}SEDs of the VHC BCU sources as TeV candidates. The dashed (dotted) line denotes the extrapolation of the best-fit {\it Fermi}-LAT spectra up to 10\,TeV for a redshift of $z = 0$ ($z = 0.5$) The shaded region indicates the possible source flux for redshifts between $0 < z \leq 0.5$. The CTA  sensitivity for 5 (50) hours of observation is shown as a gray dashed (solid) line. The sensitivity curve for either the northern or southern array is used as described in the main text. The 50 hour sensitivity for currently operating IACTs is shown as a red shaded band.}
\end{figure}

\begin{figure}[p]
\centering
\includegraphics[width=0.9\textwidth]{spectra_VHC_UCS_LP-crop.pdf}
\caption{\label{ucscta}SEDs of the VHC UCS sources as TeV candidates. The dashed (dotted) line denotes the extrapolation of the best-fit {\it Fermi}-LAT spectra up to 10\,TeV for a redshift of $z = 0$ ($z = 0.5$) The shaded region indicates the possible source flux for redshifts between $0 < z \leq 0.5$. The CTA  sensitivity for 5 (50) hours of observation is shown as a gray dashed (solid) line. The sensitivity curves for either the northern of southern array is used as described in the main text. The 50 hour sensitivity for currently operating IACTs is shown as a red shaded band. }
\end{figure}

\section{Conclusions}
\label{<6>}
Motivated to expand the sample of TeV sources,
we applied a machine learning algorithm  to variability parameters of {\it {Fermi}}-LAT blazar-like sources without firm identifications combined with new analysis of the LAT data for these sources.  Follow-up work will require additional multiwavelength studies, including finding redshifts for most of the candidates and targeted observation by IACTs.  We also recognize that this search is necessarily incomplete because of the difficulty to distinguish the blazar subclasses  by the $\gamma$-ray flux properties only. As already pointed out by \citet{3lac}, the $\gamma$-ray sources with unknown properties are generally fainter than the well-defined classes. The fainter sources offer less of the flaring information needed for the machine learning method, and so there may be HSP blazars among the sample of 3FGL sources rejected in the first step of our method. The level of incompleteness is difficult to quantify due to the very similar values of the synchrotron peaks of the three blazar subclasses. Nevertheless, some VHC HSP candidates, also thanks to the analysis of {\it {Fermi}}-LAT data, are convincing as TeV candidates and should be promising targets for currently operating IACTs, especially if the sources are located below reshifts of $\sim 0.1$. If the sources are further away, CTA should be capable of significantly detecting them.

\acknowledgments

The \textit{Fermi} LAT Collaboration acknowledges generous ongoing support
from a number of agencies and institutes that have supported both the
development and the operation of the LAT as well as scientific data analysis.
These include the National Aeronautics and Space Administration and the
Department of Energy in the United States, the Commissariat \`a l'Energie Atomique
and the Centre National de la Recherche Scientifique / Institut National de Physique
Nucl\'eaire et de Physique des Particules in France, the Agenzia Spaziale Italiana
and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education,
Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research
Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and
the K.~A.~Wallenberg Foundation, the Swedish Research Council and the
Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase
is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre
National d'\'Etudes Spatiales in France. This work performed in part under DOE
Contract DE-AC02-76SF00515.

MDM acknowledges support by the NASA {\it Fermi} Guest Investigator Program 2014 through the {\it Fermi} multi-year Large Program No. 81303 (P.I. E.~Charles) and by the NASA {\it Fermi} Guest Investigator Program 2016 through the {\it Fermi} one-year Program No. 91245 (P.I. M.~Di Mauro).

This  paper  has  gone  through  internal  review  by  the CTA  Consortium.
This research has made use of the CTA instrument response functions provided by the CTA Consortium and Observatory, see \url{http://www.cta-observatory.org/science/cta-performance/} (version prod3b-v1) for more details.
\vspace{5mm}
\facilities{{\it Fermi}, H.E.S.S., MAGIC, VERITAS, CTA}

\begin{table}
\clearpage
\begin{small}
\caption{\label{table1}Results of the {\it Fermi}-LAT analysis of the full list of BCU HSP candidates. Names of the high-confidence sources are shown in {\bf bold}. Columns: (1) Source name; (2) RA; (3) DEC; (4) $68\%$ error on the position $\theta$;  (5) Detection TS 0.1-- 300 GeV;  (6)   Photon Index; (7)   Integrated flux in the energy range 0.1-- 300 GeV;  8) $L_{HSP}$, likelihood the source is an HSP   }
\begin{tabular}{lcccccccc}
\hline
\hline
Source name  & RA   &  DEC  &  $\theta_{68\%}$	& TS  & Ph.Index   &  Flux &  $L_{HSP}$\\
\hline
            & [deg] & [deg]  & [deg]  &  &   & [$10^{-9}$ ph/cm$^2$/s] & \\
\hline
\hline
3FGL J0030.2$-$1646  &   7.59  &   -16.82   &   0.02   &   168.7 &   $1.66\pm0.08$   &   $13.4\pm3.3$ & 0.80  \\
3FGL J0039.0$-$2218  &   9.77  &   -22.32   &   0.03   &   89.3 &   $1.67\pm0.11$   &   $9.2\pm2.9$ &  0.86 \\
3FGL J0040.3+4049  &   10.09  &   40.83   &   0.03   &   75.9 &   $1.92\pm0.16$   &   $18.2\pm9.4$ & 0.87  \\
3FGL J0043.5$-$0444  &   10.88  &   -4.72   &   0.04   &   54.0 &   $1.91\pm0.15$   &   $16.2\pm7.2$ & 0.83  \\
3FGL J0043.7$-$1117  &   10.94  &   -11.31   &   0.04   &   69.4 &   $1.86\pm0.12$   &   $16.0\pm5.3$ & 0.88  \\
{\bf 3FGL J0047.9+5447}  &   12.02  &   54.81   &   0.03   &   56.7 &   $1.58\pm0.17$   &   $4.9\pm3.2$ & 0.92  \\
3FGL J0132.5$-$0802  &   23.19  &   -8.07   &   0.03   &   71.9 &   $1.87\pm0.11$   &   $16.8\pm5.4$ & 0.84  \\
{\bf 3FGL J0153.4+7114}  &   28.43  &   71.26   &   0.02   &   80.9 &   $1.82\pm0.13$   &   $17.5\pm7.3$ &  0.89 \\
3FGL J0204.2+2420  &   31.09  &   24.27   &   0.04   &   27.6 &   $1.70\pm0.16$   &   $4.7\pm2.6$ &  0.81 \\
3FGL J0305.2$-$1607  &   46.29  &   -16.14   &   0.02   &   147.6 &   $1.80\pm0.11$   &   $17.8\pm5.8$ &  0.86 \\
3FGL J0342.6$-$3006  &   55.71  &   -30.11   &   0.04   &   43.2 &   $1.96\pm0.14$   &   $12.5\pm4.7$ & 0.84  \\
3FGL J0439.6$-$3159  &   69.85  &   -32.03   &   0.03   &   119.9 &   $1.75\pm0.09$   &   $13.3\pm3.9$ & 0.81  \\
{\bf 3FGL J0506.9$-$5435}  &   76.76  &   -54.60   &   0.01   &   455.4 &   $1.50\pm0.05$   &   $14.2\pm2.2$ & 0.89  \\
3FGL J0515.5$-$0123  &   78.87  &   -1.42   &   0.04   &   45.7 &   $1.80\pm0.11$   &   $11.4\pm4.4$ & 0.85  \\
3FGL J0528.3+1815  &   82.11  &   18.27   &   0.04   &   35.7 &   $1.67\pm0.15$   &   $6.6\pm3.6$ & 0.87  \\
3FGL J0620.4+2644  &   95.17  &   26.74   &   0.02   &   92.0 &   $1.54\pm0.11$   &   $6.3\pm2.8$ & 0.85  \\
3FGL J0640.0$-$1252  &   100.01  &   -12.90   &   0.02   &   174.1 &   $1.52\pm0.09$   &   $10.3\pm3.4$ &  0.85 \\
3FGL J0646.4$-$5452  &   101.62  &   -54.92   &   0.03   &   190.3 &   $1.46\pm0.29$   &   $8.8\pm1.8$ & 0.87  \\
{\bf 3FGL J0648.1+1606}  &   102.03  &   16.09   &   0.03   &   40.1 &   $1.82\pm0.16$   &   $10.7\pm5.9$ & 0.90  \\
3FGL J0650.5+2055  &   102.64  &   20.93   &   0.02   &   206.2 &   $1.72\pm0.08$   &   $21.9\pm5.7$ & 0.82  \\
3FGL J0733.5+5153  &   113.35  &   51.86   &   0.03   &   104.3 &   $1.69\pm0.10$   &   $9.9\pm3.2$ &  0.85 \\
3FGL J0742.4$-$8133c  &   115.45  &   -81.54   &   0.05   &   32.3 &   $2.03\pm0.28$   &   $21.1\pm15.4$ &  0.88 \\
3FGL J0746.9+8511  &   117.25  &   85.22   &   0.03   &   119.0 &   $1.68\pm0.09$   &   $10.0\pm2.8$ & 0.83  \\
{\bf 3FGL J0921.0$-$2258}  &   140.24  &   -22.95   &   0.03   &   62.5 &   $1.74\pm0.14$   &   $9.4\pm4.1$ & 0.91  \\
3FGL J1040.8+1342  &   160.26  &   13.72   &   0.03   &   69.1 &   $1.71\pm0.13$   &   $8.3\pm3.4$ & 0.86  \\
3FGL J1141.2+6805  &   175.33  &   68.08   &   0.02   &   140.1 &   $1.69\pm0.09$   &   $10.9\pm2.8$ & 0.85  \\
{\bf 3FGL J1155.4$-$3417}  &   178.87  &   -34.33   &   0.02   &   147.3 &   $1.64\pm0.09$   &   $11.8\pm3.3$ & 0.92  \\
3FGL J1158.9+0818  &   179.71  &   8.31   &   0.04   &   51.5 &   $1.81\pm0.14$  &   $11.0\pm4.6$ & 0.80  \\
3FGL J1203.5$-$3925  &   180.85  &   -39.42   &   0.03   &   103.2 &   $1.70\pm0.10$   &   $13.5\pm4.5$ & 0.85  \\
3FGL J1319.6+7759  &   199.95  &   78.01   &   0.02   &   182.6 &   $1.95\pm0.8$   &   $28.3\pm5.9$ & 0.82  \\
{\bf 3FGL J1434.6+6640}  &   218.72  &   66.67   &   0.03   &   73.9 &   $1.58\pm0.12$   &   $4.4\pm1.7$ & 0.92  \\
3FGL J1446.8$-$1831  &   221.75  &   -18.51   &   0.05   &   27.9 &   $1.71\pm0.15$   &   $6.1\pm3.3$ &  0.84 \\
3FGL J1547.1$-$2801  &   236.81  &   -28.04   &   0.03   &   96.8 &   $1.78\pm0.10$   &   $19.0\pm6.1$ & 0.81  \\
3FGL J1612.4$-$3100  &   243.10  &   -30.99   &   0.02   &   495.0 &   $1.86\pm0.08$   &   $38.0\pm7.8$ & 0.81  \\
{\bf 3FGL J1714.1$-$2029}  &   258.52  &   -20.48   &   0.03   &   73.8 &   $1.44\pm0.12$   &   $5.1\pm2.3$ & 0.90  \\
{\bf 3FGL J1711.6+8846}  &   258.67  &   88.75   &   0.04   &   44.3 &   $1.83\pm0.15$   &   $8.8\pm4.2$ & 0.90  \\
3FGL J1824.4+4310  &   276.12  &   43.18   &   0.03   &   80.9 &   $1.83\pm0.15$   &   $13.6\pm5.3$ & 0.88  \\
3FGL J1841.2+2910  &   280.36  &   29.16   &   0.02   &   195.9 &   $1.80\pm0.08$   &   $29.0\pm7.1$ & 0.80  \\
3FGL J1855.1$-$6008  &   283.67  &   -60.13   &   0.06   &   21.4 &   $1.84\pm0.16$   &   $7.2\pm3.9$ &  0.84 \\
3FGL J1908.8$-$0130  &   287.20  &   -1.53   &   0.02   &   306.4 &   $1.52\pm0.21$   &   $18.6\pm2.8$ & 0.82  \\
{\bf 3FGL J1910.8+2855}  &   287.71  &   28.94   &   0.02   &   102.3 &   $1.62\pm0.10$   &   $9.8\pm3.3$ &  0.90 \\
3FGL J1939.6$-$4925  &   294.96  &   -49.47   &   0.03   &   64.6 &   $1.85\pm0.11$   &   $14.9\pm5.0$ &  0.85 \\
{\bf 3FGL J1944.1$-$4523}  &   296.11  &   -45.38   &   0.02   &   100.7 &   $1.64\pm0.10$   &   $9.4\pm3.3$ & 0.89  \\
3FGL J1959.8$-$4725  &   299.94  &   -47.43   &   0.01   &   923.8 &   $1.52\pm0.08$   &   $30.5\pm5.0$ & 0.87  \\
3FGL J2108.6$-$8619  &   316.99  &   -86.31   &   0.03   &   91.0 &   $1.65\pm0.12$   &   $10.3\pm4.0$ & 0.87  \\
3FGL J2312.9$-$6923  &   348.40  &   -69.39   &   0.04   &   35.3 &   $1.72\pm0.17$   &   $5.5\pm2.9$ &  0.86 \\
3FGL J2316.8$-$5209  &   349.28  &   -52.19   &   0.06   &   37.3 &   $1.89\pm0.16$   &   $10.8\pm5.0$ & 0.85  \\
3FGL J2347.9+5436  &   356.97  &   54.58   &   0.02   &   163.0 &   $1.79\pm0.08$   &   $24.2\pm6.4$ & 0.82  \\
\hline
\end{tabular}
\end{small}
\label{<Table 1>}
\end{table}

\begin{table}
\clearpage
\begin{small}
\caption{Same as Table~\ref{<Table 1>} for UCS$_{agn}$ HSP candidates.}
\begin{tabular}{lcccccccc}
\hline
\hline
Source name  & RA   &  DEC  &  $pos_{68}\%$	& TS  & Ph.Index   &  Flux &  $L_{HSP}$\\
\hline
            & [deg] & [deg]  & [deg]  &  &   & [$10^{-9}$ ph/cm$^2$/s] & \\
\hline
\hline
3FGL J0020.9+0323  &   5.26  &   3.36   &   0.04   &   60.7 &   $2.01\pm0.14$   &   $23.3\pm8.4$ &  0.88 \\
3FGL J0049.0+4224  &   12.26  &   42.38   &   0.04   &   37.0 &   $1.81\pm0.16$   &   $8.1\pm4.4$ &  0.82 \\
3FGL J0234.2$-$0629  &   38.59  &   -6.47   &   0.03   &   90.7 &   $1.83\pm0.11$   &   $15.6\pm4.8$ &  0.84 \\
3FGL J0251.1$-$1829  &   42.79  &   -18.50   &   0.02   &   104.3 &   $1.59\pm0.10$   &   $7.0\pm2.2$ & 0.88  \\
3FGL J0312.7$-$2222  &   48.15  &   -22.36   &   0.02   &   177.1 &   $1.84\pm0.08$   &   $22.3\pm5.2$ & 0.87  \\
3FGL J0506.9+0321  &   76.71  &   3.38   &   0.03   &   77.1 &   $1.81\pm0.12$   &   $15.3\pm6.0$ & 0.89  \\
3FGL J0524.5$-$6937  &   81.16  &   -69.61   &   0.03   &   94.1 &   $2.05\pm0.15$   &   $49.4\pm21.2$ & 0.86  \\
3FGL J0527.3+6647  &   81.86  &   66.80   &   0.03   &   51.9 &   $1.91\pm0.15$   &   $13.0\pm6.0$ & 0.83  \\
3FGL J0731.8$-$3010  &   112.96  &   -30.13   &   0.04   &   37.1 &   $1.96\pm0.17$   &   $22.4\pm11.9$ &  0.84 \\
3FGL J0813.5$-$0356  &   123.45  &   -3.95   &   0.04   &   57.0 &   $1.71\pm0.12$   &   $9.1\pm3.5$ & 0.88  \\
3FGL J0928.3$-$5255  &   142.09  &   -52.94   &   0.02   &   98.7 &   $2.09\pm0.09$   &   $88.0\pm25.5$ & 0.80  \\
3FGL J0952.8+0711  &   148.22  &   7.23   &   0.04   &   51.0 &   $1.92\pm0.15$   &   $14.0\pm6.0$ & 0.84  \\
3FGL J1057.6$-$4051  &   164.43  &   -40.87   &   0.03   &   40.2 &   $1.72\pm0.15$   &   $6.6\pm3.4$ & 0.82  \\
3FGL J1155.3$-$1112  &   178.82  &   -11.19   &   0.03   &   52.5 &   $2.03\pm0.15$   &   $21.2\pm8.9$ & 0.89  \\
3FGL J1222.7+7952  &   185.92  &   79.90   &   0.04   &   43.8 &   $2.13\pm0.22$   &   $21.1\pm12.1$ &  0.86 \\
3FGL J1225.4$-$3448  &   186.35  &   -34.75   &   0.05   &   22.3 &   $1.74\pm0.19$   &   $5.1\pm3.4$ &  0.86 \\
3FGL J1234.7$-$0437  &   188.71  &   -4.56   &   0.04   &   51.5 &   $2.01\pm0.14$   &   $23.5\pm9.7$ &  0.87 \\
3FGL J1513.3$-$3719  &   228.35  &   -37.39   &   0.03   &   54.7 &   $1.94\pm0.13$   &   $19.8\pm8.4$ &  0.87 \\
3FGL J1525.8$-$0834  &   231.53  &   -8.53   &   0.03   &   59.5 &   $1.92\pm0.12$   &   $20.0\pm7.3$ & 0.89  \\
3FGL J1528.1$-$2904  &   232.12  &   -29.11   &   0.06   &   26.3 &   $1.80\pm0.18$   &   $8.7\pm5.4$ & 0.83  \\
3FGL J1545.0$-$6641  &   236.21  &   -66.71   &   0.02   &   150.1 &   $1.60\pm0.10$   &   $11.2\pm3.9$ & 0.84  \\
{\bf 3FGL J1549.9$-$3044}  &   237.46  &   -30.75   &   0.02   &   64.3 &   $1.61\pm0.12$   &   $6.2\pm2.6$ & 0.91  \\
3FGL J1619.1+7538  &   244.78  &   75.61   &   0.02   &   107.1 &   $1.87\pm0.10$   &   $15.5\pm4.7$ & 0.88  \\
3FGL J1922.2+2313  &   290.57  &   23.25   &   0.03   &   80.8 &   $2.22\pm0.14$   &   $93.1\pm36.7$ & 0.87  \\
3FGL J2015.3$-$1431  &   303.81  &   -14.55   &   0.06   &   17.4 &   $1.81\pm0.21$   &   $5.6\pm4.2$ & 0.85  \\
3FGL J2043.6+0001  &   310.94  &   0.00   &   0.04   &   48.5 &   $2.02\pm0.15$   &   $21.5\pm8.1$ & 0.87  \\
3FGL J2053.9+2922  &   313.45  &   29.37   &   0.02   &   359.6 &   $1.77\pm0.06$   &   $46.0\pm8.4$ &  0.85 \\
{\bf 3FGL J2142.6$-$2029}  &   325.66  &   -20.50   &   0.04   &   36.1 &   $1.69\pm0.17$   &   $5.0\pm2.9$ & 0.91  \\
{\bf 3FGL J2145.5+1007}  &   326.38  &   10.13   &   0.03   &   34.1 &   $1.70\pm0.20$   &   $4.8\pm3.4$ & 0.90  \\
3FGL J2224.4+0351  &   336.12  &   3.89   &   0.05   &   29.5 &   $1.94\pm0.18$   &   $13.2\pm6.4$ & 0.89  \\
{\bf 3FGL J2300.0+4053}  &   345.06  &   40.88   &   0.03   &   52.5 &   $1.64\pm0.14$   &   $6.2\pm3.0$ &  0.90 \\
3FGL J2309.0+5428  &   347.20  &   54.41   &   0.03   &   77.1 &   $1.75\pm0.10$   &   $16.2\pm5.7$ &  0.85 \\
{\bf 3FGL J2321.6$-$1619}  &   350.40  &   -16.32   &   0.02   &   174.5 &   $1.73\pm0.08$   &   $17.3\pm4.1$ &  0.91 \\
\hline
\end{tabular}
\end{small}
\label{<Table 2>}
\end{table}
\begin{table*}
\clearpage
\begin{center}
\caption{\label{tab:zmax} Estimates of the maximum redshift values of the VHC sources for which
the sources are still detectable at $5\,\sigma$ in a 50\,hours of current IACT and 50 (5)\,hours of CTA observations,
respectively.
If no value is given, the source will not be significantly detected within the assumed observation time.
Sub- and superscript numbers give the change in the redshift values if the \emph{Fermi}-LAT spectrum
is extrapolated within 1$\,\sigma$ uncertainties.
See main text for further details and caveats.}
\begin{tabular}{l|ccc}
\hline
\hline
Source name & $\mathrm{T}_\mathrm{obs} = 50\,\mathrm{hours}$, current IACTs & $\mathrm{T}_\mathrm{obs} = 5\,\mathrm{hours}$, CTA & $\mathrm{T}_\mathrm{obs} = 50\,\mathrm{hours}$, CTA\\
\hline
\multicolumn{4}{c}{BCU HC$_{TeV} $ candidates}\\
\hline
3FGL J0047.9+5447 &  $0.18^{>0.50}_{\text{---}}$ & $0.15^{>0.50}_{\text{---}}$ & $> 0.50_{-0.38}$ \\
3FGL J1155.4$-$3417 &  $0.26^{>0.50}_{-0.17}$ & $0.35^{>0.50}_{-0.24}$ & $> 0.50$ \\
3FGL J1434.6+6640 &  $0.15^{+0.27}_{-0.13}$ & $0.14^{+0.33}_{\text{---}}$ & $> 0.50_{-0.31}$ \\
3FGL J0921.0$-$2258 &  $\text{---}$ & $\text{---}$ & $\text{---}$ \\
3FGL J0648.1+1606 &  $0.03^{+0.19}_{\text{---}}$ & $0.01^{+0.19}_{\text{---}}$ & $0.29^{>0.50}_{-0.26}$ \\
3FGL J1711.6+8846 &  $0.04^{+0.22}_{\text{---}}$ & $0.02^{+0.24}_{\text{---}}$ & $0.32^{>0.50}_{-0.29}$ \\
3FGL J1714.1$-$2029 &  $0.38^{>0.50}_{-0.28}$ & $0.47^{>0.50}_{-0.35}$ & $> 0.50$ \\
3FGL J1910.8+2855 &  $0.18^{+0.25}_{-0.15}$ & $0.16^{+0.28}_{-0.15}$ & $> 0.50_{-0.17}$ \\
3FGL J0153.4+7114 &  $\text{---}$ & $\text{---}$ & $0.43^{>0.50}_{\text{---}}$ \\
3FGL J0506.9$-$5435 &  $0.05^{>0.50}_{\text{---}}$ & $0.25^{>0.50}_{\text{---}}$ & $> 0.50_{\text{---}}$ \\
3FGL J1944.1$-$4523 &  $0.27^{>0.50}_{-0.21}$ & $0.34^{>0.50}_{-0.26}$ & $> 0.50_{-0.07}$ \\

\hline
\multicolumn{4}{c}{UCS$_{agn}$ HC$_{TeV} $ candidates}\\
\hline
3FGL J1549.9$-$3044 &  $0.13^{+0.27}_{\text{---}}$ & $0.15^{>0.50}_{-0.13}$ & $> 0.50_{-0.29}$ \\
3FGL J2142.6$-$2029 &  $0.05^{+0.28}_{\text{---}}$ & $0.07^{+0.35}_{\text{---}}$ & $0.41^{>0.50}_{-0.35}$ \\
3FGL J2321.6$-$1619 &  $0.23^{+0.20}_{-0.15}$ & $0.30^{>0.50}_{-0.20}$ & $> 0.50$ \\
3FGL J2145.5+1007 &  $0.04^{+0.33}_{\text{---}}$ & $0.02^{+0.35}_{\text{---}}$ & $0.26^{>0.50}_{\text{---}}$ \\
3FGL J2300.0+4053 &  $0.15^{+0.32}_{\text{---}}$ & $0.13^{+0.35}_{\text{---}}$ & $> 0.50_{-0.37}$ \\
\hline
\end{tabular}
\end{center}
\end{table*}