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 Schrodinger (NLS) 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 maximised 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.


Introduction
Discrete Breathers (DBs) are time-periodic and spatially-localised exact solutions which describe the motion of a nonlinear lattice, that is, a repeated arrangement of atoms. In this paper, we investigate the properties of discrete breathers on a two-dimensional honeycomb lattice, seeking conditions under which the lattice may support breather solutions.
The combination of nonlinear interactions and discreteness gives rise to breather modes. The discreteness causes gaps and cutoffs in the phonon spectrum, whilst nonlinearity allows larger-amplitude waves to have frequencies outside the phonon band.
MacKay & Aubry's work [22] established the existence of breathers in one-and higherdimensional lattice systems. Flach et al. [15] have shown that properties of breathers in the more familiar one-dimensional systems apply also to lattices in higher dimensions. As well as this analytical work, numerical methods have also been applied to DB's in higher dimensional systems. For example, Takeno [29] used lattice Green functions to determine approximations to breather solutions in one-, two-and three-dimensional lattices. Burlakov et al. [5] found breather solutions numerically on a two dimensional square lattice and in [4], Bonart et al. simulated numerically localised excitations on one-, two-and three-dimensional scalar lattices.
Although the existence of breathers does not depend on the lattice dimension, some of their properties do. Flach et al. [14] found that minimum energy threshold in order to create breathers if the lattice dimensions is equal to or greater than a critical value. This threshold energy is the positive lower energy bound attained by the breather. Strictly, only stationary breathers are necessarily time-periodic; however, MacKay and Sepulchre [21] have formulated a more precise definition of travelling breathers. Moving breathers in two-dimensional lattices were investigated in a collection of papers by Marin, Eilbeck and Russell who were motivated by the observation of dark lines formed along crystal directions in white mica [27]. In mica, potassium atoms lie in planes in which they occupy a hexagonal pattern. Numerical simulation of Marin et al. [23] exhibited moving breathers which only travelled along lattice directions. Similar results were observed in a further study of Marin et al. [24], where two-and three-dimensional lattices of various geometries were investigated. The mechanical lattice, in which each node can move horizontally and vertically is highly complex, and although attempts at a full asymptotic analysis have been made (for example, [31]), the detailed understanding of dynamics of such a system is not yet available.
Currently, there is great interest in the behaviour of honeycomb lattices, due to the development of potential applications of graphene. For example, Molina and Kivshar [25] studied the localisation and propagation of light along ribbons which have a honeycomb structure analogous to graphene. Bahat-Treidel et al. [2] have studied the propagation of a field in a photonic lattice with Kerr nonlinearity. They show that, in the honeycomb lattice, the Kerr nonlinearity produces waves with triangular symmetry. Chetverikov et al. [8] consider a system with Lennard-Jones-like interaction potentials. Using a variety of initial conditions, they use numerical simulations, to find outputs which bear strong visual similarity with results of bubble chamber experiments. A system of spherical particles in a hexagonal structure interacting with nearest neighbours via Hertzian contacts is considered by Leonard et al. [19]. They analyse the waves that spread through the system following a localised impulse. Kevrekedis et al. [17] consider interactions which include longer-range as well as nearest neighbours in a DNLS model. They find that these can stabilise and destabilise solitons. Ablowitz and Zhu [1] use perturbation theory to analyse the linear spectrum of a hexagonal lattice near its Dirac point, as well as the associated Bloch modes and envelope solutions.
Herein, we consider an electrical transmission lattice, in which a scalar quantity, for example the charge stored on a nonlinear capacitor is defined at each node, with nodes being coupled by linear inductors. This paper follows on from from previous work of Butt and Wattis [6,7] who studied discrete breathers in two-dimensional square and hexagonal electrical lattices. In such lattices, the scalar valued functions at each node can be thought of as charge, thus there is only one degree of freedom at each node. This contrasts with the models simulated by Marin et al. [23,24] where there is a vectorvalued 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 discrete breathers 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 leading-order asymptotic forms of discrete breathers 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 two-dimensional 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 nearestneighbour 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 B = {i, j}, where i = [1, 0] T and j = [0, 1] T . The position of the (m, n) node of the rectangular lattice is mi + hnj, with h = √ 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 = 6p + 1, n = odd and m = 6p + 4, n = 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.  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 use Q 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 left-facing ( 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 H = m,n;s.t. Pm,nexists 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 V (Q m,n ) = Υ ′ (Q m,n ) where we assume the potential Υ(Q m,n ) has the form Υ(Q) = 1 2 Q 2 + 1 3 aQ 3 + 1 4 bQ 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 Υ(Q) ∼ Q ν with ν < 2. 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 Thus, for left-facing nodes we have where m, n ∈ Z Z, Q m,n represents the charge at left-facing nodes andQ m,n represents the charge at right-facing nodes. The right-facing nodes in arrangement 2 are governed by

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 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 where the phase of the carrier wave ψ is given by ψ = km + lhn + ωt, where k = [k, l] T is the wavevector and ω(k) is its temporal frequency. Similarly, for left-facing nodes we seek solutions of the form where P, Q j , R 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 O(ε p e iqψ ) to mean those terms of O(ε p ) which have the coefficient e iqψ , that is, we neglect those terms which have e irψ with r = q.
Since our main calculations are only going as far as O(ε 3 ) and O(ε 4 e 0iψ ), the variables (5) are sufficient for our analysis. At higher orders of ε, we would have to include longer space scales, given by X = ε 2 m and Y = ε 2 hn, however, these make little difference to the shape of the breather, as shown in [30].
We are interested in solutions where ( F P ) = 0, equation (8) is thus an eigenvalue problem. We require the determinant of the matrix to be zero, which gives the dispersion relation 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 ω ac = 3 − |β|; and we have ω ac → 0 as k, l → 0. The surface corresponding to the positive root in (11), which clearly has larger values of ω, we denote by ω opt = 3 + |β|, 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 (k, l) ∈ T 2 = [0, 2π] × [0, 2π/h] because ω is periodic in both k and l, with period 2π in the k-direction and period 2π/h 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), (0, 2π/h), (2π/3, 2π/h) and (π/3, π/h) 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 cusp-like 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 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 .   At O(ε 2 e 0iψ ), we obtain the same equation from both (3) and (4), which are the equations for Q m,n and for Q m,n .
Note that from the ansatz, Im(G 0 ) is irrelevant, since only the combination G 0 + G * 0 ever appears in our equations. Hence, we assume Im(G 0 ) = 0, and only consider the real parts, that is, (14), and G 0 = Q 0 , but this quantity is not yet determined.

3.4.
O(ε 2 e 2iψ ): expressions for G 2 and Q 2 As previously mentioned, we aim to express all the variables G 0 , G 1 , G 2 , Q 0 , Q 1 and Q 2 in terms of F . At O (ε 2 e 2iψ ) 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) = β(2k, 2l), 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 2 = G * 2 since the the term CF 2 common to both is not necessarily real. We return to the expressions (18) in Section 5.

O(ε 2 e iψ ): the velocity profile
We now consider the governing equations at O(ε 2 e 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 Since det(M) = 0, an equation such as (19), which we write as M( G 1 Q 1 ) = d, either has no solutions, or a whole family of solutions for (G 1 , Q 1 ) T . According to the Fredholm alternative, the existence of solutions depends on d. Solutions exist only if the rhs of (19), namely d, is in the range of the matrix M, which is given by Since normals to these directions are given by the condition that d ∈Range implies n.d = 0. Note that in both the optical and the acoustic cases, (22) implies n = ( C 1 ). We also recall that the leading order quantities, P and F are related by P = CF , where both C and n = ( C 1 ) have different expressions for the acoustic and optical cases, given by (13). Using P = CF and n.d = 0 we obtain the equation 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 As expected from the standard theory of waves [32] these are simply the derivatives of the frequency with respect to the wavenumber, u = ∂ω/∂k, v = ∂ω/∂l. 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 sin 2 (lh)[cos(3k) + 2 cos(lh)] 2 + 3 sin 2 (3k) cos 2 (lh), 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 solutions of (19) are degenerate, and the one-parameter family of solutions may be written as for arbitrary G 1 . The two equations for G 1 from (19) are then identical, and are solved by Writing this as G 1 = uF Z + vF W , we have Whilst, this leaves G 1 undetermined, the quantity G 1 describes a small difference in the evolution of the left-and right-handed nodes of the honeycomb lattice.

O(ε 3 e 0iψ ): corrections to the slow mode
At O(ε 3 e 0iψ ), we obtain the equation from the substitution of (6)-(7) into both (3) and (3). We are only interested in determining the leading order terms εF , εP , which require the O(ε 2 ) terms G 0 , Q 0 , 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 G with (29), so that (32) can be simplified to G 0,τ τ = 3∇ 2 (G 0 + a|F | 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, V (Q m,n ) = V (−Q m,n ). Under this assumption, the quadratic coefficient of the force, a, is zero, leading to G 0 = Q 0 = 0 as the solution of (32). In this case, we also have G 2 = Q 2 = 0. 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, P τ τ = F τ τ = 0. In this case, G 0 , Q 0 also have no τ -dependence, and so equation (32) reduces to G 0 = Q 0 = −a|F | 2 . This solution can be substituted into equations (33)-(35) given below which again can be reduced to a single NLS equation.

Nonlinear Schrödinger equation
The final equation we need to investigate comes from terms of O(ε 3 e iψ ) which yield where the matrix M is identical to that in (8), and the rhs components are given by As in the case of the equations at O(ε 2 e iψ ), in order for this system of equations to have solutions, there is the consistency condition n.( A B ) = 0 that must be satisfied. An equation for F (Z, W, T ) can be obtained from (33)-(35) by the following procedure: (i) calculating the consistency condition n.( A B ) = 0, 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) so that V (−φ) = V (φ) and V ′ (−φ) = −V ′ (φ). In section 3.4, whilst equation (18) remains valid, we recover G 2 = Q 2 = 0, that is, there is no generation of second harmonics. Furthermore, from Section 3.7 we gain G 0 = 0 = Q 0 , which satisfies the relationship G 0 = Q 0 from Section 3.3. This means that there is no 'slow' mode, which is independent of t (corresponding to ω = 0), and the localised mode evolves only on the slower τ, T timescales.

O(ε 3 e iψ ) -derivation of NLS
We now apply the procedure (i)-(iv) from Section 3.8 and so simplify the NLS-like system (19)-(35). Taking G 1 , Q 1 as given by (28) with G = − 1 2 G and using P = CF to evaluate n.( A B ) = 0, we obtain a single equation for F , namely Thus we have completed stages (i)-(iii) of the procedure from Section 3.8. In stage (iv) we eliminate the τ -derivative terms using the travelling wave substitution F (X, Y, τ, T ) = F (Z, W, T ) with u and v representing the horizontal and vertical components of velocity found in (25)- (26). Using F τ τ = u 2 F ZZ + 2uvF ZW + v 2 F W W , we rewrite (36) as Hence, from the governing equations (3)-(4), we have found (37), which is an NLS equation in 2+1 dimensions.

The elliptic NLS equation
To make further progress towards understanding the form of possible solutions of the system (37)-(39), we make the substitution to remove the mixed derivative term. This transformation yields The NLS equation in 2+1 dimensions has two forms depending on whether the second-differential 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 localised 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, E ac and E opt . In Figure 8 we plot the sign of the ellipticity functions E(k, l) from (42) for the acoustic and optical cases. In the acoustic case, E ac ≤ 0 for almost all (k, l), there being small trefoil-shaped areas of positive ellipticity near the Dirac points. However, for the optical case, there is a wide range of wavenumbers where E opt > 0, 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 ω max = √ 3, 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 & 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 opt = −e 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 D Z < 0. 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 E = H/C 0 . Thus, upto quadratic order, m,n + (P m,n − P m+2,n ) 2 + (P m,n − P m−1,n−1 ) 2 Now our aim is to work out an expression for energy at leading order in ε, given our solution for Q, 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 and so, in the optical case, the frequency lies above the frequency of linear waves. From (6)-(7) we find Q m,n = 2εAφ(r) cos(km + lhn + ωt), Q m,n = − 2εAφ(r) cos(km + lhn + ωt + θ), where r is given by the argument of φ in (52), and using P = CF = −e iθ F since we are considering optical modes. For both left-and right-facing nodes, Q m,n and Q m,n , we have d 2 Q/dt 2 = −ω 2 Q, and dP/dt = −Q; hence, at leading order, dQ/dt = ω 2 P and P m,n = − 2εA ω φ(r) sin(km + lhn + ωt), P m,n = 2εA ω φ(r) sin(km + lhn + ωt + θ).
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.
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 maximised 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) Since u = v = 0, in this case, the breather is stationary; in addition, from (18), we find (γ = 3 and) G 2 = Q 2 = 0. Assuming G 0 is independent of τ , equation (32) can be solved by G 0 = Q 0 = −a|F | 2 , enabling us to perform stage (ii) of the process outlined in Section 3.8. Now we turn to deriving the NLS equation. From stage (i) in Section 3.8, and since C = −1, we form the equation A = B from (34)-(35). Since P = −F , stage (iii) of the calculation leads to 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 > 4 3 a 2 . In place of (45), the breather's frequency is now given by Ω = √ 6 + 3ε 2 A 2 (3b − 4a 2 )/2 √ 6, which still lies above the top of the phonon band.

O(ε 3 e 3iψ ): expression for the third harmonic
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 O(ε 3 e 3iψ ) to see if third harmonic terms are generated. At O(ε 3 e 3iψ ), we obtain similar equations to those of Section 3.4, more specifically we obtain These are the equations for general k, l; however, for stationary breathers we are only concerned with k = l = 0, in which case Q 2 = G 2 = 0, ω = √ 6 and P = −F , hence we obtain the solution H 3 = 1 8 bF 3 , R 3 = − 1 8 bF 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 0 = −a|F | 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.
Property \ Geometry Square [6] Hexagonal [7] Honeycomb Second harmonic 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 parametrise 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 cylindricallysymmetric 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 twodimensional 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 localised breather mode if BD > 0 and P K > 0. Clearly the presence of a higher derivative terms can mollify the blow-up singularity, whilst higher order nonlinearities can also help stabilise the soliton, as discussed by Kuznetsov [18]. If we pursue higher order correction terms, for example, from O(ε 5 e iψ ) terms, then terms such as ∇ 4 F and |F | 4 F occur, which may stabilise the Townes soliton provided their coefficients have the correct combinations of signs. However, such an expansion also yields terms of the form ∇ 2 (|F | 2 F ), |F | 2 ∇ 2 F and F 2 ∇ 2 F * , 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 V ′ (φ) = φ + bφ 3 + gφ 5 (that is, a = 0 and no quartic nonlinearity). We analyse the special case given by k = l = 0, so that u = v = 0, G 0 = Q 0 = Q 2 = G 2 = 0. Note that we also have u = v = 0 so that G = 0, H 0 = R 0 and we can take G = 0 so that G 1 = Q 1 = 0 = H 0 = R 0 = 0. Hence, in place of the ansatzes (6)- (7), used earlier, we use the simplified forms with H 3 = 1 8 bF 3 , and R 3 = − 1 8 bF 3 as derived in Section 5.1. Our aim is to calculate the form of the higher-order terms, namely those at O(ε 4 e iψ ) and O(ε 4 e iψ ).
Combining the results from O(εe iψ ), O(ε 3 e iψ ), O(ε 5 e iψ ), we obtain the governing equations Here, in addition to the long scales defined in (5), we have introduced even longer time and length scales given by T = ε 4 t, X = ε 3 m, and Y = ε 3 hm, and ∇ is the corresponding vector derivative with respect to X and Y .
Since ω = √ 6, the right-hand-sides of (55)-(56) must be equal. Combining this with the relation P = −F leads to |F | 4 F , since the third-derivative terms cancel. Whilst these terms generate nonzero solutions for J 1 , U 1 , such contributions do not concern us here, where our aim is to determine the properties of F . The effect of the F 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 T T from (57), to find the final governing equation As the last three terms do not appear in (53), we cannot formally determine the stability properties of the system. However, if we were to simply ignore the last three terms, (53) suggests that if g > 3 16 b 2 (and b > 0) then the combined influence of the fifth-order nonlinearities and fourth order derivatives stabilise the breather. Since we expect that the 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 discrete breathers on a two-dimensional 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 O(ε 4 e 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 localised 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 maximised 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 one-dimensional 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 & 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 [28,27].
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].