Article pubs.acs.org/JPCA
Helium Isotope Enrichment by Resonant Tunneling through Nanoporous Graphene Bilayers Salvatore Mandrà,*,†,‡ Joshua Schrier,*,§ and Michele Ceotto*,⊥ †
Department of Physics, Università degli Studi di Milano, via Celoria 16, 20133 Milano, Italy Department of Chemistry and Chemical Biology, Harvard University, 12 Oxford Street, Cambridge, Massachusetts 02138, United States § Department of Chemistry, Haverford College, Haverford, Pennsylvania 19041, United States ⊥ Department of Chemistry, Università degli Studi di Milano, via Golgi 19, 20133 Milano, Italy ‡
S Supporting Information *
ABSTRACT: Graphene is impermeable to gases, but introducing subnanometer pores can allow for selective gas separation. Because graphene is only one atom thick, tunneling can play an important role, especially for low-mass gases such as helium, and this has been proposed as a means of separating 3 He from 4He. In this paper, we consider the possibility of utilizing resonant tunneling of helium isotopes through nanoporous graphene bilayers. Using a model potential fit to previously reported DFT potential energy surfaces, we calculate the thermal rate constant as a function of interlayer separation using a recently described time-independent method for arbitrary multibarrier potentials. Resonant transmission allows for the total flux rate of 3He to remain the same as the best-known single-barrier pores but doubles the selectivity with respect to 4He when the optimal interlayer spacing of 4.6 Å is used. The high flux rate and selectivity are robust against variations of the interlayer spacing and asymmetries in the potential that may occur in experiment.
1. INTRODUCTION Two stable isotopes of helium exist. The predominant isotope on earth is 4He, which is produced as uranium α-decays into lead. Very low concentrations of 3He can be found in various terrestrial sources (remaining from the accretion of the earth) and lunar regolith (deposited by solar wind), but currently the only viable source is from the decay of tritium, as a byproduct of nuclear weapons stockpiles.1 The dwindling supply of 3He, combined with the growing use of this isotope for neutrondetection sensors for nuclear security monitoring as well as basic research and hydrocarbon explorationhas raised concerns over the long-term supply of this isotope.2 A vast quantity of 3He could potentially be separated from existing helium stockpiles, natural gas, or the atmosphere, but the energy costs of performing the separation make these exceedingly impractical with current cryogenic distillation methods that require liquefying the gases.3 Membrane-based methods of gas separation are potentially more efficient because they avoid the energy cost of liquefying the gases.4 A new research direction is the use of subnanometer sized pores introduced into graphene as a high-performance separation membrane. The theoretical simulations of gas separation5,6 and experimental synthesis7,8 of nanoporous graphene materials have been the subject of several recent © 2014 American Chemical Society
reviews. Because these membranes are only one atom thick, and the pores are so small, quantum tunneling can play a significant role in the transmission of very light elements such as helium atoms, even at room temperature.9 The mass dependence of tunneling provides a potential means of separating the two isotopes under conditions where no classical separation occurs.10 Because the concentrations of 3He are so low, it is important to have a very high selectivity for the transmission of this isotope if one is to create a practical process. Certain nitrogen-functionalized pores in graphene11,12 and two-dimensional polymers13 have been predicted to have both high flux and high selectivity for helium isotope separation. The highest isotope selectivity ratio of 19 is achieved only at 10 K and rapidly approaches no selectivity by 25 K, but at the lowest temperature the flux rate is only about 10−9 mol cm−2 s−1 under a feed pressure of 1 bar. The goal of this paper is to develop nanostructures with higher selectivity and flux rate, capable of operating at higher temperatures, which will further lower the cost of obtaining 3He from natural sources. Special Issue: Franco Gianturco Festschrift Received: March 14, 2014 Revised: May 12, 2014 Published: May 22, 2014 6457
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
All of the previous studies of helium isotope separation by tunneling have considered a single barriera single nanoporous graphene or graphene-like sheet. In this paper, we investigate resonant tunneling through double-barrier nanostructures as a way to perform isotope separation. Resonant tunneling through double barriers is a textbook problem in quantum mechanics,14 and resonant tunneling of electrons is widely used for semiconductor heterostructure electronics devices.15 Unlike single-barrier tunneling, complete transmission occurs when the incoming particles have kinetic energies matching the quasi-bound states of the potential well formed by the double barrier. These resonant transmissions occur at energies much lower than the barrier height. Due to these resonant peaks the transmission can be high, even when the temperature is so low that few atoms have sufficient kinetic energy to cross the potential energy barrier classically. What we are proposing is not to be confused with the resonant tunneling of electrons within the π-bands of a graphene sheet where the potential has been modulated to form a double barrier.16 Rather, we propose transmission of a helium atom through two aligned pores of two nanoporous graphene sheets; the double barrier arises due to the nature of the potential energy surface for performing this atom transmission. In the following section we describe a model potential based on previous ab initio potential energy calculations. From a theoretical perspective, the calculation of the thermal rate constants for resonant tunneling of atoms is numerically challenging, and we describe how to perform this calculation in section 3. In section 4 we explore the optimal separation between the two barriers and the effect of van der Waal’s interactions on the isotope separation. Finally, section 5 concludes the paper.
Figure 1. Single-barrier (blue dots) and double-barrier (red line) potentials at different distances, Δ, between the peaks. As described in section 1, each potential is divided into three parts: (I) and (III) the PES is the single-barrier potential VS(x) calculated by Hauser et al.;11 (II) the PES is given by the scaled sum of two VS(x) in a way that the minimum between the two barriers maintains the same depth as for the single barrier.
2. A MODEL POTENTIAL FOR NANOPOROUS GRAPHENE MULTILAYERS The global potential energy surface (PES) of an atom approaching a two-dimensional nanoporous surface is very complex if taken in its full dimensionality. For simplicity, we consider only the minimum energy path, where the atom approaches the center of the hole in a perpendicular direction to the graphene surface, turning this into a one-dimensional problem. This is reasonable because regions of the graphene surface without a hole have a very large potential energy barrier and are thus impermeable to helium atoms.17 Both the repulsive barrier and the attractive wells are important for efficient quantum isotope separation.12 Potential wells are important for keeping the flux high by attracting the isotopes to the membrane, and the barrier height and width are important for tunneling selection (Supporting Information). The zeropoint energy and tunneling contribute in opposite directions with respect to the temperature gradient.10 At low temperatures, tunneling favors lighter isotopes whereas the zero-point energy differences of the transition state favor the heavier ones; at high temperatures, just the opposite behavior occurs. Previously reported exploratory calculations using a model PES found the optimal conditions for helium isotope separation to be a barrier height of a few hundred wavenumbers, surrounded by two minima and operating between 10 and 20 K.12 This can be achieved by creating a two-ring defect in a graphene sheet and then nitrogen doping one side of this defect;11,12 this PES is shown as the blue dots in Figure 1. We use this single-barrier PES as a model for investigating the double-barrier system. In principle, experimental approaches can control the spacing between graphene18 and
graphene oxide19 sheets (for the purposes of hydrogen storage) to yield many possible interlayer spacings. Our goal is to explore a wide range of spacings between the layers and identify the best target. Because it is too computationally demanding to obtain a high-quality ab initio PES for each of these systems, we instead use a simplified model to expose the essential aspects of this problem. The continuous red lines in Figure 1 show the double-barrier potentials offset by different barrier distances, Δ. The potential is constructed by dividing the minimum energy path into three regions, called I, II, and III in Figure 1, defined by the changing of sign of the potential in between the two barriers. In regions I and III the PES is the same as the original single-barrier potential, VS(x). This is reasonable because the transiting atom primarily interacts with the nearest membrane. The potential in region II (between the two graphene sheets) may have minima that are greater than or equal to the minima of the single-barrier case. Molecular dynamics simulations and popular empirical dispersion corrections applied to density functional theory (e.g., by Grimme20) typically assume that dispersion behaves as a pairwise sum of atom−atom interactions. However, random phase approximation21 and quantum Monte Carlo22 calculations have shown that this simple additive model is not applicable for interactions involving extended two-dimensional systems such as the graphene sheets studied here. As we will show below, the 6458
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
inconvenient for resonant tunneling calculations.35 Hence, a time-independent approach is preferred for the calculation of the thermal rate constant36−49
isotope selectivity increases with the number of quasi-bound states within this potential, which is proportional to the depth of the potential. Considering the limiting case where the potential in region II is rescaled to have the same minima depth as the single-barrier case places a conservative lower bound on the isotope selectivity. However, this assumption is not necessary for the computational method used to compute the transmission, nor for the resonant tunneling separation effect described below. Finally, for Δ larger than the support of the single-barrier potential VS(x), the double-barrier potential is well approximated by the following equation VD(x ;Δ) ≈ VS(x − Δ/2) + VS(x + Δ/2)
k(T ) =
+∞
1 2π ℏQ R (T )
∫−∞
e − βE T ( E ) d E (2)
where QR(T) is the reactants partition function and T (E ) =
1 (2π ℏ)2 Tr[F1̂δ(E − Ĥ ) F2̂ δ(E − Ĥ )] 2
(3)
is the cumulative reaction probability evaluated between the dividing surface f1(s) = 0 and f 2(s) = 0 located at s = 0. F̂ is the flux operator and Ĥ is the Hamiltonian operator.23 To be consistent with commonly used notations and to avoid confusion, T (with no argument specified) refers to temperature, and T(E) refers to the transmission probability as a function of kinetic energy. We will always include the argument in the transmission probability. In a recent paper,35 we presented a time-independent method for calculating the thermal rate constant, k(T), for any arbitrary multibarrier potential. The method successfully treats double-barrier potentials with many quasi-bound states. It is based on a fast and robust procedure for calculating k(T) as the thermal average of the transmission probability T(E) according to eq 2. It solves the Schrödinger equation as an ordinary differential equation (ODE) with the energy as a parameter, thus avoiding the related eigenvalue problem. The main idea is to force the asymptotic eigenfunctions ψp(x) of equation
|x | ≫ Δ (1)
Small variations in the PES modify the barrier height and break the symmetry of the double barrier and can potentially ruin the resonant tunneling effect. These variations might arise from defects in the pores, thermal fluctuations of pore atoms, or angled trajectories through the double-barrier structure. To study this possibility, we also considered the case where one barrier is 20% higher than the other, as shown in Figure 2. We discuss the results in the following section.
3. A TIME-INDEPENDENT METHOD FOR RESONANT TUNNELING Time-dependent theories typically calculate thermal rates by integrating a flux autocorrelation function,23−34 but this is
⎡ ℏ2 d2 ⎤ Ĥ ψp(x) = ⎢ − + V (x)⎥ψp(x) = Eψp(x) 2 ⎣ 2m dx ⎦
(4)
at fixed eigenvalue E = p /2m to satisfy the planewave boundary conditions at distances, L, that are sufficiently far away from the scattering potential, 2
⎧ e ipx / ℏ x ≪ −L ⎪ ψp(x) = ⎨ t ipx / ℏ r −ipx / ℏ + e x ≫ +L ⎪ 2e ⎩ |t | |t |
(5)
These are obtained by taking the boundary conditions lim ψp(x) = 1
(6a)
lim ∂xψp(x) = ip /ℏ
(6b)
x →−L
x →−L
where x is in the asymptotic region. Finally, one obtains expressions for the transmission and reflection probabilities, T (E) = |t |2 = lim 4 ψp(x) + x→L
ℏ ∂xψ (x) ip p
−2
(7a)
2 ⎞−1 ⎛ 1 ℏ ⎜ R(E) = |r | = lim 1 − ⎜1 + ψp(x) − ∂xψp(x) ⎟⎟ x→L 4 ip ⎝ ⎠ 2
(7b)
where x is an asymptotic position after the barrier passage. The advantage of this method over the time-dependent approach is that it allows for a fine energy grid when the transmission probability is calculated, which is crucial for accurately integrating eq 2 when there are very narrow resonant
Figure 2. Potential energy surfaces for a model double nanoporous graphene pair of sheets. The potential is not symmetric for a more realistic modeling. Notation as in Figure 1. 6459
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
Figure 3. 3He and 4He transmission probabilities T(E) for the potential profile of Figure 1. Left panel: lin−lin scale. Right panel: log−log scale. All of the distances considered in this paper are between the two extreme distances Δ = 8.0 au and Δ = 12.0 au, shown here.
significantly the thermal rate at very low temperatures.35 The low energy portion of the log−log plot (E < 0.0001 eV) also shows that T(E) cannot be safely approximated as zero at these low energies and is different for the two isotopes. Although this is typically ignored at high temperatures, it is too important at low temperatures. For these reasons, the lower energy bound has been fixed at Emin = 10−15 au and the lower bound of the transmission probability fixed at T(E)min = 10−40 au. A further inspection of the top-left panel of Figure 3 shows that there are eight resonant peaks for the 4He transmission probability and only seven for the 3He. This difference is more evident as the graphene sheets are separated, as shown in the bottom left panel. Each of the resonant transmissions corresponds to a quasi-bound (metastable) state of the atom inside the double-barrier potential. Approximating this region as a harmonic potential, then the eigenvalues spacing is inversely proportional to the square root of the oscillator (i.e., the helium atom) mass. As a result, the 4He peaks occur at lower energies. Because the potential has a finite height, there are only a limited number of these quasi-bound states, and because the 4He energies are lower, there are more quasi-bound states for this isotope. However, in the Δ = 8.0 au case, the 3He resonant peak is the lower one because there are no quasibound states with energies less than the original asymptotic energy (explained below). Boltzmann-weighted energy integration of T(E) as in eq 2 yields the thermal rate constants, k(T). In Figure 4 we show the product of the thermal rate constant with the partition function, i.e., k(T) QR(T), to more clearly show the tunneling effects. The direct comparison with the single barrier can still be done because, for one-dimensional systems, partition functions are independent of the asymptotic shape of PES. In the same figure, the single barrier (violet left-triangles; the topmost line in both graphs) and the transition state theory (TST) (dashed black line; the bottom-most line in both graphs) rates are reported. On the upper panel the rates are for 3He passage and on the lower for 4He. Even though the thermal rate for the single barrier is typically greater than that of the double-barrier,
transmission peaks. In this case, the method allows for using a coarse grid in flat regions and a fine one in the regions around resonant peaks. This improves the numerically convergence for the energy integration in eq 2. Details about the numerical convergence are reported in the Supporting Information. A more formal presentation can be found in ref 35.
4. RESULTS AND DISCUSSION In this section, results on the transmission probability T(E) and the thermal rate constant k(T) are presented together with He isotope selection considerations. The numerical convergence of these calculations is challenging: Extra precision is required to calculate the rate at the low temperatures (5−30 K) needed for isotope separation. Integrating the contribution of the resonant peaks requires calculating T(E) for very low energies. Thermal Rate Constants. Figure 3 shows the transmission probabilities for two resonant tunneling cases, demonstrating the presence of several resonance peaks at kinetic energies below 0.022 eV. To demonstrate that our computational approach can treat any arbitrary barrier distance, here we show the results for the most closely (Δ = 8.0 au, top panels) and distantly (Δ = 12.0 au, bottom panels) spaced graphene sheets that we considered; all other calculations in this paper are intermediate between these two extremes. In the left panels of Figure 3 both axes are in linear−linear scale, and the right panels are for the same data but in log−log scale. The linear scale representation shows that each resonant peak is very narrow and has a maximum transmission of 1, in agreement with the symmetric double-barrier potential resonant tunneling rule. The narrowness of the peaks presents a numerical challenge for the energy integration of eq 2. For this purpose, we employed both a coarse-grained grid, as explained above, and a resonant peak transmission splined profile as suggested in refs 35 and 50. Thus, the energy grid spacing is of the order of 10−15 au. The log−log plots on the right panels of Figure 3 show that the chosen integration energy range includes all resonant peaks. This is important because exclusion of even one of the low-energy resonant transmission peaks can change 6460
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
Figure 4. Thermal rate constant tunneling effects. The product of k(T) QR(T) is plotted for 3He (upper panel) and 4He (lower panel). Transition state theory (TST) results in the dashed black line and the single barrier in the purple line (with the largest value in both insets). Other lines are for different barrier distances, Δ.
Figure 5. Thermal rate constant ratio, k3/k4, where k3 and k4 are the thermal rate constants for 3He and 4He, respectively. The upper panel figure shows the temperature (T) dependence of the thermal rate constant ratio; the horizontal black dashed line shows the isotopeindependent classical transition state theory (TST) value. The bottom panel shows the interbarrier distance (Δ) dependence of the thermal rate ratio; the horizontal black dashed lines show the single-barrier result at each temperatures.
order-of-magnitude discrepancies occur only at very low temperatures (T < 10 K); at higher temperatures (T > 10 K) the thermal rate constants of single and double barriers are comparable. The deviation of the thermal rate constants from the TST result show that transmission through both the singleand double-barrier graphene nanopores is greatly increased by quantum tunneling effects. This figure also shows nonlinear behavior of the thermal rate for the double-barrier potential as a function of the interbarrier distance, Δ. For 3He, shown in the top panel, the order of k(T) QR(T) as a function of Δ is 10 au > 8 au ≈ 11 au > 12 au > 9 au. For 4He, shown in the bottom panel, the order of k(T) QR(T) as a function of Δ is 9 au > 11 au > 10 au ≈ 12 au > 8 au. As we will discuss below, this is a direct consequence of the variation of the number and shape of the resonant tunneling peaks with the distance Δ and the masses of the particles. To see if the double porous graphene sheets can better select 3 He from 4He than the single sheet, we show the ratio of the thermal rate constants for the two isotopes, k3/k4, in the upper panel of Figure 5. In a classical approximation no quantum tunneling occurs, and if one also neglects the different (massdependent) zero-point energy contributions to the potential energy barrier, then the TST ratio is temperature-independent and given by the ratio of the partition functions. In contrast, for the deep-resonant tunneling regime k3/k4 varies significantly at small temperatures and gradually converges to the single-barrier ratio at large temperature (T > 15 K). This is due to the change in the number of resonant tunneling peaks for 3He versus 4He. Even if the rate ratio is the same for the single and the double barrier at larger Δ distance, the rate value is different. Interestingly, for sufficiently low temperature, the ratio for the double-barrier potential can be smaller, comparable to, or even larger than the ratio for the single barrier. This effect is
exponentially enhanced when the temperature is further lowered. The double-barrier resonant tunneling effect is so pronounced that for temperatures less than T < 15 K, an inverse isotope effect is observed for Δ = 9.0 au and Δ = 9.4 au To better appreciate the He isotope enrichment effects by resonant tunneling, the lower panel of Figure 5 shows the ratio k3/k4 for the double-barrier potential as a function of the interbarrier distance Δ. For reference, the (Δ-independent) single-barrier k3/k4 ratio is shown as horizontal dashed lines for each temperature. For temperatures ≥15 K, no significant variations in the isotope enrichment between the single- and double-barrier cases is observed. At lower temperatures, the double-barrier system offers the possibility of improved 3He selectivity only for certain values of Δ. For example, at 10 K, the double-barrier potential can select the lighter isotope with more than twice the selectivity of single-barrier potential when Δ = 8.72 au. However, a small increase in Δ can reduce the isotope selectivity below that of the single barrier. A key result of this paper is that for temperatures below 12 K, an appropriate choice of the distance Δ between the double porous graphene sheets can substantially improve selection compared to the single graphene layer. Isotope Selection versus Barriers Distances and Potential Shape. To better understand the mechanism of enhancing the isotope selection induced by the He resonant tunneling, we count the number of resonant transmission peaks in T(E), as a function of Δ. As shown in Figure 6, the number of resonant peaks for 4He is always greater-than-or-equal-to that of 3He, but the number of resonant peaks is not a monotonically increasing step function. Even though the 6461
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
potential value (see PES in Figure 1), energies of lower metastable states can also be negative. However, only the metastable states with positive energy contribute to the transmission probability T(E), because the energies of incoming helium atoms are always positive. Therefore, even if the total number of metastable states with both positive and negative energies grows by increasing the distance between the graphene sheets, states with negative energy are trapped and do not contribute to the transmission. The resonant peaks shifts to lower energies as Δ is increased, until one gets to a case with the closed channels, shown as the orange lines in Figure 7. The sudden closure of a channel is responsible for the number of peaks decreasing observed in Figure 6 for both isotopes. In conclusion, there are two opposite trends responsible for the number of resonant peaks. The positive trend is due to the increasing of anharmonicity as Δ is increased. The negative trend is due to the metastable state energy switching from positive to negative values. We stress that similar results occur if the van der Waals well between the two graphene sheets (region II) is deeper. At low temperature, the isotope separation is mainly driven by the difference of the number of resonant peaks between 3He and 4 He (as shown in Figures 5 and 6), but a deeper central well primarily increases the number of negative energy resonant peaks, which do not contribute to the transmission. Because 4 He is heavier, more negative energy peaks for 4He are present for a deeper central well which will increase the 3He/4He ratio. Asymmetric Potential Model. In any real experiment, imperfect alignment of the pores, variations in the pore geometries, and thermal fluctuations will all distort the ideal PES. When the double barrier is no longer symmetrical, the resonant transmission peaks no longer have unit transmission probabilities. To study whether these distorted potentials are still capable of isotope selection and of high flux, we examined a set of model asymmetric PES shown in Figure 2, where one barrier is scaled to lower values. The transmission probability for both isotopes is shown in Figure 8 for the smaller (upper panel) and the larger (lower panel) interbarrier distances considered. As expected, the T(E) peaks are no longer equal to unity because of the asymmetry in the potential profile. The position of the T(E) peaks changes more with increasing Δ for the heavier isotope, similar to what was observed for the symmetric potential. Does asymmetry change the thermal rate values, k(T) QR(T)? Yes. As shown in Figure 9, all of the asymmetric PESs actually have a higher rate than the corresponding symmetric PESs for both isotopes, for any interbarrier distance Δ. This is of practical importance, because the higher rates mean a faster flux of atoms across the membrane. Does asymmetry change the isotope selectivity? No. As shown in Figure 10, the thermal rate has a ratio k3/k4 similar to the symmetric one for different temperatures and distances Δ. (The ratios for the symmetric potential, previously shown in Figure 5, are indicated here as open symbols, so as to enable a direct comparison against the asymmetrical PES results.) The isotope selectivity is preservedand in some cases slightly enhancedfor the asymmetric potential. This is very encouraging for future experiments, as it means that the resonant tunneling effect is robust against these types of small fluctuations, and suggests that trajectories deviating by a small angle or displacement from the one studied here will have the same qualitative behavior. The behavior of the resonant enhancements as a function of Δ is preserved because the
Figure 6. Comparison between the number of resonant peaks for 4He and 3He. Vertical lines are in correspondence with significant k3/k4 ratio changes in Figure 5. As one can see, variations of k3/k4 and of the number of peaks are only partially correlated.
difference of the number of resonant peaks and the discontinuities of the ratio k3/k4 have a common trend, there is no direct correspondence between them. This is because k(T) is not only the sum of the resonant peaks but it also depends on the shape of the transmission probability T(E). To better understand the physical reasons why the number of peaks decrease for some Δ values, we focus on the log−log plot of T(E) for the lowest energy scattering region. The upper panel of Figure 7 reports T(E) for different values of Δ for 3He
Figure 7. Transmission probability T(E) for 3He (top panel) and 4He (bottom panel) at low scattering energies. Different curves correspond to different barrier distances Δ. At large enough distance, the lower resonant peak is suppressed because the corresponding metastable state energy is below the asymptotic one.
and the lower panel shows the 4He results. (The specific set of Δ values are different in the two insets, to better illustrate the varying number and shape of the peaks for the two isotopes.) As Δ is increased, peaks shift to lower energy due to the increment of anharmonicity caused by the van der Waals interaction between He and the graphene sheets. Because the van der Waals well minimum is below the zero asymptotic 6462
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
Figure 10. Thermal rate constant ratios and isotope selection as a function of Δ. Open symbols and lines show the symmetric and asymmetric PES results, respectively. Lower temperatures and the asymmetric potential of Figure 2 offer better isotope selection.
metastable state energies located at the bottom of the intrabarrier well are not modified. Because these are the states mainly responsible for the thermal rate values at low temperatures, the optimal value of Δ for resonant tunneling are nearly the same for both the symmetric and asymmetric potentials. Feasibility. Can these double-barrier nanostructures purify a macroscopic amount of 3He? We can obtain an upper-bound estimate of the helium flux by making two assumptions. First, we assume that any incident particle on the two-dimensional surface is transmitted through the potential associated with a pore. In reality, the pores will be spread over a larger distance, with a density determined by the synthesis procedure. Second, we assume that all trajectories are normal to the graphene membranes. Angled trajectories are neglected as they are typically reflected by the potential or, if the angle is very small, will introduce a small asymmetry in the potential that does not qualitatively change the results. Both of these assumptions increase the predicted flux rate by a multiplicative factor. As shown in Figure 5, the isotope separation for a double graphene layer is maximized when the distance between the two graphene sheets is Δ = 8.72 au. Assuming a pressure of 1 atm, the daily amounts of 3He atoms that cross a square-meter of a single-layer or a double-layer membrane at T = 10 K are almost identical and equal to ≈800 g per day. However, as shown in Figure 5, the output from the double graphene layer will be twice as pure as the single layer; i.e., the permeate will contain half as much 4He. Due to resonant tunneling a double graphene layer is more efficient than a single graphene layer for this separation. Can this double-barrier structure be made? Yes. In principle, experimental approaches can control the spacing between graphene18 and graphene oxide19 sheets to yield the optimal interlayer spacing of 8.7 au (4.6 Å) predicted above. One possible strategy is to intercalate molecules between the graphene layers. Our optimal interlayer spacing is larger than the 3.8 Å spacing achieved by intercalating methane molecules,51 but only slightly smaller than the 4.7 Å spacing achieved by intercalating tetrabutylammonium cations,52 which suggests that intermediate values can be achieved by using different types of molecules. Another possibility is to intercalate metal atoms between the graphene layers. Our optimal interlayer spacing is larger than the 3.7 Å interlayer spacing for lithium-intercalated graphene,53 but smaller than the 6.5 Å spacing for Au6 intercalated graphene,54 again suggesting that intermediate values can be achieved by intercalating different
Figure 8. Transmission probabilities for the asymmetric potential setup of Figure 2. The scale is linear for both axes. Upper panel for a distance of Δ = 8 au and lower panel for Δ = 12.0 au. The resonant peaks transmission probability is not unity because of the asymmetry in the potential.
Figure 9. k(T) QR(T) values comparison between the symmetric (SYM) and the asymmetric (ASYM) PES of Figure 4. Upper panel for 3 He and lower one for 4He. Dashed line for transition state theory (TST) values. Asymmetric potential rates are always bigger than symmetric ones at same barrier distances Δ.
asymmetry is sufficiently small that it does not affect the resonant peaks in the deep tunneling regime. Consequently, the 6463
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
■
types of metals. The pores could be simultaneously created and decorated with nitrogen by means of nitrogen ion implantation.55,56 Our analysis of the asymmetrical potential shows that typical variations will not be detrimental to performance.
ASSOCIATED CONTENT
S Supporting Information *
Numerical convergence analysis for the tunneling calculations and derivation of the separation rate. This material is available free of charge via the Internet at http://pubs.acs.org.
■
REFERENCES
(1) Wittenberg, L.; Cameron, E.; Kulcinski, G.; Ott, S.; Santarius, J.; Sviatoslavsky, G.; Sviatoslavsky, I.; Thompson, H. A Review Of 3He Resources And Acquisition For Use As Fusion Fuel. Fusion Technology 1992, 21, 2230−2253. (2) Managing Critical Isotopes: Weaknesses In DOE’s Management Of Helium-3 Delayed The Federal Response To A Critical Supply Shortage. GAO-11-472, U.S. Government Accountability Office, 2011. (3) Shea, D. A.; Morgan, D. Helium-3 Shortage: Supply, Demand, And Options For Congress; Congressional Research Service: Washington, DC, 2010 (4) Bernardo, P.; Drioli, E.; Golemme, G. Membrane Gas Separation: A Review/State Of The Art. Ind. Eng. Chem. Res. 2009, 48, 4638− 4663. (5) Jiao, Y.; Du, A.; Hankel, M.; Smith, S. C. Modelling Carbon Membranes For Gas And Isotope Separation. Phys. Chem. Chem. Phys. 2013, 15, 4832−4843. (6) Jiao, Y.; Hankel, M.; Du, A.; Smith, S. C. Porous Graphene And Nanomeshes. Graphene Chemistry: Theoretical Perspectives 2013, 129− 151. (7) Russo, P.; Hu, A.; Compagnini, G. Synthesis, Properties And Potential Applications Of Porous Graphene: A Review. Nano-Micro Lett. 2013, 5, 260−273. (8) Jiang, L.; Fan, Z. Design Of Advanced Porous Graphene Materials: From Graphene Nanomesh To 3D Architectures. Nanoscale 2014, 6, 1922−1945. (9) Schrier, J. Helium Separation Using Porous Graphene Membranes. J. Phys. Chem. Lett. 2010, 1, 2284−2287. (10) Schrier, J.; McClain, J. Thermally-Driven Isotope Separation Across Nanoporous Graphene. Chem. Phys. Lett. 2012, 521, 118−124. (11) Hauser, A. W.; Schwerdtfeger, P. Nanoporous Graphene Membranes For Efficient 3He/4He Separation. J. Phys. Chem. Lett. 2012, 3, 209−213. (12) Hauser, A. W.; Schrier, J.; Schwerdtfeger, P. Helium Tunneling Through Nitrogen-Functionalized Graphene Pores: Pressure-And Temperature-Driven Approaches To Isotope Separation. J. Phys. Chem. C 2012, 116, 10819−10827. (13) Brockway, A. M.; Schrier, J. Noble Gas Separation Using PG-ES X (X= 1, 2, 3) Nanoporous Two-Dimensional Polymers. J. Phys. Chem. C 2012, 117, 393−402. (14) Bohm, D. Quantum Theory; Prentice-Hall Physics Series; Prentice Hall: Englewood Cliffs, NJ, 1951 (15) Tsu, R.; Esaki, L. Tunneling In A Finite Superlattice. Appl. Phys. Lett. 2003, 22, 562−564. (16) Pereira, J. M.; Vasilopoulos, P.; Peeters, F. Graphene-Based Resonant-Tunneling Structures. Appl. Phys. Lett. 2007, 90, 132122− 132122. (17) Leenaerts, O.; Partoens, B.; Peeters, F. Graphene: A Perfect Nanoballoon. Appl. Phys. Lett. 2008, 93, 193107. (18) Jin, Z.; Lu, W.; O’Neill, K. J.; Parilla, P. A.; Simpson, L. J.; Kittrell, C.; Tour, J. M. Nano-Engineered Spacing In Graphene Sheets For Hydrogen Storage. Chem. Mater. 2011, 23, 923−925. (19) Kim, B. H.; Hong, W. G.; Moon, H. R.; Lee, S. M.; Kim, J. M.; Kang, S.; Jun, Y.; Kim, H. J. Investigation On The Existence Of Optimum Interlayer Distance For H2 Uptake Using Pillared-Graphene Oxide. Int. J. Hydrogen Energy 2012, 37, 14217−14222. (20) Grimme, S. Semiempirical GGA-Type Density Functional Constructed With A Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (21) Dobson, J. F.; White, A.; Rubio, A. Asymptotics Of The Dispersion Interaction: Analytic Benchmarks For Van Der Waals Energy Functionals. Phys. Rev. Lett. 2006, 96, 073201. (22) Drummond, N.; Needs, R. Van Der Waals Interactions Between Thin Metallic Wires And Layers. Phys. Rev. Lett. 2007, 99, 166401. (23) Miller, W. H.; Schwartz, S. D.; Tromp, J. W. Quantum Mechanical Rate Constants For Bimolecular Reactions. J. Chem. Phys. 1983, 79, 4889−4898. (24) Sklarz, T.; Kay, K. Semiclassical Initial Value Treatment Of Correlation Functions. J. Chem. Phys. 2004, 120, 2606−2617.
5. CONCLUSIONS We have applied a recently developed time-independent method for thermal rate calculations of multibarrier potentials to demonstrate the feasibility of isotope separation using resonant tunneling through bilayer graphene nanopores. From a fundamental perspective, this extends a basic quantum mechanical phenomenon to heavy atoms. From the perspective of chemical engineering, this is the first application of resonant tunneling to separations and it makes a significant improvement to the helium isotope separation. Considering that higher flux means lower selectivity, and vice versa,57 there is a trade-off between these two trends. However, at a working temperature of 10 K we can achieve the same 3He flux as the best single barriers but double the selectivity with respect to 4He. The high flux rate and isotope selectivity are robust against asymmetries in the potential. Moreover, the optimal separation between the two graphene sheets is intermediate to known results for molecular and metal intercalation compounds of graphene. We hope that the significant improvement in isotope separation will motivate future computational and experimental work on bilayer graphene nanopore systems. Future investigations will consider ab initio calculations of the two layer graphene sheets potential in more than one dimension, multidimensional tunneling, and the effects of different size pores.
■
Article
AUTHOR INFORMATION
Corresponding Authors
*S. Mandrà: e-mail, [email protected]. *J. Schrier: e-mail, [email protected]. *M. Ceotto: e-mail, [email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We thank Geoffrey Martin-Noble, Anahita Nourmahnad, and Kylen Solvik for careful readings of the manuscript. S.M. and M.C acknowledge the CINECA and the Regione Lombardia award under the LISA initiative, for the availability of high performance computing resources and support. S.M. also acknowledges the Air Force Office of Scientific Research (prime sponsor) and the University of California, San Diego, for the partial support to this work under the grants FA955012-1-0046 and 10323836-SUB. J.S. was supported by the Research Corporation for Scientific Advancement’s Cottrell Scholar grant and used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. 6464
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
(47) Zimmermann, T.; Vaníček, J. Three Applications Of Path Integrals: Equilibrium And Kinetic Isotope Effects, And The Temperature Dependence Of The Rate Constant Of The [1, 5] Sigmatropic Hydrogen Shift In (Z)-1, 3-Pentadiene. J. Mol. Model. 2010, 16, 1779−1787. (48) Buchowiecki, M.; Vaníček, J. Direct Evaluation Of The Temperature Dependence Of The Rate Constant Based On The Quantum Instanton Approximation. J. Chem. Phys. 2010, 132, 194106. (49) Vaníček, J. Beyond Transition State Theory: Accurate Description Of Nuclear Quantum Effects On The Rate And Equilibrium Constants Of Chemical Reactions Using Feynman Path Integrals. CHIMIA Int. J. Chem. 2011, 65, 715−719. (50) Ryaboy, V.; Lefebvre, R. Flux-Flux Correlation Function Study Of Resonance Effects In Reactive Collision. J. Chem. Phys. 1993, 99, 9547−9552. (51) Huang, Q.; Chen, X.; Lin, J.; Li, K.; Jia, Y.; Liu, J.; Guo, L.; Wang, W.; Wang, G. Preparation Of Quasi-Free-Standing Graphene With A Super Large Interlayer Distance By Methane Intercalation. J. Phys. Chem. C 2011, 115, 20538−20545. (52) Sirisaksoontorn, W.; Adenuga, A. A.; Remcho, V. T.; Lerner, M. M. Preparation And Characterization Of A Tetrabutylammonium Graphite Intercalation Compound. J. Am. Chem. Soc. 2011, 133, 12436−12438. (53) Qi, Y.; Guo, H.; Hector, L. G.; Timmons, A. Threefold Increase In The Young’s Modulus Of Graphite Negative Electrode During Lithium Intercalation. J. Electrochem. Soc. 2010, 157, A558−A566. (54) Sapkota, I.; Roundtree, M. A.; Hall, J. H.; Wang, X.-Q. Tunable Band Gap In Gold Intercalated Graphene. Phys. Chem. Chem. Phys. 2012, 14, 15991−15994. (55) Guo, B.; Liu, Q.; Chen, E.; Zhu, H.; Fang, L.; Gong, J. R. Controllable N-Doping Of Graphene. Nano Lett. 2010, 10, 4975− 4980. (56) Åhlgren, E.; Kotakoski, J.; Krasheninnikov, A. Atomistic Simulations Of The Implantation Of Low-Energy Boron And Nitrogen Ions Into Graphene. Phys. Rev. B 2011, 83, 115424. (57) Robeson, L. M. The Upper Bound Revisited. J. Membr. Sci. 2008, 320, 390−400.
(25) Pollak, E.; Martin-Fierro, E. New Coherent State Representation For The Imaginary Time Propagator With Applications To ForwardBackward Semiclassical Initial Value Representations Of Correlation Functions. J. Chem. Phys. 2007, 126, 164107. (26) Martin-Fierro, E.; Pollak, E. Forward-Backward Semiclassical Initial Value Series Representation Of Quantum Correlation Functions. J. Chem. Phys. 2006, 125, 164104. (27) Liu, J.; Miller, W. H. Real Time Correlation Function In A Single Phase Space Integral Beyond The Linearized Semiclassical Initial Value Representation. J. Chem. Phys. 2007, 126, 234110. (28) Yamamoto, T.; Wang, H.; Miller, W. H. Combining Semiclassical Time Evolution And Quantum Boltzmann Operator To Evaluate Reactive Flux Correlation Function For Thermal Rate Constants Of Complex Systems. J. Chem. Phys. 2002, 116, 7335−7349. (29) Tromp, J. W.; Miller, W. H. The Reactive Flux Correlation Function For Collinear Reactions H+ H2, Cl+ HCl And F+ H2. Faraday Discuss. Chem. Soc. 1987, 84, 441−453. (30) Bonella, S.; Monteferrante, M.; Pierleoni, C.; Ciccotti, G. Path Integral Based Calculations Of Symmetrized Time Correlation Functions. I. J. Chem. Phys. 2010, 133, 164104. (31) Shi, Q.; Geva, E. Nonradiative Electronic Relaxation Rate Constants From Approximations Based On Linearizing The PathIntegral Forward-Backward Action. J. Phys. Chem. A 2004, 108, 6109− 6116. (32) Zhu, W.; Zhao, Y. Effects Of Anharmonicity On DiffusiveControlled Symmetric Electron Transfer Rates: From The Weak To The Strong Electronic Coupling Regions. J. Chem. Phys. 2008, 129, 184111. (33) Bonella, S.; Montemayor, D.; Coker, D. F. Linearized Path Integral Approach For Calculating Nonadiabatic Time Correlation Functions. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6715−6719. (34) Ceotto, M.; Yang, S.; Miller, W. H. Quantum Reaction Rate From Higher Derivatives Of The Thermal Flux-Flux Autocorrelation Function At Time Zero. J. Chem. Phys. 2005, 122, 044109. (35) Mandrà, S.; Valleau, S.; Ceotto, M. Deep Nuclear Resonant Tunneling Thermal Rate Constant Calculations. Int. J. Quantum Chem. 2013, 113, 1722−1734. (36) Seideman, T.; Miller, W. H. Quantum Mechanical Reaction Probabilities Via A Discrete Variable Representation-Absorbing Boundary Condition Green’s Function. J. Chem. Phys. 1992, 97, 2499−2514. (37) Manthe, U.; Miller, W. H. The Cumulative Reaction Probability As Eigenvalue Problem. J. Chem. Phys. 1993, 99, 3411−3419. (38) Guo, Y.; Thompson, D. L.; Miller, W. H. Thermal And Microcanonical Rates Of Unimolecular Reactions From An Energy Diffusion Theory Approach. J. Phys. Chem. A 1999, 103, 10308− 10311. (39) Tannor, D. J.; Garashchuk, S. Semiclassical Calculation Of Chemical Reaction Dynamics Via Wavepacket Correlation Functions. Annu. Rev. Phys. Chem. 2000, 51, 553−600. (40) Garashchuk, S.; Vazhappilly, T. Wavepacket Approach To The Cumulative Reaction Probability Within The Flux Operator Formalism. J. Chem. Phys. 2009, 131, 164108. (41) van Harrevelt, R.; Nyman, G.; Manthe, U. Accurate Quantum Calculations Of The Reaction Rates For H/ D+ CH4. J. Chem. Phys. 2007, 126, 084303. (42) Wu, T.; Werner, H.-J.; Manthe, U. First-Principles Theory For The H+ CH4→ H2+ CH3 Reaction. Science 2004, 306, 2227−2229. (43) Miller, W. H.; Zhao, Y.; Ceotto, M.; Yang, S. Quantum Instanton Approximation For Thermal Rate Constants Of Chemical Reactions. J. Chem. Phys. 2003, 119, 1329−1342. (44) Ceotto, M.; Miller, W. H. Test Of The Quantum Instanton Approximation For Thermal Rate Constants For Some Collinear Reactions. J. Chem. Phys. 2004, 120, 6356−6362. (45) Ceotto, M. Vibration-Assisted Tunneling: A Semiclassical Instanton Approach. Mol. Phys. 2012, 110, 547−559. (46) Zimmermann, T.; Vaníček, J. Path Integral Evaluation Of Equilibrium Isotope Effects. J. Chem. Phys. 2009, 131, 024111. 6465
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
Helium Isotope Enrichment by Resonant Tunneling through Nanoporous Graphene Bilayers Salvatore Mandrà,*,†,‡ Joshua Schrier,*,§ and Michele Ceotto*,⊥ †
Department of Physics, Università degli Studi di Milano, via Celoria 16, 20133 Milano, Italy Department of Chemistry and Chemical Biology, Harvard University, 12 Oxford Street, Cambridge, Massachusetts 02138, United States § Department of Chemistry, Haverford College, Haverford, Pennsylvania 19041, United States ⊥ Department of Chemistry, Università degli Studi di Milano, via Golgi 19, 20133 Milano, Italy ‡
S Supporting Information *
ABSTRACT: Graphene is impermeable to gases, but introducing subnanometer pores can allow for selective gas separation. Because graphene is only one atom thick, tunneling can play an important role, especially for low-mass gases such as helium, and this has been proposed as a means of separating 3 He from 4He. In this paper, we consider the possibility of utilizing resonant tunneling of helium isotopes through nanoporous graphene bilayers. Using a model potential fit to previously reported DFT potential energy surfaces, we calculate the thermal rate constant as a function of interlayer separation using a recently described time-independent method for arbitrary multibarrier potentials. Resonant transmission allows for the total flux rate of 3He to remain the same as the best-known single-barrier pores but doubles the selectivity with respect to 4He when the optimal interlayer spacing of 4.6 Å is used. The high flux rate and selectivity are robust against variations of the interlayer spacing and asymmetries in the potential that may occur in experiment.
1. INTRODUCTION Two stable isotopes of helium exist. The predominant isotope on earth is 4He, which is produced as uranium α-decays into lead. Very low concentrations of 3He can be found in various terrestrial sources (remaining from the accretion of the earth) and lunar regolith (deposited by solar wind), but currently the only viable source is from the decay of tritium, as a byproduct of nuclear weapons stockpiles.1 The dwindling supply of 3He, combined with the growing use of this isotope for neutrondetection sensors for nuclear security monitoring as well as basic research and hydrocarbon explorationhas raised concerns over the long-term supply of this isotope.2 A vast quantity of 3He could potentially be separated from existing helium stockpiles, natural gas, or the atmosphere, but the energy costs of performing the separation make these exceedingly impractical with current cryogenic distillation methods that require liquefying the gases.3 Membrane-based methods of gas separation are potentially more efficient because they avoid the energy cost of liquefying the gases.4 A new research direction is the use of subnanometer sized pores introduced into graphene as a high-performance separation membrane. The theoretical simulations of gas separation5,6 and experimental synthesis7,8 of nanoporous graphene materials have been the subject of several recent © 2014 American Chemical Society
reviews. Because these membranes are only one atom thick, and the pores are so small, quantum tunneling can play a significant role in the transmission of very light elements such as helium atoms, even at room temperature.9 The mass dependence of tunneling provides a potential means of separating the two isotopes under conditions where no classical separation occurs.10 Because the concentrations of 3He are so low, it is important to have a very high selectivity for the transmission of this isotope if one is to create a practical process. Certain nitrogen-functionalized pores in graphene11,12 and two-dimensional polymers13 have been predicted to have both high flux and high selectivity for helium isotope separation. The highest isotope selectivity ratio of 19 is achieved only at 10 K and rapidly approaches no selectivity by 25 K, but at the lowest temperature the flux rate is only about 10−9 mol cm−2 s−1 under a feed pressure of 1 bar. The goal of this paper is to develop nanostructures with higher selectivity and flux rate, capable of operating at higher temperatures, which will further lower the cost of obtaining 3He from natural sources. Special Issue: Franco Gianturco Festschrift Received: March 14, 2014 Revised: May 12, 2014 Published: May 22, 2014 6457
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
All of the previous studies of helium isotope separation by tunneling have considered a single barriera single nanoporous graphene or graphene-like sheet. In this paper, we investigate resonant tunneling through double-barrier nanostructures as a way to perform isotope separation. Resonant tunneling through double barriers is a textbook problem in quantum mechanics,14 and resonant tunneling of electrons is widely used for semiconductor heterostructure electronics devices.15 Unlike single-barrier tunneling, complete transmission occurs when the incoming particles have kinetic energies matching the quasi-bound states of the potential well formed by the double barrier. These resonant transmissions occur at energies much lower than the barrier height. Due to these resonant peaks the transmission can be high, even when the temperature is so low that few atoms have sufficient kinetic energy to cross the potential energy barrier classically. What we are proposing is not to be confused with the resonant tunneling of electrons within the π-bands of a graphene sheet where the potential has been modulated to form a double barrier.16 Rather, we propose transmission of a helium atom through two aligned pores of two nanoporous graphene sheets; the double barrier arises due to the nature of the potential energy surface for performing this atom transmission. In the following section we describe a model potential based on previous ab initio potential energy calculations. From a theoretical perspective, the calculation of the thermal rate constants for resonant tunneling of atoms is numerically challenging, and we describe how to perform this calculation in section 3. In section 4 we explore the optimal separation between the two barriers and the effect of van der Waal’s interactions on the isotope separation. Finally, section 5 concludes the paper.
Figure 1. Single-barrier (blue dots) and double-barrier (red line) potentials at different distances, Δ, between the peaks. As described in section 1, each potential is divided into three parts: (I) and (III) the PES is the single-barrier potential VS(x) calculated by Hauser et al.;11 (II) the PES is given by the scaled sum of two VS(x) in a way that the minimum between the two barriers maintains the same depth as for the single barrier.
2. A MODEL POTENTIAL FOR NANOPOROUS GRAPHENE MULTILAYERS The global potential energy surface (PES) of an atom approaching a two-dimensional nanoporous surface is very complex if taken in its full dimensionality. For simplicity, we consider only the minimum energy path, where the atom approaches the center of the hole in a perpendicular direction to the graphene surface, turning this into a one-dimensional problem. This is reasonable because regions of the graphene surface without a hole have a very large potential energy barrier and are thus impermeable to helium atoms.17 Both the repulsive barrier and the attractive wells are important for efficient quantum isotope separation.12 Potential wells are important for keeping the flux high by attracting the isotopes to the membrane, and the barrier height and width are important for tunneling selection (Supporting Information). The zeropoint energy and tunneling contribute in opposite directions with respect to the temperature gradient.10 At low temperatures, tunneling favors lighter isotopes whereas the zero-point energy differences of the transition state favor the heavier ones; at high temperatures, just the opposite behavior occurs. Previously reported exploratory calculations using a model PES found the optimal conditions for helium isotope separation to be a barrier height of a few hundred wavenumbers, surrounded by two minima and operating between 10 and 20 K.12 This can be achieved by creating a two-ring defect in a graphene sheet and then nitrogen doping one side of this defect;11,12 this PES is shown as the blue dots in Figure 1. We use this single-barrier PES as a model for investigating the double-barrier system. In principle, experimental approaches can control the spacing between graphene18 and
graphene oxide19 sheets (for the purposes of hydrogen storage) to yield many possible interlayer spacings. Our goal is to explore a wide range of spacings between the layers and identify the best target. Because it is too computationally demanding to obtain a high-quality ab initio PES for each of these systems, we instead use a simplified model to expose the essential aspects of this problem. The continuous red lines in Figure 1 show the double-barrier potentials offset by different barrier distances, Δ. The potential is constructed by dividing the minimum energy path into three regions, called I, II, and III in Figure 1, defined by the changing of sign of the potential in between the two barriers. In regions I and III the PES is the same as the original single-barrier potential, VS(x). This is reasonable because the transiting atom primarily interacts with the nearest membrane. The potential in region II (between the two graphene sheets) may have minima that are greater than or equal to the minima of the single-barrier case. Molecular dynamics simulations and popular empirical dispersion corrections applied to density functional theory (e.g., by Grimme20) typically assume that dispersion behaves as a pairwise sum of atom−atom interactions. However, random phase approximation21 and quantum Monte Carlo22 calculations have shown that this simple additive model is not applicable for interactions involving extended two-dimensional systems such as the graphene sheets studied here. As we will show below, the 6458
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
inconvenient for resonant tunneling calculations.35 Hence, a time-independent approach is preferred for the calculation of the thermal rate constant36−49
isotope selectivity increases with the number of quasi-bound states within this potential, which is proportional to the depth of the potential. Considering the limiting case where the potential in region II is rescaled to have the same minima depth as the single-barrier case places a conservative lower bound on the isotope selectivity. However, this assumption is not necessary for the computational method used to compute the transmission, nor for the resonant tunneling separation effect described below. Finally, for Δ larger than the support of the single-barrier potential VS(x), the double-barrier potential is well approximated by the following equation VD(x ;Δ) ≈ VS(x − Δ/2) + VS(x + Δ/2)
k(T ) =
+∞
1 2π ℏQ R (T )
∫−∞
e − βE T ( E ) d E (2)
where QR(T) is the reactants partition function and T (E ) =
1 (2π ℏ)2 Tr[F1̂δ(E − Ĥ ) F2̂ δ(E − Ĥ )] 2
(3)
is the cumulative reaction probability evaluated between the dividing surface f1(s) = 0 and f 2(s) = 0 located at s = 0. F̂ is the flux operator and Ĥ is the Hamiltonian operator.23 To be consistent with commonly used notations and to avoid confusion, T (with no argument specified) refers to temperature, and T(E) refers to the transmission probability as a function of kinetic energy. We will always include the argument in the transmission probability. In a recent paper,35 we presented a time-independent method for calculating the thermal rate constant, k(T), for any arbitrary multibarrier potential. The method successfully treats double-barrier potentials with many quasi-bound states. It is based on a fast and robust procedure for calculating k(T) as the thermal average of the transmission probability T(E) according to eq 2. It solves the Schrödinger equation as an ordinary differential equation (ODE) with the energy as a parameter, thus avoiding the related eigenvalue problem. The main idea is to force the asymptotic eigenfunctions ψp(x) of equation
|x | ≫ Δ (1)
Small variations in the PES modify the barrier height and break the symmetry of the double barrier and can potentially ruin the resonant tunneling effect. These variations might arise from defects in the pores, thermal fluctuations of pore atoms, or angled trajectories through the double-barrier structure. To study this possibility, we also considered the case where one barrier is 20% higher than the other, as shown in Figure 2. We discuss the results in the following section.
3. A TIME-INDEPENDENT METHOD FOR RESONANT TUNNELING Time-dependent theories typically calculate thermal rates by integrating a flux autocorrelation function,23−34 but this is
⎡ ℏ2 d2 ⎤ Ĥ ψp(x) = ⎢ − + V (x)⎥ψp(x) = Eψp(x) 2 ⎣ 2m dx ⎦
(4)
at fixed eigenvalue E = p /2m to satisfy the planewave boundary conditions at distances, L, that are sufficiently far away from the scattering potential, 2
⎧ e ipx / ℏ x ≪ −L ⎪ ψp(x) = ⎨ t ipx / ℏ r −ipx / ℏ + e x ≫ +L ⎪ 2e ⎩ |t | |t |
(5)
These are obtained by taking the boundary conditions lim ψp(x) = 1
(6a)
lim ∂xψp(x) = ip /ℏ
(6b)
x →−L
x →−L
where x is in the asymptotic region. Finally, one obtains expressions for the transmission and reflection probabilities, T (E) = |t |2 = lim 4 ψp(x) + x→L
ℏ ∂xψ (x) ip p
−2
(7a)
2 ⎞−1 ⎛ 1 ℏ ⎜ R(E) = |r | = lim 1 − ⎜1 + ψp(x) − ∂xψp(x) ⎟⎟ x→L 4 ip ⎝ ⎠ 2
(7b)
where x is an asymptotic position after the barrier passage. The advantage of this method over the time-dependent approach is that it allows for a fine energy grid when the transmission probability is calculated, which is crucial for accurately integrating eq 2 when there are very narrow resonant
Figure 2. Potential energy surfaces for a model double nanoporous graphene pair of sheets. The potential is not symmetric for a more realistic modeling. Notation as in Figure 1. 6459
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
Figure 3. 3He and 4He transmission probabilities T(E) for the potential profile of Figure 1. Left panel: lin−lin scale. Right panel: log−log scale. All of the distances considered in this paper are between the two extreme distances Δ = 8.0 au and Δ = 12.0 au, shown here.
significantly the thermal rate at very low temperatures.35 The low energy portion of the log−log plot (E < 0.0001 eV) also shows that T(E) cannot be safely approximated as zero at these low energies and is different for the two isotopes. Although this is typically ignored at high temperatures, it is too important at low temperatures. For these reasons, the lower energy bound has been fixed at Emin = 10−15 au and the lower bound of the transmission probability fixed at T(E)min = 10−40 au. A further inspection of the top-left panel of Figure 3 shows that there are eight resonant peaks for the 4He transmission probability and only seven for the 3He. This difference is more evident as the graphene sheets are separated, as shown in the bottom left panel. Each of the resonant transmissions corresponds to a quasi-bound (metastable) state of the atom inside the double-barrier potential. Approximating this region as a harmonic potential, then the eigenvalues spacing is inversely proportional to the square root of the oscillator (i.e., the helium atom) mass. As a result, the 4He peaks occur at lower energies. Because the potential has a finite height, there are only a limited number of these quasi-bound states, and because the 4He energies are lower, there are more quasi-bound states for this isotope. However, in the Δ = 8.0 au case, the 3He resonant peak is the lower one because there are no quasibound states with energies less than the original asymptotic energy (explained below). Boltzmann-weighted energy integration of T(E) as in eq 2 yields the thermal rate constants, k(T). In Figure 4 we show the product of the thermal rate constant with the partition function, i.e., k(T) QR(T), to more clearly show the tunneling effects. The direct comparison with the single barrier can still be done because, for one-dimensional systems, partition functions are independent of the asymptotic shape of PES. In the same figure, the single barrier (violet left-triangles; the topmost line in both graphs) and the transition state theory (TST) (dashed black line; the bottom-most line in both graphs) rates are reported. On the upper panel the rates are for 3He passage and on the lower for 4He. Even though the thermal rate for the single barrier is typically greater than that of the double-barrier,
transmission peaks. In this case, the method allows for using a coarse grid in flat regions and a fine one in the regions around resonant peaks. This improves the numerically convergence for the energy integration in eq 2. Details about the numerical convergence are reported in the Supporting Information. A more formal presentation can be found in ref 35.
4. RESULTS AND DISCUSSION In this section, results on the transmission probability T(E) and the thermal rate constant k(T) are presented together with He isotope selection considerations. The numerical convergence of these calculations is challenging: Extra precision is required to calculate the rate at the low temperatures (5−30 K) needed for isotope separation. Integrating the contribution of the resonant peaks requires calculating T(E) for very low energies. Thermal Rate Constants. Figure 3 shows the transmission probabilities for two resonant tunneling cases, demonstrating the presence of several resonance peaks at kinetic energies below 0.022 eV. To demonstrate that our computational approach can treat any arbitrary barrier distance, here we show the results for the most closely (Δ = 8.0 au, top panels) and distantly (Δ = 12.0 au, bottom panels) spaced graphene sheets that we considered; all other calculations in this paper are intermediate between these two extremes. In the left panels of Figure 3 both axes are in linear−linear scale, and the right panels are for the same data but in log−log scale. The linear scale representation shows that each resonant peak is very narrow and has a maximum transmission of 1, in agreement with the symmetric double-barrier potential resonant tunneling rule. The narrowness of the peaks presents a numerical challenge for the energy integration of eq 2. For this purpose, we employed both a coarse-grained grid, as explained above, and a resonant peak transmission splined profile as suggested in refs 35 and 50. Thus, the energy grid spacing is of the order of 10−15 au. The log−log plots on the right panels of Figure 3 show that the chosen integration energy range includes all resonant peaks. This is important because exclusion of even one of the low-energy resonant transmission peaks can change 6460
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
Figure 4. Thermal rate constant tunneling effects. The product of k(T) QR(T) is plotted for 3He (upper panel) and 4He (lower panel). Transition state theory (TST) results in the dashed black line and the single barrier in the purple line (with the largest value in both insets). Other lines are for different barrier distances, Δ.
Figure 5. Thermal rate constant ratio, k3/k4, where k3 and k4 are the thermal rate constants for 3He and 4He, respectively. The upper panel figure shows the temperature (T) dependence of the thermal rate constant ratio; the horizontal black dashed line shows the isotopeindependent classical transition state theory (TST) value. The bottom panel shows the interbarrier distance (Δ) dependence of the thermal rate ratio; the horizontal black dashed lines show the single-barrier result at each temperatures.
order-of-magnitude discrepancies occur only at very low temperatures (T < 10 K); at higher temperatures (T > 10 K) the thermal rate constants of single and double barriers are comparable. The deviation of the thermal rate constants from the TST result show that transmission through both the singleand double-barrier graphene nanopores is greatly increased by quantum tunneling effects. This figure also shows nonlinear behavior of the thermal rate for the double-barrier potential as a function of the interbarrier distance, Δ. For 3He, shown in the top panel, the order of k(T) QR(T) as a function of Δ is 10 au > 8 au ≈ 11 au > 12 au > 9 au. For 4He, shown in the bottom panel, the order of k(T) QR(T) as a function of Δ is 9 au > 11 au > 10 au ≈ 12 au > 8 au. As we will discuss below, this is a direct consequence of the variation of the number and shape of the resonant tunneling peaks with the distance Δ and the masses of the particles. To see if the double porous graphene sheets can better select 3 He from 4He than the single sheet, we show the ratio of the thermal rate constants for the two isotopes, k3/k4, in the upper panel of Figure 5. In a classical approximation no quantum tunneling occurs, and if one also neglects the different (massdependent) zero-point energy contributions to the potential energy barrier, then the TST ratio is temperature-independent and given by the ratio of the partition functions. In contrast, for the deep-resonant tunneling regime k3/k4 varies significantly at small temperatures and gradually converges to the single-barrier ratio at large temperature (T > 15 K). This is due to the change in the number of resonant tunneling peaks for 3He versus 4He. Even if the rate ratio is the same for the single and the double barrier at larger Δ distance, the rate value is different. Interestingly, for sufficiently low temperature, the ratio for the double-barrier potential can be smaller, comparable to, or even larger than the ratio for the single barrier. This effect is
exponentially enhanced when the temperature is further lowered. The double-barrier resonant tunneling effect is so pronounced that for temperatures less than T < 15 K, an inverse isotope effect is observed for Δ = 9.0 au and Δ = 9.4 au To better appreciate the He isotope enrichment effects by resonant tunneling, the lower panel of Figure 5 shows the ratio k3/k4 for the double-barrier potential as a function of the interbarrier distance Δ. For reference, the (Δ-independent) single-barrier k3/k4 ratio is shown as horizontal dashed lines for each temperature. For temperatures ≥15 K, no significant variations in the isotope enrichment between the single- and double-barrier cases is observed. At lower temperatures, the double-barrier system offers the possibility of improved 3He selectivity only for certain values of Δ. For example, at 10 K, the double-barrier potential can select the lighter isotope with more than twice the selectivity of single-barrier potential when Δ = 8.72 au. However, a small increase in Δ can reduce the isotope selectivity below that of the single barrier. A key result of this paper is that for temperatures below 12 K, an appropriate choice of the distance Δ between the double porous graphene sheets can substantially improve selection compared to the single graphene layer. Isotope Selection versus Barriers Distances and Potential Shape. To better understand the mechanism of enhancing the isotope selection induced by the He resonant tunneling, we count the number of resonant transmission peaks in T(E), as a function of Δ. As shown in Figure 6, the number of resonant peaks for 4He is always greater-than-or-equal-to that of 3He, but the number of resonant peaks is not a monotonically increasing step function. Even though the 6461
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
potential value (see PES in Figure 1), energies of lower metastable states can also be negative. However, only the metastable states with positive energy contribute to the transmission probability T(E), because the energies of incoming helium atoms are always positive. Therefore, even if the total number of metastable states with both positive and negative energies grows by increasing the distance between the graphene sheets, states with negative energy are trapped and do not contribute to the transmission. The resonant peaks shifts to lower energies as Δ is increased, until one gets to a case with the closed channels, shown as the orange lines in Figure 7. The sudden closure of a channel is responsible for the number of peaks decreasing observed in Figure 6 for both isotopes. In conclusion, there are two opposite trends responsible for the number of resonant peaks. The positive trend is due to the increasing of anharmonicity as Δ is increased. The negative trend is due to the metastable state energy switching from positive to negative values. We stress that similar results occur if the van der Waals well between the two graphene sheets (region II) is deeper. At low temperature, the isotope separation is mainly driven by the difference of the number of resonant peaks between 3He and 4 He (as shown in Figures 5 and 6), but a deeper central well primarily increases the number of negative energy resonant peaks, which do not contribute to the transmission. Because 4 He is heavier, more negative energy peaks for 4He are present for a deeper central well which will increase the 3He/4He ratio. Asymmetric Potential Model. In any real experiment, imperfect alignment of the pores, variations in the pore geometries, and thermal fluctuations will all distort the ideal PES. When the double barrier is no longer symmetrical, the resonant transmission peaks no longer have unit transmission probabilities. To study whether these distorted potentials are still capable of isotope selection and of high flux, we examined a set of model asymmetric PES shown in Figure 2, where one barrier is scaled to lower values. The transmission probability for both isotopes is shown in Figure 8 for the smaller (upper panel) and the larger (lower panel) interbarrier distances considered. As expected, the T(E) peaks are no longer equal to unity because of the asymmetry in the potential profile. The position of the T(E) peaks changes more with increasing Δ for the heavier isotope, similar to what was observed for the symmetric potential. Does asymmetry change the thermal rate values, k(T) QR(T)? Yes. As shown in Figure 9, all of the asymmetric PESs actually have a higher rate than the corresponding symmetric PESs for both isotopes, for any interbarrier distance Δ. This is of practical importance, because the higher rates mean a faster flux of atoms across the membrane. Does asymmetry change the isotope selectivity? No. As shown in Figure 10, the thermal rate has a ratio k3/k4 similar to the symmetric one for different temperatures and distances Δ. (The ratios for the symmetric potential, previously shown in Figure 5, are indicated here as open symbols, so as to enable a direct comparison against the asymmetrical PES results.) The isotope selectivity is preservedand in some cases slightly enhancedfor the asymmetric potential. This is very encouraging for future experiments, as it means that the resonant tunneling effect is robust against these types of small fluctuations, and suggests that trajectories deviating by a small angle or displacement from the one studied here will have the same qualitative behavior. The behavior of the resonant enhancements as a function of Δ is preserved because the
Figure 6. Comparison between the number of resonant peaks for 4He and 3He. Vertical lines are in correspondence with significant k3/k4 ratio changes in Figure 5. As one can see, variations of k3/k4 and of the number of peaks are only partially correlated.
difference of the number of resonant peaks and the discontinuities of the ratio k3/k4 have a common trend, there is no direct correspondence between them. This is because k(T) is not only the sum of the resonant peaks but it also depends on the shape of the transmission probability T(E). To better understand the physical reasons why the number of peaks decrease for some Δ values, we focus on the log−log plot of T(E) for the lowest energy scattering region. The upper panel of Figure 7 reports T(E) for different values of Δ for 3He
Figure 7. Transmission probability T(E) for 3He (top panel) and 4He (bottom panel) at low scattering energies. Different curves correspond to different barrier distances Δ. At large enough distance, the lower resonant peak is suppressed because the corresponding metastable state energy is below the asymptotic one.
and the lower panel shows the 4He results. (The specific set of Δ values are different in the two insets, to better illustrate the varying number and shape of the peaks for the two isotopes.) As Δ is increased, peaks shift to lower energy due to the increment of anharmonicity caused by the van der Waals interaction between He and the graphene sheets. Because the van der Waals well minimum is below the zero asymptotic 6462
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
Figure 10. Thermal rate constant ratios and isotope selection as a function of Δ. Open symbols and lines show the symmetric and asymmetric PES results, respectively. Lower temperatures and the asymmetric potential of Figure 2 offer better isotope selection.
metastable state energies located at the bottom of the intrabarrier well are not modified. Because these are the states mainly responsible for the thermal rate values at low temperatures, the optimal value of Δ for resonant tunneling are nearly the same for both the symmetric and asymmetric potentials. Feasibility. Can these double-barrier nanostructures purify a macroscopic amount of 3He? We can obtain an upper-bound estimate of the helium flux by making two assumptions. First, we assume that any incident particle on the two-dimensional surface is transmitted through the potential associated with a pore. In reality, the pores will be spread over a larger distance, with a density determined by the synthesis procedure. Second, we assume that all trajectories are normal to the graphene membranes. Angled trajectories are neglected as they are typically reflected by the potential or, if the angle is very small, will introduce a small asymmetry in the potential that does not qualitatively change the results. Both of these assumptions increase the predicted flux rate by a multiplicative factor. As shown in Figure 5, the isotope separation for a double graphene layer is maximized when the distance between the two graphene sheets is Δ = 8.72 au. Assuming a pressure of 1 atm, the daily amounts of 3He atoms that cross a square-meter of a single-layer or a double-layer membrane at T = 10 K are almost identical and equal to ≈800 g per day. However, as shown in Figure 5, the output from the double graphene layer will be twice as pure as the single layer; i.e., the permeate will contain half as much 4He. Due to resonant tunneling a double graphene layer is more efficient than a single graphene layer for this separation. Can this double-barrier structure be made? Yes. In principle, experimental approaches can control the spacing between graphene18 and graphene oxide19 sheets to yield the optimal interlayer spacing of 8.7 au (4.6 Å) predicted above. One possible strategy is to intercalate molecules between the graphene layers. Our optimal interlayer spacing is larger than the 3.8 Å spacing achieved by intercalating methane molecules,51 but only slightly smaller than the 4.7 Å spacing achieved by intercalating tetrabutylammonium cations,52 which suggests that intermediate values can be achieved by using different types of molecules. Another possibility is to intercalate metal atoms between the graphene layers. Our optimal interlayer spacing is larger than the 3.7 Å interlayer spacing for lithium-intercalated graphene,53 but smaller than the 6.5 Å spacing for Au6 intercalated graphene,54 again suggesting that intermediate values can be achieved by intercalating different
Figure 8. Transmission probabilities for the asymmetric potential setup of Figure 2. The scale is linear for both axes. Upper panel for a distance of Δ = 8 au and lower panel for Δ = 12.0 au. The resonant peaks transmission probability is not unity because of the asymmetry in the potential.
Figure 9. k(T) QR(T) values comparison between the symmetric (SYM) and the asymmetric (ASYM) PES of Figure 4. Upper panel for 3 He and lower one for 4He. Dashed line for transition state theory (TST) values. Asymmetric potential rates are always bigger than symmetric ones at same barrier distances Δ.
asymmetry is sufficiently small that it does not affect the resonant peaks in the deep tunneling regime. Consequently, the 6463
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
■
types of metals. The pores could be simultaneously created and decorated with nitrogen by means of nitrogen ion implantation.55,56 Our analysis of the asymmetrical potential shows that typical variations will not be detrimental to performance.
ASSOCIATED CONTENT
S Supporting Information *
Numerical convergence analysis for the tunneling calculations and derivation of the separation rate. This material is available free of charge via the Internet at http://pubs.acs.org.
■
REFERENCES
(1) Wittenberg, L.; Cameron, E.; Kulcinski, G.; Ott, S.; Santarius, J.; Sviatoslavsky, G.; Sviatoslavsky, I.; Thompson, H. A Review Of 3He Resources And Acquisition For Use As Fusion Fuel. Fusion Technology 1992, 21, 2230−2253. (2) Managing Critical Isotopes: Weaknesses In DOE’s Management Of Helium-3 Delayed The Federal Response To A Critical Supply Shortage. GAO-11-472, U.S. Government Accountability Office, 2011. (3) Shea, D. A.; Morgan, D. Helium-3 Shortage: Supply, Demand, And Options For Congress; Congressional Research Service: Washington, DC, 2010 (4) Bernardo, P.; Drioli, E.; Golemme, G. Membrane Gas Separation: A Review/State Of The Art. Ind. Eng. Chem. Res. 2009, 48, 4638− 4663. (5) Jiao, Y.; Du, A.; Hankel, M.; Smith, S. C. Modelling Carbon Membranes For Gas And Isotope Separation. Phys. Chem. Chem. Phys. 2013, 15, 4832−4843. (6) Jiao, Y.; Hankel, M.; Du, A.; Smith, S. C. Porous Graphene And Nanomeshes. Graphene Chemistry: Theoretical Perspectives 2013, 129− 151. (7) Russo, P.; Hu, A.; Compagnini, G. Synthesis, Properties And Potential Applications Of Porous Graphene: A Review. Nano-Micro Lett. 2013, 5, 260−273. (8) Jiang, L.; Fan, Z. Design Of Advanced Porous Graphene Materials: From Graphene Nanomesh To 3D Architectures. Nanoscale 2014, 6, 1922−1945. (9) Schrier, J. Helium Separation Using Porous Graphene Membranes. J. Phys. Chem. Lett. 2010, 1, 2284−2287. (10) Schrier, J.; McClain, J. Thermally-Driven Isotope Separation Across Nanoporous Graphene. Chem. Phys. Lett. 2012, 521, 118−124. (11) Hauser, A. W.; Schwerdtfeger, P. Nanoporous Graphene Membranes For Efficient 3He/4He Separation. J. Phys. Chem. Lett. 2012, 3, 209−213. (12) Hauser, A. W.; Schrier, J.; Schwerdtfeger, P. Helium Tunneling Through Nitrogen-Functionalized Graphene Pores: Pressure-And Temperature-Driven Approaches To Isotope Separation. J. Phys. Chem. C 2012, 116, 10819−10827. (13) Brockway, A. M.; Schrier, J. Noble Gas Separation Using PG-ES X (X= 1, 2, 3) Nanoporous Two-Dimensional Polymers. J. Phys. Chem. C 2012, 117, 393−402. (14) Bohm, D. Quantum Theory; Prentice-Hall Physics Series; Prentice Hall: Englewood Cliffs, NJ, 1951 (15) Tsu, R.; Esaki, L. Tunneling In A Finite Superlattice. Appl. Phys. Lett. 2003, 22, 562−564. (16) Pereira, J. M.; Vasilopoulos, P.; Peeters, F. Graphene-Based Resonant-Tunneling Structures. Appl. Phys. Lett. 2007, 90, 132122− 132122. (17) Leenaerts, O.; Partoens, B.; Peeters, F. Graphene: A Perfect Nanoballoon. Appl. Phys. Lett. 2008, 93, 193107. (18) Jin, Z.; Lu, W.; O’Neill, K. J.; Parilla, P. A.; Simpson, L. J.; Kittrell, C.; Tour, J. M. Nano-Engineered Spacing In Graphene Sheets For Hydrogen Storage. Chem. Mater. 2011, 23, 923−925. (19) Kim, B. H.; Hong, W. G.; Moon, H. R.; Lee, S. M.; Kim, J. M.; Kang, S.; Jun, Y.; Kim, H. J. Investigation On The Existence Of Optimum Interlayer Distance For H2 Uptake Using Pillared-Graphene Oxide. Int. J. Hydrogen Energy 2012, 37, 14217−14222. (20) Grimme, S. Semiempirical GGA-Type Density Functional Constructed With A Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (21) Dobson, J. F.; White, A.; Rubio, A. Asymptotics Of The Dispersion Interaction: Analytic Benchmarks For Van Der Waals Energy Functionals. Phys. Rev. Lett. 2006, 96, 073201. (22) Drummond, N.; Needs, R. Van Der Waals Interactions Between Thin Metallic Wires And Layers. Phys. Rev. Lett. 2007, 99, 166401. (23) Miller, W. H.; Schwartz, S. D.; Tromp, J. W. Quantum Mechanical Rate Constants For Bimolecular Reactions. J. Chem. Phys. 1983, 79, 4889−4898. (24) Sklarz, T.; Kay, K. Semiclassical Initial Value Treatment Of Correlation Functions. J. Chem. Phys. 2004, 120, 2606−2617.
5. CONCLUSIONS We have applied a recently developed time-independent method for thermal rate calculations of multibarrier potentials to demonstrate the feasibility of isotope separation using resonant tunneling through bilayer graphene nanopores. From a fundamental perspective, this extends a basic quantum mechanical phenomenon to heavy atoms. From the perspective of chemical engineering, this is the first application of resonant tunneling to separations and it makes a significant improvement to the helium isotope separation. Considering that higher flux means lower selectivity, and vice versa,57 there is a trade-off between these two trends. However, at a working temperature of 10 K we can achieve the same 3He flux as the best single barriers but double the selectivity with respect to 4He. The high flux rate and isotope selectivity are robust against asymmetries in the potential. Moreover, the optimal separation between the two graphene sheets is intermediate to known results for molecular and metal intercalation compounds of graphene. We hope that the significant improvement in isotope separation will motivate future computational and experimental work on bilayer graphene nanopore systems. Future investigations will consider ab initio calculations of the two layer graphene sheets potential in more than one dimension, multidimensional tunneling, and the effects of different size pores.
■
Article
AUTHOR INFORMATION
Corresponding Authors
*S. Mandrà: e-mail, [email protected]. *J. Schrier: e-mail, [email protected]. *M. Ceotto: e-mail, [email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We thank Geoffrey Martin-Noble, Anahita Nourmahnad, and Kylen Solvik for careful readings of the manuscript. S.M. and M.C acknowledge the CINECA and the Regione Lombardia award under the LISA initiative, for the availability of high performance computing resources and support. S.M. also acknowledges the Air Force Office of Scientific Research (prime sponsor) and the University of California, San Diego, for the partial support to this work under the grants FA955012-1-0046 and 10323836-SUB. J.S. was supported by the Research Corporation for Scientific Advancement’s Cottrell Scholar grant and used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. 6464
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
The Journal of Physical Chemistry A
Article
(47) Zimmermann, T.; Vaníček, J. Three Applications Of Path Integrals: Equilibrium And Kinetic Isotope Effects, And The Temperature Dependence Of The Rate Constant Of The [1, 5] Sigmatropic Hydrogen Shift In (Z)-1, 3-Pentadiene. J. Mol. Model. 2010, 16, 1779−1787. (48) Buchowiecki, M.; Vaníček, J. Direct Evaluation Of The Temperature Dependence Of The Rate Constant Based On The Quantum Instanton Approximation. J. Chem. Phys. 2010, 132, 194106. (49) Vaníček, J. Beyond Transition State Theory: Accurate Description Of Nuclear Quantum Effects On The Rate And Equilibrium Constants Of Chemical Reactions Using Feynman Path Integrals. CHIMIA Int. J. Chem. 2011, 65, 715−719. (50) Ryaboy, V.; Lefebvre, R. Flux-Flux Correlation Function Study Of Resonance Effects In Reactive Collision. J. Chem. Phys. 1993, 99, 9547−9552. (51) Huang, Q.; Chen, X.; Lin, J.; Li, K.; Jia, Y.; Liu, J.; Guo, L.; Wang, W.; Wang, G. Preparation Of Quasi-Free-Standing Graphene With A Super Large Interlayer Distance By Methane Intercalation. J. Phys. Chem. C 2011, 115, 20538−20545. (52) Sirisaksoontorn, W.; Adenuga, A. A.; Remcho, V. T.; Lerner, M. M. Preparation And Characterization Of A Tetrabutylammonium Graphite Intercalation Compound. J. Am. Chem. Soc. 2011, 133, 12436−12438. (53) Qi, Y.; Guo, H.; Hector, L. G.; Timmons, A. Threefold Increase In The Young’s Modulus Of Graphite Negative Electrode During Lithium Intercalation. J. Electrochem. Soc. 2010, 157, A558−A566. (54) Sapkota, I.; Roundtree, M. A.; Hall, J. H.; Wang, X.-Q. Tunable Band Gap In Gold Intercalated Graphene. Phys. Chem. Chem. Phys. 2012, 14, 15991−15994. (55) Guo, B.; Liu, Q.; Chen, E.; Zhu, H.; Fang, L.; Gong, J. R. Controllable N-Doping Of Graphene. Nano Lett. 2010, 10, 4975− 4980. (56) Åhlgren, E.; Kotakoski, J.; Krasheninnikov, A. Atomistic Simulations Of The Implantation Of Low-Energy Boron And Nitrogen Ions Into Graphene. Phys. Rev. B 2011, 83, 115424. (57) Robeson, L. M. The Upper Bound Revisited. J. Membr. Sci. 2008, 320, 390−400.
(25) Pollak, E.; Martin-Fierro, E. New Coherent State Representation For The Imaginary Time Propagator With Applications To ForwardBackward Semiclassical Initial Value Representations Of Correlation Functions. J. Chem. Phys. 2007, 126, 164107. (26) Martin-Fierro, E.; Pollak, E. Forward-Backward Semiclassical Initial Value Series Representation Of Quantum Correlation Functions. J. Chem. Phys. 2006, 125, 164104. (27) Liu, J.; Miller, W. H. Real Time Correlation Function In A Single Phase Space Integral Beyond The Linearized Semiclassical Initial Value Representation. J. Chem. Phys. 2007, 126, 234110. (28) Yamamoto, T.; Wang, H.; Miller, W. H. Combining Semiclassical Time Evolution And Quantum Boltzmann Operator To Evaluate Reactive Flux Correlation Function For Thermal Rate Constants Of Complex Systems. J. Chem. Phys. 2002, 116, 7335−7349. (29) Tromp, J. W.; Miller, W. H. The Reactive Flux Correlation Function For Collinear Reactions H+ H2, Cl+ HCl And F+ H2. Faraday Discuss. Chem. Soc. 1987, 84, 441−453. (30) Bonella, S.; Monteferrante, M.; Pierleoni, C.; Ciccotti, G. Path Integral Based Calculations Of Symmetrized Time Correlation Functions. I. J. Chem. Phys. 2010, 133, 164104. (31) Shi, Q.; Geva, E. Nonradiative Electronic Relaxation Rate Constants From Approximations Based On Linearizing The PathIntegral Forward-Backward Action. J. Phys. Chem. A 2004, 108, 6109− 6116. (32) Zhu, W.; Zhao, Y. Effects Of Anharmonicity On DiffusiveControlled Symmetric Electron Transfer Rates: From The Weak To The Strong Electronic Coupling Regions. J. Chem. Phys. 2008, 129, 184111. (33) Bonella, S.; Montemayor, D.; Coker, D. F. Linearized Path Integral Approach For Calculating Nonadiabatic Time Correlation Functions. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 6715−6719. (34) Ceotto, M.; Yang, S.; Miller, W. H. Quantum Reaction Rate From Higher Derivatives Of The Thermal Flux-Flux Autocorrelation Function At Time Zero. J. Chem. Phys. 2005, 122, 044109. (35) Mandrà, S.; Valleau, S.; Ceotto, M. Deep Nuclear Resonant Tunneling Thermal Rate Constant Calculations. Int. J. Quantum Chem. 2013, 113, 1722−1734. (36) Seideman, T.; Miller, W. H. Quantum Mechanical Reaction Probabilities Via A Discrete Variable Representation-Absorbing Boundary Condition Green’s Function. J. Chem. Phys. 1992, 97, 2499−2514. (37) Manthe, U.; Miller, W. H. The Cumulative Reaction Probability As Eigenvalue Problem. J. Chem. Phys. 1993, 99, 3411−3419. (38) Guo, Y.; Thompson, D. L.; Miller, W. H. Thermal And Microcanonical Rates Of Unimolecular Reactions From An Energy Diffusion Theory Approach. J. Phys. Chem. A 1999, 103, 10308− 10311. (39) Tannor, D. J.; Garashchuk, S. Semiclassical Calculation Of Chemical Reaction Dynamics Via Wavepacket Correlation Functions. Annu. Rev. Phys. Chem. 2000, 51, 553−600. (40) Garashchuk, S.; Vazhappilly, T. Wavepacket Approach To The Cumulative Reaction Probability Within The Flux Operator Formalism. J. Chem. Phys. 2009, 131, 164108. (41) van Harrevelt, R.; Nyman, G.; Manthe, U. Accurate Quantum Calculations Of The Reaction Rates For H/ D+ CH4. J. Chem. Phys. 2007, 126, 084303. (42) Wu, T.; Werner, H.-J.; Manthe, U. First-Principles Theory For The H+ CH4→ H2+ CH3 Reaction. Science 2004, 306, 2227−2229. (43) Miller, W. H.; Zhao, Y.; Ceotto, M.; Yang, S. Quantum Instanton Approximation For Thermal Rate Constants Of Chemical Reactions. J. Chem. Phys. 2003, 119, 1329−1342. (44) Ceotto, M.; Miller, W. H. Test Of The Quantum Instanton Approximation For Thermal Rate Constants For Some Collinear Reactions. J. Chem. Phys. 2004, 120, 6356−6362. (45) Ceotto, M. Vibration-Assisted Tunneling: A Semiclassical Instanton Approach. Mol. Phys. 2012, 110, 547−559. (46) Zimmermann, T.; Vaníček, J. Path Integral Evaluation Of Equilibrium Isotope Effects. J. Chem. Phys. 2009, 131, 024111. 6465
dx.doi.org/10.1021/jp502548r | J. Phys. Chem. A 2014, 118, 6457−6465
