Isotopic Separation of Helium through Nanoporous Graphene Membranes: A Ring Polymer Molecular Dynamics Study

Microscopic-level understanding of the separation mechanism for two-dimensional (2D) membranes is an active area of research due to potential implications of this class of membranes for various technological processes. Helium (He) purification from the natural resources is of particular interest due to the shortfall in its production. In this work, we applied the ring polymer molecular dynamics (RPMD) method to graphdiyne (Gr2) and graphtriyne (Gr3) 2D membranes having variable pore sizes for the separation of He isotopes. We found that the transmission rate through Gr3 is many orders of magnitude greater than Gr2. The selectivity of either isotope at low temperatures is a consequence of a delicate balance between the zero-point energy effect and tunneling of $^4$He and $^3$He. RPMD provides an efficient approach for studying the separation of He isotopes, taking into account quantum effects of light nuclei motions at low temperatures, which classical methods fail to capture.


Introduction
The use of two-dimensional (2D) materials is an open field for many technological applications as well as for basic science. 1,2 Graphene, one of the pioneer materials that started this very active area of research, acts as an almost impermeable sheet for most atomic and molecular species due to its high electron density surrounding the aromatic rings. 3,4 Some years ago, experiments by Lozada-Hidalgo et al. 5,6 and recently by Creager et al. 7 have shown that small charged particles such as the proton and its isotopes can penetrate through a monolayer of pristine graphene. There is no complete explanation despite many theoretical calculations, [8][9][10] although some kind of chemical interaction seems to be responsible for this process. Even more recently, it has also be found that hydrogen can permeate pristine graphene. 11 However, the penetration energy barrier for other neutral atoms is significantly higher and therefore supports the widely recognized notion of graphene impermeability. 12 Introducing defects and moderate annealing into the membrane, and nanopores of different sizes, may enable graphene and other 2D layers to act as atomic and molecular sieves. [12][13][14][15][16][17][18] The use of 2D materials as a filter at the molecular level is one of the most interesting applications. 14 Besides using 2D membranes with fabricated nanopores, some compounds contain nanopores in their structures at different positions and sizes; therefore, they are more naturally suited for filtering at the molecular level.
An interesting class of 2D nanoporous material is graphyne, 19,20 which is composed of sp − sp 2 hybridized C atoms and can be considered as a graphene derivative. In graphyne, one-third of C−C bonds have been replaced by mono-and poly-acetylenic units ( -C ---C -) N .
Baughman et al., 21 theoretically proposed the first stable structures of graphyne, which was two decades later synthesized on a copper substrate. 22,23 N defines the number of acetylenic linkages and, consequently, the size of the uniformly distributed and repeating sub-nanometer triangular pores. They are termed as graph-N -yne membranes, such as graphdiyne (Gr2), graphtriyne (Gr3), etc., for N = 2 and 3, respectively. 24 The chemical and mechanical properties of these graphyne membranes have many useful features, such as they are chemically inert and stable at ambient temperatures, [25][26][27] and flexible enough to withstand deformations induced by high pressures. 28,29 The above-mentioned properties of graphynes, coupled with their unique geometrical structure, make them an excellent candidate to be utilized in gas separation and water filtration technologies. Indeed, one can find many reports on the purification of gases such as H 2 , 30 N 2 , 31 and O 2 32 from a mixture of gases and desalination and filtration of water 24,29,33 in the literature.
There is a growing worldwide demand for helium purification from its natural resources due to the shortfall in its production. 34 It has numerous industrial and scientific applications such as superconducting magnets, space rockets, arc welding, etc. In particular, the lighter isotope, 3 He, also plays a pivotal role in fundamental research, for example, in neutronscattering centers, ultracold physics and chemistry, etc. 35 The relative abundance of 3 He is low (≈1.34×10 −6 %) 36 in comparison to its heavier isotope, and its extraction from natural gas is usually done by means of expensive cryogenic distillation and pressure-swing adsorption methods. 37 An alternate and more energy-efficient process is to use the 2D porous membranes in isotopic gas separation since they usually do not involve costly liquefaction of the gases. 38  and references mentioned therein. In the permeation of H + and D + on a pristine graphene, the RPMD method also has been used, 8 and the authors point out, again, the importance of considering both ZPE and tunneling effects for isotopic separation.
Coupled with the major scientific and industrial appeal, in this paper, we examine the feasibility of separation of He isotopes using 2D Gr2 and Gr3 membranes at low temperatures (20−250 K) using the RPMD method. A significant part of the success of the RPMD method is attributed to the fact that it gives the exact quantum-mechanical rate coefficient for the transmission through a parabolic barrier, 51 which is advantageous for the present investigation. The primary objective of this study is to provide a reliable estimate of the selectivity and to show that RPMD is a necessary alternative to study these processes where quantum effects are expected to be very important. To the best of our knowledge, the RPMD method has not been applied for the He isotope separation previously and therefore presents an excellent opportunity to test the accuracy of this method.
The paper is organized as follows: in the next section (section 2), we provide the details of the RPMD approach along with the PES used in the present study. The results of RPMD rate coefficients and selectivity have been compared with earlier studies [47][48][49] in section 3.
Concluding remarks are provided in the last section (section 4). The potential energy curves of He along the transmission pathway (z− coordinates) to both membranes are illustrated in Figure 1, while the contour plots in Figure 2(a) and (b) highlight the in-pore displacements of He along x− and y− coordinates. It is quite evident from both these plots that the minimum energy path (MEP) for an effective He transmission may correspond to a straight line perpendicular to the center of the pore. For the transmission through the Gr2 membrane, the He atom lying at the center of a pore is the saddle point. The maximum potential barrier height for the MEP is 37.85 meV, which is similar to the previous studies. 35,47,48 In contrast, MEP for the He − Gr3 PEC is devoid of any potential barrier; instead, a large well (of depth 17.74 meV) appears within the vicinity of the membrane. 49 These results suggest that He penetration of larger Gr3 pores should be much easier than that of smaller Gr2 pores, analogous to the report on water permeation by Bartolomei et al. 24 Furthermore, along this MEP and He atom-membrane perpendicular distance of around 6.0 Å, the interaction potential reaches a plateau and therefore can be characterized as an asymptotic reactant site. It is interesting to note that the interaction potential steeply rises to a high value for any movement along the in-pore

Ring polymer molecular dynamics (RPMD) method
The classical Hamiltonian for the system composed of 4 He or 3 He atom under the influence of external potential arising from M carbon atoms of fixed Gr2 or Gr3 membrane can be written in the atomic unit as: where,p andr He (orr C(i) ) are the momentum and position vectors of the ith atom of mass m i . In the RPMD method, the He atom is treated as n classical replicas of the original particle, each connected with its nearest neighbor by harmonic springs. The modified ring polymer Hamiltonian has the form: where, n is the number of classical beads representing quantum He atom which are connected by a harmonic potential with force constant ω n (=β nh ). β n ≡ β/n is the reciprocal temperature of the system, β = 1/k B T . k B is the Boltzmann constant, T is the system temperature, andh is the Dirac constant (h=1 in atomic units). p (j) and r (j) are the momentum and position vectors of the jth bead in the ring polymer necklace of He atom, respectively. The ring polymer trajectory now evolves in f = 3n total degrees of freedom (in atomic Cartesian coordinates and including translation and rotational degrees of the freedom of the entire system which is a convenient method to propagate RPMD trajectories 69 ).
We perform the RPMD simulations at temperatures T below 250 K because, previously, it was observed that the maximum selectivity was obtained at low temperatures. 47, 48 We choose seven different T = 20, 30, 50, 100, 150, 200, and 250 K, to study the variation of selectivities with T and to find an optimal temperature range for isotopic separation.
Furthermore, the number of beads is a measure of the resolution of the path integral calculations, i.e., RPMD calculations scale linearly with n. 54 However, since RPMD calculations are approximately n times slower than the purely classical calculations, 54 after several trials, we have meticulously chosen the number of beads to be 128. This value of n gives a fair compromise between computational cost and accuracy (quantum effects of ZPE and tunneling) at low temperatures. The values of n=1, will correspond to a classical calculation.
Since RPMD is simply classical molecular dynamics in an extended (n-bead imaginary time path integral) phase space, the ring polymer rate coefficient can be expressed as: 51,52,69 where, Q (n) r (T ) is the n-bead path integral approximation to the quantum mechanical partition function of the reactants per unit volume, andc Here, subscript 0 and t indicates time, δ(r 0 ) is a delta function centered at r 0 , v(r 0 , p 0 ) is the velocity, and h(r t ) is the Heaviside step function. The RPMD rate coefficient in Equation 4 is not straightforward to solve numerically. Therefore, we introduce the Bennett-Chandler factorization scheme 74,75 to simplify Equation 4 that can be solved numerically without compromising its generality. 75 The method has been extensively discussed previously 52,68,69,76 and will not be repeated here. Briefly, in this approach, a reaction coordinate s(r) is defined, which monitors the progress of a reaction from the reactant (s > 0) to the product (s < 0) site. For the He transmission through 2D membranes, a reasonable reaction coordinate is s(r) =z, wherez is the z-component of the He atom centroid, which follows the minimum energy path of Figure 1. Within the Bennett-Chandler approach, the RPMD rate coefficient for a process in which reactants and products separated by a dividing surface at s ‡ (for instance, center of the pore) can be expressed as a product of two terms: Here, the first factor, k QTST (T ), is the centroid-density quantum transition-state theory W (s)'s can be computed by employing the umbrella integration procedure of Kästner and Thiel. [77][78][79] To calculate the PMF profiles, the reaction coordinate s of the He atom has been divided into 130 equally spaced windows (of width 0.05 Å) within the range from 6.00 Å (reactant site, s ∞ ) to −1.00 Å (product site). The PMFs in Equation 7 is then calculated as: with Here, N windows is the number of biasing windows placed along the reaction coordinate.
The strength of the force constant (k i ) of the harmonic biasing potential was chosen to be 2.72×10 −3 T (K) eV a −2 0 . N i is the total number of steps sampled for window i,s i and σ 2 i are the mean value and variance for the trajectory calculated for the ith window. In each umbrella sampling windows, 100 trajectories with different initial configurations were propagated for 100 ps following an initial equilibration period of 20 ps in the presence of an Andersen thermostat. 80 The ring polymer equation of motion were integrated using velocity Verlet integrator with a step size of 0.1 fs that involves alternating momentum updates and free ring polymer evolutions. 58 The second term in Equation 6, κ(t p ), is the long-time limit of a time-dependent ring polymer transmission coefficient or the ring polymer recrossing factor 63,69,70,76 and is a dynamic correction to k QTST . Typically this factor is calculated at the top of the free energy barrier on the PMF profile so as to minimize the time required to reach plateau value. 75 This factor ensures that the final k RPMD value is independent of the choice of the dividing surface. 52 The mathematical expression for κ(t p ) can be written as: 64,70 where, δ[s(r 0 )] constrains the initial configurations to the dividing surface,ṡ(r 0 ) is the velocity factor that accumulates the flux through the dividing surface, and h[s(r t )]is a Heaviside function that gathers trajectories that have crossed over to the product side of the dividing surface. h[ṡ(r 0 )] is basically a normalization factor that ensures κ(t → 0 + ) = 1. To calculate the recrossing factor, a long "parent" trajectory for He atom of length 2 ns has been carried out after an initial thermalization period of 20 ps in the presence of Andersen thermostat with its centroid pinned at the dividing surface using RATTLE algorithms. 81 After each 2 ps propagation period of the parent trajectory, 100 trajectories have been generated that have initial position of the parent trajectory, but their momenta is randomly generated from a Boltzmann distribution. These "child" trajectories are then propagated for 1 ps in the absence of thermostat and dividing surface constraints. Since the RPMD method averages points in configurational space on both sides of the potential energy barrier, modeling the effects of the tunneling and ZPE and leaving a mark on the free energy profiles, 57 we can roughly attempt to identify these properties. Comparing the PMF profiles of 4 He at a particular temperature with the corresponding 3 He's on Gr2 indicate the existence of a prominent ZPE effect; heavier isotope has a lower free energy at the TS than the lighter ones except at 20 K (see Figure 3(B)), in which quantum tunneling may play a crucial role. However, the difference between them at TS decreases with increasing temperature. For example, the difference at 30 K is 2.3 meV (see Figure 3(B)), which decreases to 0.1 meV at 200 K (see Figure 3(C)). Similarly, the 3 He and 4 He free energy comparison for Gr3 follows the same trend observed for Gr2, i.e., they decreases with increasing temperature. Within 20 − 30 K, 4 He has a smaller TS free energy than 3 He (≈1.1 meV, Figure 4(B)), and with the increase in temperature in the range 50 − 150 K, the heavier isotope progressively requires more free energy to reach TS, i.e., PMF profiles for both isotope are almost identical. Within 200 − 250 K, the free energies of 4 He are consistently larger than those obtained for 3 He (see Figure 4(C)). These observations can be notionally interpreted as follows: within 20 − 30 K, the free energy barrier is too broadened for effective tunneling, and consequently, the ZPE effect takes precedence. 55 There is a competition between them at moderate temperatures (50 − 150 K). At still more elevated temperature (200 − 250 K), the greater tunneling probability of 3 He may facilitate a smaller TS free energy compared to 4 He. Note that this analysis is conceptual because these two quantum effects cannot be rigorously separated within the RPMD formalism.
There are considerable differences observed between the PMF profiles obtained by the RPMD method with the classical ones at low temperatures (see supplementary information).
This, therefore, reinforces our earlier argument of the existence of quantum effects in these transmission processes, which the classical calculation fails to capture. For the He transmission through Gr2 membrane, the free energy barrier height in the classical calculations is always lower than those obtained by the RPMD method. The maximum difference for 3 He obtained at 100 K (≈ 10 meV) and that for 4 He obtained at 20 K (≈ 12 meV). Note that the free energy values exponentially contribute to the rate coefficients. As expected, these PMF plots tend to merge with increasing temperature as the quantum effects fade.
Comparison between the classical and RPMD method on Gr3 shows that, for both isotopes, the maximum free energy difference at the TS is found at the lowest temperature (≈ 3 meV), and this difference smoothly decreases with increasing temperature.  this study on both membranes are always close to unity (0.98 − 1.00) and are moved to the supplementary information ( Figure S5 and Figure S6). The corresponding plateau value of the transmission coefficients, κ(t → ∞), are provided in Table S1. The κ(t) value is always close to unity (0.98 − 1.00). Clearly, recrossing dynamics do not play any significant role for the He transmission on both membranes. Particularly, from 50 K and higher temperature, where κ(t) value is always 1.00. This implies that most of the He trajectories that reach the center of the pore of either Gr2 or Gr3 would overcome the free energy barrier and transport to the other side of the dividing surface. At 20 K, κ(t) for 3 He is marginally smaller than that obtained for 4 He (0.98 and 0.99 respectively). Note that the classical recrossing factor is always 1, even at low temperatures.

Thermal rate coefficient
Since the plateau of κ(t) is always close to unity, there is practically no significant difference between the values of k QTST and k RPMD on both membranes. The variation of k RPMD with temperature is plotted in Figure 5 and Figure S7 of the supplementary material, while the numerical values of k RPMD and k QTST are reported in Table S1 (see supplementary information). It is evident that k RPMD 's on Gr3 is many orders of magnitude greater than that obtained on Gr2. The maximum difference was observed at the lowest temperature (by a factor of 10 17 at 20 K), and with the increase in temperature, this difference decreases rapidly (by a factor of 10 2 at 250 K) as the k RPMD vs. T curves start to converge. This is due to the fact that the value of k RPMD on Gr3 does not change drastically with temperature (confined within 3.63×10 10 s −1 − 1.01×10 11 s −1 ). It starts to decrease slightly with temperature up to 50 K and then increases successively with increasing T . On Gr2, however, for both isotopes, k RPMD increases manifold with temperature, in agreement with the previous calculations. 35,47,48 This increment in the rate coefficient is more pronounced at the modest rise in T in the low temperature regime than at the high temperature range. For example, the increase in the value of k RPMD is in the order of 10 5 for the temperature rise from 20 K to 30 K or 30 K to 50 K. However, k RPMD increases much less than an order for the consecutive 50 K temperature jumps, starting at 150 K. The above observations can be explained by inspecting the corresponding PMF profiles.
At low temperatures (20 K − 50 K), there is practically no free energy barrier for the transmission of either isotope on Gr3. However, on Gr2, the free energy barrier is already above 60 meV at 20 K. The slowly moving He atoms do not have enough kinetic energy to overcome this barrier to reach the other side of the membrane. Therefore, it is not surprising that k RPMD 's on Gr3 is many orders of magnitude greater than those obtained on Gr2, and the actual He flux through the pores of Gr2 membranes will be extremely slow. 47,48 Similar arguments can be presented at higher temperatures, although the increased kinetic energy of the incoming He atom will contribute to smaller differences in the k RPMD values. We also point out the analogous quantum wave packet observation reported on the Gr2 and holey graphene (P7) sheet having a more diffused pore. 48 Finally, when comparing the classical and RPMD rate coefficients, it is obvious that the classical method overestimates the rate coefficient within the whole temperature regime studied in this work. This discrepancy is particularly more apparent at low temperatures, when the quantum effects dominate, and on the Gr2 membrane. For example, at 20 K and for Gr2 membrane, the classical rate coefficient is more than three orders of magnitude greater than the corresponding k RPMD .
However, they do closely follow a similar temperature dependence as obtained by the RPMD calculations.

4 He/ 3 He selectivity
The 4 He/ 3 He selectivity, defined as the ratio between the RPMD rate coefficient of the heavier He isotope to the lighter one, k RPMD ( 4 He)/k RPMD ( 3 He), is plotted as a function of temperature in Figure 6, and the corresponding values are reported in Table S1 (see supplementary   information). For Gr2 membrane, 4 He/ 3 He selectivity increases for the temperature rise from 20 K to 30 K. At 30 K, the maximum selectivity is obtained (≈ 2) favouring the heavier isotope and simultaneously indicating pronounced quantum effects as discussed previously.
However, this increased selectivity should be accompanied by a decreased permeability. 82 With the increase in temperature, this selectivity ratio progressively becomes smaller and starts to flatten out, commencing from 100 K. Although at higher temperatures (> 100 K), the transmission of the lighter isotope is favoured to some extent. We note the analogous selectivity profile was obtained on the P7 sheet. 48 Similarly, for the Gr3 membrane, the maximum 4 He/ 3 He selectivity obtained at 20 K (1.17), which gradually decreases till 50 K, and then remains almost constant (0.87 to 0.83) marginally preferring 3 He for higher tem-perature regime. The variation of the selectivity is fully consistent with the nature of their PMF profiles that inherits quantum effects. It is interesting to note that the selectivity plot for both Gr2 and Gr3 approach each other with increasing T , one from the top and another from the bottom, and at 250 K, they almost become equal (≈ 0.8 and partially in favour of 3 He). This is reasonable since as the quantum effects diminish with temperature and with fewer constraints, the lighter isotope transmission will be slightly favoured. The selectivity ratios obtained by the RPMD method show good agreement with those obtained by the quantum three-dimensional wave packet propagation (WP3D) results obtained by Gijón et al. 48 and Hernández et al. 49 for the whole temperature range studied in this work. In general, the RPMD selectivity profile closely resembles the WP3D results; however, they do overestimate and underestimate to some degree. For example, on Gr2, the RPMD selectivity ratio overestimates the quantum calculations for the whole temperature range. The maximum difference between these two methods is found at 30 K (0.2). Similarly, on Gr3, the RPMD method seems to overestimate the selectivity ratio within 20 K − 30 K and underestimate them for higher temperatures by a maximum of 0.15. Overall, the difference between RPMD and WP3D results for the selectivity ratio lies within 9−23% on Gr2 and 9−15% on Gr3 in comparison to the quantum calculation. On the other hand, the classical selectivity shows no significant change with temperature and remains almost constant (0.91 − 0.81), slightly favoring 3 He transmission. This is valid for both membranes.
As expected, with the rise in temperature, the classical selectivity follows the RPMD and WP3D ones.

Concluding remarks
In this work, we have calculated the thermal rate coefficient for the transmission of He isotopes through the pores of one atom thick graphdiyne (Gr2) and graphtriyne (Gr3) membranes using the ring polymer molecular dynamics (RPMD) method within the temperature range 20 K − 250 K. Transmission through Gr2 has a substantial free energy barrier even at 20 K, and this barrier height increases with increasing temperature. On the other hand, transmission through Gr3 can either have marginally negative (up to 30 K) or small positive (T≥ 50 K) free energy barrier. The extent of the barrier height directly impacts the calculated rates, as evident from the fact that the rate coefficient on Gr3 is at least 10 17 order of magnitude greater than on Gr2 at 20 K. The rate coefficient on Gr2 increases rapidly with temperature and starts to converge, starting from 150 K (∼ 10 8 s −1 ). However, the rate coefficients on Gr3 do not vary appreciably with temperature and remain almost constant for the whole temperature range (3.63×10 10 s −1 − 1.01×10 11 s −1 ). In general, the rate coefficient calculated for the transmission through Gr3 is always greater than the corresponding one on Gr2 over the whole temperature regime considered in this work. Moreover, we found that the recrossing dynamics have little or no effect on the final value of the rate coefficients.
The selectivity ratio, which indicates the preference of either isotope for its permeation through the membranes, has been calculated as a function of temperature.

Supporting Information Available
The following files are available free of charge.
• supplementary information file: computational details and additional results.