Elastic energy driven multivariant selection in
martensites via quantum annealing

Lara C P. dos Santos, Tian Hang, Roland Sandt, Martin Finsterbusch, Yann

Le Bouar, Robert Spatschek

► To cite this version:

Lara C P. dos Santos, Tian Hang, Roland Sandt, Martin Finsterbusch, Yann Le Bouar, et al.. Elastic energy driven multivariant selection in martensites via quantum annealing. Physical Review Research, 2024, 6 (2), pp.023076. 10.1103/physrevresearch.6.023076. hal-04650576

HAL Id: hal-04650576

https://hal.science/hal-04650576v1

Submitted on 16 Jul 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Elastic energy driven multivariant selection in martensites via quantum annealing

Lara C. P. dos Santos,¹ Tian Hang,¹ Roland Sandt,¹ Martin Finsterbusch,² Yann Le Bouar,³ and Robert Spatschek ¹‚⁴

¹Structure and Function of Materials, Institute of Energy and Climate Research IEK-2, Forschungszentrum Jülich GmbH,
52425 Jülich, Germany

²Materials Synthesis and Processing, Institute of Energy and Climate Research IEK-1, Forschungszentrum Jülich GmbH,
52425 Jülich, Germany

³Université Paris-Saclay, ONERA, CNRS, Laboratoire d'Etude des Microstructures, 92320 Châtillon, France

⁴JARA-ENERGY, 52425 Jülich, Germany

<figure><img src="image_2.png" /></figure>

(Received 30 January 2024; accepted 15 March 2024; published 22 April 2024)

We demonstrate the use of quantum annealing for the selection of multiple martensite variants in a microstructure with long-range coherency stresses and external mechanical load. The general approach is illustrated for martensites with four different variants based on the minimization of the linear elastic energy. The equilibrium variant distribution is then analyzed under application of tensile and shear strains and for different values of the considered shear and tetragonal contributions of the different martensite variants. The interface orientations between different domains of variants can be explained using the perspective of the elastic energy anisotropy for regular stripe patterns. For random grain orientations, the response to an external elastic strain is weaker and variant changes can be interpreted based on the rotated eigenstrain tensor.

DOI: 10.1103/PhysRevResearch.6.023076

I. INTRODUCTION

The formation of microstructures is critical for many applications, as it strongly impacts the properties of materials. In many cases, simulation approaches can nowadays be used to support experimental investigations related to the understanding and prediction of microstructures. Among these approaches, the phase field approach plays a central role, as it is able to predict the temporal evolution of nonequilibrium microstructures [1–5]. To this end, continuous order parameter fields are used to discriminate between different phases, grains, and, in particular, martensite variants, and in this way phase transformation kinetics can be simulated. Contrary to this approach, we recently presented an approach which complements the phase field picture in the sense that it separates discrete from continuous degrees of freedom [6]. Specific applications are stress- or strain-driven transformations in shape memory alloys or martensitic transformations in general. There, we assumed for simplicity that a grain can be in one out of two martensite variant states, which differ by their eigenstrain, i.e., their stress-free configuration. By application of external stresses and the presence of internal stresses, the grains can flip between the two variants, and these transitions are driven by a reduction of the overall elastic energy.

The prediction of martensite formation and variant selection has been widely discussed in the literature, see,

e.g., Ref. [7]. In this context, the concept of a lattice correspondence between austenite and martensite gained attention, resulting, e.g., in the well-known Kurdjumov-Sachs or Nishiyama-Wasserman orientation relationships. From a micromechanical perspective, the strain energy minimization led to an important development, which is the phenomenological theory of martensite crystallography [8,9], which is to a large extent in agreement with experimental observations. These descriptions are different but related to global energetic considerations, which use the overall energy of the system to predict the formation of microstructures.

In the present paper, we focus on the long-ranged elastic effects, which result from mismatches between different martensite variants. They lead to elastic misfits, which induce interactions not only between neighboring phases and grains but also depend on and influence the entire microstructure. As discussed in Ref. [6], the proper minimization of the elastic energy requires keeping all these interactions, as an artificial cutoff leads to improper predictions of variant selections. However, this global perspective introduces an additional degree of complexity to the overall problem, as now a huge space of microstructural configurations needs to be considered. For example, for a microstructure consisting of *N* grains and two martensite variants, we therefore have altogether $2^N$ combinations, and the energy needs to be determined for each of them to find the global ground state, which corresponds to the (low temperature) thermodynamic equilibrium. It is obvious that this configuration space is already too large for moderate grain size numbers, and therefore not only a brute force energy minimization but also the use of heuristic approaches like simulated annealing are limited and require large amounts of computing times. If we consider, e.g., a system consisting of just *N* = 20 grains, the number of

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
configurations is $2^N \simeq O(10^6)$; for a typical FFT calculation in 2D with sufficient accuracy, the single core runtime on current processors is of the order of one minute. This implies that scanning all configurations to determine the global energetic minimum requires about 17 000 core hours of computing time. The question, which arises also generally in the context of materials science, is how this optimization problem can be solved efficiently via quantum computing.

A general quantum computer with a sufficient number of qubits is not available yet, but quantum annealing (QA) technology [10–14] with several thousands of qubits and couplers is available. In contrast to general-purpose quantum computers, only specific classes of problems can be handled via QA, i.e., binary quadratic optimization problems [15]. The basic concept of QA is the initialization of its qubits in well-defined Hamiltonians, whose ground state is unique and known [16]. This Hamiltonian is then changed adiabatically under the operation at cryogenic temperatures to the desired, final problem Hamiltonian, ideally staying in the ground state throughout this transformation [16,17] and therefore enabling efficient global energy minimizations. The application of QA in materials science is still rare, however, a few publications have focused on sampling techniques via QA [18–24], phase transitions in the transverse field Ising model [25], critical phenomena in frustrated magnetic systems [26], energy calculations of defective graphene structures [27], and the automated design of metamaterials [28]. Instead, the focus of actual research concerning QA rather lies on performance tests and benchmarking of QA against conventional approaches [29–32]. Among them is the discovery that QA can be superior to conventional simulated annealing, as quantum (multiqubit) tunneling can be a more efficient path to the ground state than thermal excitations [33,34].

In our previous publication [6], we succeeded in mapping the martensite problem to a quantum annealer formulation, i.e., a binary quadratic model expressed through an Ising model or, equivalently, in terms of a quadratic unconstrained binary optimization [15], which allows us to find the true ground-state configurations even for many thousand grains in an extremely short time [6]. This approach builds upon earlier works in Refs. [35–38]. Although this conceptual progress demonstrates well the benefits of QA for materials science modeling, the chosen example was limited, in particular, by the restriction to just two variants, whereas, e.g., typical shape memory alloys or martensites in steel possess many more variants [39–42]. In the present paper, we demonstrate how the approach can be generalized to multivariant martensites and show how the resulting variant distributions are influenced by different tensile strains and random grain rotations.

## II. METHODS

The overall concept of the investigations is that we assume that in a microstructure, as depicted, e.g., in Fig. 2, each (martensite) grain is allowed to select between different variants, such that the overall elastic energy is minimized in equilibrium. As the grains are coherently connected to each other, a variant flip as sketched in Fig. 1 leads possibly to an energetically unfavorable distortion of the surrounding grains,

<figure><img src="image_3.png" /><figcaption>FIG. 1. Visualization of the variants used in the present paper. Spin $s_1$ corresponds to the shear deformation, $s_2$ to the tetragonal deformation, see Eq. (14).</figcaption></figure>

and therefore the discrete optimization of the variants is a nontrivial problem.

To achieve a description in the framework of continuum elasticity, we consider eigenstrain (stress-free strain) expressions of the form

$$ \epsilon_{ij}^{(0)}(\mathbf{r}) = \sum_{n=1}^{N} \theta_n(\mathbf{r}) \sum_{k=1}^{K} s_v \epsilon_{ij}^{(0,n,k)}. \quad (1) $$

Here, the function $\theta_n(\mathbf{r})$ is equal to one inside grain $n$ and vanishes everywhere else. The eigenstrain in each grain with index $n$ can be a superposition using $K$ spin variables $s_v = s_{n,k} = \pm 1$ with the total number of spins $N_v = N \cdot K$, which weight eigenstrain tensor components $\epsilon_{ij}^{(0,n,k)}$. Therefore, altogether we have $2^K$ martensite variants in each grain by suitable selection of the parameters $\epsilon_{ij}^{(0,n,k)}$.

To determine the thermodynamic equilibrium state (minimum elastic energy), first the elastic problem has to be solved for a given eigenstrain. This can be expressed through the minimization of the elastic energy with respect to the displacement field $u_i$ in combination with suitable mechanical boundary conditions. In the case of isotropic linear elasticity with phase and variant independent Lamé coefficient and shear modulus, the elastic energy functional reads

$$ E_{el} = \int_V \left( \frac{1}{2} \lambda (\epsilon_{kk} - \epsilon_{kk}^0)^2 + \mu (\epsilon_{ik} - \epsilon_{ik}^0)^2 \right) dV, \quad (2) $$

with the system volume $V$ and the strain tensor $\epsilon_{ik} = \frac{1}{2}(\partial_i u_k + \partial_k u_i)$. In the above expression, the Einstein sum convention is used.

Contrary to the optimization of the displacement and strain fields as continuous variables, the variant selection becomes a discrete problem on top of the preceding step. Following the Fourier transformation approach outlined in Ref. [6], we can calculate the elastic energy $E(\{s_v\})$ of the microstructure for a given spin configuration $\{s_v\}$. We note that the approach does not necessarily require us to use Fourier transformation methods but, e.g., finite element methods can also be used,
provided that the elastic interactions can be computed with sufficient accuracy, as discussed in Ref. [6].

As the expression for the elastic energy is quadratic in the framework of the linear theory of elasticity, the final expression for the elastic energy can be expressed as

$$E(\{s_v\}) = E_0 + \sum_{v=1}^{N_v} \tilde{H}_v s_v + \sum_{v, \eta=1}^{N_v} \tilde{J}_{v\eta} s_v s_\eta, \quad (3)$$

with a quadratic, symmetric matrix $\tilde{J}_{v\eta}$. The coefficients $E_0$, $\tilde{H}_v$, and $\tilde{J}_{v\eta}$ can be expressed directly in terms of the Fourier transformation solution or, alternatively, be calculated using any elastic solver, which delivers the total elastic energy $E(\{s_v\})$ for a given spin configuration $\{s_v\}$. We note that the above form of the energy holds for arbitrary values of the spin variables and does not apply only to the special case $s_v = \pm 1$. Therefore, we directly obtain

$$E_0 = E(\{s_v = 0\}), \quad (4)$$

where all spins have the value 0, which we also denote as the austenite reference state. Next, we perform calculations where all spins but one vanish, and the selected one has either the value +1 or -1. From that, we get

$$E(\{s_i = \pm 1, \text{all others } 0\}) = E_0 \pm \tilde{H}_i + \tilde{J}_{ii}. \quad (5)$$

Therefore, we can extract the diagonal element for the self-interactions as

$$\tilde{J}_{ii} = \frac{1}{2}[E(\{s_i = +1, \text{all others } 0\}) + E(\{s_i = -1, \text{all others } 0\} - 2E_0)] \quad (6)$$

and, consequently,

$$\tilde{H}_i = \frac{1}{2}[E(\{s_i = +1, \text{all others } 0\}) - E(\{s_i = -1, \text{all others } 0\})]. \quad (7)$$

Finally, to calculate the interaction coefficients $\tilde{J}_{ij}$ for $i \neq j$, we perform elastic calculations with two nonvanishing spins $s_i = s_j = 1$. With

$$E(\{s_i = s_j = 1, \text{others } 0\}) = E_0 + \tilde{H}_i + \tilde{H}_j + 2\tilde{J}_{ij} + \tilde{J}_{ii} + \tilde{J}_{jj}, \quad (8)$$

where we can directly calculate $\tilde{J}_{ij}$ using the expressions obtained above,

$$2\tilde{J}_{ij} = E(\{s_i = s_j = +1, \text{all others } 0\}) - E(\{s_i = +1, \text{all others } 0\}) - E(\{s_j = +1, \text{all others } 0\}) + E_0. \quad (9)$$

For zero stress boundary conditions, we get for the eigenstrain of type (1) the simplifications $E_0 = 0$ and $\tilde{H}_i = 0$ for all $i$.

For the use of the quantum annealer, we restrict the allowed spin values to $s_v = \pm 1$ and we also have to rewrite the expression (3) in a slightly different way, which excludes, in particular, self-interactions of the type $s_i s_i$. The classical Hamiltonian reads

$$H = H_0 + \sum_{i=1}^{N_v} H_i s_i + \sum_{i=1}^{N_v} \sum_{j=i+1}^{N_v} J_{ij} s_i s_j, \quad (10)$$

where the offset $H_0$ only affects the energy value but not the spin configuration minimizing the energy, and therefore it can be dropped for the desired Ising type of the QA problem Hamiltonian. Physically, it corresponds to the elastic energy of an austenitic microstructure as reference state; see Fig. 1. Comparison with the expression (3) above gives

$$H_0 = E_0 + \sum_{i=1}^{N_v} \tilde{J}_{ii}, \quad (11)$$

$$H_i = \tilde{H}_i, \quad (12)$$

$$J_{ij} = 2\tilde{J}_{ij} \quad \text{for } i < j. \quad (13)$$

This general formulation has the advantage that it decouples the approach for setting up the Ising coefficients for the annealer from the specific linear elastic energy solver. Nevertheless, one has to keep in mind that the formulation (3) relies on the assumption of equal elastic constants in all phases and variants and coherent boundary conditions. However, it is valid also in situation beyond the original work [6], where we used only a single spin to discriminate between the two martensite variants in each grain. Here, the same approach also applies for multiple variants using several spins per grain.

On the quantum annealer level, the determined Ising coefficients are used to define an energy landscape, where superconducting loops with clockwise or anticlockwise circulating currents define qubits with different spin states [16]. These superconducting loops interact with external fluxases, which allow us to control the energy difference and barrier height of the constructed energy landscape [16]. The starting point of the annealing process is the initialization of the system in the ground state of a known Hamiltonian $H_D \sim -\sum_i \sigma_i^x$, where $\sigma_i$ denotes the Pauli matrices, i.e., corresponding to a strong transverse magnetic field [17,43]. This Hamiltonian is turned during the annealing into the desired Ising model (10) for which an energetic minimum is sought, $\min_{\{s_i=\pm 1\}} H_p$ [15]. The Hamiltonians $H_D$ and $H_p$ do not commute [15], and the time taken by the initial Hamiltonian to reach the low-energy state is sufficiently large to establish the validity of the adiabatic theorem of quantum mechanics [44], which postulates that a system stays in its eigenstate during adiabatic changes. In contrast to classical approaches, further quantum mechanical principles like tunneling to leave metastable regions and the usage of entangled states inside QA processors (QPU) are employed [45]. Especially, if energetically close low-energy states exist, the QA process does not always determine the true ground state, therefore a suitable number of repetitions is performed and the lowest detected energy is chosen. Additionally, for larger systems, hybrid QA utilizes classical algorithms and the interplay with QA to address areas with high computational demands using a QPU coprocessor working with generic parameters for up to 11 616 spin variables on the D-Wave Advantage system [46–48].

### III. RESULTS

To illustrate the workflow, we consider a case with four variants per grain, i.e., two spins per grain. Similar to the previous work [6], we use specific cases of shear and tetragonal
distortions, which are now combined and lead to eigenstrains

$$ \epsilon_{ij}^{(0,n)} = s_{n,1} \begin{pmatrix} 0 & \epsilon_0 & 0 \\ \epsilon_0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} + s_{n,2} \begin{pmatrix} \epsilon_1 & 0 & 0 \\ 0 & -\epsilon_1 & 0 \\ 0 & 0 & \epsilon_1 \end{pmatrix}, \quad (14) $$

where n numerates the grains. For the first application, we assume that the individual tensors of the shear (involving the strain parameter $\epsilon_0$) and tetragonal deformations (associated with $\epsilon_1$) are the same in all grains; later we will also discuss the case of grain rotations, where the tensors differ from grain to grain. Figure 1 visualizes the different strain variants and illustrates the color codings used for the following plots. The grain structure is generated using a Voronoi tesselation in two dimensions, and we assume plane strain conditions, using isotropic linear theory of elasticity with periodic boundary conditions. According to the above approach, first the self-energies related to each spin are calculated, followed by a computation of spin-spin interactions, from which the Ising vector components $H_n$ and interaction parameters $J_{ij}$ are extracted.

Apart from the self-generated internal stresses due to the different variants, a homogeneous external (average) strain field $\langle \epsilon_{ij} \rangle$ can also be superimposed, which leads to an additional contribution to the magnetic field vector $H_n$, which couples the external strain $\langle \epsilon_{ij} \rangle$ to the $\mathbf{k} = 0$ mode of the Fourier transform of the eigenstrain tensor. This means that a change of the external strain leaves the coupling constants $J_{ij}$ invariant and only modifies a contribution to $H_n$, which can be treated analytically, therefore minimizing the computational cost for repeated calculations with different external boundary conditions. Figure 2 shows a sequence of equilibrium microstructures and their response to a tensile strain in the horizontal ($x$) direction for $\epsilon_1 = \epsilon_0$. Here, the computer-generated, random 2D microstructure consists of 200 grains. Each grain selects its optimum variant state ($s_{n,1}, s_{n,2}$) to minimize the total energy. In Fig. 2(a), the fixed external strain vanishes, which implies that in the ground state the volume fractions of the blue and green (respectively, red and grey) variants should be close to each other, as then the positive and negative contributions to the shear and tetragonal (eigen)strain cancel each other to minimize the elastic energy. The arrangement of the differently colored patches is then a result of the mutual spin-spin interactions. We can observe in Fig. 2(a) a structure of horizontal stripes with alternating layers of shear strain (blue-grey and red-green). Such stripes with (100) orientation have already been discussed in Ref. [6], as the formation of these structures minimizes the overall shear deformation. Additionally, tilted stripe pairs red-blue and green-grey correspond to the same tetragonal strain and annihilating shear strain, similar to the patterns for pure tetragonal eigenstrain discussed in Ref. [6].

The patterns formed by several variants can be better understood by analyzing the elastic anisotropy of the transformation. Specifically, the normal to the interface between two variants must be close to a direction that minimizes the elastic kernel of the transformation from one variant to the neighboring one. This can be demonstrated for platelike domains [49] but is also very commonly observed in complex microstructures whose anisotropy is of elastic origin, such as cuboidal

<figure><img src="image_4.png" /><figcaption>FIG. 2. Equilibrium variant distribution with uniform grain orientation. The microstructure consists of 200 grains and a tensile strain is applied in horizontal ($x$) direction. The tensile strain is (a) $\langle \epsilon_{xx} \rangle / \epsilon_1 = 0$, (b) $\langle \epsilon_{xx} \rangle / \epsilon_1 = 0.25$, (c) $\langle \epsilon_{xx} \rangle / \epsilon_1 = 0.5$, (d) $\langle \epsilon_{xx} \rangle / \epsilon_1 = 0.75$, (e) $\langle \epsilon_{xx} \rangle / \epsilon_1 = 1.0$, and (f) $\langle \epsilon_{xx} \rangle / \epsilon_1 = 1.5$. In all figures, $\epsilon_0 = \epsilon_1$ and the Poisson ratio is chosen as $\nu = 1/4$ (i.e., $\lambda = \mu$).</figcaption></figure>

microstructures in superalloys [50], Widmanstätten structures [51], or chessboard structures in Co-Pt alloys [52]. The elastic kernel can be obtained analytically in Fourier space [49] or can be calculated in real space, using a regular stripe pattern of both variants. Indeed, in this geometry, only Fourier modes perpendicular to the stripes remain, and the elastic energy is proportional to the elastic kernel. In the following, the orientation dependence of the pattern is analyzed by considering successively each pair of variants. We assume fixed vanishing average strain and equal volume fractions of the two involved variants. The angular dependence of relevant combinations of stripe pairs is shown in Fig. 3. In these pole figures, the angle 0° corresponds to vertical stripes, and positive angles to anticlockwise rotations of the pattern. From this figure, one can conclude that both green-grey and red-blue interfaces are preferred either at 0° or 90°, and this is in agreement with the observed interface orientation between these variant pairs in Fig. 2(a).
<figure><img src="image_5.png" /><figcaption>FIG. 3. Orientation dependence of the normalized elastic energy for equal volume fraction stripes consisting of pairs of variants for $\langle\epsilon_{ij}\rangle = 0$ and $\epsilon_1 = \epsilon_0$.</figcaption></figure>

Both green-red and grey-blue pairs have preferred orientations around $\pm 38^\circ$. Due to constraints of the chosen system sizes and the used periodic boundary conditions, typically deviations from the analytical prediction for infinite system sizes are expected, as discussed in detail in Ref. [6].

In the equilibrium microstructure in Fig. 2(a), we indeed find inclined interfaces between green-red and grey-blue pairs, which are in line with the theoretical expectation.

Finally, green-blue and red-grey interfaces show a very pronounced anisotropy with a strong cusp at either $17^\circ$ and
<figure><img src="image_6.png" /><figcaption>FIG. 4. 2 × 2 supercell representation of Fig. 2 to emphasize the microstructure patterns. The dashed yellow lines mark the boundaries or the simulation with periodic boundary conditions. In all figures, $\epsilon_1 = \epsilon_0$ and $\nu = 1/4$.</figcaption></figure>

−62° (green-blue) or 62° and −17° (red-grey). Due to this strong orientation dependence, these interfaces either have to appear at these preferred angles or they are strongly suppressed. In fact, the corresponding interface lengths in Fig. 2(a) are much shorter than the other interfaces.

We note for the subsequent panels of Fig. 2 that the elaborated anisotropy of the interactions is not affected by a homogeneous external strain, and therefore the polar plots in Fig. 3 also remain valid in these cases. In these pictures, the applied tensile external strain in the horizontal direction is incrementally increased. This favors the variants with a positive tetragonal spin value, and consequently the fraction of green and grey grains increases. They form more and more pronounced bands in the (100) directions, as explained above, partly interrupted by disappearing inclined red-blue bands with unfavorable tetragonal strain contribution. Finally, only grains with $s_{n,2} = +1$ remain, and the volume fractions with positive and negative shear are essentially the same, as much as the irregular grain sizes permit.

It is useful to emphasize the periodicity of the patterns by looking at a 2 × 2 supercell of the simulation domain, as shown in Fig. 4. The first striking observation is that the pattern in Fig. 4(c) for $\langle\epsilon_{xx}\rangle/\epsilon_0 = 0.5$ differs significantly from the other cases. These latter cases show both horizontal and tilted stripes, as analyzed above. There, grey-green and red-blue interfaces appear at 0° or 90°, and green-blue as well

as red-grey interfaces play at most a minor role, as discussed above. The red-green and grey-blue interfaces are roughly at ±45° as a compromise between the preferred analytical predictions for infinite systems and the constraints by the periodic boundary conditions.

The special case $\langle\epsilon_{xx}\rangle/\epsilon_0 = 0.5$ shows a remarkable topological change, as now interfaces between green-blue, green-red, grey-blue, and grey-red appear, whereas the previously important red-blue and green-grey interfaces are absent. Moreover, the newly appearing interfaces are at their preferred orientations. For example, the marked 24° angle inside the red lozenge is the difference between the 62° minimum energy orientation for the red-grey interface and the 38° orientation for the red-green pattern. This leads to the geometry sketched in Fig. 4(c), from which the overall periodic pattern is constructed. According to the equilibrium angles, this demands that the joint volume fraction of the diamond-shaped green and grey domains is $\sin(80°)/[\sin(24°) + \sin(80°)] \approx 0.7$, and therefore this exceptional pattern is only observed close to a specific external strain value. We note that this chessboard structure is very similar to the one observed in Co-Pt and (CuAu)_{1-x}Pt_x systems [52]. However, the symmetry is lower in the present case because the green and grey regions are not tetragonal.

The same microstructure as in Fig. 2 is subjected to a shear strain $\langle\epsilon_{xy}\rangle$ in the (110) direction, and the results of
<figure><img src="image_7.png" /><figcaption>FIG. 5. Equilibrium variant distribution with uniform grain orientation. The microstructure consist of 200 grains and an average shear strain $\langle\epsilon_{xy}\rangle$ is applied. The strain is (a) $\langle\epsilon_{xy}\rangle/\epsilon_0 = 0$, (b) $\langle\epsilon_{xy}\rangle/\epsilon_0 = 0.25$, (c) $\langle\epsilon_{xy}\rangle/\epsilon_0 = 0.5$, (d) $\langle\epsilon_{xy}\rangle/\epsilon_0 = 0.75$, (e) $\langle\epsilon_{xy}\rangle/\epsilon_0 = 1.0$, and (f) $\langle\epsilon_{xy}\rangle/\epsilon_0 = 1.5$. In all figures, $\epsilon_0 = \epsilon_1$, and the chosen Poisson ratio is $\nu = 1/4$ (i.e., $\lambda = \mu$).</figcaption></figure>

the energy minimization are shown in Fig. 5. The starting configuration without external strain Fig. 5(a)] is therefore identical to Fig. 2(a). An increase of the shear strain favors variants with positive $s_{n,1}$, as they accommodate the given external strain to a large amount. Therefore, for large positive shear strain $\langle\epsilon_{xy}\rangle/\epsilon_0$ values, only grey and blue patches remain. As they differ by the tetragonal eigenstrain only, they align as inclined bands to minimize the total elastic energy, as discussed in detail in Ref. [6]. In particular, grey-green and red-blue interfaces are at 0° or 90°, as predicted in Fig. 3. Similarly, blue-green and red-grey interfaces have orientations close to the expected ones (±17°, ±62°). For blue-grey and red-green interfaces, the periodicity constraint enforces orientations close to 45°.

An interesting case is the three-variant configuration in Fig. 5(d) (apart from a negligible fraction of the red variant). A 2 × 2 supercell of the equilibrium microstructure is shown in Fig. 6. The overall 45° alignment is imposed by periodic boundary conditions. The green domain slightly changes the orientation of the blue-grey interfaces so that they are rotated closer to the optimal orientation 90° – 38° = 52°. This configuration would not be favorable at –45° (instead of the 45°) because the blue-grey interfaces would be rotated away from the optimal orientation. Instead, for –45°, red domains appear in the blue platelets [see Fig. 5(e)].

For the present four-variant setup, the response of the microstructure does not only depend on the ratio $\langle\epsilon_{ij}\rangle/\epsilon_0$ but

<figure><img src="image_8.png" /><figcaption>FIG. 6. 2 × 2 supercell of the simulation domain for external shear $\langle\epsilon_{xy}\rangle/\epsilon_0 = 0.75$ from Fig. 5(d) with $\epsilon_1 = \epsilon_0$. Apart from a negligible amount of the red variant, three phases are present here. The sketch shows the geometrical explanation of the equilibrium pattern.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. Equilibrium variant distribution with uniform grain orientation for $\langle\epsilon_{ij}\rangle = 0$. The microstructure consists of 200 grains and the dimensionless parameter $\epsilon_1/\epsilon_0$ is varied. The ratios are (a) $\epsilon_1/\epsilon_0 = 0.2$, (b) $\epsilon_1/\epsilon_0 = 1$, and (c) $\epsilon_1/\epsilon_0 = 5$. In all figures, the chosen Poisson ratio is $\nu = 1/4$ (i.e., $\lambda = \mu$).</figcaption></figure>

also on the relative relevance of shear and tetragonal deformations, as expressed through the dimensionless parameter $\epsilon_1/\epsilon_0$. Whereas in the above investigations we have fixed this ratio to $\epsilon_1/\epsilon_0 = 1$, Fig. 7 shows the influence of this parameter for fixed external strain. Here, we note that a change of $\epsilon_0 \rightarrow a\epsilon_0$ and $\epsilon_1 \rightarrow b\epsilon_1$ does not require us to repeat the entire computations for the interaction parameters $J_{ij}$ and the external magnetic field $H_i$. Instead, it follows from the theory of linear elasticity that the coupling constants $J_{ij}$ are quadratic in the rescaling factors $a, b$ and $H_i$ scales linearly with them. In detail, for a spin pair $(i, j)$, where both of them refer to shear transformations, the coupling constants are rescaled according to $J_{ij} \rightarrow a^2 J_{ij}$, and similarly for two tetragonal spins as $J_{ij} \rightarrow b^2 J_{ij}$. For mixed interactions between a tetragonal and shear degree of freedom, the rescaling obeys $J_{ij} \rightarrow abJ_{ij}$. The magnetic terms $H_n$ scale linearly with either $a$ or $b$, depending on the spin type.

The starting case (center panel in Fig. 7) is again the already discussed reference situation with $\epsilon_1 = \epsilon_0$. In all three pictures, the boundary conditions are fixed to vanishing average strain, $\langle\epsilon_{ij}\rangle = 0$. The right panel shows the equilibrium microstructure for $\epsilon_1/\epsilon_0 = 5$, which means that the weight of the tetragonal deformations is significantly larger than that of the shear. Here, the variants select essentially diagonal bands, as discussed in Ref. [6], in order to match the periodic boundary conditions. The blue and green variants have canceling tetragonal and shear eigenstrains, such that for equal volume fractions the zero strain boundary conditions are met best. We can also expect that the opposite color pair (red-grey) can lead to an energetically equivalent solution, and therefore a weak, locally noisy pattern, which mixes these two configurations, may appear. We expect this effect to be more pronounced in larger systems. One may expect that the opposite case with $\epsilon_1/\epsilon_0 = 0.2$, which is dominated by the shear eigenstrain, should similarly lead to bands with alternating shear in the $\langle 100 \rangle$ direction, as discussed for the two-variant case in Ref. [6]. However, as can be seen from the left panel of Fig. 7, this is not the case. Instead, each horizontal stripe consists of blue and grey (or red and green) patches, which share the same shear eigenstrain but have opposite tetragonal

strain. The reason for this asymmetry—as compared to the tetragonal dominated case discussed above—is that a small, but nonvanishing tetragonal eigenstrain also leads to a tangential mismatch strain at the interfaces between the $\langle 100 \rangle$ bands. Therefore, a purely blue horizontal band would lead to a contraction in the tangential direction, whereas the compensating green horizontal band would expand in the same direction, hence leading to large mismatch coherency stresses. Therefore, it is energetically favorable to generate $\langle 100 \rangle$ structures with an additional zigzag substructure, such that all four variants are present with equal volume fraction and only the irregularity of the microstructure adds some minor noise.

It is again instructive to look at the angular dependence of the elastic energy for infinite stripe patterns (again for zero mean strain boundary conditions), as it depends on the ratio $\epsilon_1/\epsilon_0$, see Fig. 8. Here, one can clearly see that the interaction becomes very anisotropic for green-grey interfaces at low ratios $\epsilon_1/\epsilon_0$ and essentially isotropic at high values, with minimum energetic cost at 0° and 90°. For the green-red interfaces, the trend is opposite. Green-red (and similarly grey-blue) interfaces have preferred interface orientations $\pm 38^\circ$. The green-blue and similarly grey-red interfaces have the strongest dependence on $\epsilon_1/\epsilon_0$ and, therefore, the overall energy minimization clearly favors the optimal orientations for these interfaces. From this perspective, Fig. 7(a) is particularly interesting, and its zigzag structure becomes obvious in the $2 \times 2$ supercell representation in Fig. 9. The white lines emphasize the structure and illustrate that for this microstructure all interface orientations are optimal and not affected by periodicity constraints. The last example concerns a case where the grains are oriented randomly, i.e., rotated around the $[001]$ axis, as depicted in the left panel of Fig. 10. Like in the previous publication [6], the response to a uniaxial strain $\langle\epsilon_{xx}\rangle$ is less pronounced than in the cases with uniform grain orientation above. The center panel shows the equilibrium microstructure for vanishing mean strain, whereas the right panel is for a rather high tensile strain $\langle\epsilon_{xx}\rangle/\epsilon_1 = 2.5$.

For a clockwise rotation of a grain by an angle $\alpha$ around the $[001]$ axis, the relevant components of the eigenstrain tensor
<figure><img src="image_10.png" /><figcaption>FIG. 8. Polar representation of the orientation dependence of the elastic energy (scaled to the maximum value) for mismatching pairs of variants. The calculations are for fixed strain boundary conditions $\langle \epsilon_{ij} \rangle = 0$.</figcaption></figure>

become

$$ \epsilon_{xx}^{(0)'} = \cos^2 \alpha \epsilon_{xx}^{(0)} + 2 \cos \alpha \sin \alpha \epsilon_{xy}^{(0)} + \sin^2 \alpha \epsilon_{yy}^{(0)}, \quad (15) $$

$$ \epsilon_{xy}^{(0)'} = -\sin \alpha \cos \alpha \epsilon_{xx}^{(0)} + (\cos^2 \alpha - \sin^2 \alpha) \epsilon_{xy}^{(0)} + \sin \alpha \cos \alpha \epsilon_{yy}^{(0)}, \quad (16) $$

$$ \epsilon_{yy}^{(0)'} = \sin^2 \alpha \epsilon_{xx}^{(0)} - 2 \cos \alpha \sin \alpha \epsilon_{xy}^{(0)} + \cos^2 \alpha \epsilon_{yy}^{(0)}. \quad (17) $$

Therefore, for the chosen example, a rotation by $\pi/4$ converts the tetragonal deformation to a shear transformation (and vice versa), apart from the eigenstrain contribution in the $z$ direction, which only affects the energetic cost but not the strain distribution. Therefore, we expect for a tensile load in the

[100] direction, that grains with orientation close to 0 or $\pi$ (the response is $\pi$-periodic), have a tetragonal spin preferentially being in state +1; for a rotation near 90°, the tetragonal contribution prefers the $s = -1$ state. For both of these orientations (and nearby values), the shear contributions are indifferent, as they have only a vanishing or small contribution $\epsilon_{xx}^{(0)'}$, which is sensitive to the tensile strain. As explained above, for rotations by approximately $\pi/4$ (or $3\pi/4$), the role of shear and tetragonal deformations is exchanged, but the original tetragonal transformation has an energy penalty due to the contribution in the $z$ direction for the used plane strain setup.

For some selected grains in Fig. 10, the change of the variant states is discussed in the following. Grain 1 has an
<figure><img src="image_11.png" /><figcaption>FIG. 9. 2 × 2 supercell of the periodic solution Fig. 7(a) with $\epsilon_1/\epsilon_0 = 0.2$. The structure's geometry is sketched and highlights the interface orientations between the variants.</figcaption></figure>

orientation close to 0 (or $\pi$). In the zero strain state, it is marked in red, hence the shear spin equals –1 and the tetragonal component is also –1. Under stronger tensile load it turns grey, hence both spins flipped their sign, leading to a stronger alignment with the external field. Grain 2 has an orientation close to 90°, is originally red (shear: –1, tetragonal –1) and becomes blue in the right panel. In this case, the dominant effect is to have negative tetragonal spin values, as discussed above. Finally, for grain 3, which has a mixed orientation, we effectively have a transition from blue to grey, where the positive shear spin is relevant due to the orientation near $\pi/4$, whereas the flipping tetragonal spin component is less relevant.

IV. CONCLUSIONS

From a materials science perspective, quantum quantum is only starting to influence the opportunities and methodologies of modeling approaches. It can be expected that, in the future, general purpose quantum computing will have the potential to strongly accelerate materials science related simulations and to enable completely different possibilities. To our knowledge, sufficiently large quantum computers are not yet available, and also suitable algorithms have not yet been developed. Nevertheless, QA, as a specific type of adiabatic quantum computing, is already commercially available with large numbers of qubits today. The type of problems which can be addressed is limited, as it is required to map the materials science related problem of interest to the minimization of an Ising Hamiltonian.

In the present paper, we have demonstrated that the selection of multiple martensite variants can be mapped to an Ising problem for specific conditions and have demonstrated how the energy minimization can be performed efficiently via QA. The superiority of QA compared to classical computing has been demonstrated in our previous publication [6]. The multiplicity of the domains together with the long-range, anisotropic, and shape-dependant character of the elastic interactions make the determination of the ground state a particularly difficult problem, and therefore QA opens completely different possibilities. The central result of the present publication is the generalization to multivariant systems and the analysis of the arising microstructures. Also, we demonstrate how the proper coefficients of the Ising model can be obtained from general linear elastic calculations.

To illustrate the approach, we use a four variant per grain case, where each grain is represented by two Ising spins. We use a linear combination of shear and tetragonal deformations to demonstrate the methodology. As a result, the equilibrium microstructure depends on the applied external strains (or stresses), as well as on the strength of the distortions induced by the different variants. The observed patterns can largely be explained by computations of the anisotropy of lamellar arrangements of the different variants, showing that the

<figure><img src="image_12.png" /><figcaption>FIG. 10. Equilibrium variant distribution with random grain orientation. (a) Grain orientation map corresponding to the microstructure. The grain rotation angle in the color bar is given in radian (modulo $\pi$ due to symmetry). The rotation axis is along the [001] direction. The microstructure consists of 200 grains and tensile strain in horizontal direction is (b) $\langle\epsilon_{xx}\rangle/\epsilon_1 = 0$ and (c) $\langle\epsilon_{xx}\rangle/\epsilon_1 = 2.5$. In all panels, $\epsilon_0 = \epsilon_1$, and the chosen Poisson ratio is $\nu = 1/4$ (i.e., $\lambda = \mu$). The grains have a random orientation based on a uniform distribution.</figcaption></figure>
mismatch of neighboring domains dominates the overall elastic energy.

From a fundamental perspective, it is important to demonstrate whether QA is superior to classical global minimization strategies such as simulated annealing. The present example of long-range elastic interactions between misfitting martensite variants falls into the class of random spin glass models, and therefore it is plausible that—at least for small system sizes—the minimization directly benefits from multiqubit tunneling [33,34]. From the perspective of the investigated application, suitable scenarios to study this effect could be constructed through specific arrangements of grains such that the underlying energy landscape has controlled false minima, similar to Ref. [34]. However, from a materials science perspective, large numbers of grains are usually more representative and desirable, and therefore it is less likely that the tunneling advantage can be detected in a straightforward way in the investigated examples. Moreover, beyond the order of 100 grains, the use of hybrid approaches is currently required and used in the present paper, which changes the minimization dynamics. The results indicate that the ground-state microstructural configuration is not always found, therefore requiring repeated simulations, but, in general, the observed

structures and interfaces are in agreement with the expected true minimum energy interface orientations, as demonstrated in this paper.

For many applications in science, the formation of microstructures is critical and efficient simulation approaches to support experimental findings are highly desired. Our developed QA approach presents an efficient technological opportunity for the determination of equilibrium microstructures with long-range elastic interactions, where multiple martensite variants need to be considered.

## ACKNOWLEDGMENTS

This work was funded by the German Federal Ministry of Education and Research (BMBF) via the project ALANO (Grant No. 13XP0396B) and the Helmholtz project ZeD-aBase. Open access was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—No. 491111487. The authors gratefully acknowledge the Jülich Supercomputing Centre [53] for funding this project by providing computing time on the D-Wave Advantage System JUPSI through the Jülich UNified Infrastructure for Quantum computing (JUNIQ).

[1] A. Karma and W. J. Rappel, Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics, Phys. Rev. E **53**, R3017 (1996).

[2] A. Karma and W. J. Rappel, Quantitative phase-field modeling of dendritic growth in two and three dimensions, Phys. Rev. E **57**, 4323 (1998).

[3] G. Boussinot and E. A. Brener, Achieving realistic interface kinetics in phase-field models with a diffusional contrast, Phys. Rev. E **89**, 060402 (2014).

[4] K. Wang, G. Boussinot, C. Hüter, E. A. Brener, and R. Spatschek, Modeling of dendritic growth using a quantitative nondiagonal phase field model, Phys. Rev. Mater. **4**, 033802 (2020).

[5] A. Finel, Y. Le Bouar, B. Dabas, B. Appolaire, Y. Yamada, and T. Mohri, Sharp phase field method, Phys. Rev. Lett. **121**, 025501 (2018).

[6] R. Sandt, Y. Le Bouar, and R. Spatschek, Quantum annealing for microstructure equilibration with long-range elastic interactions, Sci. Rep. **13**, 6036 (2023).

[7] B. Yue, F. Hong, N. Hirao, R. Vasin, H.-R. Wenk, B. Chen, and H.-K. Mao, A simple variant selection in stress-driven martensitic transformation, Proc. Natl. Acad. Sci. **116**, 14905 (2019).

[8] M. Wechsler, D. Lieberman, and T. Read, On the theory of the formation of martensite, Trans. AIME **197**, 1503 (1953).

[9] J. S. Bowles and J. K. Mackenzie, The crystallography of martensite transformations I, Acta Metall. **2**, 129 (1954).

[10] A. B. Finnila, M. A. Gomez, C. Sebenik, C. Stenson, and J. D. Doll, Quantum annealing: A new method for minimizing multidimensional functions, Chem. Phys. Lett. **219**, 343 (1994).

[11] J. Brooke, D. Bitko, T. F. Rosenbaum, and G. Aeppli, Quantum annealing of a disordered magnet, Science **284**, 779 (1999).

[12] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse Ising model, Phys. Rev. E **58**, 5355 (1998).

[13] S. Morita and H. Nishimori, Mathematical foundation of quantum annealing, J. Math. Phys. **49**, 125210 (2008).

[14] A. Rajak, S. Suzuki, A. Dutta, and B. K. Chakrabarti, Quantum annealing: An overview, Philos. Trans. R. Soc. A **381**, 20210417 (2023).

[15] R. H. Warren, Mathematical methods for a quantum annealing computer, JAAM **3**, 82 (2018).

[16] M. W. Johnson, M. H. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, Quantum annealing with manufactured spins, Nature (London) **473**, 194 (2011).

[17] S. Boixo, T. F. Rønnow, S. V. Isakov, Z. Wang, D. Wecker, D. A. Lidar, J. M. Martinis, and M. Troyer, Evidence for quantum annealing with more than one hundred qubits, Nat. Phys. **10**, 218 (2014).

[18] R. Sandt and R. Spatschek, Efficient low temperature Monte Carlo sampling using quantum annealing, Sci. Rep. **13**, 6754 (2023).

[19] J. Nelson, M. Vuffray, A. Y. Lokhov, T. Albash, and C. Coffrin, High-quality thermal Gibbs sampling with quantum annealing hardware, Phys. Rev. Appl. **17**, 044046 (2022).

[20] M. Vuffray, C. Coffrin, Y. A. Kharkov, and A. Y. Lokhov, Programmable quantum annealers as noisy Gibbs samplers, PRX Quantum **3**, 020317 (2022).

[21] S. Mandrà, Z. Zhu, and H. G. Katzgraber, Exponentially biased ground-state sampling of quantum annealing machines with transverse-field driving Hamiltonians, Phys. Rev. Lett. **118**, 070502 (2017).
[22] M. S. Könz, G. Mazzola, A. J. Ochoa, H. G. Katzgraber, and M. Troyer, Uncertain fate of fair sampling in quantum annealing, Phys. Rev. A **100**, 030303(R) (2019).

[23] S. Mukherjee and B. K. Chakrabarti, On the question of ergodicity in quantum spin glass phase and its role in quantum annealing, J. Phys. Soc. Jpn. **88**, 061004 (2019).

[24] M. Yamamoto, M. Ohzeki, and K. Tanaka, Fair sampling by simulated annealing on quantum annealer, J. Phys. Soc. Jpn. **89**, 025002 (2020).

[25] R. Harris, Y. Sato, A. J. Berkley, M. Reis, F. Altomare, M. H. Amin, K. Boothby, P. Bunyk, C. Deng, C. Enderud, S. Huang, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, T. Lanting, R. Li, T. Medina, R. Molavi, R. Neufeld, T. Oh, I. Pavlov, I. Perminov, G. Poulin-Lamarre, C. Rich, A. Smirnov, L. Swenson, N. Tsai, M. Volkmann, J. Whittaker, and J. Yao, Phase transitions in a programmable quantum spin glass simulator, Science **361**, 162 (2018).

[26] P. Kairys, A. D. King, I. Ozfidan, K. Boothby, J. Raymond, A. Banerjee, and T. S. Humble, Simulating the Shastry-Sutherland Ising model using quantum annealing, PRX Quantum **1**, 020320 (2020).

[27] B. Camino, J. Buckeridge, P. A. Warburton, V. Kendon, and S. M. Woodley, Quantum computing and materials science: A practical guide to applying quantum annealing to the configurational analysis of materials, J. Appl. Phys. **133**, 221102 (2023).

[28] K. Kitai, J. Guo, S. Ju, S. Tanaka, K. Tsuda, J. Shiomi, and R. Tamura, Designing metamaterials with quantum annealing and factorization machines, Phys. Rev. Res. **2**, 013319 (2020).

[29] M. Jünger, E. Lobe, P. Mutzel, G. Reinelt, F. Rendl, G. Rinaldi, and T. Stollenwerk, Quantum annealing versus digital computing: An experimental comparison, ACM J. Exp. Algorithms **26**, 1 (2021).

[30] O. Parekh, J. Wendt, L. Shulenburger, A. Landahl, J. Moussa, and J. Aidun, Benchmarking adiabatic quantum optimization for complex network analysis, Technical Report Sandia National Lab. (April 2015), https://www.osti.gov/biblio/1459086.

[31] B. Yan and N. A. Sinitsyn, Analytical solution for nonadiabatic quantum annealing to arbitrary Ising spin Hamiltonian, Nat. Commun. **13**, 2212 (2022).

[32] A. D. King, J. Raymond, T. Lanting, R. Harris, A. Zucca, F. Altomare, A. J. Berkley, K. Boothby, S. Ejtemaee, C. Enderud, E. Hoskinson, S. Huang, E. Ladizinsky, A. J. R. MacDonald, G. Marsden, R. Molavi, T. Oh, G. Poulin-Lamarre, M. Reis, C. Rich, Y. Sato, N. Tsai, M. Volkmann, J. D. Whittaker, J. Yao, A. W. Sandvik, and M. H. Amin, Quantum critical dynamics in a 5,000-qubit programmable spin glass, Nature (London) **617**, 61 (2023).

[33] S. Mukherjee and B. K. Chakrabarti, Multivariable optimization: Quantum annealing and computation, Eur. Phys. J.: Spec. Top. **224**, 17 (2015).

[34] S. Boixo, V. N. Smelyanskiy, A. Shabani, S. V. Isakov, M. Dykman, V. S. Denchev, M. H. Amin, A. Y. Smirnov, M. Mohseni, and H. Neven, Computational multiqubit tunnelling in programmable quantum annealers, Nat. Commun. **7**, 10327 (2016).

[35] D. Sherrington, A simple spin glass perspective on martensitic shape-memory alloys, J. Phys.: Condens. Matter **20**, 304213 (2008).

[36] S. Kartha, T. Castán, J. A. Krumhansl, and J. P. Sethna, Spin-glass nature of tweed precursors in martensitic transformations, Phys. Rev. Lett. **67**, 3630 (1991).

[37] J. P. Sethna, S. Kartha, T. Castán, and J. A. Krumhansl, Tweed in martensites: A potential new spin glass, Phys. Scr. **1992**, 214 (1992).

[38] R. Vasseur and T. Lookman, Spin models for ferroelastics: Towards a spin glass description of strain glass, Solid State Phenom. **172–174**, 1078 (2011).

[39] E. Polatidis, M. Šmíd, I. Kuběna, W.-N. Hsu, G. Laplanche, and H. Van Swygenhoven, Deformation mechanisms in a superelastic NiTi alloy: An in-situ high resolution digital image correlation study, Mater. Des. **191**, 108622 (2020).

[40] Z. Wang, J. Chen, R. Kocich, I. P. Dolbnya, S. Tardif, L. Kunčická, J.-S. Micha, K. Liogas, O. V. Magdysyuk, I. Szurman, and A. M. Korsunsky, Grain structure engineering of NiTi shape memory alloys by intensive plastic deformation, ACS Appl. Mater. Interfaces **14**, 31396 (2022).

[41] L. Mañosa and A. Planes, Mechanocaloric effects in shape memory alloys, Philos. Trans. R. Soc. A **374**, 20150310 (2016).

[42] K. Otsuka and X. Ren, Physical metallurgy of Ti-Ni-based shape memory alloys, Prog. Mater. Sci. **50**, 511 (2005).

[43] T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer, Defining and detecting quantum speedup, Science **345**, 420 (2014).

[44] A. Lucas, Ising formulations of many NP problems, Front. Phys. **2**, 5 (2014).

[45] T. Lanting, A. J. Przybysz, A. Y. Smirnov, F. M. Spedalieri, M. H. Amin, A. J. Berkley, R. Harris, F. Altomare, S. Boixo, P. Bunyk, N. Dickson, C. Enderud, J. P. Hilton, E. Hoskinson, M. W. Johnson, E. Ladizinsky, N. Ladizinsky, R. Neufeld, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, S. Uchaikin, A. B. Wilson, and G. Rose, Entanglement in a quantum annealing processor, Phys. Rev. X **4**, 021041 (2014).

[46] J. Raymond, R. Stevanovic, W. Bernoudy, K. Boothby, C. C. McGeoch, A. J. Berkley, P. Farré, and A. D. King, Hybrid quantum annealing for larger-than-QPU lattice-structured problems, ACM Trans. Quantum Comput. **4**, 1 (2023).

[47] J. J. Berwald, The mathematics of quantum-enabled applications on the d-wave quantum computer, Not. Am. Math. Soc. **66**, 832 (2019).

[48] D-Wave Leap quantum cloud service, https://cloud.dwavesys.com.

[49] A. G. Khachaturyan, Theory of Structural Transformations in Solids (Dover Publications, Inc., Mineola, New York, 2008).

[50] Y. Wang, D. Banerjee, C. C. Su, and A. G. Khachaturyan, Field kinetic model and computer simulation of precipitation of $L_{12}$ ordered intermetallics from f.c.c. solid solution, Acta Mater. **46**, 2983 (1998).

[51] H. Lebbad, B. Appolaire, Y. Le Bouar, and A. Finel, Insights into the selection mechanism of Widmanstätten growth by phase-field calculations, Acta Mater. **217**, 117148 (2021).

[52] Y. Le Bouar, A. Loiseau, and A. G. Khachaturyan, Origin of chessboard-like structures in decomposing alloys. Theoretical model and computer simulation, Acta Mater. **46**, 2777 (1998).

[53] https://www.fz-juelich.de/ias/jsc.