## I. Type-I Penrose P3 system

## Target system

One of the main objectives of this study is to confirm that the E-VAE model can be used as an efficient computational approach to search for energetically stabilized spin states on a complicated spin-ice system. For this, we implement a Penrose P3 dipolar spin-ice system which is referred to as a Type-I system in this study (Fig. 2A).

This system has a quasi-crystalline pattern composed of seven different types of unit vertices [25]; Each magnetic moment is centered on the side of two types of rhombuses composed of Penrose tiles, and its direction is strongly constrained to be parallel with the side. Translational symmetry is not available on this system [24], thus the difficulty of searching for the true ground state increases exponentially with the size of the system. Nevertheless, Shi et al. have proposed a method to build up the ground state for this system through elaborated logical steps [25]. (See Fig. S1 in Supplementary Information, SI.) They found that the proposed ground state consists of two different spin groups, the skeleton and flippable parts. The skeleton part is a quasi-one-dimensional rigid part, and it has a unique (up to time-reversal symmetry) ground state with long-range ordering. The flippable part is composed of topologically induced emergent frustrated vertices or spins [25, 29]. These are usually surrounded by the skeleton part, and they can be flipped to other degenerate spin configurations with identical energy. To demonstrate the utility of E-VAE model in searching for the ground states of spin-ice systems, we set this Type-I system as our first target system and investigate whether the E-VAE model can generate the proposed ground state properly.

As shown in the sub-figure of Fig. 2A, the dipole-dipole interaction model is considered only between the nearest neighbors as indicated by Eq. (1),

\(\epsilon =\frac{E}{ND}=-\frac{1}{N}\sum _{⟨ij⟩}\frac{3\left({\widehat{m}}_{i}\bullet {\overrightarrow{r}}_{ij}\right)\left({\widehat{m}}_{j}\bullet {\overrightarrow{r}}_{ij}\right)-{\widehat{m}}_{i}\bullet {\widehat{m}}_{j}{\left|{\overrightarrow{r}}_{ij}\right|}^{2}}{{\left|{\overrightarrow{r}}_{ij}\right|}^{5}},\) | (1) |

where \(\epsilon\) is a unitless energy density parameter, \(E\) is the total energy of the system, \(D\) is the dipole-dipole interaction strength written in energy unit, \({\widehat{m}}_{i}\) is a magnetic moment located at i-th spin site, and \({\overrightarrow{r}}_{ij}\) is the unitless displacement vector between the dipole moments \({\widehat{m}}_{i}\) and \({\widehat{m}}_{j}\). The \(⟨ij⟩\) index pair under the summation refers to the nearest neighbor pairs. This short-range interaction scheme is adopted to exactly follow the logical steps that Shi et al. used to obtain the ground state by tiling vertices.

## Dataset generation

To secure the dataset composed of numerous metastable spin states, we independently perform multiple simulated annealing processes on the system and gather the final local energy minimum states. The total dataset is composed of 50,000 spin states, and it is divided into three sub-datasets: training, validation, and test datasets with 30,000, 10,000, and 10,000 data, respectively. The details about the simulated annealing process and data generation are given in the Methods section.

A sample of the simulated annealing results is shown in Fig. 2B. Here, we first show that simulated annealing is not an appropriate method to search for the ground state of complex spin-ice systems. To demonstrate this clearly, we separated the spins in this sample into two groups which are composed of the spins located on the skeleton (Fig. 2C) and flippable (Fig. 2D) parts of the proposed ground state (Fig. S1A in SI). It is clearly shown that the two colors (red and blue) form several domains, and the domains are intermixed in most regions in Fig. 2C. Considering that the magnetic moments in the skeleton part in the proposed ground state make a specific long-range ordering, the intermixed domains indicate that the simulated annealing method which is one of the representative conventional optimization methods is not enough to investigate the true ground state due to the enormous complexity and the existence of numerous metastable states of this system. In addition, there are lots of vertices in their excited states as shown in Fig. 2D, and it clearly indicates that the spin state shown in Fig. 2B is not in the lowest energy state.

We investigated the energy histogram for the 30,000 spin states in our training dataset (simulated annealing results) as shown in Fig. 2E. The \({\epsilon }_{G}\), which is the energy value calculated using the Eq. (1) and the spin configuration of the proposed ground state by Shi et al., is significantly lower than the energy values in the histogram. This fact qualitatively supports that the simulated annealing method has limitations for solving this complex spin-ice system.

## Training results

We train our E-VAE model using the metastable spin states in the training dataset and investigate the energy distributions of generated states from the trained E-VAE models for each \(\gamma\) value as shown in Fig. 2F. As \(\gamma\) increases, the energy distribution shifts to a lower energy region. Considering the Hamiltonian loss, \({L}_{\mathcal{H}}\), is related to the energy of the generated state and \(\gamma\) is the controlling parameter for the independent variation of the magnitude of\({L}_{\mathcal{H}}\), the behavior of energy distribution shown in Fig. 2F indicates that the \({L}_{\mathcal{H}}\) is properly minimized during the training process of E-VAE model. (The behaviors of each loss term during the training process are given in Fig. S2 and Note 1 in SI) In addition, the energy distribution collapses into a sharp peak when \(\gamma\) exceeds a certain value; in Fig. 2F, it is approximately \(3.0 \tilde 4.0\). In other words, the collapsed E-VAE model generates a sufficiently lower energy spin state in the Type-I system compared with the simulated annealing results. Hence, we confirm that the E-VAE model can be utilized as an innovative numerical method significantly outperforming conventional methods in estimating the ground state of the complex spin-ice system.

The generated spin state from the collapsed E-VAE model trained with \(\gamma =5.0\) is shown in Fig. 2G. Comparing with the skeleton and flippable parts of the simulated annealing result shown in Fig. 2C and D, there is a perfectly ordered skeleton structure with the 5-fold rotational symmetry which is implied in the frame structure of the system, and all vertices in the flippable part are in one of their degenerate ground state configurations. Thus, the state found by the E-VAE model is one of the possible configurations of the spin state proposed by Shi et al., clearly supporting that it is a candidate for the true lowest energy state of the Type-I system.

In addition, the E-VAE reveals a surprising fact: the spins highlighted as the black arrows in Fig. 2G at several specific five-fold vertices, though known as representative flippable vertices [25], should actually be classified as the skeleton part, and not the flippable part. One can notice that there is an anti-clockwise flow formed by the black and red spins, as indicated by the black dotted circle in Fig. 2G, which can be found repeatedly throughout the spin states at the specific five-fold vortex sites. If the five-fold vertices are truly independent flippable vertices, then the anti-clockwise and clockwise flows should appear randomly. Considering the spin state shown in Fig. 2G is a significantly low energy spin state generated by the trained E-VAE model, we suspect that the spins highlighted as the black arrows are not flippable but rigid. In fact, all possible configurations are checked for the rigidity of these sites (see Fig. S3 in SI), and we confirm that the spins highlighted as the black arrows cannot be flipped until the entire skeleton structure is inversed. Consequently, we propose a new skeleton-flippable configuration including the newly revealed additional skeleton parts as shown in Fig. 2G.

## II. Type-II Penrose P3 system

## Target system

To search for the ground states of unexplored spin-ice systems, we use a different type of aperiodic Penrose P3 pattern to implement complicated spin-ice systems which are referred to as Type-II systems in this study. We construct three Type-II systems including 640, 1195, and 2150 spins; Fig. 3A shows the frame structure of a Type-II system with 640 spins. Note that the frame structure of the system is clearly distinct from the Type-I system as shown in the sub-figures in Fig. 3A. The ground states of Type-II systems have never been investigated before to the best of our knowledge. The short-range interaction scheme shown in Eq. (1) is also considered in Type-II systems to generate datasets and to define the Hamiltonian losses for each system.

## Proposing the ground state

We first apply the E-VAE model to the Type-II system with \(N=640\). The same network structure shown in Fig. 1 is also used to implement the E-VAE model for this case. We generate datasets and train the E-VAE model using the same strategy as for the Type-I system.

After the training process, we obtain a spin state from the trained E-VAE model. The energy value of the spin state is significantly lower than the values in the energy distribution of the simulated-annealing-generated spin states (Fig. 3B). Figure 3C shows the spin state obtained from the trained E-VAE model, and we propose it as a candidate of the ground state of the Type-II system. Similar to the Type-I system, the Type-II system also shows a robust skeleton part, although it has a different quasi one-dimensional ordering. We also confirm the existence of the flippable part: each vertex can be replaced by a different spin configuration with identical energy. Although only applied to two different spin-ice systems, we believe our E-VAE model has the potential to be utilized to estimate the ground state of various complex spin systems.

## System size dependency

As discussed earlier, the difficulty of determining the ground state of the Penrose P3 spin-ice systems using conventional methods increases exponentially with the system size, owing to the system being composed of aperiodic patterns with no translational symmetry in the spin network structures. On the contrary, we speculate that the difficulty of searching for the ground state using the E-VAE model may not be greatly affected by the size of the system owing to an efficient grouping process of the encoder network structure.

The E-VAE model is designed to include an encoder and decoder structure as shown in Fig. 1. The encoder network compresses the given input data into a single latent code. Through this encoding process, it is expected that several essential features from the vast amount of information implied in the input data are extracted and encoded into the latent code. Therefore, each component in the code represents collective characteristics (usually called high-level features in the deep learning field), combining several simple characteristics of the input data. In the case of the spin-ice systems considered in this study, each component in the latent codes can be connected to various collective states composed of multiple spins. In other words, the encoder part in E-VAE will play the role of reducing the total degree of freedom of the raw data through an efficient grouping process.

In order to validate this claim, we apply E-VAE model to two other Type-II systems with different \(N\) (\(N=1195\) and \(2150\)). The spin states found by the trained E-VAE models for each of these \(N\) cases are shown in Fig. S4 in SI, while the ground state energies are shown in Fig. 4.

For all cases, the trained E-VAE model discovers a significantly lower energy state than the energy distribution of the simulated-annealing-generated dataset used to train each model. Considering that standard deviation value of each distribution, \({\sigma }_{\text{S}\text{A}}\), represents the energy difference between the spin states that can be reasonably obtained through the simulated annealing process, it is surprising that the energy value of each spin state obtained from the trained E-VAE model is several times of \({\sigma }_{\text{S}\text{A}}\) away from the mean value of each distribution, \({\mu }_{\text{S}\text{A}}\).

To quantitatively show the remarkable performance of our method, we investigate the \(\varDelta \epsilon /{\sigma }_{\text{S}\text{A}}\) which is a measure of the energy difference between \({\epsilon }_{\text{E}-\text{V}\text{A}\text{E}}\) and \({\mu }_{\text{S}\text{A}}\) (\(\varDelta \epsilon /{\sigma }_{\text{S}\text{A}}\) graph in Fig. 4). Note that \(\varDelta \epsilon /{\sigma }_{\text{S}\text{A}}\) increases with the system size (\(N\)). For instance, in the case of \(N=2150\), \(\varDelta \epsilon\) is as large as \(13{\sigma }_{\text{S}\text{A}}\). Following empirical rule, the probability (\(P\)) of obtaining a spin state with the energy value of \({\epsilon }_{\text{E}-\text{V}\text{A}\text{E}}\) through the simulated annealing process is infinitesimal (\({P\approx 10}^{-38}\) for \(13{\sigma }_{\text{S}\text{A}}\)). Therefore, we believe that our E-VAE method significantly outperforms the simulated annealing method especially as the system size increases, demonstrating that this method transcends the limitations of conventional methods for the problem of searching the ground state of complicated spin-ice systems.