Designing electronic properties of two-dimensional crystals through optimization of deformations

With their intrinsic two-dimensionality, 'electronic membranes' are easily pulled or pinched by atomic-scale tips [1–3], can be made to conform to the substrate topography [4–6], can be inflated as balloons [7], can be stretched [8] or bent [9], crumpled on demand [10], and so on. Hence, two-dimensional (2D) crystals are an excellent case of (and opportunity for) correlation between electronic behavior and shape with tremendous implications in bridging soft and hard condensed matter. For example, if a physical property is sensitive to the state of deformation of the system it can be used to monitor its shape, strain, etc; conversely, the shape variables can be manipulated so that the physical quantity in question behaves in a desired way, has a certain magnitude, or a particularly useful spatial profile. In addition, the fact that some of these 2D electronic membranes can be easily, and non-detrimentally, embedded in living tissues, organs or plants [11, 12], brings the tantalizing prospect of using them in bioelectronics. The method to be discussed next can be a valuable tool there, in the cases where the systemʼs functionality is determined by the shape or deformation state of the membrane.

To be specific—but by no means implying a limitation in scope—consider the problem of strained graphene. It is well-established that a mechanically-strained graphene sheet is very resilient [1] even in polycrystaline form [13, 14], and has altered electronic transport properties. In particular, and among other features, it exhibits an unconventional contribution in the electron–phonon coupling leading to the emergence of so-called pseudomagnetic fields (PMFs) [15–18]. These fields appear naturally in the effective (low-energy) description of the electronic problem in deformed graphene, and are a consequence of the peculiar lattice structure. Briefly, the celebrated Weyl–Dirac equation that captures most of the electronic phenomenology of graphene ($\mathcal{H}={{v}_{F}}\;{\bf p}\cdot {\boldsymbol{ \sigma }}$ for one of the K points in the Brillouin zone) is corrected in the presence of lattice deformations in a way that amounts to substituting ${\bf p}\to {\bf p}+e{\bf{\mathcal{A}}}$, where $\bf{\mathcal{A}}$ encodes all the details of the deformation and how it perturbs the electronic hopping amplitudes (defined below) [15]. As a result, even though $\bf{\mathcal{A}}$ is not a magnetic vector potential, the actual dynamics has the same characteristics and the Dirac electrons in graphene react to static and non-uniform lattice deformations as though they were under the influence of an effective magnetic field, with all the consequences that a magnetic field brings to electronic motion, except that time-reversal symmetry is not broken and, thus, $\bf{\mathcal{A}}$ will have an opposite sign for the effective Hamiltonian at the time-reversal transformed $K^{\prime}$ point. One such consequence is the modification of the electronic energy spectrum with the development of local Landau levels for certain lattice deformations. This has been recently confirmed by local scanning tunneling spectroscopy (STS) measurements on nanometer-scale graphene nano-blisters which revealed Landau level resonances associated with PMFs in the range 300–600 T [19, 20]. Equally interesting space-dependent Fermi velocities have also been reported in recent experiments on strained graphene [21], bringing this other theoretical prediction [22–24] and implication of non-uniform deformation fields closer to reality.

The possibilities associated with these discoveries and the confirmation of the drastic impact that moderate lattice deformations can have in grapheneʼs electronic spectrum have spurred researchers to investigate deformation modes allowing a degree of control over PMFs that can be tailored for specific ends, such as electronic confinement, guiding, and so on. This is a concept known as strain engineering or straintronics [6, 25–28].

Since the electronic dynamics can be straightforwardly determined once a space-dependent (pseudo)magnetic field $B({\boldsymbol{X}} )$ is prescribed, and since much is already known about the behavior of Dirac electrons in graphene under the influence of magnetic field profiles such as barriers, wells, channels, and so on, it is natural to approach this strain-engineering problem from the perspective of seeking which deformation fields applied to the carbon lattice lead to that prescribed PMF profile. As will be clear in subsequent sections, the solution is not unique. If not for anything else, this should be clear from the fact that there is a 'gauge' freedom in selecting the vector potential $\bf{\mathcal{A}}$ from ${\boldsymbol{B}} =\nabla \times {\bf{\mathcal{A}}}$. The simplest of such problems is to determine which displacement fields lead to a strictly uniform (space-independent) B. The first notable theoretical investigation along these lines was that of [26], who restricted their analysis to deformations in the plane. In this regime the PMFs are linear in the displacement field, allowing one to calculate an in-plane deformation field giving rise to any given PMF. In particular, to generate a constant or mostly constant PMF requires a characteristic deformation with three-fold symmetry (see figure 1), and the magnitude of the resulting PMF depends explicitly on the relative orientation of the deformation field and the underlying graphene lattice. This particular strain configuration has been recently explored in experiments on 'artificial graphene' [29].

Zoom In Zoom Out Reset image size

Extending deformations to three dimensions introduces nonlinearity into the strain field, and such simple solutions are no longer available. Continuum mechanical theoretical investigations in simple geometries such as one-dimensional bending [30, 31] or radially-symmetric bubbles [32], and atomistic simulations of graphene sheets adhered to nanoscale patterned substrates [33–35] are examples of forward problems: calculating the PMF associated with a certain deformation. But such approaches are unlikely to solve the inverse problem of finding the deformation mode required to produce a given PMF. In addition, given the current surge of experimental interest in deliberately inducing non-uniform strains in graphene, various possible routes are being explored [6, 20, 36, 37]. To be experimentally relevant, an attempt to effectively tackle the crucial inverse problem should be generic enough to encompass such diverse means to experimentally generate the desired strain fields.

This report presents a general-purpose framework which may be used to solve such inverse problems in graphene. In particular, for a given target PMF and experimental configuration, the method aims to find the optimal deformation control that, when applied to the graphene sheet, produces the desired PMF. Desired PMF refers to any specified space dependence of $B({\boldsymbol{X}} )$. Deformation control is the name for the geometric and mechanical parameters of the experiment that may be varied to change the deformation field. In the 2D example of [26] the deformation control is the displacement field applied to the outer boundary. In the case of graphene adhered to a patterned substrate the shape of the substrate performs that role. In this particular setup, which we will use extensively as an illustrative example in this report, the graphene sheet is assumed to have been transferred onto a patterned substrate, and forced to conform to its shape by combined hydrostatic pressure and adhesion forces (see figure 2). The aim in this case is to find the substrate pattern and pressure (the two deformation controls for this example) for which the deformed graphene sheet exhibits a desired target PMF. But we underline that the approach is straightforwardly applicable to any other target quantity with a known dependence on the strain field.

Zoom In Zoom Out Reset image size

We begin by summarizing the elastic properties of graphene and the elastic plate equations that govern its deformation when considered as a continuum elastic membrane. We then discuss the coupling of local deformations to the electronic degrees of freedom by means of the PMF and the optimization framework that forms the basis of our solution method. This is followed by a summary of the numerical algorithm used to solve the problem and, finally, as an example calculation, we present the computed substrate shapes that generate various predefined PMFs in an overlaid graphene sheet, and discuss the versatility of our framework for application in different experimental and theoretical scenarios well beyond the example calculations shown here for illustration. For completeness, various technical considerations and details are included as appendices to the main text.

2.1. Grapheneʼs elastic parametrization

The deformation of graphene is modeled using the equations from continuum elasticity. This formulation is chosen for its applicability across a wide range of lengthscales. In spite of recent developments based on discrete differential geometry to directly relate atomistic configurations with electronic properties of the type we envision [24], an atomistic approach to the elastic relaxation problem becomes easily unfeasible at scales of a few nanometers due to the intrinsically more numerically demanding nature of inverse problems.

The deformation of a graphene sheet is thus described in terms of its deviation from a flat 2D surface. The point ${\boldsymbol{X}} =(X,Y,0)$ is transformed to ${\boldsymbol{r}} =({{x}_{1}},{{x}_{2}},{{x}_{3}})=(x,y,z)$ in three-dimensional (3D) space, where ${{x}_{\alpha }}={{X}_{\alpha }}+{{v}_{\alpha }}({{X}_{1}},{{X}_{2}})$, and $z=w(X,Y)$. The deformation measures which describe the sheetʼs local stretching and bending are respectively the strain and curvature tensors ${{\varepsilon }_{\alpha \beta }}$ and ${{\rho }_{\alpha \beta }}$. Since these are complicated to write in terms of the displacement components, in practice simplified forms are used (for completeness, a detailed discussion of the form of ${{\varepsilon }_{\alpha \beta }}$ and ${{\rho }_{\alpha \beta }}$ is included in appendix A). The most common simplification is perhaps the von Kármán approximation, which uses the expressions

$$\begin{eqnarray}&&{{\varepsilon }_{\alpha \beta }}=\frac{1}{2}\left( \frac{\partial {{v}_{\alpha }}}{\partial {{X}_{\beta }}}+\frac{\partial {{v}_{\beta }}}{\partial {{X}_{\alpha }}}+\frac{\partial w}{\partial {{X}_{\alpha }}}\frac{\partial w}{\partial {{X}_{\beta }}} \right),\qquad {{\rho }_{\alpha \beta }}=\frac{{{\partial }^{2}}w}{\partial {{X}_{\alpha }}\partial {{X}_{\beta }}}.\end{eqnarray} \tag{ 1 }$$

The stress and moment resultants are assumed to be isotropic and linear in the strain and curvature tensors:

$$\begin{eqnarray}&&{{N}_{\alpha \beta }}=C\;{{\mathbb{A}}_{\alpha \beta \gamma \delta }}\ {{\varepsilon }_{\gamma \delta }},\qquad {{M}_{\alpha \beta }}=D\;{{\mathbb{A}}_{\alpha \beta \gamma \delta }}\ {{\rho }_{\gamma \delta }},\end{eqnarray} \tag{ 2 }$$

where we have used Einstein's summation convention, and

$$\begin{eqnarray}&&{{\mathbb{A}}_{\alpha \beta \gamma \delta }}=\left( \frac{1-\nu }{2} \right)({{\delta }_{\alpha \gamma }}{{\delta }_{\beta \delta }}+{{\delta }_{\alpha \delta }}{{\delta }_{\beta \gamma }})+\nu {{\delta }_{\alpha \beta }}{{\delta }_{\gamma \delta }},\end{eqnarray} \tag{ 3 }$$

ν being the Poisson ratio and C, D the stretching and bending moduli, respectively. To calculate the stretching modulus C, we use the results of [38] which, in our notation, are

$$\begin{eqnarray}&&C=358.1\;{\rm N}\;{{{\rm m}}^{-1}},\qquad C\nu =60.4\;{\rm N}\;{{{\rm m}}^{-1}}\qquad \Rightarrow \qquad \nu =0.169.\end{eqnarray} \tag{ 4 }$$

This value of ν agrees with the experimentally-measured Poisson ratio in graphite [39]. The value used for the bending modulus was that of [40], calculated ab initio as $D=1.46\;{\rm eV}=2.34\times {{10}^{-19}}\;{\rm N}\;{\rm m}$.

Note that the moduli C and D are independent in our formulation of the elastic response of graphene, which is treated as a purely 2D sheet. This means that thickness is not a parameter in our modeling. We emphasize this aspect because grapheneʼs elasticity is often modeled by treating it as a 3D material which is thin in one dimension, i.e. a conventional elastic thin plate. In those cases the stiffness and bending moduli are often written in terms of the 3D Youngʼs modulus E and the thickness h:

$$\begin{eqnarray}&&C=\frac{Eh}{1-{{\nu }^{2}}},\qquad D=\frac{E{{h}^{3}}}{12(1-{{\nu }^{2}})},\end{eqnarray} \tag{ 5 }$$

where a typical value $h\approx 0.3\;{\rm nm}$ for grapheneʼs effective thickness [41] is used. For example such expressions have been used to cite grapheneʼs Youngʼs modulus as of the order of $1\;{\rm TPa}$ [1, 42, 43]. While this may be useful to convey the scale and exceptional strength of graphene, the same numbers lead to an inaccurate value for the bending modulus D. Treating graphene as a continuous 3D elastic object might be a convenient approximation, but in keeping with grapheneʼs 2D nature, we retain the parameters C and D as our main quantities here rather than express them in terms of Youngʼs modulus E.

2.2. Equilibrium conditions

In addition to the kinematic and constitutive equations for a sheet of graphene, one must establish the equations of force balance to close the system. These are typically found by minimizing the potential energy functional consisting of two terms: ${{\mathcal{E}}^{{\rm elast}}}$, the stored elastic energy, and ${{\mathcal{E}}^{{\rm ext}}}$, the potential energy associated with external forces applied to the sheet. The latter may be surface tractions or adhesive forces (for simplicity we neglect any forces explicitly applied to the edge of the graphene sheet). The two energy terms are given by

$$\begin{eqnarray}&&{{\mathcal{E}}^{{\rm elast}}}[w,{{v}_{1}},{{v}_{2}}]={{\int }_{\Omega }}\left( \frac{1}{2}{{\varepsilon }_{\alpha \beta }}{{N}_{\alpha \beta }}+\frac{1}{2}{{\rho }_{\alpha \beta }}{{M}_{\alpha \beta }} \right)\;\;{{{\rm d}}^{2}}{\boldsymbol{X}} ,\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{{\mathcal{E}}^{{\rm ext}}}[w,{{v}_{1}},{{v}_{2}};{{\lambda }_{i}}]={{\int }_{\Omega }}V[w,{{v}_{1}},{{v}_{2}};{{\lambda }_{i}}]\;\;{{{\rm d}}^{2}}{\boldsymbol{X}} .\end{eqnarray} \tag{ 6b }$$

For now the energy density of external forces, $V[w,{{v}_{1}},{{v}_{2}};{{\lambda }_{i}}]$, is left unspecified. However, we do note that it is in this term that the influence of the control variables ${{\lambda }_{i}}$ is encoded; these may include, for instance, a parametrization of an underlying substrate, or the components of a surface traction field. The specific example corresponding to figure 2 will be presented in detail later.

In the standard variational formulation of the problem the potential energy is minimized by setting its first variation to zero, giving us three weak form equations for v₁, v₂, and w. However, this requires some regularity in the behavior of the transverse displacement w: its first derivative must be continuous (C¹). Choosing C¹ elements in an arbitrary triangular discretization is not trivial, however. To overcome this difficulty, we use a mixed variational principle [44], based on the work of Herrmann and Miyoshi [see 45–48], which involves treating the moment tensor ${{M}_{\alpha \beta }}$ as a separate variable. This allows us to treat the variables as continuous, and affine over each triangular element. For a graphene sheet with clamped conditions at the boundary, the six weak form equations that result are (see section A.2 for a detailed derivation):

$$\begin{eqnarray}&&{{\int }_{\Omega }}\left[ \frac{\partial {{\tilde{v}}_{1}}}{\partial X}{{N}_{11}}+\frac{\partial {{\tilde{v}}_{1}}}{\partial Y}{{N}_{12}}+{{\tilde{v}}_{1}}{{V}_{{{v}_{1}}}}[w,{{v}_{1}},{{v}_{2}};{{\lambda }_{i}}] \right]\;\;{{{\rm d}}^{2}}{\boldsymbol{X}} =0,\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{{\int }_{\Omega }}\left[ \frac{\partial {{\tilde{v}}_{2}}}{\partial X}{{N}_{12}}+\frac{\partial {{\tilde{v}}_{2}}}{\partial Y}{{N}_{22}}+{{\tilde{v}}_{2}}{{V}_{{{v}_{2}}}}[w,{{v}_{1}},{{v}_{2}};{{\lambda }_{i}}] \right]\;\;{{{\rm d}}^{2}}{\boldsymbol{X}} =0,\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&\begin{array}{l} {{\int }_{\Omega }}\left\{ \frac{\partial \tilde{w}}{\partial X}\left( -\frac{\partial {{M}_{11}}}{\partial X}-\frac{\partial {{M}_{12}}}{\partial Y}+{{N}_{11}}\frac{\partial w}{\partial X}+{{N}_{12}}\frac{\partial w}{\partial Y} \right) \right. \\ \qquad +\frac{\partial \tilde{w}}{\partial Y}\left( -\frac{\partial {{M}_{12}}}{\partial X}-\frac{\partial {{M}_{22}}}{\partial Y}+{{N}_{12}}\frac{\partial w}{\partial X}+{{N}_{22}}\frac{\partial w}{\partial Y} \right) \\ \qquad +\left. \tilde{w}{{V}_{w}}[w,{{v}_{1}},{{v}_{2}};{{\lambda }_{i}}] \right\}\;{{{\rm d}}^{2}}{\boldsymbol{X}} =0, \\ \end{array}\end{eqnarray} \tag{ 7c }$$

$$\begin{eqnarray}&&{{\int }_{\Omega }}\left[ \frac{1}{D(1-{{\nu }^{2}})}({{M}_{11}}-\nu {{M}_{22}}){{\tilde{M}}_{11}}+\frac{\partial w}{\partial X}\frac{\partial {{\tilde{M}}_{11}}}{\partial X} \right]\;{{{\rm d}}^{2}}{\boldsymbol{X}} =0,\end{eqnarray} \tag{ 7d }$$

$$\begin{eqnarray}&&{{\int }_{\Omega }}\left[ \frac{1}{D(1-\nu )}{{M}_{12}}{{\tilde{M}}_{12}}+\frac{1}{2}\frac{\partial w}{\partial X}\frac{\partial {{\tilde{M}}_{12}}}{\partial Y}+\frac{1}{2}\frac{\partial w}{\partial Y}\frac{\partial {{\tilde{M}}_{12}}}{\partial X} \right]\;{{{\rm d}}^{2}}{\boldsymbol{X}} =0,\end{eqnarray} \tag{ 7e }$$

$$\begin{eqnarray}&&{{\int }_{\Omega }}\left[ \frac{1}{D(1-{{\nu }^{2}})}({{M}_{22}}-\nu {{M}_{11}}){{\tilde{M}}_{22}}+\frac{\partial w}{\partial Y}\frac{\partial {{\tilde{M}}_{22}}}{\partial Y} \right]\;{{{\rm d}}^{2}}{\boldsymbol{X}} =0.\end{eqnarray} \tag{ 7f }$$

As will be clear shortly, these equilibrium equations provide the physical constraints to the optimization procedure. Our task is to seek a set of control parameters (substrate topography, boundary shape, etc) that, upon solution of the variational problem to find the equilibrium configuration of the elastic medium, yields a PMF distribution that best approaches the prescribed target. The implementation of this optimization is done numerically. We have chosen to use piecewise affine finite elements combined with a patch recovery method in our calculations for their simplicity and ease of implementation (see appendix E). But it should be noted that the method allows higher-order elements to be used, as long as one ensures that those formulations are stable and solvable.

2.3. Coupling deformations to electrons

To the weak-form elastic equilibrium equations we must add an equation linking the strain field to the generated PMF, $B({\boldsymbol{X}} )$. This is because we wish to find the deformation field that best approximates $B({\boldsymbol{X}} )$ to a desired target, say $\hat{B}({\boldsymbol{X}} )$. The origin of this PMF that appears in the low-energy effective Hamiltonian of deformed graphene is the local modification of the electronic hopping amplitudes, t, between neighboring carbon atoms brought about by the space dependent deformation of the crystal lattice. The hopping is constant in the perfect crystal: t₀ = 2.7 eV. But, since t depends strongly on the inter-atomic distance, any local change caused by a deformation leads to perturbations to this equilibrium value and, hence, more generically, $t\left( {{{\boldsymbol{X}} }_{i}},{{{\boldsymbol{X}} }_{i}}+{\boldsymbol{n}} \right)={{t}_{0}}+\delta t\left( {{{\boldsymbol{X}} }_{i}},{{{\boldsymbol{X}} }_{i}}+{\boldsymbol{n}} \right)$. The presence of $\delta t$, which is a relatively small perturbation to t₀ in practical situations, adds a correction to the low-energy Dirac-like Hamiltonian that emerges from a tight-binding description of the electronic hopping among p_z bands of adjacent carbon atoms. The effective Hamiltonian around the point ${\boldsymbol{K}} =(4\pi /3\sqrt{3}a,\;0)$ in the first Brillouin zone has the form [15, 16]

$$\begin{eqnarray}&&\mathcal{H}={{v}_{F}}\;{\boldsymbol{ \sigma }} \cdot \left( {\boldsymbol{p}} +e \right),{\bf{\mathcal{A}}}\end{eqnarray} \tag{ 8 }$$

where ${\boldsymbol{ \sigma }}$ is a vector of Pauli matrices, and $\hbar {{v}_{F}}=3{{t}_{0}}a/2$, with a = 1.42 Å the carbon–carbon distance in equilibrium. For deformations on scales that are large compared to a, the components of the pseudomagnetic vector potential ${\bf{\mathcal{A}}}={{\mathcal{A}}_{x}}{{{\boldsymbol{e}} }_{x}}+{{\mathcal{A}}_{y}}{{{\boldsymbol{e}} }_{y}}$ are explicitly given by [16]

$$\begin{eqnarray}&&{{\mathcal{A}}_{x}}({\boldsymbol{X}} )-{\rm i}{{\mathcal{A}}_{y}}({\boldsymbol{X}} )\simeq -\frac{\hbar c}{2ea}\left( {{\varepsilon }_{xx}}-{{\varepsilon }_{yy}}+2{\rm i}\;{{\varepsilon }_{xy}} \right),\end{eqnarray} \tag{ 9 }$$

where $c=-\partial {\rm log} t(r)/\partial {\rm log} r{{|}_{r=a}}$ (see appendix B). For static deformations, a value $c\approx 3.37$ captures the changes in various physical properties arising from strain-induced modifications of the π bands in agreement with first-principles calculations [49–54]. We note, however, that the effective low-energy Hamiltonian (8) contains only the leading order corrections arising from non-uniform deformations; further expanding in higher orders of smallness in the strain magnitude and the momentum with respect to ${\boldsymbol{K}}$ leads to terms that introduce, for example, Fermi surface anisotropy [54, 55] and space-dependent v_F [23, 24, 56, 57]. For simplicity, since we want to tailor only the PMF distribution as illustration of the method, and to keep the focus on the optimization framework rather than the details of the different levels of approximation for the effective strain-dependent Hamiltonian, we shall focus the subsequent analysis on the Hamiltonian (8). But it should be clear that, as far as the optimization procedure is concerned (which does not take into account the energy of the electronic system), the particular form of $\mathcal{H}$ is only relevant in order to identify the target quantity that we wish to optimize and its expression in terms of the strain components, as in (9). If instead of the PMF we were interested in optimizing, for example, towards a desired space-dependence of the Fermi velocity [23] or that of the deformation potential [16], the method requires only the specification of its functional dependence on strain.

Finally, the PMF itself, B, being defined as the 2D curl of $\bf{\mathcal{A}}$, reads:

$$\begin{eqnarray}&&B({\boldsymbol{X}} )=\frac{\hbar c}{ae}\left[ \frac{\partial }{\partial Y}\left( \frac{{{\varepsilon }_{11}}-{{\varepsilon }_{22}}}{2} \right)+\frac{\partial {{\varepsilon }_{12}}}{\partial X} \right].\end{eqnarray} \tag{ 10 }$$

As noted above, by virtue of our choice of piecewise affine finite elements for the numerical interpolation, the six variables v₁, v₂, w, M₁₁, M₁₂, and M₂₂ are treated as continuous, and affine over each triangular element. As a consequence, the strain field will be discontinuous and constant in each triangular element, leaving the PMF (10) undefined within this interpolation scheme. To overcome this we use the technique of patch recovery [58] detailed in appendix E. In brief, this is a mechanism that uses the discontinuous strain data ${{\varepsilon }_{\alpha \beta }}$ to recover a strain field $\varepsilon _{\alpha \beta }^{{\rm rec}}$ of the same type as the primary variables: continuous and affine over each element. The derivative of $\varepsilon _{\alpha \beta }^{{\rm rec}}$ is well-defined, and thus so is the PMF if it is calculated using this recovered strain field.

Solution of the weak form variational (equilibrium) equations constitutes the forward problem: in other words, given a set of control variables (here chosen to be substrate shape, encoded in the external potential V), what deformation and PMF do these conditions impose on the graphene sheet? This report is aimed at answering the corresponding inverse problem: what are the control variables that will give rise to a desired PMF?

This inverse question is posed as an optimization problem, where an integral is minimized subject to the weak form equations written explicitly in equations (7a )–(7f ). If we let $\hat{B}(X,Y)$ be the desired PMF in Lagrangian coordinates, we then seek to minimize the functional

$$\begin{eqnarray}&&\mathcal{I}[w,{{v}_{1}},{{v}_{2}}]={{\int }_{\Omega }}{{\left( B[w,{{v}_{1}},{{v}_{2}}]-\hat{B}(X,Y) \right)}^{2}}\;{{{\rm d}}^{2}}{\boldsymbol{X}} ,\end{eqnarray} \tag{ 11 }$$

to find a PMF, B, which is (ideally everywhere) as close as possible to the prescribed $\hat{B}(X,Y)$ (the reader will note once more at this stage that the quantity $B[w,{{v}_{1}},{{v}_{2}}]$, which is here associated with the PMF, can be replaced by any other of interest, as long as its dependence on the strain or deformation field can be specified; the scope of applicability of this method extends, therefore, well beyond the PMF example chosen here for definiteness).

This sort of optimization problem, however, is typically mathematically ill-posed, in the sense that there are infinitely-many solutions to such a minimization and, in order to find a solution which also satisfies the weak form equations, the numerical method often yields a solution which is not smooth. To counter this phenomenon, one must add to the minimization integral, $\mathcal{I}$, a regularization term which penalizes high spatial variations in the control variables ${{\lambda }_{i}}$:

$$\begin{eqnarray}&&\mathcal{I}[w,{{v}_{1}},{{v}_{2}}]={{\int }_{\Omega }}{{\left( B[w,{{v}_{1}},{{v}_{2}}]-\hat{B}(X,Y) \right)}^{2}}\;{{{\rm d}}^{2}}{\boldsymbol{X}} +\eta \;{{\mathcal{I}}^{{\rm reg}}}[{{\lambda }_{i}}],\end{eqnarray} \tag{ 12 }$$

where η is a tunable parameter. The precise form of ${{\mathcal{I}}^{{\rm reg}}}[{{\lambda }_{i}}]$ will depend on what the control variables ${{\lambda }_{i}}$ represent; for the specific example of substrate shape optimization a typical form will be discussed below. Thus the full problem is to minimize, by varying the six state variables v₁, v₂, w, M₁₁, M₁₂, M₂₂ and control variables ${{\lambda }_{i}}$, the objective function (12) subject to the six equations (7a )–(7f ), solved for all admissible variations $(\tilde{\cdot })$. In these expressions, B is given by (10), ${{N}_{\alpha \beta }}$ is given by (2), and ${{\varepsilon }_{\alpha \beta }}$ is given by (1). This problem is an example of a PDE-constrained optimization. For technical details regarding well-posedness and solution methods for such problems the reader is referred to [59] or Borzì and Schulz [60]. This general procedure has also been applied to shape optimization in elastic plates experiencing differential growth fields [61].

Finally, in solving this problem numerically, all equations are adimensionalized in such a way that most variables are $O(1)$ to ensure good numerical behavior (details described in appendix C).

We wish to apply the previously-developed general theory to a specific example, to illustrate its practical implementation, and its utility in the problem at hand. The example we have in mind is of a graphene sheet forced to conform to a certain substrate shape by pressure and adhesive forces, as depicted in figure 2. The elevation of the substrate is denoted $z=\hat{z}(x,y)$ (we use coordinates x, y as Eulerian coordinates rather than the Lagrangian X, Y used in the definition of the graphene deformation). For definiteness, we assume that we are searching for target PMFs $\hat{B}$ of typical scale ${{B}_{0}}=10\;{\rm T}$. We further assume that the domain Ω, representing the shape of the computational domain, has typical dimensions $L=100\;\overset{\circ}{A}$.

4.1. The external forces

As discussed above, the external potential term $V[{{v}_{1}},{{v}_{2}},w;{{\lambda }_{i}}]$ will have two components: the work done due to hydrostatic pressure, $pw$, and the adhesion energy ${{V}_{{\rm adh}}}$ between the graphene flake and the substrate. To find the adhesion energy, consider the graphene sheet as a collection of atoms, interacting with a field ${{V}_{p}}(x,y,z)$ in ${{\mathbb{R}}^{3}}$-space, such that a particle $\;{\rm d}S$ of the plate, located at $(X+{{v}_{1}},Y+{{v}_{2}},w)$ in Eulerian coordinates, contributes ${{V}_{p}}(X+{{v}_{1}},Y+{{v}_{2}},w)\;{\rm d}S$ to the adhesion energy. Then ${{V}_{{\rm adh}}}[w,{{v}_{1}},{{v}_{2}};{{\lambda }_{i}}]={{V}_{p}}(X+{{v}_{1}},Y+{{v}_{2}},w)$.

Given a substrate shape $z=\hat{z}(x,y)$, we could determine the adhesion potential V_p at every point in 3D space. But this will be time-intensive in general, and for optimization problems could prove prohibitively expensive. As an alternative, assume that the gradient of the substrate is small, so that we can approximate ${{V}_{p}}(x,y,z)=J(z-\hat{z}(x,y))$, where J(s) is some one-dimensional adhesion potential, such as the Lennard-Jones potential between surfaces [62],

$$\begin{eqnarray}&&J(s)=\frac{{{J}_{0}}}{3}\left[ {{\left( \frac{{{s}^{*}}}{s} \right)}^{8}}-4{{\left( \frac{{{s}^{*}}}{s} \right)}^{2}} \right],\end{eqnarray} \tag{ 13 }$$

where s* is the adhesion well position, or the distance from the substrate at which a particle is in equilibrium. Thus,

$$\begin{eqnarray}&&V[w,{{v}_{1}},{{v}_{2}};{{\lambda }_{i}}]=pw+J\left( w-\hat{z}(X+{{v}_{1}},Y+{{v}_{2}}) \right),\end{eqnarray} \tag{ 14 }$$

with the control variables ${{\lambda }_{i}}$ encoded in $\hat{z}$. As a representative value for s*, we use the value that [63] give for C–Cu, namely $2.243\;\overset{\circ}{A}$. Similarly we use ${{J}_{0}}=0.45\;{\rm J}\;{{{\rm m}}^{-2}}$ as a representative value, from the investigation of [7] into the adhesion strength between graphene and SiO$_{2}$. We select a typical value for the hydrostatic pressure as $p=100\;{\rm bar}={{10}^{7}}\;{\rm Pa}$.

4.2. Parametrization of the substrate topography

In writing an expression for the substrate geometry, a naïve approach would be to use the same type of discretization as for the graphene sheet itself. The domain Ω is discretized into a collection of triangles, and each of the state variables (w, ${{v}_{\alpha }}$, ${{M}_{\alpha \beta }}$) is posited to be continuous and affine over each triangle. This allows each variable to be described entirely in terms of its values at the nodes of the triangulation, and is the main advantage in writing the state equations in weak form³ . However, this approach will not work when it comes to describing the shape of the substrate, $z=\hat{z}(x,y)$. Since the graphene sheet can move in a lateral direction, the triangulations of the substrate and the sheet itself will not remain in registration. Therefore, for a given nodal point the distance measure $w-\hat{z}(X+\varepsilon {{v}_{1}},Y+\varepsilon {{v}_{2}})$ will not vary smoothly as ${{v}_{\alpha }}$ are varied. The alternative, which we will follow here, is to construct a smooth shape for the substrate.

For the numerical experiments in this article, we assume that the substrate is patterned periodically in the two horizontal directions, and set the repeating 2D unit cell to be a rhombus. If we introduce the two coordinates ${{\xi }_{1}}=y/\sqrt{3}+x$ and ${{\xi }_{2}}=y/\sqrt{3}-x$ the unit cell corresponds to the domain

$$\begin{eqnarray}&&\Omega =\left\{ ({{\xi }_{1}},{{\xi }_{2}})\ :\ 0\leqslant {{\xi }_{1}}\lt 1,\ 0\leqslant {{\xi }_{2}}\lt 1 \right\}.\end{eqnarray} \tag{ 15 }$$

The topography of the substrate, $\hat{z}(x,y)$, can then be resolved into a truncated Fourier expansion with the period of the unit cell Ω. Since the (finite) set of expansion coefficients determines the overall topography, they play the role of the control variables ${{\lambda }_{i}}$: varying the topography of the substrate is, therefore, achieved by varying these expansion constants (refer to appendix D for particulars of this approach). The periodicity of the substrate places limits on the patterns of PMF that can be sought. If we integrate the PMF over the unit rhombus, we find

$$\begin{eqnarray}&&{{\int }_{\Omega }}B\;\;{{{\rm d}}^{2}}{\boldsymbol{X}} ={{\oint }_{\partial \Omega }}\left[ \left( \frac{{{\varepsilon }_{11}}-{{\varepsilon }_{22}}}{2} \right){{n}_{2}}+{{\varepsilon }_{12}}{{n}_{1}} \right]\;{\rm d}s=0,\end{eqnarray} \tag{ 16 }$$

because the strain fields generated by the periodic substrate will also be periodic. Thus (if we limit ourselves to periodic substrates) it is impossible to generate PMFs whose integral over the unit cell is non-zero—in particular this rules out the generation of strictly constant non-zero PMFs by periodic deformations. It should be emphasized, however, that our general procedure is applicable to arbitrary domains, geometries, and target PMFs, which require different parametric expansions of the substrate shape in place of the Fourier series employed here. The periodic choice is used by us simply out of convenience, precisely for its straightforward Fourier expansion that allows the description of an arbitrarily patterned substrate using mathematically simple trigonometric functions.

As discussed earlier, in order to avoid convergence towards solutions that are ill-behaved during the numerical optimization, a regularization term, ${{\mathcal{I}}^{{\rm reg}}}$, is added to the objective integral, as per equation (12) (in this case such ill-behaved solutions could be, for example, substrate profiles with discontinuities or sharp topographical features). We choose it to be

$$\begin{eqnarray}&&{{\mathcal{I}}^{{\rm reg}}}=\frac{1}{{\rm Area}(\Omega )}{{\int }_{\Omega }}{{\left| \nabla \hat{z} \right|}^{2}}{{{\rm d}}^{2}}{\boldsymbol{X}} ,\end{eqnarray} \tag{ 17 }$$

which is simple to calculate using the orthogonality of the basis functions in the Fourier expansion over Ω. This expression measures the fineness of spatial variation in the substrate: ${{\mathcal{I}}^{{\rm reg}}}$ is larger for substrate profiles with smaller wavelengths. The minimization of (12) thus leads to the penalizing of such rough profiles that would be unrealistic in view of the finite feature resolution associated with any experimental approach to substrate patterning. Moreover, if a particular experimental implementation is to be carried out, the regularization term can be further refined or adapted to reflect the specific geometric, fabrication, or other constraints.

For illustration we initially seek PMFs with a typical strength of $10\;{\rm T}$. The unit cell is chosen to have edges of length $L={{10}^{-8}}\;{\rm m}$ and we apply a pressure of $100\;{\rm bar}$. These parameters might seem at the threshold of current experimental applicability, but they are chosen for their numerical tractability—since if s* is too small, standard numerical algorithms will iterate over trial configurations with negative graphene–substrate separations s, a highly non-smooth problem. To overcome this issue, one would need to carefully design algorithms in which negative separations were avoided. But, for the proof-of-principle calculations reported here, we choose this acceptable compromise.

We obtain good results for the numerical parameters K = 2, $\eta ={{10}^{-9}}$, and for an isometric mesh of 800 triangles in the unit cell. For illustration in this report we have chosen the following four target PMF patterns to optimize for:

• (i)  
  An alternating pattern of equilateral triangles,
• (ii)  
  A kagome-like pattern of zero-valued hexagons interspersed with alternating triangular regions,
• (iii)  
  Alternating strips of positive and negative fields,
• (iv)  
  A rhombus of negative PMF surrounded by a field of positive PMF.

Each of these target fields (shown in the first column of figure 3) adheres to the condition that its integral over the unit cell must be zero. In the second column we see the PMF attained by the optimization code, and in the third column we see the substrate topography that produced such a field. It is clear that the converged solutions reproduce with very good accuracy the spatial dependence of the target PMF in all four cases, including the rapid sign changes imposed by the target field which, being necessarily smooth in the solution, are still quite sharp, with the sign change occurring in a very short length scale. This documents how this optimization strategy is able to capture all the features, global and detailed, of the target PMF. From an experimental point of view, the power of this method is clear: by providing an accurate solution to the inverse problem, it allows one to specify in detail what substrate pattern and topography leads to a PMF of given magnitude and space dependence. This, in principle, provides all the experimental information needed to fabricate the corresponding structures.

Zoom In Zoom Out Reset image size

It is important to recall that, as pointed out earlier, the solution to the pseudomagnetic inverse problem might not be unique. This is because the objective function (12) is defined in terms of the field B, whose relation to the deformation field expressed in equation (10) allows for a large 'gauge' freedom. As a result, more than one set of control parameters within a certain parameter range might be simultaneously compatible with the target field within a desired accuracy and obey the elastic plate equilibrium equations. On the other hand, having found a set of parameters that optimizes the induced PMF to the target sought is not a guarantee that such set will remain optimal upon finite changes of an external variable, or a constant scaling of the target function. This latter aspect is best illustrated with a specific case for our example system. The panels in the last column of figure 3 show the suitable substrate profiles for a pressure of 100 bar and PMFs with an amplitude of 10 T. The equivalent calculations for a target PMF magnitude of 100 T, a ten-fold increase, yield the results shown in figure 4. It is clear that the optimal substrate topographies that guarantee the same degree of proximity between the induced and target PMF as before are markedly different from those in figure 3. This means that in some experimental setups, such as the one sketched in figure 2, the proximity of the induced PMF to the target might need to be compromised in favor of a having a fixed set of control parameters suitable in a range of PMF field amplitudes (i.e., a single substrate profile able to generate acceptable PMF profiles of different amplitude). However, the power of the method and its experimental practicality in allowing a direct, one-step route from PMF design to substrate fabrication, should largely compensate such compromises, when they are unavoidable.

Zoom In Zoom Out Reset image size

There are many ways in which one can experimentally control the deformation of a graphene sheet [64], with each choice leading to a different set of control variables. Most obviously, the graphene sheet may be manipulated directly, whether by substrate topology (such as nanopillars) [6, 33], a distribution of attached structures like nanotubes [27], or a nano-manipulation of substrate adhesion properties [20]. The corresponding control variables would be the configurations of the nanostructures, including their shapes, their positions in relation to the graphene sheet, and their height. Edge actuation, where the control variables are the displacements applied to the edge, is another deformation mode [26, 30]. In experiments, it may be preferable to apply these edge displacements indirectly, by applying electromechanical forces to the electrodes attached to a graphene flake [65]. The position and shape of the electrodes form the control variables in this case.

A further class of deformations in graphene are the inflation of bubbles by suspending graphene over a particularly-shaped cavity, whether by hydrostatic pressure or electromechanical forces [2, 9]. Since the forcing in these examples is global, it is the shape of the cavity that provides the variation in the calculated strain field, and as such the control variables in an optimization calculation would be a parametrization of the cavity shape. Such inflation problems can be coupled with local deformation, in the form of a point deformation due to a STM tip [2, 3], providing additional control variables of tip position and strength and allowing a greater ability to achieve desired strain fields and consequential electronic properties.

Finally, it should be noted that, despite our focus in this report on optimizing the control parameters for a target PMF (which constitutes the core of the strain-engineering concept in graphene), this optimization framework can be rather easily extended to other target quantities by replacing the objective function in equation (12) by the relevant measure of 'distance' for that problem, and specifying its dependence on the strain or displacement field analogously to the specification in equation (10). Each of the three main components of the procedure—namely the objective function that is minimized, the state equations that form the constraints (here being the elastic plate equations), and the control variables that are open to experimental variation—can be changed to answer different questions of interest.

As simple examples, we suggest that one may wish to minimize or maximize the degree of rippling obtained in the edge actuation of a suspended graphene sheet [8], or that the resonant frequency of a graphene flake suspended over a cavity [66, 67] may be optimized by varying the cavity shape according to the principles outlined in this exposition.

We presented a general-purpose framework suitable to answer the following inverse problem in graphene: which set of external control parameters (substrate topography, sample shape, load distribution, etc) guarantees that the resulting equilibrium state of graphene exhibits a PMF that varies in space in a prescribed way? The ability to answer this question in general, given only a potential experimental setup and the target field profile, is paramount towards fulfilling the vision of tailored transport and other electronic properties in graphene by strain-engineering. This concept calls for expedited ways to answer the question above. The method presented here relies on a PDE-constrained optimization strategy to minimize the generic objective function (12) which penalizes significant deviations between the induced and target PMFs. It thus affords a one-step route from PMF design to experimental implementation, is unbiased and general enough to accommodate a multitude of experimental parameters and conditions that can be envisaged to produce the desired deformations in the graphene lattice, and always ensures compliance with the constraints imposed by elasticity theory and the equilibrium conditions of graphene treated as a continuous elastic medium. We trust that it can be an important tool in designing or guiding experimentally realistic conditions for strain-engineered graphene devices and beyond—the versatility to define, in principle, any target function for other physical quantities entails a broad applicability.

The authors acknowledge support from the NRF-CRP grant 'Novel 2D materials with tailored properties: beyond graphene' (R-144-000-295-281), as well as productive discussions with A H Castro Neto and L Mahadevan.

Here we recapitulate the results of sections 2.1 and 2.2 with further discussion on their validity and applicability.

A.1. Definitions and assumptions

The deformation of a graphene sheet is described in terms of its deviation from a flat 2D surface. A point ${\boldsymbol{X}}$ in an undeformed flat surface is defined by its coordinates $({{X}_{1}},{{X}_{2}},0)=(X,Y,0)$, for all (X, Y) belonging to some set Ω that defines the physical domain. Under a deformation the point ${\boldsymbol{X}} =(X,Y,0)$ is transformed to ${\boldsymbol{r}} =({{x}_{1}},{{x}_{2}},{{x}_{3}})=(x,y,z)$ in 3D space, where

$$\begin{eqnarray}\begin{array}{rcl} x & = & X+{{v}_{1}}(X,Y), \\ y & = & Y+{{v}_{2}}(X,Y), \\ z & = & w(X,Y), \\ \end{array}\end{eqnarray} \tag{ A.1 }$$

are the expressions for the new coordinates in terms of the in-plane displacements ${{v}_{\alpha }}$, and the vertical deflection w.⁴ From these expressions one can define the base vectors ${{{\boldsymbol{r}} }_{,\alpha }}=({{x}_{,\alpha }},{{y}_{,\alpha }},{{z}_{,\alpha }})$ and the metric tensor ${{g}_{\alpha \beta }}={{{\boldsymbol{r}} }_{,\alpha }}\cdot {{{\boldsymbol{r}} }_{,\beta }}$ of the deformed surface. Then the true strain tensor is defined in terms of the difference of this metric tensor from its original value ${{\delta }_{\alpha \beta }}$ (Kroneckerʼs delta, or the identity tensor, from our choice of Cartesian coordinates):

$$\begin{eqnarray}&&\varepsilon _{\alpha \beta }^{{\rm true}}=\frac{1}{2}({{g}_{\alpha \beta }}-{{\delta }_{\alpha \beta }}).\end{eqnarray} \tag{ A.2 }$$

In terms of displacements, this becomes

$$\begin{eqnarray}&&\varepsilon _{\alpha \beta }^{{\rm true}}=\frac{1}{2}\left( {{v}_{\alpha ,\beta }}+{{v}_{\beta ,\alpha }}+{{w}_{,\alpha }}{{w}_{,\beta }}+{{v}_{\gamma ,\alpha }}{{v}_{\gamma ,\beta }} \right).\end{eqnarray} \tag{ A.3 }$$

The second deformation measure of the surface is the curvature tensor, defined by ${{\rho }_{\alpha \beta }}={{{\boldsymbol{r}} }_{,\alpha \beta }}\cdot {\boldsymbol{n}}$, where ${\boldsymbol{n}}$ is the unit normal vector to the deformed surface. We do not derive the result here (for details in the case of a curved elastic shell see [68] or [69]) but the full expression for the curvature tensor in Cartesian coordinates is

$$\begin{eqnarray}&&\begin{array}{l} \rho _{\alpha \beta }^{{\rm true}}={{g}^{-1/2}}{{w}_{,\alpha \beta }}\left( 1+{{v}_{\lambda ,\lambda }}+\frac{1}{2}{{v}_{\lambda ,\lambda }}{{v}_{\mu ,\mu }}-\frac{1}{2}{{v}_{\lambda ,\mu }}{{v}_{\mu ,\lambda }} \right) \\ -{{g}^{-1/2}}{{v}_{\kappa ,\alpha \beta }}{{w}_{,\lambda }}\left( {{\delta }_{\kappa \lambda }}+{{\delta }_{\kappa \lambda }}{{v}_{\mu ,\mu }}-{{v}_{\lambda ,\kappa }} \right), \\ \end{array}\end{eqnarray} \tag{ A.4 }$$

where $g={\rm det} \;{{g}_{\alpha \beta }}$.

Expressions (A.3) and (A.4) are far too unwieldy for most purposes. Based on assumptions regarding the relative sizes of the displacement components and the length scale of deformations, the strain and curvature tensors are simplified. We chose the von Kármán approximation for its simplicity and its capacity to model moderate deflections. This simplification uses the expressions

$$\begin{eqnarray}&&{{\varepsilon }_{\alpha \beta }}=\frac{1}{2}\left( {{v}_{\alpha ,\beta }}+{{v}_{\beta ,\alpha }}+{{w}_{,\alpha }}{{w}_{,\beta }} \right),\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&{{\rho }_{\alpha \beta }}={{w}_{,\alpha \beta }}.\end{eqnarray} \tag{ A.6 }$$

Denote the corrections to these approximations by $\varepsilon _{\alpha \beta }^{{\rm corr}}=\varepsilon _{\alpha \beta }^{{\rm true}}-{{\varepsilon }_{\alpha \beta }}$, and similarly for $\rho _{\alpha \beta }^{{\rm corr}}$. As an a posteriori check on the validity of our solutions, we can verify that the approximations are close to the true values, or

$$\begin{eqnarray}&&\varepsilon _{\alpha \beta }^{{\rm corr}}\ll {{\varepsilon }_{\alpha \beta }},\qquad \rho _{\alpha \beta }^{{\rm corr}}\ll {{\rho }_{\alpha \beta }}.\end{eqnarray} \tag{ A.7 }$$

The stress and moment resultants are assumed to be isotropic and linear in the strain and curvature tensors:

$$\begin{eqnarray}&&{{N}_{\alpha \beta }}=C\;{{\mathbb{A}}_{\alpha \beta \gamma \delta }}\ {{\varepsilon }_{\gamma \delta }},\qquad {{M}_{\alpha \beta }}=D\;{{\mathbb{B}}_{\alpha \beta \gamma \delta }}\ {{\rho }_{\gamma \delta }},\end{eqnarray} \tag{ A.8 }$$

where

$$\begin{eqnarray}&&{{\mathbb{A}}_{\alpha \beta \gamma \delta }}=\left( \frac{1-\nu }{2} \right)({{\delta }_{\alpha \gamma }}{{\delta }_{\beta \delta }}+{{\delta }_{\alpha \delta }}{{\delta }_{\beta \gamma }})+\nu {{\delta }_{\alpha \beta }}{{\delta }_{\gamma \delta }},\end{eqnarray} \tag{ A.9 }$$

$$\begin{eqnarray}&&{{\mathbb{B}}_{\alpha \beta \gamma \delta }}=\left( \frac{1-\sigma }{2} \right)({{\delta }_{\alpha \gamma }}{{\delta }_{\beta \delta }}+{{\delta }_{\alpha \delta }}{{\delta }_{\beta \gamma }})+\sigma {{\delta }_{\alpha \beta }}{{\delta }_{\gamma \delta }},\end{eqnarray} \tag{ A.10 }$$

ν being the Poisson ratio and C, D the stretching and bending moduli, respectively. We have defined σ to be the analog of the Poisson ratio for bending deformations; if D_G is the Gaussian bending rigidity in the Helfrich free energy for the bending of a membrane [70], then $\sigma =1+{{D}_{G}}/D$. To calculate the stretching modulus C, we use the results of [38], who fitted a polynomial stress–strain relation to ab initio calculations up to strains of 50%. For simplicity, we will assume a linear stress–strain relationship, which is valid only up to strains of around 10%. The linear terms of [38] are, in our notation,

$$\begin{eqnarray}\left. \begin{array}{ccccccccccccccc} C & = & 358.1\;{\rm N}\;{{{\rm m}}^{-1}} \\ C\nu & = & 60.4\;{\rm N}\;{{{\rm m}}^{-1}} \\ \end{array} \right\}\quad \Rightarrow \quad \nu =0.169.\end{eqnarray} \tag{ A.11 }$$

This value of ν agrees with the experimentally-measured Poisson ratio in graphite [39]. The value we chose for the bending modulus was that of [40], calculated ab initio as $D=1.46\;{\rm eV}=2.34\times {{10}^{-19}}\;{\rm N}\;{\rm m}$. We have found only two investigations into the value of D_G (and hence σ) in graphene; the calculations of [71] lead to $\sigma =-0.056$, whereas the numerical study of [72] gives a value of $\sigma =0.565$. In the absence of consensus, in our calculations σ is set to be equal to the Poisson ratio: $\sigma =\nu =0.169$, and thus ${{\mathbb{B}}_{\alpha \beta \gamma \delta }}={{\mathbb{A}}_{\alpha \beta \gamma \delta }}$.

The constitutive equations for macroscopic materials are usually derived from full 3D isotropic elasticity expressions in the limit that the plate thickness is small. In the most rigorous treatments this analysis leads to limits on the validity of simplifying expressions such as (A.6) in terms of the relative sizes of stored elastic energy, applied surface tractions, and plate thickness [73–75]. Using such analyses, $\sigma =\nu$ and the stiffness and bending moduli may be written in terms of the 3D Youngʼs modulus E and the thickness h:

$$\begin{eqnarray}&&C=\frac{Eh}{1-{{\nu }^{2}}},\qquad D=\frac{E{{h}^{3}}}{12(1-{{\nu }^{2}})}.\end{eqnarray} \tag{ A.12 }$$

Using a typical value $h\approx 0.3\;{\rm nm}$ for graphene thickness [41], this has been used to cite grapheneʼs Youngʼs modulus as of the order of $1\;{\rm TPa}$ [1, 42, 43]. While this may be useful to convey the scale and exceptional strength of graphene, the same numbers lead to an inaccurate value for the bending modulus D. Treating graphene as a continuous 3D elastic object is a convenient approximation, so for definiteness we keep the 2D parameters C and D as our main quantities here rather than express them in terms of Youngʼs modulus E. A rigorous justification of the plate equations used to model graphene deformations is beyond the scope of this paper, and would involve a detailed analysis of the stored energy involved together with applied surface tractions. For the purposes of this paper it is enough to ensure that the strains and curvatures are within reasonable limits ($|{{\varepsilon }_{\alpha \beta }}|\lt 0.1$, $|{{\rho }_{\alpha \beta }}|\lt {{h}^{-1}}$) and that the corrections to the strains and curvatures are small.

A.2. Weak form equations

In section 2.2 we stated that equations (6a )–(6f ) could be derived from a minimization of the energy integrals (5a )–(5b ). In this section we justify this claim.

Recall that the stored energy in the plate was given by

$$\begin{eqnarray}&&\mathcal{E}={{\int }_{\Omega }}\left( \frac{1}{2}{{\varepsilon }_{\alpha \beta }}{{N}_{\alpha \beta }}+\frac{1}{2}{{\rho }_{\alpha \beta }}{{M}_{\alpha \beta }}+V[w,{{v}_{1}},{{v}_{2}};{{\lambda }_{i}}] \right)\;{{{\rm d}}^{2}}{\boldsymbol{X}} .\end{eqnarray} \tag{ A.13 }$$

To help understand the mixed variational principles that we rely on, let us consider a simplified problem of purely transverse deflections of an elastic plate subject to hydrostatic pressure:

$$\begin{eqnarray}&&{{\mathcal{E}}^{{\rm T}}}={{\int }_{\Omega }}\left( \frac{1}{2}{{\rho }_{\alpha \beta }}{{M}_{\alpha \beta }}+pw \right)\;{{{\rm d}}^{2}}{\boldsymbol{X}} ={{\int }_{\Omega }}\left( \frac{D}{2}{{\mathbb{B}}_{\alpha \beta \gamma \delta }}{{w}_{,\alpha \beta }}{{w}_{,\gamma \delta }}+pw \right)\;{{{\rm d}}^{2}}{\boldsymbol{X}} .\end{eqnarray} \tag{ A.14 }$$

Assume that the boundary of this plate is formed of three disjoint regions: $\partial \Omega ={{\Gamma }_{c}}\cup {{\Gamma }_{s}}\cup {{\Gamma }_{f}}$, with clamped conditions along Γ_c, simply supported conditions along Γ_s, and free conditions along Γ_f.

The standard variational approach is to minimize ${{\mathcal{E}}^{{\rm T}}}$ over the space of all twice-differentiable w satisfying w = 0 on ${{\Gamma }_{c}}\cup {{\Gamma }_{s}}$ and ${{\partial }_{n}}w=0$ on Γ_c. The first variation of ${{\mathcal{E}}^{{\rm T}}}$ is

$$\begin{eqnarray}&&\delta {{\mathcal{E}}^{{\rm T}}}={{\int }_{\Omega }}\left( D{{\mathbb{B}}_{\alpha \beta \gamma \delta }}{{w}_{,\alpha \beta }}{{\tilde{w}}_{,\gamma \delta }}+p\tilde{w} \right)\;{{{\rm d}}^{2}}{\boldsymbol{X}} ,\end{eqnarray} \tag{ A.15 }$$

where $\tilde{w}=\delta w$ is the variation in w. The weak solution is then the twice-differentiable function $w(X,Y)$ that satisfies $\delta {{\mathcal{E}}^{{\rm T}}}=0$ for each twice differentiable variation $\tilde{w}$ satisfying $\tilde{w}=0$ on ${{\Gamma }_{c}}\cup {{\Gamma }_{s}}$ and ${{\partial }_{n}}\tilde{w}=0$ on Γ_c. To find the strong form equation and boundary conditions to which this weak formulation corresponds, assume that w is four-times differentiable and integrate (A.15) twice by parts:

$$\begin{eqnarray}\begin{array}{rcl} \delta {{\mathcal{E}}^{{\rm T}}} & = & {{\int }_{\Omega }}\tilde{w}\left( D{{\mathbb{B}}_{\alpha \beta \gamma \delta }}{{w}_{,\alpha \beta \gamma \delta }}+p \right){{{\rm d}}^{2}}{\boldsymbol{X}} +{{\oint }_{\partial \Omega }}\left( {{\tilde{w}}_{,\alpha }}{{M}_{\alpha \beta }}{{n}_{\beta }}-\tilde{w}{{M}_{\alpha \beta ,\alpha }}{{n}_{\beta }} \right){\rm d}s \\ {} & = & {{\int }_{\Omega }}\tilde{w}\left( D{{\mathbb{B}}_{\alpha \beta \gamma \delta }}{{w}_{,\alpha \beta \gamma \delta }}+p \right){{{\rm d}}^{2}}{\boldsymbol{X}} +{{\oint }_{\partial \Omega }}\left[ ({{\partial }_{n}}\tilde{w}{{n}_{\alpha }}+{{\partial }_{t}}\tilde{w}{{t}_{\alpha }}){{M}_{\alpha \beta }}{{n}_{\beta }}-\tilde{w}{{M}_{\alpha \beta ,\alpha }}{{n}_{\beta }} \right]{\rm d}s \\ {} & = & {{\int }_{\Omega }}\tilde{w}\left( D{{\mathbb{B}}_{\alpha \beta \gamma \delta }}{{w}_{,\alpha \beta \gamma \delta }}+p \right){{{\rm d}}^{2}}{\boldsymbol{X}} +{{\int }_{{{\Gamma }_{s}}\cup {{\Gamma }_{f}}}}\frac{\partial \tilde{w}}{\partial n}{{M}_{\alpha \beta }}{{n}_{\alpha }}{{n}_{\beta }}\;{\rm d}s \\ {} & {} & -{{\int }_{{{\Gamma }_{f}}}}\tilde{w}\left[ \frac{\partial }{\partial s}({{M}_{\alpha \beta }}{{t}_{\alpha }}{{n}_{\beta }})+{{M}_{\alpha \beta ,\alpha }}{{n}_{\beta }} \right]\;{\rm d}s, \\ \end{array}\end{eqnarray} \tag{ A.16 }$$

using the boundary conditions for $\tilde{w}$. Setting this to zero for each admissible variation $\tilde{w}$, we find that the governing equation is

$$\begin{eqnarray}&&D{{\nabla }^{4}}w+p=0,\end{eqnarray} \tag{ A.17 }$$

with boundary conditions

$$\begin{eqnarray}&&w=\frac{\partial w}{\partial n}=0\qquad {\rm on}\;{{\Gamma }_{c}},\end{eqnarray} \tag{ A.18 }$$

$$\begin{eqnarray}&&w={{M}_{\alpha \beta }}{{n}_{\alpha }}{{n}_{\beta }}=0\qquad {\rm on}\;{{\Gamma }_{s}},\end{eqnarray} \tag{ A.19 }$$

$$\begin{eqnarray}&&{{M}_{\alpha \beta }}{{n}_{\alpha }}{{n}_{\beta }}=\frac{\partial }{\partial s}({{M}_{\alpha \beta }}{{t}_{\alpha }}{{n}_{\beta }})+{{M}_{\alpha \beta ,\alpha }}{{n}_{\beta }}=0\qquad {\rm on}\;{{\Gamma }_{f}}.\end{eqnarray} \tag{ A.20 }$$

This is the Euler–Lagrange equation associated with the minimization of (A.14).

A.2.1. Mixed variational principles

In a standard variational principle, the weak form equations are found by minimizing the energy functional. In a typical mixed variational principle, a dual variable is selected and a new functional is introduced. For the simple plate bending problem above the dual variable is usually selected to be the bending moment tensor ${{M}_{\alpha \beta }}$.⁵ The variational functional is a version of the Hellinger–Reissner principle [44, 47], given by

$$\begin{eqnarray}&&\mathcal{H}[w,{{M}_{\alpha \beta }}]={{\int }_{\Omega }}\left[ -\frac{1}{2D}{{\mathbb{E}}_{\alpha \beta \gamma \delta }}{{M}_{\alpha \beta }}{{M}_{\gamma \delta }}+{{w}_{,\alpha \beta }}{{M}_{\alpha \beta }}+pw \right]\;{{{\rm d}}^{2}}{\boldsymbol{X}} .\end{eqnarray} \tag{ A.21 }$$

Here

$$\begin{eqnarray}&&{{\mathbb{E}}_{\alpha \beta \gamma \delta }}=\frac{1}{(1-{{\sigma }^{2}})}\left[ \frac{(1+\sigma )}{2}({{\delta }_{\alpha \gamma }}{{\delta }_{\beta \delta }}+{{\delta }_{\alpha \delta }}{{\delta }_{\beta \gamma }})-\sigma {{\delta }_{\alpha \beta }}{{\delta }_{\gamma \delta }} \right]\end{eqnarray} \tag{ A.22 }$$

is the inverse of ${{\mathbb{B}}_{\alpha \beta \gamma \delta }}$. The weak form equations are derived from this principle by finding the stationary value of $\mathcal{H}$ over all functions ${{M}_{\alpha \beta }}$ and w satisfying the aforementioned conditions on Γ_c and Γ_s. Note that this stationary value of $\mathcal{H}$ will be neither a minimum nor a maximum; it is for this reason that these methods are often called saddlepoint methods. The weak form equations are found by setting the variation $\delta \mathcal{H}$ to zero, where

$$\begin{eqnarray}&&\delta \mathcal{H}={{\int }_{\Omega }}\left[ -\frac{1}{D}{{\mathbb{E}}_{\alpha \beta \gamma \delta }}{{\tilde{M}}_{\alpha \beta }}{{M}_{\gamma \delta }}+{{w}_{,\alpha \beta }}{{\tilde{M}}_{\alpha \beta }}+{{\tilde{w}}_{,\alpha \beta }}{{M}_{\alpha \beta }}+p\tilde{w} \right]\;{{{\rm d}}^{2}}{\boldsymbol{X}} .\end{eqnarray} \tag{ A.23 }$$

Performing integration by parts one may once more recover the strong formulation (A.17) with the correct boundary conditions.

However, the formulation (A.21) still requires a certain regularity of the deflection $\tilde{w}$; broadly speaking the square of its second derivative must be integrable. Meanwhile, the only regularity required of ${{M}_{\alpha \beta }}$ is that its square must be integrable. One of the main advantages of the mixed variational approach is that it allows regularity to be transferred from the displacement to the moment. On integrating (A.21) by parts we obtain

$$\begin{eqnarray}&&\begin{array}{l} \mathcal{H}[w,{{M}_{\alpha \beta }}] \\ \qquad ={{\int }_{\Omega }}\left[ -\frac{1}{2D}{{\mathbb{E}}_{\alpha \beta \gamma \delta }}{{M}_{\alpha \beta }}{{M}_{\gamma \delta }}-{{w}_{,\alpha }}{{M}_{\alpha \beta ,\beta }}+pw \right]\;{{{\rm d}}^{2}}{\boldsymbol{X}} +{{\oint }_{\partial \Omega }}{{M}_{\alpha \beta }}{{w}_{,\alpha }}{{n}_{\beta }}\;\;{\rm d}s. \\ \end{array}\end{eqnarray} \tag{ A.24 }$$

This functional is one that can be minimized over the space of all w and ${{M}_{\alpha \beta }}$ whose first derivatives are square-integrable. However, such w are unable to account for zero normal-derivatives on the boundary, so we encode that information directly in (A.24): since $w={{w}_{n}}=0$ on Γ_c, the boundary integral in (A.24) is zero over Γ_c, and hence

$$\begin{eqnarray}&&\begin{array}{l} \mathcal{H}[w,{{M}_{\alpha \beta }}] \\ \quad ={{\int }_{\Omega }}\left[ -\frac{1}{2D}{{\mathbb{E}}_{\alpha \beta \gamma \delta }}{{M}_{\alpha \beta }}{{M}_{\gamma \delta }}-{{w}_{,\alpha }}{{M}_{\alpha \beta ,\beta }}+pw \right]\;{{{\rm d}}^{2}}{\boldsymbol{X}} +{{\int }_{{{\Gamma }_{s}}\cup {{\Gamma }_{f}}}}{{M}_{\alpha \beta }}{{w}_{,\alpha }}{{n}_{\beta }}\;\;{\rm d}s. \\ \end{array}\end{eqnarray} \tag{ A.25 }$$

The weak form equations are obtained by finding the stationary value of (A.25) over the space of admissible functions satisfying w = 0 on ${{\Gamma }_{s}}\cup {{\Gamma }_{f}}$; in other words

$$\begin{eqnarray}&&\begin{array}{l} {{\int }_{\Omega }}\left[ -\frac{1}{D}{{\mathbb{E}}_{\alpha \beta \gamma \delta }}{{\tilde{M}}_{\alpha \beta }}{{M}_{\gamma \delta }}-{{\tilde{w}}_{,\alpha }}{{M}_{\alpha \beta ,\beta }}-{{w}_{,\alpha }}{{\tilde{M}}_{\alpha \beta ,\beta }}+p\tilde{w} \right]\;{{{\rm d}}^{2}}{\boldsymbol{X}} \\ \qquad +{{\int }_{{{\Gamma }_{s}}\cup {{\Gamma }_{f}}}}({{\tilde{w}}_{,\alpha }}{{M}_{\alpha \beta }}+{{w}_{,\alpha }}{{\tilde{M}}_{\alpha \beta }}){{n}_{\beta }}\;\;{\rm d}s=\;0 \\ \end{array}\end{eqnarray} \tag{ A.26 }$$

for all trial functions $\tilde{w}$ and ${{\tilde{M}}_{\alpha \beta }}$ satisfying these conditions. Again, these weak form equations lead to the same strong form (A.17) together with appropriate boundary conditions.

This elementary exposition has omitted technical details regarding the regularity of the solutions; for a more rigorous consideration the reader is referred to [44, 47], and [48].

A.2.2. Application to nonlinear plate bending

The application of these mixed variational principles to nonlinear plates was first analyzed by [45] and [46]. They were interested in developing numerical methods to study the buckling of compressed plates. This meant that their boundary conditions were ones of applied force, which allowed them to use an Airy stress function approach, leading to coupled fourth-order differential equations. We are unable to use these equations directly as our boundary conditions are ones of zero-displacement, which is difficult to express in terms of the stress function. Instead, we simply add the in-plane stored elastic energy to the variational formulation (A.25), together with an arbitrary external potential. Writing ${{N}_{\alpha \beta }}=C{{\mathbb{A}}_{\alpha \beta \gamma \delta }}{{\varepsilon }_{\gamma \delta }}$ for simplicity, where

$$\begin{eqnarray}&&{{\varepsilon }_{\alpha \beta }}=\frac{1}{2}\left( {{v}_{\alpha ,\beta }}+{{v}_{\beta ,\alpha }}+{{w}_{,\alpha }}{{w}_{,\beta }} \right),\end{eqnarray} \tag{ A.27 }$$

the mixed variational functional we use is

$$\begin{eqnarray}\begin{array}{rcl} \mathcal{H}[{{v}_{\kappa }},w,{{M}_{\kappa \lambda }}] & = & {{\int }_{\Omega }}\left( \frac{1}{2}{{N}_{\alpha \beta }}{{\varepsilon }_{\alpha \beta }}-\frac{1}{2D}{{\mathbb{E}}_{\alpha \beta \gamma \delta }}{{M}_{\alpha \beta }}{{M}_{\gamma \delta }} \right)\;{{{\rm d}}^{2}}{\boldsymbol{X}} \\ {} & {} & +{{\int }_{\Omega }}\left( -{{w}_{,\alpha }}{{M}_{\alpha \beta ,\beta }}+V[{{v}_{\kappa }},w] \right)\;{{{\rm d}}^{2}}{\boldsymbol{X}} +{{\int }_{{{\Gamma }_{s}}\cup {{\Gamma }_{f}}}}{{w}_{,\alpha }}{{M}_{\alpha \beta }}{{n}_{\beta }}\;\;{\rm d}s. \\ \end{array}\end{eqnarray} \tag{ A.28 }$$

According to the discussion above, in order to derive the weak form equations we should find the stationary value of $\mathcal{H}$ over all admissible ${{v}_{\alpha }}$, w, ${{M}_{\alpha \beta }}$ satisfying ${{v}_{\alpha }}=w=0$ on ${{\Gamma }_{c}}\cup {{\Gamma }_{s}}$. The first variation of $\mathcal{H}$ can be straightforwardly derived along the same lines discussed above, whereupon one obtains

$$\begin{eqnarray}\begin{array}{rcl} \delta \mathcal{H} & = & {{\int }_{\Omega }}\left( {{\tilde{v}}_{\alpha ,\beta }}{{N}_{\alpha \beta }}+{{\tilde{v}}_{\alpha }}{{V}_{{{v}_{\alpha }}}} \right)\;{{{\rm d}}^{2}}{\boldsymbol{X}} \\ {} & {} & +{{\int }_{\Omega }}\left( -{{\tilde{M}}_{\alpha \beta }}\frac{1}{D}{{\mathbb{E}}_{\alpha \beta \gamma \delta }}{{M}_{\gamma \delta }}-{{\tilde{M}}_{\alpha \beta ,\beta }}{{w}_{,\alpha }} \right)\;{{{\rm d}}^{2}}{\boldsymbol{X}} +{{\int }_{{{\Gamma }_{s}}\cup {{\Gamma }_{f}}}}{{\tilde{M}}_{\alpha \beta }}{{w}_{,\alpha }}{{n}_{\beta }}\;\;{\rm d}s \\ {} & {} & +{{\int }_{\Omega }}\left[ {{\tilde{w}}_{,\alpha }}({{N}_{\alpha \beta }}{{w}_{,\beta }}-{{M}_{\alpha \beta ,\beta }})+\tilde{w}{{V}_{w}} \right]\;{{{\rm d}}^{2}}{\boldsymbol{X}} +{{\int }_{{{\Gamma }_{s}}\cup {{\Gamma }_{f}}}}{{\tilde{w}}_{,\alpha }}{{M}_{\alpha \beta }}{{n}_{\beta }}\;\;{\rm d}s. \\ \end{array}\end{eqnarray} \tag{ A.29 }$$

These equations then lead naturally to the weak form equations (6a )–(6f ) on assuming that the entire boundary is clamped (${{\Gamma }_{s}}={{\Gamma }_{f}}=\varnothing$), that $\sigma =\nu$, and on writing out the equations for the six components $\tilde{w}$, ${{\tilde{v}}_{\alpha }}$ and ${{\tilde{M}}_{\alpha \beta }}$ explicitly.

The equations hold for all continuous integrable variations ${{\tilde{v}}_{\alpha }}$, $\tilde{w}$, ${{\tilde{M}}_{\alpha \beta }}$ that satisfy ${{\tilde{v}}_{\alpha }}=\tilde{w}=0$ on the boundary.

To the six weak-form elastic equations derived in appendix A.2 we must add an equation linking the strain field to the generated PMF, $B({\boldsymbol{X}} )$. This is because we wish to find the deformation field that best approximates $B({\boldsymbol{X}} )$ to a desired (target), $\hat{B}({\boldsymbol{X}} )$. The origin of this PMF that appears in the low-energy effective Hamiltonian of deformed graphene is the local modification of the electronic hopping amplitudes between neighboring carbon atoms brought about by the space dependent deformation of the crystal lattice.

A single orbital nearest-neighbor tight-binding model for the $\pi$ bands derived from electronic hopping among p_z orbitals of neighboring carbons has been extremely successful in describing the behavior of electrons in graphene, and their response to various kinds of external perturbations and fields [17]. The Hamiltonian that reflects this physics is given by

$$\begin{eqnarray}&&H=-\mathop{\sum }\limits_{i,{\boldsymbol{n}} }t\left( {{{\boldsymbol{X}} }_{i}},{{{\boldsymbol{X}} }_{i}}+{\boldsymbol{n}} \right)\;a_{{{{\boldsymbol{X}} }_{i}}}^{\dagger }b_{{{{\boldsymbol{X}} }_{i}}+{\boldsymbol{n}} }^{{}}+{\rm H}.{\rm c}..\end{eqnarray} \tag{ B.1 }$$

The bipartite nature of the honeycomb lattice is evident in this expression by the explicit distinction between the lattice sites belonging to sub-lattice A or B. The second-quantized operator ${{a}_{{{{\boldsymbol{X}} }_{i}}}}$(${{b}_{{{{\boldsymbol{X}} }_{i}}}}$) destroys an electron in a p_z orbital that belongs to a carbon atom located on site A(B) of the unit cell placed at ${\boldsymbol{X}}$_i. The parameter $t\left( {{{\boldsymbol{X}} }_{i}},{{{\boldsymbol{X}} }_{i}}+{\boldsymbol{n}} \right)$ is the hopping amplitude between two neighboring $\pi$ orbitals, and ${\boldsymbol{n}}$ runs over the three unit cells containing a B atom neighboring the A atom from the unit cell at ${\boldsymbol{X}}$_i. The hopping amplitude is constant in the perfect crystal: $t\left( {{{\boldsymbol{X}} }_{i}},{{{\boldsymbol{X}} }_{i}}+{\boldsymbol{n}} \right)={{t}_{0}}=2.7$ eV. But, since t depends strongly on the inter-atomic distance, any local change caused by a deformation leads to perturbations to this equilibrium value and, hence, more generically, $t\left( {{{\boldsymbol{X}} }_{i}},{{{\boldsymbol{X}} }_{i}}+{\boldsymbol{n}} \right)={{t}_{0}}+\delta t\left( {{{\boldsymbol{X}} }_{i}},{{{\boldsymbol{X}} }_{i}}+{\boldsymbol{n}} \right)$. The presence of $\delta t$, which is a relatively small perturbation to t₀ in practical situations, adds a correction to the low-energy Dirac-like Hamiltonian that emerges from (B.1) so that the effective Hamiltonian around the point ${\boldsymbol{K}} =(4\pi /3\sqrt{3}a,\;0)$ in the first Brillouin zone has the form [15, 16]

$$\begin{eqnarray}&&\mathcal{H}={{v}_{F}}\;{\boldsymbol{ \sigma }} \cdot \left( {\boldsymbol{p}} +e \right){\bf{\mathcal{A}}},\end{eqnarray} \tag{ B.2 }$$

where ${\boldsymbol{ \sigma }}$ is a vector of Pauli matrices, and $\hbar {{v}_{F}}=3{{t}_{0}}a/2$, with a = 1.42 Å the carbon–carbon distance in equilibrium. For deformations on scales that are large compared to a, the curvature-induced tilting of neighboring p_z orbitals can be neglected⁶ . In this situation the hopping amplitude t depends only on the distance between neighboring atoms, and we straightforwardly obtain the components of the vector potential ${\bf{\mathcal{A}}}={{\mathcal{A}}_{x}}{{{\boldsymbol{e}} }_{x}}+{{\mathcal{A}}_{y}}{{{\boldsymbol{e}} }_{y}}$ by expanding t to linear order in the deformation tensor. Choosing the coordinate system so that ${{{\boldsymbol{e}} }_{x}}$ is along the zig-zag direction of the honeycomb lattice one obtains [16]

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{A}}_{x}}({\boldsymbol{X}} )-{\rm i}{{\mathcal{A}}_{y}}({\boldsymbol{X}} ) & \equiv & -\frac{1}{e{{v}_{F}}}\mathop{\sum }\limits_{{\boldsymbol{n}} }\delta t\left( {\boldsymbol{X}} ,{\boldsymbol{X}} +{\boldsymbol{n}} \right){{{\rm e}}^{{\rm i}{\boldsymbol{K}} \cdot {\boldsymbol{n}} }} \\ {} & \simeq & -\frac{\hbar c}{2ea}\left( {{\varepsilon }_{xx}}-{{\varepsilon }_{yy}}+2{\rm i}\;{{\varepsilon }_{xy}} \right), \\ \end{array}\end{eqnarray} \tag{ B.3 }$$

where $c=-\partial {\rm log} t(r)/\partial {\rm log} r{{|}_{r=a}}$. For static deformations, a value $c\approx 3.37$ captures the changes in various physical properties arising from strain-induced modifications of the $\pi$ bands in agreement with first-principles calculations [49–54]. Finally, the PMF B, being defined as the 2D curl of $\bf{\mathcal{A}}$, reads:

$$\begin{eqnarray}&&B({\boldsymbol{X}} )=\frac{\hbar c}{ae}\left[ \frac{\partial }{\partial Y}\left( \frac{{{\varepsilon }_{11}}-{{\varepsilon }_{22}}}{2} \right)+\frac{\partial {{\varepsilon }_{12}}}{\partial X} \right].\end{eqnarray} \tag{ B.4 }$$

As noted in the previous section, by virtue of our choice of piecewise affine finite elements for the numerical interpolation, the six variables v₁, v₂, w, M₁₁, M₁₂, and M₂₂ are treated as continuous, and affine over each triangular element. As a consequence, from (A.4) the strain field will be discontinuous and constant in each triangular element. As the strain components ${{\varepsilon }_{\alpha \beta }}$ are discontinuous under this approximation, the PMF (B.4) using this scheme is undefined. To overcome this difficulty, we use the technique of patch recovery. For details of the technique, first described by [58], we refer the reader to appendix E. In brief, this is a mechanism that uses the discontinuous strain data ${{\varepsilon }_{\alpha \beta }}$ to recover a strain field $\varepsilon _{\alpha \beta }^{{\rm rec}}$ of the same type as the primary variables: continuous and affine over each element. The derivative of $\varepsilon _{\alpha \beta }^{{\rm rec}}$ is well-defined, and thus so is the PMF if it is calculated using this recovered strain field:

$$\begin{eqnarray}&&B({\boldsymbol{X}} )=\frac{\hbar c}{ae}\left[ \frac{\partial }{\partial Y}\left( \frac{\varepsilon _{11}^{{\rm rec}}-\varepsilon _{22}^{{\rm rec}}}{2} \right)+\frac{\partial \varepsilon _{12}^{{\rm rec}}}{\partial X} \right].\end{eqnarray} \tag{ B.5 }$$

In solving the optimization problem of section 3 numerically, the first step is to non-dimensionalize the system of equations in such a way that most variables are $O(1)$ to ensure good numerical behavior. To accomplish this we choose the following scalings, where an overbar represents the dimensionless quantity. Set

$$\begin{eqnarray}&&\varepsilon =\frac{ae{{B}_{0}}L}{\hbar c}\end{eqnarray} \tag{ C.1 }$$

to be the typical scaling of the strain field, then

$$\begin{eqnarray}&&(X,Y)=L(\bar{X},\bar{Y}),\quad {{v}_{\alpha }}=L\varepsilon {{\bar{v}}_{\alpha }},\quad w=L\sqrt{\varepsilon }\bar{w},\end{eqnarray} \tag{ C.2 }$$

$$\begin{eqnarray}&&{{\varepsilon }_{\alpha \beta }}=\varepsilon {{\bar{\varepsilon }}_{\alpha \beta }},\quad \varepsilon _{\alpha \beta }^{{\rm rec}}=\varepsilon \bar{\varepsilon }_{\alpha \beta }^{\,{\rm rec}},\quad {{N}_{\alpha \beta }}=C\varepsilon {{\bar{N}}_{\alpha \beta }},\quad {{M}_{\alpha \beta }}=\frac{D\sqrt{\varepsilon }}{L}{{\bar{M}}_{\alpha \beta }},\end{eqnarray} \tag{ C.3 }$$

$$\begin{eqnarray}&&V=C{{\varepsilon }^{2}}\bar{V},\quad B={{B}_{0}}\bar{B},\quad \hat{B}={{B}_{0}}\overline{{\hat{B}}}.\end{eqnarray} \tag{ C.4 }$$

All these constants are previously-defined, with the exception of L, representing the typical size of the domain Ω, and B₀, the typical magnitude of the target PMF $\hat{B}$. Under these scalings the equations exhibit only one dimensionless parameter, namely the dimensionless bending stiffness $\kappa$:

$$\begin{eqnarray}&&\kappa =\frac{D}{C{{L}^{2}}\varepsilon }=\frac{D\hbar c}{C{{L}^{3}}ae{{B}_{0}}}.\end{eqnarray} \tag{ C.5 }$$

For completeness we will summarize the minimization problem in its dimensionless form:

$$\begin{eqnarray}&&{\rm Minimize}\qquad \bar{\mathcal{I}}={{\int }_{{\bar{\Omega }}}}{{\left( \bar{B}-\overline{{\hat{B}}}(\bar{X},\bar{Y}) \right)}^{2}}\;{{{\rm d}}^{2}}{\boldsymbol{\bar{X}}} +\eta \;{{\bar{\mathcal{I}}}^{{\rm reg}}}[{{\lambda }_{i}}],\end{eqnarray} \tag{ C.6 }$$

subject to the six equations

$$\begin{eqnarray}&&{{\int }_{{\bar{\Omega }}}}\left[ \frac{\partial {{\tilde{v}}_{1}}}{\partial \bar{X}}{{{\bar{N}}}_{11}}+\frac{\partial {{\tilde{v}}_{1}}}{\partial \bar{Y}}{{{\bar{N}}}_{12}}+{{\tilde{v}}_{1}}{{{\bar{V}}}_{{{{\bar{v}}}_{1}}}}[\bar{w},{{{\bar{v}}}_{1}},{{{\bar{v}}}_{2}};{{\lambda }_{i}}] \right]\;{{{\rm d}}^{2}}{\boldsymbol{\bar{X}}} =0,\end{eqnarray} \tag{ C.7 }$$

$$\begin{eqnarray}&&{{\int }_{{\bar{\Omega }}}}\left[ \frac{\partial {{\tilde{v}}_{2}}}{\partial \bar{X}}{{{\bar{N}}}_{12}}+\frac{\partial {{\tilde{v}}_{2}}}{\partial \bar{Y}}{{{\bar{N}}}_{22}}+{{\tilde{v}}_{2}}{{{\bar{V}}}_{{{{\bar{v}}}_{2}}}}[\bar{w},{{{\bar{v}}}_{1}},{{{\bar{v}}}_{2}};{{\lambda }_{i}}] \right]\;{{{\rm d}}^{2}}{\boldsymbol{\bar{X}}} =0,\end{eqnarray} \tag{ C.8 }$$

$$\begin{eqnarray}&&\begin{array}{l} {{\int }_{{\bar{\Omega }}}}\left\{ \frac{\partial \tilde{w}}{\partial \bar{X}}\left[ -\kappa \left( \frac{\partial {{{\bar{M}}}_{11}}}{\partial \bar{X}}+\frac{\partial {{{\bar{M}}}_{12}}}{\partial \bar{Y}} \right)+{{{\bar{N}}}_{11}}\frac{\partial \bar{w}}{\partial \bar{X}}+{{{\bar{N}}}_{12}}\frac{\partial \bar{w}}{\partial \bar{Y}} \right] \right. \\ \quad +\frac{\partial \tilde{w}}{\partial \bar{Y}}\left[ -\kappa \left( \frac{\partial {{{\bar{M}}}_{12}}}{\partial \bar{X}}+\frac{\partial {{{\bar{M}}}_{22}}}{\partial \bar{Y}} \right)+{{{\bar{N}}}_{12}}\frac{\partial \bar{w}}{\partial \bar{X}}+{{{\bar{N}}}_{22}}\frac{\partial \bar{w}}{\partial \bar{Y}} \right] \\ \quad +\left. \tilde{w}{{{\bar{V}}}_{{\bar{w}}}}[\bar{w},{{{\bar{v}}}_{1}},{{{\bar{v}}}_{2}};{{\lambda }_{i}}] \right\}{{{\rm d}}^{2}}{\boldsymbol{\bar{X}}} =0, \\ \end{array}\end{eqnarray} \tag{ C.9 }$$

$$\begin{eqnarray}&&{{\int }_{{\bar{\Omega }}}}\left[ \frac{1}{(1-{{\sigma }^{2}})}({{{\bar{M}}}_{11}}-\sigma {{{\bar{M}}}_{22}}){{\tilde{M}}_{11}}+\frac{\partial \bar{w}}{\partial \bar{X}}\frac{\partial {{\tilde{M}}_{11}}}{\partial \bar{X}} \right]\;{{{\rm d}}^{2}}{\boldsymbol{\bar{X}}} =0,\end{eqnarray} \tag{ C.10 }$$

$$\begin{eqnarray}&&{{\int }_{{\bar{\Omega }}}}\left[ \frac{1}{(1-\sigma )}{{{\bar{M}}}_{12}}{{\tilde{M}}_{12}}+\frac{1}{2}\frac{\partial \bar{w}}{\partial \bar{X}}\frac{\partial {{\tilde{M}}_{12}}}{\partial \bar{Y}}+\frac{1}{2}\frac{\partial \bar{w}}{\partial \bar{Y}}\frac{\partial {{\tilde{M}}_{12}}}{\partial \bar{X}} \right]\;{{{\rm d}}^{2}}{\boldsymbol{\bar{X}}} =0,\end{eqnarray} \tag{ C.11 }$$

$$\begin{eqnarray}&&{{\int }_{{\bar{\Omega }}}}\left[ \frac{1}{(1-{{\sigma }^{2}})}({{{\bar{M}}}_{22}}-\sigma {{{\bar{M}}}_{11}}){{\tilde{M}}_{22}}+\frac{\partial \bar{w}}{\partial \bar{Y}}\frac{\partial {{\tilde{M}}_{22}}}{\partial \bar{Y}} \right]\;{{{\rm d}}^{2}}{\boldsymbol{\bar{X}}} =0,\end{eqnarray} \tag{ C.12 }$$

together with the additional definitions

$$\begin{eqnarray}&&{{\bar{N}}_{11}}={{\bar{\varepsilon }}_{11}}+\nu {{\bar{\varepsilon }}_{22}},\qquad {{\bar{N}}_{12}}=(1-\nu ){{\bar{\varepsilon }}_{12}},\qquad {{\bar{N}}_{22}}=\nu {{\bar{\varepsilon }}_{11}}+{{\bar{\varepsilon }}_{22}},\end{eqnarray} \tag{ C.13 }$$

$$\begin{eqnarray}&&{{\bar{\varepsilon }}_{11}}=\frac{\partial {{{\bar{v}}}_{1}}}{\partial \bar{X}}+\frac{1}{2}{{\left( \frac{\partial \bar{w}}{\partial \bar{X}} \right)}^{2}},\quad {{\bar{\varepsilon }}_{22}}=\frac{\partial {{{\bar{v}}}_{2}}}{\partial \bar{Y}}+\frac{1}{2}{{\left( \frac{\partial \bar{w}}{\partial \bar{Y}} \right)}^{2}},\end{eqnarray} \tag{ C.14 }$$

$$\begin{eqnarray}&&{{\bar{\varepsilon }}_{12}}=\frac{1}{2}\left( \frac{\partial {{{\bar{v}}}_{1}}}{\partial \bar{Y}}+\frac{\partial {{{\bar{v}}}_{2}}}{\partial \bar{X}}+\frac{\partial \bar{w}}{\partial \bar{X}}\frac{\partial \bar{w}}{\partial \bar{Y}} \right),\end{eqnarray} \tag{ C.15 }$$

$$\begin{eqnarray}&&\bar{B}=\frac{\partial }{\partial \bar{Y}}\left( \frac{\bar{\varepsilon }_{11}^{\,{\rm rec}}-\bar{\varepsilon }_{22}^{{\rm rec}}}{2} \right)+\frac{\partial \bar{\varepsilon }_{12}^{\,{\rm rec}}}{\partial \bar{X}},\end{eqnarray} \tag{ C.16 }$$

and $\bar{\varepsilon }_{\alpha \beta }^{{\rm rec}}$ obtained from ${{\bar{\varepsilon }}_{\alpha \beta }}$ by strain recovery.

To make the expression (14) for the substrate–graphene interaction dimensionless, set

$$\begin{eqnarray}&&(x,y)=L(\bar{x},\bar{y}),\quad s=L\sqrt{\varepsilon }\bar{s},\quad \hat{z}(x,y)=L\sqrt{\varepsilon }\overline{{\hat{z}}}(\bar{x},\bar{y}),\quad J(s)={{J}_{0}}\bar{J}(\bar{s}),\end{eqnarray} \tag{ C.17 }$$

so that

$$\begin{eqnarray}&&\bar{V}[\bar{w},{{\bar{v}}_{1}},{{\bar{v}}_{2}};{{\lambda }_{i}}]=\bar{p}\bar{w}+{{\bar{J}}_{0}}\bar{J}\left( \bar{w}-\overline{{\hat{z}}}(\bar{X}+\varepsilon {{{\bar{v}}}_{1}},\bar{Y}+\varepsilon {{{\bar{v}}}_{2}}) \right),\end{eqnarray} \tag{ C.18 }$$

where

$$\begin{eqnarray}&&\bar{J}(\bar{s})=\frac{1}{3}\left[ {{\left( \frac{{{{\bar{s}}}^{*}}}{{\bar{s}}} \right)}^{8}}-4{{\left( \frac{{{{\bar{s}}}^{*}}}{{\bar{s}}} \right)}^{2}} \right],\end{eqnarray} \tag{ C.19 }$$

and the three dimensionless parameters are given by

$$\begin{eqnarray}&&\bar{p}=\frac{pL}{C{{\varepsilon }^{3/2}}},\qquad {{\bar{J}}_{0}}=\frac{{{J}_{0}}}{C{{\varepsilon }^{2}}},\qquad {{\bar{s}}^{*}}=\frac{{{s}^{*}}}{L\sqrt{\varepsilon }}.\end{eqnarray} \tag{ C.20 }$$

As a representative value for s^*, we use the value that [63] give for C–Cu, namely $2.243\;\overset{\circ}{A}$. Similarly we use ${{J}_{0}}=0.45\;{\rm J}\;{{{\rm m}}^{-2}}$ as a representative value, from the investigation of [7] into the adhesion strength between graphene and SiO$_{2}$. We select a typical value for the hydrostatic pressure as $p=100\;{\rm bar}={{10}^{7}}\;{\rm Pa}$. These values, for a lengthscale $L={{10}^{-8}}\;{\rm m}$ and a target PMF scale ${{B}_{0}}=10\;{\rm T}$, give $\bar{p}=0.545$, ${{\bar{J}}_{0}}=30.66$ and ${{\bar{s}}^{*}}=0.280$.

The derivatives of the potentials appearing in the dimensionless weak form equations are

$$\begin{eqnarray}&&{{\bar{V}}_{{\bar{w}}}}[\bar{w},{{\bar{v}}_{1}},{{\bar{v}}_{2}};{{\lambda }_{i}}]=\bar{p}+{{\bar{J}}_{0}}\bar{J}^{\prime} (\bar{w}-\overline{{\hat{z}}}(\bar{X}+\varepsilon {{\bar{v}}_{1}},\bar{Y}+\varepsilon {{\bar{v}}_{2}})),\end{eqnarray} \tag{ C.21 }$$

$$\begin{eqnarray}&&{{\bar{V}}_{{{{\bar{v}}}_{\alpha }}}}[\bar{w},{{\bar{v}}_{1}},{{\bar{v}}_{2}};{{\lambda }_{i}}]=-\varepsilon {{\bar{J}}_{0}}\bar{J}^{\prime} (\bar{w}-\overline{{\hat{z}}}(\bar{X}+\varepsilon {{\bar{v}}_{1}},\bar{Y}+\varepsilon {{\bar{v}}_{2}})){{\left. \frac{\partial \overline{{\hat{z}}}}{\partial {{{\bar{x}}}_{\alpha }}} \right|}_{(\bar{X}+\varepsilon {{{\bar{v}}}_{1}},\bar{Y}+\varepsilon {{{\bar{v}}}_{2}})}},\end{eqnarray} \tag{ C.22 }$$

where

$$\begin{eqnarray}&&\bar{J}^{\prime} (\bar{s})=\frac{8}{3{{{\bar{s}}}^{*}}}\left[ {{\left( \frac{{{{\bar{s}}}^{*}}}{{\bar{s}}} \right)}^{3}}-{{\left( \frac{{{{\bar{s}}}^{*}}}{{\bar{s}}} \right)}^{9}} \right].\end{eqnarray} \tag{ C.23 }$$

In section 4.2 and subsequently, all variables are assumed to be dimensionless, and overbars are omitted for clarity.

For the numerical experiments in this article, we assume that the substrate is patterned periodically in the two horizontal directions, and set the repeating unit cell to be the rhombus depicted in figure D1 . If we introduce the two coordinates

$$\begin{eqnarray}&&{{\xi }_{1}}=\frac{y}{\sqrt{3}}+x,\qquad {{\xi }_{2}}=\frac{y}{\sqrt{3}}-x,\end{eqnarray} \tag{ D.1 }$$

the unit cell corresponds to the domain

$$\begin{eqnarray}&&\Omega =\left\{ ({{\xi }_{1}},{{\xi }_{2}})\ :\ 0\leqslant {{\xi }_{1}}\lt 1,\ 0\leqslant {{\xi }_{2}}\lt 1 \right\}.\end{eqnarray} \tag{ D.2 }$$

The topography of the substrate, $\hat{z}(x,y)$, can then be resolved in terms of a sum of functions which are periodic on the unit cell Ω. Such functions take one of the four following forms:

$$\begin{eqnarray}&&f_{kl}^{1}={\rm cos} (2\pi k{{\xi }_{1}}){\rm cos} (2\pi l{{\xi }_{2}}),\end{eqnarray} \tag{ D.3 }$$

$$\begin{eqnarray}&&f_{kl}^{2}={\rm cos} (2\pi k{{\xi }_{1}}){\rm sin} (2\pi l{{\xi }_{2}}),\end{eqnarray} \tag{ D.4 }$$

$$\begin{eqnarray}&&f_{kl}^{3}={\rm sin} (2\pi k{{\xi }_{1}}){\rm cos} (2\pi l{{\xi }_{2}}),\end{eqnarray} \tag{ D.5 }$$

$$\begin{eqnarray}&&f_{kl}^{4}={\rm sin} (2\pi k{{\xi }_{1}}){\rm sin} (2\pi l{{\xi }_{2}}),\end{eqnarray} \tag{ D.6 }$$

and so we choose a truncated expansion in terms of these as follows:

$$\begin{eqnarray}&&\hat{z}(x,y)=\mathop{\sum }\limits_{k=0}^{K}\mathop{\sum }\limits_{l=0}^{K}\left( {{\alpha }_{kl}}f_{kl}^{1}+{{\beta }_{kl}}f_{kl}^{2}+{{\gamma }_{kl}}f_{kl}^{3}+{{\delta }_{kl}}f_{kl}^{4} \right),\end{eqnarray} \tag{ D.7 }$$

with the constants ${{\alpha }_{kl}}$, ${{\beta }_{kl}}$, ${{\gamma }_{kl}}$, ${{\delta }_{kl}}$ playing the role of the control variables ${{\lambda }_{i}}$: varying the topography of the substrate is, therefore, achieved by varying these $4{{(K+1)}^{2}}$ constants. For convenience we set

$$\begin{eqnarray}&&{{\beta }_{k0}}={{\gamma }_{0k}}={{\delta }_{k0}}={{\delta }_{0k}}=0\end{eqnarray} \tag{ D.8 }$$

for each k and, since rigid vertical displacements of the substrate do not affect the objective function, we further set ${{\alpha }_{00}}=0$.

Zoom In Zoom Out Reset image size

In order to avoid convergence towards solutions that are ill-behaved (in this case those could be, for example, substrate profiles with discontinuities or sharp topographical features) during the numerical optimization, a regularization term, ${{\mathcal{I}}^{{\rm reg}}}$, is added to the objective integral, as per equation (12). We choose it to be

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{I}}^{{\rm reg}}} & = & \frac{1}{{\rm Area}(\Omega )}{{\int }_{\Omega }}{{\left| \nabla \hat{z} \right|}^{2}}{{{\rm d}}^{2}}{\boldsymbol{X}} \\ {} & = & \frac{4}{3}\int _{0}^{1}\int _{0}^{1}\left[ {{\left( \frac{\partial \hat{z}}{\partial {{\xi }_{1}}} \right)}^{2}}+{{\left( \frac{\partial \hat{z}}{\partial {{\xi }_{2}}} \right)}^{2}}-\frac{\partial \hat{z}}{\partial {{\xi }_{1}}}\frac{\partial \hat{z}}{\partial {{\xi }_{2}}} \right]\;\;{\rm d}{{\xi }_{1}}\;{\rm d}{{\xi }_{2}}, \\ \end{array}\end{eqnarray} \tag{ D.9 }$$

which is simple to calculate using the orthogonality of the basis functions f_klⁿ over Ω.

The domain Ω for the graphene sheet will also be the unit rhombus, with periodic boundary conditions applied to all six state variables. However, we set the displacement components ${{v}_{1}}={{v}_{2}}=0$ at the corner points to disallow arbitrarily-large horizontal rigid displacements. This is a reasonable constraint on account of the 2D periodicity of the substrate. Due to the geometry of the unit rhombus, we can set the triangulation to be a regular isometric grid.

As noted in appendix A, we choose a finite element discretization for our six variables v₁, v₂, w, M₁₁, M₁₂, M₂₂ that approximates these quantities with functions that are continuous across the domain Ω, and affine over each triangular element in the discretization (see figure E1 (a) for a representation of such a function). Thus the quantities can be parametrized by their values at each nodal point of the triangulation. Differentiating such a function leads to a discontinuous function, which is constant on each triangular element, as shown in figure E1(b).

Zoom In Zoom Out Reset image size

The piecewise constant function is a less accurate approximation than the continuous piecewise affine function, and this led to the patch recovery method [58], which reconstructs an accurate continuous piecewise affine representation of a quantity calculated as a piecewise constant function. The canonical example where this recovery method becomes relevant is in elasticity with piecewise affine displacements leading to a piecewise constant stress field. The original purpose of the patch recovery method was to find a better approximation to the stress field calculated from a displacement-based finite element method. In this article a piecewise affine displacement field leads to a piecewise constant approximation ${{\varepsilon }_{\alpha \beta }}$ to the strain field, whereas we require a differentiable approximation. By using the patch recovery method we recover a piecewise affine strain field, $\varepsilon _{\alpha \beta }^{{\rm rec}}$, which we are able to differentiate to find the PMF B according to the prescription in equation (B.5).

To illustrate the patch recovery method, consider a triangulation of the domain Ω which defines the spatial extent of the medium with triangles $k=1,\ldots ,{{N}_{t}}$ and nodes $i=1,\ldots ,{{N}_{p}}$. We have a function f, constant on each element (so $f(X,Y)={{f}_{k}}$ if ${\boldsymbol{X}}$ is in triangle k), from which we want to recover a piecewise affine function ${{f}^{{\rm rec}}}$ (defined by its values $f_{i}^{{\rm rec}}$ at each node ${{{\boldsymbol{X}} }_{i}}$ of the triangulation).

The strength of the patch recovery method is that the nodal values $f_{i}^{{\rm rec}}$ are calculated individually in turn, rather than in a global optimization over all values at once. For each nodal point i, we identify the patch, which (for triangular elements with a piecewise affine target) is the set of all elements that contain the node i as a vertex, as displayed in figure E2 . The key step in the process is to fit a function $f_{i}^{{\rm fit}}(X,Y)$ to the patch for node i that is of the same order as the proposed target function. So, in this case, we need to fit an affine function $f_{i}^{{\rm fit}}(X,Y)=a+bX+cY$ to the patch. We use the values f_k, evaluated at the centroids $(X_{k}^{c},Y_{k}^{c})$ of the elements, to calculate the parameters $(a,b,c)={{{\boldsymbol{a}} }^{T}}$ through a least-squares optimization. As noted in [58], ${\boldsymbol{a}}$ is thus found by solving the system

$$\begin{eqnarray}&&\left[ \mathop{\sum }\limits_{k\in {\rm patch}(i)}{{{\boldsymbol{p}} }_{k}}\;{\boldsymbol{p}} _{k}^{T} \right]\;{\boldsymbol{a}} =\mathop{\sum }\limits_{k\in {\rm patch}(i)}{{f}_{k}}\;{{{\boldsymbol{p}} }_{k}},\end{eqnarray} \tag{ E.1 }$$

where ${{{\boldsymbol{p}} }_{k}}={{(1,X_{k}^{c},Y_{k}^{c})}^{T}}$. Having found $f_{i}^{{\rm fit}}(X,Y)$, the nodal value of the recovered function ${{f}^{{\rm rec}}}$ is simply $f_{i}^{{\rm rec}}=f_{i}^{{\rm fit}}({{X}_{i}},{{Y}_{i}})$, the fit function evaluated at the nodal point. At the domain boundaries there will usually be too few elements in the patch for the system (E.1) to be well-conditioned. Instead, we would follow [58] and find the boundary nodal values of ${{f}^{{\rm rec}}}$ by using the interior patches, and average over all the calculated values. This consideration does not apply for periodic boundary conditions, since in that case we can treat the entire domain as being of infinite extent, and all points are interior points.

Zoom In Zoom Out Reset image size