Commensurate-incommensurate phase transition and a network of domain walls in bilayer graphene with a biaxially stretched layer
Irina V. Lebedeva Andrey M. Popov




liv_ira@hotmail.com
Nano-Bio Spectroscopy Group and ETSF, Universidad del País Vasco, CFM CSIC-UPV/EHU, 20018 San Sebastián, Spain
[Present address: ]CIC nanoGUNE,
San Sebastián 20018, Spain

popov-isan@mail.ru
Institute for Spectroscopy of Russian Academy of Sciences, Troitsk, Moscow 108840, Russia



The two-chain Frenkel-Kontorova model is applied for an analytical description of the energy and structure of the network of domain walls in bilayer graphene. Using this approach, the commensurate-incommensurate phase transition upon biaxial stretching of one of the graphene layers is considered. We demonstrate that formation of the equilateral triangular network of domain walls becomes energetically favourable above the critical relative biaxial elongation of the bottom layer of 3.0· 10^-3. It is shown that the optimal period of the triangular network of domain walls is inversely proportional to the difference between the biaxial elongation of the bottom layer and the critical elongation as long as it is much greater than the width of domain walls. Quantitative estimates of the contribution of a single dislocation node to the system energy and the period of the network of domain walls are obtained. Experimental measurements of the period could help to verify the energy of the fully incommensurate state (such as obtained by relative rotation of the layers) with respect to the commensurate one. 

Commensurate-incommensurate phase transition and a network of domain walls in bilayer graphene with a biaxially stretched layer
    Andrey M. Popov
    Received: date / Accepted: date
===============================================================================================================================





§ INTRODUCTION

Stacking dislocations arising from variation in stacking of graphene layers were recently predicted for bilayer graphene <cit.>. Since then, characteristic patterns consisting of commensurate domains separated by stacking dislocations, which play the role of incommensurate domain walls, have been observed in numerous experiments <cit.>. It has been demonstrated that stacking dislocations affect electronic <cit.> and optical <cit.> response of graphene and, similar to in-plane defects, they should be taken into account upon development of nanoelectronic and nanoelectromechanical devices.

In stacking dislocations, one of the layers in a bilayer is slightly stretched and the other one slightly compressed and/or there is a different shear strain in the layers so that the stacking changes from the one ground-state stacking to another one <cit.>. The variation of the stacking is mostly localized in narrow strips with the width much smaller than the size of commensurate domains where the stacking is close to the ground-state one. These narrow strips are referred to as domain walls <cit.> or boundaries between commensurate domains <cit.>. 

In previous theoretical works <cit.>, it was assumed that domain walls do not cross as long as the distance between them is large and can be treated as isolated. It was shown that formation of domain walls becomes energetically favourable above some critical uniaxial elongation of the bottom layer <cit.>. At small elongations, the layers are maintained commensurate by the interlayer interaction. However, at some critical elongation, the interlayer interaction energy can no longer compensate the large elastic energy. As a result, formation of the first domain wall takes place. The distance between adjacent domain walls decreases continuously as the elongation of the bottom layer is increased. Therefore, there is a commensurate-incommensurate phase transition of the second order, which can be characterized by the inverse distance between domain walls as the order parameter <cit.>. So far such a phase transition has been studied for bilayer graphene  <cit.>, bilayer boron nitride <cit.> and graphene-boron nitride heterostructure <cit.> only upon uniaxial stretching and crossing of domain walls was neglected.

Nevertheless, some experimental images of bilayer graphene obtained by transmission electron microscopy <cit.> and near-field infrared nanoscopy <cit.> revealed the presence of  triangular networks of intercrossed domain walls. Such networks are associated with global interlayer biaxial and rotation strains, i.e. different strains within the graphene layers and relative rotation of the layers <cit.>. Therefore, deformation of one of the layers and/or its rotation with respect to the other should promote formation of networks of domain walls. In the present paper, we develop an analytical approach for description of networks of domain walls based on the two-chain Frenkel-Kontorova <cit.>. While diverse patterns of domain walls and loading conditions can be studied using this model, here we limit our consideration to the case when one graphene layer is biaxially stretched with respect to the other. The application of a tensile strain favors formation of tensile domain walls perpendicular to the Burgers vector <cit.>. Correspondingly, a regular network of tensile domain walls separating commensurate domains with the shape of equilateral triangles (Fig. <ref>) should arise upon application of a biaxial strain. We show that similar to uniaxial stretching, the commensurate-incommensurate phase transition takes place upon increasing the biaxial elongation of one of the graphene layers.  It should be noted that the commensurate-incommensurate phase transition considered here is very different from the crossover from the incommensurate state with commensurate domains separated by domain walls to the fully incommensurate state observed for graphene on boron nitride upon rotation of the graphene layer <cit.>. 





In the following, we review the theory of domain walls in bilayer graphene, develop the model for bilayer graphene with a biaxailly stretched bottom layer and apply it to estimate the critical biaxial elongation and the period of the network of domain walls. Finally we discuss the accuracy of our model and summarize our conclusions.




§ STRUCTURE AND ENERGY OF DOMAIN WALLS

To understand the structure of domain walls in bilayer graphene and the reasons for formation of the network of domain walls, it is necessary to consider the potential surface of interlayer interaction energy of graphene layers, i.e. the dependence of the interlayer interaction energy on the relative in-plane displacement of the layers. 
This potential energy surface is determined by hexagonal symmetry of a single graphene layer  <cit.> (Fig. <ref>). The minima of the potential energy surface are located in vertices of hexagons, the maxima in the centers of the hexagons and the barriers for transition between the minima in the middle of the sides of the hexagons. The maxima correspond to the AA stacking in which all atoms of one layer are on top of the atoms of the second one. The minima are observed for the AB stacking obtained from the AA stacking by the shift one of the layers by one bond length l in any armchair direction. The barriers for transition between the minima correspond to the saddle-point (SP) stacking, which can be obtained from the AA stacking by the shift of one of the layers by 3l/2 in any armchair direction. 

The structure and energy of domain walls depend only on changes of the interlayer interaction energy for different stackings and not on its absolute values. Therefore, in Fig. <ref> and below, we consider the interlayer interaction energy relative to the energy of the AB stacking, i.e. set it to zero at the minima of the potential energy surface.    

Domain walls in bilayer graphene correspond to the relative displacement of the layers in the armchair direction through the SP stacking <cit.>, which is the minimum energy path between adjacent minima AB. The dislocations are partial since the Burgers vector b⃗ equal in magnitude to the bond length l is smaller than the lattice constant a_0 = l √(3). In the bilayer with a biaxially stretched bottom layer, a triangular network of domain walls should be formed with six domain walls corresponding to six armchair directions merging at each dislocation  node (Fig. <ref>). It is clear, therefore, that dislocation nodes correspond to the transition through the AA stacking (Fig. <ref>). 









To study domain walls in bilayer graphene we use the two-chain Frenkel-Kontorova model  <cit.>, which takes into account structural relaxation of both of the layers and was applied previously to study domain walls in double-walled carbon nanotubes <cit.>,  bilayer graphene  <cit.>,  bilayer boron nitride <cit.> and graphene-boron nitride heterostructure <cit.>. The out-of-plane buckling of graphene layers <cit.> is neglected. This is justified when the bilayer is supported, for example, on misoriented graphene layers <cit.>, boron nitride <cit.>, etc. In such cases, the two-chain Frenkel-Kontorova model parameterized on the basis of first-principles calculations <cit.> gives the width of domain walls in close agreement with the experimental data <cit.>. The account of out-of-plane buckling in suspended graphene can be performed using more advanced continuum models <cit.> or atomistic approaches <cit.>.

To describe a straight and isolated domain wall in free bilayer graphene within the two-chain Frenkel-Kontorova model, we consider two chains of particles connected by harmonic springs (Fig. <ref>). Each particle corresponds to a ribbon of graphene parallel to the domain wall, of unit length and of width equal to the bond length l. The number of particles in each chain N≫ 1. The elastic constant of the springs is k_s. The equilibrium distance between particles in isolated chains is l. Then the elastic energy of the chains can be computed as 

    W_el =    k_s l^2/2∑_n=-N/2^N/2{(v_n-v_n-1)^2 
        + (v_n+u_n-v_n-1-u_n-1)^2},

where v_n is the displacement of the n-th particle  in the chain corresponding to the bottom layer relative to the position p_n=nl in the isolated chain and u_n is the relative displacement of the n-th particles in the chains corresponding the upper and bottom layers. Both of these quantities are in the units of the bond length l.

The interaction energy of the chains is given by

    W_int   =  ∑_n=-N/2^N/2V(u_n)l,

where V(u) is the interlayer interaction energy of bilayer graphene per unit area relative to that in the commensurate state with the AB stacking. 

In the continuum approximation and in the limit N→∞, the total energy of the system, W = W_el+W_int, can be written as

    W =     ∫_-∞^+∞{k_s l^2/2(dv/dn)^2  + k_s l^2/2(dv/dn+du/dn)^2
       +V(u)l}dn.


Introducing the variable a = v + u/2, the above equation can be presented in the form

    W =    ∫_-∞^+∞{ k_s l^2(da/dn)^2 + k_s l^2/4(du/dn)^2
       +V(u)l}dn.


The energy described by Eq. (<ref>) is minimized when da/dn is set to zero and 0 ≤ u(n) ≤ 1 corresponds to the relative displacement of atoms of the layers between two adjacent AB minima through the SP stacking (Figs. <ref> and <ref>). In the present paper we limit our consideration to tensile domain walls. In this case, i.e. when the Burgers vector is perpendicular to the domain wall, k_sl=k/(1-ν^2), where k is the elastic constant of graphene under uniaxial stress (related to the Young's modulus Y and the thickness of graphene layers h as k=Yh) and ν is the Poisson's ratio. Then Eq. (<ref>) is reduced to

    Δ W =    ∫_-∞^+∞{k l^2/4(1-ν^2)|du/dx|^2  + V(u)}dx,

where Δ W=W  is the energy associated with creation of a single domain wall per unit length and coordinate x=nl corresponds to the direction perpendicular to the domain wall and parallel to the Burgers vector (i.e. along the armchair direction).

The relative displacement u(x) that minimizes the formation energy of domain walls in Eq. (<ref>) is determined by the Euler-Lagrange equation δΔ W/ δ u=0, which upon integration gives

    k l^2/4(1-ν^2)|du/dx|^2 = V(u).


It was shown in the previous papers <cit.> that the potential energy surface of bilayer graphene can be described with high accuracy by the expression containing only the first Fourier harmonics:

    V(u_x,u_y) =    V_0(3/2+cos(2k_0u_x-2π/3) 
       -2cos(k_0u_x-π/3)cos(k_0u_y√(3))),
 
where k_0 = 2π/3 and u_x and u_y correspond to relative displacements of the layers in the armchair and zigzag directions, respectively, in units of the bond length l. Here, as mentioned above, the interlayer interaction energy is given relative to that in the AB stacking.  

For u_y=0 and 0 ≤ u_x=u(x) ≤ 1, the dependence of the interlayer interaction energy along the straight line between two adjacent minima AB, which also corresponds to the minimum energy path, can be written as

    V(u) = V_max(2cos(k_0 u + 2π/3) +1)^2,

where V_max=V_0/2 is the barrier to relative sliding of graphene layers through the SP stacking. 

From Eqs. (<ref>) and (<ref>), it follows that the domain wall is described by a soliton:

    u(x)    = 1/2+
       3/πarctan[1/√(3)tanh(2π x √(V_max(1-ν^2)/3kl^2))].


The dependence u(x) has a nearly constant slope close to the center of the domain wall at x=0, which is roughly equal to |d u/d x|_x=0.
Therefore, the dislocation width corresponding to the characteristic width of domain walls can be defined as <cit.>

    l_D = |d u/d x|^-1_x=0 =l/2√(k/V_max(1-ν^2)).


The dislocation energy per unit length of domain walls based on Eqs. (<ref>), (<ref>) and (<ref>) becomes

    Δ W = √(k l^2/(1-ν^2))   ∫_0^1 √(V(u))du
       =√(k l^2 V_max/(1-ν^2))(3√(3)/π-1).


The use of the continuum approximation in Eqs. (<ref>) – (<ref>) is justified if l_D≫ l. In the present paper, we use the following parameters for the model obtained by density functional theory (DFT) calculations <cit.> using the the second version of the van der
Waals density functional (vdW-DF2) <cit.>: l = 1.430 Å, k = 331 ± 1 J/m^2, ν = 0.174 ± 0.002 and V_max = 1.61 meV/atom (in meV per atom of the upper/adsorbed layer). The latter value for the barrier is in agreement with the data deduced from shear mode frequencies of bilayer and few-layer graphene and graphite <cit.> and within the range of other DFT results <cit.>. Using these parameters, we estimate that for tensile domain walls in bilayer graphene, l_D =13.4 nm and this is indeed much greater than the bond length. This result is close to the corresponding experimental data of 10.1 ± 1.4 nm (Ref. Alden2013) and 11 nm (Ref. Yankowitz2014).

Let us now discuss the effect of uniaxial stretching of the bottom layer on the formation energy of domain walls. For that we consider relative elongations of the layers in the case when one domain wall is created in the system of finite length L in the direction of the Burgers vector <cit.>.  
According to the two-chain Frenkel-Kontorova model, atoms of both of the layers are displaced relative to their positions in the initial commensurate system upon creation of a domain wall <cit.>. One domain wall corresponds to the relative displacement of the layers equal to the bond length l. Therefore, if the layers are free and composed of the same material, the resulting relative elongations of the bottom and upper layers with one domain wall are ϵ_0= l/(2L) and -ϵ_0, respectively. This means that Eqs. (<ref>) and (<ref>) describe the energy of the finite bilayer with one domain wall and relative elongation of the bottom layer ϵ_0 with respect to the commensurate system (the bilayer with no domain walls) with zero elongation. Analogously, these equations also describe the energy of the bilayer with one domain wall and relative elongation of the bottom layer ϵ with respect to the commensurate system with the relative elongation ϵ-ϵ_0. Thus, to calculate the energy of the bilayer with one domain wall  with respect to the commensurate system with the same relative elongation ϵ of the bottom layer, it is needed to substract from Eqs. (<ref>) and (<ref>) the elastic energy coming from the increase of the relative elongation of the commensurate bilayer from ϵ-ϵ_0 to ϵ. While the formation energy of domain walls in the case when the relative elongation of the bottom layer is greater by the extra elongation ϵ_0 compared to the commensurate bilayer (see Eqs. (<ref>) and (<ref>)) is always positive, the formation energy of domain walls in the system with the fixed length of the bottom layer can be negative for sufficiently large ϵ (above the critical value) <cit.>.



§ TRIANGULAR NETWORK OF DOMAIN WALLS

For a biaxially stretched bottom layer, we consider the triangular network with domain walls in all zigzag directions (Fig. <ref>). We assume that commensurate domains are large compared to the width of domain walls so that the side of the triangles L_T≫ l_D. In this case, Eqs. (<ref>) and (<ref>) describe the energy of domain walls in the bilayer with the relative biaxial elongation of the bottom layer ϵ with respect to the commensurate system with the relative  biaxial elongation ϵ-ϵ_0, where ϵ_0=√(3)l/(2L_T). Correspondingly, the energy of the bilayer with the network of domain walls relative to the commensurate bilayer at the same relative elongation ϵ of the bottom layer can be written as

    Δ W_T = -Δ W_el + Δ W_dw+ Δ W_dn,

where Δ W_el is the increase of the elastic energy of the commensurate bilayer due to the extra elongation ϵ_0, Δ W_dw and Δ W_dn are the contributions of domain walls (edges of triangles) and dislocation nodes (vertices of triangles, see Fig. <ref>) into the formation energy of the network of domain walls, respectively. Below we give these relative energies per unit area of the overlap between the layers of the bilayer. 

The increase of the elastic energy of the commensurate bilayer due to the extra elongation ϵ_0 is simply given by

    Δ W_el =   2k/(1-ν)(ϵ^2-(ϵ-ϵ_0)^2)
       =2√(3)kϵ/(1-ν)l/L_T-3k/2(1-ν)(l/L_T)^2.


The contribution of domain walls is Δ W_dw=3L_TΔ W/(2S_T), where Δ W is the energy per unit length of domain walls given by Eq. (<ref>) and S_T=√(3)L_T^2/4 is the area of one commensurate domain. Therefore,

    Δ W_dw = 2l/L_T√(3k V_max/(1-ν^2))(3√(3)/π-1).


The contribution of dislocation nodes should be positive and proportional to the density of nodes. It, therefore, should depend quadratically on the period of the network of domain walls Δ W_dn∝(l/L_T)^2. 

From this assumption and Eqs. (<ref>)–(<ref>), the expression for the energy of the bilayer with the  triangular network of domain walls relative to the commensurate system can be written as 

    Δ W_T= Al/L_T +Bl^2/L^2_T,

where 

    A=-2√(3)(kϵ/(1-ν)-√(k V_max/(1-ν^2))(3√(3)/π-1))

and B is determined by the contribution Δ W_dn of dislocation nodes  and the term of the elastic energy Δ W_el  quadratic in l/L_T (see Eq. (<ref>)). Both of these contributions to the total energy are positive and B is positive.

It should be noted that limit L_T→∞ corresponds to the commensurate system. Accordingly, the energy Δ W_T of the bilayer with the network of domain walls with respect to the energy of the commensurate system given by Eq. (<ref>) in this limit tends to zero.

The optimal period of the network of domain walls can be found from the conditions ∂Δ W_T/∂ L_T=0 and ∂^2 Δ W_T/∂ L_T^2≥0. It follows from these conditions that when A is positive, the optimal period of the network tends to infinity,  i.e. the ground state corresponds to the commensurate system. When A is negative, the ground state corresponds to the bilayer with the network of domain walls characterized by the period

    L_0 = -2B/Al.


From Eqs. (<ref>) and (<ref>), it can be checked that in the latter case, the relative energy of the bilayer with the optimal network of domain walls is given by 

    Δ W_0 = -A^2/4B.

Since B is positive (see Eqs. (<ref>)), the structure with the triangular network of domain walls is indeed more energetically favourable than the commensurate system for negative A.

The critical relative biaxial elongation ϵ_c above which the network of domain walls corresponds to the ground state of the system is determined by the condition A=0. Using Eq. (<ref>), this gives

    ϵ_c=√(V_max (1-ν)/k(1+ν))(3√(3)/π-1)=(1-ν)Δ W/kl.

Note that this critical elongation ϵ_c does not depend on the contribution  Δ W_dn of dislocation nodes to the relative energy of the bilayer with the network of domain walls having a quadratic dependence on l/L_T.


Using the parameters obtained by the DFT calculations <cit.> (see Sec. II), we estimate that the triangular network of domain walls corresponds to the ground state of bilayer graphene with a biaxially stretched bottom layer at relative elongations  ϵ > ϵ_c=3.0· 10^-3. For the bilayer with a uniaxially stretched layer, the critical elongation for becoming domain walls energetically favourable is given by ϵ_c0=Δ W/kl and is about 3.6· 10^-3 (Ref. Lebedeva2016). It is seen that the ratio of the critical biaxial and uniaxial elongations is determined by the Poisson's ratio: ϵ_c/ϵ_c0=1-ν, i.e. the critical value for biaxial stretching is smaller than in the case of uniaxial stretching by only 20%. This may seem surprising as the elastic energy is about twice greater upon introduction of the biaxial elongation compared to the same uniaxial elongation. The density of dislocations, however, is also different in two cases and the critical elongation turns out to be similar. 

From Eq. (<ref>), it follows that above the critical elongation, the optimal period L_0 of the network of domain walls is proportional to L_0∝(ϵ-ϵ_c)^-1.  The inverse quantity L_0^-1 can be considered as the order parameter of the phase transition. 
From Eq. (<ref>), it can be concluded that  the relative energy of the bilayer with the optimal network of domain walls changes as Δ W_0 ∝ -l^2/L_0^2 ∝-(ϵ-ϵ_c)^2 above the critical elongation. This demonstrates that the commensurate-incommensurate phase transition in bilayer graphene with a biaxially stretched bottom layer is of the second order. 

The predicted dependence of the optimal period L_0 of the network of domain walls L_0∝(ϵ-ϵ_c)^-1 is different from the logarithmic law  L_0∝ln[(ϵ-ϵ_c)^-1] for the optimal separation between adjacent domain walls in the standard Frenkel-Kontorova model for a one-dimensional chain on a fixed substrate <cit.>. This difference is related to the contributions to the energy of the system with domain walls that have a non-linear dependence on the inverse distance between the domain walls. In our case, such contributions come from the extra elongation of the bottom layer that is needed to fulfil the boundary conditions (see Eq. (<ref>)) and the presence of dislocation nodes. Both of these contributions are positive and depend quadratically the period of the network of domain walls. In the standard Frenkel-Kontorova model these contributions are absent. Instead the repulsive interaction between the domain walls that depends exponentially on their separation is taken into account <cit.>. We assume that this exponential energy term can be neglected as compared to the quadratic terms for large periods of the network of domain walls, L ≫ l_D. 

Note that the quadratic contribution to the energy should be present in any two-layer system where there are boundary conditions for the substrate layer but it is not completely fixed (even in one-dimensional systems, bilayers with a uniaxially stretched bottom layer and no dislocation nodes <cit.>, etc.). Therefore, in such systems, the separation between domain walls should also change inversely proportional to the difference between the elongation of the substrate layer and critical elongation upon the commensurate-incommensurate phase transition. 

To get quantitative estimates of the optimal period of the network of domain walls, we assume that the nodes have the shapes of hexagons with the side l_D corresponding to the width of domain walls (Fig. <ref>). We also suppose that the layers within the nodes are uniformly stretched with the biaxial strain ± l/(2l_D). Accordingly, the relative displacement between the layers in these regions changes linearly and corresponds to the AB stacking at the vertices of the dislocation nodes and AA stacking at the centers (Fig. <ref>). This means that the average interlayer interaction energy within the nodes is given by 

    V_in = ∫_hex V(u_x,u_y)du_xdu_y/∫_hexdu_xdu_y=3/2V_0=3V_max,

where the integration is performed over one hexagon on the potential energy surface with vertices at the AB stacking and center at the AA stacking  (Fig. <ref>) and we use Eq. (<ref>) for the interlayer interaction energy V(u_x,u_y). It is clear that V_in is equal simply to the average of the interlayer interaction energy over the whole potential energy surface. Therefore, in our model, the structure of the layers within the dislocation nodes can be referred to as fully incommensurate. A similar state is achieved when the layers are rotated with respect to each other by an arbitrary angle that does not correspond to any Moiré pattern <cit.> (and even in the structures corresponding to Moiré patterns, the average interlayer interaction energy is only slightly different from V_in as shown in Ref. Yakobson).

Considering that the layers are uniformly stretched within the nodes and the interlayer interaction energy is equal to its average over the potential energy surface, the contribution of dislocation nodes to the relative energy of the bilayer with domain walls can be written as 

    Δ W_dn = 3(l_D/L_T)^2(V_in + k/2(1-ν)(l/l_D)^2).


Taking into account Eqs. (<ref>) and (<ref>), this contribution can be presented as

    Δ W_dn = (l/L_T)^2 3(5+2ν)k/4(1-ν^2).


Then the parameter B in Eq. (<ref>) takes the form

    B= 3(7+4ν)k/4(1-ν^2).


The dependence of the optimal period L_0 of the network of domain walls on the biaxial elongation estimated using Eq. (<ref>) with the parameter B from Eq. (<ref>)  is shown in Fig. <ref>. The contribution of a single dislocation node to the relative energy of the bilayer with the triangular network of domain walls is w_dn=Δ W_dn√(3)L_T^2/2=151 eV.  This is the same as the contribution of a domain wall of the length of 0.14  μm.



Our estimates are valid as long as L_0 ≫ l_D =13.4 nm. From  Eqs. (<ref>), (<ref>), (<ref>) and (<ref>)  it follows that this condition is satisfied for 

    ϵ≪ϵ_c+√(3)/2(7+4ν)√(V_max (1-ν)/k(1+ν)),

which corresponds to ϵ≪ 3.3· 10^-2.



§ ACCURACY OF THE MODEL

Let us now discuss the accuracy of our model for dislocation nodes in which we assume that the layers are just uniformly stretched and fully incommensurate. A similar approximation can be also considered for domain walls, which formation energy is known from the exact analytical solution (see Eq. (<ref>)). For simple estimates, we can suppose that the layers are completely commensurate within the commensurate domains, while in a domain wall, the layers are uniformly stretched with the uniaxial tensile strain ± l/(2l_D) and the relative displacement between them changes linearly across the domain wall and lies on the straight line between two adjacent minima AB on the potential energy surface (Fig. <ref>). In this case, the formation energy of domain walls given by Eq. (<ref>) can be written as

    Δ W =  k l^2/4(1-ν^2)l_D  + V_avl_D,

where V_av is the average interlayer interaction energy along the line between two adjacent minima. 

From Eq. (<ref>), it follows that

    V_av =  ∫_0^1 V(u)du= V_max(3-9√(3)/2π).


Optimization of the dislocation width l_D within this simple model using Eq. (<ref>) and the condition ∂Δ W/∂ l_D=0, gives 

    l_D=l/2√((3-9√(3)/2π)^-1k/V_max(1-ν^2))

and 

    Δ W =√((3-9√(3)/2π)k l^2 V_max/(1-ν^2)).

The latter value exceeds the result of the two-chain Frenkel-Kontorova model corresponding Eq.  (<ref>) only by 10%. 

The above value is the minimal error of the model in which the layers are uniformly stretched and incommensurate within domain walls. However, in the case of dislocation nodes formed by crossing domain walls, their size is determined by the dislocation width l_D from Eq.  (<ref>). Using a similar assumption for domain walls, gives the formation energy 

    Δ W =1/2(4-9√(3)/2π)√(k l^2 V_max/(1-ν^2)).

This is 16% greater than the result given by Eq.  (<ref>). Still this error is acceptable given that the barrier to relative sliding of the layers is not known with a high precision (the DFT values for the barrier lie in the range of 0.5 – 2.1 meV/atom <cit.>, although the interval can be reduced if the experimental interlayer distance is used <cit.>, and the estimates of the barrier from 
the experimental measurements of the shear mode frequencies <cit.> and dislocation width <cit.> correspond to 1.7 meV/atom and 2.4 meV/atom, respectively). It can be expected that our model used to analyze the triangular network of domain walls in bilayer graphene with a biaxially stretched bottom layer has a similar accuracy. To improve the accuracy of the model, a non-uniform strain distribution within dislocation nodes should be taken into account.




§ CONCLUSIONS AND DISCUSSION

The two-chain Frenkel-Kontorova model has been applied to analyze the commensurate-incommensurate phase transition in bilayer graphene with a biaxially stretched bottom layer. It is found that the critical biaxial elongation above which the triangular network of domain walls corresponds to the ground state of the system does not depend on the structure and energetics of dislocation nodes. The critical biaxial elongation is estimated to be about 3.0· 10^-3 and its difference with the critical elongation in the case of the uniaxially stretched bottom layer is determined purely by the Poisson's ratio. It is also predicted that above the critical elongation, the optimal period of the network of domain walls changes inversely proportional to the difference between the biaxial elongation of the bottom layer and the critical elongation. Considering the inverse period of the network as the order parameter, it is proved that the commensurate-incommensurate phase transition is of the second order. 

To get quantitative estimates of the period of the triangular network of domain walls and energy of dislocation nodes, it is assumed that in the nodes, the layers are uniformly stretched and fully incommensurate. The estimated contribution of a single dislocation node to the total energy of the system is found to be about 150 eV. This quantity can be overestimated in our model of dislocation nodes by no more than 20%, as demonstrated by the example of the formation energy of domain walls. 

In our estimates, we use the approximation of the potential energy surface of bilayer graphene by the first Fourier harmonics. This approximation has been justified by a number of DFT calculations <cit.>. It was also applied to estimate the barrier V_max to relative sliding of graphene layers from measurements of the shear mode frequency of graphene layers in bilayer and few-layer graphene and graphite and the value of 1.7 meV/atom was deduced <cit.>, similar to the most reliable DFT results <cit.>. The close value of 2.4 meV/atom was obtained from measurements of the dislocation width in few-layer graphene <cit.>. 
Experimental measurements of the period of the triangular network of domain walls, e.g., by  transmission electron microscopy <cit.>, scanning tunneling microscopy <cit.> or near-field infrared nanoscopy <cit.>, would allow to get an experimental estimate for the energy of transition of bilayer graphene from the commensurate state to the fully incommensurate one and to further check the accuracy of the approximation of the potential energy surface.

The approach for description of networks of domain walls in bilayer graphene developed in the present paper on the basis of the two-chain Frenkel-Kontorova model  has allowed us to get an analytical solution for the case of biaxial elongation of the bottom layer, i.e. when the commensurate domains have the shape of equilateral triangles and all the domain walls are tensile. 
In the experiments  <cit.>, the triangular domains are not always equilateral and the domain walls are not equivalent. The corresponding  images can be obtained when strains applied to the bottom layer are not equal along different axes and/or there is a shear strain applied. Relative rotation of the layers and bending of the bilayer also affect the local structure of domain walls and the global structure of the network of domain walls. The current analytical approach based on the two-chain Frenkel-Kontorova model can be straightforwardly extended to describe formation of networks of domain walls for such types of the external load and they will be considered in subsequent papers.  

Other diverse patterns of domain walls have been observed experimentally including irregular triangular networks with curved domain walls  <cit.>, L-shape domain walls and closed-loop circles  <cit.>, etc. A non-uniform distribution of strains in the layers, presence of defects in graphene layers and substrate as well as existence of barriers to formation and transformation of domain walls <cit.> can be considered as possible reasons for formation of such patterns and require further investigation.



§ ACKNOWLEDGMENTS

AMP acknowledges the Russian Foundation for Basic Research (Grant 18-02-00985).


