Revealing the Innate Subnanometer Porous Structure of Carbon Nanomembranes with Molecular Dynamics Simulations and Highly-Charged Ion Spectroscopy

TL;DR

This paper uses simulations and a new ion spectroscopy method to reveal the tiny porous structure of carbon nanomembranes.

Contribution

The study introduces a novel method combining molecular dynamics and highly-charged ion spectroscopy to characterize CNM structure.

Findings

CNMs likely contain under-coordinated carbon with subnanometer pores.

The porous structure is reactive in air but stabilized by hydrogen and oxygen groups.

Simulated ion charge exchange data matches experimental results.

Abstract

Carbon nanomembranes (CNMs) are nanometer-thin disordered
carbon
materials that are suitable for a range of applications, from energy
generation and storage through to water filtration. The structure–property
relationships of these nanomembranes are challenging to study using
traditional experimental characterization techniques, primarily due
to the radiation sensitivity of the free-standing membrane. Highly
charged ion spectroscopy is a novel characterization method that is
able to infer structural details of the carbon nanomembrane without
concern about induced damage affecting the measurements. Here we employ
molecular dynamics simulations to produce candidate structural models
of terphenylthiol-based CNMs with varying degrees of nanoscale porosity
and compare predicted ion charge exchange data and tensile moduli
to experiment. The results suggest that the in-vacuum CNM composition
likely comprises a significant fraction of under-coordinated carbon,
with an open subnanometer porous structure. Such a carbon network
would be reactive in the atmosphere and would be presumably stabilized
by hydrogen and oxygen groups under atmospheric conditions.

Figures

Exit charge state spectra of 90 keV Xe10+ after transmission through the CNM before (a) and after (b) contaminant removal through ion irradiation. (c) A fluence of ∼3 × 1011 ion impacts per cm2 is necessary to reach an equilibrium mean charge state of q out ∼ 2.

Schematic diagram of the specular reflection imposed on atoms that enter the exclusion cylinders during region-restricted dynamics.

Planar top-down view of a selection of exclusion cylinder-enforced carbon nanomembrane structures with atoms colored according to their coordination: sp = blue, sp 2 = orange, sp 3 = purple, and single-coordinated atoms are shown in green.

Predicted average tensile moduli of the CNM models as a function of the number of exclusion cylinders (panel a), where the color of each point represents the annealing time to which each structure has been subjected. Each data point has been averaged over both in-plane directions, and error bars denote one standard deviation from the mean. The shaded green region indicates the reported experimental tensile modulus range. Panel b presents the same tensile data as a function of sp 2/sp for 1, 16, and 64 ps simulated annealed structures, where the enforced cylinder counts are indicated by the different markers and colors.

Simulated exit charge state for Xe ions incident on various hole-enforced nanomembrane structures. Structures that have been annealed for 0, 16, and 64 ps are shown in dark blue shades, and the 9 ps structure is shown in orange. It should be noted that since there is only one unique CNM structure for each configuration, there is no statistical uncertainty evaluation in the calculated exit charge state values.

Experimental (left half) and simulated (right half) charge exchange spectra through carbon nanomembranes for 72 keV Xe8+, 135 keV Xe15+, and 180 keV Xe20+ from left to right. For the simulations, we used a 150 exclusion cylinder membrane annealed for 9 ps shown far right in a top-down view, with atoms colored according to coordination, as per Figure . Here, the left half of the image is rendered using a bond-based visualization and the right half uses a van der Waals representation with an atom radius of 1.7 Å. Note that the scattering angles are different in experiment and simulation, details are given in the text.

Pore area distributions shown in blue-shaded regions (left y-axis) for event-varied (a), force-varied (b), and the 4 ps annealed region-restricted simulations (panel c), with the total pore count shown in solid red triangles (right y-axis). Mean pore area is indicated by horizontal blue lines for each distribution. Note that panels a) and b) have different y-axis scales compared to panel c).

Full text

Introduction

Nanometer-thin carbon nanomembranes (CNMs) represent an intriguing class of thin-film materials with great potential in various nanotechnologies such as fuel cells,? energy storage, ?,? photoactive electronics, ?,? nanofiltration, and water desalination. ?−? ? ? ? They are formed by irradiating self-assembled monolayers (SAMs) of organic molecules (e.g., biphenylthiol or terphenylthiol) grown on gold substrates with low-energy electrons, resulting in a cross-linked soft film approximately one nanometer thick. ?−? ? ? ? ? ? ? ? CNMs are generally considered to be noncrystalline soft carbon films,? with high selectivity for gas and water permeation. ?,? Unlike graphite and related fullerenes, CNMs exhibit a lower tensile modulus between 6–12 GPa and ultimate tensile strengths between 400–700 MPa, depending on the precursor SAM. ?,?

In all applications, the structural characterization of the CNMs is crucial, as to a large extent the structure defines their functional properties. However, the characterization of freestanding CNMs is significantly more challenging in comparison to typical 2D materials, such as graphene and MoS_2_. ?,? For example, imaging CNMs via high-resolution transmission electron microscopy (HRTEM) induces graphitization due to the high electron beam energies, and hence hinders structural investigation.? Scanning tunneling or atomic force microscopy (STM, AFM) can provide structural information only for membranes on solid substrates. Only recently, successful studies of CNMs at low temperatures were reported using noncontact AFM.? As such, there is a pressing need for characterization methods that can suitably probe the disordered structure of freestanding CNMs.

An alternative to these traditional experimental methods are keV beams of highly-charged ions (HCIs), which have recently been introduced as a characterization technique to provide structural information on CNMs.? There are two key aspects to this method. Even though ions may modify the local site after impact, each ion carries with it information regarding the pristine sample state that it encountered, which is then interrogated after transmission at a detector. The second is that total ion fluences used are extremely low, <1000 ions/μm^2^; therefore, it is unlikely that ions impact regions that have already experienced previous ion transmission. HCIs are generated by removing several electrons from an atom. Two key quantities define the HCI, the kinetic energy and the potential energy, which for the latter equals the sum of the binding energies of all of the removed electrons. Upon approaching a surface, HCIs begin to neutralize via a two-center Auger–Meitner process known as interatomic Coulombic decay (ICD). ?−? ? ICD is strongly distance-dependent, going as ∼1/R ^6^, with R being the distance between the HCI and the material atoms. ?,? Hence, the charge exchange process is sensitive to the material thickness and subnanoscale structure.

HCI transmission spectroscopy of CNMs can be used to study the structure of the membranes by analyzing the outgoing charge states and scattering patterns. High charge exchange (and large scattering angles) gives insight regarding the thickness and density of the membranes. Less charge exchange (and small scattering angles) gives access to regions of lower density, and even atomic-scale voids, in the sample. As with most spectroscopy methods, a model is required to infer sample properties from the measured spectra. In the case of HCI experiments, a model of the target’s atomistic structure is necessary to link the experiment to HCI transmission simulations.

The disordered structure of CNMs poses a challenge not only for experimental characterization but also for theoretical descriptions, due in part to the limited experimental data. Molecular dynamics (MD) simulations can be useful for elucidating structure–property relationships at the nanoscale and have been successfully used to study carbon materials such as tetrahedral amorphous carbon,? carbon foams,? and fullerenes.? Recently, MD simulations have also been used to explore the SAM-to-CNM conversion process. ?−? ? Ehrens et al.? reported the first full MD simulations of carbon-only SAM-to-CNM conversion by imparting momentum impulses to the SAM and removing random molecules, considering biphenyl, terphenyl, and naphthalene thiol precursor SAMs.? The authors report that momentum transfer simulations are able to produce cross-linked CNM structures with some degree of porosity; however, the predicted tensile moduli of the CNMs are significantly higher than experimental values. While these prior efforts were able to produce promising CNM structures, they could not recover the experimental tensile modulus.

To address these limitations, we employ two distinct MD simulation approaches to produce candidate terphenylthiol precursor CNM structures. The first method uses no spatial information on the precursor SAM. Instead, a predefined hole structure is enforced via modified particle dynamics and used in conjunction with a liquid-quench and annealing simulation protocol to produce porous CNM structures with various degrees of porosity and carbon network development. With this approach, CNM structures were created with varying degrees of subnanometer porosity and sp ^2^ fractions ranging from 40% to 90%. The second approach is an extension of the momentum transfer simulation technique reported by Ehrens et al.,? where discrete momenta are imparted to the carbon-only terphenylthiol SAM structure.? We herein refer to the former as the “exclusion cylinder” method and the latter as “momentum transfer” simulations.

The modeled CNM structures were then used as input for HCI transmission simulations. These calculations were performed using the time-dependent potential (TDPot) method,? which was modified to be compatible with non-2D target materials. The modified TDPot description of HCI neutralization is based on the ICD model and is able to predict angle-resolved exit charge state distributions in agreement with previous experiments of highly charged Xe transmission through graphene.? The simulated angle-resolved charge exchange spectra were then directly compared with the experiment, thereby placing further constraints on the true atomistic structure of the terphenyl-based CNMs. Our combined MD simulation and HCI spectroscopy approach is a unique tool that can be useful for characterizing other materials in the broader class of thin, short-range ordered organic nanosheets.

Methodology

Experimental Methods

HCI transmission experiments were performed at the ion spectrometer at TU Wien and are discussed in detail in refs ? and ?. A Dreebit Dresden EBIS-A ?,? was used, capable of producing Xe^ q+^ ions in charge states up to q = 44 at a fixed acceleration potential of 9 kV, resulting in ions with kinetic energies of E kin = 9 × q keV. A Wien filter for charge-state selection, a set of deflection and focusing electro-optical components, and a pair of slits guide and shape the beam toward the terphenylthiol-based CNM samples. The CNM samples used were prepared according to refs ? and ?.

The HCI experiments were performed under ultrahigh vacuum conditions (<1 × 10^–9^ mbar) in transmission geometry, i.e. ions impinged normal to the CNM surface, and were detected after transmission on a 2D position-sensitive multichannel plate (MCP) detector. A slit and a pair of deflection plates with a field orthogonal to the slit were located between the sample and the MCP to separate distinct charge states for analysis. This arrangement permits the study of the scattering angle ϕ on the horizontal axis and the charge state of ions along the vertical axis of the 2D MCP spectra.

Contaminant Removal by Electronic Excitation

Standard surface cleaning techniques cannot be easily applied to very thin samples such as CNMs, as high-fluence sputtering can lead not only to contaminant removal but also to the destruction of the sample itself. Heating the sample is also only possible up to a limit, as the CNM undergoes a phase transition to graphene at temperatures above ∼900 K. ?,?,? To overcome these challenges, we apply slow HCI irradiation and spectroscopy at different incident Xenon charge states. With low charge states (≲q = 10) membranes were cleaned by means of soft desorption induced by electronic transitions. ?,? Higher charge states were then used to probe the angle-dependent charge exchange patterns and assess the cleanliness of the CNM.?

HCI neutralization and the accompanying potential energy deposition induce electronic excitations into the sample. It has been shown that small electronic excitations from low-energy electron irradiation and ultraviolet photons can lead to contaminant removal from single-layer graphene. ?−? ? For HCIs, the potential energy deposited to the sample can be fine-tuned by selecting the incident charge state, with energies ranging from several eV up to hundreds of keV per ion. We navigate a fine line between contaminant removal at low charge states and damage to the sample at high charge states. Depending on the electronic properties of the probed material, nanopore formation can occur upon HCI impact, i.e., pores with diameters up to several nanometers can form in the sample surface at the ion impact position.? Many materials, including CNMs, exhibit a potential energy threshold under which there is no HCI-induced structural damage. ?−? ? We conclude that damage by Xe ions with q = 10 and a kinetic energy of 90 keV is sufficiently small as long as the applied fluence is less than 10^12^ ions/cm^2^, which corresponds to a fluence of ≲10^–3^ ions per carbon atom. Note that conventional kinetic sputtering by momentum transfer from the ion to the carbon atoms is limited to only one or two atoms per incident ion since an extended collision cascade that typically drives sputtering is not possible in a freestanding nanometer-thin membrane.

The exit charge state q out distribution was monitored as a function of the applied fluence. Figure(a) shows the exit charge state spectrum before cleaning, consisting mostly of neutrals and ions with q out = 1 and 2, all scattered up to large angles. This is a characteristic charge exchange pattern of thick materials where multiple scattering leads to large scattering angles and ions can fully neutralize.? Panel (b) shows the exit charge state after sample preparation, where all charge states from incident q = 10 through to neutrals are visible, and the majority were scattered to scattering angles ϕ < 0.2°. The appearance of more charge states and the overall decrease in scattering angles during sample cleaning are consistent with sample thinning, meaning contaminants were removed from the surface.

Exit charge state spectra of 90 keV Xe10+ after transmission through the CNM before (a) and after (b) contaminant removal through ion irradiation. (c) A fluence of ∼3 × 1011 ion impacts per cm2 is necessary to reach an equilibrium mean charge state of q out ∼ 2.

Figurec presents the mean charge state as a function of the applied fluence for the lower half and broader distribution (neutral Xe to Xe^4+^ in blue), and the upper half and smaller distribution (Xe^5+^ to Xe^9+^ in red). Xe^10+^ ion data were omitted from this analysis as it includes ions transmitted through micrometer-scale cracks.? For the low charge state distribution (blue in (c)), a continuous increase in the mean exit charge state with fluence can be observed before it plateaus at ∼3 × 10^11^ ions/cm^2^ (∼0.003 ions/nm^2^). These data support the idea that the membrane is not significantly damaged during ion irradiation and that only physisorbed atomic and molecular species were removed. If this were not the case, ions would constantly remove atoms from the material, making it thinner and leading to a monotonic increase in the mean exit charge state. Thus, after ∼3 × 10^11^ ions/cm^2^ the sample is deemed “clean”. Note that this fluence corresponds to 3.6 × 10^–5^ ions per carbon atom in terphenyl-based CNMs.

Applying this cleaning procedure, all spectra exhibit a bimodal charge state distribution, which was also observed in previous studies that used nitro biphenyl thiol-based CNMs. ?,?,? Additional experimental data on nitro biphenyl thiol CNM cleaning, as well as a comparison of nitro biphenyl thiol and terphenylthiol membranes (from different suppliers), can be found in SI. All results presented herein were achieved using CNMs cleaned according to the description above.

Computational Methods

Structural Models of Carbon Nanomembranes via Molecular Dynamics

Simulations

MD simulations were performed using a combination of LAMMPS? and an in-house MD simulation code, both of which used the carbon EDIP potential to describe particle interactions.? The exclusion cylinder capability, which is detailed in the following section, has been implemented in the stand-alone carbon EDIP MD code and has been made available for use under the GPL 3.0 license.? Structures were visualized using the OVITO software package,? and analyses were performed using in-house codes.

Exclusion Cylinder Molecular Dynamics Simulations

The exclusion cylinder MD simulation method was developed as a means of studying possible CNM structures by enforcing a distribution of voids during the structural evolution of the carbon network. Specifically, we explored the effect of increasing the areal hole density for a uniformly distributed set of holes with a fixed number of carbon atoms and unit cell dimensions (and therefore initial membrane thickness). It should be emphasized that these simulations do not aim to model the SAM-to-CNM formation process, and that we instead use a “means-to-an-end” approach to produce candidate CNM structures with various pore arrangements.

For each exclusion cylinder configuration, N _ p _ cylinders were generated and placed about the xy-plane within the unit cell. Positions of the circles were then moved using a simple hard-sphere-like dynamics algorithm which terminated when no circles were within the specified interpore distance cutoff. This cutoff was varied based on the cylinder count to maintain a roughly homogeneous distribution of exclusion cylinders. These cylinders then defined the regions that atoms were not permitted to enter during certain stages of dynamics. Nine unique configurations were used: 0, 10, 30, 50, 70, 90, 110, 130, and 150 cylinders, which correspond to an exclusion cylinder density of up to 0.84 nm^–2^. The exclusion cylinders extend infinitely into the principal z-coordinate, i.e., perpendicular to the membrane surface, and obey the periodic boundary conditions in x and y. Periods of dynamics that incorporate the exclusion cylinders are referred to as region-restricted dynamics. It should be noted that the exclusion cylinder algorithm treats the carbon atoms as point particles with zero size; as such the diameter of the exclusion cylinders is not directly comparable to any experimental measure of porosity.

Each simulation begins with the same number of carbon atoms (all other elements were omitted) and unit cell dimensions, i.e., 14040 atoms in a fully periodic (135 × 135 × 12) Å^3^ unit cell, which corresponds to the experimental terphenylthiol SAM density of ∼1.3 g/cm^3^. Initial atom coordinates were randomly assigned within the allowed region for each configuration, so that no atoms were initially located within the exclusion cylinder locations. In other words, each configuration of enforced exclusion cylinders begins from a unique starting configuration of randomly placed carbon atoms within the allowed volume. A mean exclusion cylinder radius of 2.5 Å was used, sampled from a Gaussian distribution with σ = 0.5. This value was selected based on recent gas permeation experiments ?,? that show molecules with kinetic diameters under 3.5 Å are able to permeate through terphenyl-based CNMs, as well as prior ion selectivity experiments? and water transport data. ?,?,?,? An integration time step of 0.35 fs was used for these simulations.

To prohibit atoms from entering the exclusion regions during region-restricted dynamics, the following routine was performed for all atoms at each integration time step. After the forces and updated positions have been calculated, we check if x,y coordinates of any atoms lie within the radius of any exclusion cylinders. If so, the velocity is reflected as

$$[eqn]$$

where v _ x _ and v _ y _ are the velocity components in x and y, respectively, and θ is the angle between v _ x _ and the x-coordinate of the unit cell. Figure schematically illustrates the specular reflection in several steps. This specular reflection conserves total energy, and hence, the standard equations of motion may be used.

Schematic diagram of the specular reflection imposed on atoms that enter the exclusion cylinders during region-restricted dynamics.

The overall simulation protocol for each hole configuration was as follows. Starting with the carbon atoms randomly packed in the unit cell, each hole configuration was subjected to region-restricted dynamics with full 3D periodic boundary conditions for (i) 3 ps in the constant particle, volume, and energy (NVE) ensemble, and (ii) an instant quench to a temperature of 300 K which was held for 3 ps using the Bussi thermostat.? Periodic boundary conditions in the z-direction were then switched off, and a vacuum gap was used to create a thin slab with a surface normal to z. This was followed by (iii) 3 ps of simulation thermostated to 300 K to allow the surface to relax and (iv) a linear increase in thermostat temperature up to 3000 K over 10 ps, i.e., at a heating rate of 270 K/ps, and then further maintained at 3000 K. Nine structural configurations were branched off after 0, 1, 4, 9, 16, 25, 36, 49, and 64 ps of simulation time. Once branched off, the structures are considered to be independent models. For each of these simulation branches, the exclusion cylinder restrictions were then switched off, and atoms were permitted to move freely about the unit cell. It is important to note that the annealing does not have any physical analog to the experimental preparation of CNMs. As all molecular simulations are limited by the time scales that are accessible to them, higher temperatures are used to accelerate the evolution of the system so that many carbon atom configurations can be studied in a computationally tractable manner.

Each branched frame was then immediately quenched to 300 K, and the pressure was maintained in the x and y unit cell directions at 1 atm via the Nosé–Hoover barostat? for a further 300 ps. This additional simulation time after quenching allows the structure to relieve excessive internal stresses before tensile deformation simulations. With the structure relaxed, the coordination fractions were calculated by using a bond length cutoff of 1.95 Å. The resultant bonding topology of the CNMs then consisted of a major single bond network and several disconnected small chains and isolated atoms. All carbon atoms not connected to the dominant topology of the membrane were removed. The cell dimensions in x and y remained within a few Ångströms of their initial values during this relaxation process. A total of 81 unique candidate CNM structures were produced, spanning the high-temperature simulation time and the exclusion cylinder count parameter space.

Tensile moduli for each CNM structure were predicted by deforming the x or y unit cell direction at an engineering strain rate of 5 × 10^–3^ ps^–1^. An integration time step of 0.1 ps was used for these tensile deformation simulations. During the deformation, the unstrained cell dimension was allowed to vary to maintain 1 atm of pressure. Stress–strain curves (Figures S5 in the SI) were then used to calculate the tensile modulus by applying a linear fit up to a strain of 0.02, where the slope of the fit gives the tensile modulus.? Tensile moduli were calculated for both x and y deformation directions for each CNM structure, and hence, the predicted modulus was taken as the mean of these two values.

Momentum Transfer Simulations

An alternative MD simulation approach was used to explore the formation of pores in terphenylthiol SAMs without the explicit enforcement of holes. The momentum transfer method is a refinement of the MD simulations reported by Ehrens et al.,? and begins with the assumption that the initial correlations in the SAM may partially carry through to the resultant CNM. Specifically, these new simulations employ periodic boundary conditions in the lateral (x and y) membrane dimensions. Details regarding the method can be found in ref ?; however, the process is briefly summarized below.

Each simulation begins with the same ∼135 Å square unit cell. Unlike the region-restricted simulations, the initial coordinates of the 14040 atoms were arranged as a terphenyl SAM (with all elements other than carbon removed) on a gold substrate, oriented such that the SAM surface was normal to the z-coordinate. Gold was modeled using a Lennard-Jones wall, employing the C–Au parameters from ref ?. Figure S8 in the SI presents an image of the initial SAM unit cell.

Irradiation of the SAM was modeled by a series of discrete momentum transfer events, where select groups of atoms were imparted additional momentum to mimic collective electron impacts. To this end, the SAM unit cell was divided into 1444 square regions, resulting in approximately 10 carbon atoms per region. A single momentum transfer event imparts a primary downward force to all atoms within a randomly selected region, i.e. along the −z direction. Regions that neighbor this primary region were subjected to a secondary planar force that was directed away from the primary region. For these simulations, a primary force of 50 eV/Å was used. Secondary force magnitudes and the number of force events were varied. This applied force was exerted for a single time step only, which was then followed by a relaxation period of regular NVE dynamics for a time of t _ r _. For a given structure, a loop of N _ e _ impact events was performed, followed by a period of regular NVE dynamics for 200 fs, and a further 35 ps of thermostated NVT dynamics at a temperature of 300 K, using the Langevin thermostat.? To mimic a cooling gradient through the gold surface, the thermostat is applied only in a defined region near the minimum of the Lennard-Jones potential. Some example results of this method are shown in Figure S8 in SI.

Pore Detection

A pore detection algorithm based on image processing techniques was developed to characterize the pore distribution in the CNM structures. To quantify the pores that would be “visible” to incident ions, xy coordinates of all atoms were binned into a spatial grid of pixels with a bin width of ∼0.25 Å along the x and y axes. Each pixel was then marked as either occupied or vacant, depending on if xy coordinates of any atom were present within the pixel or not. This array of pixels was then converted into a bitmap for further processing. At this stage, pores were visible as white areas and dense regions were shown as black areas. Due to the relatively small bin widths, white noise was present in the dense regions.

Two morphological operations from the OpenCV image processing library were used to remove noise pixels and determine the pore edges: “erosion” and “dilation”.? This combined erosion and dilation is also referred to as ‘opening’. In both operations the CNM pore image was convolved with a 17 px × 17 px elliptical kernel S, with a central anchor point. The size and shape of this kernel were tuned by inspection, such that noise no longer affected the final pore counts. Erosion was applied first to de-noise the dense regions by scanning S over the image and comparing the kernel pixels to overlapping pixels in the image. If they match completely, the image pixel under the anchor point retains its value, otherwise the pixel is eroded (set to zero). This process eliminates noise pixels but simultaneously thickens the edges of pores artificially. To correct for this edge thickening, a dilation operation was applied where S was once more scanned over the image, but this time the anchor point value was set to one if at least one element covered by S was non-zero. With the bitmap de-noised, contour detection was performed using the Suzuki and Abe method,? taking the periodic xy boundary conditions into account. The pixel count within the identified closed contours was then used to determine the pore area, with the diameter estimated from a circle of equivalent area. Figure S6 in the SI illustrates this process for a sample CNM structure. Total CNM porosity was then taken as the percentage ratio of pore area to surface area of the unit cell.

Highly Charged Ion Transmission Simulations

Charge exchange and scattering of HCIs transmitted through CNMs were modeled using the time-dependent potential (TDPot) approach of Wilhelm and Grande.? The original TDPot method was developed to study the stopping and scattering of HCIs transmitted through 2D materials, e.g., graphene, and as such had to be adapted for targets of finite thickness. To this end, parameters that control the ICD rate ?,? were tuned by comparing the simulated ion charge state distributions to benchmark experimental data for single-layer, bilayer, and trilayer graphene (SLG, BLG, and TLG, respectively). Figure S2 in the SI presents the mean exit charge states comparison between simulation and the experimental fit functions for SLG, BLG, and TLG.?

Simulations of 72 keV Xe^8+^, 135 keV Xe^15+^, and 180 keV Xe^20+^ ions transmitted normal to the CNM surfaces were performed. For each set of ion charge and kinetic energy parameters, 150000 individual ion trajectories starting at random x and y points within the unit cell were used to ensure sufficient sampling of the CNM. Exit charge state of the ions and deflection angle data were used to create exit charge state histograms to compare with the experimental data.

Results and Discussion

A total of 81 CNM structures were produced using the region-restricted MD simulation approach. Nine exclusion cylinder arrangements were used with a mean radius of 2.5 Å and a varying hole count from 0 through to 150 cylinders (which corresponds to a hole density of 0.84 holes per nm^2^). As discussed in the Computational Methods section, the term “hole” refers to the enforced exclusion cylinders that were used as part of the simulation method, and so the terms “hole count”, “hole density”, and “hole diameters” refer only to this exclusion cylinder construction and not the resultant membrane.

Figure presents a top-down view of all candidate CNM structures produced using the exclusion cylinder method with atoms colored according to their coordination: sp in blue, sp ^2^ in orange, and sp ^3^ in purple. All exclusion cylinder configurations began with roughly similar coordination fractions: 50% sp, 40% sp ^2^, with trace amounts of sp ^3^. During simulation, the coordination fractions progress similarly. There is a slight increase of sp carbon up until 9 ps, followed by a rapid decrease and plateauing below 15%. Note that here we use the term “sp” imprecisely, as 2-fold coordinated edge carbon atoms are in fact sp ^2^ radicals in bonding character, and only linear carbon chains would be true sp. However, for the sake of simplicity, we do not make this distinction in our analysis. For sp ^2^ carbon, there is a continuous increase at three distinct rates during simulation, and sp ^3^ carbon increases slightly from ∼2%, up to 4%. Graphs of the coordination fraction as a function of simulation time are shown in Figure S3 in SI.

Planar top-down view of a selection of exclusion cylinder-enforced carbon nanomembrane structures with atoms colored according to their coordination: sp = blue, sp 2 = orange, sp 3 = purple, and single-coordinated atoms are shown in green.

Tensile moduli for all structures shown in Figure were predicted from the linear regime of the CNM stress–strain response, where the tensile modulus is given by the slope of the linear fit. Then the CNM tensile modulus as a function of the number of exclusion cylinders is presented in Figurea, where the color of each line represents the simulation time of each of the structures. The shaded green region denotes the typical experimental tensile modulus of terphenyl-based CNMs, which has been reported to be 5–12 GPa.? Two overall trends are apparent: (i) for each hole configuration the predicted modulus increases with simulation time and (ii) the spread of tensile modulus values due to annealing generally decreases with an increase in exclusion cylinder count. Figureb presents a sample of the tensile data for low, medium, and high annealing times (1, 16, and 64 ps, respectively) as a function of the sp ^2^/sp ratio. Tensile modulus appears to be strongly correlated to the total sp ^2^ fraction of the membrane as the modulus generally increases with increasing sp ^2^ percentage, up to some threshold value for high hole densities. Structures with low hole counts, up to about 31 holes, do not exhibit any plateau in the tensile modulus.

Predicted average tensile moduli of the CNM models as a function of the number of exclusion cylinders (panel a), where the color of each point represents the annealing time to which each structure has been subjected. Each data point has been averaged over both in-plane directions, and error bars denote one standard deviation from the mean. The shaded green region indicates the reported experimental tensile modulus range. Panel b presents the same tensile data as a function of sp 2/sp for 1, 16, and 64 ps simulated annealed structures, where the enforced cylinder counts are indicated by the different markers and colors.

These data suggest that there are two ways of achieving a relatively (compared to that of pure graphite/graphene) low tensile modulus that is comparable to experimental CNM values. One way to achieve this is to have a more disordered structure, comprising significant amounts of under-coordinated carbon ≳30% sp. However, carbon structures with high sp ^2^ content can also result in a low modulus given a sufficiently high hole density, which then suppresses the in-plane stiffness of the membrane. Visual inspection of the highly annealed CNMs with 150 holes reveals that at high hole densities there is not enough space between pores for graphitic-like carbon to form lateral layers with the basal plane normal to the surface. Instead, the carbon structure tends to form single- or few-layered fullerenes that wrap around the exclusion regions with the basal surfaces normal to the void. In the extreme case, highly annealed membranes with high hole density resemble a forest of carbon nanotubes, interconnected with rows of sp ^3^ carbon. These results are consistent with previous theoretical reports of the transverse elastic constants for carbon nanotubes. Popov et al.? calculated the elastic moduli of single-walled carbon nanotubes of various diameters using force-constant lattice dynamics, and reported that nanotubes with a diameter of about 10 Å have a transverse tensile modulus of ∼10 GPa, which decreases for increasing tube diameters.

It is worth noting that a consequence of a membrane structure with large amounts of under-coordinated carbon is that it would be chemically reactive and that dangling bonds would most likely be stabilized by hydrogen and possibly water or other oxygen-bearing groups when exposed to atmospheric conditions. Our results therefore suggest that a nonreactive CNM with a tensile modulus of around 10 GPa could only exist if the CNM resembled a cross-linked forest of nanotubes. Such a structure is not compatible with the experimental Raman spectra data of terphenyl-based CNMs, which do not exhibit strong graphitic marker peaks, namely the D and G bands.?

HCI transmission experiments were performed using terphenylthiol-based CNMs with 72 keV Xe^8+^, 135 keV Xe^15+^, and 180 keV Xe^20+^ incident ion beams. TDPot simulations were then conducted for a range of exclusion cylinder CNM structures shown in Figure, yielding angle-resolved charge exchange spectra for incident Xe^ q+^ ions, which can be directly compared to the experimental data.

Simulations of 72 keV Xe^8+^, 135 keV Xe^15+^, and 180 keV Xe^20+^ incident on the CNM model structures were performed. Figure presents xenon exit charge state histograms for the 0, 30, 90, and 150 cylinder CNM structures for selected simulation times. Exit charge state distributions for structures annealed for 0, 16, and 64 ps are shown in darkening shades of blue to illustrate the progressive effect of annealing. Here, we highlight the structures annealed for 9 ps, as after this time structures begin to exhibit a bimodal exit charge state distribution. Lower values of exit charge (i.e., more neutral) are attributed to ions that pass through dense material and undergo significant neutralization, whereas the higher charge states are attributed to ions that pass through less-dense, or even open (relative to the ion ICD interaction radius), regions of the CNM. The transition from a single skewed distribution to a bimodal distribution occurs when the model CNM structure develops significant layering, which then presents distinct single and bilayer regions to incoming ions. Additionally, structures annealed for more than 9 ps are in general much stiffer than the expected experimental tensile modulus.

Simulated exit charge state for Xe ions incident on various hole-enforced nanomembrane structures. Structures that have been annealed for 0, 16, and 64 ps are shown in dark blue shades, and the 9 ps structure is shown in orange. It should be noted that since there is only one unique CNM structure for each configuration, there is no statistical uncertainty evaluation in the calculated exit charge state values.

The incident charge state is also prominent in many of the spectra shown in Figure. This particular charge state increases with increasing exclusion cylinder count, which is due to the increase of open regions within the CNM. These results show that the relative peak height of the incident charge state is sensitive to subnanoscale porosity and that the ratio between this and the mean of the lower charge state distribution could yield information regarding the total porosity of the CNM. Unfortunately, it is currently not possible to extract such detail from experimental charge exchange spectra, as there is always an overwhelming contribution of incident charge states due to microscale cracks in the samples. However, careful examination of high exit charge states may still yield valuable information regarding the porosity of the CNMs. In fact, Figure shows that the simulation captures the higher charge state distribution well, which is direct evidence of an innate subnanometer porosity.

Experimental (left half) and simulated (right half) charge exchange spectra through carbon nanomembranes for 72 keV Xe8+, 135 keV Xe15+, and 180 keV Xe20+ from left to right. For the simulations, we used a 150 exclusion cylinder membrane annealed for 9 ps shown far right in a top-down view, with atoms colored according to coordination, as per Figure . Here, the left half of the image is rendered using a bond-based visualization and the right half uses a van der Waals representation with an atom radius of 1.7 Å. Note that the scattering angles are different in experiment and simulation, details are given in the text.

Figure presents a comparison between the experimental charge exchange spectra and the TDPot simulation result using the best-fitting CNM structure, which was the 150 exclusion cylinder CNM annealed for 9 ps. A top-down view of the CNM with the usual coordination-based color scheme with two types of visual representation is also shown: a bond-based visualization on the left and a van der Waals representation with a carbon radius of 1.7 Å on the right. In general, the angle-resolved charge state spectra agree well with the experimental data. Both exhibit a distinct bi-modal charge exchange distribution where high charge states scatter within a narrow band of angles (<0.1°), and low exit charge states are significantly broader, scattering up to 0.8°. However, some discrepancies are present. In general, the simulations exhibit less overall neutralization compared to the experiment. For Xe^20+^, there is a more pronounced separation between the upper and lower charge state distributions in the experiment compared to the simulations.

These data in combination with the predicted tensile modulus suggest that low simulation time CNMs with higher numbers of enforced cylinders may be good candidate models for terphenylthiol-based CNM. The low charge state distribution can be linked to ions transmitting through densified parts of the CNM. The densification of the membrane and minimal carbon loss (≲ 5%) already imply a certain porosity of the resulting membrane, and ion charge exchange in these parts of the membrane is therefore indirect evidence of subnanometer porosity. Our TDPot simulations in comparison to angle-resolved experimental data (cf. Figure) also enables a quantification of which potential structure yields the best match to the densified part of the membrane. The slight discrepancy with experiment can be put into context by examining the key the assumptions in modeling process. The charge exchange simulations are dependent on two inputs: the CNM structure modeling and the description of the ion-surface interaction.

Key aspects of CNMs, such as the poorly understood reaction pathway from pre-cursor to CNM, lack of long-range order, the nano lengthscale of their features, and multi-element precursor composition, all pose significant challenges to MD simulation efforts. No current method is able to model the entire conversion process from the well-known and multi-element monolayer structure, through to the final, less-defined CNM. However, our method explores a significant fraction of the possible phase space of CNM structures within the confines of current experimental data, e.g. tensile modulus and density. The candidate CNM structural models that we have produced and highlighted as a best-fit to experimental data capture the key experimental aspects of terphenylthiol-based CNMs. Despite this, it remains difficult to determine if our models completely describe the actual CNM samples that were used in the charge exchange experiments, given the complexity of the SAM-to-CNM conversion process, and the wet-chemical transfer process involved in sample preparation.

The TDPot method, which is used to model the interaction between an ion and solid matter, utilizes benchmarked nuclear scattering potentials which are screened according to the time-dependent charge state of the projectile ion.? The rate of charge state neutralization is mediated by a distance-dependent function. The long-range part of this function is the main source of uncertainty within the model itself. The short-range behavior of the charge neutralization function was fit against experimental charge exchange data for single-, bi-, and tri-layer graphene. Therefore, some uncertainty remains regarding the long range amplitude and shape of the charge neutralization function. It may be possible to refine this charge neutralization description in the future by using experimental charge exchange spectra of HCIs on materials with nanoscale pores and well-defined, experimentally verified, atomistic structures. However, such experiments lie beyond the scope of this work.

While the exclusion cylinder-enforced MD simulations were able to provide insights regarding the possible atomistic structure of CNMs, the question remains whether it is possible to form a porous carbon network without the constraint of forced cylinder geometry. To explore this idea, additional MD simulations were performed using the momentum transfer method, where localized force events and light annealing are used to induce pore formation in a carbon-only terphenyl SAM. Here, two simulation parameters were varied to systematically study the formation of pores: (i) the so-called secondary force applied to the SAM at each impact event (from 300 to 550 eV/Å), and (ii) the number of impact events (30 through 50). The former will be referred to as “force-varied” and the latter as “event-varied”. The resultant pore structure of these simulations was then analyzed as well as the hole-enforced structures.

Figure presents the pore analysis data for the event-varied simulations (panel a), force-varied simulations (panel b), and exclusion cylinder 10 ps annealed CNMs (panel c). The distribution of pore areas is shown by the shaded blue violin-shaped regions. For context, a pore area of 30 Å^2^ would correspond to a pore radius of about 3 Å, assuming a circular pore geometry. Absolute pore counts are also shown in the same plots marked by the solid red triangles, with values indicated on the right vertical axis. For both event-varied and force-varied simulations, pore area distributions appear as either skewed normal distributions or bimodal with a dominant distribution of smaller pores with an area between 5 and 80 Å^2^, and a minor distribution of larger pores above ∼100 Å^2^. The overall range in pore areas is, in general, smaller for the hole-enforced simulations (panel c), which is intuitive given that hole geometry was contrived to an extent.

Pore area distributions shown in blue-shaded regions (left y-axis) for event-varied (a), force-varied (b), and the 4 ps annealed region-restricted simulations (panel c), with the total pore count shown in solid red triangles (right y-axis). Mean pore area is indicated by horizontal blue lines for each distribution. Note that panels a) and b) have different y-axis scales compared to panel c).

Considering the event-varied data in Figurea, a somewhat counterintuitive result can be observed. Here, the pore count reduces from roughly 60 and plateaus at 30 for an increasing number of force events, with 40 events seemingly being a threshold for producing fewer but larger pores. Visual inspection of the trajectories revealed that below the threshold event count of about 40 events, pores were primarily formed due to the overall low carbon atom density in the membrane, resulting in many smaller pores as the carbon structure anneals and develops. Regions near event sites tended to densify as carbon atoms were pushed away from the locus of the force event. Hence after a sufficient number of force events, the remaining material had enough carbon atoms to feed the structural evolution during annealing, resulting in fewer small pores. Figureb indicates that pore count also tends to decrease with increasing force magnitude. Pore area distributions, however, do not seem to follow any obvious patterns, with distributions similar to those seen in Figurea.

In the case of the enforced cylinder CNM structures at a simulation time of 9 ps, pore area distributions shown in Figurec have a lower spread compared to those in panels a) and b), with pore areas ranging from a few square Ångström up to about 300 Å^2^. For low cylinder count CNMs, the mean pore area is approximately 30 Å for all cylinder counts. When assuming a circular pore shape, this corresponds to a pore diameter of 5.5 Å (measured from carbon nuclei centers, and therefore excluding the effective electron cloud size), similar to the exclusion cylinder radius. It is interesting that all exclusion cylinder counts used result in similar final pore counts after relaxation without the exclusion region constraints.

The momentum transfer simulations in general produced broader and more complex pore area distributions, when compared to the constrained region-restricted CNM structures. These results suggest that there may be multiple pore creation mechanisms that occur during electron irradiation of the SAMs. As the conversion process from SAM to CNM is not yet fully understood, it is difficult to determine what form the pore distribution should take. Unfortunately, information regarding the pore size distributions cannot yet be extracted from the HCI charge exchange data. A deeper understanding of the SAM-to-CNM conversion process is still required to understand the pore area distribution further within the limitations set by the experimental data. Future work should explore this aspect of CNM structure and its influence on the resultant structure–property relationships.

Conclusion

The atomistic structure of CNMs was explored by using a combined experimental and simulation approach. HCI charge exchange spectra for cleaned CNMs exhibit a characteristic distribution of angle-dependent charge exchange patterns. Ions exiting the membrane in high charge states are not significantly scattered, whereas ions in low outgoing charge states are deflected more. The high charge state distribution can be clearly linked to ions passing through thinned and void regions in the membrane, whereas ions in the lower distribution pass through dense material regions.? The lower distribution is sensitive to the specific material (areal) density, and under the assumption that only a negligible amount of carbon is lost during CNM cross-linking, inhomogeneities in membrane density are an additional independent indicator of subnanometer porosity.

Two independent MD simulation techniques were used to produce candidate models of the CNM. One method explored possible carbon network configurations for differing degrees of porosity, while the other elucidated a possible mechanism of pore formation due to momentum impulses that might occur during precursor SAM irradiation. We link the candidate atomistic models from MD to the observed charge exchange spectra in HCI transmission with the help of a model based on classical equations of motion of the ions propagating through a given atomic structure. A quantitative comparison between the experimental and simulated ion transmission spectra for the different MD structures reveals that a subnanometer porous structure with a significant fraction of under-coordinated carbon is likely. This is further corroborated by comparing the predicted tensile moduli to the literature data. The MD simulation techniques that we have developed can be used to study other porous systems beyond that of pure carbon, both for structure modeling, and subsequent application-specific process modeling, e.g. water transport through membranes.

The CNM structures that best fit the experimental touchpoints would likely be chemically reactive at ambient conditions and could be passivated by hydrogen, as well as water and other oxygen-bearing groups, after the CNM has been removed from vacuum. This aspect of CNM production has seldom been discussed in the literature. We suggest that structure passivation is a crucial aspect of the CNM, stabilizing the reactive carbon network that is left behind after electron irradiation of the precursor SAM, and may even play a role in the functional properties of the CNM. Our study is the first to give quantitative insight into the atomistic structure of CNMs from coinciding experiment and simulation, paving the way to a comprehensive understanding of the mechanical and chemical properties of organic short-range ordered membranes.

Supplementary Material