domain: phys.qphy chim spi

text:
Roland Sandt
email: r.sandt@fz-juelich.de

Yann Le Bouar

Robert Spatschek

Quantum annealing for microstructure equilibration with long-range elastic interactions

We demonstrate the use and benefits of quantum annealing approaches for the determination of equilibrated microstructures in shape memory alloys and other materials with long-range elastic interaction between coherent grains and their different martensite variants and phases. After a one dimensional illustration of the general approach, which requires to formulate the energy of the system in terms of an Ising Hamiltonian, we use distant dependent elastic interactions between grains to predict the variant selection for different transformation eigenstrains. The results and performance of the computations are compared to classical algorithms, demonstrating that the new approach can lead to a significant acceleration of the simulations. Beyond a discretization using simple cuboidal elements, also a direct representation of arbitrary microstructures is possible, allowing fast simulations with currently up to several thousand grains.

Introduction

The modeling of microstructures is an important approach to the understanding, improvement and development of new materials for various applications. However, as mechanisms at different length and time scales are intimately linked, such descriptions and model implementations are typically challenging and require massive computational resources. Although phase field simulation approaches -the most prominent method for predicting microstructure evolution -strongly benefit from developments like the thin interface limit [START_REF] Karma | Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics[END_REF][START_REF] Karma | Quantitative phase-field modeling of dendritic growth in two and three dimension[END_REF] , nondiagonal phase field models [START_REF] Boussinot | Achieving realistic interface kinetics in phase-field models with a diffusional contrast[END_REF][START_REF] Wang | Modeling of dendritic growth using a quantitative nondiagonal phase field model[END_REF] and sharp phase field approaches [START_REF] Finel | Sharp phase field method[END_REF] , simulations containing large microstructural domains to obtain predictions with a certain statistical significance are rare, strongly limited by the available (super-)computer resources and their associated costs and energy consumption. Despite the enormous progress in this research field and the extended use of parallel computers and graphics cards for the simulations, limitations of the computational techniques remain a serious thread to the basic scientific progress and applied research.

One of the striking questions, which arises at the horizon of materials science modeling is how quantum computing will potentially change the simulation landscape in the future. However, at present a general-purpose quantum computer of sufficient size is not yet available. In the meantime, a technology known as quantum annealing (QA) [START_REF] Finnila | Quantum annealing: a new method for minimizing multidimensional functions[END_REF][START_REF] Brooke | Quantum Annealing of a Disordered Magnet[END_REF][START_REF] Kadowaki | Quantum annealing in the transverse Ising model[END_REF][START_REF] Morita | Mathematical foundation of quantum annealing[END_REF][START_REF] Rajak | Quantum annealing: an overview[END_REF] has emerged and is available on several sites worldwide. The use of such machines differs significantly from traditional gate based computers and therefore currently only specific problems can be handled by quantum annealers [START_REF] Warren | Mathematical methods for a quantum annealing computer[END_REF] . The concept of a quantum annealer is that its qubits are initialized in a well defined state which is described by a Hamiltonian with a unique ground state [START_REF] Johnson | Quantum annealing with manufactured spins[END_REF] . During the operation at cryogenic temperatures, this Hamiltonian is changed adiabatically such that the ground state converts into the one of the final, desired Hamiltonian [START_REF] Johnson | Quantum annealing with manufactured spins[END_REF][START_REF] Boixo | Evidence for quantum annealing with more than one hundred qubits[END_REF] , and therefore allows to perform global energy minimization computations efficiently. The structure of these Hamiltonians is a binary quadratic model, which can be expressed in terms of a quadratic unconstrained binary optimization or equivalently through an Ising model [START_REF] Warren | Mathematical methods for a quantum annealing computer[END_REF] . Due to this specific structure, so far, materials science related applications of this technology are still rare. Instead, actual research focuses mainly on the benchmarking and performance tests of quantum annealing compared to classical approaches [START_REF] Jünger | Quantum Annealing versus Digital Computing: An Experimental Comparison[END_REF][START_REF] Parekh | Benchmarking Adiabatic Quantum Optimization for Complex Network Analysis[END_REF][START_REF] Yan | Analytical solution for nonadiabatic quantum annealing to arbitrary Ising spin Hamiltonian[END_REF] . Some first applications in the field of biology and traffic research in the sense of optimization problems have been developed recently. Here, quantum annealing enables the efficient analysis of transcription factors in gene expression with combined machine learning algorithms [START_REF] Li | Quantum annealing versus classical machine learning applied to a simplified computational biology problem[END_REF] , identification of conformations of lattice protein models [START_REF] Perdomo-Ortiz | Finding low-energy conformations of lattice protein models by quantum annealing[END_REF] and their folding [START_REF] Irbäck | Folding lattice proteins with quantum annealing[END_REF] , detection of tree cover in aerial images [START_REF] Boyda | Deploying a quantum annealing processor to detect tree cover in aerial imagery of california[END_REF] , real-world traffic flow optimization problems [START_REF] Neukart | Traffic Flow Optimization Using a Quantum Annealer[END_REF] or control of automated guided vehicles [START_REF] Ohzeki | Control of Automated Guided Vehicles Without Collision by Qantum Annealer and Digital Devices[END_REF] . However, the usage of quantum annealing in materials science is not widespread and few publications correspond to phase transitions in the transverse field Ising model [START_REF] Harris | Phase transitions in a programmable quantum spin glass simulator[END_REF] , the investigation of critical phenomena in frustrated magnets via the Shastry-Sutherland Ising model [START_REF] Kairys | Simulating the Shastry-Sutherland Ising Model Using Quantum Annealing[END_REF] , Monte-Carlo sampling [START_REF] Sandt | Efficient low temperature Monte Carlo sampling using quantum annealing[END_REF] and the automated materials design of metamaterials [START_REF] Kitai | Designing metamaterials with quantum annealing and factorization machines[END_REF] . The purpose of the present paper is therefore to demonstrate that this novel technology can indeed lead to completely new possibilities beyond the existing and commonly used descriptions for the modeling of microstructures.

In order to be as explicit and illustrative as possible, we demonstrate here the case of coherent solid state transformations involving austenitic and martensitic phases, where the latter are allowed to appear in different variants. Such transitions play a role for shape memory alloys like NiTi, which can be deformed easily at low temperatures, but heating to higher temperatures lets the material return to its previous, trained shape [START_REF] Otsuka | Physical metallurgy of Ti-Ni-based shape memory alloys[END_REF] . The modelling and mapping of shape memory alloys to spin glass systems was previously established in several studies [START_REF] Sherrington | A simple spin glass perspective on martensitic shape-memory alloys[END_REF][START_REF] Kartha | Spin-Glass Nature of Tweed Precursors in Martensitic Transformations[END_REF][START_REF] Sethna | Tweed in Martensites: A Potential New Spin Glass[END_REF][START_REF] Vasseur | Spin Models for Ferroelastics: Towards a Spin Glass Description of Strain Glass[END_REF] and can here be exploited for QA applications. In the following we will mainly stick to the terminology of the shape memory alloys but emphasize that similar approaches can be used to model e.g. the transformation and deformation behavior in steels, ferroelastic materials, as well as phase changes in solid electrolytes for rechargeable batteries. The particular aspect that plays a central role here are the anisotropic long-range elastic interactions, which are common for solid state transformations [START_REF] Khachaturyan | Theory of Structural Transformations in Solids[END_REF] , and therefore the ground state configuration does not only depend on phase concentrations and fractions, but also on the detailed microstructural arrangement of phases and grains. In a typical phase field simulation [START_REF] Chen | Phase-field models for microstructure evolution[END_REF] , the microstructural evolution is solved together with the relaxation of the mechanical deformations in the spirit of a continuum description, which leads to very long simulation times. Here, we show that the separation of the discrete degrees of freedom for the variant distribution of martensitic phases from the continuous development of the microstructure and the QA treatment allow to drastically increase the performance of the computations and therefore to simulate significantly larger, application relevant systems compared to existing approaches.

Results

One-dimensional model. For a simplified 1D model we consider only a "martensitic" phase which is assumed to exist in two different variants. Hence the microstructure consists of a line of grains of these variants, as depicted in the inset of Fig. 1(a). To be explicit, we assume that both variants have a stress free strain (eigenstrain), which leads to a shear deformation relative to the austenitic mother phase, and denote these variants by state variables s i = ±1. As in the end we will map the description to a one-dimensional Ising model, we also use here the terminology of "spins" with two possible alignments in the spirit of a magnetic model. As each of the variants leads to a shearing of the cell, we get an overall stress free deformation of this line (compared to the shear strain free austenite), depending on the spin configuration. We assume that all grains have the same height d, the same elastic constants, and opposite shear eigenstrain ±ε 0 . As one can readily see from the inset of Fig. 1(a), the stress free equilibrium position of the top grain x 0 depends only on the number of variants N + with orientation s i = +1 and N -with s i = -1, but not on the individual arrangement, which is a particularity of the simplified 1D model and the chosen eigenstrain. Hence, for a fixed number N = N + + N -spins in a row, the macroscopic stress free strain is ε¯ = (N + -N -)ε 0 /N, which leads to x 0 = Ndε ¯. If an external deformation is enforced, i.e. x ̸ = x 0 the elastic energy is F el = µ eff (xx 0 ) 2 with an effective shear modulus µ eff . Obviously, the elastic energy is minimized if the spin configuration is such that x = x 0 , which implies (N + -N -) min = x/ε 0 d, up to the point of saturation, where all spins are aligned. This expression serves as reference for the comparison with the numerical minimization approaches below. We note that we neglected at this stage the discrete nature of the variants, which means that the integer value N + -N -should be as close as possible to the continuum value (N + -N -) min above. Although the energy in the simple 1D model does not depend on the arrangement of the variants but only on the total numbers N + and N -= N -N + , we formulate the model here on the level of the individual "spins" s i for the later extension towards higher dimensions and the use of the quantum annealer. Hence we get N + -N -= ∑ i s i . Inserting this into the elastic energy expression yields

F el = µ eff ε 0 2 d 2 ∑ i, j s i s j -2µ eff xε 0 d ∑ i s i + µ eff x 2
, where the summations run over all spins. For the implementation on a quantum annealer, we need to bring this to the Ising form of a Hamiltonian H with

H = ∑ i h i s i + ∑ i< j J i j s i s j + H 0 , (1)
where the first term corresponds to the interaction with an external magnetic field h i and the second term to a spin-spin interaction, which favors ferromagnetic (antiferromagnetic) ordering in case that the coupling constant J i j is negative (positive). The last, spin-independent term H 0 is only an irrelevant additive constant. From the comparison of the two above expressions we identify h i = -2µ eff xε 0 d and J i j = 2µ eff ε 2 0 d 2 . First, we note that the external deformation is here analogous to the magnetic field in the Ising description. Second, the spin-spin interaction term J i j is positive, hence favoring "antiferromagnetic ordering". Also, this term is independent of the spin numbers i, j, which means that this interaction does not depend on the distance between the grains. In other words, the elastic interaction depends only on the averaged "magnetization" N + -N -, which implies a mean field interaction.

The goal of the formulation is to minimize the elastic energy and to find the optimal spin or variant configuration {s i }. To this end, we use three different numerical approaches (see methods section), and the results are compared to the analytical solution above: First, a brute force approach iterates over all spin configurations to find the energetic minimum exactly, second we use simulated annealing as probabilistic ground state finder, and finally the quantum annealing approach. Fig. 1(a Comparison between the results obtained by numerical minimization (solid lines) versus the analytical theory for an infinite and continuous system (dashed line). For large displacement, all "spins" are aligned and therefore the "magnetization" saturates. The inset shows a sketch of the one-dimensional arrangement of martensite variants s i = +1 (red) and s i = -1 (green). The bottom row is fixed to position x = 0, whereas the top grain has a mean position x 0 in the stress free state. If an additional external stress or strain is applied, the top layer is moved to position x, and the entire microstructure is sheared to the dashed configuration. (b) Elapsed computation time as a function of the number of grains. Different methods and algorithms are compared. Dashed parts of the QA curve belong to the regime of chain breaks. For large system sizes, only the hybrid quantum annealing approach remains feasible, showing an almost constant computing time need for less than 1000 spin variables (inset).

the resulting "magnetization" (N + -N -)/N, i.e. the average variant orientation, as function of applied displacement x/dNε 0 , which corresponds to the magnetic field in the Ising model. As expected, the results agree with the analytical theory up to the aforementioned discretization effect, which becomes less pronounced for large grain numbers. For high displacements saturation sets in when all variants are de-twinned, which means that all spins are either in the state +1 or -1. We note that for the investigated number of spins all used algorithms lead to the same energy minimum, which confirms that also the probabilistic approaches indeed find the global minimum states. Fig. 1(b) shows the required computation time for the different methods and algorithms as a function of the number of grains N. All conventional algorithm implementations are based on single core computations without parallelization and are mainly shown for a qualitative comparison, as the focus of the investigations is on the quantum annealing approach. For the latter, we use quantum processing unit (QPU) implementations up to the highest possible number of spins (typically N ≈ 170 for the Pegasus architecture [START_REF] Boothby | Next-Generation Topology of D-Wave Quantum Processors[END_REF] of a D-Wave machine). The brute force approach, where iterations over all spin configurations are run, has the highest computation time. Even at small spin systems of around N ≈ 40 the elapsed user time was too large for practical applications due to the simulation time scaling ∼ O(2 N ). The pure quantum annealing method produces the fastest results and ends up with an almost constant elapsed QPU access time. Overall, the computations for N ≈ 150 are roughly three orders of magnitude faster than for the other classical approaches. Beyond around N ≈ 50 spins, so called chain breaks [START_REF] Grant | Benchmarking embedded chain breaking in quantum annealing[END_REF] occur occasionally. They result from the need to encode strongly coupled spins as a single logical spin. Ideally, these spins should represent the same state as the individual spins, but in practise this identity can be violated. To avoid this issue and to simulate even larger systems, which are essential for higher dimensional modeling in the following sections, hybrid classical and quantum annealing approaches can be used, which combine pure QA with conventional minimization approaches [START_REF] Raymond | Hybrid quantum annealing for larger-than-QPU lattice-structured problems[END_REF] . The numerical results in Fig. 1(b) show an increase of the computation time of the hybrid solver compared to the pure QA, but the relative acceleration compared to the classical algorithms becomes even more striking. For the hybrid solver, the elapsed computation time is essentially independent of the number of spin variables and increases only beyond 10 3 grains to several seconds. Altogether, the hybrid QA is clearly the fastest approach for large grain numbers and is therefore used in the following two-dimensional simulations.

Transformations in higher dimensions. For the determination of the linear elastic energy beyond one dimension, we consider coherent precipitates of different variants which form inside the matrix. In this way, the entire material can be considered to consist of small entities (in the following denoted as grains), which can be in one of the different martensitic states.

3/12

The simplest possible (cartesian) discretization is to use small cuboidal grains with edge length a. All grains are assumed to be coherent (the elastic displacements and tractions are continuous at the interfaces between the grains), and we use homogeneous elasticity, i.e. we ignore differences in the elastic constants between the different phases or variants. This has the consequence that the elastic energy reduces to combinations of pairwise interactions between all grains [START_REF] Johnson | Elastic interactions and stability of misfitting cuboidal inhomogeneities[END_REF] .

For demonstrational purposes we perform here two-dimensional simulations in a plane strain setup, but a transfer to three dimensions is straightforward. In particular, the annealer part does not depend on the dimensionality of the description. The qualitatively new aspect beyond 1D is the appearance of distance and orientation dependent "spin-spin" interactions, which decay only slowly with the distance between the grains, and therefore leads to fully populated matrices J i j . As it turns out that an accurate determination of the elastic interaction energy is essential for a precise prediction of the equilibrium microstructure, we use Fourier transformation approaches with periodic boundary conditions as outlined in the methods section. As boundary conditions, we use either vanishing average stress in the periodic volume V , ⟨σ i j ⟩ = 1

V σ i j (r) dr = 0, or, similarly to the 1D description, a given average strain ⟨ε i j ⟩. We employ in the following for simplicity isotropic elasticity, which is e.g. described by the Lamé coefficient λ and the shear modulus µ, i.e. the stress-strain relationship reads

σ i j = 2µ(ε i j -ε (0) i j ) + λ δ i j (ε kk -ε (0) kk ),
where implicit summation over repeated indices is used. The position dependent eigenstrain ε (0) i j (r) is known for a given microstructure with fixed phase dependent stress free strains (relative to the austenitic mother phase), ε (0) i j (r) = θ (r)ε 0 i j , where the indicator function θ is zero in the austenite and either +1 or -1 in the two considered martensite variants. For a given microstructure, the elastic energy can then be computed in reciprocal space, as shown in the methods section. For the formulation as Ising model we discretize our microstructure using small non-overlapping cuboidal grains as discussed above and assign a "spin" s i to each of them like before, such that the indicator field becomes a superposition θ (r) = ∑ i s i θ i (r), where θ i equals one inside the corresponding square and is zero outside. Therefore, the elastic energy decomposes into pairwise interactions (for i ̸ = j) and self-energy terms (for i = j)

E i, j = s i s j 1 2V dr dr ′ B(r -r ′ )θ i (r)θ j (r ′ ), (2)
where the integral kernel B(r) is defined through the inverse of the elastic Green's function. Hence, it is sufficient to perform the Fourier transform energy calculations for all pairs of the same martensite variant s i = s j = 1 on the discrete lattice sites in the volume V ; for periodic boundary conditions and identical grain shapes, it is sufficient to calculate the elastic interaction energy between a reference grain and all the other grains, due to translation invariance. In case of fixed average strain boundary conditions, an additional homogeneous term appears (see methods section), contributing both to the spin-spin interaction J i j as well as to the magnetic field term h i , which is absent for zero average stress boundary conditions. The resulting fully populated matrix of coupling constants with both positive and negative entries has similarities to spin glass systems with random couplings, which have been investigated in the literature with conventional approaches, see e.g. [START_REF] Ray | Sherrington-Kirkpatrick model in a transverse field: Absence of replica symmetry breaking due to quantum fluctuations[END_REF] .

For the simplest case that the eigenstrain is purely dilatational and isotropic the Bitter-Crum theorem applies and the total energy depends only on the volume fraction of the martensite variant, where no interaction between the grains is present and only a self energy term remains [START_REF] Fratzl | Modeling of Phase Separation in Alloys with Coherent Elastic Misfit[END_REF] .

For a nontrivial elastic interaction and the link to the previous 1D description, we consider a shear transformation strain with ε x 0 y = ε 0 yx = ε 0 , where all other components vanish. In this case, we obtain a distance and orientation dependent interaction as depicted in Fig. 2(a), which is computed here for the case of vanishing average stress, ⟨σ ik ⟩ = 0. Here and in the following parts the Poisson ratio is chosen as ν = 1/4 (i.e. λ = µ). The interaction energy is obtained by subtracting the elastic self energies E self for each of the two (isolated) martensite grains inside the austenitic matrix from the total elastic energy E el of the two-grain arrangement, i.e. E int = E el -2E self , to normalize the interaction energy such that it decays to zero for large grain separations. For short distances, a transition between attraction and repulsion is found for the ⟨100⟩ direction, whereas a purely repulsive interaction results for the diagonal ⟨110⟩ directions. Due to the periodic boundary conditions, the result depends on the system size V = L x × L y , as the grains also interact with their periodic images, hence r ≪ L x ,L y is required to observe the decay of the interaction.

We note that in two dimensions the interaction energy decays asymptotically as r -2 , whereas in three dimensions it scales as r -3 in large systems, which follows from the elastic Green's function [START_REF] Landau | Theory of Elasticity[END_REF] . For the quantum annealer implementation, the interaction energies are needed only for the discrete lattice points (symbols on the curves). Although the decay of the elastic interaction may suggest to cut it off beyond a certain distance in real space, it turns out that such an approach is inappropriate, as it leads in the end to invalid equilibrium microstructures, and it is therefore essential to keep all interaction terms J i j with high precision to avoid spurious effects. We note that the formulation on the quantum annealer does not depend on the dimensionality, therefore the scaling plot in Fig. 1(b) applies here as well.

Based on the calculation of the elastic interactions, we obtain from the Ising model implementation on the quantum annealer with hybrid solver stripe patterns in ⟨100⟩ direction as equilibrium structures. These patterns are irregular in the sense that

E int /λa 3 ε 0 2 r/a [100] [010] <110> Figure 2.
Interaction energies of two grains of equal variant type (s i = s j ). Interaction energies in the case of (a) shear eigenstrain and vanishing average stress and (b) tetragonal eigenstrain. The interaction energy per length is given in units of λ a 3 ε 2 0 , and the computations were done using a system size of L x /a = L y /a = 50, where a is the edge length of the grains. At distance r/a = 0 the grains touch each other. The symbols on the continuous curves indicate the information for the interaction at discrete lattice sites, which is actually used in the annealer simulations. the widths of the stripes are not uniform. This is in analogy to the 1D model, which was discussed above, where we found that the arrangement of the two variants is not determined. This coincidence, which is physically expected, is nontrivial from the model formulation, as i) in the 1D model we had a distance independent interaction in the discretized model, where here the interaction is significantly more complex, but adds up to the same effective descriptions for the periodic arrangement; ii) a rotation of the pattern by 90 degree is possible and sometimes obtained from the optimal configuration due to the discrete rotational symmetry; iii) the fixing of the average stress compared to the given average strain in the 1D formulation can lead to unequal distributions of the different variants. In particular, for the presently considered absence of an external strain (implying a vanishing magnetic field in the Ising terminology), there is no constraint of the sort ⟨s i ⟩ = 0 for the average spin alignment. All stripe configurations are energetically equivalent, which includes the possibility of a single variant configuration. These results therefore confirm simultaneously the accuracy of the elastic interaction calculation with the pairwise decomposition as well as the ability of the quantum annealer to identify the true ground state configurations.

As next example, we use a tetragonal eigenstrain with the only nonvanishing components ε 0 xx = -ε 0 yy = ε 0 zz = ε 0 . First, we consider again the situation with vanishing average stress, ⟨σ i j ⟩ = 0. The corresponding interaction energy is shown in Fig. 2(b) for ν = 1/4. In this case, the equilibrium microstructure is trivial and consists of a single variant, as in this case the elastic energy is zero for the periodic system. Therefore, the situation differs from the previous shear transformation, where also lamellar arrangements with both variants lead to stress free situations. The reason is that any interface between two variants leads to a mismatch between adjacent variants for the tetragonal transformation, and therefore such a situation is energetically unfavorable here. A change of boundary conditions to vanishing average strain, ⟨ε i j ⟩ = 0, alters the situation, since then arrangements with equal amounts of both variants are preferred, as this lowers the volumic part of the elastic energy. In this case, we find regular inclined stripes as equilibrium pattern, as shown in Fig. 3(a). Again, the solution is not unique; in particular, due to translation invariance, the annealer returns also configurations where the stripes are shifted. Also, a switch of the sign of the inclination angle φ (see definition in the figure) leads to energetically equivalent solutions. However, we do not find ground state configurations which lead to different (absolute) inclination angles or strip widths or even irregular variations of the latter, contrary to the shear transformation case before.

The reason for the observed ground state morphologies is a combination of continuum elasticity effects, the granular structure of the material and constraints induced by periodic boundary conditions. Fig. 3(b) shows the computed elastic energy for different numbers of regular arrangements of stripes in the periodic system as function of the inclination angle φ . Here we see a pronounced influence of the grain size, as the elastic energy of configurations with regular stripe pairs with a discretization by 50 × 50 grains (squares in the figure) is higher than for corresponding cases with very fine grains, where discretization effects do not play a role anymore (smooth curves). The oscillating nature is due to the periodic boundary conditions, as improper choices of the inclination angle lead to discontinuities of the stripe patterns at the boundaries, which is energetically unfavorable. Therefore, continuous patterns correspond to the stationary points of the curves. For specific angles, the curves for 5/12 The width of the stripes is uniform, consisting of grains with configuration s i = +1 (red) and s i = -1 (green). (b) Elastic energy of stripe patterns with different inclination angles φ . The solid curves correspond to smooth stripes (the grain size a/L x , a/L y ≪ 1 is negligible) and show a pronounced stationary point for inclinations for which the pattern repeats periodically without kinks at the boundaries. The squares correspond to situations with the same number of stripes, where the system is discretized by 50 × 50 quadratic grains, leading to pronounced aliasing effects, and the resulting elastic energy is higher than for the smooth stripes. This shifts the energetic minimum for 6 stripe pairs at φ ≈ 40 • to a lower angle φ ≈ 33 • with 3 stripe pairs. The infinite system size limit for smooth stripes is depicted as black dotted curve. three and six stripe pairs meet at local minima, which is a consequence of the scale invariance of linear elasticity. From the smooth, continuum limit curves one would conclude that an angle of about φ ≈ 40 • should lead to the energetically lowest configuration (absolute minimum of the smooth red curve). Moreover, in the limit of infinite systems, where periodic boundary conditions do not play a role anymore, an analytical treatment is possible, leading to the energy expression E ∞ el = V B(n)/2 for equal volume fraction of the two variants with

B(n) = 4µ λ + 2µ ε 2 0 [(3λ + 2µ) -2(3λ + 2µ)n 2 + 4(λ + µ)n 4 ]
with n = cos φ . Energy minimization gives the optimal angle φ = cos -1 5/8 ≈ 37.8 • , see Fig. 3(b) (minimum of the black dotted curve). However, these predictions disagree with the finding from the quantum annealer, which favors a configuration with three stripe pairs at a lower angle of φ ≈ 33 • . This observation can be understood by consideration of the granular structure of the patterns investigated here, as the microstructure in the annealer simulations consists of 50 × 50 square grains. First, the explicit appearance of the length scale a breaks the scale invariance of the periodic pattern, and therefore the minima of the energy curves belonging to the discrete microstructures (squares in Fig. 3(b)) do not coincide anymore at the local minima. Additionally, with increasing inclination antialiasing effects of the patterns become more relevant, and therefore the energy curves show an increasing disagreement with the continuum limit curves. As a result, the energetic minimum in the discrete microstructure indeed shifts toward a configuration with three stripe pairs at φ ≈ 33 • (absolute minimum of the blue squares in Fig. 3(b)), which is in agreement with the prediction of the quantum annealer. Consequently, details of the granular structure can change the energetics compared to a full continuum approximation, especially since many local minima of the elastic energy are located close to each other.

Variant selection in realistic microstructures. The approach presented above is not limited to mutually interacting cuboidal grains, but can also be applied to realistic microstructures. To illustrate the procedures, we have generated a microstructure consisting of N = 400 grains using a Voronoi tesselation [START_REF] Rycroft | Voro++: A three-dimensional voronoi cell library in c++[END_REF] . Each grain is allowed to take one out of two martensite variants with the tetragonal eigenstrain tensor, and we pre-compute all mutual elastic interactions between them. We note that contrary to the case with the cuboidal grains in a periodic array here we cannot exploit translational invarince due to the different shapes of the grains, and hence these elastic interaction energy calculations scale here as O(N 2 ) instead of O(N) before, although we still use periodic boundary conditions. Additionally, we consider now arbitrary given external strains ⟨ε i j ⟩, which leads to the appearance of a "magnetic" term like in the one dimensional description. With that, we can predict the equilibrium

6/12 (a) (b) (c) (d) (e) (f)
⟨ε xx ⟩/ε 0 = 1.3.
variant distribution within the microstructure using the hybrid quantum annealer, and this step is typically executed within a few seconds of runtime.

Examples for the equilibrium microstructures are shown in Fig. 4 as function of the externally applied strain ⟨ε xx ⟩, whereas the other average strain components vanish. The observed microstructures are indeed similar to what we have found before using the square discretization, although here the band widths and orientation deviate from the previous case due to microstructural details and the smaller number of grains, and these effects can be explained using an analysis similar to the one done for Fig. 3(b). We note that in these microstructures all grains have the same orientation, and therefore the application of a tensile strain strongly favors the selection of the grain variant s i = +1 (for a compressive situation we observe the opposite behavior), and we find a full alignment of all variants in the last snapshot.

Additionally, we have performed the same analysis for grains with uniformly distributed random orientation, which implies a rotation of the local transformation strain tensor, see Fig. 5 for the grain orientations and for the variant selection. Here, also the equilibrated spatial distribution of the variants appears to be irregular. Application of a tensile strain again favors the "alignment" of the variant, but this time even for high strains not all grains select the same variant, which is due to the local rotation. In fact, a grain, which is rotated by 90 • with respect to the straining direction has a preference to be in variant state s i = -1, as then the direction of expansion is aligned with the external tensile strain. This can be clearly seen e.g. in Fig. 5(c) for the highest tensile strain in x direction, where the remaining patches with "spin" s i = -1 correspond to the grains with orientation close to π/2 (or 3π/2). We emphasize that for a given microstructure (shapes of all grains) the mutual elastic grain-grain interactions have to be computed only once. As mentioned before, this step has to be done with high precision, and consequently this is the step which demands the highest amount of computing time. After that, all changes of the external boundary conditions affect only the k = 0 mode contributing to the interactions J i j and h i , and these terms can be calculated analytically (see methods section). As each hybrid quantum annealing calculation typically requires only a few seconds, the entire microstructural change during mechanical loading can be calculated extremely fast.

Discussion

The central result of the present paper is the shown optimization of microstructures via quantum annealing, exhibiting a clear performance advantage of the novel approach compared to conventional energy minimization strategies. The brute force approach is not recommended, whereas optimized simulated annealing algorithms produce good results. However, quantum annealing represents the by far fastest method for optimization problems, particularly for systems with high numbers of grains (spins) and non-vanishing coupling constants and biases, and allows the determination of ground state configurations for system sizes, which are not accessible for the classical algorithms on reasonable computing timescales.

For a system consisting of N grains, we need to compute O(N 2 ) spin-spin interactions. These elastic energy calculations have to be done with high accuracy, and therefore they dominate the overall computing time. After that, O(2 N ) spin configurations need to be compared to identify the equilibrium configuration. For the conventional algorithm, this combinatorial step dominates the total computational effort already for low values of N. In contrast, with the quantum annealer or its hybrid variant the computation time for the minimization of the Ising energy expression is completely negligible compared to the elastic interaction energy computations. Hence, we have demonstrated that QA is able to drastically optimize the search for microstructural equilibrium states in solid phases with long-range elastic interactions. Already today, the usage of hybrid quantum annealing enables the computation of microstructures with several thousand grains which all interact with each other, which is essential for a realistic modeling of microstructures inside various materials.

For many application relevant investigations, it is critical to understand whether and how models can be formulated that they are suitable for quantum computing. We have demonstrated this here for the particular case of long-range elastic interactions. Extensions towards the consideration of interfacial energy, multiple martensite variants, anisotropic elasticity, orientation relationships between grains and phases, and different spatial dimensions are obvious, as they do not conceptually influence the presented strategy of formulating the problem in terms of an Ising model. Inhomogeneous elasticity and the proximity to free surfaces can effectively lead to many-body interactions, for which perturbative extensions or the introduction of product spin variables are promising directions [START_REF] Johnson | Elastic interactions and stability of misfitting cuboidal inhomogeneities[END_REF][START_REF] Spatschek | Scale bridging description of coherent phase equilibria in the presence of surfaces and interfaces[END_REF] . Beyond the purely elastic effects, further potential applications comprise phase changes in multi-phase solid state batteries, phase transformations in high strength steels or other materials like ferroelectrics. Overall, the separation of continuous and discrete degrees of freedom and the quantum treatment of the latter may also be beneficial for hybrid phase field and quantum annealing descriptions which combine a variant selection with a grain morphology evolution in an efficient way to drastically reduce computing time demands of existing approaches for application relevant sample sizes.

Methods

Quantum annealing. Like general purpose quantum computers, a quantum annealer is built from qubits, which here store and process information using superconducting loops. A clockwise or anticlockwise circulating current in such a loop represents different spin states [START_REF] Johnson | Quantum annealing with manufactured spins[END_REF] . In each qubit superconducting loops interact with external flux biases, which allows to construct an 8/12 energy landscape, where the fluxes influence barrier height and energy di fference [START_REF] Johnson | Quantum annealing with manufactured spins[END_REF] . At the start of the co mputation, the system is initialized in the ground state of a known Hamiltonian H 0 ∼ -∑ i σ i x with Pauli matrices σ i , i.e. a strong transverse magnetic field 13, 4 3 . During the annealing process, the Hamiltonian is turned into the desired one based on an Ising model 11 H p = ∑ i h i s i + ∑ i< j J i j s i s j with spin states s i = ±1, bias h i and coupling strength J i j between qubits i and j, for which an energetic minimum is sought, min {s i =±1} H p . Both Hamiltonians do not commute [START_REF] Warren | Mathematical methods for a quantum annealing computer[END_REF] , and the time of the initial Hamiltonian to adopt the low energy state is sufficiently large to ensure the validity of the adiabatic theorem of quantum mechanics [START_REF] Lucas | Ising formulations of many NP problems[END_REF] , which states that a system remains in its eigenstate, if changes occur adiabatically. Notice that the quantum annealing employs tunneling to leave metastable regions, contrary to the simulated annealing [START_REF] Finnila | Quantum annealing: a new method for minimizing multidimensional functions[END_REF] . Another important quantum mechanical principle in quantum annealing is the entanglement and the usage of entangled states inside quantum annealing processors (QPU) [START_REF] Lanting | Entanglement in a Quantum Annealing Processor[END_REF] .

As in practise this approach does not always deliver the lowest energy state, especially if energetically close low energy states exist, a suitable number of repetitions is made and the configuration with the lowest detected energy is t aken. If the Ising problems do not match the architecture of the QPU, a subgraph of coupled qubits, know as chains, cover one variable of the problem in the so called minor embedding [START_REF] Raymond | Hybrid quantum annealing for larger-than-QPU lattice-structured problems[END_REF][START_REF] Choi | Minor-embedding in adiabatic quantum computation: II. Minor-universal graph design[END_REF] . Additionally, for huge systems hybrid quantum annealing exploits classical algorithms and the interplay with quantum annealing in areas of high computational demands using a QPU coprocessor working with generic parameters for up to 11616 spin variables on the D-Wave Advantage system [START_REF] Raymond | Hybrid quantum annealing for larger-than-QPU lattice-structured problems[END_REF][START_REF] Berwald | The Mathematics of Quantum-Enabled Applications on the D-Wave Quantum Computer[END_REF] . In practise, the D-Wave framework Leap 48 allows a direct formulation in terms of a problem Ising Hamiltonian.

Brute force minimization. For N spins we compute the energy of all 2 N possible configurations to determine the minimum. This deterministic approach delivers the true ground state energy but has a high computational effort.

Simulated annealing. For this probabilistic approach 49 a random initial configuration is chosen. A new candidate configuration, which we generate here by a single spin flip, is accepted if its energy is lower than the previous value. If the energy is higher by an amount ∆E, the configuration is accepted with a probability given by the Boltzmann factor exp(-∆E/T ), in order not to get stuck in local energy minima. During the simulation, the temperature T is reduced according to a specific cooling strategy, in order to converge towards an energetic minimum at the end of the simulation. As our main goal is not to maximize the performance of the (classical) algorithms but rather to demonstrate the general scaling behavior, we refrain from a detailed discussion of the parameter optimization of the probabilistic simulated annealing approach. This includes in particular the use of suitable stopping criteria when no further reduction of the energy occur, as well as the use of problem adapted cooling strategies. For the simulated annealing approach we use single spin flips trials in each iteration, and the temperature T is decreased each time by ∆T /µ eff ε 2 0 d 2 = 10 -6 , which delivers a good performance for large system sizes. The simulations are stopped after a fixed number of 10 7 steps, which is optimized for the largest considered spin system with N = 150 in Fig. 1(b), leading to a scaling of the computation time ∼ N 2 due to the calculation of the interaction energy.

Elasticity. We solve the elastic problem of a multi-grain setup with homogeneous linear elasticity, i.e. all variants and phases are assumed to have the same elastic constants. Also, coherent interfaces are assumed, which means continuity of displacements at the interfaces. The martensite variants have different stress free strains (or eigenstrains) compared to the mother austenite phase, hence the stress-strain relation reads for general linear elasticity σ i j = λ i jkl (ε klε (0) kl ), where ε (0) kl (r) is the local stress free strain tensor and λ i jkl the elastic tensor. We determine the elastic equilibrium configuration, which obeys the condition ∂ σ i j /∂ x j = 0 in bulk domains and the continuity of normal stresses at interfaces, using Fourier transformation approaches [START_REF] Khachaturyan | Theory of Structural Transformations in Solids[END_REF] . From that, the elastic energy can be computed in reciprocal space as 32

E el = V 2 ∑ k̸ =0 B(n)| θ (k)| 2 (3)
for a periodic system with vanishing average stress as boundary condition, where θ (k) is the Fourier transform of the indicator field θ (r) and B(n) with n = k/k equals B(n) = σ 0 i j ε 0 i jn i σ 0 i j Ω jk σ 0 kl n l with σ 0 i j = λ i jkl ε 0 kl . Here, Ω i j (n) is the normalized Green tensor for displacements, which is defined through its inverse as Ω -1 ik = λ i jkl n j n l . The summation in Eq. ( 3) is over discrete vectors due to the periodic boundary conditions in real space. The summation is in principle infinite, and can be efficiently computed using the decoration technique [START_REF] Rodney | Phase field methods and dislocations[END_REF] . For average strain boundary conditions, i.e. a prescribed value of ⟨ε i j ⟩, an additional homogeneous (k = 0) contribution appears in Eq. (3), which reads E hom = V λ i jkl (⟨ε i j ⟩ -⟨ε (0) i j ⟩)(⟨ε kl ⟩ -⟨ε (0) kl ⟩)/2, which can be calculated analytically.

3 Figure 1 .

 31 Figure 1. Results of the one-dimensional model comparing different numerical and analytical methods. (a) Mean variant orientation (N + -N -)/N as function of the displacement x/dNε 0 .Comparison between the results obtained by numerical minimization (solid lines) versus the analytical theory for an infinite and continuous system (dashed line). For large displacement, all "spins" are aligned and therefore the "magnetization" saturates. The inset shows a sketch of the one-dimensional arrangement of martensite variants s i = +1 (red) and s i = -1 (green). The bottom row is fixed to position x = 0, whereas the top grain has a mean position x 0 in the stress free state. If an additional external stress or strain is applied, the top layer is moved to position x, and the entire microstructure is sheared to the dashed configuration. (b) Elapsed computation time as a function of the number of grains. Different methods and algorithms are compared. Dashed parts of the QA curve belong to the regime of chain breaks. For large system sizes, only the hybrid quantum annealing approach remains feasible, showing an almost constant computing time need for less than 1000 spin variables (inset).

Figure 3 .

 3 Figure 3.Resulting stripe patterns for tetragonal eigenstrain. (a) Equilibrium structure with three stripe pairs (counted along the horizontal axis) in a system consisting of 50 × 50 cuboidal grains. A vanishing mean strain, ⟨ε i j ⟩ = 0, is imposed. The width of the stripes is uniform, consisting of grains with configuration s i = +1 (red) and s i = -1 (green). (b) Elastic energy of stripe patterns with different inclination angles φ . The solid curves correspond to smooth stripes (the grain size a/L x , a/L y ≪ 1 is negligible) and show a pronounced stationary point for inclinations for which the pattern repeats periodically without kinks at the boundaries. The squares correspond to situations with the same number of stripes, where the system is discretized by 50 × 50 quadratic grains, leading to pronounced aliasing effects, and the resulting elastic energy is higher than for the smooth stripes. This shifts the energetic minimum for 6 stripe pairs at φ ≈ 40 • to a lower angle φ ≈ 33 • with 3 stripe pairs. The infinite system size limit for smooth stripes is depicted as black dotted curve.

Figure 4 .

 4 Figure 4. Resulting equilibrium variant distribution with uniform grain orientation. The microstructures consist of 400 grains and tensile strain is applied in horizontal (x) direction. Red (green) grains correspond to variant s i = +1 (s i = -1). The tensile strain is (a) ⟨ε xx ⟩/ε 0 = 0, (b) ⟨ε xx ⟩/ε 0 = 0.1, (c) ⟨ε xx ⟩/ε 0 = 0.5, (d) ⟨ε xx ⟩/ε 0 = 0.9, (e) ⟨ε xx ⟩/ε 0 = 1.1 and (f) ⟨ε xx ⟩/ε 0 = 1.3.

7 / 12 Figure 5 .

 7125 Figure 5. Resulting equilibrium variant distribution with random grain orientation. (a) Grain orientation map corresponding to the microstructures. In the color bar the grain rotation angle is given in radian (modulo π due to symmetry). The rotation axis is along the [001] direction. The microstructures consist of 400 grains and tensile strain is applied in horizontal (x) direction. The grains have a random orientation, which is the same for all cases, based on a uniform distribution. The tensile strain in horizontal direction is (b) ⟨ε xx ⟩/ε 0 = 0 and (c) ⟨ε xx ⟩/ε 0 = 2.1. Red (green) grains correspond to variant s i = +1 (s i = -1).

Acknowledgements

This research was funded by the German Federal Ministry of Education and Research (BMBF) via the project Meet HiEnD 3 and the Helmholtz project ZeDaBase. Open access was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -491111487. The authors gratefully acknowledge the Jülich Supercomputing Centre (https://www.fzjuelich.de/ias/jsc) for funding this project by providing computing time on the D-Wave Advantage TM System JUPSI through the Jülich UNified Infrastructure for Quantum computing (JUNIQ).

Data Availability

Data that was obtained during this project will be made available by the corresponding author upon request.

9/12

Author contributions statement R.Sa., Y.L.B. and R.Sp. contributed to analytical and numerical calculations, methodology, visualization and analysis. All authors contributed to the writing of the manuscript.

Competing Interests

The authors declare no competing financial or non-financial interests.

12/12


Jeoffrey Tourneur

Loïc Joanny

Loïc Perrin

Stéphane Paul

Bruno Fabre

Efficient and Highly Stable 3D Printed NiFe and NiCo Bifunctional Electrodes for Practical HER and OER

Keywords: Water splitting, Hydrogen Evolution reaction (HER), Oxygen Evolution Reaction (OER), Alkaline Electrolysis, 3D Printing, Electrocatalysis

Considered like a fully renewable and clean energy carrier with the highest energy density, dihydrogen (H2) constitutes the best alternative to fossil fuels for insuring the sustainability of energy. This technologically and industrially relevant energy vector can be produced by water electrolysis that is now widely recognized as an eco-friendly, scalable and carbon-free route. In that context, alkaline water electrolysis is an affordable technology because it allows the use of transition metals as electrocatalysts for both hydrogen (HER) and oxygen (OER) evolution reactions. Here, combining catalytically active transition metal alloys (NiFe and NiCo) with a 3D printing technique, namely selective laser melting (SLM), enables access to ca. 25 cm 2 microstructured cylindrical electrodes efficiently promoting alkaline HER and OER. NiCo is found to be the most efficient electrode for HER, with an overpotential of 210 mV at 10 mA cm -2 . It is also the most stable electrode when studied in operation during prolonged electrolysis with a potential change of only 40 mV after 140 h electrolysis at 50 mA cm -2 (1.25 A). Regarding OER, NiFe shows the highest catalytic activity with an overpotential of 300 mV at 10 mA cm -2 and the greatest stability in operation with an electrolysis-induced potential change of only 30 mV. Further surface characterization techniques (scanning electron microscopy SEM, energy dispersive X-ray spectroscopy EDS, contact angle measurements, and analysis of electrogenerated gas bubbles) are used to monitor the electrolysis-induced chemical changes and to get valuable information on the bubble dynamics. Interestingly, electrolysis is found to be beneficial for enhancing the hydrophilicity/aerophobicity of the electrodes and thus facilitating the detachment of the H2 or O2 bubbles.

Introduction

In the context of an increase in worldwide energy consumption of about 60 % from 2010 to 2040 [START_REF] Lewis | Powering the Planet: Chemical Challenges in Solar Energy Utilization[END_REF] , it becomes crucial to develop large-scale economically viable technological processes to produce and store energy. Facing traditional fossil fuels responsible for climate change and global warming, dihydrogen H2 constitutes a clean energy carrier and sustainable alternative for meeting future energy demands [START_REF]The Future of Hydrogen[END_REF] . Nowadays, H2 can be used either as an energy vector owing to its high mass energy density (123 MJ kg -1 ), or as a chemical reagent for Haber-Bosch process [START_REF] Ghahraloud | Hydrogen Production through Thermal Decomposition of Hydrogen Sulfide: Modification of the Sulfur Recovery Unit To Produce Ultrapure Hydrogen[END_REF] . As an energy vector, it can be consumed into a fuel cell to produce electricity for mobility [START_REF] Acar | The Potential Role of Hydrogen as a Sustainable Transportation Fuel to Combat Global Warming[END_REF] , or to store the excess energy produced by intermittent and renewable energy sources [START_REF] Tourneur | Molecular and Material Engineering of Photocathodes Derivatized with Polyoxometalate-Supported {Mo3S4} HER Catalysts[END_REF] . Around 95% of the global H2 production still relies on methane steam reforming which leads to the emission of greenhouse gases such as carbon dioxide and carbon monoxide [START_REF]A National Vision of Amrica's Trasition to a Hydrogen Economy -to 2030 and Beyond[END_REF][START_REF]A Hydrogen Strategy for a Climate Neutral Europe[END_REF] . The generation of H2 by water electrolysis is now widely recognized as a clean and carbon-free route [START_REF] Terlouw | Large-Scale Hydrogen Production via Water Electrolysis: A Techno-Economic and Environmental Assessment[END_REF][START_REF] Kang | Mesoporous Multimetallic Nanospheres with Exposed Highly Entropic Alloy Sites[END_REF][START_REF] Jiang | Nanoarchitectonics of Metallene Materials for Electrocatalysis[END_REF] .

However, this process is still costly and cannot become affordable as long as it uses scarce and noble metals such as Pt, Pd, Ir, Ru as electrode materials [START_REF] Cao | Hierarchical NiMo Alloy Microtubes on Nickel Foam as an Efficient Electrocatalyst for Hydrogen Evolution Reaction[END_REF] . These metals are principally employed for PEM (Proton Exchange Membrane) electrolysis in acidic medium [START_REF] Lewis | Powering the Planet: Chemical Challenges in Solar Energy Utilization[END_REF][START_REF] Shiva Kumar | Hydrogen Production by PEM Water Electrolysis -A Review[END_REF] . This type of electrolyzer yields the highest current densities for H2 production but is prone to corrosion owing to the harsh oxidizing conditions in the anodic compartment. Research continues to explore alkaline water (AWE) and anion exchange membrane (AEM) electrolyzers because these technologies allow the use of transition metals as electrocatalysts for both hydrogen (HER) and oxygen (OER) evolution reactions [START_REF] Ehlers | Affordable Green Hydrogen from Alkaline Water Electrolysis: Key Research Needs from an Industrial Perspective[END_REF][START_REF] Yu | Clean and Affordable Hydrogen Fuel from Alkaline Water Splitting: Past, Recent Progress, and Future Prospects[END_REF] . Transition metals are generally widespread and quite cheap, in addition to offer a high electrocatalytic efficiency competing with that of Pt group metals [START_REF] Trasatti | Work Function, Electronegativity, and Electrochemical Behaviour of Metals[END_REF][START_REF] Roger | Earth-Abundant Catalysts for Electrochemical and Photoelectrochemical Water Splitting[END_REF] . A benchmarking study by Jaramillo and co-workers has enabled to evaluate and standardize the HER and OER electrocatalytic activity of numerous electrodeposited Ni-based binary alloys (incorporating Co, Cr, Mo, Fe, La or W) [START_REF] Mccrory | Benchmarking Hydrogen Evolving Reaction and Oxygen Evolving Reaction Electrocatalysts for Solar Water Splitting Devices[END_REF][START_REF] Mccrory | Benchmarking Heterogeneous Electrocatalysts for the Oxygen Evolution Reaction[END_REF] . Other groups have also shown that ternary alloys, such as NiCoFe [START_REF] Barati Darband | Facile Electrodeposition of Ternary Ni-Fe-Co Alloy Nanostructure as a Binder Free, Cost-Effective and Durable Electrocatalyst for High-Performance Overall Water Splitting[END_REF][START_REF] Bachvarov | Electrodeposited Ni-FeCo and NiFeCoP Alloy Cathodes for Hydrogen Evolution Reaction in Alkaline Medium[END_REF] , NiCoMo [START_REF] Jeghan | Mo-Incorporated Three-Dimensional Hierarchical Ternary Nickel-Cobalt-Molybdenum Layer Double Hydroxide for High-Efficiency Water Splitting[END_REF] , NiCrMo [START_REF] Nady | Novel Ni-Cr-Based Alloys as Hydrogen Fuel Sources through Alkaline Water Electrolytes[END_REF] , NiFeMo [START_REF] Hsieh | Alloy Inverse-Opals on Ni Foam as Outstanding Bifunctional Catalysts for Electrolytic Water Splitting of Ultra-Low Cell Voltages at High Current Densities[END_REF] and NiFeCr [START_REF] Yang | Highly Active Trimetallic NiFeCr Layered Double Hydroxide Electrocatalysts for Oxygen Evolution Reaction[END_REF] , could be excellent candidates for electrocatalyzing these two half-reactions of interest. In these alloys, each metal plays a key role with a positive impact on electrocatalysis. For instance, Mo shows a high affinity to adsorb protons while Ni can dissociate water [START_REF] Nairan | NiMo Solid Solution Nanowire Array Electrodes for Highly Efficient Hydrogen Evolution Reaction[END_REF] . The precise role of Fe is still in discussion, it could either act as an active site for the reaction, or activate it [START_REF] Bao | Suppressor, and Promoter-The Roles of Ni and Fe on Oxygen Evolution Reaction Activity and Stability of NiFe Alloy Thin Films in Alkaline Media[END_REF][START_REF] Feng | Fe-Based Electrocatalysts for Oxygen Evolution Reaction: Progress and Perspectives[END_REF] . Co alloyed with Ni yields a bifunctional material able to electrocatalyze both alkaline HER and OER. Concerning Cr, its presence enables the catalytic alloy to be effectively protected against corrosion [START_REF] Babu | Recent Developments on Cr-Based Electrocatalysts for the Oxygen Evolution Reaction in Alkaline Media[END_REF][START_REF] Peng | Stabilizing the Unstable: Chromium Coating on NiMo Electrode for Enhanced Stability in Intermittent Water Electrolysis[END_REF] . These catalytic materials are usually deposited on conducting (metal foam or carbon) substrates by either electrochemical or hydrothermal methods from the corresponding metal salts [START_REF] Sobha Jayakrishnan | Electrodeposition: The Versatile Technique for Nanomaterials[END_REF][START_REF] Feng | Hydrothermal and Solvothermal Syntheses[END_REF] . Although these simple and fast methods allow the catalyst structure to be controlled at the nanometer scale, their industrial scale-up is however not so easy.

In that context, 3D printing has appeared as an emerging technique to produce a new kind of metal components with complex geometries [START_REF] Ngo | Additive Manufacturing (3D Printing): A Review of Materials, Methods, Applications and Challenges[END_REF][START_REF] Debroy | Technological and Economic Issues in Metal Printing and Their Solutions[END_REF] . 3D printing offers numerous advantages such as the opportunities of preparing metallic alloys with customized composition, novel and complex electrode shapes and porous systems with high resolution and high control on the surface microstructure. Among the different 3D printing methods, the SLM (Selective Laser Melting) technique is particularly suitable to print different types of metals. The principle of SLM is based on the local melting of metal by a high-power laser beam scanning over a metal powder bed in a controlled manner (Figure 1). Thus, metal pieces with complex geometries can be fabricated layer-by-layer with a resolution between two layers of 20 µm. Some groups, such as Pumera et al [START_REF] Ambrosi | Helical 3D-Printed Metal Electrodes as Custom-Shaped 3D Platform for Electrochemical Devices[END_REF][START_REF] Ambrosi | Multimaterial 3D-Printed Water Electrolyzer with Earth-Abundant Electrodeposited Catalysts[END_REF][START_REF] Browne | 3D Printing for Electrochemical Energy Applications[END_REF] and Esposito et al [START_REF] Esposito | Membraneless Electrolyzers for Low-Cost Hydrogen Production in a Renewable Energy Future[END_REF][START_REF] Bui | 3D-Printed Electrodes for Membraneless Water Electrolysis[END_REF] , have already used this additive manufacturing technology to print electrodes or electrolyzers with original geometries. SLM has a current resolution of about 100 µm and can routinely printed 10 cm-sized metal pieces. This technology is already used for aeronautics industry and its development towards the energy sector remains still to be explored. Indeed, although the 3D printing methods hold great promises to fabricate practical catalytic electrodes [START_REF] Lee | 3D Printing for Electrocatalytic Applications[END_REF] , their use to elaborate metallic alloys specifically designed for electrochemical purposes has not been yet explored to the best of our knowledge [START_REF] Lee | 3D Printing for Electrocatalytic Applications[END_REF][START_REF] Liu | CFD-DEM Modelling of Mixing of Granular Materials in Multiple Jets Fluidized Beds[END_REF] .

In this work, ca. 25 cm 2 cylindrical NiFe and NiCo electrodes have been printed by SLM from the mixture of the corresponding metal powders. Their electrocatalytic performance for both HER and OER in alkaline medium has been examined and compared with that of reference Ni cylinder. Interestingly, it is also demonstrated that the introduction of cylindrical micropatterns on the electrode surface was an efficient strategy to facilitate the gas bubbles release. Moreover, complementary surface characterization techniques have been used to investigate the chemical changes of the electrodes caused by their prolonged cathodic or anodic polarization.

Experimental Section

Chemicals

Spherical powders for 3D printing of Ni (purity 99.99 %, 15-35 µm granulometry) was purchased from Hoganas. Spherical powders for 3D printing of Fe and Co (purity 99.99 %, 15-53 µm granulometry) were purchased from Neyco. KOH (>85 %, analysis grade) and absolute ethanol (RSE grade) were purchased from Sigma-Aldrich and Carlo Erbatron, respectively. Ultra-pure water (Veolia Water-STI) with a resistivity of 18.2 MΩ cm was used to prepare the aqueous solutions.

3D printing

All electrodes were designed using the software Auto Desk Fusion 360. Magics software was used for implementation on the platform using block and contour supports. A Turbula T2F station was used to mix the atomic powders in a powder bed of NiFe and NiCo. Electrodes were produced by SLS France company using the SLM technique in a Tech-Gine TS120 under nitrogen atmosphere. Each layer of melt metal was 20 µm-thick. The so fabricated electrodes were then removed from the support, scraped and sand, then cleaned under pressurized air. Photographs of these electrodes are shown in Figure S1 (Supporting Information).

Electrochemistry

The electrochemical measurements were performed at room temperature on a HTDS PARSTAT PMC-200 multichannel electrochemical workstation equipped with an AC frequency response analyzer (FRA). The 3D printed electrode was directly used as the working electrode (geometric surface area: 25 cm 2 ). A graphite rod or a double stabilized anode (DSA, model AO2002, titanium foil covered by mixed ruthenium and iridium oxides, ECS Tarn) was used as the counter electrode for OER or HER, respectively. The reference electrode was either Hg/HgO from Origalys for chronopotentiometry or KCl-saturated Ag/AgCl for cyclic voltammetry. For comparison, a pure Ni electrode (25 cm 2 geometric surface area) was considered. Before use, the working electrode was sonicated for 10 min in ultra-pure water and absolute ethanol to remove the oxide layer and was dried at 70°C for 3 h. An electrical contact was established with an ergot on the electrode and masked with an insulating epoxy resin (Loctite EA9460, Ellswoth Additives). The electrode was then placed into an oven at 90°C overnight to cure the resin.

A home-made H-type electrochemical macrocell with a constant headspace volume was used to carry out long-term electrolytic tests at high current intensities with ca. 25 cm 2 electrodes (Figure S2). The macrocell consisted of two compartments with a total electrolyte capacity of 1.4 L, in order to neglect the volume consumed during the water electrolysis. The headspace volume was 1.25 L per half-cell, so that volume changes caused by the produced gas (H2 and O2) could be neglected. This large headspace volume was chosen to avoid overpressure during GC analysis. The cell was sealed by a flat and circular ground joint, closed by a stainless steel flange. The junction was made by a male-female ground joint sealed with Teflon tape. Septa and epoxy resin (Loctite, EA3423, Ellswoth Additives) were used to seal the openings for gas measurement and introduction of the working electrode.

All cyclic voltammograms (CVs) were recorded at 5 mV s -1 and the Ohmic drop was corrected using the cell resistance determined at high frequencies (typically 100 kHz). The corresponding Tafel slopes were determined from the linear overpotential vs. log (current density) (η vs. log j) plots. The current densities were normalized by the geometric surface area of the 3D printed electrodes determined with Auto Desk Fusion 360 software and all the measured potentials were converted versus the reversible hydrogen electrode (RHE) according to the equations given below: In an alkaline medium (1.0 M KOH, pH = 14 ± 0.1) using a Hg/HgO reference electrode (filled with 0.1 M KOH solution):

ERHE = EHg/HgO + 0.059 pH + E 0 Hg/HgO with E 0 Hg/HgO = +0.14 V (1)
In an alkaline medium (1.0 M KOH, pH = 14 ± 0.1) using an Ag/AgCl reference electrode (filled with saturated KCl solution):

ERHE = EAg/AgCl + 0.059 pH + E 0 Ag/AgCl with E 0 Ag/AgCl = +0.20 V (2)
where EHg/HgO and EAg/AgCl are the applied potential values versus the selected reference electrode. Potentials are given with 5 mV precision and the reference electrodes were controlled every six months [START_REF] Kawashima | Accurate Potentials of Hg/HgO Electrodes: Practical Parameters for Reporting Alkaline Water Electrolysis Overpotentials[END_REF] . Electrochemical impedance spectroscopy (EIS) measurements were performed at an applied overpotential of 0.1 V for HER and 0.3 V for OER over a frequency range from 100 kHz to 100 mHz (50 frequencies, logarithmic distribution) at 10 mV amplitude. The long-term stability of the catalytic electrodes was examined using chronopotentiometry at an applied ±50 mA cm -2 current density. EIS data were analyzed using the Zview software.

Characterizations

Scanning electron microscopy (SEM) images and EDS (Energy Dispersive X-Ray Spectroscopy) elemental mapping were performed using a IT-300 JEOL microscope and an X-Max 50 mm 2 EDS detector (Oxford Instruments).

Profilometry measurements were performed using a Bruker Dektak X profilometer with a Stylus B type tip (12.5 µm curvature radius) and tilt angle correction.

ICP-AES (Inductive Coupled Plasma Atomic Emission

Spectroscopy) measurements were performed using a 725 ES Varian station with an argon plasma.

Static water contact angles were measured with an automated Krüss Easy Drop goniometer. At least three small droplets of 2.0 μL of ultra-pure water were dispensed, and the contact angles were determined using a Tangent 2 fitting model.

Gas chromatography (GC) measurements for the detection and quantification of the gases produced by electrolysis were performed with a micro-GC Agilent 990 equipped with a TCD detector. The column was heated at 145°C, the injector at 90°C, with Alphagaz 2 Smartop argon (Air Liquide) as the carrier gas at 29 psi. Sampling time was 50 ms and the results are obtained in vol. %. For the determination of the Faradaic efficiency (FE), the apparatus was calibrated with a gas mixture containing N2 (3%, v/v), CH4 (5%, v/v), CO (10%, v/v), H2 (10%, v/v), CO2 (20%, v/v) and He (52%, v/v) (from GazDetect, UN1954) and air from the laboratory. A H-type cell perfectly closed was used and degassed under Ar for 10 min prior to analysis. A current density of ±10 mA cm -2 was applied during 1.5 h and one measurement was realized every 10 min. The current density was selected considering the theoretical gas production to not exceed 10% of the headspace volume of the cell. The FE was determined using the gas perfect law and the pressure value measured by a Vernier sensor. It was expressed versus a reference material (Pt for HER and DSA for OER) studied under the same conditions. We have checked that a 1 cm 2 Pt sheet or DSA grid immersed in ca. 100 mL electrochemical cell and polarized at ± 10 mA yielded FEs of 100% for HER and OER, respectively.

Bubble dynamics

Pictures of the electrode covered with hydrogen bubbles were acquired at 35 frames per second with a ThorLabs CS165MU/M camera equipped with a MVL35M23 objective. Screenshots from the videos were processed using the freeware Image J to estimate the average surface coverage and the size of the bubbles.

Results and Discussion

Characterization of 3D printed binary alloys

The morphological characteristics and composition of 25 cm 2 3D printed materials were evaluated by profilometry, SEM and EDS. First, profilometry measurements revealed rough surfaces with an average roughness of 5.0 ± 1.3 µm irrespective of the composition of the binary alloy. It is worth emphasizing that high roughness is a main characteristic of the SLM technique [START_REF] Koutiri | Influence of SLM Process Parameters on the Surface Finish, Porosity Rate and Fatigue Behavior of as-Built Inconel 625 Parts[END_REF] . Dimensions of the cylindrical patterns engraved on the electrode surface ranged from 200 to 300 µm for their height, 700 to 800 µm and 400 to 500 µm for their vertical and horizontal elongations, respectively (Figure S3). The chemical composition and its homogeneity were determined by SEM-EDS at the micrometer scale (Figure 2). NiCo showed a homogeneous distribution of Ni and Co elements (Figure S4). The Co content was estimated at 10% from EDS. In contrast, Ni and Fe were not homogeneously distributed in the NiFe alloy and two phases were clearly evidenced: i) one phase showing Fe microparticles of 20-30 µm diameter (red large areas in Figure 2) and ii) a second one characterized by µm-sized Fe particles in a Ni matrix (red dots in the green area). The average Fe content was around 6%.

3.2

Hydrogen Evolution Reaction

Electrocatalytic activity and stability

The electrocatalytic activity of the 25 cm 2 3D printed electrodes was evaluated in alkaline medium (1.0 M KOH, pH 14) both for HER and OER. Based on the electrochemical surface area of these electrodes, it is worth emphasizing that current intensities as high as 1 A could be delivered.

Before the HER process, a reversible redox wave corresponding to the formation and oxidation of nickel hydrides [START_REF] Hall | The Electrochemistry of Metallic Nickel: Oxides, Hydroxides, Hydrides and Alkaline Hydrogen Evolution[END_REF] was observed within the potential range 0.0-0.3 V vs RHE for the three electrodes (Figure S5). This wave was followed by the HER one starting below ca. -0.1 V vs RHE (Figure 3a). NiCo exhibited the highest HER activity, as supported by the lowest iR-corrected overpotential of 210 mV measured at 10 mA cm -2 (250 mA) (Table 1). This value was 20 and 60 mV lower than those determined for Ni and NiFe, respectively. Comparable Tafel slopes between 89 (Ni) and 103 mV per decade (NiFe) were determined, in relatively close agreement with those reported for Ni-based HER cathodes studied in alkaline medium [START_REF] Hsieh | Alloy Inverse-Opals on Ni Foam as Outstanding Bifunctional Catalysts for Electrolytic Water Splitting of Ultra-Low Cell Voltages at High Current Densities[END_REF][START_REF] Nairan | NiMo Solid Solution Nanowire Array Electrodes for Highly Efficient Hydrogen Evolution Reaction[END_REF][START_REF] Ambrosi | Multimaterial 3D-Printed Water Electrolyzer with Earth-Abundant Electrodeposited Catalysts[END_REF][START_REF] Zhang | NiFe-Based Nanostructures on Nickel Foam as Highly Efficiently Electrocatalysts for Oxygen and Hydrogen Evolution Reactions[END_REF] . a HER overpotential at -10 mA cm -2 . b After 140 h electrolysis at -50 mA cm -2 . c At -0.1 V vs RHE.

In order to gain further insights on the electron transfer kinetics, EIS measurements were performed at an applied overpotential of 0.1 V and in the frequency range 10 5 -10 -1 Hz. The data are presented as plots of the imaginary part (Z″) vs the real part (Z′) of the complex impedance (Z = Z′ + iZ″) as a function of frequency f (Nyquist diagram). All Nyquist plots are characterized by a semicircle and can be reasonably fitted by assuming the equivalent electrical circuit comprising a resistance (the cell resistance Rs) in series with the charge transfer resistance (Rct) and a constant phase element (CPE) in parallel (Figure 4a). For this circuit model, the diameter of the semicircle on the Z′ axis corresponds to Rct. Moreover, CPE is responsible for the deformation of the semi-circle in the diagram due to a frequency dispersion of the capacitance and accounts for surface inhomogeneities such as roughness or structuration [START_REF] Pajkossy | Impedance Spectroscopy at Interfaces of Metals and Aqueous Solutions -Surface Roughness, CPE and Related Issues[END_REF][START_REF] Jorcin | CPE Analysis by Local Electrochemical Impedance Spectroscopy[END_REF][START_REF] Hirschorn | Determination of Effective Capacitance and Film Thickness from Constant-Phase-Element Parameters[END_REF] . Consistent with CV data, the Rct of the NiCo electrode (215 mΩ cm -2 ) was found to be smaller than that calculated for Ni (280 mΩ cm -2 ) and NiFe (454 mΩ cm -2 ) demonstrating faster HER charge transfer kinetics (Table 1). The long-term stability of 3D printed cathodes was examined at a constant current density of -50 mA cm -2 over six days. The corresponding chronopotentiometric curves showed that NiCo was stabilized at a less negative potential (namely -0.33 V vs RHE) with respect to the two other electrodes (Figure 4b). It was also the most robust electrode as the potential change was only 40 mV after 140 h electrolysis, while a potential change of 100 mV was monitored for Ni and NiFe.

FE for HER was about 90% vs reference Pt for the three electrodes (Figure S6 and Table 1). This value lower than 100% could be reasonably explained by the large-scale electrodes and the large headspace volume of the cell leading to unavoidable loss of produced gas. Moreover, it is noticeworthy that the temperature of the electrolytic solution increased by 7°C after electrolysis which could also account for a non-quantitative FE for HER.

Bubble dynamics

Among the factors influencing the yield of an electrolyzer, the efficiency of the release of H2 and O2 bubbles produced at the cathode and anode, respectively, represents a critical issue for high-performance electrocatalysis. This problem is often minimized in laboratory-scale electrochemical cells by using small electrode surfaces of less than 1 cm 2 . This constitutes however a major issue in industry, as bubble films formed at large scale (a few dozen to a few hundred tens cm 2 ) electrodes can be more than 1 cm thick, leading to high Ohmic resistances and mass transfer barriers of the electrode reactions, and consequently resulting in a strong enhancement in the electrical power consumption of the electrolyzer [START_REF] Kempler | Effects of Bubbles on the Electrochemical Behavior of Hydrogen-Evolving Si Microwire Arrays Oriented against Gravity[END_REF][START_REF] Zhao | Gas Bubbles in Electrochemical Gas Evolution Reactions[END_REF] . Consequently, the control of both the nucleation and the flow of gas bubbles is a crucial issue to optimize electrolysis process and reduce costs. In that context, the N. S. Lewis' group produced some valuable data about this topic [START_REF] Kempler | Effects of Bubbles on the Electrochemical Behavior of Hydrogen-Evolving Si Microwire Arrays Oriented against Gravity[END_REF][START_REF] Kempler | Optical and Electrochemical Effects of H 2 and O 2 Bubbles at Upward-Facing Si Photoelectrodes[END_REF] . The electrode microstructuring, in the form of e.g. conic arrays [START_REF] Ma | Directional and Continuous Transport of Gas Bubbles on Superaerophilic Geometry-Gradient Surfaces in Aqueous Environments[END_REF] , canals [START_REF] Bintein | Grooves Accelerate Dew Shedding[END_REF] , or microwires [START_REF] Kempler | Optical and Electrochemical Effects of H 2 and O 2 Bubbles at Upward-Facing Si Photoelectrodes[END_REF] , and the surface functionalization [START_REF] Xie | Engineered Switchable-Wettability Surfaces for Multi-Path Directional Transportation of Droplets and Subaqueous Bubbles[END_REF][START_REF] Bae | Superaerophobic Polyethyleneimine Hydrogels for Improving Electrochemical Hydrogen Production by Promoting Bubble Detachment[END_REF] have been demonstrated to be promising approaches to control the aerophobicity of the surface and thus the bubble release [START_REF] Xu | High-Performance Water Electrolysis System with Double Nanostructured Superaerophobic Electrodes[END_REF] .

In this study, photographs have been captured from videos recorded during the electrolysis of 3D printed cathodes. Two current densities were selected, namely -1 and -30 mA cm -2 , in order to observe the nucleation phase of bubbles and their detachment from the electrode surface when a high density of bubbles was electrogenerated. Effects of the long-term electrolysis on the efficiency of the bubbles release were also investigated. Globally, the characteristics of H2 bubbles electrogenerated at -1 mA cm -2 on freshly prepared electrodes were not dependent on the nature of the electrode (Figures 5 andS7). The electrode surface was covered by about 5% of H2 bubbles and the average radius of one bubble was 55 ± 15 µm (determined from the image processing by the ImageJ freeware, Figure S8). After a longterm electrolysis at -50 mA cm -2 for 140 h, the photographs showed a vertically directed bubble cloud with few bubbles remaining on the electrode surface. This apparent facilitated release of H2 bubbles caused by electrolysis was thought to be ascribed to the enhanced hydrophilicity/aerophobicity of the electrodes. Indeed, the contact angles measured with KOH 1M for the three electrodes were dramatically decreased from ca. 110° to 30° ± 5° after electrolysis. A plausible reason for the enhanced aerophobicity of the electrodes could be the formation of nickel hydrides due to the hydrogen diffusion into Ni, which is commonly observed for Ni electrodes studied for HER in alkaline medium [START_REF] Hall | The Electrochemistry of Metallic Nickel: Oxides, Hydroxides, Hydrides and Alkaline Hydrogen Evolution[END_REF] . At -30 mA cm -2 , approximately 80% of the micropatterns acted as coalescence sites of H2 bubbles (Figure S9, movie 1). The mean bubble departure radius was estimated to be 120 ± 30 µm. After electrolysis, only 10-15% of the patterns were still covered by bubbles with unchanged departure radius (movie 2).

Post-electrolysis characterizations

Different surface characterization techniques have been used to monitor the potential morphology and composition changes of the electrodes caused by 140 h electrolysis at -50 mA cm -2 . First, post-electrolysis profilometry analysis confirmed that cylindrical micropatterns geometry and dimensions remained unchanged. Some lower atomic number phases were now visible on NiCo and NiFe to a lesser extent which were evidenced by the darker areas on the backscattered electron SEM images (Figure 6). The EDS analysis showed that the metallic composition of the three electrodes was not significantly altered after electrolysis (Table S1). It was however noticed an increase in the O content and the presence of potassium varying from 2 at. % (Ni and NiFe) to 15% (NiCo). These results are consistent with the adsorption of potassium carbonates on the electrode surface and are not surprising in sofar as such species were systematically detected for Ni-based HER and OER catalysts polarized in alkaline KOH medium [START_REF] Mccrory | Benchmarking Hydrogen Evolving Reaction and Oxygen Evolving Reaction Electrocatalysts for Solar Water Splitting Devices[END_REF] . Although all the detected elements were homogeneously distributed on the electrode surface in the case of Ni and NiFe, two distinct phases could be distinguished for NiCo (Figures S10-S12). The Ni and Co contents of the majority phase were relatively close to those before electrolysis while the second one which was poor in metals and rich in K and O was probably composed of potassium carbonates. It is obvious that these two chemically different phases were responsible for the bright and dark areas, respectively, in the backscattered electron SEM images (Figure 6).

Finally, the ICP-AES analysis of the electrolytic solution showed the absence of released metals after electrolysis.

3.3

Oxygen Evolution Reaction

Electrocatalytic activity and stability

The same protocol as that used for HER was followed for measuring the electrocatalytic activity of the three 3D printed electrodes for OER. Before OER, a scan rate-dependent reversible redox wave (Ni 3+ /Ni 2+ ) was first observed for all anodes between 1.2 and 1.5 V vs RHE, corresponding to the oxidation of Ni hydroxide to NiOOH [START_REF] Hall | The Electrochemistry of Metallic Nickel: Oxides, Hydroxides, Hydrides and Alkaline Hydrogen Evolution[END_REF] which is the catalytically active species for OER (Figure S5) [START_REF] Trotochaud | Solution-Cast Metal Oxide Thin Film Electrocatalysts for Oxygen Evolution[END_REF][START_REF] Bediako | Structure-Activity Correlations in a Nickel-Borate Oxygen Evolution Catalyst[END_REF] . An additional anodic peak was also visible at ca. 1.2 V for NiCo which could be also ascribed to the conversion of Co hydroxide to its oxyhydroxide form [START_REF] Sivanantham | Hierarchical NiCo2S4 Nanowire Arrays Supported on Ni Foam: An Efficient and Durable Bifunctional Electrocatalyst for Oxygen and Hydrogen Evolution Reactions[END_REF] . This transition metal hydroxide/oxyhydroxide redox process was followed by the OER signal displaying a much larger current density (Figure 7a). NiFe showed an overpotential for OER of 300 mV at 10 mA cm -2 which was lower than those measured for Ni and NiCo (~350 mV) (Table 2). The higher efficiency of NiFe for OER was confirmed by its Tafel slope of 33 mV per decade (Figure 7b) which was smaller than those of NiCo (40 mV per decade) and Ni (63 mV per decade). It is worth noting that such values were in perfect line with those determined for other Ni-based OER electrocatalysts deposited on nickel foam or stainless steel [START_REF] Ambrosi | Multimaterial 3D-Printed Water Electrolyzer with Earth-Abundant Electrodeposited Catalysts[END_REF][START_REF] Zhang | NiFe-Based Nanostructures on Nickel Foam as Highly Efficiently Electrocatalysts for Oxygen and Hydrogen Evolution Reactions[END_REF] . Faster OER charge transfer kinetics of NiFe was also assessed by the Rct value of 6.0 mΩ cm -2 extracted from Nyquist plots which was much smaller than those determined for Ni (20.0 mΩ cm -2 ) and NiCo (87.2 mΩ cm -2 ) (Figure 8a). Additionally, NiFe showed a relatively good catalytic stability under operation. Long-term electrolysis at +50 mA cm -2 for 140 h revealed that NiFe was the most robust anode as the potential change was only 30 mV after 140 h electrolysis, while potential changes of 90 and 100 mV were monitored for Ni and NiCo, respectively (Figure 8b). FE for OER was 100% vs reference DSA for the three electrodes (Figure S13). Nevertheless, it can be noticed that the FE of DSA was only 83% vs the theoretical FE calculated from the consumed electrical charge. This discrepancy can be reasonably explained by the large headspace volume and the increase in the electrolyte temperature during electrolysis, as mentioned before for HER.

Bubble dynamics

Globally, the surface coverage of O2 bubbles electrogenerated at +1 mA cm -2 was lower than that of H2 bubbles (Figures 9 and S14) which is quite logical given that two moles of H2 are produced for one mole of O2 in the overall water splitting. Furthermore, the mean bubble nucleation and departure radii were not significantly different from those estimated for HER (Figure S15). After a long-term electrolysis at +50 mA cm -2 for 140 h, the bubbles were more easily released from the electrode surface, as a result of the increased hydrophilicity/aerophobicity of the electrodes. Indeed, similar to HER, the contact angles measured with KOH 1M for the three electrodes were decreased from ca. 110° to 30° ± 5° after electrolysis. Such a trend can be explained by the formation of metal hydroxides/oxyhydroxides during the electrolysis [START_REF] Wei | Insight into the Boosted Electrocatalytic Oxygen Evolution Performance of Highly Hydrophilic Nickel-Iron Hydroxide[END_REF] . At +30 mA cm -2 , very few bubbles were observed on the electrode surface either before or after electrolysis (Figure S16, movies 3 and 4). Moreover, a careful examination of bubbles photographs showed that the flow of released O2 bubbles was higher for NiFe and NiCo than for Ni at both +1 and +30 mA cm -2 .

Post-electrolysis characterizations

As for HER, the profilometry analysis of the three electrodes after 140 h electrolysis at +50 mA cm -2 showed no significant changes in the engraved micropatterns geometry and dimensions. A slight decrease in the content of Ni was observed for Ni and NiFe whereas the metallic composition of NiCo remained stable (Table S2).

Conclusion

To sum-up, ca. 25 cm 2 cylindrical Ni, NiFe and NiCo electrodes were fabricated by 3D printing technique and their electrocatalytic performance for alkaline HER and OER was evaluated. Nearly cylindrical micropillars engraved on the electrode surface acted as preferential sites for the nucleation of the electrogenerated gas bubbles. Regarding HER, NiCo was found to be the most efficient and the most stable electrode when studied in operation during electrolysis at -50 mA cm -2 (-1.25 A) over 140 h. Compared with Ni and NiFe, this result was somewhat surprising if one considers that the chemical composition of this material was the most strongly modified after electrolysis with a higher proportion of adsorbed potassium carbonates in particular. For the other half-reaction, NiFe showed a superior OER activity and stability in operation. Furthermore, the chemical composition of the bimetallic alloys was not significantly changed after 140 h electrolysis at +50 mA cm -2 . Interestingly, electrolysis was found to be beneficial for enhancing the hydrophilicity/aerophobicity of the electrodes and thus facilitating the detachment of the H2 or O2 bubbles.

This work constitutes a first promising step towards practical and industrially relevant functional electrodes able to deliver current intensities higher than 1 A for electrolytic water splitting. By taking the advantage of the versatility of the SLM technique, other Ni-based bi-and trimetallic alloys with different geometry and structuring will be envisaged in order to enhance the catalytic performance and stability in operation. In that context, the resolution of the engraved micropatterns is still low and needs to be increased in the future to reach the micrometer scale. This is a necessary step to finely control the flow and the detachment of the electrogenerated gas bubbles, and consequently the efficiency of the electrolyzer. Work in this stimulating direction is currently in progress.

ASSOCIATED CONTENT

Supporting Information. Supplementary profilometry curves, SEM-EDS mapping, photographs and analysis of the electrogenerated gas bubbles, GC chromatograms and additional cyclic voltammograms. This material is available free of charge via the Internet at http://pubs.acs.org. ANRT (Agence Nationale de la Recherche Technologique) and H2X Ecosystems are acknowledged to fund this work through the CIFRE program.

Notes

Due to the industrial work, the SLM process laser parameters and Auto Desk Fusion instructions for the patterning of the electrodes are protected and under confidentiality for submitting further patent.

ACKNOWLEDGMENT

Muriel Escadeillas (CRMPO, Univ. Rennes) is fully acknowledged for ICP-AES measurements performed at CRMPO platform. Axel Dubois and Hugo Schaal (CEO, SLS France) are fully acknowledged for their contribution to the fabrication of the 3D printed electrodes.

ABBREVIATIONS

AEM, Alkaline Exchange Membrane; AWE, Alkaline Water Electrolysis; CPE, Constant Phase Element; DSA, Double Stabilized Anode; EDS, Energy Dispersive X-Ray Spectroscopy; EIS, Electrochemical Impedance Spectroscopy; FE, Faradaic Efficiency; HER, Hydrogen Evolution Reaction; ICP-AES, Inductive Coupled Plasma Atomic Emission Spectroscopy; KPI, Key Performance Indicator; OER, Oxygen Evolution Reaction; PEM, Proton exchange Membrane; SEM, Scanning Electron Microscopy; SLM, Selective Laser Melting.

Figure 1 .

 1 Figure 1. Schematic illustration of the metal additive manufacturing with a metal powder bed employing SLM (Selective Laser Melting) technology. Reproduced with permission from reference 32. Copyright 2018 Elsevier.

Figure 2 .

 2 Figure 2. EDS elemental mapping of NiFe showing the inhomogeneous distribution of Ni (blue) and Fe (red).

Figure 3 .

 3 Figure 3. a) Ohmic drop-corrected cyclic voltammograms of Ni, NiFe and NiCo at 5 mV s -1 in KOH 1M. b) Corresponding Tafel plots.

Figure 4 .

 4 Figure 4. a) Nyquist plots at -0.1 V vs RHE for Ni, NiFe and NiCo from 100 kHz to 100 mHz and equivalent electrical circuit used for fitting (inset). b) Chronopotentiometric curves at -50 mA cm -2 for Ni, NiFe and NiCo in KOH 1M.

Figure 5 .

 5 Figure 5. Photographs of the best performing NiCo electrode for HER at -1 mA cm -2 in KOH 1M before (a) and after (b) longterm electrolysis at -50 mA cm -2 for 140 h.

Figure 6 .

 6 Figure 6. Large-scale backscattered electron SEM images of Ni (a), NiFe (b) and NiCo (c) after electrolysis at -50 mA cm -2 for 140 h.

Figure 7 .

 7 Figure 7. a) Ohmic drop-corrected cyclic voltammograms of Ni, NiFe, and NiCo at 5 mV s -1 in KOH 1M. b) Corresponding Tafel plots.

Figure 8 .

 8 Figure 8. a) Nyquist plots at 1.5 V vs RHE for Ni, NiFe and NiCo from 100 kHz to 100 mHz. b) Chronopotentiometric curves at +50 mA cm -2 for Ni, NiFe and NiCo in KOH 1M.

Figure 9 .

 9 Figure 9. Photographs of the best performing NiFe electrode for OER at +1 mA cm -2 in KOH 1M before (a) and after (b) longterm electrolysis at +50 mA cm -2 for 140 h.

  The analysis by ICP-AES of the electrolytic solution after electrolysis revealed that metallic Ni was released from the Ni and NiFe electrodes at the level of traces, i.e 60-200 ppm. Consistent with these findings, chemically different phases were clearly evidenced by backscattered electron SEM for Ni (Figure10). The darker areas corresponded to phases richer in O and poorer in K in the EDS map (FigureS17) and therefore could be reasonably ascribed to Ni oxides. In contrast, all the globally homogeneously distributed on the NiFe and NiCo surfaces (FiguresS18 and S19).

Figure 10 .

 10 Figure 10. Large-scale backscattered electron SEM images of Ni (a), NiFe (b) and NiCo (c) after electrolysis at +50 mA cm -2 for 140 h.

Table 1 .

 1 Key Performance Indicators (KPIs) Measured for Ni, NiFe and NiCo for Hydrogen Evolution Reaction.

	Cathode	Ni	NiFe NiCo
	η -10 / mV a	230	270	210
	Tafel slope / mV dec -1	89	103	97
	Faradaic efficiency / %	88	92	90
	Electrolysis-induced potential	100	100	40
	change / mV b
	Rs / mΩ cm -2 c	7.2	7.2	9.2
	Rct / mΩ cm -2 c	280	454	215

Table 2 .

 2 KPIs Measured for Ni, NiFe and NiCo for Oxygen Evolution Reaction.

	Anode	Ni	NiFe NiCo
	η +10 / mV a	350	300	360
	Tafel slope / mV dec -1	63	33	40
	Faradaic efficiency / %	100	100	100
	Electrolysis-induced potential	90	30	100
	change / mV b
	Rs / mΩ cm -2 c	6.0	7.2	7.2
	Rct / mΩ cm -2 c	20.0	6.0	87.2

a OER overpotential at +10 mA cm -2 . b After 140 h electrolysis at +50 mA cm -2 . c At 1.5 V vs RHE.

Funding Sources

AUTHOR INFORMATION

Corresponding Author * Loïc Perrin: loic.perrin@h2x-ecosystems.com * Bruno Fabre: bruno.fabre@univ-rennes.fr

Author Contributions