\section{Introduction}
\label{sec:intro}
Lossy compression is necessary to support the current infrastructure of modern media applications. However it introduces a wide range of spatial and temporal artifacts that negatively impact the end user quality of experience~\cite{9035388}. Statistical techniques for artifact removal as a post-processing step has been a long studied research topic~\cite{SHEN19982}.  More recently, Deep Neural Networks (DNNs) have been used in different areas of the compression pipeline, from end-to-end learnt compression~\cite{toderici2017full}, to purely post-processing~\cite{zhang2021video} and mixed post/pre-processing~\cite{guleryuz2021sandwiched}. The use of DNNs for compressed video enhancement has high potential for its ability to exploit perceptually salient metrics in the training process through the use of generative adversarial networks (GANs), something which is difficult to exploit with traditional statistical methods.

The use of transfer learning for crafting Deep Feature Quality Metrics (DFQMs) has been an emerging area of interest in the Deep Learning Community. Image DFQMs have been shown have a high correlation with human based subjective scores~\cite{ding2020image}. The underlying features used to generate these scores are the internal representations of high performing pre-trained image classification networks e.g. VGG~\cite{simonyan2014very}. LPIPS~\cite{zhang2018unreasonable} has been shown to outperform traditional quality metrics such as $\ell_{2}$, SSIM and FSIM. FloLPIPS~\cite{danier2022flolpips}, an extension of LPIPS that emphasizes areas of motion in a patch, has outperformed other traditional metrics including VMAF and MS-SIM.

\begin{figure}
    \centering
    \includegraphics[width=\columnwidth]{imgs/OverAll.png}

    \caption{\em Proposed W-GAN setup for compressed video enhancement. The generator take in motion information, the current patch and the patch from the nearest I-frame. The critic network includes selected extracted outputs EfficientNetB3 to assist with restoration.}
    \label{fig:Overall}
    \vspace{-1.8em}
\end{figure}

While these DFQMs are used in wide range of applications for assessing algorithmic performance and also in the loss-functions for training networks, the features that are used to generate these new metrics lack explainability. This is a characteristic hallmark of deep neural networks, however there are some approaches that can be used to help segmenting deep features and their importance in the current application.
\vspace{-1em}
\subsubsection*{Contributions of this paper}
This paper outlines the use of perceptually relevant deep features in the creation of a neural critic for a Wasserstein-GAN (W-GAN) with gradient penalty~\cite{gulrajani2017improved}. This neural critic is then used for training a compressed video enhancement network using data that is compressed with the H.265 codec. The following contributions are made:
\begin{itemize}
    \item A method for evaluating deep features of a pre-trained image classification network to isolate feature blocks that have perceptually relevant features.
    \item The design of a critic network that uses these selected feature blocks as a primary input and an associated compressed video enhancement network that is trained adversarially with the critic.
\end{itemize}

\vspace{-1em}
\subsection{Related Work}
Perceptually relevant DNN approaches have been recently used within the compression pipeline. Chen et al.~\cite{https://doi.org/10.48550/arxiv.2007.02711} investigates the use of loss functions that include perceptual metrics (VMAF, SSIM or MS-SIM) alongside MSE for training an end-to-end DNN image codec. Their results show more than 20\% coding gain with VMAF in the loss function when compared to a pure MSE loss. There has been sufficient work that show DFQMs~\cite{ding2020image, zhang2018unreasonable, danier2022flolpips} are more closely related to human subjective scores. This suggests that using the features that these metrics are dependent on has the potential to yield better results beyond using hand crafted metrics like VMAF. In our proposed model that focuses on post-processing for videos, we utilize extracted prior deep features from EfficientNetB3 directly in the neural critic of our W-GAN for training. We do not optimize for any other traditional metrics or DFQMs directly.

Mohammadi et al.~\cite{mohammadi2022perceptual} presents an investigation of the perceptual impact of changing the optimization metric of the loss function for a DNN image codec. Results indicate that training with the DFQM DISTS~\cite{ding2020image} and MS-SSIM yields highest perceptual gain. Images that were compressed with DISTS were selected with the highest winning frequency of 21\% when compared to a total of 9 different optimization schemes. In~\cite{mohammadi2022perceptual}, the entire set of deep features is used for the DFQM. In our work we instead identify the subset of features which are most relevant and deploy them in a more efficient adversarial training approach.

Ma et al.~\cite{https://doi.org/10.48550/arxiv.2011.09190} presents a post-processing GAN based approach that uses a mixture of multiple objective metrics. Their loss function incorporates the distance between complete deep features of a degraded-reference pair directly. In our approach, we do not use a difference between deep feature pairs. Instead, we improve performance by only using a subset of features that are perceptually relevant. These features are given to our critic network. This allows the critic to learn distinguishing features between deep representations of clean and degraded patches.

\vspace{-0.9em}
\section{Deep Feature Selection}

\begin{figure}
    \centering
    \includegraphics[width=0.9\columnwidth]{imgs/SelectedLayersBar.png}
    \caption{\em Averaged weights per layer in RDM analysis. Based on this measure of importance, layers $0$ and $3$ were selected as the top two most significant internal representations. This follows the notion that layers closer to the input of classification networks encode higher perceptual information.}
    \label{fig:barplot}
    \vspace{-1.8em}
\end{figure}

DFQMs such as LPIPS, FloLPIPS and DISTS are all dependent on the internal representations of pre-trained classification networks. The deep representation of an image is generated by extracting the output of intermediate layers in the network. Features are usually collected at the final output of the convolution filter bank operation of that layer.

These features can be represented by a set of rank 3 arrays whose height and width depends on the size of the input image. The number of channels in the feature array is dependent on the number of convolution filters used at that layer and can be represented by the layer depth. For an arbitrary image input $\tilde{I}$, the set of generated features $F$ can be represented by $F_{z}(\tilde{I}) \in F(\tilde{I})$, where $z$ represents the depth at which a feature is accessed.

Historic DFQMs ($Q$)~\cite{ding2020image, zhang2018unreasonable, danier2022flolpips} are based on comparing the internal representations of a reference ($x$) and degraded ($y$) patch, $Q(F(x), F(y))$. A key point is that the $F_z$ are extracted from {\em frozen} layers of a network trained for classification. Hence not all of the features will have significant relevance for perceptual criteria. Therefore,
in our proposed method, our critic utilizes only a subset of perceptually relevant features $F^{\alpha}$ from a single input $\hat{x}$: $F^{\alpha}(\hat{x}) \subset F(\hat{x})$. Discovering which subset is perceptually relevant is one of our main contributions. The intuition is that there is a set of compression invariant features that can be used by the W-GAN in restoration of a given patch. In addition, by selecting this subset we reduce the dimensionality of Q() and thus improve computational performance of the critic network by allowing to focus on a subset of key features.

A dataset of 99 reference $512\times512$ patches was used to generate deep feature representations from the classification network EfficientNetB3~\cite{tan2019efficientnet} pre-trained on ImageNet. This network was chosen for its relatively small number of parameters and increased performance over other classificaton networks. The 99 patches were compressed using the H.265 codecs at 6 linearly spaced constant rate factor (CRF) points in the operating range. This resulted in 99 reference and 594 degraded full feature representations.

Each $F$ for a single patch contains rank-3 arrays that total $\approx$ 3 million elements. Traditional statistical techniques for dimensionality reduction proves intractable for these high element multi-dimensional feature arrays. To overcome this limitation, we take inspiration from neuroscience~\cite{KRIEGESKORTE2013401, SupervisedMay}. In~\cite{SupervisedMay}, inferior temporal (IT) response to images presented to humans were compared to responses from layers in a pre-trained classification network when presented with those same images. Their goal was to discover which CNN layers were able to best model the human response. To do that, the authors exploited a dissimilarity measure \cite{KRIEGESKORTE2013401} for the CNN layer outputs which was then used to predict the same dissimilarity response from the human stimulus. In our case we are interested in discovering the specific layers in our representation (EfficientNetB3) which best enable the critic network to discriminate between degraded and clean video. In a sense the stimulus here is the input images and the response is the output of our CNN layers. Hence we follow the same strategy as the neuroscience case.

We first define Representational Dissimilarity Matrices (RDMs) for responses $F_z$ from layer $z$ and ${\bf F}$ from all layers, as follows.
\begin{align}
    R^{Y}(i,j) & = d\bigl({\bf F}(y^{i}), {\bf F}(y^{j})\bigr) \\
    R^{{}_{c}X_{z}}(i,j) & = d\bigl(F_{z}(x^{i}_{c}), F_{z}(x^{j}_{c})\bigr)
\end{align}
where $d(\dot, \cdot)$ is the euclidean distance, $y_k$ is the $k$th patch in our clean dataset and $x^k_c$ is the $k$ th patch in the datset degraded through compression with CRF $c$. Hence $F_z(x^i_c)$ is the response of the $z$th layer to the degraded patch $i$. $R^{Y}(i,j)$ is therefore a measurement of the auto-dissimilarity between all the clean patches in the dataset, while $R^{{}_{c}X_{z}}(i,j)$ is the auto-dissimilarity across the degraded patches.

By building a linear predictor for $R^{Y}(i,j)$ based on linear combination of elements across layers in $R^{{}_{c}X_{z}}(i,j)$, we can select layers based on significance for prediction. The linear predictor is defined through weights ${}_{c}\beta_{z}$ applied to each layer/degradation element in $R^{{}_{c}X_{z}}$  as follows.
\begin{equation}
\hat{R}^{Y}({}_{c}\beta_{z}) =
     \sum_{\forall z}\sum_{\forall c}{}_{c}\beta_{z}R^{{}_{c}X_{z}},
\end{equation}
We estimate ${}_{c}\beta_{z}$ by minimising the cosine distance between $\hat{R}^{Y}({}_{c}\beta_{z})$ and $R^{Y}$. We apply the additional constraint of non-negative weights with the BFGS algorithm\cite{9781118723203}.

By averaging the estimated weights across degradation level $c$ (6 in our case) we estimate the relative importance of each layer. This is shown in Figure~\ref{fig:barplot}. Layers 0 and 3 of EfficientNetB3 are therefore chosen for our critic network as they clearly encode most perceptual information.

\begin{figure}[]
\centering
\subfloat[\em Critic. An input patch is processed by EfficientNetB3 and only layers 0 and 3 (purple) are used in producing output scores. Computation stops at layer 3, other layers are included for clarity.]{%
  \includegraphics[clip,width=\columnwidth]{imgs/CriticNetwork.png}%
}

\subfloat[\em{
Generator. Motion magnitude and direction $m_t$ of the current frame is concatenated with the current RGB patch $f_t$. The corresponding patch from the nearest intra-coded frame is fed into a CNN block (blue). The resulting features are concatenated with the center of the UNet (red). These features are then fed into an attention module (CBAM~\cite{woo2018cbam}). The restoration parameter $d_t$ guides the enhancement network and $\Delta_{t}$ is a weighting factor that scales the features extracted from the I-frame.
}]{%
  \includegraphics[clip,width=\columnwidth]{imgs/Generator.png}%
}

\caption{{\em Proposed Critic (a) and Generator (b) networks.}}
\label{fig:networks}
\end{figure}

\section{Network Architecture}
 \vspace{-1em}

\textbf{Critic}: Our critic network takes a patch as input and produces a feature block that represents the patch degradation level. This is adapted from the PatchGAN architecture~\cite{isola2017image}. The main contribution of our critic network is the introduction of the selected EfficientNetB3 layers as a primary input. This is shown in Figure~\ref{fig:networks}(a).

The extracted features are processed in parallel by CNN Blocks. The input frame is also processed directly to capture information that EfficientNetB3 may have missed. These are then concatenated and processed further to produce the feature block.
There is no activation function in the final layer. This follows the requirement of a W-GAN for the critic output to be unbounded.

To compare our proposed critic with optimal layers ($C_{opt}$), we also train with critic networks that:
\begin{enumerate}
    \item Uses the output from layers 5 and 6 of EfficientNetB3 ($C_{5,6}$).
    \item Does not use the EfficientNetB3 module, and hence the entire critic is fully trainable ($C_{t}$).
\end{enumerate}

\noindent\textbf{Generator}: Our generator (Figure~\ref{fig:networks}(b)) uses a UNet structure as the main
network. Motion vector information (magnitude and direction) are concatenated with the input RGB patch. This is then fed into a UNet with skip connections. Features from the nearest Intra-coded frame (I-frame) are extracted by a CNN block and is concatenated with the center of the UNet. It is assumed that the nearest I-frame has a high amount of texture related to the current scene. The features extracted from the I-frame are scaled by the distance parameter $\Delta_{t}$, where $0 \leq \Delta_{t} \leq 1$. This allows the I-frame features to have a greater weight if they are closer to the current patch. $\Delta_{t}$ is calculated as : $\Delta_{t} = e^{-0.02(\gamma)}$, where $\gamma$ is the absolute distance between the frame index of input patch and the nearest I-frame.

The user-defined restoration parameter $d_{t}$ specifies the strength of the restoration and is proportional to the amount of degradation present in the input patch. This value is fed into all stages of the UNet to guide the restoration.

\begin{table}[]
\caption{\em The impact of varying $\Lambda$ for $C_{opt}$ compared against other critics and VBM4D.}
\label{tab:Results}
\resizebox{\columnwidth}{!}{%
\begin{tabular}{@{}ccccccccc@{}}
\cmidrule(l){2-9}
                                 & \textbf{Lambda} & \textbf{PSNR$_Y$}     & \textbf{PSNR$_{CBCR}$}  & \textbf{SSIM}       & \textbf{VMAF}        & \textbf{LPIPS}      & \textbf{FID}        & \textbf{KID}           \\ \midrule
\multirow{6}{*}{\textbf{C$_{opt}$}} & 0.0             & 26.92                & 33.24                & 0.80                & 46.82                & 0.45                & 0.44                & 3.3e-4                 \\ \cmidrule(l){2-9}
                                 & 1.0             & 26.07                & 34.82                & 0.79                & 43.67                & 0.45                & 0.44                & 3.5e-04                \\ \cmidrule(l){2-9}
                                 & 10.0            & 27.98                & 35.06                & 0.79                & {\ul \textit{50.43}} & 0.44                & 0.45                & 3.6e-04                \\ \cmidrule(l){2-9}
                                 & 100.0           & 28.25                & 36.87                & 0.81                & 48.72                & 0.43                & 0.40                & 3.8e-04                \\ \cmidrule(l){2-9}
                                 & 1000.0          & 29.11                & 37.75                & 0.83                & 47.28                & {\ul \textit{0.41}} & {\ul \textit{0.36}} & {\ul \textit{3.2e-04}} \\ \cmidrule(l){2-9}
                                 & 10000.0         & 29.43                & 38.12                & 0.83                & 47.28                & 0.42                & 0.37                & 3.7e-04                \\ \midrule
\textbf{C$_{5,6}$}                & 1000.0          & 29.64                & 38.41                & {\ul \textit{0.83}} & 45.91                & 0.43                & 0.39                & 3.8e-04                \\ \midrule
\textbf{C$_t$}                  & 1000.0          & 29.06                & 37.65                & 0.82                & 49.50                & 0.44                & 0.40                & 3.5e-04                \\ \midrule
\textbf{VBM4D}                   & -               & {\ul \textit{30.86}} & {\ul \textit{39.21}} & 0.82                & 41.09                & 0.43                & 1.2                 & 1.2e-02                \\ \midrule
\textbf{Degraded}                & -               &                      &                      &                     &                      &                     &                     &                        \\ \bottomrule
\end{tabular}}
\end{table}

\begin{figure*}[]
    \centering
    \includegraphics[width=0.95\textwidth]{imgs/Results.png}
    \caption{\em Generated patches. VBM4D (b) has the highest PSNR but suffers from oversmoothening. This is a characteristic of optimzing with respect to MSE alone. Our proposed method $C_{opt}$ with $\Lambda$ = 10 yields a high detailed image but the color does not match the original material. $C_{5,6}$ has the highest PSNR of all critic approaches, but still has a loss of detail. $C_{t}$ has higher DFQM scores than $C_{5,6}$, however there is a consistent streaking effect overlaid on degraded parts of the image. $C_{opt}$ with $\Lambda$ = 1000 has the best performing DFQM metrics.}
    \label{fig:results}
\end{figure*}

\section{Training}

A dataset with 2649 degraded/reference pairs was created by degrading clips using the CRF parameter on an H.265 codec such that the resulting PSNR$_{Y}$ of the degraded video was $\leq 38$dB. More details about the dataset can be found in the supplementary material~\cite{suppMaterial}. All training patches were 512x512.

The degradation parameter $d_{t}$, which controls the strength of the restoration network, was set to be the LPIPS score between the degraded input patch and the reference patch. Our motion magnitude and direction ($m_{t}$) was calculated using the DeepFlow algorithm~\cite{6751282} and the Luma channel of the previous patch, current patch and next patch.

Our generator was initially trained using the squared $l_{2}$ norm loss for 25 epochs. The outputs from this initially trained generator was then used to train the critic network for 10 epochs. This type of training is adapted from the NoGAN approach~\cite{salmona2022deoldify}. The Adam optimizier with learning rate 1e-04 was used in both cases.
The generator and critic network were then trained adversarially for 50 epochs using Adam optimizers with an initial learning rate of 1e-04 that decayed to 1e-06. Our adversarial training routine follows~\cite{gulrajani2017improved}, with the loss function for both the critic and generator as follows.
\begin{align}
    \mathcal{L}_{critic} & = - \biggl[ C(y_t) - C(G(f_t)) \biggr] + \lambda\Gamma \\
    \mathcal{L}_{gen} & = -C(G(f_t)) + \Lambda\bigl(l_2(\hat{f}_t, y_t)\bigr)
\end{align}
Where $C$, $G$ are the critic and generator network respectively. $\lambda$, $\Gamma$ are the gradient penalty weight and gradient penalty function (as in~\cite{gulrajani2017improved}). $\Lambda$ is a weight that allows the generator to strike a compromise between MSE ($l_2(\cdot, \cdot)$), and perceptual criteria. This is important
because generative networks like these do have the tendency to synthesise or hallucinate content which {\em looks like} picture content but is not constrained by the input pixel texture.
The MSE weight $\Lambda$ therefore allows the output to be constrained to be related more closely to the input content.

\section{Results}

Table~\ref{tab:Results} shows results of performance of the various systems tested using a unique test set of 501 reference/degraded clips. In addition to the usual performance metrics of VMAF and PSNR, we report Frechet Inception Distance (FID) and Kernel Inception Distance (KID)~\cite{heusel2017gans}. These latter metrics express distance in terms of the output of features from an Inception network. They have recently increased in importance for perceptual evaluation of picture quality as used~\cite{mentzer2022neural}.

In order to select the optimal value of $\Lambda$, the table reports performance for $C_{opt}$ under the effect of  varying $\Lambda$ on a log scale. $\Lambda = 1000$
 is chosen as optimal given that it minimises KID and FID.
Our proposed critic $C_{opt}$ (using perceptual layers) is the best performer with up to 10\% (FID) and 15\% (KID) average improvement  compared to other critics ($C_{5,6}$, $C_{t}$).

$\Lambda$ is used to weight the effect of the MSE loss function when training the generator. With lower $\Lambda$ the output content has a high level of detail but the color of the generated patches does not match the source material. This effect can be seen in Figure~\ref{fig:results}(e) and is corroborated by the low PSNR$_{Y/CBCR}$ scores. As $\Lambda$ increases the resulting PSNR$_{Y/CBCR}$ also increases as MSE is inversely proportional to PSNR.
The other critic configurations $C_{5,6}$ and $C_{t}$ were also trained with $\Lambda$ = 1000. Patches generated from training with $C_{5,6}$ contain less detail than the other critic configurations, however PSNR$_{Y/CBCR}$ was the best performing among critics. This effect can be seen in Figure~\ref{fig:results}(f).  Training with $C_{t}$ yielded DFQMs that were better performing than $C_{5,6}$, however there are visible streaking effects over areas of high distortion. This can be seen in Figure~\ref{fig:results}(g).
Considering the performance of more traditional statistical approaches, we find our optimal critic to outperform VBM4D in DFQMs. Patches generated from VBM4D suffered from the ``over-smoothing'' effect that is characteristic of MSE optimized algorithms. The visual impact of this is shown in Figure~\ref{fig:results}(b). This validates the need for perceptual criteria.

\section{Discussion \& Conclusion}
Drawing upon the principles of DFQMs, such as LPIPS, we developed a critic that employs layers from a pre-trained classification network to enhance critic proficiency. Although our model displays lower performance in conventional metrics such as VMAF, PSNR, and SSIM, it generates output that exhibit a close resemblance to the source content and this is reflected in performance increases in LPIPS, FID and KID (Shown in Table~\ref{tab:Results}) as well as visual evidence.

Our proposed critic $C_{opt}$, outperforms the other critics which indicates that the use of the selected layers from EfficientNetB3 exploits latent information that is invariant to degradation from compression. The superior perceptual performance achieved by our proposed model, as compared to VBM4D, provides evidence of the capactiy of the critic to focus on highly informative perceptual features beyond the conventional distance metrics.

A drawback of pure adversarial training for compression enhancement purposes is the introduction of hallucination artifacts, which can cause the appearance of plausible image features that are not present at all in the source content. To mitigate this, we use a weighted loss function using $\Lambda$ to combine both perceptual critic and MSE.
Lower values of $\Lambda$ imply more emphasis on perceptual loss and hence result in pictures with greater detail.

In future, we plan to conduct a subjective study to further validate our approach. We also plan to measure the importance of different layers of the critic and the generator. Data, full network architectures and models are available at~\cite{suppMaterial}.
\clearpage
\bibliographystyle{IEEEbib}


\section{\label{sec1:introduction}Introduction}
Materials based on III-V semiconductor compounds are attracting much attention in science and engineering due to their diverse applications in fields such as optoelectronics \cite{Soref1993Silicon,CardanoBook}. One of the main goals of basic and applied research is to tailor materials' optical properties to a specific application \cite{Hepp2022Room,Fuchs2018,Mokkapati2009IIIV,Dimroth2016Four,Mitchell2011Four,PhilippsS2018Chapter}. One of the most critical fundamental properties in this respect is the bandgap, both in terms of size and type (direct or indirect). For example, optical telecommunication applications require materials with direct bandgaps in the range of 0.80--0.95 eV \cite{Hepp2022Room,Fuchs2018,Mokkapati2009IIIV}, while solar cell applications require a range of 0.5--2.0 eV \cite{Dimroth2016Four, Mitchell2011Four, PhilippsS2018Chapter}. Composition engineering, i.e., changing the relative composition of group 13 and 15 elements in ternary III-V compounds, is one of the most important approaches to adjusting the bandgap  \cite{Beyer2015Metastable,Vurgaftman2001Band,Beyer2017Local,Ludewig2016Movpe,Liebich2011Laser,Supplie2018Metalorganic,Stringfellow2019Fundamental,Volz2009Movpe,Feifel2017Movpe,Kunert2004Movpe,Volz2004Specific,Veletas2019Bismuth,Wegele2016Interface,Hepp2019Movpe}. Systematic application of strain such as mechanical strain (e.g., external pressure \cite{Van1990Pressure,Potter1956Indirect,Alekseev2020Effect,Katiyar2020Breaking}, mechanical bending of nanowire \cite{Lim2021Strain,Signorello2014Inducing,Signorello2013Tuning}) or strain due to lattice mismatch (e.g., core-shell mismatch in nanowires \cite{Balaghi2019Widely,Gronqvist2009Strain,Hetzl2016,Montazeri2010,Skold2005}) on a system are alternative strategies to tailor the bandgap. Combining composition and strain engineering, the bandgap can be tuned over a wide range of values, and direct or indirect semiconductors can be designed. In thin-layer heteroepitaxy, choosing the substrate-layer combination with minimum lattice mismatch is often desirable to minimize the strain effect from the substrate. However, in practice, perfect lattice matching is rarely possible. In such cases, not only the composition but the effect of inherent strain from the substrate also substantially affects the active layer's bandgap \cite{Beyer2015Metastable,Vurgaftman2001Band,Beyer2017Local,Ludewig2016Movpe,Liebich2011Laser,Supplie2018Metalorganic,Stringfellow2019Fundamental,Volz2009Movpe,Feifel2017Movpe,Kunert2004Movpe,Volz2004Specific,Veletas2019Bismuth,Wegele2016Interface,Hepp2019Movpe,Bir1974Symmetry,Sun2007Physics,Tao2020,Tsutsui2019,Fang2011}. Therefore, one requires a complete knowledge of the material-specific dependence of the bandgap on composition and strain to optimize the proper choice of materials. However, exploring the vast chemical space of all possible combinations of III- and V-elements with variation in composition and strain is experimentally not feasible. Additionally, growing a new material is often challenging because of thermodynamic or kinetic limitations, such as phase separation or surface roughening, in addition to the demanding task of optimizing the growth conditions \cite{Beyer2015Metastable,Liebich2011Laser,Volz2009Movpe,Volz2004Specific,Wegele2016Interface,Hepp2019Movpe}. This makes an experimental screening approach of vast compound and strain spaces unrealistic. We thus aim to develop a reliable and predictive theoretical approach.

Two major theoretical approaches that have been developed to analyze strain effects on the bandgap of III-V materials are (semi-)empirical methods and ab initio approaches. Although (semi-)empirical methods such as k$\cdot$p theory \cite{Sun2007Physics,Bannow2017} and tight-binding methods \cite{Sun2007Physics,Anderson1991Optimized,Tan2016Transferable,Jancu1998} are computationally efficient, they rely on empirical parameters which require system-specific experimental input data. This strongly limits the predictive ability of these methods for new materials. Additionally, in case of a large mismatch in atomic sizes, the ternary material shows local strain effects, severely affecting the bandgap \cite{Rosenow2018Ab}. These local strain effects, however, can not be included in empirical approaches and, hence, are neglected. Then again, ab initio approaches such as density functional theory (DFT) \cite{Bannow2017,Mondal2022,Kim2010Towards,Tran2009Accurate,Jiang2013Band,Rehman2016Electronic,Koller2012Improving} allow for the calculation of electronic properties from electronic structure theory and are thus predictive if accurate density functionals are used. The local relaxation of the atomic positions also allows to properly include and investigate the effect of the local strain on the electronic properties in these approaches. Additionally, recent advancements in the modeling strategies of alloy systems using quasirandom supercells \cite{Wang1998,Medeiros2015,Popescu2010,Rubel2014,Popescu2012} allow for electronic properties calculations in the ab initio approaches, even for diluted and disordered materials. An alternative to ab initio approaches is the use of GW-based methods, which are nevertheless too computationally demanding for screening approaches as intended here \cite{Rosenow2018Ab,Mondal2022,Kim2010Towards,Hinuma2014Band}.

In a previous study, we established a computational protocol for predictive modeling based on DFT for binary III-V compounds over a wide range of strain values \cite{Mondal2022}. In this study, we are extending this approach to ternary III-V compounds, which then allow the combination of strain and composition to fully explore bandgap design. For ternary systems, only the effects of composition variations on the bandgap in unstrained materials have been studied \cite{Bannow2017,Anderson1991Optimized,Rosenow2018Ab,Wang1998,Medeiros2015,Popescu2010,Rubel2014,Popescu2012}. For strained materials, a suitable theoretical framework is still lacking. We present here a predictive first-principle protocol for a complete mapping of the mutual correlation of composition, strain, and bandgap in ternary III-V semiconductor systems. The goal is to provide guidelines for assessing and identifying the most promising target materials for experimental investigations in the future.

We start by describing the computational methods in \cref{sec2:computationaldetails}. Next, we describe the protocol for determining the nature of the bandgap from supercell calculations using GaAsP as an example in \cref{sec3:protocolfordeterminingbandgapnature}. We further present the composition-strain-bandgap correlation results for different ternary III-V semiconductors in \cref{sec4:results}. We start with GaAsP, an experimentally well-studied and promising candidate for LEDs, detectors, and Si-based multi-junction solar cells \cite{Craford1973,Henning1983,Tanaka1994,Sato2002,Geisz2002,Lang2013,Hayashi1994,Grassman2016}. The results for the GaAsN compound, a promising laser-active material \cite{Rosenow2018Ab,Weyers1992,Kunert2006,Zhao2004}, are presented next. To show the general applicability of our approach, we then show selected results for (i) GaPSb, a candidate for vertical cavity emitting surface laser \cite{Loualiche1998,Shimomura1996,Nakajima2000,Russell2016,Jou1988}, (ii) GaAsSb, a material for tandem solar cell application \cite{Cherng1984,Jen1998}, (iii) GaPBi, a promising material for near-infrared photonic device application on Si \cite{Christian2016,Christian2015}, and (iv) GaAsBi, another material discussed for near- and mid-infrared photonic device application \cite{Sweeney2011,Wang2013,Cooke2006}. We then discuss the comparison of our computations with experiments in \cref{sec5:discussion}, underlining the accuracy and predictive capability of our computational approach.

\section{\label{sec2:computationaldetails}Computational details}
The calculations were performed with DFT-based approaches as implemented in the Vienna Ab initio Simulation Package (VASP 5.4.4) \cite{Kresse1993Ab,Kresse1994Ab,Kresse1996Efficient,Kresse1996Efficiency}, using plane wave basis sets in conjunction with the projector-augmented wave (PAW) approach \cite{Kresse1999From,Blochl1994Projector}. The ternary materials were generated using the special quasirandom structures (SQS) approach \cite{Zunger1990Special} with a supercell of size $6 \times 6 \times 6$. The SQS cells were generated using the Alloy Theoretic Automated Toolkit (ATAT) \cite{VandeWalle2002,VandeWalle2009,VanDeWalle2013}. Geometry optimization of the supercells was performed using the PBE functional \cite{Perdew1996Generalized}, including the dispersion-correction method DFT-D3 with an improved damping function \cite{Grimme2010Consistent,Grimme2011Effect}. The basis set energy cut-off was set to 450 eV. The electronic energy convergence criteria of $10^{-6}$ eV and the force convergence of $10^{-2}$ eV\AA $^{-1}$ were used. The reciprocal space was sampled at the $\Gamma$-point only, given the large supercells used \cite{Monkhorst1976}. The meta-GGA functional TB09 \cite{Tran2009Accurate} was used to calculate the electronic properties (bandgaps and band structures). The effects of spin-orbit coupling were considered in the TB09 calculations. For the meta-GGA calculations, the energy cut-off of the basis set and the convergence criterion for the electronic energy were lowered to 350 eV and $10^{-4}$ eV, respectively, to minimize the computational costs. Structure optimizations were carried out by consecutive volume and position optimization until convergence was reached. This setup was previously used to generate bandgaps in excellent agreement with experimental data \cite{Rosenow2018Ab}.

We considered the materials investigated here all feature the zincblende type structure only. Moreover, [100] crystal direction is the most common choice of substrate orientation and growth direction in epitaxy. Therefore, we modeled the strain application along [100] directions only. The isotropic strain was modeled by increasing (decreasing) all the lattice parameters of the unstrained structure by the same amount. In this case, only the atomic positions of the strained structure were optimized, keeping the volume fixed. For biaxial strain, the in-plane lattice parameters were kept fixed and relaxed in the out-of-plane direction. No structural phase transition is assumed under strain application. More details on the strain modeling can be found in Ref.~\cite{Mondal2022}. In the following, we indicate the tensile strain with a `$+$' sign and the compressive strain with a `$-$' sign.

The calculations were performed at discrete points in composition-strain space (\cref{fig:figS6}). The discrete DFT calculated bandgap values were then interpolated to create the final images in \cref{fig:fig2,fig:fig3,fig:fig4,fig:fig5,fig:fig6}. Noticeably, for the systems we addressed in this article, the variations of bandgap values with concentration and strain are mostly non-monotonic (see \cref{fig:figS6}). This resulted in non-smooth interpolation in \cref{fig:fig2,fig:fig3,fig:fig4,fig:fig5,fig:fig6}. It is to be stressed that the origin of the non-smooth patterns is neither an interpolation artifact nor a deficiency of our DFT protocol. This solely originated because of the non-monotonic variation of the bandgap values (in the composition-strain space) of the SQS cells that we used to calculate bandgaps. A choice of positive smoothening during interpolation (e.g., bivariate B-spline, gaussian filtering) could mitigate the problem but significantly reduced the interpolated bandgap values and, thus, avoided. Further detail of the interpolation procedures can be found in \cref{secSVI:interpolationdetails}. Moreover, the nature of bandgaps can solely be deduced from the enclosed regions by the direct-indirect transition lines and thus require no interpolation.

\section{\label{sec3:protocolfordeterminingbandgapnature}Protocol for determining bandgap nature}
\begin{figure}[hbtp]
\centering
  \subfloat[]{\label{fig:fig1a}\includegraphics[width=3.4in]{figure1a.png}}
  \hspace{.2in}
 \subfloat[]{\label{fig:fig1b}\includegraphics[width=3.4in]{figure1b.png}}
  \caption{\label{fig:fig1}Variation of the bandgap under isotropic compressive strain for GaAs$_{0.963}$P$_{0.037}$. The $\Gamma$-, L-, and X-BSW of the folded supercell conduction band are given in parentheses. The vertical lines in (a) separate regions where the CBM changes character. In (b), the strain resolution is increased to determine the point of direct-indirect transition more accurately, indicated by the red circle (where the highest BSW changes from $\Gamma$ to L).}
\end{figure}
Supercell calculations, as required for modeling ternary semiconductors, lead to the folding of band structures \cite{Ku2010,Yang2018}. Although the size of the bandgap can be well extracted from the folded band structure, determining the bandgap's nature requires the primitive Bloch character of the bands to be known, which gets mixed up in the supercell eigenstates. With the band unfolding method, one projects these supercell eigenstates on the eigenstates of a suitable reference primitive cell. This requires the calculation of Bloch spectral weights (BSW), which measure the fraction of the primitive Bloch character in a supercell eigenstate. The result is an effective band structure (EBS) \cite{Wang1998,Medeiros2015,Popescu2010,Rubel2014,Popescu2012}. The spectral weights, $w_{n, {\bf K}}\left( {\bf k} \right)$, can be calculated from the plane wave coefficients as described in Ref.~\cite{Popescu2012}:
\begin{equation}
    w_{n, {\bf K}}\left( {\bf k}_j \right) = \sum _{\bf g} | C_{n, {\bf K}} \left( {\bf g} + {\bf G}_j\right) |^2
\end{equation}
Where $n$ is the band index, ${\bf g}$ and ${\bf G}_j$ are the reciprocal lattice vectors of the primitive and supercell, respectively. The index $j$ accounts for the series of primitive vectors, ${\bf k}_j = {\bf K} + {\bf G}_j$. The code `fold2Bloch' from Ref.~\cite{Rubel2014} was used to calculate the spectral weights.

In our previous study on binary III-V systems \cite{Mondal2022}, we have shown that the valence band maxima (VBM) always remain at the $\Gamma$-point, and only the conduction band minima (CBM) change their position in reciprocal space under strain. We have also shown that the CBM occurs only at the $\Gamma$-, L-, and (near) X-point in the band structure under strain. Therefore, it is sufficient to trace the conduction band (CB) at these points to determine the nature of the bandgap. As in our previous study on binary systems, we focus here on analyzing ternary III-V compounds with zincblende structures. For these structures in the ($6 \times 6 \times 6$) supercell dimensions chosen here, the $\Gamma$-, L- and X-point of the primitive band structure fold to the $\Gamma$-point in the supercell \cite{Popescu2010,Rubel2014,Popescu2012}. Therefore, it is sufficient to calculate the BSWs of solely the CB at the $\Gamma$-point in the supercell calculation to determine the nature of the bandgap. Consequently, we performed the supercell calculations by sampling the reciprocal space only at the $\Gamma$-point and unfolded the CB.

Figure~\ref{fig:fig1} shows the steps for determining the bandgap nature from supercell calculations more clearly. Figure~\ref{fig:fig1a} shows the bandgap variation for GaAsP with 3.7\% P concentration under isotropic compressive strain. The $\Gamma$-, L-, and X-BSWs of the folded supercell CB are given in parentheses. This shows 100\% $\Gamma$-BSW for the unstrained structure in line with the direct bandgap. With increasing strain, the $\Gamma$-BSW decreases (first number in brackets), and the L-BSW increases (second number in brackets). After a certain amount of strain, the L-character of the CB dominates. The bandgap becomes indirect in nature. Further compressing the system then leads to a transition of the CB character from L to X. In \cref{fig:fig1b}, we show calculations with increased resolution in strain to accurately determine the transition to L corresponding to the sought point of direct-indirect transition (DIT) at $-1.4$\% strain. We define the last strained structure with bandgap of direct nature before the transition to the indirect bandgap as the transition point (the red circle in \cref{fig:fig1b}). In the supporting information (\cref{fig:figS1}), we have given the EBSs of GaAs$_{0.963}$P$_{0.037}$ for different strain values. These confirm our analyses.

If the difference in BSW between different points in k-space is large, the nature of the bandgap can be unanimously determined. However, close to the transition points, in some cases, the differences are more subtle (\cref{fig:figS2}). We, therefore, set a cut-off criterion of 20\% BSW. If the $\Gamma$-BSW is larger than the cut-off criterion, then the direct transition has a finite probability even if the L- or X-BSW is larger than the $\Gamma$-BSW. In such cases, the bandgap is called `partially direct'. This defines a `region of uncertainty in the bandgap nature'. We chose the 20\% cut-off criterion because this produces results that agree best when compared to the experiments for several systems. For the GaPBi system, however, the 10\% BSW cut-off criterion produces the best agreement.

In some systems such as GaAsN, the band originating from the added nitrogen atoms, the so-called `defect N state' \cite{Rubel2014,Goodrich2019,Wei1996,Wu2004,Misiewicz2005}, is strongly dispersed under strain (see \cref{fig:figS3b,fig:figS4b}). Therefore, we set another cut-off criterion of 20\% BSW as a minimum limit for a (defect) eigenstate to be considered an eigenstate (\cref{fig:figS5}). If none of the CBs satisfy the cut-off criterion, we use the lowest CB for determining the bandgap nature.

\section{\label{sec4:results}Results}
In this section, we present the bandgaps calculated for different materials and determine their nature according to the above protocol. We mapped the bandgaps in terms of their size and nature for various strained ternary III-V compounds. We start with two important ternary III-V semiconductor materials, GaAsP and GaAsN. Then we show selected data for the material systems GaPSb, GaAsSb, GaPBi, and GaAsBi.
\subsection{\label{subsec41:gaasp}G\lowercase{a}A\lowercase{s}P}
\begin{figure}
    \centering
    \includegraphics[width=3.4in]{figure2.png}
    \caption{Isotropic strain for GaAsP. The variation of bandgap magnitudes (E$_{\text{g}}$) and type as a function of composition and strain. The dashed black horizontal line indicates unstrained GaAsP. The black circles are the calculated DIT points. The `DIRECT' and `INDIRECT' correspond to the enclosed region with the nature of bandgap direct and indirect, respectively. The hatched pattern region is the `uncertainty region' (see \cref{sec3:protocolfordeterminingbandgapnature}).}
    \label{fig:fig2}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=3.4in]{figure3.png}
    \caption{Biaxial strain for GaAsP. The variation of bandgap magnitudes (E$_{\text{g}}$) and type as a function of composition and strain. The dashed black horizontal line indicates unstrained GaAsP. The black circles are the calculated DIT points. The DIT points are fitted with 5th order polynomial. The `DIRECT' and `INDIRECT' correspond to the enclosed region with the nature of bandgap direct and indirect, respectively. The hatched pattern region is the `uncertainty region' (see \cref{sec3:protocolfordeterminingbandgapnature}). Solid black lines indicate the substrate lines under `epitaxial growth' model.}
    \label{fig:fig3}
\end{figure}
For the case of isotropic strain, \cref{fig:fig2} shows the bandgap as a function of composition ($x=$ 0--100\% in GaAs$_{1-x}$P$_{x}$) from 5\% tensile to 5\% compressive strain. The bandgap value varies between 0.32--2.42 eV in the strain regime investigated. For the same amount of P concentration, the bandgap primarily increases in moving from tensile to compressive. Furthermore, the figure shows that in going from compressive to tensile strain, the DIT occurs at a higher concentration of P atoms. The dashed horizontal line marks the data corresponding to the unstrained structures for different fractions of P. The intersection of this line with the DIT line shows at which percentage of phosphorous contribution the unstrained structure shows a DIT. This transition occurs at $x = 37$\%. Here, the bandgap shows a value of 1.96 eV. The terms `DIRECT' and `INDIRECT' in the figure correspond to the area where the bandgap is direct and indirect, respectively. Due to the similarity with commonly used phase diagrams, we call this representation a ``bandgap phase diagram". This and the following figures thus provide a 2D representation of the bandgap phase diagram for the ternary materials.

For the biaxial strain regime, \cref{fig:fig3} shows the bandgap phase diagram for GaAsP as a function of composition from 5\% tensile to 5\% compressive strain. The value of the bandgap varies in a range of 0.82--2.42 eV. For the same amount of P atoms, the bandgap reaches a maximum around the unstrained structure and gets smaller for tensile as well as compressive strain. This is different from the isotropic strain case. For unstrained GaP, the bandgap value is 2.36 eV. The nature of the bandgap also shows a different trend compared to \cref{fig:fig2}. The range of strain around the unstrained structure where a direct bandgap is found gets smaller for higher amounts of P. This is in line with GaAs (0\% P) being a direct and GaP (100\% P) being an indirect semiconductor. The largest amount of P concentration where a direct semiconductor is found is 39--40\% P in the unstrained structure. This is similar to the previous experimental result (45\% P) \cite{Vurgaftman2001Band,Capizzi1981}.

One of the most common approaches to experimentally realize biaxial strain in III-V semiconductors is epitaxial growth. As pointed out in Ref.~\cite{Mondal2022}, biaxial strain can be used to model epitaxial growth. We thus investigate the effect of different substrates in our bandgap phase diagram (\cref{fig:fig3}), where each solid line corresponds to one substrate: GaAs, GaP, InP, or Si. These solid lines indicate how much biaxial strain would develop in the GaAsP system as the respective value of \%P when grown on the respective substrates. The (substrate) strains are calculated
according to \cref{eqn:equation2}.
\begin{equation}
    \text{Substrate strain} (\%) = \frac{a_{sub}-a}{a}\times 100 \label{eqn:equation2}
\end{equation}
Where $a_{sub}$ is the equilibrium lattice parameters of the substrates, and $a$ is the lattice parameters of unstrained GaAsP systems at their respective P concentrations. E.g. for 100\% P, the strain on the GaP substrate is zero, while growing GaAs (0\% P) on GaP would result in 3.8\% in-plane compressive strain. This, of course, neglects defect formation and strain relaxations and assumes perfect epitaxial growth. Clearly, by choosing different substrates, the nature and the size of the bandgap can be tuned over a wide range. We refer to the next section for a comparison of our calculations to experimental data.

\subsection{\label{subsec42:gaasn}G\lowercase{a}A\lowercase{s}N}
\begin{figure}
    \centering
    \includegraphics[width=3.4in]{figure4.png}
    \caption{Isotropic strain for GaAsN (upto 12\% N). The variation of bandgap magnitudes (E$_{\text{g}}$) and type as a function of composition and strain. The dashed black horizontal line indicates unstrained GaAsN. The black circles are the calculated DIT points. Beyond 7\% N, the DIT is outside the investigated strain regime. 10 SQS cells are used for each configuration and strain point. The bandgaps plotted are the average bandgaps. The `DIRECT' and `INDIRECT' correspond to the enclosed region with the nature of bandgap direct and indirect, respectively.}
    \label{fig:fig4}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=3.4in]{figure5.png}
    \caption{Biaxial strain for GaAsN (upto 12\% N). The variation of bandgap magnitudes (E$_{\text{g}}$) and type as a function of composition and strain. The dashed black horizontal line indicates unstrained GaAsN. 10 SQS cells are used for each configuration and strain point. The bandgaps plotted are the average bandgaps. All the bandgaps are direct in nature. Solid black lines indicate the substrate lines under the `epitaxial growth' model.}
    \label{fig:fig5}
\end{figure}
As the next material, we investigate GaAsN. First, we show results for isotropic strain, which results in the bandgap phase diagram shown in \cref{fig:fig4}. The results are markedly different from GaAsP, and the data set is much more limited. In agreement with previous data \cite{Rosenow2018Ab}, we found that the size of the bandgap strongly depends on the distribution of N atoms in the supercell, even when using the SQS approach. We thus used 10 SQS cells for each data point in the figure and averaged the resulting bandgaps. This results in an error bar for the DIT points, which is rather large for medium amounts of nitrogen atoms due to the formation of small clusters and chains. Calculations were only possible for up to 12\% N atoms. For higher concentration and/or high compressive strain, our chosen supercell is not large enough to avoid the unphysical electronic interaction of N atoms with their images in the periodic boundary condition approach. This effect has already been discussed in Ref.~\cite{Rosenow2018Ab}. For the strain and composition regions where computation was possible, an indirect gap is only found for low values of \%N and rather large compressive strain values. The EBSs for selected \%N and strain values are shown in \cref{fig:figS3,fig:figS4}.

For biaxial strain in GaAsN, the data are shown in \cref{fig:fig5}. In contrast to GaAsP, in this case, the bandgap gets smaller with the increasing amount of nitrogen in the system, from 1.47 eV for the unstrained case of GaAs to 0.10 eV for the highly strained systems with a large number of N atoms. All bandgaps computed are direct. Epitaxial growth on GaAs is reasonably possible for moderate strain values and results in a variation of bandgap from 1.47--0.45 eV. For GaP and Si substrates, a large strain would be exerted on the system, and mostly lower bandgap values are found.

\subsection{\label{subsec43:gapsbgaassbgapbigaasbi}G\lowercase{a}PS\lowercase{b}, G\lowercase{a}A\lowercase{s}S\lowercase{b}, G\lowercase{a}PB\lowercase{i}, G\lowercase{a}A\lowercase{s}B\lowercase{i}}
\begin{figure*}
\centering
    \subfloat[GaPSb]{\label{fig:fig6a}\includegraphics[width=3.4in]{figure6a.png}}
    \hspace{0.2in}
    \subfloat[GaAsSb]{\label{fig:fig6b}\includegraphics[width=3.4in]{figure6b.png}} \\
    \subfloat[GaPBi]{\label{fig:fig6c}\includegraphics[width=3.4in]{figure6c.png}}
    \hspace{0.2in}
    \subfloat[GaAsBi]{\label{fig:fig6d}\includegraphics[width=3.4in]{figure6d.png}}
  \caption{\label{fig:fig6}Bandgap phase diagram for ternary III-V semiconductors GaEY (E= P, As; Y = Sb, Bi) under biaxial strain. The bandgap magnitudes (E$_{\text{g}}$) are shown in the color bar. The dashed black horizontal line indicates unstrained structures. The black circles are the calculated DIT points. The `DIRECT' and `INDIRECT' correspond to the enclosed region with the nature of bandgap being direct and indirect, respectively. The hatched pattern region is the ‘uncertainty region’ (see \cref{sec3:protocolfordeterminingbandgapnature}). Solid black lines indicate the substrate lines under the `epitaxial growth' model.}
\end{figure*}
The approach outlined here can be extended to other combinations of constituents in III-V semiconductor materials. Exemplarily, we present the bandgap phase diagrams for four other ternary compounds in \cref{fig:fig6}. Since epitaxial growth is the most interesting experimental realization method for these compounds, we only present the data for biaxial strain.

For all compounds investigated, we find a bandgap range of 0.0--2.36 eV, with the largest values found for the host materials GaP and GaAs in the unstrained case. The alloys with Sb could be investigated over the full range of 0--100\% Sb in GaPSb and GaAsSb (\cref{fig:fig6a,fig:fig6b}). We find a DIT for unstrained GaPSb at 30\% Sb concentration (\cref{fig:fig6a}). For this compound, the DIT is shown as a region, including the uncertainty in determining the nature of the bandgap, as outlined in \cref{sec3:protocolfordeterminingbandgapnature}. With the increase in the Sb fraction, the strain at which the DITs take place increases. For GaAsSb and GaAsBi, the bandgap is direct throughout the range investigated (\cref{fig:fig6b,fig:fig6d}), while it is indirect for GaPBi (\cref{fig:fig6c}). Notably, we only investigate the bismides up to a fraction of 15\% Bi. The reason is that, similar to GaAsN (\cref{fig:figS4}), for structures with large Bi content, the strongly dispersed bands decrease the reliability in the determination of bandgap nature.  Additionally, GaPBi and GaAsBi become metallic for higher Bi fractions. Our DFT protocol based on the TB09 functional produces less reliable bandgap values in metallic systems. Although we find no transition within 15\% Bi, it can not be excluded that the DIT appears at higher percentages of bismuth.

Again, we indicate the strain values associated with different typical substrates for epitaxial growth by solid black lines in the figures. The data show that deviating from the substrate-layer lattice-matching condition quickly leads to high strain, and defects are highly likely to occur during growth. Also, Si can be used as a substrate for GaPSb and GaPBi epitaxial growth if the Sb or Bi content is not too large. The epitaxial growth of the respective GaAs-based materials (GaAsP, GaAsN, GaAsSb, GaAsBi) will give rise to high strain on Si-substrate throughout the whole composition region. A noticeable change in the slope in substrate lines is found close to 10\% Bi and 40\% Sb concentration in GaAsBi and GaAsSb, respectively. Although we did not find any structural phase transition in those regions, the origin of the change in the slope is not clear to us yet.

From the above discussions, it becomes clear that bandgap phase diagrams can be a valuable aid in deciding which substrates are good choices for targeting a specific bandgap size and nature for a given ternary material. And vice versa, which material to grow for a specific application and a given substrate? We will discuss this further in the next section.

\section{\label{sec5:discussion}Discussion}
All data were derived from DFT computations to this point. In \cref{tab:table1}, we now compare our calculated bandgaps with experimental data from measurements on heteroepitaxial layer structures. The GaAsP/GaAs samples were grown by low-pressure hydride vapor phase epitaxy (LP-HVPE). Further details can be found in Ref.~\cite{Stromberg2020}. The remaining samples were grown by metalorganic vapor phase epitaxy (MOVPE). The details of the growth characteristics of the MOVPE samples can be found in the Refs.~\cite{Hepp2019Movpe,Jou1988,Nattermann2017,Ludewig2013,Ludewig2017,Lewis2012,Masnadi-Shirazi2014,Ludewig2014GrowthMOVPE}.  Experimentally, the layer thickness and bandgaps of the MOVPE samples were determined using X-ray diffraction and room-temperature photoluminescence (RT-PL), respectively. Except for GaPSb samples from Ref.~\cite{Jou1988}, in which cases, the PL were measured at 10K.

The comparison of the experimental bandgaps with our computed results shows impressive agreement. The deviation is determined with respect to the root-mean-square-deviation (RMSD) from all available experimental samples and is around 0.1 eV for most structures; for GaAsP, it is even more accurate. This confirms our previous findings on unstrained structures that the DFT protocol we developed gives excellent agreement to experimental bandgaps \cite{Beyer2017Local,Bannow2017,Rosenow2018Ab,Mondal2022}. In this study, we show that it is also applicable to compound semiconductors under strain. For samples with very small layer thickness, the matching of experiment and computation is less good. This can be observed for GaAsSb/GaAs samples with a deviation of ca. 0.4 eV. We attribute this to the 2D quantum confinement effect which is found for thin samples. This confinement effect leads to an increase in the bandgap with respect to thicker samples \cite{Cipriano2020}. This effect is not captured in our computational model as the calculations were performed for 3D periodic strained structures. The large deviation for GaAsN/GaP and GaAsN/Si samples can not be explained by this effect alone, though. An additional effect here is the strong dependency of the bandgap on the distribution of N atoms which has been found for GaAsN before \cite{Rosenow2018Ab}. The dependency is further enhanced under strain in those samples. No RT-PL was observed for GaPSb/Si and GaPBi/GaP samples. This is consistent with our findings that those materials show indirect bandgaps (\cref{fig:fig6a,fig:fig6c}). Experimental measurements of magnitude of the indirect bandgaps are not available yet.
\begin{table*}
\caption{\label{tab:table1}Comparison of the calculated bandgaps for investigated ternary III-V semiconductors under biaxial strain with the experiments. The experimental data are for the heteroepitaxial layer structures, and the bandgaps are determined from photoluminescence (PL) measurements. The `I' in the brackets indicate the indirect nature of the bandgap. The remaining bandgaps are direct. (RT-PL: room temperature PL)}
\begin{ruledtabular}
\begin{tabular}{cccccccc}
 & & & &\multicolumn{2}{c}{Bandgap (eV)} & & \\ \nmidrule{5-6}
 System & Substrate & $x(\%)$ & Layer thickness (nm) & Calculated & Experiment & Deviation (eV) & RMSD (eV) \\ \hline\Tstrut\Bstrut
 \multirow{3}{*}{GaAs$_{1-x}$P$_x$} & \multirow{3}{*}{GaAs} & 18.0 & 9000 & 1.65 & 1.66 & 0.01 & \multirow{3}{*}{0.01} \\
 & & 25.0 & 12250 & 1.72 & 1.72 & 0.00 &  \\
 \cite{Stromberg2020}& & 28.0 & 13000 & 1.75 & 1.76 & 0.01 &  \\ \hline\Tstrut\Bstrut

 \multirow{5}{*}{GaP$_{1-x}$Sb$_x$\footnote{For GaPSb samples, no specific thicknesses were reported in the reference.}} & GaP & 14.0 & $-$ & 1.66(I) & 1.61(I) & $-0.05$ & \multirow{5}{*}{0.10} \\ \nmidrule{2-7}
 & \multirow{3}{*}{GaAs} & 29.0 & $-$ & 1.33(I) & 1.39(I) & 0.06 & \\
 & & 32.0 & $-$ & 1.30 & 1.31 & 0.01 & \\
 \cite{Jou1988} & & 37.0 & $-$ & 1.24 & 1.33 & 0.09 &  \\ \nmidrule{2-7}
 & GaSb & 93.0 & $-$ & 0.56 & 0.74 & 0.18 &  \\ \hline\Tstrut\Bstrut

\multirow{2}{*}{GaAs$_{1-x}$Sb$_x$} & \multirow{2}{*}{GaAs} & 5.5 & 46.3 & 1.22 & 1.34 & 0.12 & \multirow{2}{*}{0.13} \\
 & & 7.0 & 51.2 & 1.17 & 1.31 & 0.14 &  \\ \hline\Tstrut\Bstrut

 \multirow{8}{*}{GaAs$_{1-x}$Bi$_x$} & \multirow{8}{*}{GaAs} & 0.9 & 75.0 & 1.28 & 1.33 & 0.05 & \multirow{8}{*}{0.10} \\
 & & 1.9 & 67.0 & 1.19 & 1.26 & 0.07 &  \\
 & & 2.9 & 60.0 & 1.10 & 1.20 & 0.10 &  \\
 & & 3.2 & 59.0 & 1.08 & 1.18 & 0.10 &  \\
 & & 3.8 & 54.0 & 1.04 & 1.14 & 0.10 &  \\
 \cite{Nattermann2017,Ludewig2013,Lewis2012} & & 4.8 & 25.0 & 0.98 & 1.11 & 0.13 &  \\
 \cite{Masnadi-Shirazi2014,Ludewig2014GrowthMOVPE,Mohmad2011TheGaAs1xBix}& & 5.3 & 50.0 & 0.95 & 1.07 & 0.12 &  \\
 & & 6.0 & 25.0 & 0.91 & 1.04 & 0.13 &  \\ \hline\Tstrut\Bstrut

 \multirow{5}{*}{GaAs$_{1-x}$N$_x$} & \multirow{5}{*}{GaAs} & 1.2 & 6.3 & 1.20 & 1.25 & 0.05 & \multirow{5}{*}{0.12} \\
 & & 2.0 & 17.0 & 1.10 & 1.16 & 0.06 &  \\
 & & 2.3 & 7.0 & 1.06 & 1.17 & 0.11 &  \\
 \cite{Vurgaftman2001Band,Rosenow2018Ab,Ludewig2017}& & 2.9 & 7.0 & 1.00 & 1.11 & 0.11 &  \\
 & & 5.0 & 4.0 & 0.82 & 1.01 & 0.19 &  \\ \hline\Tstrut\Bstrut
 \multirow{2}{*}{GaAs$_{1-x}$Sb$_x$} & \multirow{2}{*}{GaAs} & 27.8 & 3.7 & 0.67 & 1.10 & 0.43 & \multirow{2}{*}{0.42} \\
 & & 28.0 & 4.1 & 0.66 & 1.07 & 0.41 &  \\ \hline\Tstrut\Bstrut
 \multirow{5}{*}{GaAs$_{1-x}$N$_x$} & GaP & 4.9 & 6.0 & 0.66 & 1.18 & 0.52 & \multirow{5}{*}{0.66} \\ \nmidrule{2-7}
 & \multirow{4}{*}{Si} & 6.9 & 5.5 & 0.47 & 1.21 & 0.74 & \\
 & & 8.9 & 5.5 & 0.46 & 1.17 & 0.71 & \\
 \cite{Vurgaftman2001Band,Rosenow2018Ab,Ludewig2017}& & 9.5 & 6.0 & 0.45 & 1.11 & 0.66 &  \\
 & & 10.9 & 5.4 & 0.45 & 1.11 & 0.66 &  \\ \hline\Tstrut\Bstrut
 GaP$_{1-x}$Sb$_x$ & Si & $\leq 17.5$ & 7--9 & 1.91--1.36(I) & \multicolumn{3}{c}{\textrm{No RT-PL observed}} \\ \hline\Tstrut\Bstrut
 GaP$_{1-x}$Bi$_x$ & \multirow{2}{*}{GaP} & \multirow{2}{*}{$\leq 12.0$} & \multirow{2}{*}{17--73} & \multirow{2}{*}{2.03--1.18(I)} & \multicolumn{3}{c}{\multirow{2}{*}{No RT-PL observed}} \\
 \cite{Nattermann2017} & & & & &\multicolumn{3}{c}{}  \\
\end{tabular}
\end{ruledtabular}
\end{table*}

The excellent and consistent agreement between the experiment and calculated bandgaps (both in magnitude and in nature) suggests that we are able to quantitatively predict the bandgap over a wide range of compounds, compositions, and strain regions. However, as discussed above, the effect of 2D confinement is also crucial for relatively thin quantum well heterostructures and, hence, needs further investigation.
\begin{figure*}
\centering
  \subfloat[]{\label{fig:fig7a}\includegraphics[width=7in]{figure7a.png}}
  \hfill
  \subfloat[]{\label{fig:fig7b}\includegraphics[width=7in]{figure7b.png}}
  \hfill
  \subfloat[]{\label{fig:fig7c}\includegraphics[width=7in]{figure7c.png}}
  \caption{Proposals on how the bandgap phase diagram of biaxially strained GaAsP can be used in designing optoelectronic devices. (a) Defines the bound of composition region for creating a quantum-well heterostructure with direct bandgap GaAsP on GaAs substrate. (b) Choosing the different composition regions appropriately to make a multi-junction heterostructure with successive direct and indirect cells on the GaAs substrate. (c) In the vicinity of the transition point, the bandgap properties of the GaAsP epilayer on the GaP substrate can be changed by appropriately varying the composition.}
\end{figure*}

\begin{figure*}\ContinuedFloat
\centering
  \subfloat[]{\label{fig:fig7d}\includegraphics[width=7in]{figure7d.png}}
  \caption{\label{fig:fig7}\textit{(Continued)} Proposals on how the bandgap phase diagram of biaxially strained GaAsP can be used in designing optoelectronic devices. (d) Depending on the choice of substrate, GaAs or Si, the particular composition indicated by the vertical line can be made direct or indirect bandgap, respectively.}
\end{figure*}

Finally, based on the bandgap phase diagram, we propose several design strategies to optimize the selection of material combinations for achieving specific optical applications and new design principles for devices (\cref{fig:fig7}).

In \cref{fig:fig7a}, we propose a quantum-well heterostructure (QWH) composed of biaxially strained GaAsP on GaAs substrate. As the QW layers are made out of a single material with varied composition only, the epitaxial growth can be performed efficiently. The bandgap phase diagram shows the areas in compositional phase space where a direct bandgap in GaAs$_{1-x}$P$_x$ can be achieved ($x <$ 34\%). For $x >$ 35\%, the bandgaps are indirect and hence, are inappropriate for the heterostructure.

Figure~\ref{fig:fig7b} shows an efficient approach for the monolithic integration of multiple QWH to construct multi-junction photovoltaics. In this case, the QWHs are separated by thin indirect bandgap layers of the same material as QWH but only with a different composition. This would make the integration approach efficient, as no sample transfer is required during growth.

In \cref{fig:fig7c}, we propose a device with a gradual change in the bandgap properties. The concept utilizes the continuous transition in the nature of bandgap with alloy concentration in the vicinity of the DIT region. At the amount of P chosen here ($x =$ 15--35\%), we grow the GaAsP epitaxial layer on GaP with P-concentration continuously changing from direct to indirect bandgap region or vice versa. This way, changes in the bandgap magnitude, as well as the nature of the bandgap, are possible. Note that the concentration gradient can be implemented both in the horizontal and vertical directions.

Figure~\ref{fig:fig7d} shows another application of this concept. By appropriately choosing the substrate, we can tune the epitaxial layer (here: GaAsP) to show either a direct or indirect bandgap. Depending on the substrate, GaAs, or Si, the particular composition indicated by the vertical line will show direct or indirect bandgap, respectively.

\section{\label{sec6:summary}Summary}
Using density functional theory and the concept of band unfolding, we developed a first-principle computational protocol for the exhaustive mapping of the bandgap magnitude and type over a wide range of composition and strain values for several ternary III-V semiconductors. We constructed the composition-strain-bandgap relationship, the bandgap phase diagram, for several ternary III-V semiconductors: GaAsP, GaAsN, GaPSb, GaAsSb, GaPBi, and GaAsBi. We showed that this way of mapping the effect of strain can be used to choose application-specific best-suited material systems and hence is highly beneficial to device design. In addition, we developed an efficient approach based on Bloch spectral density for determining the nature of bandgap from supercell calculation. Notably, our computational protocol can be generalized to explore the vast chemical space of III-V materials with all other possible combinations of III- and V-elements. The comparison to experimental bandgap data underlines the accuracy of the computational approach chosen. This approach will be extended to more complex materials in the future.

\begin{acknowledgments}
We thank the German Research Foundation (DFG) in the framework of the Research Training Group ``Functionalization of Semiconductors" (GRK 1782) for funding this project and to HRZ Marburg, GOETHE-CSC Frankfurt, ZIH Dresden, and HLR Stuttgart for providing the computational resources.
\end{acknowledgments}

\subsection*{\label{datavailabilitystatement}DATA AVAILABILITY STATEMENT}
The density functional theory calculations data are openly available in the NOMAD repository (\url{DOI: https://doi.org/10.17172/NOMAD/2023.02.27-1}). The interactive bandgap phase diagrams are available on GitHub (\url{https://bmondal94.github.io/Bandgap-Phase-Diagram/}, last accessed 24.02.2023).

\subsection*{\label{ordisids}ORCID \lowercase{i}D\lowercase{s}}
\noindent
Badal Mondal \\ \hspace*{10pt} \url{https://orcid.org/0000-0002-0522-1254}\\
Ralf Tonner-Zech \\ \hspace*{10pt} \url{https://orcid.org/0000-0002-6759-8559}\\
Thilo Hepp \\ \hspace*{10pt} \url{https://orcid.org/0000-0002-0782-987X}\\
Kerstin Volz \\ \hspace*{10pt} \url{https://orcid.org/0000-0002-4456-5439}

\pagebreak
\onecolumngrid
\subsection*{\label{supportinginformation}SUPPORTING INFORMATION}
The Supporting Information (SI) is available free of charge at \url{https://doi.org/....} The SI contains:
\begin{itemize}
    \item Effective bandstructures of GaAsP, GaAsN, and GaAsN under selected isotropic strain values.
    \item Determining uncertainty in bandgap nature close to direct-indirect transition region.
    \item The DFT calculated bandgap values for all the systems described in the manuscript.
    \item The details of the interpolation procedure.
\end{itemize}
\twocolumngrid

\bibliographystyle{apsrev4-2}

\section{\label{secSI:ebsgaasp}Effective bandstructures of G\lowercase{a}A\lowercase{s}$_{0.963}$P$_{0.037}$ under isotropic strain}
\begin{figure}[htb]
\centering
  \subfloat[Strain = 0.0\%]{\label{fig:figS1a}\includegraphics[width=3.4in]{figureS1a.png}}
  \hspace{0.2in}
  \subfloat[Strain = --1.3\%]{\label{fig:figS1b}\includegraphics[width=3.4in]{figureS1b.png}} \\
  \subfloat[Strain = --1.5\%]{\label{fig:figS1c}\includegraphics[width=3.4in]{figureS1c.png}}
  \hspace{0.2in}
  \subfloat[Strain = --2.0\%]{\label{fig:figS1d}\includegraphics[width=3.4in]{figureS1d.png}}
  \caption{Effective bandstructures of GaAs$_{0.963}$P$_{0.037}$ under isotropic strain. The bandstructures were plotted using `BandUP' \cite{Medeiros2015Sup,Medeiros2014Sup}. Figures~\ref{fig:figS1a} and \ref{fig:figS1b} are direct bandgaps. Figure~\ref{fig:figS1c} is the direct-indirect transition point. Figure~\ref{fig:figS1d} is indirect bandgap.The `+' and `--' strains correspond to the tensile and compressive strain, respectively. For simplicity, only the L$\rightarrow \Gamma \rightarrow $ X high symmetry path is shown.}
  \label{fig:figS1}
\end{figure}
\vspace*{\fill}
\pagebreak
\vspace*{\fill}
\section{\label{secSII:bandgapnature}Determining the nature of bandgap using Bloch spectral weights}
\begin{figure}[h]
    \centering
    \includegraphics[width=7in]{figureS2.png}
    \caption{Variation of the bandgap (E$_{\text{g}}$) under biaxial tensile strain for GaAs$_{0.63}$P$_{0.37}$. The Bloch spectral weights (BSW) are given for the conduction band at the $\Gamma$ and L or X point [$\Gamma$ : (L or X)]. The vertical lines separate regions where the CBM changes character. As we are only interested in the direct and indirect character of the bandgap, here, we do not explicitly show the L and X BSWs.}
    \label{fig:figS2}
\end{figure}
\vspace*{\fill}
\pagebreak
\vspace*{\fill}
\section{\label{secSIII:ebsgaasn}Effective bandstructures of G\lowercase{a}A\lowercase{s}$_{0.995}$N$_{0.005}$ under isotropic strain}
\begin{figure}[h]
\centering
  \subfloat[Strain = 0.0\%]{\label{fig:figS3a}\includegraphics[width=3.4in]{figureS3a.png}}\\
  \subfloat[Strain = --4.5\%]{\label{fig:figS3b}\includegraphics[width=3.4in]{figureS3b.png}}
  \hspace{0.2in}
  \subfloat[Strain = +4.5\%]{\label{fig:figS3c}\includegraphics[width=3.4in]{figureS3c.png}}

  \caption{Effective bandstructures of GaAs$_{0.995}$P$_{0.005}$ under isotropic strain. The bandstructures were plotted using `BandUP' \cite{Medeiros2015Sup,Medeiros2014Sup}. Figures~\ref{fig:figS3a} and \ref{fig:figS3c} are direct bandgaps. In Fig.~\ref{fig:figS3b} the Bloch spectral weight of the ‘N-defect state’ at the conduction band $\Gamma$-point is very small, indicated by red circle. Therefore, Fig.~\ref{fig:figS3b} is categorized as an indirect bandgap. The `+' and `--' strains correspond to the tensile and compressive strain, respectively. For simplicity, only the L$\rightarrow \Gamma \rightarrow $ X high symmetry path is shown.}
  \label{fig:figS3}
\end{figure}
\vspace*{\fill}
\clearpage
\vspace*{\fill}
\begin{figure}
\centering
  \subfloat[Strain = 0.0\%]{\label{fig:figS4a}\includegraphics[width=3.4in]{figureS4a.png}}
  \hspace{0.2in}
  \subfloat[Strain = --5.0\%]{\label{fig:figS4b}\includegraphics[width=3.4in]{figureS4b.png}}
  \caption{Effective bandstructures of GaAs$_{0.889}$P$_{0.111}$ under isotropic strain. The bandstructures were plotted using `BandUP' \cite{Medeiros2015Sup,Medeiros2014Sup}. The bands are strongly dispersed for high N concentrations. The `+' and `--' strains correspond to the tensile and compressive strain, respectively. For simplicity, only the L$\rightarrow \Gamma \rightarrow $ X high symmetry path is shown.}
  \label{fig:figS4}
\end{figure}

\vspace*{\fill}
\pagebreak
\vspace*{\fill}
\section{\label{secSIV:bandredefinition}Conduction band redefinition based on minimum cut-off Bloch spectral weights criterion}
\begin{figure}[h]
    \centering
    \includegraphics[width=7in]{figureS5.png}
    \caption{Variation of the bandgap (E$_{\text{g}}$) under biaxial compressive strain for GaAs$_{0.995}$N$_{0.005}$. The folded supercell conduction bands' $\Gamma$-, L-, and X-BSWs are given at each point. The blue and orange lines are the E$_{\text{g}}$ with and without redefinition of the CBM. 20\% BSW is the cut-off criterion for redefining an eigenstate.}
    \label{fig:figS5}
\end{figure}
\vspace*{\fill}
\pagebreak
\vspace*{\fill}
\section{\label{secSV:dftvalues}The density functional theory bandgap values (without interpolation)}

\begin{figure}[h!]
\centering
  \subfloat{\label{fig:figS6_}\includegraphics[scale=0.6]{figureS6_legends.png}}\\
  \setcounter{subfigure}{0}
  \subfloat[Isotropic strain for GaAsP]{%
    \includegraphics[width=3.4in]{figureS6a_1.png}
    \hspace{0.2in}
    \includegraphics[width=3.4in]{figureS6a_2.png}
    \label{fig:figS6a}
   } \\ %
   \subfloat[Biaxial strain for GaAsP]{%
    \includegraphics[width=3.4in]{figureS6b_1.png}
    \hspace{0.2in}
    \includegraphics[width=3.4in]{figureS6b_2.png}
    \label{fig:figS6b}
   }

\end{figure}
\vspace*{\fill}
\clearpage
\vspace*{\fill}
\addtocounter{figure}{1}
\begin{figure}[h!]\ContinuedFloat
\centering
   \subfloat[Isotropic strain for GaAsN]{%
    \includegraphics[width=3.4in]{figureS6c_1.png}
    \hspace{0.2in}
    \includegraphics[width=3.4in]{figureS6c_2.png}
    \label{fig:figS6c}
   } \\%
   \subfloat[Biaxial strain for GaAsN]{%
    \includegraphics[width=3.4in]{figureS6d_1.png}
    \hspace{0.2in}
    \includegraphics[width=3.4in]{figureS6d_2.png}
    \label{fig:figS6d}
   } \\ %
   \subfloat[Biaxial strain for GaAsSb]{%
    \includegraphics[width=3.4in]{figureS6e_1.png}
    \hspace{0.2in}
    \includegraphics[width=3.4in]{figureS6e_2.png}
    \label{fig:figS6e}
   }

\end{figure}
\addtocounter{figure}{1}
\begin{figure}[h!]\ContinuedFloat
\centering
  \subfloat[Biaxial strain for GaAsBi]{%
    \includegraphics[width=3.4in]{figureS6f_1.png}
    \hspace{0.2in}
    \includegraphics[width=3.4in]{figureS6f_2.png}
    \label{fig:figS6f}
   } \\%
   \subfloat[Biaxial strain for GaPSb]{%
    \includegraphics[width=3.4in]{figureS6g_1.png}
    \hspace{0.2in}
    \includegraphics[width=3.4in]{figureS6g_2.png}
    \label{fig:figS6g}
   } \\ %
   \subfloat[Biaxial strain for GaPBi]{%
    \includegraphics[width=3.4in]{figureS6h_1.png}
    \hspace{0.2in}
    \includegraphics[width=3.4in]{figureS6h_2.png}
    \label{fig:figS6h}
   } %

  \caption{Effect of strain on the bandgap of GaE$_{1-x}$Y$x$ (E=As, P; Y=P, N, Sb, Bi). The left column shows the variation of bandgap magnitudes (E$_{\text{g}}$) as a function of strain at different compositions for the systems. The composition ($x\%$) legends are shown at the top. In the right column, we project the E$_{\text{g}}$ on the composition-strain space. The E$_{\text{g}}$ from these figures is interpolated to create Figs. 2--6 of the manuscript. We chose the same color bar scale limit for all the systems for consistency. The black crosses in (e) and (f) indicate the bandgap calculations that did not converge.}
  \label{fig:figS6}
\end{figure}

\clearpage

In Fig.~\ref{fig:figS6}, for the binary compounds (GaAs, GaP, and GaSb), we used the bandgap values from the `III-V binary semiconductors strain study' dataset of Ref.~\cite{Mondal2022Sup}. The strain axes were sampled primarily on the regular interval of 0.5\%. For GaAsP and GaPSb systems, sampling resolutions were increased close to the direct-indirect transition regions, as explained in Sec. III of the manuscript. For binaries, the intervals were not regular in the reference dataset. The composition axes were sampled in non-regular intervals. Because of the finite supercell size, the compositions could not be sampled in regular intervals. Note, however, the choice of compositions was arbitrary. Also, some of the DFT calculations at high strain could not be converged to the required accuracy and hence, are excluded. For GaAsN under isotropic strain because of the effective conduction band redefinition (Fig.~\ref{fig:figS5}), as explained in Sec. III of the manuscript, the sudden rapid increase in bandgap values under compressive strain are noticeable (Fig.~\ref{fig:figS6c}). The sudden increase in band energies due to band redefinition are also visible for compressively strained GaPSb at 25\% Sb concentration (Fig.~\ref{fig:figS6g}).
\vspace*{\fill}
\section{\label{secSVI:interpolationdetails}Interpolation details}
The E$_{\text{g}}$ from the right column figures of Fig.~\ref{fig:figS6} were interpolated over a fine grid to create Figs. 2--6 of the manuscript. We chose 0.1 grid resolution in both composition and strain for all the systems. We used Clough-Tocher piecewise cubic, C1 smooth, curvature-minimizing 2D interpolator from python \verb|scipy.interpolate| class \cite{SciPyNMeth2020}. The interpolation scheme follows:
\begin{itemize}
    \item[Step 1:] The input data are first triangulated with \verb|Qhull| \cite{Www.qhull.org} (Fig.~\ref{fig:figS7b}).
    \item[Step 2:] The interpolant then constructs a piecewise cubic interpolating B{\'e}zier polynomial on each triangle using the Clough-Tocher scheme \cite{Alfeld1984,Farin1986} (Fig.~\ref{fig:figS7c}). The interpolant is guaranteed to be continuously differentiable. The interpolant gradients are determined such that the interpolating surface curvature is approximately minimized \cite{Nielson1983,Renka1984}.
    \item[Step 3:] In the end, nearest-neighbor extrapolation is performed on the remaining points that can not be covered via \verb|Qhull| (Fig.~\ref{fig:figS7d}). Note that mostly the boundary points beyond the $x$ and $y$ limits of input points can not be triangulated under \verb|Qhull|. As long as those points are not far away from the available input points, the nearest-neighbor extrapolations to those points are acceptable.
    \item[Step 4:] Finally, bicubic interpolation, as implemented in python \verb|matplotlib.pyplot.imshow| class \cite{Hunter:2007,thomas_a_caswell_2022_7084615} is applied to smoothen the grid image (Fig.~\ref{fig:figS7e}).
\end{itemize}
Note both the input data and interpolation grid are rescaled to unit square during interpolation to eliminate the effect of incommensurable units and/or large differences between the $x$ \& $y$ scales.

\begin{figure}
\centering
  \subfloat[]{\label{fig:figS7a}\includegraphics[width=3.4in]{figureS7a.png}}
  \hspace{0.2in}
  \subfloat[]{\label{fig:figS7b}\includegraphics[width=3.4in]{figureS7b.png}} \\
  \subfloat[]{\label{fig:figS7c}\includegraphics[width=3.4in]{figureS7c.png}}
  \hspace{0.2in}
  \subfloat[]{\label{fig:figS7d}\includegraphics[width=3.4in]{figureS7d.png}} \\
  \subfloat[]{\label{fig:figS7e}\includegraphics[width=3.4in]{figureS7e.png}}
  \caption{The Clough-Tocher piecewise cubic, C1 smooth, curvature-minimizing 2D interpolation scheme. The color dots are the input points. The black dots indicate the interpolation grid. (a) Input data. (b) Delaunay triangulation, indicated by black triangles. (c) Nearest-neighbor interpolation points, indicated by the black stars. These points could not be triangulated. (d) Nearest-neighbor interpolation. (e) Bicubic interpolation over (d).}
  \label{fig:figS7}
\end{figure}
\vspace*{\fill}
\clearpage