## 2 Calculation Methods  
### Elastic Constants  
In this work, by using the large-scale atomic/molecular massively parallel simulator (LAMMPS) [35], we performed the molecular statics (MS) and molecular dynamics (MD) simulations for structural relaxations and the calculations of the elastic constant tensor \(C\). The Fe-Fe interaction is described by the embedded atom method (EAM) potential developed by [PERSON] [36], and the Fe-H interaction is fitted by [PERSON] et al. [37]. \(C\) is defined by [PERSON]'s law:  
\[\sigma=C\epsilon \tag{1}\]  
where \(\sigma\) and \(\epsilon\) are the stress and strain tensors, respectively. Each component of \(\sigma\) is calculated according to the following formula [38]:  
\[\sigma_{ij}=-\frac{1}{V}\left(\sum_{a}^{N}m_{a}v_{a_{i}}v_{a_{j}}+\sum_{a}^{N ^{\prime}}r_{a_{i}}f_{a_{j}}\right) \tag{2}\]  
\(V\) is the volume;\(f_{a}\) is the force applied to atom a by other atoms. Respectively, \(m_{a}\), \(v_{a}\), and \(r_{a}\) are the mass, velocity, and position vector of atom a. Different orientations are denoted by \(i,j\); \(N\) is the total number of atoms in the system; and \(N^{\prime}\) is the number considering the periodic image atoms. For cubic structures, \(C\) has only three independent components: \(C_{11}\), \(C_{12}\), and \(C_{44}\). If a deformation is applied to a single direction, it satisfies:  
\[C_{ij}=\sigma_{i}/\varepsilon_{j} \tag{3}\]where \(i,j=1,2,4\), which, respectively, represent the stretching along the \(x,y\) axis, and the angle change between the \(y\), \(z\) directions. In this work, \(\varepsilon_{j}=0.00001\); for accuracy, we rotate the \(x,y\), and \(z\) axes, and then average the results when calculating \(C_{ij}\). The influence of interstitial hydrogen atoms on the elastic parameters is often decomposed into a linear combination of various factors. For example, [PERSON] et al. considered the elastic parameters of \(\alpha\)-Fe single crystals under high \(c_{H}\) as effects of the chemical interactions and volume expansion [27]. [PERSON] et al. decomposed the increase of the Young's modulus under hydrogen-charging conditions into the hydrogen-induced changes in cohesion energy and the effect of stresses [39]. Similarly to them, we write the elastic constant \(C_{ij}\) of the structure \(q\) as follows:  
\[ \begin{array}{c}C_{ij}^{q}(c_{H})=\langle C_{ij}^{p}(c_{H})\rangle+\Delta C_ {ij}^{dis,q}(c_{H})+\Delta C_{ij}^{def,q}\big{(}c_{H},c_{def}\big{)}\\ =C_{ij}(0)+\Delta C_{ij}^{coh}(c_{H})+\Delta C_{ij}^{vol}(\Delta V(c_{H}))\\ \qquad\qquad+\Delta C_{ij}^{dis,q}(c_{H})+\Delta C_{ij}^{def,q}\big{(}c_{H},c_ {def}\big{)}\end{array}  \tag{4}\]  
in which \(p\) is the structure number. \(\langle C_{ij}^{p}(c_{H})\rangle\) is the average value of \(C_{ij}\) for all structures with hydrogen concentration of \(c_{H}\), representing the elastic constant of the solid solution of hydrogen perfectly uniformly distributed in \(\alpha\)-Fe. It is decomposed into the elastic constant \(C_{ij}(0)\) of pure \(\alpha\)-Fe single crystal; the change \(\Delta C_{ij}^{coh}(c_{H})\) caused by the cohesion of iron and hydrogen; and the change \(\Delta C_{ij}^{col}(\Delta V(c_{H}))\) caused by the volume change of \(\Delta V(c_{H})\). \(\Delta C_{ij}^{dis,q}(c_{H})\) is the deviation of the elastic constant of structure q from the average result, which is caused by the deviation of hydrogen from the ideal uniform distribution in microscopic distribution. \(\Delta C_{ij}^{def,q}\big{(}c_{H},c_{def}\big{)}\) is the change in elastic constant caused by the concentration of other defects (\(c_{def}\)) except for the interstitial hydrogen atoms. In this work, as \(\Delta C_{ij}^{def,q}\big{(}c_{H},c_{def}\big{)}\) is not considered for the time being, Equation (4) can be simplified as:  
\[ \begin{array}{c}C_{ij}^{q}(c_{H})=\langle C_{ij}^{p}(c_{H})\rangle\ +\Delta C_{ij}^{dis,q}(c_{H})\\ =\ C_{ij}(0)+\Delta C_{ij}^{coh}(c_{H})+\Delta C_{ij}^{vol}(\Delta V(c_{H}))+ \Delta C_{ij}^{dis,q}(c_{H})\end{array}