Predicting synthesizable multi-functional edge reconstructions in two-dimensional transition metal dichalcogenides

Two-dimensional (2D) transition metal dichalcogenides (TMDCs) have attracted tremendous interest as functional materials due to their exceptionally diverse and tunable properties, especially in their edges. In addition to the conventional armchair and zigzag edges common to hexagonal 2D materials, more complex edge reconstructions can be realized through careful control over the synthesis conditions. However, the whole family of synthesizable, reconstructed edges remains poorly studied. Here, we develop a computational approach integrating ensemble-generation, force-relaxation, and electronic-structure calculations to systematically and efficiently discover additional reconstructed edges and screen their functional properties. Using MoS2 as a model system, we screened hundreds of edge-reconstruction to discover over 160 reconstructed edges to be more stable than the conventional ones. More excitingly, we discovered nine new synthesizable reconstructred edges with record thermodynamic stability, in addition to successfully reproducing three recently synthesized edges. We also find our predicted reconstructed edges to have multi-functional properties—they show near optimal hydrogen evolution activity over the conventional edges, ideal for catalyzing hydrogen-evolution reaction (HER) and also exhibit half-metallicity with a broad variation in magnetic moments, making them uniquely suitable for nanospintronic applications. Our work reveals the existence of a wide family of synthesizable, reconstructed edges in 2D TMDCs and opens a new materials-by-design paradigm of ‘intrinsic’ edge engineering multifunctionality in 2D materials.


INTRODUCTION
Two-dimensional (2D) transition metal dichalcogenides (TMDCs) have been shown to exhibit exceptional electronic, magnetic, optical, mechanical, and catalytic properties [1][2][3][4][5] . The discontinuity represented by the one-dimensional (1D) edges of 2D TMDCs are analogous to the surfaces of bulk materials, which can impart new properties as well as functionalities [6][7][8][9] . For example, 2D MoS 2 monolayer is a direct gap semiconductor with a bandgap of 1.8 eV, while the 1D MoS 2 zigzag (ZZ) edge is metallic 10,11 . Moreover, it is the edges of 2D TMDCs that are responsible for high catalytic activities, while the basal planes are chemically inert [12][13][14][15] . Therefore, understanding the properties as well as the synthesis conditions of the edges in 2D TMDCs is critical to the design of new functional materials.
The lower symmetry of TMDCs compared with graphene, which normally exhibits ZZ or armchair (AC) edges 16 , enables them to present a greater diversity of edge structures and compositions. For example, using molecular beam epitaxy (MBE) under conditions of high Mo flux and atomic-resolution electron microscopy, two types of new reconstructed Mo-terminated edges have been synthesized and imaged. Attributed to the presence of the reconstructed edges, nanoporous MoS 2 films were found to have excellent activity for hydrogen evolution reaction (HER) and are magnetic up to 400 K 17 . Moreover, recent in situ heating experiments using scanning transmission electron microscopy (STEM) to track the edge evolution in monolayer Mo 1−x W x Se 2 (x = 0.05) flakes demonstrated that by varying the local chemical environment, a variety of reconstructed nonstoichiometric edges can be formed 18 . Using first-principles density functional theory (DFT), we have shown that these reconstructed MoSe 2 edges can have near optimal HER activity over conventional stoichiometric edges 19 . These studies demonstrated that the edges of TMDCs can be tuned at the nanoscale level, in some cases with high fidelity, giving rise to reconstructions that provide improved functional properties.
Despite these aforementioned successes, the whole family of synthesizable reconstructed edges in 2D TMDCs remains largely unknown. This presents a unique opportunity to computationally screen and discover additional functional reconstructed TMDC edges; in this article, we lay out such a computational approach. Starting with configuration ensemble generations, we screen for stable edges using a computationally affordable force-field method. The obtained stable edges are then further refined with DFT based electronic-structure calculations to generate phase diagrams and screen for their functional properties. Our study thus provides a comprehensive yet affordable computational scheme for predicting synthesizable functional edges of 2D materials and provides useful guidelines to experimental researchers. Many studies have investigated tuning the catalytic, electronic and magnetic properties of TMDCs by external doping in an Edisonian approach 20-23 , but the merit of our study is to discover a family of 'intrinsically' (based on the metal/chalcogen ratio) tunable materials with widely varying functional properties. Success of our study opens a new materials-by-design paradigm to search and discover other families of 2D 'intrinsic' edge reconstructions with specific multi-functionalities, with potential edge polytypism, for a wide range of nanoscale applications.

RESULTS AND DISCUSSION
Workflow for predicting synthesizable functional edge reconstructions in 2D TMDCs We use a nanoribbon model to study the edge reconstructions in 2D TMDCs. As shown in Fig. 1, for a given type (MoS 2 , MoSe 2 , WS 2 , WSe 2 …) and phase (2H, 1T, 1T′…) of TMDC, we set a region of width 'w' to look for all possible edge reconstructions in a nanoribbon of length 'L'. The total width of the nanoribbon in our study is seven pairs of Mo and S rows (~18.96 Å). We then generate n m edge configurations by tiling the underlying triangular lattice, where m and n correspond to the number of rows in the reconstructed region (i.e., m = w) and number of possible occupancies considered for each row, respectively. For example, in our work with MoS 2 , we set the width of the reconstructed region to be 4 rows, and consider 5 possible occupancies (blank, Mo, S up, S down, and S up and down) for each row, to obtain 5 4 (625) edge configurations. After the full set of edge configurations are obtained, we use a cost-effective forcefield method which has been benchmarked and validated against DFT to screen for stable edges. To date, several force-field methods have been developed for TMDCs [24][25][26] , and here we choose a recently reported reactive force-field (ReaxFF) model for MoS 2 27 . In principle, the structures can be also used to fit a neuralnetwork based surrogate model, as has been shown for other materials 28 . The stable edges are then further refined with explicit DFT calculations to generate phase diagrams and provide guidelines to experimentalists for synthesis. At last, the synthesizable edges are chosen to study the functionalities, including their catalytic activities (e.g., HER), half-metallicity and magnetic variations, with potentially more properties waiting to be discovered, such as topologically protected (magnetic) edgestate for dissipationless quantum nanotransport. In the following sections, we'll apply this workflow to a model system MoS 2 .
Screening for stable reconstructed edges To screen for stable reconstructed edges, we would like to use force fields or surrogate (neural-network-based) methods, as geometry optimizations (and simulated annealing if necessary) of thousands of edges at the DFT level of theory would be computationally prohibitive. Herein we choose a ReaxFF reactive potential, which can adequately describe the thermodynamic and structural properties of MoS 2 sheets, guided by extensive DFT simulations 27 . To validate the reliability of ReaxFF in identifying stable edges, we compare the formation energies from ReaxFF with DFT for all the generated 625 edge configurations without geometry optimizations, on both 2H and 1T MoS 2 phases. The convex hull plots (formation enthalpy as a function of Mo percent) are shown in Fig. 2a, b, d and e. Both DFT and ReaxFF calculations show that there exists a large quantity of edge configurations which have negative formation enthalpies. Some were found to be even more stable than the conventional unreconstructed edges (as indicated by the cyan dots in the figures), as they fall below the lines that connect the pure elements and the unreconstructed edges. Here, the unreconstructed edge refers to the ZZMo edge. As shown in Fig. 2c and f, with the exception of data points with unreasonably high formation enthalpies (indicating unstable edges with unphysical atomic geometries), there is a strong linear correlation between the DFT and ReaxFF energies, indicating that the ReaxFF potential can identify stable reconstructed edges.
Next, we use the ReaxFF potential to perform geometry optimizations of the generated 625 edge configurations to screen for stable edges. Figure 3 shows the convex hull plots after geometry optimizations with ReaxFF for 2H and 1T MoS 2 . One can see that most of the edges become energetically favorable (i.e., have negative formation enthalpies) after geometry optimizations. In addition, 161 2H and 168 1T MoS 2 edges fall below the lines connecting the unreconstructed edges and the pure elements. This indicates a significant portion of reconstructed edges for each phase are more stable than the conventional unreconstructed edges. The left and right slopes of each data point can give an idea about the stability range as a function of chemical potential, but we'll further generate phase diagrams to provide more accurate information on the synthesis conditions where these edges can be experimentally realized.
Screening for synthesizable reconstructed edges The lowest data points of each stoichiometry shown in Fig. 3 correspond to the edges which are most likely synthesized under specific chemical environments. Thus we choose these edge structures for 2H MoS 2 to perform further DFT calculations. The Workflow for predicting synthesizable multi-functional edge reconstructions in 2D TMDCs. Starting with configuration ensemble generations for a given size of the TMDC nanoribbon in a particular phase, stable edges are screened by using a computationally affordable force field method. The energies of stable edges are then further refined with DFT based electronic-structure calculations to generate phase diagrams to identify synthesizable edge reconstructions, whose functional properties can subsequently be computed.
reason why we choose the phase of 2H in this work is that it is more stable than 1T, and its known edges are active for catalysis [12][13][14][15] . Figure 4 shows the DFT refined edge structures for the lowest data points of each stoichiometry. The edges are labeled based on the orientation with respect to the hexagonal lattice of MoS 2 (e.g., Mo-oriented ZZMo), and the outermost termination group (e.g., S-terminated (-S) and Mo-terminated (-Mo)). Since our work considers the reconstructions on one side of the nanoribbons, the edges are all derived from the Mooriented unreconstructed edge (ZZMo). We have successfully reproduced 3 Mo-oriented edges found from previous experimental studies: ZZMo, ZZMo-S, and ZZMo-S2 17,18 , which validates our approach. These edges currently serve as the basis for most studies of physical and chemical properties of MoS 2 edges. More excitingly, we have also predicted nine new reconstructed edges with record stability not found in the literature (highlighted in orange). This discovery reveals that there is a large undiscovered set of experimentally realizable edges in the structure space which is being currently overlooked.
Notably, two of the predicted reconstructed edges contain grain boundaries (GBs): ZZMo-GB-Mo and ZZMo-GB-S. GBs in 2D materials connect two regions with different orientations. Various GBs have been observed in MoS 2 monolayers, such as the ones made of 4 | 4, 4 | 6, 4 | 8, 5 | 7, or 6 | 8 ring pairs [29][30][31] . GBs can affect the properties of MoS 2 monolayers significantly. For example, recent study shows that GBs in MoS 2 exhibit superior HER catalytic activity, combined with long-term stability 32 . The new predicted ZZMo-GB-Mo edge contains octagon rings, and this type of GB is firstly discovered. Meanwhile, ZZMo-GB-S contains a previously observed 6 | 8 boundary, which switches the orientation of the edge from Mo-oriented to S-oriented.
To better understand the synthesis conditions of the predicted reconstructed edges, we generated the phase diagram with the formation energy as a function of the sulfur chemical potential    5). As one can see, the most stable edges from the Mo-rich limit to S-rich limit are ZZMo-Mo2, ZZMo-MoS2, ZZMo-S, ZZMo-S2, and ZZMoS4. Therefore, these edges can be synthesized by tuning the chemical potential via controlled MBE growth 17 or via controlling heating and electron beam irradiation under STEM 18 . Our results on ZZMo-S and ZZMo-S2 are in agreement with previous studies 18 . Besides the five edges, the other seven edges can also be formed, as they have extremely close formation energies to the five edges under various chemical potential windows. For example, at the sulfur chemical potential of −3.9 eV, the formation energies of ZZMo-S2, ZZMo-S3, and ZZMo-S4 are within 0.1 eV/Å, and could likely coexist. In addition, it is reported that the ZZMo edge (indicated by the red dash line in Fig. 5) can be formed under Mo-rich conditions 17 . ZZMo-Mo2S, ZZMo-MoS, and ZZMo-GB-Mo have lower formation energies than ZZMo under Mo-rich conditions, so it is anticipated that they can also be formed under the same Mo-rich conditions.
For comparison, we have also plotted the often-discussed Soriented edge with 50% S coverage 33 in Fig. 5. It is a Mo-rich edge with the same stoichiometry as the ZZMo-MoS edge, so it is parallel to ZZMo-MoS in the figure. Compared with the Mo-rich Mo-oriented edges, the S-orientated edge with 50% S coverage is more stable than ZZMo-Mo2S, ZZMo-MoS, ZZMo-GB-Mo, but is less stable than ZZMo-Mo when µ S is from −7.0 to −6.4 eV. However, it becomes the most stable edge among all the studied edges in this work when µ S increases to the range of −6.4 to −6.1 eV. As µ S further increases, the stoichiometric ZZMo-MoS2 becomes the most stable edge, followed by three S-rich Moorientated edges ZZMo-S, ZZMo-S2, and ZZMo-S4.
One point of interest to note is that ZZMo is not the most stable stoichiometric edge (33.3% Mo). The most stable stoichiometric edge we discovered is the reconstructed ZZMo-MoS2 edge, which is 0.63 eV/supercell lower in energy than the unreconstructed ZZMo edge, and more importantly, 0.35 eV/supercell lower than the currently recognized most stable stoichiometric edge 34 . This energetically strongly preferred stoichiometric ZZMo-MoS2 edge should be considered in future studies of MoS 2 edges. It may have interesting consequences in growth. For example, if growth continues from this edge, it could propagate along a different angle, and then a GB which bends the overall sheet structure can be obtained. The metastable edges like ZZMo might be favored due to kinetic stabilization, and they could also be more catalytically active 35 . Thus we select three metastable edges to study to simply demonstrate the viability of the metastable reconstructions. Beyond these, there are hundreds of other metastable edges which are worthy of investigation in the future. Supplementary Figure 1 shows the structures of four metastable edges. ZZMo-S2Mo is another metastable stoichiometric edge, and contains seven-member rings, which often appear in GBs of MoS 2 monolayers. ZZMo-S and ZZMo-GB contain 31.5% Mo (same  stoichiometry as ZZMo-S). ZZMo-S has similar configuration as the observed Mo-Klein edge 17 , while ZZMo-GB contains a 6 | 8 boundary similar as ZZMo-GB-S, but with less S at the edge. It is expected that these metastable edges can be kinetically stabilized and show promising catalytic activities.
One thing to note is that µ Mo and µ S need to remain smaller than µ Mo (bulk) and µ S (bulk), respectively, to avoid the decomposition of MoS 2 to bulk Mo or S. Our computed µ S (bulk) ≈ E S (bulk) = −4.2 eV, and µ Mo (bulk) ≈ E Mo (bulk) = −11.5 eV. By assuming µ Mo + 2µ S = ε at equilibrium, where ε is the energy per formula unit of MoS 2 monolayer, the upper and lower bounds of µ S were determined to be −4.2 and −5.4 eV(as indicated by the two vertical dashed lines in Fig. 5). This is in good agreement with previous studies 36 . Although some of the edges are predicted to be stabilized under conditions beyond this limit, these edges have been observed experimentally due to some kinetic reasons 17,18 .
HER activity on predicted synthesizable reconstructed edges One major promising application for the predicted synthesizable reconstructed edges are for HER, part of the water splitting reaction [37][38][39] . Electrochemical HER using renewable energy can provide a sustainable supply of fuel for future societies with hydrogen as the key energy carrier 40 . Since H adsorption is the first step in HER, and the Gibbs free energy for hydrogen adsorption (ΔG H ) is a well-established descriptor for the HER 41-43 , we screen for HER activity via ΔG H . Figure 6a shows the computed ΔG H for the fifteen predicted edges (including both stable and metastable edges and the metastable edges are marked with asterisk). Since a close-to-zero value of ΔG H suggests a good HER catalyst, we predict six edges (ZZMo-S, ZZMo, ZZMo-MoS2, ZZMo-S, ZZMo-Mo2, and ZZMo-GB-Mo) to be good HER catalysts (the absolute values of ΔG H for these edges are within 0.2 eV of the optimal point). Our results on the existing edges (ZZMo, ZZMo-S, ZZMo-S2) are in agreement with previous studies 14,17 . Notably, we discovered four new synthesizable reconstructed edges could be better HER catalysts than the existing active edges. Figure 6b shows the structures of H adsorption at the six most active edges. For the S-terminated edges (ZZMo-S, ZZMo-MoS2, and ZZMo-S), H is adsorbed at the S sites, while for the Mo-terminated edges (ZZMo, ZZMo-Mo2, and ZZMo-GB-Mo), H is adsorbed at the Mo sites. Thus, H always prefers to adsorb at the under-coordinated sites, as these sites are more active 44 . We have considered H coverage effects on ΔG H for all the edges, and the optimal H coverages (when ΔG H is most close-to-zero) of the six most active edges are 50%, 100%, 50%, 25%, 25%, and 50%, respectively. H coverage is defined by n H /n site , where n H is the number of adsorbed H, and n site is total number of equivalent adsorption site.
Our previous study shows that the optimal HER performance cannot be obtained in the pristine ZZMo and ZZS (or Se) edges 19 . Indeed, here we further prove that the reconstructed edges could achieve better HER performance. Specifically, we find that ΔG H is 0.10 eV for ZZMo, while it decreased to 0.03 eV for ZZMo-MoS2, and −0.01 eV for ZZMo-S. Since ZZMo-MoS2 is the most stable stoichiometric edge discovered so far, it is highly expected that this predicted reconstructed edge can be synthesized and applied to HER in the near future. While the energetic difference is not dramatic, a tiny change in ΔG H near zero could result in a difference of several orders of magnitude improvement in exchange current density 45 . One thing to note is that our calculations assume that the catalyst is charge neutral for computational simplicity, but in reality, the catalyst is usually charged to match its Fermi level with the applied electrode potential. Previous study has shown that due to the charge induced change in the occupation of electronic states, the charge on 2D materials can have a strong impact on the electrochemical reactions 43 . Therefore, while our current approach looks at the intrinsic HER activities, we hope to further explore HER activities in operating conditions on the reconstructed edges by including the charge effects in the future. Moreover, the chemical interactions between the reaction intermediate adsorbed on the catalyst and the surrounding solvent molecules are neglected in our study. These interactions can be strong for reactions in aqueous solution when the intermediate contains a highly charged atom that can form hydrogen bond(s) with water, for example electrochemical CO 2 reduction (CO 2 R) 46 . Since the adsorbed intermediate H at the studied edges is nearly charge neutral, it is expected that these interactions won't be very strong, which we also hope to further explore in the future. Table 1 lists the electronic and magnetic properties of the nanoribbons containing the predicted synthesizable reconstructed edges. It shows that edge reconstructions can effectively modify the properties of the nanoribbons. Specifically speaking, nanoribbons with reconstructed edges could become halfmetallic. Figure 7 shows the density of states (DOS) of the halfmetallic nanoribbon with the unreconstructed ZZS edge on one side and the reconstructed ZZMo-MoS3 edge on the other side. The spin up channel has a gap of 0.38 eV, while the spin down channel is metallic. Our spin density and projected DOS plots show that the half-metallicity is mainly from the reconstructed side. The reconstruction-induced peculiar half-metallic nature which allows electrons of only one spin direction to move can be used in spin-based electronics as proposed for half-metallic graphene nanoribbons 47,48 . Meanwhile, the magnetic moments of the nanoribbons vary significantly due to the edge reconstructions. It was recently reported that by integrating multiple magnetic units into a single piece of MoS 2 nanoribbon, magnetic devices for information storage can be achieved 49 . Here we provide many synthesizable reconstructed edges with different local magnetic moments, which can further improve the density of information storage.

Further applications in nanodevices
To further understand potential applications of reconstructed edges in nanodevices, we examine the electronic and magnetic properties of MoS 2 nanoribbons with all the generated 625 edge configurations. Interestingly, we find that 41 2H MoS 2 nanoribbons become half-metallic due to the edge reconstruction, representing a significant proportion of the generated structures. Supplementary Fig. 2 shows the DOS for one of the half-metallic nanoribbons (the geometry is also shown). One can see that near the Fermi level, one spin channel is continuous, while other spin channel has a relatively wide gap of~0.44 eV. Meanwhile for the 1T MoS 2 nanoribbons, we find they are all metallic, consistent with the fact that 1T MoS 2 monolayer itself is metallic. Regarding their magnetic properties, Fig. 8 shows the histograms of the magnetic moments for the 625 2H and 1T MoS 2 nanoribbons. The pristine ZZ MoS 2 nanoribbons possess magnetic moment of~1.1 µ B /supercell, in agreement with previous studies 11 . However, with edge reconstructions the magnetic moment can be dramatically changed. The range of the magnetic moment for 2H and 1T MoS 2 is 0.3~6.2 µ B /supercell, and 0.0~4.8 µ B /supercell, respectively. Most of the 2H nanoribbons are in the range of 0.5~1.5 µB/supercell, while most of the 1T nanoribbons are in 0.0~1.0 µ B /supercell. Therefore, the magnetic moment of a nanoribbon can be increased by 4~6 times by edge reconstructions. These results show that reconstructed edges have great potential applications in nanodevices due to their tunable electronic and magnetic properties.
Given the potential for nanoribbons to host topologically protected states, in the future, it is worthy to investigate how edge reconstructions modify the topology of the electronic-bands and related properties of MoS 2 nanoribbons as recently reported for the carbon nanoribbons 50 . It is known that pristine MoS 2 nanoribbons demonstrate the quantum valley Hall effect, while some other TMDC nanoribbons present a topological insulator phase 51 . Moreover, 1D topological superconductivity can be realized by placing a TMDC monolayer with M terminated ZZ edge on top of a superconducting substrate 52 . From these revelations one is led to conceive edge reconstructions would affect the topology. Given that topology also plays a role in catalysis [53][54][55] , there are naturally interesting unexplained implications for (say) HER performance. Applying machine-learning algorithms on our workflow-generated dataset to screen the functional properties of more on-lattice and off-lattice reconstructions, even with larger wavelength reconstructions than we could compute using current DFT-approaches, is another interesting direction.

DISCUSSION
In conclusion, we have developed a computational approach to rapidly and efficiently discover more synthesizable functional edges in the family of 2D TMDCs. Using MoS 2 as an example, we screened 625 edge configurations for 2H and 1T MoS 2 phases, and predicted stable edges to guide the experimental synthesis. Subsequently we studied the functional properties of these edges. We discover many of these intrinsically tunable edge reconstructions to be near-optimal for HER. More excitingly, they exhibit halfmetallicity and broad variation in local magnetic moments (0.3~6.2 µ B /supercell) making them suitable in nanodevices for applications in spintronics and information storage. In the future, the effects of edge reconstructions on the topological properties of TMDC nanoribbons are worthy of further investigations. Combining our current computational approaches to other neural-network and/or evolutionalry structure-search algorithms can lead to accelerating discoveries of more functional edgereconstructions in 2D materials. Our work thus opens a new materials-by-design paradigm, to 'intrinsically' edge engineer multi-functionalities in 2D TMDCs for a wide range of technological applications.

Details of DFT and force field simulations
Spin-polarized DFT calculations were performed using the Vienna ab initio simulation package (VASP) 56,57 . Electron exchange-correlation was represented by the functional of Perdew, Burke and Ernzerhof (PBE) of generalized gradient approximation (GGA) 58 . The ion-electron interaction was described with the projector augmented wave (PAW) method 59 . The plane-wave cutoff is set to 400 eV, and a conjugate gradient method is

Formation energy calculations
The formation energy per atom was calculated by where E tot (Mo x S y ) is the total energy of the Mo x S y nanoribbon, and μ Mo and μ S are the chemical potentials of Mo and S. To determine whether the edges can be formed and are stable against decomposition into bulk Mo or bulk S, we constructed convex hull plots in Figs. 2 and 3. The chemical potentials used there are µ Mo in body-centered cubic (bcc) Mo and µ S in orthorhombic alpha-S, which were obtained from the computed energies per atom at T = 0 K by neglecting the entropic terms for the condensed states (i.e., µ S (bulk) ≈ E S (bulk) = −4.2 eV, and µ Mo (bulk) ≈ E Mo (bulk) = −11.5 eV). To generate the phase diagram (formation energy as a function of chemical potential), it was assumed that µ Mo + 2µ S = ε, where ε is the energy per formula unit of MoS 2 monolayer 61 . Then, the formation energy changed to H adsorption energy calculations The differential H adsorption energy ΔE H was calculated by where E(Mo x S y + nH) and E[Mo x S y + (n − 1)H] represent the total energy of the Mo x S y nanoribbon with n and (n -1) adsorbed H atoms, respectively, and E(H 2 ) represents the total energy of a gas phase H 2 molecule. A negative value of ΔE H suggests favorable absorption. The Gibbs free energy of H adsorption ΔG H was obtained by where ΔE ZPE is the difference in zero-point energy between the adsorbed H and H in the gas phase H 2 molecule, and ΔS H is the entropy difference between the adsorbed H and ½ H 2 in the gas phase under standard conditions. The zero-point energy was calculated by summing vibrational frequencies ω over normal modes ν: E ZPE ¼ 1 2 P hω. The entropy of the free H 2 molecule at 298.15 K and 1 atm was taken from the NIST database.

DATA AVAILABILITY
Data reported in this paper is available from the corresponding author upon reasonable request.

CODE AVAILABILITY
Code used to generate the edge configurations in this paper is available from the corresponding author upon reasonable request.