{'fraction_non_alphanumeric': 0.03807910839160839, 'fraction_numerical': 0.016144012237762236, 'mean_word_length': 4.735390882030393, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dirac plasmons in bipartite lattices of metallic nanoparticles

We study theoretically"graphene-like"plasmonic metamaterials constituted by two-dimensional arrays of metallic nanoparticles, including perfect honeycomb structures with and without inversion symmetry, as well as generic bipartite lattices. The dipolar interactions between localised surface plasmons in different nanoparticles gives rise to collective plasmons that extend over the whole lattice. We study the band structure of collective plasmons and unveil its tunability with the orientation of the dipole moments associated with the localised surface plasmons. Depending on the dipole orientation, we identify a phase diagram of gapless or gapped phases in the collective plasmon dispersion. We show that the gapless phases in the phase diagram are characterised by collective plasmons behaving as massless chiral Dirac particles, in analogy with electrons in graphene. When the inversion symmetry of the honeycomb structure is broken, collective plasmons are described as gapped chiral Dirac modes with an energy-dependent Berry phase. We further relax the geometric symmetry of the honeycomb structure by analysing generic bipartite hexagonal lattices. In this case we study the evolution of the phase diagram and unveil the emergence of a sequence of topological phase transitions when one hexagonal sublattice is progressively shifted with respect to the other.


Introduction
The way materials interact with light has been of interest, both scientifically and artistically, for millennia. The first lenses and mirrors were created to bend the trajectory of light, leading to a huge variety of practical applications. Among them, focusing light to small regions of space has been extensively used in classical optics to image tiny objects. However, the diffraction limit has seriously hindered our ability to observe microscopic structures with dimensions less than the wavelength of the detecting light [1]. This limitation has been overcome with the use of plasmonic nanostructures [2,3], such as isolated metallic nanoparticles [4]. When illuminated by an external radiation, the electrons in the nanoparticle oscillate collectively, forming a localised surface plasmon (LSP) resonance [5,6]. The evanescent field at the surface of the nanoparticle associated with the LSP [1,7] produces strong optical field enhancements in the subwavelength region. This phenomenon overcomes the diffraction limit and allows for resolution at the molecular level [8].
Within the field of plasmonics, single or few nanostructures have been the cynosure so far. However, the attention is rapidly shifting to metamaterials constituted by ordered arrays of nanostructures. The periodic arrangement of the nanostructures in the array, as well as the interactions between the LSPs on each nanostructure can be exploited in a variety of novel artificial materials exhibiting properties far beyond those seen in nature, leading to fascinating applications such as electromagnetic invisibility cloaking [9][10][11], perfect lensing [12,13] and slow light [14]. The interaction between LSPs on individual nanostructures generates collective plasmonic modes that extend over the whole array [15][16][17][18][19][20][21][22][23][24][25][26][27]. As a consequence the collective plasmons can exhibit a variety of properties that crucially depend on the lattice structure of the metamaterial, and on the microscopic interactions between LSP resonances.
Among the infinite possible nanostructure arrangements, one structure that has gained considerable attention in the condensed matter community within the last decade is the honeycomb structure of graphene (see figure 1), a two-dimensional (2D) monolayer of carbon atoms [28]. The honeycomb structure can be regarded as a special bipartite lattice formed by two inequivalent hexagonal sublattices, denoted by A and B, with equal nearest neighbour separation for all carbon-carbon bonds. The hopping of electrons between neighbouring atoms in graphene produces a spectrum characterised by the presence of massless fermionic Dirac quasiparticles [29][30][31][32]. These pseudo-relativistic Dirac fermions have an associated chirality which leads to several of the remarkable properties of graphene, such as a nontrivial Berry phase of π responsible for the anomalous quantum Hall effect [30,31], as well as the suppression of electronic back-scattering from smooth scatterers [33]. The latter property is ultimately responsible for the very high electron mobility of graphene samples.
The majority of the unique electronic properties of graphene stems from the honeycomb structure alone. This feature has been exploited in a variety of "artificial graphene" systems [34][35][36][37][38][39][40][41][42][43][44][45][46] that share the honeycomb lattice structure, albeit in different physical contexts. Among them, in a recent publication [44], it has been theoretically demonstrated that the properties exhibited by electrons in graphene can be replicated within the context of plasmonic metamaterials. In fact, a 2D honeycomb array of metallic nanoparticles has been shown to support collective plasmons (CPs) that behave as unprecedented chiral massless Dirac bosons. Moreover, the band structure of these collective modes was found to be fully tunable with the polarisation of the incident light that triggers the CPs, a feature which is not available Figure 1. (a) Honeycomb lattice with lattice constant a, consisting of two sublattices, labelled A and B, connected by the nearest neighbour vectors in real graphene devices. Honeycomb plasmonic arrays of metallic nanoparticles thus open interesting perspectives for the realisation of ultrathin tunable metamaterials where the electromagnetic radiation can be transported effectively by chiral pseudorelativistic Dirac modes.
In this paper we further explore the optical response of bipartite arrays of metallic nanoparticles beyond the artificial graphene plasmonic lattice. Firstly, in section 2, we exploit the tunability of the CP band structure and analyse in detail the phase diagram of gapless and gapped phases that emerge while tilting the orientation of the dipole moments associated with LSPs. We show that there are several topologically disconnected regions of the phase diagram which support CPs behaving as massless Dirac quasiparticles. In section 3 we then extend the analysis of the plasmonic array by considering a honeycomb structure with broken inversion symmetry in which the two inequivalent sublattices are made of nanoparticles with different LSP frequencies. Such an approach is relevant for the description of bipartite arrays where the nanoparticles in the two sublattices are either made of different metals, or they are intentionally designed to have different size or shape. In this case we highlight the emergence of gapped chiral Dirac plasmons characterised by an energy-dependent Berry phase. In the last part of the paper (section 4) we explore the interesting scenario of bipartite lattices that break the angular symmetry of the honeycomb array and study the evolution of the aforementioned phase diagram. We show that a sequence of topological phase transitions occurs in the phase diagram while progressively shifting one sublattice with respect to the other. This analysis highlights the existence of massless Dirac plasmons in a vast family of bipartite lattices beyond the honeycomb one, but at the same time it identifies critical values of lattice displacements that lead to the progressive disappearance of massless Dirac plasmons into gapped modes. Our conclusions are finally presented in section 5. A short description of the videos accompanying our paper is given in Appendix A.

Perfect honeycomb lattice of identical nanoparticles
In this section we set the scene for the rest of the paper by revisiting the properties of an ensemble of identical spherical metallic nanoparticles forming a 2D honeycomb struture [44], and conducting a thorough analysis of the resultant CPs.

Model
The honeycomb array has a lattice constant a (see figure 1) and is embedded in an (effective) medium with dielectric constant m . Each individual spherical nanoparticle of radius r supports an LSP resonance, which can be triggered by an oscillating external electric field, with wavelength λ r. Under such a condition, the collective electronic excitation associated with the LSP can be modelled as a point dipole with a natural frequency of oscillation ω 0 that depends on the size, shape and structure of the nanoparticles, as well as on their environment [6]. The dipolar LSPs have a dipole moment p = −eN e h(R)p, where e is the electronic charge, N e is the number of valence electrons in each nanoparticle, h is the displacement field associated with the electronic centre-of-mass motion that corresponds to the LSP at position R, and p is the dipole unit vector. The LSP can be regarded as a bosonic mode, in particular when the size of the nanoparticle is small enough such that quantum size effects are important [47][48][49][50][51][52][53][54].
The nature of the coupling between LSPs in different nanoparticles depends on their size and inter-particle distance. The interactions between the dipolar LSPs lead to collective plasmonic modes that extend over the whole system. Under the condition that r a/3, we consider each dipole to interact with its neighbours through dipole-dipole interactions and only take into account the near field generated by each dipole [17,21]. Such a quasistatic approximation is reliable if λ a and has been shown to qualitatively reproduce the results of more sophisticated simulations which include retardation effects [26]. Within this approximation, the interaction between two dipoles p and p , located at R and R respectively, reads with n = (R − R )/|R − R |. As the CPs will be triggered by an external radiation of a given polarisation, we assume that all the dipoles point along the same direction and parametrise their orientations asp = sin θ(sin ϕx − cos ϕŷ) + cos θẑ, where θ is the angle betweenp andẑ and ϕ is the angle between the projection ofp in the xy-plane and e 1 (see figure 1). Hence, the Hamiltonian reduces to a system of coupled oscillators, where the non-interacting term H 0 reads as [50,52] while the dipole-dipole interaction term is given by [44] H int = (eN e ) 2 m a 3 Here, s = A, B is the sublattice index, h s (R s ) is the displacement field associated with the LSP at position R s of the lattice, Π s (R s ) is the conjugated momentum to h s (R s ), M = N e m e is the mass of the displaced valence electrons (with m e being the electron mass), and the vectors e j (j = 1, 2, 3), defined in figure 1, connect the A and B sublattices. Finally, the coefficients appearing in the Hamiltonian (4) depend on the polarisation angles (θ, ϕ) and hence represent tunable interaction parameters that crucially influence the CP band structure. Note that in equation (4) we only consider the interaction between nearest neighbours, as the effect of interactions beyond nearest neighbours has been shown not to qualitatively change the plasmonic spectrum [44]. The analogy between our plasmonic structure and the electronic properties of graphene becomes apparent by introducing the bosonic ladder operators which annihilate an LSP on a nanoparticle located at position R belonging to the A or B sublattice, respectively. These operators obey the commutation relations , and give access to the CP dispersion and to the nature of the CP quantum states. The bosonic operators above can be converted to momentum space through a R = N −1/2 q exp(iq · R)a q and b R = N −1/2 q exp(iq · R)b q , with N the number of unit cells of the honeycomb lattice. Thus the components of the Hamiltonian (2), equations (3) and (4), become respectively. Here we introduced Ω = ω 0 (r/a) 3 (1 + 2 m )/6 m ω 0 , while the information on the LSP polarisation is encoded in the function

Exact diagonalisation
To find the normal modes of the system we introduce the new bosonic operators and impose that the Hamiltonian (2) (with H 0 and H int given in equations (8) and (9), respectively) is diagonal in this new basis, The Heisenberg equation of motion β ± q , H = ω ± q β ± q then leads to the eigenvalue problem This procedure yields the CP dispersion and the coefficients in the coherent superposition (11) with It is evident by the form of f q in equation (10) that the dispersion of the two CP branches ω ± q crucially depends on the dipole orientation (θ, ϕ). In what follows, we will explore several relevant cases that highlight the tunability of the CP band structures with the LSP polarisation.
Since Ω ω 0 , the CP dispersion in (14) can be approximated as and the coefficients in equations (16) and (17) reduce to cosh ϑ ± q 1 and sinh ϑ ± q 0, respectively, so that the Bogoliubov operators in (11) take the simpler form Incorporating these expansions in the Hamiltonian (12), we find that the system is approximately described by the Hamiltonian The expression above demonstrates that the non-resonant terms in the interaction Hamiltonian H int (cf. the third and fourth terms in the right-hand side of equation (9)) are irrelevant for the description of CPs in the honeycomb array of nanoparticles.

Dirac-like collective plasmons
When the LSP polarisation points normal to the plane (θ = 0), the interaction parameters (5) are all equal to 1. This case is directly analogous to the electronic tight-binding problem in pristine graphene, where the hopping parameters between nearest neighbour carbon atoms are all equal. In fact, the spectrum given in (14) presents two inequivalent Dirac cones [44], similar to those present in the electronic band structure of graphene [29,32]. These Dirac cones occur at a frequency ω ± q = ω 0 and are centred at the K and K points located at ±K D = 4π 3 √ 3a (±1, 0) in the first Brillouin zone (see figure 1(b)). Close to the Dirac points, the function f q defined in (10) expands as f q − 3a 2 (±k x + ik y ), where k is the wavevector away from the Dirac point (q = ±K D + k, where |k| |K D |), such that the CP dispersion (14) is linear and forms a Dirac cone, ω ± q ω 0 ± v|k|, with group velocity v = 3Ωa/2. ‡ We hence find with equation (20) that the system is described close to the Dirac points by the effective Hamiltonian Here, 1 corresponds to the 4 × 4 identity matrix, τ z to the Pauli matrix acting on the valley space (K/K ), σ = (σ x , σ y ) is the vector of Pauli matrices acting on the sublattice space (A/B), and a k,K(K ) and b k,K(K ) annihilate plasmons with wavevector k in the vicinity of the K (K ) valley in the A and B sublattices, respectively. Up to a global energy shift of ω 0 , equation (21) corresponds to a massless Dirac Hamiltonian for the CPs and shows that the corresponding spinor eigenstates represent Dirac-like massless bosonic excitations that present similar effects to those of electrons in graphene [44].

Phase diagram for arbitrary localised surface plasmon polarisations
We now depart from the specific case of LSPs polarised perpendicular to the plane (θ = 0), where there exists a one-to-one correspondence between graphene and our metasurface of metallic nanoparticles. Instead we explore the rich phase diagram of gapless and gapped CP band structures that emerges when tilting the LSP polarisation.
For an arbitrary polarisation of the LSPs, we can determine if the CP dispersion is gapless by imposing |f q | = 0 in (14). Together with equation (10), this leads to the condition [44] 0 for having gapless plasmonic modes. This allows us to produce a phase diagram in figure 2(a) which shows the polar and azimuthal polarisation angles θ and ϕ, respectively, for which the CP band structure is gapless (white regions) or gapped  (correspondingly labelled) in figure 2(a). In figure 2(a), the black solid, dashed and dotted lines indicate the angles where one of the nearest-neighbour coupling strengths C j defined in (5) equals zero, reducing the system to a collection of non-interacting 1D chains. This condition renders the system equivalent to waveguides with a dispersion that is translationally invariant along one direction [44]. As can be seen from figure  2(a), there are also points where two of these lines intersect, signalling that two nearestneighbour bonds are "cut" and the system can effectively be described as isolated dimers, leading to flat bands with no wavevector dependence. The CP dispersions that arise upon continuously tuning the polarisation angle along these lines is shown in video tunable CP dispersion.mp4. In what follows, we predict Dirac physics in these gapless regions, even in pockets of the phase diagram topologically disconnected from the "trivial" graphene case for which θ = 0 (see section 2.3). To explore this, we expand the Hamiltonian in these regions. We have investigated three different gapless positions: (i) (θ, ϕ) = (π/2, π/2), (ii) (θ, ϕ) = (π/4, π), and (iii) (θ, ϕ) = (π/6, 0), as indicated in figure 2(a). In all three cases there are gapless modes with two cones centred at the two inequivalent Dirac points ±K D . In cases (i) and (ii), these are located at ±K D = 2 √ 3a arccos ( 2 5 )(±1, 0), and for case (iii), ±K D = 2 √ 3a arccos (− 2 13 )(±1, 0). Close to the two inequivalent Dirac points, the function f q of equation (10) These expansions show that the CP dispersion (14) forms elliptical cones in the vicinity of the Dirac points above, since the magnitude of the k x and k y components of f q are not equal, unlike in the purely out-of-plane polarisation case (see section 2.3). This is illustrated in figure 3 for case (i). Moreover, by expanding equation (12) in the vicinity of the Dirac points, an effective Hamiltonian can be identified for each case that adequately describes the CPs and, up to a global energy shift of ω 0 , corresponds to a massless Dirac Hamiltonian. Thus, Dirac-like physics can be recovered in other gapless regions away from the case where the polarisation points normal to the plane (θ = 0).
The emergence of a gap in the CP spectrum occurs by the annihilation of Dirac points with opposite Berry phases while varying the polarisation direction. This is illustrated in figure 4 which shows the position of the Dirac points when continuously varying the polarisation angle from ϕ π/10 to ϕ π/5 while keeping θ = π/2 constant (purely in-plane polarisation). These two extrema are denoted by A and B in the inset of the figure, which reproduces part of the phase diagram of figure 2(a). As can be seen from figure 4, the Dirac points move in momentum space as ϕ is increased (see the arrows in the figure that indicate the direction of the motion). At the angles where a gap opens, the two inequivalent Dirac points merge and annihilate each other. Hence in the gapped regions there are no longer two distinct valleys created by the two sublattices. At the angle labelled A in the inset of figure 4 (θ = π/2, ϕ π/10), the gap in the spectrum closes and the two Dirac points in each Brillouin zone appear at the same position, for example two at point (i), two at point (ii) and two at point (iii) in figure 4. At the angle labelled B in the inset of the figure (θ = π/2, ϕ π/5) the two 180 π, where the dispersion for a given ϕ either changes from gapped to gapless or vice versa when θ is increased. In figure 5 we investigate for ϕ = 0 the dispersion at θ values just below θ c ( figure 5(a)), exactly at θ = θ c ( figure 5(b)), and just above θ c (figure 5(c)). It is evident from the figure that two inequivalent Dirac cones annihilate at q = 0 when θ = θ c , subsequently opening a gap and forming two paraboloids for θ > θ c . Indeed, this becomes transparent when expanding (10) in the vicinity of q = 0. Introducing θ = θ c + δθ with 0 δθ θ c , we find Thus, any finite positive δθ leads to the formation of a gap in the CP spectrum at q = 0 (see equation (14)). In the vicinity of this point, the bands are parabolic with two inequivalent effective masses whose ratio is controlled by δθ (see figure 5(c)). In this case, the gap in the spectrum is due to the non-vanishing f q at q = 0 (see equation (23)) and does not require breaking of the symmetry between the two sublattices (see section 3). The phenomenon described in figures 4 and 5 is reminiscent of Dirac point merging in graphene that has been predicted to occur with large mechanical deformations of the lattice [55,56]. While these deformations seem to be impossible to reach experimentally in real graphene, other systems such as cold atoms in optical lattices [41] or microwaves in artificial deformed honeycomb structures [42,43] have been realised and Dirac point merging has been observed. The experimental feasibility of our proposal to observe this phenomenon only requires the external light polarisation to be varied.

Honeycomb structure with broken inversion symmetry
In this section we consider a perfect honeycomb array of spherical metallic nanoparticles with broken inversion symmetry by envisaging that the LSPs on the A and B sublattices have inequivalent resonance frequencies, ω A and ω B , respectively. This could be realised experimentally by either manufacturing the two sublattices out of different materials or by constructing them with different sizes [6]. For our quasistatic approximation of point-like interacting dipoles to hold we still assume the radii (r A and r B ) of the particles in the two sublattices to be much smaller than λ. The analysis in section 2 applies to this case as well with only minor changes: (i) the natural LSP frequency ω 0 is replaced by ω A and ω B in the two sublattices, and (ii) the interaction coefficient Ω is replaced byΩ = (ω A ω B /ω 2 0 ) 1/2 (r A r B /r 2 ) 3/2 Ω. As a result, the dispersions of the two CP branches now read where f q is defined in equation (10). It is convenient to express the two LSP frequencies as ω A =ω + δω/2 and ω B =ω − δω/2, withω the average frequency and δω the difference in frequency. § In figure 6 we plot the dispersion (25) in the special case where the polarisation points perpendicular to the 2D plane (θ = 0), for δω/ω = 0 (figure 6(a)), δω/ω = 0.004 ( figure 6(b)) and δω/ω = 0.04 (figure 6(c)). As can be seen from the figure, any finite difference in the LSP frequencies (δω = 0) introduces an asymmetry into the system that corresponds to a gap of size δω opening in the spectrum while the extrema of the two bands still occur at the K and K points in the Brillouin zone. Notice also that the symmetry-breaking term δω affects the bandwidth of each plasmonic subband, defined as W ± = max {ω ± q }−min {ω ± q }, as shown in figure 7. In particular, increasing δω leads to an algebraic decrease of both bandwidths. These effects of the symmetry-breaking term δω are further examplified in the video broken inversion symmetry.mp4. In the vicinity of the Dirac points, the system is effectively described by the Hamiltonian for gapped quasiparticles with eigenvalues ω ± k =ω ± ṽ 2 |k| 2 + (δω/2) 2 , whereṽ = 3Ωa/2. In contrast to the analysis in section 2.4, here the dispersion is gapped for any value of wavevector and the quasiparticles acquire a finite effective mass m * = δω/2ṽ 2 due to the inversion symmetry breaking term δω = 0. This scenario is reminiscent of the electronic dispersion in hexagonal boron nitride [57] and in transition metal dichalcogenides, such as MoS 2 [58][59][60].
Besides the global frequency shiftω, the effective Hamiltonian (26) in each valley evolves continuously between a purely out-of-plane Zeeman term ( δω/2)σ z (wheñ v|k| δω) and a 2D massless Dirac Hamiltonian − ṽσ · k (whenṽ|k| δω). The eigenstates of the Hamiltonian in the K valley thus correspond to the normalised vector in the Bloch sphere V ± k = ±[ṽ 2 |k| 2 + (δω/2) 2 ] −1/2 (−ṽk x , −ṽk y , δω/2). This allows us to calculate the Berry phase of the pseudo-spin s = 1/2 Dirac quasiparticles described by (26) as φ B = sΩ B , with Ω B the solid angle enclosed by the Bloch-sphere vector V ± k while the electronic state |k is transported anticlockwise in a closed loop around the Dirac point in the 2D wavevector space [37,61,62]. Depending on the ratiõ v|k|/δω, and thus on the quasiparticle energy, the solid angle Ω B changes (see the zoom in figure 6(b)). This ranges from Ω B = 0 at the band gap edges (ω =ω ± δω/2 where V ± k=0 = ±ẑ) to Ω B 2π ifṽ|k| δω/2 (where V ± k rotates in the x-y plane), the latter case being analogous to that found for electrons in monolayer graphene [30,31].
Interestingly, we find that CPs in honeycomb arrays with broken inversion symmetry are naturally described as gapped Dirac quasiparticles with an energydependent Berry phase (defined up to integer multiples of 2π) φ B (ω) = π |ω −ω| − δω/2 |ω −ω| (27) in the K valley if |ω −ω| δω/2. A similar analysis yields an equal and opposite Berry phase for the K valley.
The explicit connection between pseudo spin and orbital degrees of freedom evident in the effective Hamiltonian (26) is at the very heart of several properties exhibited by quasiparticles in real and artificial graphene systems, such as the suppression of elastic backscattering from smooth impurity potentials [33], weak antilocalisation [63,64] as well as the anomalous quantum Hall effect [30,31]. In a forthcoming publication [65] we will analyse in detail the consequences of the energydependent Berry phase (27) in quantum transport of gapped Dirac quasiparticles.

Deformed bipartite lattices
The four-component wavefunction of the quasiparticles, along with the effective massless Dirac Hamiltonian in real and artificial graphene stem from the bipartite nature of the 2D lattice as well as from the time-reversal symmetric and parityinvariant nature of the system [66]. In particular, the honeycomb structure is a special case of the bipartite hexagonal lattice boasting the additional invariance under rotation by integer multiples of 2π/3. In this section we explore the band structure and the stability of the phase diagram of gapless and gapped spectra for CPs in generic hexagonal bipartite lattices obtained by rigidly shifting one sublattice with respect to the other (see figure 8). Our analysis highlights the possibility to design a variety of tunable plasmonic metamaterials with the desired polarisation-dependent optical response, along with the additional benefit of supporting Dirac CPs.
We investigate the properties of 2D arrays of identical spherical metallic nanoparticles with LSP frequency ω 0 sharing the same hexagonal Bravais lattice as the honeycomb structure considered so far, but with an arbitrary position of the second basis nanoparticle. Thus the lattice can still be considered as made of two hexagonal sublattices A and B, but with the B sublattice shifted by a vector such that the new nearest-neighbour vectors now become e j = e j − (see figures 1(a) and 8). This rigid shift breaks the three-fold rotational symmetry of the original honeycomb case. The main questions then to be addressed concern the very existence of Dirac CPs and the fate of the phase diagram while progressively distorting the original lattice.
The mathematical procedure to obtain the Hamiltonian and the band structure for CPs in the general bipartite hexagonal system is a straightforward extension of that presented in section 2. The replacement of the new nearest neighbour vectors e j in the Hamiltonian (4) leads to the dispersion relation with Here, C j = (a/|e j |) 3 (1 − 3 sin 2 θ cos 2 ϕ j ) is a generalisation of C j defined in equation (5), with ϕ j = ϕ + arctan(e j,x /e j,y ) the angle between e j and the projection of the LSP polarisationp into the x-y plane. Unlike in the honeycomb case, the lengths of the nearest neighbour vectors are not all equivalent, hence the extra term (a/|e j |) 3 in C j . An equivalent analysis to that in section 2.4 for the new bipartite lattices yields the resulting phase diagrams which identify the domains of LSP polarisation leading to the existence of Dirac CP quasiparticles. From equation (29) we obtain the condition for gapless plasmonic dispersions at an arbitrary polarisation, Figure 9 shows the resulting phase diagrams in the (θ, ϕ) parameter space for five different lattices, where the cream and orange regions indicate gappless and gapped CP spectra, respectively. While figure 9(a) corresponds to the undeformed honeycomb lattice for which = (0, 0) (see also figure 2(a)), figures 9(b)-(e) are obtained by progressively displacing the B sublattice in the −ŷ direction (see figure 8). The phase diagram in figure 9(b), obtained for = (0, −0.0893a) is topologically equivalent to the one of the undeformed lattice (figure 9(a)). Despite the rather small displacement (less than 10 %), the broken discrete rotational symmetry of the lattice induces an evident bulging of the phase diagram compared to the undeformed case. For larger deformations, in figures 9(b) to 9(e) we observe that the topology of the phase diagram changes, as indicated by the arrows in the figure. For example, when going from = (0, −0.0893a) to = (0, −0.1433a) (figure 9(b) to 9(c)), a compact gapless domain splits into two disconnected pockets separated by a gapped phase. In contrast, when going from = (0, −0.1433a) to = (0, −0.1437a) (figure 9(c) to 9(d)), the topological phase transition occurs by annihilation of a gapless domain into a gapped one. Finally, when going from = (0, −0.1437a) to = (0, −0.1967a) (figure 9(d) to 9(e)), the transition occurs by merging three topologically disconnected gapless domains into one. Even larger deformations lead to a phase diagram that is mostly dominated by gapped phases.
The five phase diagrams shown in figure 9 are snapshots from video deformed bipartite lattice.mp4, which shows a large succession of phase diagrams corresponding to lattice deformations from = (0, 0) to = (0, −0.2a). As is clear from the video and can be inferred from panels (c) and (d) in figure 9, the spectrum for the polarisation θ = 0, for which the analogy between our CPs and electrons in graphene holds [44], is gapped for deformations larger than (0, −0.1435a). Whilst we are not discussing elastic deformations of real graphene membranes [67][68][69] it is insightful to note that the literature on strained graphene suggests that a ca. 20 % change in one of the nearest-neighbour hopping lengths is needed to merge the two Dirac points and open a gap in the spectrum [55,56,70]. Figure 10 further illuminates the qualitatively distinct physical regimes realisable in different lattices for the same dipole polarisation. For θ = 0, figure 10(a) shows the CP band gap ∆ in units of the coupling Ω as a function of the position of the B sublattice (white regions indicate a gapless CP dispersion). For a perfect honeycomb array, the B sublattice is located at the center of the triangle (black dot) defined by the three neighbouring nanoparticles in sublattice A (depicted as blue circles in the figure). The grey areas indicate regions for which the interparticle distance is smaller than 3r and where the dipole-dipole approximation (1) breaks down such that higher multipoles need to be taken into account. As can be seen from the figure, the gapless dispersion is robust in a sizeable region around the point = (0, 0) (see dot in the figure). While shifting the B sublattice towards the A one, a non-vanishing gap forms, whose magnitude progressively increases.
To emphasise that indeed the gapless domain in figure 10(a) supports Dirac CPs, we focus on the particular displacement = (0, a/4) indicated by a cross in the figure.
Note that the discussions hitherto have considered nearest-neighbour interactions only, an approach shown to qualitatively capture the relevant physics for the perfect honeycomb structure [44]. One must however be cautious in wielding such a model for too large values of where further neighbour terms in the interaction become increasingly important.
For such a bipartite hexagonal lattice, the two inequivalent Dirac points are located at ±K D = ± 2 √ 3a (arccos (−   √  2197 250 ), 0), close to which the function f q defined in equation (29) expands as f q − 96a 125 (± 20101 2197 k x + ik y ). Thus, close to the Dirac points the CP dispersion describes massless Dirac CPs with an angular-dependent group velocity, as shown in figure 10(b).

Conclusion
We have studied the tunable plasmonic response of an artificial analogue of graphene constituted by a 2D hexagonal bipartite lattice of metallic nanoparticles. We have analysed the band structure of collective plasmons stemming from the near field dipolar coupling between localised surface plasmons in individual nanoparticles. In the case of a perfect honeycomb array, we have explored in detail the rich phase diagram of gapless and gapped phases emerging from tilting the polarisation of the dipole moments. We have shown that each topologically disconnected gapless phase supports collective plasmons that are effectively described by a massless Dirac Hamiltonian. The emergence of gapped phases from gapless ones occurs via the coalescence of Dirac cones while progressively tilting the dipole orientation, in analogy to what would occur for electrons in real graphene under extreme elastic strain. For perfect honeycomb arrays with two inequivalent sublattices we unveil that collective plasmons are effectively described as gapped Dirac particles with an energy-dependent Berry phase.
We also studied the plasmonic properties of bipartite hexagonal arrays obtained from the honeycomb structure by rigidly shifting one sublattice with respect to the other. This system is still highly tunable with the orientation of the dipoles and supports gapless Dirac plasmons as well as gapped phases. The phase diagram undergoes significant changes for lattice displacements of a few percent, with a series of topological phase transitions that split and annihilate gapless Dirac phases. Our analysis shows that a large family of bipartite arrays of metallic nanoparticles can support Dirac plasmonic modes. These could efficiently transport highly confined radiation through the 2D metamaterial. While gapless Dirac phases exist in generic hexagonal bipartite lattices, they occupy a progressively small portion of the phase diagram for increasing deformations. This analysis will be important for designing the properties of novel plasmonic metamaterials based on bipartite lattices, e.g., by assembling dimers in ordered structures or by means of self-assembled materials with less symmetry than the honeycomb structure.
The collective plasmons in our work constitute the building block for the analysis of the optical properties of Dirac metamaterials. This would involve a quantum treatment of the coupling between photons and collective plasmons resulting in tunable plasmon polaritons that may inherit some of the effective Dirac-like properties studied in the present work. A preliminary work in this direction [71] has shown that plasmon polaritons in a simple 3D cubic array of metallic nanoparticles present a band structure that is tunable with the polarisation of light, leading to birefringence which is purely due to interaction effects between dipolar localised surface plasmons.
A crucial question to be addressed in the quest of designing arrays of metallic nanoparticles supporting Dirac-like collective plasmons is the inherent damping that these collective modes suffer from. Localised surface plasmons in individual metallic nanoparticles are subject to three main sources of dissipation [6]: absorption (Ohmic) losses, that do not depend on the size of the nanoparticle, radiation losses that scale as the volume of the nanoparticle, and Landau damping, a purely quantum-mechanical effect whose associated decay rate is inversely proportional to the nanoparticle radius. The two latter damping mechanisms suggest that there exists an optimal nanoparticle radius for which the total damping rate is minimised, estimated in reference [44] to be around 8 nm for silver nanoparticles. The interaction between the localised surface plasmons in each nanoparticle forming the bipartite lattice may change the picture above. Indeed, recent theoretical results on metallic nanoparticle dimers [72] suggest that while Landau damping is weakly influenced by the dipolar interaction, radiation damping strongly depends of the wavelength of the collective mode. Further work will be needed in order to investigate this important issue.