\documentclass[review,3p]{elsarticle}
\usepackage{graphicx}\usepackage{color}
\usepackage{subfigure}\usepackage{epstopdf}\usepackage{upgreek}
\usepackage{amssymb,amsmath}
\usepackage{booktabs}
\usepackage{mathtools}

\journal{Materials Science and Technology}

\bibliographystyle{elsarticle-num}
%\bibliographystyle{elsarticle-harv}
%\bibliographystyle{elsarticle-num-names}
\biboptions{sort&compress}

\begin{document}

\begin{frontmatter}
\title{Hydrogen and the tetragonality of FeNi (tetrataenite)}

 \author[a]{Jae Hoon Jang} \author[b]{H. K. D. H. Bhadeshia}
\address[a]{Korea Institute of Materials Science,  Changwon, 51508, Republic of Korea}
\address[b]{Queen Mary University of London,  School of Engineering and Materials Science, United Kingdom}

\date{\today}

%\begin{graphicalabstract}
%\includegraphics[width=1.0\textwidth]{abstract.jpg}
%\end{graphicalabstract}

\begin{abstract}
Calculations have been conducted to show that it may be feasible to enhance the tetragonality ($c/a$ ratio) of tetrataenite ($L1_{0}, P4/mmm$) by charging it with hydrogen, a process that can be achieved at ambient temperature, well below the ordering temperature of the FeNi compound. Although the saturation magnetisation is slightly compromised, the magnetocrystalline anisotropy should increase significantly. It is now worth doing experiments to establish or disprove the results, because the energy of inserting a single hydrogen atom into the location that enhances tetragonality is slightly greater than that in a sub-optimal site in the context of increasing the $c/a$ ratio. Changes in tetragonality as a function of the location of the hydrogen atom can be explained in terms of the relative strengths of the Fe-H and Ni-H bonds.
\end{abstract}

\begin{keyword} Tetrataenite \sep ordered FeNi \sep tetragonality \sep magnetic coercivity \sep hydrogen additions \sep first principles calculations

\end{keyword}

\end{frontmatter}

\section{Tetrataenite}
%:
The metal-rich phases found in iron, stony iron, and  as inclusions in chondrite  meteorites, consist of $\upalpha$ (A2; $Im\overline 3m$), $\upgamma$ (A1; $Fm\overline 3m$), $\upgamma''$ (L1$_0$, $P4/mmm$) and $\upgamma'$ (L1$_2$; $Pm\overline 3m$) \cite{Albertsen:1978,Clarke:1980}. In the absence of impacts, the meteorites cool from the red-hot state ($>1500$\,K) at 1-1000\,K\,per million years \cite{Lyons:2019}. All the phases are ferromagnetic below their respective Curie temperatures, but tetrataenite is particularly interesting because of its high magnetocrystalline anisotropy~\cite{pauleve1968magnetization}. High magnetocrysalline anisotropy can primarily be considered to enhance coercivity, which also is influenced by the microstructure~\cite{coey2010magnetism}. In the perfectly ordered state, the magnetic moment is 1.625\,$\mu_{\rm B}\,{\rm atom ^{-1}}$, with iron contributing 2.6-2.63\,$\mu_{\rm B}$ and nickel 0.6-0.62\,$\mu_{\rm B}$ \cite{Guenzburger:2005,Lewis:2014}. The moment per iron atom is greater than that found for any of the three common allotropic forms of pure iron \cite[pp.30-39,][]{Bhadeshia:2021},\cite{Steinle:2004}. The first-principles calculated (0\,K) direction of magnetisation of a single crystal of tetrataenite is along the four-fold axis of symmetry  [001] \cite{Guenzburger:2005}.

Tetrataenite has the chemical composition close to FeNi and orders at  $T_{\rm ord}=321\pm 2$\,$^\circ$C \cite{Pauleve:1962,Swartzendruber:1991} where diffusion is imperceptibly slow. It is therefore impossible to replicate the compound found in meteorites, on Earth, though valiant efforts have produced samples by difficult means and in an imperfectly ordered state in the form of minute particles~\cite{scorzelli1997study,lima2001new}. It is likely, however, that progress may be made in the future so that permanent magnets can be made that don't quite compete with rare-earth varieties but have a favourable position in the cost versus $(-BH)_{\rm max}$ plot \cite{Mandal:2023}.

A low-symmetry crystal structure that is ferromagnetic will tend to have greater magnetocrystalline anisotropy. The tetragonaility of $\upgamma''$ as a face-centred tetragonal structure is rather small, with lattice parameters of $a = 3.568\pm0.0018$\,\AA\  and $c = 3.5901 \pm0.0015$\,\AA\ as measured using synchrotron X-rays [15, 16]. However, the primitive cell of $\upgamma''$  is body-centred tetragonal, with $a_{\rm bct} = a/\sqrt 2$ and  $c/a_{\rm bct} \approx \sqrt2$, resulting in a large magnetocrystalline anisotropy \cite{woodgate2023revisiting}. If something could be done to enhance the tetragonality, that would inevitably lead to a greater magnetocrystalline anisotropy. The best way to do this might be to introduce small but misfitting-atoms at $T<T_{\rm ord}$, into the irregular-octahedral interstices that are slightly longer along the $c$-axis.\footnote{The interstices are not strictly octahedral, but the term {\em irregular} is henceforth omitted for brevity. The same applies to the tetrahedral interstices.}   Single-crystal moduli $C_{ij}$ ($P4/mmm$) have been calculated as a function of the long-range order parameter $\eta$ \cite{Tian:2019}, from which the compliances $S_{ij}$ have been calculated here in a manner consistent with the $4/mmm$ tetragonal symmetry, as in \cite{Ballato:1995}:

\begin{center}
\begin{tabular}{cccccccccccccc}
GPa: & $C_{11}$ & $C_{12}$ & $C_{13}$ & $C_{33}$ & $C_{44}$ & $C_{66}$ \\[3pt]
$\eta=1$ & 226 & 173 & 148 & 268 & 131 & 151  \\[3pt]
$\eta=0$ & 207 & 168 & & & 119 & \\[3pt] \hline \\[3pt]
1/GPa: & $S_{11}$ & $S_{12}$ & $S_{13}$ & $S_{33}$ & $S_{44}$ & $S_{66}$ \\[3pt]
$\eta=1$& 0.011504 & 0.007246 & 0.002351 & 0.006324 & 0.007647 & 0.006608\\[3pt]
$\eta=0$ & 0.017483 & -0.007820 & & & 0.008430 & \\[3pt]\
\end{tabular}\end{center}

Using the data, it is straightforward to show  that the $\upgamma''$ single crystal is much stiffer along [001]$_{\upgamma''}$ than along $\langle100\rangle_{\upgamma''}$.\footnote{The modulus along a specific direction is given by $E^{-1}\{{\bf n}\}=S_{11}(n_1^4+n_2^4) +S_{33}n_3^4 +(2S_{13}+S_{44})(1-n_3^2)n_3^2+ (2S_{12}+S_{66})n_1^2n_2^2, $ \cite{Gorodtsov:2021}. The $n_i$ are the components of a unit vector {\bf n}. Using this, for $\upgamma''$, $E_{100}=E_{010}=87$\,GPa, whereas $E_{001}=158$\,GPa.} Substitution of misfitting atoms into the octahedral interstices should then lead to a reduction in the $c/a$ ratio. This is confirmed by first-principles calculations of tetrataenite that has been nitrided \cite{Matar:2019}, perhaps leading to the general conclusion that any small atom that causes dilation  would reduce the tetragonality. An elasticity analysis does not account for subtleties based on the local environment around the any interstitial atoms in an ordered $\upgamma''$ as discussed below in the context of hydrogen alloying.

The alloying of tetrataenite with hydrogen has never previously been examined, but the huge advantage is that it can be incorporated into $\upgamma''$  by electrochemical charging under {\em ambient} conditions, i.e. well below the ordering temperature. Hydrogen is known to expand ordinary austenite \cite{Ulmer:1993}, i.e., it's effective size is greater than that of the interstice, and therefore might be expected to reduce the $\upgamma''$ $c/a$ ratio given that the $a$ parameter is smaller than $c$.   The conclusion reached on the basis of $S_{ij}$ that the tetragonality of $\upgamma''$ should decrease when hydrogen is introduced into the lattice neglects H-metal interactions that are dependent on the specific environment around each hydrogen atom, which will not be identical for similarly-shaped interstices in $\upgamma''$-$P4/mmm$. Assuming that there is an iron atom at the origin of the unit cell, an octahedral interstice at $[ \frac{1}{2}\, 0 \,0]$ has nickel atoms at $[ \frac{1}{2}\, 0\,\pm \frac{1}{2}]$, whereas the corresponding interstice at $[ \frac{1}{2}\, \frac{1}{2}\, \frac{1}{2}]$ has iron atoms at $[ \frac{1}{2}\, \frac{1}{2}\,0]$ and $[ \frac{1}{2}\, \frac{1}{2}\,1]$. This kind of subtlety is absent in the elastic analysis presented above.  It is necessary therefore to conduct first-principles calculations to study  changes in the $c/a$ that could in principle promote magnetocrystalline anisotropy. Insertion of hydrogen into tetrahedral interstices, which are less favoured for hydrogen occupation in paramagnetic austenite in pure iron  \cite{Song:2013b}, has also been examined out of interest.


\section{Calculations}
First-principles calculations were conducted using  the Vienna Ab-initio Simulation Package \cite{Hafner:2008,Kresse:1994}  with the projector augmented-wave formalism. The $\upgamma''$ ($P4/mmm$) was implemented using a conventional unit cell with two Fe atoms and two Ni atoms, as shown in Fig.~\ref{fig01} which includes the locations of the octahedral-interstices. Unlike a disordered structure, the two kinds of octahedral interstices at $z=0,1$ and $z=\frac{1}{2}$  are not equivalent. Those at $z=0,1$ are surrounded on the plane normal to $z$ by four near-neighbour iron atoms, whereas those  at $z=\frac{1}{2}$ have four near-neighbour nickel atoms. The two sets of interstices are, therefore, not equivalent, and the consequences will become apparent in the calculation results. The tetrahedral interstices are located at all equivalent  [1/4, 1/4, 1/4] positions, with each such interstice surrounded by two Fe atoms and two Ni atoms.
\begin{figure}[ht!]
\centering
\subfigure[]{
\includegraphics[width=0.45\textwidth]{tetrataenite_structure.eps}}\\
\subfigure[]{
\includegraphics[width=0.75\textwidth]{sites.eps}}
\caption{(a) Projection of the conventional tetragonal unit-cell of FeNi containing two Fe and two Ni atoms on to the (100) plane, with the fractional coordinates along the [001] $c$-axis   marked. The origin of the cell-axes is set to an iron atom. The small circles represent the octahedral interstices. (b) Configurations of octahedral interstices with fractional $z$ coordinates indicated; metal atoms omitted  for clarity. If each configuration is occupied completely by hydrogen atoms, the net contribution to the unit cell would be one hydrogen atom because the sites may be shared with adjacent cells.} \label{fig01}
\end{figure}


The total energies were obtained using the Perdew-Burke-Ernzerhof  realisation of the generalised gradient approximation for the exchange-correlation potential \cite{Perdew:1996}. A plane-wave basis set with a cut-off energy of 720 eV was used in all calculations. Integrations over the three-dimensional Brillouin zone were performed using the tetrahedron method over an $11\times11\times11$ Monkhorst-Pack mesh for the conventional unit cell, and the total energy is converged to an accuracy of $10^{-8}$ eV. All structures were assumed to be in the ferromagnetic state. The initial lattice parameters were adapted from experimental values ($a=3.568$\,\AA, $c=3.5901$\,\AA) \cite{Abdu:2002}. Internal atomic relaxation and geometry optimisation were achieved when the change in total energy between two ionic steps was smaller than $10^{-3}$\,eV\,atom$^{-1}$. Spin polarisation was considered for all calculations, and the convergence of these computational parameters was checked carefully. All thermal contributions were neglected in the present calculations. The formation energy of each structure was calculated using the following formula:

\begin{equation}
\Delta E = E\{{\rm FeNi}+n{\rm H}\} - E\{{\rm FeNi}\} - n\times\frac{1}{2} E\{{\rm H}_2\},
\end{equation}
where $E\{{\rm FeNi}+n{\rm H}\}$ and $E\{{\rm FeNi}\}$ are the total energies of FeNi + $n$H and FeNi structure, respectively. $E\{{\rm H}_2\}$ is the total energy of hydrogen in its gaseous state.

Table~\ref{tab01}  summarises the changes in lattice parameters and formation energies from the first-principles calculations. For FeNi, the calculated lattice parameters (representing 0\,K) are about 0.2\% smaller than the experimentally measured ambient-temperature values in both the $a$ and $c$ directions. The $c/a$ ratio closely matches the experimental values.

The location of the hydrogen atoms for the purposes of calculation are designated A, B, C and D in Figure~\ref{fig01}b, terminology consistently used in Table~\ref{tab01} an throughout to express the variety of configurations of octahedral interstices studied, together with their outcomes.

The insertion of a single hydrogen atom at a mole fraction $x=\frac{1}{5}$, in octahedral interstices in a manner consistent with either the configuration A or B yields the most favoured location in terms of $\Delta E$. The $c/a$ ratio is reduced because the positive strain $\epsilon_{\rm H-Fe}$ along Fe-H-Fe which is parallel to the $a$ parameter  is much greater than along Ni-H-Ni which is parallel to the $c$-axis.

The $\Delta E$ is greater  when both A and B  are occupied with $x=\frac{1}{3}$.  The $c/a$ ratio is only slightly smaller than for pure FeNi  because $\epsilon_{\rm H-Ni}$ along the Ni-H-Ni orientation becomes similar to $\epsilon_{\rm H-Fe}$ (Table~\ref{tab01}, Figure~\ref{fig02}).  It follows that placing hydrogen atoms at the equivalent A or B sites, in any combination, does not increase the tetragonality of $\upgamma''$. This is consistent with the Fe-H-Fe chains parallel to the $a$-axes on the (100) planes being  more compliant to the insertion of hydrogen at $z=0,1$.   Adding a third H-atom (A+B+C) uniformly expands the lattice but with only a small increase in $c/a$.

\begin{table}[h!]
\centering
\begin{tabular}{lcccccccc}
\toprule
& & & & & \multicolumn{2}{c}{Nearest neighbours} \\
& $\Delta E$ / eV & $a$ / \AA & $c$ / \AA  & $c/a$ & H-Ni\,/\,\AA  & H-Fe\,/\,\AA  & $\epsilon_{\rm H-Ni}$ & $\epsilon_{\rm H-Fe}$  \\
\midrule
2FeNi (experimental) & - & 3.568\phantom{0} & 3.5901 & {\bf 1.0062} & &   \\
2FeNi (calculated) & - & 3.5587 & 3.581\phantom{0} & 1.0063 & {\em 1.779} & {\em 1.790}  \\[5pt]
2FeNi + H (tetrahedral) & 0.347 & 3.6133 & 3.7288 & 1.0320 & 1.613  & 1.687  & &   \\[5pt]
2FeNi + H\{A\} & 0.139 & 3.6248 & 3.6124 & 0.9966 & $1.806\parallel c$  & $1.812\parallel a$ & 0.009 & 0.019\\
2FeNi + 2H\{A+B\} & 0.271 & 3.6598 & 3.681 & 1.0057& $1.841\parallel c$  & $1.830\parallel a$  & 0.029 & 0.029 \\
2FeNi + 3H\{A+B+C\} & 0.302 & 3.6996 & 3.7232 & 1.0064 & $1.850\parallel a$ & $1.862\parallel c$  & 0.040 & 0.040\\[5pt]
2FeNi + H \{C\} & 0.178 & 3.5584 & 3.7604 & 1.0568 & $1.779\parallel a$  & $1.880\parallel c$ & 0.000 & 0.050 \\
2FeNi + 2H \{C+D\} & 0.364 & 3.5642 & 3.9302 & 1.1027 & $1.782\parallel a$  & $1.965\parallel c$ & 0.002 & 0.098\\[5pt]
2FeNi + 3H \{A+C+D\} & 0.386 & 3.6478 & 3.8705 & 1.0611& $1.824\parallel a$ & $1.935\parallel c$  & 0.025 & 0.081\\
2FeNi + 4H\{A+B+C+D\} & 0.327 & 3.7236 & 3.8117 & 1.0237& $1.862\parallel a$ & $1.906\parallel c$ & 0.047 & 0.065 \\
\bottomrule
\end{tabular}
\caption{Calculated lattice parameters and formation energies at 0\,K for configurations with one hydrogen atom added at octahedral ($z=0 $ or 1/2) or tetrahedral ($x=1/4$, $y=1/4$, $z=1/4$) interstitial positions.  The designations A,B,C,D refer to Figure~\ref{fig01}b. The unit cell contains 2Fe and 2Ni atoms so the designations indicate the correct ratios of Fe, Ni and H atoms in a cell. The number of octahedral interstices per cell is the same as the number of metal atoms, i.e., 4. There are 8 tetrahedral interstices per metal atom. The strains along Fe-H-Fe and Ni-H-Ni are represented by $\epsilon_{\rm H-Fe}$ and $\epsilon_{\rm H-Ni}$, respectively, stated to three significant figures. The italicised numbers in row two represent the distance between an octahedral interstice and a metal atom, since the 2FeNi does not contain hydrogen.}
\label{tab01}
\end{table}

\begin{figure}[!ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{calculations.eps}
\caption{Changes in formation energy and  $c/a$ ratio from first-principles calculations as a function of the location and concentration of  hydrogen, with the latter quantities defined by the designations A,B,C,D in Figure~\ref{fig01}b in the case of the octahedral interstices. } \label{fig02}
\end{center}
\end{figure}

Placing a single hydrogen atom in an octahedral interstice at `C' (or its equivalent `D')  leads to a dramatic increase in tetragonality, with $\Delta E=0.178$\,eV being slightly larger than the most favourable location `A' ($\Delta E=0.139\,$eV), Table~\ref{tab01}. At those locations, the Fe-H-Fe chains are parallel to the $c$-axis with $\epsilon_{\rm H-Fe}>\epsilon_{\rm H-Ni}$, leading to the increased $c/a$ ratio. Two hydrogen atoms at C+D maximises the  tetragonality. The general conclusion from all the data presented in Figure~\ref{fig02} is that locating hydrogen atoms in configurations C,D where it lies in the nickel-atom planes, maximises tetragonality, with hydrogen concentrations in excess of $x=\frac{1}{3}$ causing reductions in $c/a$ because atoms would necessarily have to locate on the A,B sites.

It is useful to understand why the Fe-H-Fe chains are more flexible relative to Ni-H-Ni.  Figure~\ref{fig03} shows the calculated total and partial local density of states (DOS) for the hydrogen atom, and the nearest Fe and Ni atoms. At both configurations, the H-$s$ orbital commonly overlaps with the $s$ and $p$ states of the nearest Fe and Ni atoms in the energy range $-12$\,eV to $-7$\,eV. In the Fe-H bonding, the contribution of the $p$-state is relatively large, while in the Ni-H bonding, the $s$-state at lower energy levels plays a more significant role. Moreover, the peak height of the DOS for Ni-H bonding is greater than for the Fe-H bonding, and the energy value is lower, indicating a stronger bond. This is confirmed by the  maps in Figure~\ref{fig04}  which show the charge density contour plots in the plane normal to the $b$-axis intersecting the Fe, Ni, and H atoms for configurations `A' and `C'. In both cases, the charge density of H overlaps more with Ni than with Fe, and consequently, the chain of atoms in the Fe-H-Fe configuration is extended to a significantly larger spacing compared to the Ni-H-Ni chain on the (002) plane.

 \begin{figure}[ht!]
\centering
\subfigure[]{
\includegraphics[height=0.6\textheight]{DOSa.eps}}
\subfigure[]{
\includegraphics[height=0.6\textheight]{DOSb.eps}}
\caption{The calculated total and partial local density of states  for Fe, Ni and H on configurations (a) `A' and (b) `C'. Fe and Ni refer to the atoms in the nearest vicinity to H atoms. The spin down DOS values are factored by $-1$, and the Fermi levels ($E_{\rm F}$) are set to zero. For comparison, the density corresponding to the $s$ and $p$-states of Fe and Ni atoms was factored by 10.}
\label{fig03}
\end{figure}

\begin{figure}[ht!]
\centering
\subfigure[]{
\includegraphics[height=0.35\textheight]{chargea3.eps}}
\subfigure[]{
\includegraphics[height=0.35\textheight]{chargeb3.eps}}
\caption{Charge density contour plots in the plane normal to the $b$-axis, cutting the Fe, Ni and H atoms for configurations (a) `A' and (b) `C'. The horizontal axis represents [100] direction, while the vertical axis is parallel to [001].}
\label{fig04}
\end{figure}


It is noted that insertion of a single H-atom into a tetrahedral interstice leads to the largest $\Delta E$ for $x=\frac{1}{5}$ and further calculations for larger concentration, not listed here,  do not improve the situation. Tetrahedral interstices are therefore not considered further.


The presence of hydrogen in $\upgamma''$ does reduce the net magnetic moment per atom which in turn should reduce the saturation magnetisation of the material (Table~\ref{tab02}). However, the reduction of total $\mu_{\rm B}$  for configurations `C' and `C+D' which have large $c/a$ ratios might be tolerated if the magnetocrystalline anisotropy is enhanced to produce more-permanent magnets.


Several studies in the past have investigated the magnetocrystalline anisotropy energy in tetrataenite using both first-principles calculations and experimental methods \cite{izardar2020interplay,lewis2014magnete,edstrom2014electronic}. Defining the energy difference $E_{[100]}-E_{[001]}$, where $E_{[100]}$ and $E_{[001]}$ are the total energies, including spin-orbit coupling, we obtained with magnetisation aligned along the [100] and [001] directions, respectively. This calculation is performed for both the structure without hydrogen and in the of hydrogen in the `A' and `C' configurations, using optimised lattice parameters. 

Since the magnetic easy axis is [001], the anisotropy should be positive \cite{izardar2020interplay,neel1964magnetic}. For FeNi without hydrogen, it is calculated to be 0.64 $\upmu$eV/f.u., consistent  76-128 $\upmu$eV/f.u. reported in previous studies\cite{neel1964magnetic,kojima2014fe}. The corresponding estimates for configurations `A' and `C' are 221.0 and 133.9 $\upmu$eV/f.u., respectively, both significantly greater than without hydrogen, suggesting that the presence of hydrogen enhance the magnetocrystalline anisotropy in both configurations.

\begin{table}[hpb]
\centering
\begin{tabular}{lcccccccc}
\toprule

  & Total & Fe & Ni & H \\
&  $\mu_{\textrm{B}}$& $\mu_{\textrm{B}}$ & $\mu_{\textrm{B}}$ & $\mu_{\textrm{B}}$ \\
\midrule
2FeNi & 6.583 & 2.662 & 0.629 & - \\[5pt]
2FeNi + H (tetrahedral)   & 6.426 & 2.629 & 0.593 & $-0.018$\\[5pt]
2FeNi + H \{A\}  & 6.404 & 2.547 & 0.665 & $-0.019$\\
2FeNi + 2H \{A+B\}   & 6.008 & 2.342 & 0.680 & $-0.017$\\
2FeNi + 3H \{A+B+C\}  & 5.385 & 2.221 & 0.491 & $-0.013$ \\[5pt]
2FeNi + H \{C\} & 6.364 & 2.659 & 0.528 & $-0.010$ \\
2FeNi + 2H \{C+D\}  & 6.071 & 2.682 & 0.359 & $-0.006$\\[5pt]
2FeNi + 3H \{A+C+D\}  & 5.478 & 2.469 & 0.289 & $-0.012$\\
2FeNi + 4H \{A+B+C+D\} & 4.906 & 2.27 & 0.209 & $-0.013$ \\
\bottomrule
\end{tabular}
\caption{Calculated total magnetic moment and magnetic moments per atom ($\mu_{\textrm{B}}$).}
\label{tab02}
\end{table}

\section{Conclusions}

It is possible that the alloying of tetrataenite with hydrogen, for example by electrochemical charging, can lead to significant increases in the tetragonality of its unit cell, and hence in its magnetocrystalline anisotropy and coercivity. The optimum location of the hydrogen atoms in this context is in the octahedral interstices at $\frac{1}{2}z$ where it is surrounded on the plane (002) by four near neighbour nickel atoms. This is because the chain of atoms Fe-H-Fe normal to that plane (parallel to $c$) is extended to a significantly larger spacing than any Ni-H-Ni chain on the plane (002) because the Ni-H bonding is significantly stronger than Fe-H bonds. The $c/a$ ratio and {magnetocrystalline anisotropy therefore increases though the saturation magnetisation is reduced by a small amount.

Although the energy $\Delta E$ for the introduction of a single hydrogen atom on (002) is about 28\%\ greater than on (001), it still is reasonably small (0.178) so it may be feasible for hydrogen atoms to reside there during electrochemical-charging.

Insertion of hydrogen into the octahedral sites on (001) where it is surrounded by iron atoms as nearest neighbours is not useful in altering $c/a$ since the expansion of the Fe-H-Fe chains on that plane lead to increases in the $a$-parameter. Nevertheless, the MAE calculation indicate an enhancement in magnetocrystalline anisotropy.

\section*{Declaration of Competing Interest }
The author declare that they have no known competing financial interests or personal relationship that could have appeared influence the work reported in this paper.

\section*{Acknowledgement}
This work was funded by the Fundamental R\&D Program of the Korea Institute of Materials Science (KIMS), Grant No.~PNK9630. This work was also supported by the Industrial Strategic Technology Development Program (No.~20010942) funded by the Ministry of Trade, Industry \& Energy (MOTIE, Korea).


\bibliography{reference}

\end{document} 