## Abstract

We extend the three-point XPFC model of Seymour & Provatas (Seymour & Provatas 2016 *Phys. Rev. B* **93**, 035447 (doi:10.1103/PhysRevB.93.035447)) to two components to capture chemical vapour deposition-grown graphene, and adapt a previous two-point XPFC model of Greenwood *et al.* (Greenwood *et al.* 2011 *Phys. Rev. B* **84**, 064104 (doi:10.1103/PhysRevB.84.064104)) into a simple model of two-component graphene. The equilibrium properties of these models are examined and the two models are compared and contrasted. The first model is used to study the possible roles of hydrogen in graphene grain boundaries. The second model is used to study the role of hydrogen in the dendritic growth morphologies of graphene. The latter results are compared with new experiments.

This article is part of the theme issue ‘From atomistic interfaces to dendritic patterns’.

## 1. Introduction

Graphene is a two-dimensional carbon allotrope with a honeycomb structure. Its many potential applications are well documented [1,2]. Defects in graphene are known to impact its functionality in these and other applications [3]. Understanding how to control defects during graphene growth is thus crucial [4].

Graphene can be grown via chemical vapour deposition (CVD) on nickel [5] and copper substrates [6]. The partial pressures of H_{2} and CH_{4} gases streamed into the system lead to adsorption of H and C onto the metal surface. Depending on these pressures, the surface carbon can nucleate dendritic graphene flakes or graphene islands, with a coverage of approximately one monolayer thick. Coexistence can be found between graphene islands and disordered regions of C and H [4,7]. It is observed that hydrogen gas is necessary to produce dendritic graphene morphologies. There have been experimental predictions made about hydrogen’s role in the CVD process [8,9], but there are still aspects about its interactions with carbon in graphene growth that are unknown.

While experimentally the carbon is easy to image, imaging hydrogen in CVD graphene growths [4] is difficult. Moreover, this topic is difficult to do theoretically due to the many length and time scales involved in graphene microstructure.

First-principle calculations have been used to observe the role of hydrogen in CVD growths on copper using methane and hydrogen gases, and to study graphene morphologies on various metal substrates [10]. Molecular dynamics simulations have been used to examine CVD growths on nickel surfaces to predict carbon structures that form [11] and to calculate the energy of defect structures [12]. While these methodologies are very valuable on the atomic scale, their limited length and time scales preclude observations of extended polycrystalline domains and dendrites. Traditional phase-field models have recently been developed to simulate larger scale epitaxial graphene microstructure on copper [13]. However, these presuppose specific sharp interface conditions at the graphene interface and do not model atomic-scale defect structures.

The *phase-field crystal* (PFC) method has shown promise for simulating nano–microstructure evolution on diffusive time scales, while simultaneously incorporating several important phenomena emergent at the atomic scale [14]. The first phase-field crystal (PFC) model to capture honeycomb symmetry [15] does not allow for coexistence of the solid phase with a disordered phase, as required to model CVD graphene. The authors in [16] developed a two-component PFC model where two triangular lattices interact to a suite of interesting two-dimensional phases/defect structures. However, this approach only allows graphene-disordered coexistence for equal amounts of *A*/*B* components, and represents the graphene lattice with equal substitutional constituents of *A* and *B*, both conditions being unsuitable to study CVD graphene. Recently, a single-component structural PFC (XPFC) model was developed that employs a rotationally invariant three-point correlation, and models crystallization of graphene from a vapour phase over a wide range of densities. It produces polycrystalline grain boundaries with 5–7 defect structures [17,18].

This paper extends the three-point XPFC model of Seymour *et al.* [17] to two components to capture some features of CVD graphene growth, which can be represented as a carbon lattice with some fraction of adsorbed hydrogen. We will also adapt a previous two-point XPFC model of [19] into a simple model of two-component graphene. The equilibrium properties of these models are examined. The first model is then used to study the role of the hydrogen component in grain boundary dissolution, while the second model is used to study the role of hydrogen in dendritic graphene morphologies, and our results are compared with new experiments.

## 2. Three-point XPFC models for graphene symmetry

The present work extends the XPFC model of [17] to a two-component mixture, where one component favours a solid phase and the second component favours a disordered phase. The model is motivated by the process of CVD-grown graphene, where the primary component *A* (called ‘carbon’) and the secondary species *B* (called ‘hydrogen’) effectively mix both on graphene islands and in the disordered regions between them [20]. The XPFC methodology of [17] employs a novel, rotationally invariant, three-point correlation function that features an angular dependence that stabilizes honeycomb (and other) crystal structures, and allows their growth from a disordered state over a wide range of temperatures and densities. The two-point correlation function in [17] is based on hard-sphere interactions; however, it is straightforward to replace this by the broader Gaussians of the original XPFC approach [19]. We begin this section by briefly reviewing the pure material version of the model, and then extend it to two components.

### (a) Review of single-component model

Crystallization in the XPFC model is defined in terms of the dimensionless atomic number density , where *ρ* is the density, is a reference density in the disordered state of the system and ** r**=(

*r*,

*θ*) as we are in two dimensions. The Helmholtz free energy functional of this system is given by 2.1where

*F*

_{ideal}is the ideal free energy of cDFT expanded to fourth order, and given by 2.2The excess free energy,

*Φ*

_{excess}, is given by 2.3where

*C*

_{2}is a two-particle correlation function, and are components of

*C*

_{3}, the three-particle correlation function. The ideal free energy is an ‘ideal-gas like contribution’ that describes non-interacting particles. The excess free energy describes interactions between particles and is responsible for the emergence of structure during crystallization.

The term *C*_{2} is a repulsive term that sets the crystal lattice spacing. Its form is given by
2.4where *R* and *r*_{0} set the magnitude and length scale of the repulsion, and circ(*r*) denotes the indicator function on 0<*r*<1. To model graphene, we use *R*=6 and *r*_{0}=1.22604, while using a different choice of these parameters leads to crystals with triangular and square symmetries [17].

Three-point interactions are used here to select an angular dependence between nearest neighbours. The set the angular dependence and depend on this XPFC model’s temperature parameter, *X*, which is a measure of the strength of interaction between bonds. The components of the three-particle correlation function are split into a radial and angular dependence according to , where the radial component is given by
2.5where *X* sets the temperature, while *a*_{0} further sets the lattice spacing. The angular component is given by
2.6where *θ* is the angular separation and *m* is a parameter that controls which crystal bond angles are favoured. For threefold graphene crystals *m*=3. Here, we set *a*_{0}=1, which implies that distances are measured in units of *a*_{0}. Seymour & Provatas [17] investigated the properties of this model. It supports equilibrium phase coexistence between a disordered (vapour) phase and graphene, square or triangular phases. It was also shown that graphene grain boundaries feature experimentally relevant 5–7 defect structures, and bulk graphene responds self-consistently elastically [17].

### (b) Two-component model

We extend the model of §2a for two components by introducing the field *n*_{A} to describe the primary component (carbon), and a field *n*_{B} to describe the secondary component (hydrogen). They measure density relative to reference densities and . The total free energy is
2.7where
2.8and
2.9with *C*_{2} given by equation (2.4) and *C*_{s} by components in equation (2.5) and equation (2.6). The impurity species is assumed in this work to behave like an ideal ‘gas’, and its free energy is given by
2.10The above form assumes that the carbon favours a graphene structure and the hydrogen a disordered state. The interaction on a two-dimensional surface between species *A* and *B* is given by
2.11where *a* and *b* are constants that influence the type of interactions. The interaction kernel, *C*_{int}, is given by
2.12where *B* controls the magnitude of interaction and *α* adjusts the interface energy between the fields. In this study, *B*=2 and *α*=0.5. The function *C*_{int}(** r**) is peaked at

*k*=0. In this context, using

*b*=0 mostly couples the average of

*n*

_{A}to the local average of

*n*

_{B}, while the case

*b*≠0 is a minimal term that couples the ordering amplitude of

*n*

_{A}to the local average of

*n*

_{B}.

### (c) Equilibrium properties

Unlike previous XPFC models, rather than transform , where *n* is the local total density and *c* is the local *B* concentration, we will continue to work in the variables *n*_{A} and *n*_{B} as this representation was found to be convenient for investigating problems that consider inter-species interactions at the atomic level, such as CVD growth from a disordered phase. In this representation, a disordered phase has (*n*_{A},*n*_{B})=(〈*n*_{A}〉,〈*n*_{B}〉), where 〈*n*_{A}〉 and 〈*n*_{B}〉 are the average values of *n*_{A} and *n*_{B}. An ordered phase is one in which *n*_{A} and *n*_{B} assume different average values, plus some periodic part that has the periodicity of the lattice (i.e. graphene density was expanded in a Fourier series with reciprocal lattice vectors of four magnitudes). The amplitude of each periodic mode of the solid-state expansion is nonlinearly dependent on average densities and the parameters of the model.

We construct a phase diagram in (〈*n*_{A}〉,〈*n*_{B}〉) space. This amounts to extremizing the total free energy re-expressed in terms of average densities, and performing a common plane construction—the two-component analogue of the common tangent construction. This is done by a complex numerical algorithm we developed, which returns the convex hull over the region of coexistence defined by the free energy densities of the disordered and graphene phases. This is followed by the numerical construction of coexistence tie lines that terminate at two sets of equilibrium coexistence densities, and , where the superscripts o and d refer to ordered and disordered phases, respectively. The details of this general procedure are beyond the scope of this paper, but are found in [18].

Figure 1 shows two example phase diagrams at the model temperature *X*^{−1}=0.4. One is for the case *a*=*b*=1 and the other for the case *a*=1 and *b*=0. For the case of *b*=1, varying the average of the *n*_{B} field can fully dissolve the ordered phase of the *A*–*B* mixture as *n*_{B} couples to the crystal amplitudes of *n*_{A}. For the *b*=0 case, 〈*n*_{B}〉 plays a more passive role, coupling to the average of *n*_{A} and leading only to changes in the volume fraction of the ordered versus disordered phases. It is not yet known if either of these two scenarios of our model correspond to the experimental condition of CVD-grown graphene. While it is plausible to model CVD graphene as an effective mixture, because hydrogen chemically bonds with carbon in the ordered phase through adsorption [20–22], it is not known to what extent H atoms cover C atoms of graphene, and how this is affected by the pressure (density) of the H_{2} and CH_{4} gases. By the interaction term in equation (2.11), it is possible to model a wide range of mixing conditions. In the remainder of this section, the case of *b*=0 is considered. The main results of this subsection remain unchanged with the system described by *a*=*b*=1.

### (d) Dynamics and the evolution of graphene grain boundaries

The time evolution of *n*_{A} and *n*_{B} is simulated using model B type dynamics [23]
2.13where *t* is time and *M*_{r} is the ratio of mobilities of *A* and *B* atoms. The above equations should also have stochastic noise terms appended to each field which satisfy appropriate fluctuation dissipation statistics [24]. These simulate thermal fluctuations in each field. For simplicity we set *M*_{r}=1 here, and neglect noise in this work. Numerical simulations of equations (2.13) were conducted on a uniform square grid with d*t*=0.0002 and d*x*=0.2.

Figure 2 shows a graphene polycrystal grown with the present two-component XPFC model. The temperature was set to *X*^{−1}=0.4, and *n*_{A}, *n*_{B} were initialized uniformly with Gaussian fluctuations of zero mean and unit variance. The averages (〈*n*_{B}〉−〈*n*_{A}〉) were chosen in the solid region of the phase diagram corresponding to *a*=1 and *b*=0 in figure 1. Figure 2*a* shows the graphene structure set by *n*_{A} (carbon), while figure 2*b* shows that *n*_{B} (hydrogen) is largely disordered, with larger accumulations in the grain boundaries, consistent with experimental observations [20,21], which suggest that hydrogen chemically mixes (binds) with carbon within the bulk of graphene, but has a higher propensity for binding on defects. The larger hydrogen density in the grain boundaries is a consequence of how the *n*_{A} and *n*_{B} fields are coupled together in the free energy. This simulation begins with the nucleation of different crystal orientations from the disordered phase, followed by growth and impingement leading to grain boundaries. The insets of figure 2*a* zoom in on one region that emphasizes the bulk graphene lattice and another that shows a 5–7 grain boundary defect, which contains a ring of five atoms adjacent to one of seven atoms [3].

### (e) Dissolution of grain boundaries in the XPFC model

In CVD growth, CH_{4} and H_{2} gases are introduced in a low pressure chamber. Depending on the relative amount of each gas and the overall growth time, the surface graphene coverage ranges from isolated dendrites to full polycrystalline coverage of graphene grains [13,25,26]. In experiments we performed, it was observed that after partial coverage of the surface by graphene flakes was achieved during the above process, shutting off the CH_{4} gas and cooling the system while continuing to stream H_{2} through the system leads to ‘holes’ occurring within the bulk of some graphene grains [25,26]. We also observed that significantly increasing the H_{2} pressure in this situation leads to the disappearance of graphene flakes. These observations suggest that a high enough concentration of hydrogen may interfere with graphene growth and stability. This is also observed in other work [13]. It would be instructive to investigate further whether decreasing carbon density (pressure) at fixed H_{2} pressure causes dissolution of graphene at grain boundaries and leaves behind elevated regions of hydrogen there, due to H saturation and because of the preference of H to bind preferentially at defects [21]. This will be investigated experimentally in a future work. Here, we investigate the consequences of hydrogen on graphene grain boundaries in the context of our XPFC model.

We studied the role of changing carbon density in grain boundaries of 100% surface coverage of graphene grown by the present XPFC model [25,26]. The initial condition for these simulations is a graphene polycrystal state such as the one in figure 2. After this state is reached, the density of hydrogen (*n*_{B}) is held constant while the density of carbon (*n*_{A}) is decreased to mimic the experimental scenario described above. The average carbon density is reduced uniformly throughout the sample using a uniform flux source of the form appended to the *n*_{A} equation, where *f* controls the flux rate and is the desired target *A* (carbon) density. Here, we used *f*=1.95 and is typically chosen in the coexistence region of figure 1 (*b*=0 case), which determines the final fractions/concentrations of graphene and disordered phases.

Figure 3 shows the result of one of these simulations for the system corresponding to *X*^{−1}=0.4, *a*=1, *b*=0. Here the initial (〈*n*_{A}〉,〈*n*_{B}〉)=(1,0), and , while 〈*n*_{B}〉=0 is maintained. Figure 3*a*,*c* shows the simulation after 20 000 numerical time steps (units of *dt*) and figure 3*b*,*d* shows the density fields after 35 000 simulation time steps. As time proceeds, grain boundaries start to become dissolved because the density 〈*n*_{A}〉 moves towards the disordered state. At late times, the dissolved volume fraction will reach its equilibrium volume fraction. Based on the tie lines of the *a*=1, *b*=0 phase diagram in figure 1, increasing the concentration of hydrogen in the simulation will lead to larger regions of disordered phase (gaps between graphene flakes). The choice of *b*=0 is a minimal coupling that gives dissolution of grain boundaries as *n*_{A} is lowered (for fixed 〈*n*_{B}〉), consistent with some experimental observations. The more general coupling of *n*_{A} to *n*_{B} in equation (2.11) may also provide other paths to dissolution. For example, the case *b*≠0 may lead to complete dissolution of grain boundaries only by increasing *n*_{B} (at fixed 〈*n*_{A}〉). A precise choice of couplings for this system requires further experimental evidence.

## 3. Two-point XPFC model for graphene symmetries

We studied the growth morphology of individual graphene crystals experimentally as a function of hydrogen in the CVD process [25,26]. We found that depending on the flow of CH_{4} and H_{2} gases, we can either produce full coverage graphene or isolated dendritic graphene flakes. As has been mentioned previously, there is no established method of imaging the hydrogen in order to elucidate its role in affecting graphene morphology. However, from some of our observations, and that of others [13], it appears plausible that hydrogen interferes with the diffusion of carbon adatoms on a metal surface. This can thus have a destabilizing effect on the graphene interface due to the Mullins–Sekerka instability [27]. This section investigates this mechanism in the context of XPFC modelling, and compares the results with our experimental observations.

The model of §2 allows for graphene growth and coexistence with a disordered phase over a wide range of carbon/hydrogen densities. However, it is presently limited to small systems as it is computationally intensive due to high *k*-modes generated by its correlation functions. It is possible to soften these functions using Gaussian peaks in *C*_{2}(*r*) and *C*_{r}(*r*), or to filter out the high *k*-modes with a low-pass filter. Rather than do so in this work, we opt instead to investigate dendritic growth using the simpler two-component XPFC model of [19], which stabilizes graphene by inverting the sign on the cubic term in the ideal free energy. This is the two-component XPFC analogue of the model introduced by Hirvonen *et al.* [28]. The sign change is phenomenological—although it may potentially be argued to come from the low-wavelength contribution of a three-point correlation function [29]. This method of stabilizing graphene is less robust than the method of [17], and also leads to some incorrect elastic responses [18]. However, it will prove efficient for examining interface instability and dendritic evolution of graphene as a function of hydrogen.

The simplified mixture model follows [19] and introduces a total mass density field , where is a total reference density. A local concentration of *B* is defined as *c*=*n*_{B}/(*n*_{A}+*n*_{B}), where *n*_{A} and *n*_{B} represent the respective densities of carbon (*A*) and hydrogen (*B*). The dimensionless free energy of this model is given by
3.1The only change in the above model from that in [19] is the change of sign of the cubic term, as was done in [28]. As in [19], the entropy of mixing is defined as , where *c*_{r} is a reference concentration. The parameters *η*, *χ* and *ω* are introduced in order to fit the free energy away from and *c*_{r}. The parameter *α* is included to control the energy scale of concentration fluctuations. Here we set *c*_{r}=0.5, *η*=1.4, *χ*=1, *α*=1 and *ω*=0.2. The effective two-point correlation function, *C*^{n}_{eff}, is given by
3.2The term is the two-particle correlation function of *A* atoms, and the interpolation function *X*_{1}(*c*) is a phenomenological function [19] that modulates the strength of the *A* (carbon) ordering based on the local *c* (hydrogen concentration). As hydrogen favours a disordered state, we do not include for simplicity. We use for a triangular structure, given by the usual form, , where |*k*| is the norm of the reciprocal lattice vector ** k**. Following [19], we set

*d*=

*π*/2,

*k*

_{p}=2

*π*and

*β*=0.8. The parameter

*σ*is the XPFC model temperature. The other parameters control the peak location and width.

Figure 4 shows the numerically constructed phase diagram of this simple XPFC graphene model in the *σ*–*n*_{0} and *σ*–*c*_{0} planes, where *n*_{0} and *c*_{0} are averages of *n* and *c*, respectively. Details of this phase diagram construction are found in [25,26]. The slight inflection changes in the lines of coexistence at higher values of *σ* are probably due to numerical errors in the phase diagram construction. Figure 4 qualitatively corresponds to a small region around 〈*n*_{A}〉≈0 in the *a*=*b*=1 phase diagram of figure 1. Physically, this corresponds to the situation where carbon density is held close to some reference, while 〈*n*_{B}〉 is allowed to vary. Thus, increasing/decreasing 〈*n*_{B}〉 in the model of §2b corresponds to concomitantly increasing/decreasing concentration *c* and *n*_{0} in the XPFC model in this section. Figure 4 has stable graphene states at low hydrogen concentration and, conversely, disordered states at high values of hydrogen concentration. This is consistent with our experimental observations, where it was found that, at a fixed carbon density (fixed CH_{4} pressure), the growth rate of graphene dendrites increased/decreased when decreasing/increasing the hydrogen density (H_{2} pressure). Figure 4 also predicts that the disordered phase occurs at increasing temperatures (*σ* values).

Figure 5 shows a polycrystalline graphene sample simulated with the simple XPFC model. The dynamics of *n* and *c* were simulated by equations (2.13). Initial averages of *c* and *n* were set uniformly in the solid part of the phase diagram, and perturbed with random Gaussian fluctuations. Lighter regions correspond to higher density (concentration), and conversely for darker regions. Figure 5*a* reveals that the maxima in *n* form a honeycomb lattice in bulk regions, and chains of defects along grain boundaries of different orientations. The *c* field (figure 5*b*) has small variations about a smoother local mean. These variations are much smaller variations than in the *n* field. These variations in *c* have been analysed in other PFC models [29] and have been found to be of higher order than the smooth part of the *c*-field. The concentration attains above-average concentration values at grain boundaries. This simplified binary model also generates 5–7 defect structures at grain boundaries [25,26]. Moreover, analogously to the model of §2, when the hydrogen concentration was increased (thus decreasing the fraction of carbon in the system), the grain boundaries dissolve as shown in [25,26].

### (a) Summary of experimental dendritic growth

We grew dendritic graphene flakes experimentally using CVD with CH_{4} and H_{2} gases. Experiments were performed with the goal of obtaining isolated graphene dendrites. It is found that in order to produce dendritic morphology, the presence of hydrogen is required, otherwise only polygonal morphologies are produced. In the experiments reported in this work, the flow of carbon-carrying methane was kept at CH_{4}=1 sccm (standard cubic centimetres per minute), while the flow of hydrogen was varied from 2 to 80 sccm. The trends reported below on the role of hydrogen remained true in each trial and only a subset of experimental data obtained is presented. The complete data and analysis can be found in [25,26].

Figure 6 shows the results from trials at H_{2}=2, 20 and 40 sccm. The flow rates of CH_{4} and H_{2}, growth times and ‘*dendricity*’ (a parameter that qualitatively characterizes the degree of protrusion of dendritic branching) are shown in table 1. The dendricity is defined as the ratio between the inner and outer radius of a graphene flake. The inner radius is the radius of an approximate circle before any dendritic arms start to grow and the outer radius is the radius of an approximate circle that fits the entire graphene flake. For each H_{2} flow rate 10 optical microscopy images were taken of different flakes in the sample and the dendricity of each was calculated. As the growth time observed depended on the H_{2} flow rate, the dendricity parameter was normalized by the growth time. The averaged results of dendricity/min are reported in table 1.

Figure 6 reveals that a larger flow rate of H_{2} leads to fewer and slower dendritic flakes. Table 1 also shows that dendricity decreases with increased hydrogen flow, signalling a slower growth rate and increased stability of the graphene interface to perturbations. Conversely, a slower flow rate of H_{2} leads to faster crystal growth rates and more dendritic side branching. These results suggest that adsorbed hydrogen influences the driving force for graphene crystallization on the copper surface. These results are also consistent thermodynamically with the phase diagram in figure 4. Below we also examine the graphene growth kinetically using our simple XPFC model.

### (b) Dendritic growth simulated by the simplified XPFC model

To better understand the above dendritic growth experiments, we simulated dendritic growth using the simple XPFC model at different average hydrogen (*B*) concentrations. Interface instability occurs by the well-known Mullins–Sekerka instability mechanism [27], according to which the growth rate and curvature of dendritic branches is inversely proportional to the thermodynamic driving force. It is noted that experimental microstructures on the scale of ∼100 μm are difficult to do with PFC simulations without impractical use of supercomputing, because XPFC theories resolve structure at the atomic scale. To overcome this problem, we thus simulate crystallization from the disordered state at sufficiently large thermodynamic driving forces to give small dendritic structures that develop on short time scales. Our simulations here are on the scale of ∼10^{−1} μm. Based on traditional phase-field studies of solidification, we do not expect the trends reported here to change due to this difference. Larger length scales may also be possible using complex *amplitude* equations derived from the present XPFC models.

We simulated crystal growth using the simple XPFC model at different concentration (hydrogen) values and fixed temperature in order to compare our simulations with our experiments. Simulations were done with the XPFC temperature set to *σ*=0.2, and seeded with a small circular initial seed surrounded by a disordered state. The seed density *n* was initialized with an inverted triangular mode expansion, with the mode amplitude evaluated at the graphene phase boundary for *σ*=0.2. The average of the density *n*_{0} was set to the coexistence value in each phase. The initial concentration in each phase was set to *c*_{0}, the system average. Three values of *c*_{0} were investigated: *c*_{0}=0.13, 0.15, and 0.17. These values of *c*_{0}, time simulated and the simulated dendricity of dendrites obtained are listed in table 2. The dendritic morphologies of the last computed time step for each *c*_{0} are shown in figure 7. The data show that a lower hydrogen concentration (*c*_{0}) leads to higher dendricity (pointier and longer branches), consistent with our CVD experiments. Furthermore, the data of table 2 show that the lower *c*_{0} (i.e. higher driving force) leads to faster rates of crystallization, also consistent with experimental observations. Figure 7 also shows a concentration ‘halo’ around the dendrite, predicting that experiments should observe larger concentrations of hydrogen within inter-dendrite branches. Results are also consistent with other theoretical work [13].

## 4. Conclusion

We introduced two new PFC models to simulate crystallization of two-component solids with graphene symmetry. The first is two-component extension of the three-point XPFC model of [17], which treats carbon as an ordering component and hydrogen as an ideal gas. In this model, carbon (*A*) can interact with the local average of hydrogen (*B*), either weakly enough to cause hydrogen to limit graphene growth or segregate to defects and dissolve grain boundaries, or strongly enough to completely block the ordering of graphene. Its equilibrium properties were reviewed, showing it is able to emulate crystallization of graphene by CVD over a wide range of carbon and hydrogen densities. Simulations and experimental observations are used to provide directions for future experiments to probe the role of hydrogen in graphene grain boundaries. The second model examined is a simple adaption of a recent two-point XPFC alloy model that stabilizes graphene. This model is less robust or realistic than the first model but is computationally less stiff, thus allowing us to simulate crystallization on larger system sizes and longer time scales, required to prove dendritic growth. We used this simple XPFC model to show how decreasing/increasing hydrogen levels in CVD enhances/reduces the growth speed and dendricity of graphene flakes, results that are consistent with new experiments we reported.

## Data accessibility

This article has no additional data.

## Authors' contributions

M.S. developed all three-point correlation XPFC models, and the general methodology for computing two-component phased diagrams. K.L.M.E. carried out simulations of this model, developed the simplified XPFC model, and carried out the experimental work discussed in this paper with help from M.L. and M.H. N.P. designed and supervised this study and co-wrote the manuscript. All the authors read and approved the manuscript.

## Competing interests

The authors declare that they have no competing interests.

## Funding

We thank the National Science and Engineering Research Council of Canada and the Canada Research Chairs Program for funding and Compute Canada for high-performance computing resources.

## Footnotes

One contribution of 16 to a theme issue ‘From atomistic interfaces to dendritic patterns’.

- Accepted September 25, 2017.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.