\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

