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

One of the enticing features common to most of the two-dimensional electronic systems that are currently at the forefront of materials science research is the ability to easily introduce a combination of planar deformations and bending in the system. Since the electronic properties are ultimately determined by the details of atomic orbital overlap, such mechanical manipulations translate into modified electronic properties. Here, we present a general-purpose optimization framework for tailoring physical properties of two-dimensional electronic systems by manipulating the state of local strain, allowing a one-step route from their design to experimental implementation. A definite example, chosen for its relevance in light of current experiments in graphene nanostructures, is the optimization of the experimental parameters that generate a prescribed spatial profile of pseudomagnetic fields in graphene. But the method is general enough to accommodate a multitude of possible experimental parameters and conditions whereby deformations can be imparted to the graphene lattice, and complies, by design, with graphene's elastic equilibrium and elastic compatibility constraints. As a result, it efficiently answers the inverse problem of determining the optimal values of a set of external or control parameters that result in a graphene deformation whose associated pseudomagnetic field profile best matches a prescribed target. The ability to address this inverse problem in an expedited way is one key step for practical implementations of the concept of two-dimensional systems with electronic properties strain-engineered to order. The general-purpose nature of this calculation strategy means that it can be easily applied to the optimization of other relevant physical quantities which directly depend on the local strain field, not just in graphene but in other two-dimensional electronic membranes.


Introduction
With their intrinsic two-dimensionality, "electronic membranes" are easily pulled or pinched by atomic-scale tips [1,2,3], can be made to conform to the substrate topography [4,5,6], can be inflated as balloons [7], can be stretched [8] or bent [9], crumpled on demand [10], and so on. Hence, two-dimensional crystals are an excellent case (and opportunity) of 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 twodimensional 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 (PMF) [15,16,17,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 (H = v F p·σ 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 p → p+eA, where A encodes all the details of the deformation and how it perturbs the electronic hopping amplitudes (defined below) [15]. As a result, even though 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, A will have an opposite sign for the effective Hamiltonian at the time-reversal transformed K 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 spacedependent Fermi velocities have also been reported in recent experiments on strained graphene [21], bringing this other theoretical prediction [22,23,24] and implication of (a) (b) Figure 1.
The two dimensional solution of Guinea et al. [26]. (a) An initially circular and isotropic graphene sheet is deformed to a rounded triangular shape. (b) The magnitude and direction of local stretch is indicated by the ellipses, which are the images under the deformation of small circles in the undeformed sheet.
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 [25,26,6,27,28].
Since the electronic dynamics can be straightforwardly determined once a spacedependent (pseudo)magnetic field B(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 strainengineering 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 A from B = ∇×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 Guinea et al. [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 3-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].
Extending deformations to three dimensions introduces nonlinearity into the strain pattern width pattern depth pressure Figure 2. Experimental situation to be modeled and taken throughout this report as a practical example of our method to solve the inverse strain-engineering problem in graphene. Graphene is placed on a patterned substrate with which it interacts via Van der Waals forces. Hydrostatic pressure and substrate profile are the two control parameters here, and the former is used to control the degree of conformation of graphene the substrate.
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,34,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(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 Guinea et al. [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.
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.

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 two-dimensional surface. The point X = (X, Y, 0) is transformed to r = (x 1 , x 2 , x 3 ) = (x, y, z) in three-dimensional space, where x α = X α + v α (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 ε αβ and ρ αβ . 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 ε αβ and ρ αβ is included in Appendix A). The most common simplification is perhaps the von Kármán approximation, which uses the expressions The stress and moment resultants are assumed to be isotropic and linear in the strain and curvature tensors: where we have used Einstein's summation convention, and ν being the Poisson ratio and C, D the stretching and bending moduli, respectively. To calculate the stretching modulus C, we use the results of Wei et al. [38] which, in our notation, are This value of ν agrees with the experimentally-measured Poisson ratio in graphite [39]. The value used for the bending modulus was that of Kudin et al. [40], calculated ab initio as D = 1.46 eV = 2.34 × 10 −19 N 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 two-dimensional 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 three-dimensional 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 three-dimensional Young's modulus E and the thickness h: where a typical value h ≈ 0.3 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 TPa [42,43,1]. 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 two-dimensional nature, we retain the parameters C and D as our main quantities here rather than express them in terms of Young's modulus E.

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: E elast , the stored elastic energy, and E 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 For now the energy density of external forces, V [w, v 1 , v 2 ; λ i ], is left unspecified. However, we do note that it is in this term that the influence of the control variables λ 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 1 , v 2 , and w. However, this requires some regularity in the behavior of the transverse displacement w: its first derivative must be continuous (C 1 ). Choosing C 1 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,46,47,48], which involves treating the moment tensor M αβ 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 Appendix A.2 for a detailed derivation): 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.

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(X). This is because we wish to find the deformation field that best approximates B(X) to a desired target, sayB(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 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(X i , X i + n) = t 0 + δt(X i , X i + n). The presence of δt, which is a relatively small perturbation to t 0 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 K = (4π/3 √ 3a, 0) in the first Brillouin zone has the form [15,16] where σ is a vector of Pauli matrices, and v F = 3t 0 a/2, with a = 1.42Å the carboncarbon distance in equilibrium. For deformations on scales that are large compared to a, the components of the pseudomagnetic vector potential A = A x e x + A y e y are explicitly given by [16] A where c = −∂ log t(r)/∂ log r| r=a (see Appendix B). For static deformations, a value c ≈ 3.37 captures the changes in various physical properties arising from straininduced modifications of the π bands in agreement with first-principles calculations [49,50,51,52,53,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 K leads to terms that introduce, for example, Fermi surface anisotropy [55,54] 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 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 pseudomagnetic field itself, B, being defined as the 2D curl of A, reads: As noted above, by virtue of our choice of piecewise affine finite elements for the numerical interpolation, the six variables v 1 , v 2 , w, M 11 , M 12 , and M 22 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 ε αβ to recover a strain field ε rec αβ of the same type as the primary variables: continuous and affine over each element. The derivative of ε rec αβ is well-defined, and thus so is the PMF if it is calculated using this recovered strain field.

Optimization
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 Eqs. (7a-7f). If we letB(X, Y ) be the desired PMF in Lagrangian coordinates, we then seek to minimize the functional to find a PMF, B, which is (ideally everywhere) as close as possible to the prescribed 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, I, a regularization term which penalizes high spatial variations in the control variables λ i : where η is a tunable parameter. The precise form of I reg [λ i ] will depend on what the control variables λ 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 1 , v 2 , w, M 11 , M 12 , M 22 and control variables λ i , the objective function (12) subject to the six equations (7a-7f), solved for all admissible variations (·). In these expressions, B is given by (10), N αβ is given by (2), and ε αβ 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 Tröltzsch [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).

Practical application: optimizing substrate shapes
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 =ẑ(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 PMFsB of typical scale B 0 = 10 T. We further assume that the domain Ω, representing the shape of the computational domain, has typical dimensions L = 100Å.

The external forces
As discussed above, the external potential term V [v 1 , v 2 , w; λ i ] will have two components: the work done due to hydrostatic pressure, pw, and the adhesion energy V 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 R 3 -space, such that a particle dS of the plate, located at ( Given a substrate shape z =ẑ(x, y), we could determine the adhesion potential V p at every point in three-dimensional 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 −ẑ(x, y)), where J(s) is some one-dimensional adhesion potential, such as the Lennard-Jones potential between surfaces [62], where s * is the adhesion well position, or the distance from the substrate at which a particle is in equilibrium. Thus, with the control variables λ i encoded inẑ. As a representative value for s * , we use the value that Xu and Buehler [63] give for C-Cu, namely 2.243Å. Similarly we use J 0 = 0.45 J m −2 as a representative value, from the investigation of Koenig et al. [7] into the adhesion strength between graphene and SiO 2 . We select a typical value for the hydrostatic pressure as p = 100 bar = 10 7 Pa.

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 α , M αβ ) 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 =ẑ(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 −ẑ(X + εv 1 , Y + εv 2 ) will not vary smoothly as v α 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 ξ 1 = y √ 3 + x and ξ 2 = y √ 3 − x the unit cell corresponds to the domain The topography of the substrate,ẑ(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 λ 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 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 nonzero -in particular this rules out the generation of strictly constant nonzero PMFs by periodic deformations. It should be emphasized, however, that our general procedure is applicable to arbitrary domains, geometries, § In this formulation the variations· are, for each nodal point i = 1, . . . , N p , the piecewise affine functions which take the value 1 at nodal point i and zero at each other nodal point. For a given weak form equation this provides N p equations for each of the N p unknown values of the function at the nodal points. 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 illbehaved during the numerical optimization, a regularization term, I 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 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: I 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.

Results and discussion
For illustration we initially seek PMFs with a typical strength of 10 T. The unit cell is chosen to have edges of length L = 10 −8 m and we apply a pressure of 100 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 graphenesubstrate 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, η = 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:   Results of the numerical optimization to make the PMF approach a prescribed spatial pattern with typical magnitude of 100 T. Column 1: the target PMF patterns (color scale in Tesla, on the left). Column 2: calculated PMF (left-hand color scale). Column 3: substrate topography associated with the PMF shown on its left (color scale in units of L = 10 nm, on the right). The unit cell is displayed in each image, and has edges of length L.
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.
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.
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) [33,6], a distribution of attached structures like nanotubes Three-dimensional plots of the optimized graphene topographies due to the substrate shapes in Figures 3 and 4, together with the resulting PMFs, are shown in Appendix F. [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 [30,26]. 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 [9,2]. 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 [3,2], 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.

Summary
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 pseudomagnetic field 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.

Appendix A. Graphene as an elastic continuum
Here we recapitulate the results of sections 2.1 and 2.2 with further discussion on their validity and applicability.

Appendix A.1. Definitions and assumptions
The deformation of a graphene sheet is described in terms of its deviation from a flat two-dimensional surface. A point 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 X = (X, Y, 0) is transformed to are the expressions for the new coordinates in terms of the in-plane displacements v α , and the vertical deflection w ¶. From these expressions one can define the base vectors r ,α = (x ,α , y ,α , z ,α ) and the metric tensor g αβ = r ,α · r ,β of the deformed surface. Then the true strain tensor is defined in terms of the difference of this metric tensor from its original value δ αβ (Kronecker's delta, or the identity tensor, from our choice of Cartesian coordinates): In terms of displacements, this becomes The second deformation measure of the surface is the curvature tensor, defined by ρ αβ = r ,αβ · n, where 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 Koiter [68] or Niordson [69]) but the full expression for the curvature tensor in Cartesian coordinates is where g = det g αβ . 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 ¶ Our convention is to have all Greek indices ∈ {1, 2}. The coordinate system of the undeformed sheet is chosen to be Cartesian, and we will make extensive use of Einstein's summation convention throughout this report. Subscripts following a comma denote partial differentiation with respect to that coordinate. 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 Denote the corrections to these approximations by ε corr αβ = ε true αβ − ε αβ , and similarly for ρ 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 The stress and moment resultants are assumed to be isotropic and linear in the strain and curvature tensors: ν 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 σ = 1 + D G /D. To calculate the stretching modulus C, we use the results of Wei et al. [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 Wei et al. [38] are, in our notation, This value of ν agrees with the experimentally-measured Poisson ratio in graphite [39]. The value we chose for the bending modulus was that of Kudin et al. [40], calculated ab initio as D = 1.46 eV = 2.34 × 10 −19 N m. We have found only two investigations into the value of D G (and hence σ) in graphene; the calculations of Wei et al. [71] lead to σ = −0.056, whereas the numerical study of Koskinen and Kit [72] gives a value of σ = 0.565. In the absence of consensus, in our calculations σ is set to be equal to the Poisson ratio: σ = ν = 0.169, and thus B αβγδ = A αβγδ . The constitutive equations for macroscopic materials are usually derived from full three-dimensional 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,74,75]. Using such analyses, σ = ν and the stiffness and bending moduli may be written in terms of the three-dimensional Young's modulus E and the thickness h: .
Using a typical value h ≈ 0.3 nm for graphene thickness [41], this has been used to cite graphene's Young's modulus as of the order of 1 TPa [42,43,1]. 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 twodimensional 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 (|ε αβ | < 0.1, |ρ αβ | < h −1 ) and that the corrections to the strains and curvatures are small.

Appendix A.2. Weak form equations
In section 2.2 we stated that equations (7a-7f) could be derived from a minimization of the energy integrals (6a-6b). In this section we justify this claim.
Recall that the stored energy in the plate was given by 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: Assume that the boundary of this plate is formed of three disjoint regions: ∂Ω = Γ c ∪ Γ s ∪ Γ f , with clamped conditions along Γ c , simply supported conditions along Γ s , and free conditions along Γ f . The standard variational approach is to minimize E t over the space of all twicedifferentiable w satisfying w = 0 on Γ c ∪ Γ s and ∂ n w = 0 on Γ c . The first variation of E t is δE t = Ω (DB αβγδ w ,αβw,γδ + pw) d 2 X, (A.15) wherew = δw is the variation in w. The weak solution is then the twice-differentiable function w(X, Y ) that satisfies δE t = 0 for each twice differentiable variationw satisfying w = 0 on Γ c ∪ Γ s and ∂ nw = 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: using the boundary conditions forw. Setting this to zero for each admissible variatioñ w, we find that the governing equation is with boundary conditions This is the Euler-Lagrange equation associated with the minimization of (A.14).
Appendix 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 αβ + . The variational functional is a version of the Hellinger-Reissner principle [44,47], given by Here is the inverse of B αβγδ . The weak form equations are derived from this principle by finding the stationary value of H over all functions M αβ and w satisfying the aforementioned conditions on Γ c and Γ s . Note that this stationary value of H will be neither a minimum nor a maximum; it is for this reason that these methods are often + Though Reinhart [46] and others have used the curvature tensor ρ αβ in place of M αβ , this merely results in a rearrangement of the governing equations, since one is a linear combination of the components of the other.
called saddlepoint methods. The weak form equations are found by setting the variation δH to zero, where 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 deflectioñ w; broadly speaking the square of its second derivative must be integrable. Meanwhile, the only regularity required of M αβ 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 This functional is one that can be minimized over the space of all w and M αβ 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 The weak form equations are obtained by finding the stationary value of (A.25) over the space of admissible functions satisfying w = 0 on Γ s ∪ Γ f ; in other words for all trial functionsw andM αβ 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 Arnold [44], Blum and Rannacher [47], and Oukit and Pierre [48].
Appendix A.2.2. Application to nonlinear plate bending The application of these mixed variational principles to nonlinear plates was first analyzed by Miyoshi [45] and Reinhart [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 αβ = CA αβγδ ε γδ for simplicity, where the mixed variational functional we use is According to the discussion above, in order to derive the weak form equations we should find the stationary value of H over all admissible v α , w, M αβ satisfying v α = w = 0 on Γ c ∪ Γ s . The first variation of H can be straightforwardly derived along the same lines discussed above, whereupon one obtains These equations then lead naturally to the weak form equations (7a-7f) on assuming that the entire boundary is clamped (Γ s = Γ f = ∅), that σ = ν, and on writing out the equations for the six componentsw,ṽ α andM αβ explicitly. The equations hold for all continuous integrable variationsṽ α ,w,M αβ that satisfỹ v α =w = 0 on the boundary.

Appendix B. Coupling deformations to electrons: PMFs
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(X). This is because we wish to find the deformation field that best approximates B(X) to a desired (target),B(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 π 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 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 secondquantized operator a X i (b 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 X i . The parameter t(X i , X i + n) is the hopping amplitude between two neighboring π orbitals, and n runs over the three unit cells containing a B atom neighboring the A atom from the unit cell at X i . The hopping amplitude is constant in the perfect crystal: t(X i , X i + n) = 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(X i , X i + n) = t 0 + δt(X i , X i + n). The presence of δt, which is a relatively small perturbation to t 0 in practical situations, adds a correction to the low-energy Diraclike Hamiltonian that emerges from (B.1) so that the effective Hamiltonian around the point K = (4π/3 √ 3a, 0) in the first Brillouin zone has the form [15,16] where σ is a vector of Pauli matrices, and v F = 3t 0 a/2, with a = 1.42Å the carboncarbon 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 A = A x e x + A y e y by expanding t to linear order in the deformation tensor. Choosing the coordinate system so that e x is along the zig-zag direction of the honeycomb lattice one obtains [16] A where c = −∂ log t(r)/∂ log r| r=a . For static deformations, a value c ≈ 3.37 captures the changes in various physical properties arising from strain-induced modifications of the π bands in agreement with first-principles calculations [49,50,51,52,53,54]. Finally, the pseudomagnetic field B, being defined as the 2D curl of A, reads: (B.4) * Note, however, that this is not a restriction on the applicability of the method. The assumption of small deviations from the planar configuration is for convenience and definiteness only. A full parametrization of the hopping modifications including curvature-induced re-hybridization would be dealt with in precisely the same way, because the only ingredient that is needed is the dependence of the PMF B on the strain components. The central and only requirement is the ability to explicitly specify this dependence, as done in equation (10) under the stated conditions.
As noted in the previous section, by virtue of our choice of piecewise affine finite elements for the numerical interpolation, the six variables v 1 , v 2 , w, M 11 , M 12 , and M 22 are treated as continuous, and affine over each triangular element. As a consequence, from (A.5) the strain field will be discontinuous and constant in each triangular element. As the strain components ε αβ 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 Zienkiewicz and Zhu [58], we refer the reader to Appendix E. In brief, this is a mechanism that uses the discontinuous strain data ε αβ to recover a strain field ε rec αβ of the same type as the primary variables: continuous and affine over each element. The derivative of ε rec αβ is well-defined, and thus so is the PMF if it is calculated using this recovered strain field:

Appendix C. Nondimensionalization
In solving the optimization problem of section 3 numerically, the first step is to nondimensionalize 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 to be the typical scaling of the strain field, then All these constants are previously-defined, with the exception of L, representing the typical size of the domain Ω, and B 0 , the typical magnitude of the target PMFB. Under these scalings the equations exhibit only one dimensionless parameter, namely the dimensionless bending stiffness κ: For completeness we will summarize the minimization problem in its dimensionless form:  and the three dimensionless parameters are given bȳ As a representative value for s * , we use the value that Xu and Buehler [63] give for C-Cu, namely 2.243Å. Similarly we use J 0 = 0.45 J m −2 as a representative value, from the investigation of Koenig et al. [7] into the adhesion strength between graphene and SiO 2 . We select a typical value for the hydrostatic pressure as p = 100 bar = 10 7 Pa.
The derivatives of the potentials appearing in the dimensionless weak form equations arē (C. 23) In section 4.2 and subsequently, all variables are assumed to be dimensionless, and overbars are omitted for clarity.

Appendix D. Parametrization of the substrate topography
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 the unit cell corresponds to the domain The topography of the substrate,ẑ(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: and so we choose a truncated expansion in terms of these as follows: with the constants α kl , β kl , γ kl , δ kl playing the role of the control variables λ i : varying the topography of the substrate is, therefore, achieved by varying these 4(K + 1) 2 constants. For convenience we set for each k and, since rigid vertical displacements of the substrate do not affect the objective function, we further set α 00 = 0. 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, I reg , is added to the objective integral, as per equation (12). We choose it to be I reg = 1 Area(Ω) Ω |∇ẑ| 2 d 2 X which is simple to calculate using the orthogonality of the basis functions f n 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 two-dimensional periodicity of the substrate. Due to the geometry of the unit rhombus, we can set the triangulation to be a regular isometric grid.

Appendix E. Strain recovery
As noted in Appendix A, we choose a finite element discretization for our six variables v 1 , v 2 , w, M 11 , M 12 , M 22 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).
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 displacementbased finite element method. In this article a piecewise affine displacement field leads to a piecewise constant approximation ε αβ to the strain field, whereas we require a differentiable approximation. By using the patch recovery method we recover a piecewise affine strain field, ε 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, . . . , N t and nodes i = 1, . . . , N p . We have a function f , constant on each element (so f (X, Y ) = f k if X is in triangle k), from which we want to recover a piecewise affine function f rec (defined by its values f rec i at each node X i of the triangulation). The strength of the patch recovery method is that the nodal values f rec i 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 fit i (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 fit i (X, Y ) = a + bX + cY to the patch. We use the values f k , evaluated at the centroids (X c k , Y c k ) of the elements, to calculate the parameters (a, b, c) = a T through a least-squares optimization. As noted in Zienkiewicz and Zhu [58], a is thus found by solving the system   k∈patch(i) where p k = (1, X c k , Y c k ) T . Having found f fit i (X, Y ), the nodal value of the recovered function f rec is simply f rec i = f fit i (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 Zienkiewicz and Zhu [58] and find the boundary nodal values of f 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.
Appendix F. Three-dimensional plots Figure F1 shows three-dimensional visualizations of deformed graphene sheets corresponding to the four solutions of Figure 3, where the target PMF value was 10 T. Vertical scales are exaggerated for clarity.
The corresponding visualizations for B = 100 T (corresponding to Figure 4) are shown in Figure F2. In this case the vertical scale is not exaggerated. Figure F1. Three-dimensional plots of graphene sheets deformed by the four substrates in Figure 3, colored with the resultant pseudomagnetic fields. Vertical scales in some of the plots are exaggerated for clarity, by factors of 2, 3, 1 and 3 respectively.