Discrete breathers in honeycomb Fermi–Pasta–Ulam lattices

We consider the two-dimensional Fermi–Pasta–Ulam lattice with hexagonal honeycomb symmetry, which is a Hamiltonian system describing the evolution of a scalar-valued quantity subject to nearest neighbour interactions. Using multiple-scale analysis we reduce the governing lattice equations to a nonlinear Schrödinger equation coupled to a second equation for an accompanying slow mode. Two cases in which the latter equation can be solved and so the system decoupled are considered in more detail: firstly, in the case of a symmetric potential, we derive the form of moving breathers. We find an ellipticity criterion for the wavenumbers of the carrier wave, together with asymptotic estimates for the breather energy. The minimum energy threshold depends on the wavenumber of the breather. We find that this threshold is locally maximized by stationary breathers. Secondly, for an asymmetric potential we find stationary breathers, which, even with a quadratic nonlinearity generate no second harmonic component in the breather. Plots of all our findings show clear hexagonal symmetry as we would expect from our lattice structure. Finally, we compare the properties of stationary breathers in the square, triangular and honeycomb lattices.

et al [23,24], where there is a vector-valued function at each node, the in-plane, horizontal and vertical displacements. In [6] the lattice considered has C 4 rotational symmetry, that is, rotations through any multiple of π 2 radians maps the lattice onto itself. In [7], a hexagonal lattice with C 6 rotational symmetry is analysed, here, a rotation through an angle which is a multiple of π 3 maps the lattice onto itself. Although this lattice was formed of tessellating triangles, it is the rotational symmetry that gives the hexagonal lattice its name. In both cases, the method of multiple scales was applied, leading to an approximation for small amplitude breathers and their properties. Asymptotic estimates for breather energies were found, confirming the existence of minimum threshold energies obtained by Flach [14]. Numerical simulations showed that there was no restriction on the allowed direction of travel. This result contrasts with the behaviour of the mechanical lattice analysed by Marin et al [23], who find breathers which only travel along lattice directions.
This model we consider is simplified, in that only weak nonlinearities are considered, and includes no onsite potential. We investigate the behaviour of DBs on the two-dimensional honeycomb lattice shown in figure 1. The lattice possesses C 3 rotational symmetry, being made up of tessellating hexagons, in which rotation through any angle of a multiple of π 2 3 leaves the lattice invariant. Our aim is to investigate the combined leading-order effects of nonlinear nearest-neighbour interactions and the honeycomb geometry, by finding leadingorder asymptotic forms of DBs in this lattice. This complements previous studies of square and hexagonal lattices [6,7]. Numerical studies of Marin et al [23,24] required the use of an onsite potential as well as nonlinear nearest-neighbour interactions to general breathers in 2D lattices. One aim of the current work is to provide parameter regimes and initial conditions where breathers may exist in a system with only nonlinear nearest-neighbour interactions. We follow a similar analytic procedure to that of [6,7], using the method of multiple scales to obtain a system of equations from which we derive a nonlinear Schrödinger (NLS) equation that allows us to determine approximate small amplitude breathers. However, the analysis of the honeycomb lattice is significantly more complicated than the square or hexagonal cases due to the geometry of the latticeʼs interconnections which mean that there are two distinct types of node, which we call left-facing and right-facing. The analysis is similar to that of diatomic lattices, as it supports two types of mode which can be termed 'acoustic' and 'optical'.
In section 2 we derive the governing equations and the Hamiltonian structure behind them. We use the method of multiple-scales in section 3 to determine approximations to small amplitude breathers. Taking the amplitude, ε, as our small parameter, we form a power series expansion, equating terms at each order in ε and each harmonic of a fundamental linear mode. A dispersion relation is found in section 3.2, the plot of which shows some of the symmetry properties that we expect to find in the lattice. We find, in section 3.8, a reduction of the governing lattice equations to an NLS equation in two special cases. The first case, analysed in section 4, is where the interaction potential is symmetric, and along with the NLS equation, we find an ellipticity criterion for moving breathers; this means that only certain combinations of wavenumbers may produce a moving breather. The second special case, investigated in section 5, covers asymmetric potentials, but is restricted to stationary breathers. A relationship between the coefficient of the quadratic and cubic nonlinearities is derived in this case. We conclude our findings in section 6 with a summary of the results derived, and suggestions for further study.

Geometry of the lattice
We consider the nodes of the hexagonal honeycomb lattice as lying on a subset of a rectangular lattice. First, we introduce orthonormal basis vectors The position of the (m, n) node of the rectangular lattice is 3 so that the resulting hexagons are regular. In order to specify the honeycomb lattice, we retain only those nodes (m, n) for which + m n is an even integer and omit = + m p 6 1, = n odd and = + = m p n 6 4, even. In figure 1 the filled circles denote the nodes retained in the honeycomb lattice which satisfy these relations, and the open circles show all the remaining nodes in the underlying rectangular lattice.
To derive governing equations of this lattice we introduce vectors , as shown on the right hand side of figure 2, to describe the two configurations by which nodes are connected to nearest neighbours. The honeycomb lattice is composed of two distinct arrangements of connecting nearest neighbour nodes, shown in figure 3. We refer to arrangement 1 as left-facing nodes, since they are connected to a nearest neighbour horizontally to the left. Arrangement 2 will be referred to as right-facing nodes. When looking at figure 3 we note that each node is connected to three nodes of the opposite arrangement.

Derivation of the governing equations
In the application we consider here, at every node there is a nonlinear capacitor, and between adjacent nodes, a linear inductor. We denote the voltage across the capacitor (m, n) by V m n , and the total charge stored on this capacitor by Q m n , . Finally, the current in the direction of the vector e i , through the inductor immediately to the right of (m, n) is denoted by I m n , , and J m n , and K m n , represent currents in the direction of the vectors e j and e k , respectively. This configuration is illustrated in figure 2, which shows an enlarged view of the lattice with the relevant currents indicated.
We derive separate governing equations for the two arrangements of nodes, and only later aim to reconcile the two into a single description. Our aim is to find equations for the variable, Q m n , at each node of the lattice. To enable equations to be derived, we need to make the distinction between left-and right-facing nodes. We useQ m n , for left-facing nodes, that is, arrangement 1 in figure 3(a), andQ m n , for right-facing nodes, namely arrangement 2 in figure 3(b). We use Q m n , for a general node, in practice, this will be either one of the leftfacing (Q m n , ) or the right-facing (Q m n , ) nodes. A derivation from Kirchoffʼs laws has been given in [7]. Here we simply quote the Hamiltonian and note that P m n , and Q m n , are the conjugate momentum and displacement variables of the system. The charge-voltage relationship is given by where we assume the potential Υ Q ( ) 4 . Since our analysis is based on small amplitude nonlinear expansions, we assume that there is a Taylor series of Υ Q ( ); we do not consider potentials of the form Υ In section 4.3 we use the Hamiltonian (1) to find the energy of small amplitude breathers.
The lattice equations are obtained by eliminating P m n , from the equations

Asymptotic analysis
We now aim to find an approximate analytic solution to the equations (3) and (4) by applying the method of multiple scales. We first rescale the current variables m n , and t, introducing the new variables 2 with ε ≪ 1 being the amplitude of the breather, the variables X Y , will be treated as continuous real variables.
We require different ansatzes for the right-and left-facing nodes, therefore we analyse each type of node individually using (3) and (4). For right-facing nodes we seek solutions of the form ... c.c., where P Q R , , j j are all functions of τ X Y T ( , , , ). After substituting the ansatzes (6) and (7), into the relevant right-and left-facing lattice equations (3) and (4), we equate the coefficients of each harmonic frequency at each order of ε to find two sets of equations, which we analyse in order below. We use the slightly unusual notation  ε ψ ( e ) p q i to mean those terms of  ε ( ) p which have the coefficient ψ e q i , that is, we neglect those terms which have ψ e r i with ≠ r q. Since our main calculations are only going as far as  ε ( ) 3 and  ε ψ ( e ) 4 0i , the variables (5) are sufficient for our analysis. At higher orders of ε, we would have to include longer space scales, given by ε = X m 2 and ε = Y hn 2 , however, these make little difference to the shape of the breather, as shown in [30].
The first order we investigate is  ε ψ ( e ) i , whence we obtain and β* is its complex conjugate. We write β β = | | θ − e i , the magnitude being We are interested in solutions where ≠ 0 ( ) The dispersion relation describes the dependence of the temporal frequency of the wave on the wavenumbers (k, l). The negative square root in (11) leads to an 'acoustic' branch, or surface in ω k l ( , , ) space with lower frequencies, which we denote by ω β = − | | 3 ac ; and we have ω → 0 ac as → k l , 0. The surface corresponding to the positive root in (11), which clearly has larger values of ω, we denote by ω β = + | | 3 opt , and we describe this surface as the 'optical' branch. The acoustic branch accounts for frequencies in the range ω ⩽ ⩽ 0 3 , whilst the optical branch satisfies ω ⩽ ⩽ 3 6. The plot of ω against k and l along with the contour plot is shown in figure 4. We have the dispersion relation (11) for the two coupled systems (3) and (4). We consider k and l such that because ω is periodic in both k and l, with period π 2 in the k-direction and period π h 2 in l-direction. The locations in (k,l) space of the minima of ω ac and the maxima of ω opt coincide are seen in figure 4 as the circles in the centres of the hexagonal shapes in the contour plot. These points are at (0, 0), π (2 3, 0), π h (0, 2 ), π π h (2 3, 2 ) and π π h ( 3, ) etc. Points where the two surfaces meet are also evident in figure 4 as the centres of the triangles surrounding the hexagonal shapes, at these points where ω = 3 . The ω k l ( , ) dispersion surfaces have cusplike singularities at these points, which can be denoted by k 1 ,...,k 6 , where π π π π π π π π π π By comparing (11) with (10), we observe that these points occur where β = 0. In graphene, these wavevectors are known as Dirac points [26]. figure 4 also illustrates the hexagonal symmetry of the lattice. figure 5 shows the magnitude and argument of β as function of (k, l).
Note the presence of sizable plateaus where arg β π ≈ ± ( ) 0, 2 3. However, equation (8) remains unsolved. Since det M ( ) = 0, solutions can be written as P = CF, where, for ω ac and ω opt , we have ac i opt ac i respectively, the latter expression arising from β Γ = θ − e i . These expressions for C ac , C opt will be used in later calculations, where we find expressions for the functions G 2 , G 1 , G 0 , Q 2 , Q 1 and Q , 0 in terms of F.
: relationship between G 0 and Q 0 At  ε ψ ( e ) 2 0i , we obtain the same equation from both (3) and (4), which are the equations for Q m n , and for Q m n , . 0 Note that from the ansatz, Im G ( ) 0 is irrelevant, since only the combination + * G G 0 0 ever appears in our equations. Hence, we assume = Im G ( ) 0 0 , and only consider the real parts, that is, (14), and = G Q 0 0 , but this quantity is not yet determined.
: expressions for G 2 and Q 2 As previously mentioned, we aim to express all the variables G G G Q Q , , , , by substituting the ansatzes (6)- (7) into (3)- (4), we obtain where γ* is the complex conjugate of γ and Note that if we think of β and γ as being functions of (k, l), they are related by γ β = k l k l ( , ) (2 , 2 ), compare (9) and (17). Solving the linear system (15)-(16) for G 2 and Q 2 as functions of F, we find Whilst the bottom term in the vector is the complex conjugate of the top, we do not have = * Q G 2 2 since the the term CF 2 common to both is not necessarily real. We return to the expression (18) in section 5.
We now consider the governing equations at  ε ψ ( e ) 2 i , which can be written as where M is the matrix given in (8), and β k , β l are the partial derivatives of β with respect to k, l respectively, namely  (19), namely d, is in the range of the matrix M, which is given by ac i opt i Since normals to these directions are given by This equation implies that F (and hence P as well) is a travelling wave. We write the horizontal and vertical velocity components are found to be β β l l As expected from the standard theory of waves [32] these are simply the derivatives of the frequency with respect to the wavenumber, Since we have different expressions for ω ac and ω opt , equations (25)-(26) also generates different formulae for u ac and u opt (and v ac and v opt ). Figure 6 shows plots of the horizontal and vertical components of the velocity as functions of the wavenumbers (k, l). Note that at the Dirac points, where β = 0, the singularity is removable since the numerators in (25)- (26) are also zero.
The overall speed, c, is given by which is plotted in figure 7. Whilst the above calculations, (25)- (27), are for the acoustic mode, similar calculations for the optical mode produce similar results. All these plots show periodic behaviour, however, the hexagonal symmetry of the system only becomes clear in the total speed, the plots of the velocities u v , have a more complicated, although complimentary form. The velocities u v , both show sensitive dependence on wave vector (k, l). At the wavevectors k 1 ,..., k 2 , found in (12), both the components of velocity are zero.
The above calculation gives the condition on the rhs of (19) for solutions to exist. However, the quantities G 1 , Q 1 remain unknown. The solution of (19) is degenerate, and the one-parameter family of solutions may be written as  for arbitrary G 1 . The two equations forĜ 1 from (19) are then identical, and are solved by Whilst, this leaves G 1 undetermined, the quantityĜ 1 describes a small difference in the evolution of the left-and right-handed nodes of the honeycomb lattice.
3.6. O ε 3 e 0iψ À Á : corrections to the slow mode At  ε ψ ( e ) 3 0i , we obtain the equation from the substitution of (6)-(7) into both (3) and (4). We are only interested in determining the leading order terms εF, εP, which require the  ε ( ) 2 terms G Q , 0 0 , so we do not pursue the determination of H R , 0 0 , which are  ε ( ) 3 correction terms and provide only a small difference between the right-facing and left-facing nodes.
In general, we cannot solve (32) to find G 0 and Q 0 , but there are two special cases when we can do so. In the cases considered later, eitherˆ= H 0 or it is not relevant to our calculations. We also choose = −Ĝ G 1 2 with (29), so that (32) can be simplified to 2 , which is similar to the previously derived results for the square and hexagonal lattices, see equations (2.23) of [6] and (2.23) of [7].
The first special case, analysed in section 4, is when the interaction potential is symmetric, that is, . Under this assumption, the quadratic coefficient of the force, a, is zero, leading to as the solution of (32). In this case, we also have . This means that the system is governed by equations (33)-(35) given below, which can be reduced to a single NLS equation.
In section 5, we consider the second case, ω ω = max where ω max represents the maxima of ω. In this case the breather is stationary, since (25)-(26) yield u = 0 and v = 0, then the system as a whole has no τ-dependence, that is, = = 2 . This solution can be substituted into equations (33)-(35) given below which again can be reduced to a single NLS equation.

NLS equation
The final equation we need to investigate comes from terms of  ε ψ ( e ) 1 1 where the matrix M is identical to that in (8), and the rhs components are given by 4e e e e e e e 3 9  , that is, + = CA B 0, (ii) substitute in expressions for G 2 , Q 2 , G 0 , Q 0 , (iii) making the substitution P = CF, (iv) transforming to travelling wave coordinates by (24).
However, in general, this equation will still be coupled to G 0 , through (32); and carrying out this procedure in general leads to extremely lengthy expressions.

Summary
We have derived a multiple scales asymptotic expansion for envelope solutions of the scalar two-dimensional honeycomb lattice. After finding the usual expressions for the frequency, and group velocity, we have found a coupled system of PDEs for the shape of the envelope given by (32) and (33)-(35). Whilst we cannot, in general, solve this resulting system of equations, there are two special cases in which the system reduces to a single NLS equation. The general case shares some similarities with the Davey-Stewartson system of equations [10] obtained in fluid mechanics. The remainder of this paper is not directed to a general analysis of equations (32) and (33)-(35), rather we consider two special cases in more detail.
These two cases are considered in more detail in sections 4 and 5 respectively, and in both cases the procedure (i)-(iv) leads to substantially simpler expressions than the general case. In both special cases we find additional criteria which, if not satisfied, mean the the lattice cannot support breather solutions.

The symmetric potential (a = 0) and moving breathers
In this section we consider the simplified case where a = 0 in (3) and (4)  , that is, there is no generation of second harmonics. Furthermore, from section 3.7 we gain = = G Q 0 0 0 , which satisfies the relationship = G Q 0 0 from section 3.3. This means that there is no 'slow' mode, which is independent of t (corresponding to ω = 0), and the localized mode evolves only on the slower τ T , timescales.

The elliptic NLS equation
To make further progress towards understanding the form of possible solutions of the system (37)-(39), we make the substitution The NLS equation in 2+1 dimensions has two forms depending on whether the seconddifferential operator part of the equation is elliptic or hyperbolic. We are only interested in elliptic systems (where the coefficients of ξξ F and ζζ F have have the same sign) as our aim is to find solutions which are localized in both spatial dimensions. We therefore define the ellipticity as where expressions for D W , D Z and D M are given in (38)-(39). Since we have two expressions for C, one for the acoustic mode and the other for the optical mode, as given in (13), we have different expressions for D Z , D W , D M in the two cases, and two expressions for the ellipticity,  ac and  opt .
In figure 8 we plot the sign of the ellipticity functions  k l ( , ) from (42) for the acoustic and optical cases. In the acoustic case,  ⩽ 0 ac for almost all (k, l), there being small trefoilshaped areas of positive ellipticity near the Dirac points. However, for the optical case, there is a wide range of wavenumbers where  > 0 opt , as shown by the white hexagonal areas in the right panel of figure 8. Note that the optical case also shows small trefoil-shaped areas of positive ellipticity near the Dirac points. Breathers corresponding to these wavenumbers are expected to be unstable as their frequencies will coincide with those of linear waves. Whilst the dispersion relation in this diatomic system does not have a gap-the frequency spectrum ω k l ( , )) includes all values from zero to ω = 3 max , the form of the breathers near the Dirac points are expected to be similar to those of gap solitons in other diatomic systems, where the dispersion relation has gaps. One-dimensional FPU problems have been studied by Livi et al [20] and James and Noble [16]. As in one-dimensional diatomic systems, the breathers corresponding to the optical domain including = k l ( , ) (0, 0) have frequencies which lie above the optical band, and so are expected to be long-lived.
We now focus the optical case, where = = − θ C C e opt i , and (41) simplifies to In order for bright breathers to exist, there is a second criterion to be satisfied, namely that the coefficients of the nonlinear term and the spatial derivative must have the same sign. Since the nonlinearity is negative, we require For the optical mode, this condition is satisfied for all (k,l).

Asymptotic estimates for breather energy
The total electrical energy in the honeycomb lattice is conserved. This quantity is related to the Hamiltonian (1) by =Ẽ H C 0 . Thus, upto quadratic order, Now our aim is to work out an expression for energy at leading order in ε, given our solution forQ, Q in terms of F. Since we are only interested in leading order approximation to the energy, the dependence of the solution for F given by (52) on T can be ignored, as the dependence on ω dominates. However, in passing we note that from (52) that the combined frequency of the breather mode is given by 2 2 and so, in the optical case, the frequency lies above the frequency of linear waves. From The equation for the energy is given by (44), here we show the calculation of the onsite (Q m n , , Q m n , ) part of this, the calculation of the interaction energy (due to P m n , ,P m n , ) can be found in a similar way and gives an identical final expression.
We replace the double sum over (m, n) in (44) by an integral over (X, Y)-space using (5), and transform into an integral over (Z,W)-space using (24). The Jacobian required to then make the transformation from (Z, W) to ξ η ( , ) coordinates using (40) is . The final stages use since r is slowly varying in m n , . From equation (48), we note that to leading order, the energy of the breather is independent of the breather amplitude, A. Thus, no matter how small the breather amplitude, there is a minimum energy required to create it. The reason for this threshhold energy is that as the amplitude reduces, the width of the breather increases, in such a way that the total energy remains constant. This property confirms the observations of Flach et al in [14].
However, the threshold energy is dependent on the wavenumbers k and l, therefore moving breathers have different threshold energies. Figure 9 shows that the energy threshold is locally maximized at = = k l 0, that is, for static breathers. Moving breathers require less energy to form. An alternative viewpoint is that as breathers lose energy, they start moving, and accelerate, to the maximum speed, where the ellipticity constraint is only just satisfied. It is also clear from (48) that the energy is closely related to the ellipticity constraint. Finally, we note that the breathers predicted near the Dirac points, which have frequencies lying in the linear spectrum, have much higher energies than the out-of phase optical breathers whose frequencies lie above the top of the linear spectrum.

Static breathers in an asymmetric potential
In this section we examine the more general case for which ≠ a 0 in (3) and (4). Recall that at the end of section 3 we obtained a system of two coupled equations for G 0 and F, namely (32) and the equation that can be derived by following the procedure (i)-(iv) in section 3.8. It is only possible to reduce this system to a single solvable equation when a = 0 (as analysed above in section 4) and when G 0 is independent of τ, which we discuss in this section. The only example where the system becomes independent of τ is the case = = k l 0, on the optical branch. Under these conditions, we have, from (9), (11), (13), (20), (25), (26), (24), (30) For bright breather solutions to exist, we require the coefficients of the nonlinearity and the spatial diffusion terms to have same signs, that is, > b a 4 3 2 . In place of (45), the breatherʼs frequency is now given by Ω ε , which still lies above the top of the phonon band. As noted above, the second harmonic terms, G 2 and Q 2 are both zero for the case of stationary breathers, that is the optical mode with = = k l 0. Hence we extend the expansion of section 3 to consider the terms at  ε ψ ( e ) 3 3i to see if third harmonic terms are generated. At  ε ψ ( e ) 3 3i , we obtain similar equations to those of section 3.4, more specifically we obtain 3 . Thus we find that the honeycomb lattice generates third harmonics but not second harmonics in the stationary breather.

Comparison with other lattice geometries
In earlier papers [6,7] we have carried out similar calculations on the square and hexagonal lattices, where the derivations are considerably simpler. In all cases we have = − | | G a F 0 2 ; however, other properties of the lattices differ, depending on the geometries concerned. In table 1 we compare the results of the honeycomb lattice analysed here with corresponding results for the square and hexagonal lattices analysed earlier.
The absence of any second harmonic is a property shared with the square lattice. Whilst the hexagonal lattice generates no third harmonic, it does generate a second harmonic. Furthermore, the inequality relating the coefficients of nonlinear terms is identical for the honeycomb lattice and the square lattice, whilst different for the hexagonal. The possibly surprising result from this table is that, at least as far as stationary breathers are concerned, the honeycomb lattice has more in common with the square lattice than the hexagonal lattice.

Stability of the breather
The solution for F is a one-parameter family, which we parametrize by the amplitude, A, as where ϕ r ( ) is the function which solves the elliptic problem  ϕ ϕ ϕ = − 2 3 in two dimensions. This elliptic problem is known to have solutions, and the cylindrically-symmetric solution we write as ϕ r ( ). Solutions such as (52) are known as Townes soliton solutions [9] of the 2D NLS. These solutions are known to be unstable in the 2D NLS, with subcritical initial conditions suffering from dispersion, leading to the amplitude converging to zero everywhere through the initial data spreading out; whilst supercritical initial conditions blow up, with the energy being focused to a single point, where the amplitude diverges. However, arbitrarily small structural perturbations to (43) can stabilize the Townes soliton. For example, results proven by Davydova et al [11] for the equation demonstrate the stability of a localized breather mode if > BD 0 and > PK 0. Clearly the presence of a higher derivative terms can mollify the blow-up singularity, whilst higher order nonlinearities can also help stabilize the soliton, as discussed by Kuznetsov [18].
If we pursue higher order correction terms, for example, from  ε ψ ( e ) 5 i terms, then terms such as  F 4 and | | F F 4 occur, which may stabilize the Townes soliton provided their coefficients have the correct combinations of signs. However, such an expansion also yields terms of the form  | | F F ( ) 2 2 ,  | | F F 2 2 and  * F F 2 2 , and the effect of such structural perturbations of (43) has, to our knowledge, not yet been determined. Fibich and Papanicolaou [12,13] have also addressed this problem, though their results do not yet extend to these nonlinear derivative terms.
To illustrate this, let us consider the case of stationary breathers on a symmetric lattice, that is, we take the nearest-neighbour restoring force to be ϕ ϕ 5 (that is, a = 0 and no quartic nonlinearity). We analyse the special case given by = = k l 0, so that = = u v 0, . Note that we also haveˆ=ˆ= u v 0 so thatˆ= G 0,   Here, in addition to the long scales defined in (5), we have introduced even longer time and length scales given by , and  is the corresponding vector derivative with respect toX andỸ .
Since ω = 6 , the right-hand-sides of (55)-(56) must be equal. Combining this with the relation = − P F leads to since the third-derivative terms cancel. Whilst these terms generate nonzero solutions for J U , 1 1 , such contributions do not concern us here, where our aim is to determine the properties of F. The effect of theF T term is to change the timescale slightly, and the   F . term rescales the space scale X, hence we will neglect these terms.
Applying ω∂ 2i T to the leading order form of (57), which is (50) in the case a = 0, yields which we use to eliminate F TT from (57), to find the final governing equation second derivatives of cubic nonlinearities can be bound by some combination of fifth order nonlinearities and fourth order derivatives of F, it is reasonable to assume that for sufficiently large g, the breather will be stable.

Conclusions
We have investigated the properties of DBs on a 2D honeycomb lattice. After applying Kirchoffʼs laws to the electrical lattice, we derived a governing set of equations for the case of nonlinear capacitors at nodes, and nodes being connected by linear inductors. Using multiple scales asymptotic methods, we reduced the governing equations to a single NLS equation from which we can determine the properties and conditions under which small-amplitude breathers may exist. There are two cases in which an NLS equation can be obtained. We analysed each case in more detail in sections 4 and 5.
The analysis of the honeycomb is more complicated than either the square or the triangular lattices, due to the necessity of treating the two types of node, which means that a diatomic analysis must be carried out. This leads to extra complications at the level of determining the evolution of the 'slow mode' at  ε ( e ) 4 0 . Part of the extra complexity is that in order to derive the leading order F and P terms in (6) and (7), it is necessary to simultaneously find the first correction terms G 1 and Q 1 .
The first special case we considered ( §4) was that of a symmetric potential in which the terms G 0 , Q 0 , Q 2 and G 2 are all zero. From this we were able to obtain an ellipticity condition, for the wavenumbers (k, l), to ensure we obtained solutions which were localized in both spatial directions. A minimum threshold energy to create breathers was also found. This confirmed the observations of Flach et al [14]. The ellipticity condition, breather energy and dispersion relation, obtained in sections 3.2 and 4, were plotted. The breather energy is maximized for these stationary breathers. For other wave vectors, moving breathers are created, with lower energies.
The second case we analysed was the case of asymmetric potentials. Here we only considered the specific wavenumber = = k l 0, and the optical branch which guarantees stationary breathers. However, this enabled us to describe the behaviour of a range of nonlinearity parameters (a b , ) for which stationary breathers may exist. We find no second harmonic term in the expansion in this case, but there is a third harmonic. These properties show a close similarity between the square lattice and the honeycomb, quite distinct from the hexagonal lattice.
It is natural to consider the relationship between the honeycomb system studied here and a one-dimensional systems. We note that the two-dimensional systems studied previously [6,7] both had a dispersion relation with a single branch that described modes with optical and acoustic characters, as in one-dimensional (monatomic) systems. However, in onedimensional diatomic systems, the dispersion relation has two branches, one optical and one acoustic. Such systems have been studied by Livi et al [20] and James and Noble [16], amongst others. In the latter paper, the authors derive the dispersion relation, showing it to have two branches: an optical and an acoustic form, as in the honeycomb lattice. However, in the one-dimensional diatomic lattice, the two branches are distinct, and do not meet at Dirac points; rather, there is a gap between the two branches, in which breathers may exist with frequencies which are not coincident with any linear wave. In the honeycomb lattice, the two branches meet at the Dirac points, and so there is no gap. However, from figure 8 we see that the honeycomb lattice still supports breather solutions near the Dirac points. Whilst it would be interesting to investigate these solutions further, we expect them to be unstable, since their frequencies coincide with those of linear waves, allowing energy interchange with phonons which could lead to the breatherʼs decay. In contrast, we now consider the larger white regions in figure 8 (right), corresponding to wavenumbers for which optical breathers exist. We expect these breathers to have frequencies above the top of the optical band, and so will have no linear wave with the same frequency. These waves, however, may still be unstable, due to other effects, such as the breatherʼs motion over lattice sites being resonant with linear modes. Such losses may still allow the breather to travel long distances before decaying, and so be relevant in applications such as the explanation of tracks in mica via quodons as suggested by Russell and Eilbeck [27,28].
In this paper we have only looked at a scalar-valued quantities at each node that is, only one degree of freedom. In future works we aim to analyse the stability of these breather solutions, and find approximate solutions to the vector-valued honeycomb lattice similar to the lattices Marin et al studied in [23].