Introduction

The prediction of equilibrium crystal structures has long been a challenge in materials science. Atomic bonding, which is determined by a variety of electronic and steric factors, governs crystal structure formation1. Among the various approaches developed, phenomenological models have been powerful tools for crystal structure prediction2. Unlike physical models that rely on costly first-principles calculations, phenomenological models are built on empirical data and excel at reproducing known results and predicting unknown phenomena3,4.

Prominent phenomenological models, such as the Hume–Rothery rules, predict the solubility and crystal structures of metallic elements based on their atomic size and electronegativity differences5. Large electronegativity differences favor compound formation, whereas similar atomic sizes and electronegativities with small valence differences increase solubility. While these rules were primarily developed for binary metallic alloys, they may not fully capture the compositional asymmetry and charge redistribution effects observed in ternary systems. Miedema’s rules6, another key model, predicted alloy formation by analyzing enthalpy changes through work function and electron density differences, favoring alloys with large work function differences and small electron density differences. Similarly, this model does not yet explicitly incorporate the electrostatic stabilization that can be important for predicting stable configurations in ternary compounds.

Visualization methods, such as Villar’s method7 and the Mooser–Pearson plot8 provide insights into the phase stability. Villar’s method uses atomic size, valence, and electronegativity differences to map the AB compound stability, whereas the Mooser–Pearson plot predicts AX compound configurations based on electronegativity and atomic size, distinguishing octahedral from tetrahedral coordination. Although both approaches have been effective for binary systems, they do not fully capture the composition-sensitive bonding variations inherent to multicomponent compounds.

Pauling’s rules stand out among bond polarity-based models9, defining cation (Rule 1) and anion (Rule 2) coordination numbers, cation–cation repulsion (Rule 3), and highlighting the separation of cations into distinct crystallographic planes (Rule 4), which are critical for understanding cation-eutaxy ABX ternary systems. In particular, Rule 4 is crucial for understanding the cation-eutaxy ABX ternary system, as high-valence, low-coordination cations tend to avoid polyhedral element sharing due to electrostatic repulsion. This drives the formation of layered structures with distinct cationic planes, such as LiMnO210, where strong cation repulsion stabilizes the layered configuration and facilitates Li-ion diffusion, a key feature in battery applications.

Cation-eutaxy ABX ternary structures are defined as layered structures in which cations preferentially occupy eutactic positions, forming a well-ordered arrangement. This structural motif, driven by strong cation-cation repulsion, has garnered significant attention as a platform for designing novel layered materials11. Topochemical etching selectively removes A-site cations from ABX ternary compounds to produce metastable BX binary layered materials12,13,14, thus expanding the van der Waals (vdW) library. Fewer than 100 BX binary materials have been synthesized using this method15, there is immense potential for growth in this field.

Our previous work advanced this approach by partially retaining the A-site cations to maintain structural stability, leading to layered A1-xBX ternary materials16. By exploiting the A-site cation vacancies created during topochemical etching, we demonstrated memristive behavior in a class of III‒V-layered semiconductors driven by a vacancy diffusion mechanism. These properties, which are absent in bulk materials, highlight the transformative potential of cationic eutaxy systems. These approaches are still in the early stages of investigation, making the identification of suitable ternary compounds an urgent research priority. Furthermore, as the study of multicomponent materials continues to gain attention, there is an increasing need for predictive models that, along with Pauling’s fourth rule, can accurately assess the phase stability of ternary systems.

In this study, we present a phenomenological model to predict the stability of cation-eutaxy ABX ternary compounds. This model categorizes cationic eutaxy compounds within [A‒III‒V] chemical systems by assessing the covalent nature of the III‒V bonds through key atomic properties, including electronegativity, ionic radius, and charge. Based on computational high-throughput screening and the Materials Project (MP) database16, we validated our proposed model using 35 known cation-eutaxy ABX compounds and nine previously unidentified candidates. Our model demonstrated high predictive accuracy by incorporating the atomic size, electronegativity, and charge. Furthermore, the experimental realization of the previously unexplored compounds, K2In2P3 and Na2In2As3, validated the model and provided key insights for guiding the discovery of layered materials and their potential applications.

Classification of cation-eutaxy structures

We propose a topological framework to classify ternary A‒III‒V compounds based on the dimensionality of the interconnected III‒V polyhedra within their crystal structures (Fig. 1a). III‒V compounds can be categorized as zero-dimensional (0D), one-dimensional (1D), two-dimensional (2D), or three-dimensional (3D) III‒V networks. Representative compounds are illustrated based on the crystallographic information files from the MP database17: 0D (Ca14AlSb11, \(I{4}_{1}/{acd}\)), 1D (Ca5(AlSb3)2, \({Pbam}\)), 2D (CaIn2P2, \(P{6}_{3}/{mmc}\)), and 3D (MgB9N, \(R\bar{3}m\)) (see Supplementary Data 1).

Fig. 1: Topological classification of ternary A‒III‒V compounds based on bonding characteristics.
figure 1

a Structure classification of A‒III‒V compounds according to the dimensionality of III‒V network. Compounds with 0D and 1D III‒V networks incorporate a 3D A‒V networks, whereas those with 2D III‒V networks adopt alternating stacked structures with 2D A‒V networks. In contrast, ternary compounds containing 3D III‒V networks include 0D or 1D A‒V networks. The structural conditions and representative structural schematics of cation-eutaxy structures are summarized below. b Phenomenological model for classifying III‒V network dimensionality in ternary compounds using descriptors relative \({E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\) (α) and relative \({E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\) (β). Bonding characteristics are quantified by Econf,avg and Econf,diff, displayed in a modified van Arkel-Ketelaar triangle (top left inset). Ternary compounds with 0D, 1D, and 3D III‒V networks are marked as blue circles, black inverted triangles, and green diamonds, respectively. Known cation-eutaxy compounds (2D III‒V network) are represented by red squares, while unknown cation-eutaxy candidates are shown as orange stars. Regions with dense populations of specific dimensionality are shaded in corresponding background colors. Known cation-eutaxy compounds are labeled with numerical indices and their chemical formulas listed in the bottom right inset, while unknown candidates are directly labeled with their chemical formulas in the top right inset. c Histogram showing the composition of element A and element III in ternary compounds. Weighted compositions, calculated based on stoichiometric ratios, are represented by solid bars for element A and hatched bars for element III.

Full size image

In particular, cation-eutaxy ternary compounds exhibit 2D III–V networks within their structures. These compounds are characterized by a unique cation ordering, in which cations A and III occupy separate planes to minimize electrostatic repulsion. The structural conditions defining cation-eutaxy A‒III‒V compounds are as follows (bottom, Fig. 1a): (1) The ionicity of the A‒V bond: cation A must exhibit stronger ionic character in the A‒V bond than in the III‒V bond. This requires cation A to have a lower electronegativity than cation III, facilitating ionic bonding with V and strengthening the polar covalent nature of the III‒V bonds. (2) Two-dimensional polyhedral connectivity: the III-V polyhedra are interconnected, similar to 2D crystals. This connectivity involves primary III–V bonding and polyanionic interactions, resulting in a layered [III–V] framework rather than a rod-like structure. (3) Distinct cation planes: cations A and III are positioned in separate planes, with the [III‒V] layers sandwiched between layers composed exclusively of cation A. (4) Absence of interlayer primary bonding: along the z-axis, there are no primary bonds between the III‒V polyhedra, ensuring the structural separation of the layers. These criteria allow for the effective classification of cation-eutaxy structures. Unlike conventional structures defined by specific space groups or prototypes (Supplementary Fig. 1), cation-eutaxy structures encompass various configurations across various space groups and prototypes.

In a previous study, we employed a structure-matching algorithm to identify 35 thermodynamically stable (energy above the convex hull (Ehull) = 0 eV/atom) cation-eutaxy A‒III‒V compounds from the MP database16. These structures were further grouped into 15 prototypes (Supplementary Fig. 2 and Supplementary Table 1). By applying an ion substitution algorithm, we identified nine additional quasi-stable (Ehull < 0.2 eV/atom) cation-eutaxy A‒III‒V compounds that were not previously listed in the MP database16. This computational screening yielded 44 cation-eutaxy A‒III‒V compounds (Supplementary Table 1). In this study, we propose a phenomenological model and validate its accuracy using these datasets.

In ionic compounds, increasing the bond covalency often leads to the formation of more discrete structures18. Low covalency typically favors 3D close-packed structures, whereas higher covalency promotes the formation of low-dimensional architectures. Applying this principle to ternary A‒III‒V compounds, we hypothesized that increasing the covalency of the III‒V bonds would enhance the formation of cation-eutaxy structures. To test this hypothesis, we developed models to classify cation-eutaxy compounds in [A‒III‒V] chemical systems by evaluating the III‒V bond covalency using atomic properties such as electronegativity, ionic size, and charge. Unlike in binary systems, where bonding characteristics remain relatively constant, ternary systems exhibit composition-dependent variations due to local charge redistribution. To account for this effect, we introduced stoichiometry-weighted descriptors that systematically incorporate composition-dependent bonding variations. This analysis provides a robust framework for predicting cation-eutaxy structures and offers valuable insights for discovering layered materials through topochemical etching.

Van Arkel-Ketelaar triangle

The van Arkel-Ketelaar triangle traditionally evaluates bonding characteristics using the average electronegativity and electronegativity difference between two elements, which classify bonding as metallic, covalent, or ionic19. However, substituting traditional electronegativity scales with the configuration energy (Econf) provides a more accurate description of the bonding characteristics20. Econf is defined as the average ionization energy of the valence electrons in the ground states21 (Supplementary Fig. 3). Building on this, we utilize a modified van Arkel-Ketelaar triangle based on the average Econf (Econf,avg) and Econf difference (Econf,diff)22. Econf,avg represents the transition from metallic to covalent bonding, whereas Econf,diff reflects the balance between ionic and metallic/covalent bonding (top-left inset of Fig. 1b).

We hypothesized that cation-eutaxy structures were more likely to form when the III‒V network exhibited greater covalency than the A‒V network. To quantitatively analyze the relative Econf of the A‒V and III‒V bonds in ternary compounds, we introduce a approach that incorporates stoichiometry, an aspect not previously considered in traditional van Arkel-Ketelaar triangles. This addresses the limitations of defining bond characteristics based solely on the Econf of individual elements, which fails to account for the compositional influence within ternary systems.

To capture the bond characteristics more comprehensively, we propose using the stoichiometry-weighted configuration energy, \({E}_{{{\rm{conf}}}}^{{{\rm{w}}}}\), as a quantitative descriptor of bonding in ternary compounds. Specifically, the weighted average Econf is denoted as \({E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\) and the weighted difference is represented as \({E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\). These metrics were formulated for A‒V and III‒V bonds in AlIIImVn, where l, m, and n are the stoichiometric coefficients of cations A, III, and anion V, respectively. The expressions for \({E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\) and \({E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\) for (A‒V) and (III‒V) are as follows:

$${E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\left({{\rm{A}}}-{{\rm{V}}}\right)=\frac{(l\times {E}_{{{\rm{conf}}}}\left({{\rm{A}}}\right)+n\times {E}_{{{\rm{conf}}}}({{\rm{V}}}))}{(l+m+n)},$$
(1)
$${E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\left({{\rm{III}}}-{{\rm{V}}}\right)=\frac{(m\times {E}_{{{\rm{conf}}}}\left({{\rm{III}}}\right)+n\times {E}_{{{\rm{conf}}}}({{\rm{V}}}))}{(l+m+n)},$$
(2)
$${E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\left({{\rm{A}}}-{{\rm{V}}}\right)=\frac{n\times l\times ({E}_{{{\rm{conf}}}}\left({{\rm{V}}}\right)-{E}_{{{\rm{conf}}}}({{\rm{A}}}))}{({l}^{2}+{m}^{2}+{n}^{2})},{{\rm{and}}}$$
(3)
$${E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\left({{\rm{III}}}-{{\rm{V}}}\right)=\frac{n\times m\times ({E}_{{{\rm{conf}}}}\left({{\rm{V}}}\right)-{E}_{{{\rm{conf}}}}({{\rm{III}}}))}{({l}^{2}+{m}^{2}+{n}^{2})}$$
(4)

Accordingly, we define two descriptors based on the modified van Arkel-Ketelaar triangle: relative \({E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\) (α) and relative \({E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\) (β). α is calculated as the \({E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\) of (III‒V) relative to (A‒V) in compounds AlIIImVn. α is calculated as the ratio of \({E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\) for III‒V to A‒V bonds:

$$\alpha=\frac{{E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}({{\rm{III}}}{-}{{\rm{V}}})}{{E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}({{\rm{A}}}{-}{{\rm{V}}})}=\frac{(m\times {E}_{{{\rm{conf}}}}({{\mbox{III}}})+n\times {E}_{{{\rm{conf}}}}({\mbox{V}}))}{(l\times {E}_{{{\rm{conf}}}}({{\mbox{A}}})+n\times {E}_{{{\rm{conf}}}}({\mbox{V}}))}$$
(5)

Similarly, β is calculated as:

$$\beta=\frac{{E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}({{\rm{III}}}{-}{{\rm{V}}})}{{E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}({{\rm{A}}}{-}{{\rm{V}}})}=\frac{n\times m\times ({E}_{{{\rm{conf}}}}({{\mbox{V}}})-{E}_{{{\rm{conf}}}}({\mbox{III}}))}{n\times l\times ({E}_{{{\rm{conf}}}}({{\mbox{V}}})-{E}_{{{\rm{conf}}}}({\mbox{A}}))}$$
(6)

To prevent unphysical results, such as negative β values or suppressed variation across compositions, we applied a weighting scheme to \({E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\) that incorporates the product of stoichiometric coefficients and normalizes by their squared sum (\({l}^{2}+{m}^{2}+{n}^{2}\)). This formulation was essential for resolving composition-dependent bonding characteristics in ternary systems (see Methods and Supplementary Fig. 4).

Changes in α reflect the transition in the metallic and covalent nature of III–V bonds compared to A–V bonds. An increase in this value indicates a stronger covalent character of the III–V bonds. Similarly, changes in β represent the ionic and metallic/covalent nature shift. A decrease in this value suggests that the III–V bonds exhibit a more pronounced metallic/covalent character than the A–V bonds.

Changes in α and β play a significant role in determining the dimensionality of the III–V network (Supplementary Data 1). Specifically, an increase in α drives the network transition from 0D to 3D, while a decrease in β shifts the network from 3D to 0D (Fig. 1b and Supplementary Fig. 5). These findings suggest a linear relationship between α and β for the dimensionality of the network. When applied to 44 ternary A–III–V compounds, the α and β descriptors revealed a distinct clustering of cation-eutaxy structures within a specific parameter space. Notably, most cation-eutaxy structures were found to satisfy the criteria of α > 1 and β <1, highlighting their strong covalent nature in III–V bonds relative to A–V bonds. This trend indicates enhanced covalency in III–V bonds when α is high and β is low. A larger α and smaller β correspond to stronger covalent bonding in III–V relative to A–V, thereby promoting the formation of discrete, layered III–V networks. These findings are consistent with the van Arkel-Ketelaar triangle, where increased covalency is associated with a higher Econf,avg and a smaller Econf,diff. Such separation was not observed when unweighted descriptors were used (Supplementary Fig. 6a), indicating that stoichiometry-weighted definitions of α and β are essential for capturing the correlation between bonding character and structural dimensionality. Notably, the model demonstrated its robustness by accurately classifying all 35 previously known cation-eutaxy compounds and 9 newly identified ones, underscoring its potential as a reliable tool for predicting and categorizing cation-eutaxy structures.

In addition, the dimensionality of the III–V network was strongly influenced by the weighted compositions of A and III (Fig. 1c and Supplementary Data 1). First, as the composition of A decreases and that of III increases, the network transitions from 0D to 3D. Specifically, 0D III–V networks are observed when the weighted composition of A is significantly higher than that of III. In contrast, 3D networks emerge when the weighted composition of III far exceeds that of A. Second, a decrease in the weighted composition of A (i.e., low l) and a corresponding increase in the composition of III (i.e., high m) lead to an increase in α, which enhances the covalent nature of III–V bonds and drives the transition of the network from 0D to 2D. Similarly, the trend of increasing β aligns with the transition from 0D to 2D, as shown in Fig. 1b. However, when β > 1, the network shifts to 3D, consistent with the results presented in Fig. 1b. These trends highlight how the compositional balance between A and III dictates the dimensionality of the III–V network. The descriptors, α, and β, provide critical insights into how variations in electronegativity and composition influence the dimensionality of III–V networks, offering valuable guidance for designing and optimizing cation-eutaxy structures.

Effective ionic radius ratio

Next, we investigated how the ionic radius influences the bonding characteristics and whether it can serve as a criterion for classifying cation-eutaxy ternary systems. According to Pauling’s first rule, the ionic radius ratio (\(\frac{{r}_{{{\rm{c}}}}}{{r}_{a}}\)) between the cation and anion in a polyhedron determines the coordination number, which in turn affects bond polarization23,24. Specifically, a smaller \(\frac{{r}_{{{\rm{c}}}}}{{r}_{{{\rm{a}}}}}\) ratio corresponds to a lower coordination number, leading to an increase in bond polarization (top right inset of Fig. 2a). Given that bond polarization is proportional to bond covalency, the ionic radius ratio provides a means to evaluate the covalency of III–V bonds. Based on this principle, we propose a phenomenological model that utilizes the ionic radius ratio to assess the III–V bond covalency and effectively classify cation-eutaxy compounds.

Fig. 2: Phenomenological model for understanding cation-eutaxy structures in ternary A‒III‒V compounds based on ionic radius and ionic potential.
figure 2

a Model for classifying the dimensionality of III‒V network in ternary compounds using weighted ionic radius ratios (\(\frac{{{{\rm{r}}}}_{{{\rm{A}}}}^{{{\rm{w}}}}}{{{{\rm{r}}}}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) and \(\frac{{{{\rm{r}}}}_{{{\rm{III}}}}^{{{\rm{w}}}}}{{{{\rm{r}}}}_{{{\rm{V}}}}^{{{\rm{w}}}}}\)). Schematic representation of polyhedral structures and corresponding electron densities in a 2D projection (top right inset). As the ionic radius ratio of the cation to the anion increases, the coordination number increases, leading to a reduction in the electron cloud overlap of ions and a decrease in polarization. Compounds with 0D, 1D, and 3D III‒V networks are marked as blue circles, black inverted triangles, and green diamonds, respectively. Known cation-eutaxy compounds (2D III‒V network) are represented by red squares, while unknown cation-eutaxy candidates are shown as orange stars. Regions densely populated by compounds of specific dimensionalities are highlighted with matching background colors. Cation-eutaxy compounds are labeled numerically, with known compounds shown in red and unknown candidates in orange. The corresponding chemical formulas for each labeled compound are summarized in the bottom right inset. b Distribution of cation-eutaxy compounds by anion type. Nitride compounds (red) are the least common, with two examples, whereas phosphide (green), arsenide (blue), and antimonide (purple) are more prevalent, comprising 13, 17, and 12 compounds, respectively. c Histogram categorizing group I and II contributions to cation-eutaxy compounds, grouped by the anion V. d Relationship between ionic potential \({\varPhi }_{{{\rm{A}}}}\) and the ratio \(\frac{{\varPhi }_{{{\rm{V}}}}}{{\varPhi }_{{{\rm{III}}}}}\). Nitrides, phosphides, arsenides, and antimonides are represented in red, green, blue, and purple, respectively. The bottom right inset schematically illustrates the interplay between the surface charge density of III‒V layer and the ionic radius of the anion V.

Full size image

First, as the ionic radius ratio (\(\frac{{r}_{{{\rm{III}}}}}{{r}_{{{\rm{V}}}}}\)) decreased, the covalency of the III–V bonds increased, creating favorable conditions for the formation of cation-eutaxy structures. Furthermore, considering that cation-eutaxy ternary systems consist of alternating layers, it is favorable for the A‒V network to adopt a 2D network. Consequently, a decrease in the ionic radius ratio (\(\frac{{r}_{{{\rm{A}}}}}{{r}_{{{\rm{V}}}}}\)) between cations A and V further promotes the formation of cation-eutaxy structures.

To validate this hypothesis, we also introduced weighted ionic radius ratios (\(\frac{{r}_{{{\rm{A}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) and \(\frac{{r}_{{{\rm{III}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\)) as descriptors, considering the stoichiometric coefficients of cations A, III, and V in the ternary compounds. The ionic radii of A, III, and V are denoted rA, rIII, and rV, respectively25. The weighted ionic radii of A (\({r}_{{{\rm{A}}}}^{{{\rm{w}}}}\)), III (\({r}_{{{\rm{III}}}}^{{{\rm{w}}}}\)), and V (\({r}_{{{\rm{V}}}}^{{{\rm{w}}}}\)) in AlIIImVn were calculated as follows:

$${r}_{{{\rm{A}}}}^{{{\rm{w}}}}=\frac{l\times {r}_{{{\rm{A}}}}}{(l+m+n)},$$
(7)
$${r}_{{{\rm{III}}}}^{{{\rm{w}}}}=\frac{m\times {r}_{{{\rm{III}}}}}{(l+m+n)},{{\rm{and}}}$$
(8)
$${r}_{{{\rm{V}}}}^{{{\rm{w}}}}=\frac{n\times {r}_{{{\rm{A}}}}}{(l+m+n)}$$
(9)

The weighted ionic radius ratios are then determined as:

$$\frac{{r}_{{{\rm{A}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}=\frac{l\times {r}_{{{\rm{A}}}}}{n\times {r}_{{{\rm{V}}}}}$$
(10)
$$\frac{{r}_{{{\rm{III}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}=\frac{m\times {r}_{{{\rm{III}}}}}{n\times {r}_{{{\rm{V}}}}}$$
(11)

The analysis showed that as \(\frac{{r}_{{{\rm{A}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) decreased, the A–V network transitioned to a 2D structure. Similarly, smaller \(\frac{{r}_{{{\rm{III}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) values led to a reduction in the dimensionality of the III–V network, shifting from 3D to 2D, 1D, and eventually to 0D (Fig. 2a and Supplementary Data 1). This indicated an inverse relationship between \(\frac{{r}_{{{\rm{A}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) and \(\frac{{r}_{{{\rm{III}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) in terms of dimensionality (Fig. 2a). Overall, the dimensionality of the III–V network can be expressed as \(\frac{{r}_{{{\rm{A}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\propto 1/\frac{{r}_{{{\rm{III}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\), with clear boundaries between regions corresponding to 0D, 1D, and 2D structures (Fig. 2a and Supplementary Fig. 7). This classification was not achievable when unweighted ionic radius ratios were used (Supplementary Fig. 6b), underscoring the essential role of stoichiometry-weighted descriptors in resolving the dimensionality trends. Notably, 2D III–V networks were predominantly observed at relatively small values for both \(\frac{{r}_{{{\rm{A}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) and \(\frac{{r}_{{{\rm{III}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) (Fig. 2a). This finding suggests that both the III–V, and A–V networks exhibit high covalency, which is a necessary condition for attaining cation-eutaxy structures.

Analysis based on the ionic radius provides valuable insights into the selection of anions (V) necessary for the formation of cation-eutaxy structures. Specifically, we observed a striking disparity in the number of cation-eutaxy structures among the different anions: only 2 nitrides, compared to 13 phosphides, 17 arsenides, and 12 antimonides (Fig. 2b). This scarcity of nitrides can be attributed to the exceptionally small ionic radius of nitrogen, which significantly reduces the covalency of the III‒V bonds, thereby increasing their ionicity. Furthermore, nitrogen’s high Econf appeared to contribute to the reduced covalency of the III‒N bonds (Supplementary Fig. 8).

Next, we examined the role of charge in influencing the electrostatic interactions in the cation-eutaxy structure. Our analysis revealed a correlation between the ionic radius of anion V and the preference for cation A. Specifically, smaller radii of V favor the formation of cation-eutaxy structures with group II cations as A, whereas larger radii of V are associated with a preference for group I cations as A (Fig. 2c). For example, in cation-eutaxy nitrides, all the A cations belong to Group II. However, as the ionic radius of V increases from phosphorus to Sb, the proportion of group I cations in A increases, culminating in antimonides, where A predominantly consists of group I cations. This relationship can be attributed to the balance between the electrostatic repulsion within the 2D III‒V networks and the shielding effect the A cations provide in the cation-eutaxic A‒III‒V structures (bottom-right inset of Fig. 2d). In 2D III‒V networks, the anion V is positioned on the surface and carries surface charges, whereas cation A is situated between the III‒V layers, maintaining structural stability. An analysis of the ionic potential (Φ), which represents the charge density at the surface of an ion, demonstrates that the balance between electrostatic repulsion across III‒V layers and the shielding effect of A is critical for forming cation-eutaxy structures.

The ionic potential (Φi) of an ion can be defined as \(\frac{{w}_{i}\times {n}_{i}}{{r}_{i}}\), where wi represents the content of ion i, ni is the charge, and ri is the ionic radius26. We note that this heuristic is associated with the effective potential of an ion’s electrostatic influence, rather than a direct physical potential. As ΦV increases and ΦIII decreases, the high charge density of anion V becomes concentrated on the surface of the layer, leading to an increase in the surface charge density of the III‒V layers. Consequently, a higher \(\frac{{\varPhi }_{{{\rm{V}}}}}{{\varPhi }_{{{\rm{III}}}}}\) ratio results in a greater negative charge on the layer surface, amplifying the electrostatic repulsion between the III‒V layers. Conversely, an increase in ΦA enhances the shielding effect of cation A, counteracting this repulsion27. As shown in Fig. 2d, a clear proportional relationship exists between \(\frac{{\varPhi }_{{{\rm{V}}}}}{{\varPhi }_{{{\rm{III}}}}}\) and ΦA in cation-eutaxy A‒III‒V structures (Supplementary Table 1). This relationship underscores the critical role of sufficient shielding by cation A in stabilizing the III–V layers when the surface charge density is high.

Unknown cation-eutaxy within A2In2V3 prototype system

In the previous analysis, cation-eutaxy structures were classified based on atomic features, and surface charge considerations determined the preference for the specific cations A. Building on this, we integrated electronegativity, ionic radius, and charge into a unified classification framework.

We focused on the A2III2V3 prototype structures (top left inset of Fig. 3a) within the [A–In–V] chemical space, where indium serves as a group-III element. Among the 15 prototype structures identified across the 44 predicted cation-eutaxy compounds, the A₂III₂V₃ type was the most frequently encountered, with In-based compositions comprising the largest subset. We therefore selected the A₂In₂V₃ system as a representative case for validating the predictive power of the model through experimental synthesis. The results revealed the following key findings by plotting cation-eutaxy A2In2V3 compounds using the relative \({E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\) (β), and the weighted ionic radius ratio of In over V (\(\frac{{r}_{{{\rm{In}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\)). First, we identified clear boundaries between cation-eutaxy and non-cation-eutaxy structures at β = 0.42 and \(\frac{{r}_{{{\rm{In}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) = 0.26 (Fig. 3a and Supplementary Table 2). In A2In2V3 cation-eutaxy structures, charge neutrality requires cation A to be an alkali metal. However, for A2In2N3, the high \(\frac{{\varPhi }_{{{\rm{N}}}}}{{\varPhi }_{{{\rm{In}}}}}\) ratio results in significant electrostatic repulsion between In2N3 layers, necessitating a cation A with high ΦA, such as those found in alkaline earth metals, to provide sufficient shielding. As a result, A2In2N3 structures with alkali metal cations are structurally unstable and classified as non-cation-eutaxy structures. Third, our analysis confirms that lower β and \(\frac{{r}_{{{\rm{In}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) values favor the formation of cation-eutaxy A2In2V3. These results highlight the effectiveness of using an integrated approach combining electronegativity, ionic radius, and charge to classify cation-eutaxy structures and provide strong evidence for the feasibility and accuracy of the proposed phenomenological model.

Fig. 3: Discovery of unknown K2In2P3 and Na2In2As3 compounds.
figure 3

a Combined phenomenological model for [A2–In2–V3] chemical systems. The background color highlights distinct regions of predicted structural behavior: the blue area corresponds to non-cation-eutaxy compounds, while the brown area encompasses the predicted cation-eutaxy region, including both known and unknown candidates. The top left inset shows a simulated image of a representative cation-eutaxy A2In2V3 structure, generated from a crystallographic information file. b Thermodynamic stability of cation-eutaxy compounds within this family, with Ehull values of Li2In2Sb3, Na2In2Sb3, K2In2Sb3, Na2In2As3, K2In2As3, and K2In2P3 calculated as 0.197, 0.010, 0.000, 0.043, 0.000, and 0.022 eV/atom, respectively. The variations in Ehull due to ion substitution of anion V and cation A are highlighted with brown and green arrows, respectively. c, d HAADF-STEM images and corresponding SAED patterns (insets) of K2In2P3 (c) and Na2In2As3 (d), viewed along the [010] zone axis. The electron diffraction patterns of (002) and (100) are indexed by yellow and blue circles, respectively. e and f ABF-STEM images of K2In2P3 (e) and Na2In2As3 (f) along the [010] zone axis. g, h Enlarged ABF-STEM images of K2In2P3 (g) and Na2In2As3 (h) showing agreement with crystallographic information file-based structural models.

Full size image

We further investigated structural stability by calculating the formation enthalpies (Fig. 3b). For example, in the transitions from Na2In2Sb3 (Ehull = 0.010 eV/atom) to Na2In2As3 (Ehull = 0.043 eV/atom) and from K2In2As3 (Ehull = 0.0 eV/atom) to K2In2P3 (Ehull = 0.022 eV/atom), we observed a trend of increasing Ehull as the size of anion V decreased. Similarly, when K2In2Sb3 (Ehull = 0.0 eV/atom) and K2In2As3 (Ehull = 0.0 eV/atom) were replaced by Li2In2Sb3 (Ehull = 0.197 eV/atom) and Na2In2As3 (Ehull = 0.043 eV/atom), respectively, a smaller cation A was associated with a higher Ehull. This increase in the formation enthalpy indicates reduced thermodynamic stability, adversely affecting the formation of the cation-eutaxy structure. These trends align with the relationships depicted in Fig. 3a, where β and \(\frac{{r}_{{{\rm{In}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) are inversely proportional to stability. For instance, K2In2As3 is classified as a stable, known cation-eutaxy structure, whereas Na2In2As3 is metastable, and Li2In2As3 is a non-cation-eutaxy structure. Similarly, K2In2As3 has a stable cation-eutaxy structure, K2In2P3 is metastable, and K2In2N3 has a non-cation-eutaxy structure. In conclusion, substituting cations A and V with smaller ions reduces the likelihood of the formation of cation-eutaxy structures.

To rigorously assess the reliability of our model, we selected K₂In₂P₃ and Na₂In₂As₃, which were identified as candidate cation-eutaxy compounds with β and \(\frac{{r}_{{{\rm{In}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) values near the classification threshold. This selection allowed us to evaluate the model’s accuracy in distinguishing cation-eutaxy from non-cation-eutaxy compounds. Plate-shaped K2In2P3 and Na2In2As3, synthesized via solid-state methods, were characterized using scanning electron microscopy (SEM), which confirmed their morphologies (Supplementary Fig. 9). The elemental composition of K₂In₂P₃ and Na₂In₂As₃ was confirmed through energy-dispersive spectroscopy (EDS) analysis (Supplementary Fig. 9). X-ray diffraction (XRD) analysis further demonstrated that these compounds possessed crystal structures identical to those predicted by the simulated structures (Supplementary Fig. 10). The calculated XRD patterns were derived from ion substitution algorithm-predicted structures and relaxed using generalized gradient approximation (GGA)28. To resolve the atomic structures precisely, scanning transmission electron microscopy (STEM) analysis was performed on both unknown compounds. High-angle annular dark-field STEM (HAADF-STEM) imaging using focused ion beam (FIB)-prepared specimens revealed distinct atomic arrangements of the In2P3 and In2As3 layers, consistent with cation-eutaxy A2III2V3 structures (Fig. 3c, d). Additionally, selected-area electron diffraction (SAED) patterns confirmed that K2In2P3 and Na2In2As3 exhibited monoclinic P-lattice symmetry (insets in Fig. 3c, d). To ensure the accurate detection of light elements such as K and Na, annular bright-field STEM (ABF-STEM) analysis was conducted, which unequivocally confirmed that the synthesized compounds aligned perfectly with the atomic arrangements predicted for cation-eutaxy A2III2V3 structures (Fig. 3e–h).

Furthermore, we attempted to synthesize Li₂In₂P₃, Li₂In₂As₃, and Na₂In₂P₃, which exhibited β and \(\frac{{r}_{{{\rm{In}}}}^{{{\rm{w}}}}}{{r}_{{{\rm{V}}}}^{{{\rm{w}}}}}\) values similar to K₂In₂P₃ and Na₂In₂As₃ but were classified as non-cation-eutaxy compounds (Supplementary Fig. 11). XRD analysis revealed that the synthesized materials consisted of zinc-blende III–V phases along with either a ternary phase with a non-cation-eutaxy structure or highly reactive binary phases, depending on the chemical system. This outcome strongly supports the validity of the descriptor threshold established for distinguishing cation-eutaxy structures within the A₂In₂V₃ family, reinforcing the predictive power of our model.

Our phenomenological model proved to be effective in classifying cation-eutaxy structures, as its predictions were consistent with experimental observations and aligned with the results of the ion substitution algorithm. This model introduces a threshold for distinguishing cation-eutaxy from non-cation-eutaxy structures, which cannot be derived using the ion substitution algorithm alone. Materials falling outside this threshold can be excluded from consideration as cationic-eutaxy compounds, thus streamlining future high-throughput screening efforts. The complementary use of high-throughput screening and a phenomenological model offers a highly efficient pathway for identifying promising cation-eutaxy structures.

In this study, we developed a phenomenological model to classify cation-eutaxy compounds based on the following key parameters: electronegativity, ionic radius, and charge. Our model demonstrated a high level of agreement with the results obtained from high-throughput screening methods that combined ion substitution algorithms and formation energy calculations.

This study provides insight into the atomic properties that govern the formation of cation-eutaxy structures in ternary compounds, deepening our understanding of the fundamental factors influencing these arrangements. Moreover, the approach can be extended to other prototype structures suitable for topochemical etching, offering a framework for the synthesis of additional families of 2D materials. The current model is particularly effective for compounds composed of s- and p-block elements, but may be less applicable to systems involving transition metals, where d-orbital interactions and crystal field effects play a significant role in structural stability (Supplementary Fig. 12 and Supplementary Table 3). In addition, because the model does not explicitly incorporate thermodynamic parameters, it does not provide direct information on thermodynamic stability. Despite these limitations, the integration of atomic-scale understanding with predictive modeling contributes to the rational design of cation-eutaxy structures and their potential applications.

Methods

Computational screening for cation-eutaxy structures

To classify the III–V network dimensionality in ternary compounds, we employed a structure matching algorithm implemented in the Pymatgen Python library17. A dataset of 35 cation-eutaxy compounds was obtained from the Materials Project database29 and analyzed using the structure matching algorithm with a fractional length tolerance of 0.2 and site tolerance of 0.3, yielding 19 unique prototype structures.

To predict unknown cation-eutaxy structures, we employed an ionic substitution algorithm30 leveraging a data-driven statistical model trained on known compounds. This model assigns probabilities to element substitutions based on oxidation states and chemical environments, enabling targeted predictions of stable ternary compounds. The substitutions were validated based on a probability threshold of 1×10⁻⁸ and common oxidation states defined in the Smact library31. This methodology, originally trained on pnictide compounds, was implemented using Pymatgen’s substitution framework32, further enhancing the prediction accuracy of potential cation-eutaxy structures.

The structural refinement and validation of the predicted compounds were performed using density functional theory (DFT) calculations within the Vienna Ab initio Simulation Package (VASP) framework. Energy minimization and lattice relaxation were conducted using projector-augmented wave (PAW) potentials33,34,35, with the PBE exchange-correlation functional applied following the calculation settings specified in the Materials Project28.

The stability of each predicted cation-eutaxy structure was assessed via Ehull calculations using the PhaseDiagram API from Pymatgen32,36,37. A metastability tolerance of 200 meV/atom above the convex hull was applied, ensuring that only thermodynamically reasonable structures were considered.

All calculations were performed using publicly available open-source tools, ensuring reproducibility. Computational workflows, including structure-matching algorithms, ionic substitution models, and phase stability analyses, are fully documented and accessible via the Pymatgen and Smact libraries.

Descriptor formulation and weighting scheme

To quantitatively capture bonding characteristics in ternary A–III–V compounds, we introduced stoichiometry-weighted descriptors based on Econf and ionic radius. The stoichiometric ratios of constituent elements were systematically incorporated into the descriptor definitions to reflect composition-dependent variations in bonding behavior.

For the \({E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\) and \({r}^{{{\rm{w}}}}\), each elemental property was scaled by its stoichiometric coefficient normalized by the total number of atoms in the formula unit (\(l+m+n\)). For example, in a compound with the composition AlIIImVn, the weighted configuration energy and weighted ionic radius of cation A are defined as

$${E}_{{{\rm{conf}}}}^{{{\rm{w}}}}\left({{\rm{A}}}\right)=\frac{{{\rm{l}}}}{\left(l+m+n\right)}{E}_{{{\rm{conf}}}}\left({{\rm{A}}}\right),\,{r}_{{{\rm{A}}}}^{{{\rm{w}}}}=\frac{{{\rm{l}}}}{\left(l+m+n\right)}{r}_{{{\rm{A}}}}$$

Based on these values, \({E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}\) for A–V and III–V bonds were calculated as follows:

$${E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}({{\rm{A}}}-{{\rm{V}}})=\frac{\,l{\times E}_{{{\rm{conf}}}}\left({{\rm{A}}}\right)+{n\times E}_{{{\rm{conf}}}}\left({{\rm{V}}}\right)}{\left(l+m+n\right)},\\ {E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}({{\rm{III}}}-{{\rm{V}}})=\frac{\,m{\times E}_{{{\rm{conf}}}}\left({{\rm{III}}}\right)+{n\times E}_{{{\rm{conf}}}}\left({{\rm{V}}}\right)}{\left(l+m+n\right)}$$

In contrast, \({E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\) required a different scheme to resolve issues such as negative descriptor values and poor resolution across compositions. To address this, we incorporated the product of relevant stoichiometric coefficients and normalized by the sum of squared coefficients (\({l}^{2}+{m}^{2}+{n}^{2}\)). For instance, \({E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\) for A–V and III–V bonds were expressed as:

$${E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\left({{\rm{A}}}-{{\rm{V}}}\right)=\frac{\,n\times l{\times }({E}_{{{\rm{conf}}}}\left({{\rm{V}}}\right)-{E}_{{{\rm{conf}}}}\left({{\rm{A}}}\right))}{\left({l}^{2}+{m}^{2}+{n}^{2}\right)},\\ {E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}\left({{\rm{III}}}-{{\rm{V}}}\right)=\frac{\,n\times m{\times }({E}_{{{\rm{conf}}}}\left({{\rm{V}}}\right)-{E}_{{{\rm{conf}}}}\left({{\rm{III}}}\right))}{\left({l}^{2}+{m}^{2}+{n}^{2}\right)}$$

Final bonding descriptors were then defined as relative metrics:

$$\alpha=\frac{{E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}({{\rm{III}}}{-}{{\rm{V}}})}{{E}_{{{\rm{conf}}},{{\rm{avg}}}}^{{{\rm{w}}}}({{\rm{A}}}{-}{{\rm{V}}})},\,\beta=\frac{{E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}({{\rm{III}}}{-}{{\rm{V}}})}{{E}_{{{\rm{conf}}},{{\rm{diff}}}}^{{{\rm{w}}}}({{\rm{A}}}{-}{{\rm{V}}})}$$

This combined weighting strategy enabled consistent prediction and classification of III–V network dimensionality across a broad range of ternary compositions, by capturing subtle variations in bonding asymmetry and charge distribution that are not accessible using unweighted descriptors.

Preparation of unknown K2In2P3

K2In2P3 was synthesized using a solid-state method in an inert atmosphere to prevent the alkali metals from reacting with moisture and oxygen. A mixture of 0.195 g of bulk potassium (Thermo Fisher Scientific, chunks in mineral oil, 98%), 0.574 g of indium shot (Thermo Fisher Scientific, teardrops (4 mm), 99.99%), and 0.232 g of red phosphorus lumps (Thermo Fisher Scientific, lumps, 99.999 + %) was placed in an alumina crucible, which was then sealed inside a quartz ampoule. The ampoule was heated in a programmable furnace to 750 °C for 3 h, held at 750 °C for 20 h, cooled to 500 °C in 50 h, maintained at 500 °C for 100 h, and finally cooled to room temperature in 100 h. After the reaction, the quartz ampoule was broken under argon atmosphere to retrieve the crystals, which were then washed with methanol to remove the byproducts.

Preparation of unknown Na2In2As3

Na2In2As3 was synthesized by reacting stoichiometric amounts of bulk sodium (Thermo Fisher Scientific, ingot, 99.8%), indium shots (Thermo Fisher Scientific, teardrop (4 mm), 99.99%), and arsenic powder (Thermo Fisher Scientific, lump, 99.999%) in a 2:2:3 molar ratio. The materials were placed in an alumina crucible and subsequently loaded into a quartz ampoule inside an Ar-filled glove box. The ampoule was sealed under an argon atmosphere using the double-ampoule technique to prevent oxidation, with the smaller ampoule inside the larger one. The sealed ampoule was heated in a box furnace at a rate of 100 °C/h to 1000 °C and held at this temperature for 4 h. Afterward, the temperature decreased from 1000 to 500 °C at a rate of 5 °C/h, and the sample was maintained at 500 °C for 100 h to facilitate crystal growth. The ampoule was then allowed to cool naturally in a furnace. The resulting Na2In2As3 crystals exhibited a powdery morphology, with a typical particle size of approximately 0.1 mm and a gray color.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.