Topological structures, spontaneous symmetry breaking and energy spectra in dipole hexagonal lattices

The interplay between the special triangular/hexagonal two dimensional lattice and the long range dipole–dipole interaction gives rise to topological defects, specifically the vortex, formed by a particular arrangement of the interacting classic dipoles. The nature of such vortices has been traditionally explained on the basis of numerical evidence. Here we propose the emerging formation of vortices as the natural minimum energy configuration of interacting (in-plane) two-dimensional dipoles based on the mechanism of spontaneous symmetry breaking. As opposed to the quantal case, where spin textures such as skyrmions or bimerons occur due to non-linearities in their Hamiltonian, it is still possible to witness classic topological structures due only to the nature of the dipole–dipole force. We shall present other (new) topological structures for the in-plane honeycomb lattice, as well as for two-dimensional out-of-plane dipoles. These structures will prove to be essential in the minimum energy configurations for three-dimensional simple hexagonal and hexagonal-closed-packed structures, whose energies in the bulk are obtained for the first time.


Motivation
Two classic identical dipoles of magnitude of µ , m u and m v localized, respectively, at positions r u and r v , with dipole moments given as interact with each other giving rise to the interaction energy where r uv is the vector between the two dipoles u and v, r uv = ||r uv || is their separation distance (constant C E is either µ 0 4π or 1 4πǫ 0 , magnetic or electric, respectively). Energy (2) will be given in units of C E µ 2 /a 3 throughout the present work.
The terms in (2) possess a well-defined physical meaning. On the one hand, the first term is essentially the classical counterpart of the Heisenberg exchange interaction, which plays a paramount role in magnetism. On the other hand, the second term favors the alignment between the two dipoles. In a particular configuration of dipoles in equilibrium, there is a non-trivial interplay between the two contributions, which, depending on the particular lattice, does not resemble the usual Heisenberg interaction. Furthermore, the Heisenberg Hamiltonian is computed for nearest-neighbors. What will make our study specially challenging is that in the case of interaction (2), one will consider all pairs, as opposed to just neighbors as in the usual Heisenberg case, and the coupling constant will not be a constant at all, but shall possess a long-range nature (going as 1/r 3 ).
Emerging fields of research such as chiral magnetism and topological spintronics are developed through the study of topologically nontrivial spin textures. These structures are found in a variety of magnetic materials [25][26][27][28][29][30][31][32][33][34] , and are potentially interesting for information processing and data storage. As opposed to the classic interaction energy functional form (2), the Hamiltonian for quantum spins, running over nearest-neighbor, next-nearestneighbor, and next-next-nearest-neighbor sites, contains additional non-linear terms that are responsible for the appearance of magnetic skyrmions in 3D or magnetic bimerons in 2D 35 . Also, these topological magnetized spin textures carry integer topological charges Q. In essence, it is the existence of non-linearities in the quantum Hamiltonian that lead to well-defined topological structures.
In spite of the stark contrast existing between the nature and form of the interactions in quantum and classical systems, one may wonder if it is still possible to observe special topological structures in classic systems of interacting dipoles. The answer is affirmative in systems with hexagonal symmetry, either in two or three dimensions. To be more specific, we shall focus on the minimum energy configuration of systems of dipoles in the thermodynamic limit, as opposed to the quantum case. The classical extremal energy configurations of a magnetic dipole system (either minimum or maximum) is one of equilibrium in which no torque should act on any given dipole. Now, the total energy of a system of dipoles in terms of (2) can also be cast in the form of the Hamiltonian The quadratic form (3) will be particularly suitable in the next section.

Luttinger-Tisza and energy decomposion methods
We can obtain the spectrum of a system of identical dipoles using the Luttinger-Tisza method 17 under the assumption that the minimum energy configuration exhibits translational symmetry. If T(i) denotes the points generated from i with discrete translations belonging to the T symmetry group the mentioned symmetry corresponds to m i = m i ′ for all i ′ ∈ T(i) . The system can be split into identical cells and thus limit the summation to one single cell. Therefore, the energy per dipole can be expressed as The energy per dipole finally reads as The method involves the sotution of an eigenvalue problem of the nd dimensional matrix Â (d is the dimension of the dipole), with k being the eigenvalues and x k the corresponding (orthogonal) eigenvectors of the system, with �x k � = √ n . The final expression for the energy reads as where a k denotes the components of m in the base {x k }. Two conditions must be satisfied for i = 1 . . . n , namely Within the framework of the Luttinger-Tisza method, these two constraints are known as the strong and the weak conditions, respectively. We can thus obtain the minimum energy per dipole as E min = 1/2 min µ 2 , where min denotes the smallest eigenvalue of Â . There exists an alternative to the previous Luttinger-Tisza method, derived by us [36][37][38] . For a certain set of positions {r 1 , r 2 , . . . , r M } corresponding to M identical dipoles (all of them with the same magnitude µ ), one resorts to a minimization procedure for its total dipole-dipole interaction energy. We chose the well-known simulated annealing 39 (SA) method. This Monte Carlo method is able to scape from local minima in the space of M(2M) variables corresponding to 2D(3D). We have checked that even for high values of the total dimensionality of the problem, it returns precise values for the optimal energy.
For a finite sample of N dipoles it is then easy to infer the minimum energy configuration. Once it is known, one needs to "grow" the plaquette in the most symmetric way, that is, equally in all possible directions. In the same spirit of the liquid drop model 40 , which considers the total nuclear energy as a sum of different contributions, such as the volume or the surface term, we shall correspond each inferred total energy E N to a series of functionally different terms depending on N, which fixed coefficients to be obtained later on.
Therefore, the total energy is decomposed as The term ∝ N is the volume contribution, in other words, the asymptotic energy per dipole in the thermodynamic limit. That is, V E ≡ lim N→∞ E N (N)/N . The second term in (9) corresponds to the surface contribution. Similarly, the third term in (9) takes into account those dipoles on the boundary ∂� . Finally there is an independent contribution. By adjusting the numerically obtained total dipole-dipole plaquette energy to the functional form (9), which employs the so-called Levenberg-Marquardt 41 non-linear regression, we do obtain remarkable results. In this fashion, we are able to infer the value of lim N→∞ E N (N)/N with great precision.
The previous overall optimal configuration − sample growing scheme, constitutes a easy-to-implement alternative to Luttinger's because (1) no periodic boundary conditions are imposed, nor (2) the computation of Ewald sums (given in terms of involved lattice sums 42 ) are required for obtaining the minimum energy per dipole in the limit N → ∞.
Also, as a general rule in any two-dimensional periodic configuration of dipoles (not necessarily optimal ones), we have discovered that every global rotation of all dipoles by a phase θ is followed by a linear dependency of the total interaction energy on the cosine of twice the phase. That is, E(θ) = a + b cos(2θ) . Heuristically, this relation can be obtained as a generalization of the 1D case. Suppose that all dipoles along an infinite straight line are rotated clockwise by an angle θ . The total (minimum) energy per dipole in the thermodynamic limit will be given by E(θ) = (1 − 3 cos 2 θ)ζ (3) , where ζ(3) is the so called Apéry constant. By using the relation between the cosine of an angle and half that angle, we obtain E(θ) = −1/2 ζ(3) − 3/2 ζ(3) cos(2θ).

Two-dimensional structures
Hexagonal lattice: out-of-plane dipoles. Let us first consider an infinite hexagonal layer where dipoles can freely rotate out-of-plane, yet always being fixed in space. The first state that one could think of is the maximum energy state, which corresponds to all dipoles either point upwards or downwards. This instance was first considered by Danilov et al. 43 , where a direct calculation of the total energy per particle in terms of lattice sums was mentioned to be divergent, which is not. In point of fact, we do obtain using our energy decomposition method TÂm . Interestingly, the minimum energy configuration when we grow the crystal of dipoles in hexagonal plaquettes corresponds to a special arrangement of dipoles as shown in Fig. 1. Crosses represent dipoles pointing inwards, whereas dots correspond to dipoles coming out of the plane. We call this configuration the "flag" state, due to its resemblance with the South African flag. Notice the special arrangement of the dipoles, that makes this state posses a point symmetry of 120 • . This state is certainly unique from the topological point of view. Also, we notice that it is highly sensitive to the boundary ∂� . The transition from the flag state Fig. 1a to a new state whose boundary is square-like Fig. 1b is also shown. This configuration possesses a higher energy, namely, −0.883517 ± 1.852 × 10 −10 . When the lattice is further distorted to form a parallelogram, we recover 1/6th of the original state in Fig. 1c. It is obvious that the shape of the boundary affects the total energy per dipole all the way to the thermodynamic limit, a feature which is not blurred by finite-size effects. We currently can only offer energetic numeric arguments in order to explain the existence of the flag state, but no physical ones. Notice that the high sensitivity of the flag state to the boundary conditions is an inherently typical feature of topological states in the quantum sense. In other words, in the quantum realm, it is very difficult to rigorously prove how the boundary ∂� affects the robustness of the corresponding enclosed state, and that is no different classically.
All that is known to happen is that the flag state is not topologically equivalent to the other two configurations. There is a clear symmetry reduction from going to the flag state to the others, an effect that has energetic consequences only in the square-like confined one, the one in Fig. 1b. It is likely that first and second configurations (a) and (b) in Fig. 1 tend to the last one (c) there in the bulk.
Hexagonal lattice: in-plane dipoles. The Luttinger-Tisza method was applied to a system of twodimensional dipole moments with identical scalar strength located at the sites of an infinite rhombic lattice with an arbitrary rhombicity angle by Brankov and Danchev 44 . This two-dimensional cell choice allowed to study the simple square lattice and the hexagonal one in a very compact way. Their study was an extension of two previous works by Belobrov et al. 45,46 that considered finite clusters of dipoles, highlighting the existence in the hexagonal case of a "macro-vortex". Brankov and Danchev 44 extended their original work and provided a better numerical accuracy to the computation of configuration energies. Although further work on dipole-dipole interacting systems followed 47-55 inspired by the Luttinger and Tisza seminal work 17 , Refs. [44][45][46] constitute the first consistent treatment of simple square lattice of dipoles and, specifically, the hexagonal lattice one. We shall provide an extended description of the hexagonal case in the present work, with an increased precision of the energy per dipole in the thermodynamic limit by several orders of magnitude. Figure 2 depicts the set of possible states obtained by either employing the Luttinger-Tisza method, or using numerical computations based on symmetry arguments. Each configuration is shown, along with their energy per dipole energies, given with exact decimal digits. Incidentally, the antiferromagnetic states of energies − 2.04746 and 2.96697 are connected by a global rotation in the sense of the E(θ) = a + b cos(2θ).
Let us consider the case of the minimum possible energy configuration given by the vortex. The Luttinger-Tisza method provides an energy given by the lattice sum where L 2 stands for the points in the hexagonal lattice. There are means to speed-up the computation of lattice sums such as (10) by using Ewald sums.
What we would like to stress out is that we shall compute the very same value for (10) without the burden of even considering the Luttinger-Tisza method at all, with all the corresponding formalism. We shall, then, show how our energy decomposition method works. Incidentally, we consider a (10) Figure 1. Depiction of the minimum energy configuration for perpendicular dipoles in a single layer, as we distort the lattice on which we grow the dipole crystal. From an hexagonal shape (a), whose ground state is the flag state with energy -0.919515, we obtain another minimum energy configuration where the shape is square-like, as shown in (b) www.nature.com/scientificreports/ slight deviation of (9) of the form N V E + N 1/2 S E + C , where the boundary term is not accounted.
Recall that we are only interested in the bulk value V E , and this omission shall not imply any variation in the final result whatsoever. Once we have obtained several hexagonal clusters of increasing size N their corresponding energies (according to Fig. 2a or 2b), the concomitant non-linear fit returns the values V E = −2.758545 ± 8.223 · 10 −9 , S = 1.22402 ± 7.662 × 10 −6 , C = 3.30666 ± 0.001497 . Extraordinary precise values for E 0 (10) as well as for the other energies are thus obtained in the same fashion without employing the tools required for the Luttinger-Tisza method. Let us focus on the description of the vortex state. First of all, it is not expected to have a peculiar configuration of in-plane dipoles as the minimum energy state, having a considerable energy of − 2.758545. As opposed to the flag state, the vortex state is very much robust and, as far as our numerical computations are concerned, does not depend on the boundary ∂� of the system. The vortex state is depicted in Fig. 3b in great detail. Let us analyze the vortex state from the symmetry point of view. The union of three spirals (red, green and blue) has a third-order axis (rotation on 120 • ) perpendicular to the plane of the lattice. This is a point group of symmetry 3. The symmetry group of the hexagonal lattice contains identical transformation, 2-, 3-and 6th order rotations and six reflections. The identical transformation leads to one spiral, shown in the inset of Fig. 3, and rotations lead to the spiral configuration that forms the vortex state.
What is remarkable about the vortex state from the symmetry perspective alone, is that the partition of the hexagonal lattice into spirals is not unique. One can construct the partition on two spirals with the second-order axis (the point group of symmetry 2, rotation of 180 • as in Fig. 3b. Furthermore, one can also obtain the partition of the hexagonal lattice on the central point with 6 spirals with the six-order axis (rotation on 60 • , the point group 6), which constitutes the vortex in Fig. 3a. Of course, all these groups are subgroups of the symmetry group of the hexagonal lattice. The vortex structure and the onset of spontaneous symmetry breaking. Despite the previous description, the question still remains: why a (macro)vortex? The works 45,46 managed to answer the question for finite clusters of dipoles, based on surface magnetic charge arguments. However, when dealing with the corresponding N → ∞ , a more intuitive reason arises: "The structure of a dipole system should be sensitive to the symmetry of the nearest-neighbor environment". However, this claim is rather vague.
One possible explanation could be the following. Let us assume that the minimum energy configuration for low temperatures is the (in-plane) ferromagnetic one in the bulk. The original Hamiltonian (3) has an O(2) dipole rotation symmetry, whilst the direction of the ferromagnetic order is not fixed for N → ∞ . Things, however, change near the edges in finite systems. There there is a clear ferromagnetic direction, and different regions or domains merge from the boundary to the interior, in a sort of inwards dipole "freezing".
The previous argument somehow implicitly assumes the six-fold symmetry of the vortex state as inferred from the lattice. There is no reason why ferromagnetic domains should not merge forming a triangular shape, for instance. Again, it is not quite convincing. In essence, the existence of vortices formed by magnetic dipoles have been suggested mainly by numerical simulations. Analytic Taylor expansion of the dipole-dipole interaction has shown 56 that it is able to obtain a local representation in the plane such that the lattice symmetry induces a magnetic anisotropy. This approximation, although purely mathematical in essence, is somehow capable of predicting the existence of vortices in the ground state magnetization.
Based on the previous facts, we shall pursue an explanation solely based on physical and symmetry arguments. To such an end, we shall first define the chirality of a set of N dipoles {m i } with respect to an arbitrary, perpendicular axis o as where r i denotes the position of the dipole m i with respect to an axis o. The chirality is obviously a quantity independent of the position of the axis. In the case of the vortex state, we obtain the following result with V 1 = 1.73205 ± 1.833 × 10 −7 , V 2 = 7.79359 ± 2.491 × 10 −4 , and V 3 = 11.4368 ± 0.08078 . Therefore, we shall define the chirality of the system of dipoles in the thermodynamic limit as ξ ≡ lim N→∞ ξ N /N 3 , which is equal to ±1.73205 , depending on the vortex pointing clockwise or counter-clockwise.
The magnetization of the vortex state is per se rather intriguing. The three embedded spirals conform an overall state with a nonzero total dipole moment M. Since the number of dipoles N goes as 1 + 3n(n − 1) (formula for hexagonal centered numbers) for each added layer n, and M is easily found to be (2n − 1)µ , we thus have the explicit dependence Therefore the magnetization per dipole of this state goes as 1/ √ N , as opposed to the expected 1/N behavior in the thermodynamic limit. Notice that in the rest of states depicted in Fig. 2, M is either 0 or o(1/N). The (11)  www.nature.com/scientificreports/ magnetization per dipole could be a potential physical magnitude to be used in order to link the minimum energy configuration of the hexagonal lattice, that is, the vortex, with the rest of the spectrum. However, it does not change continuously and remains practically zero for all practical purposes. That is why we shall consider the aforementioned chirality ξ . In the language of phase transitions, this is going to be our order parameter in the Landau sense, although at zero temperature. In order to explain the appearance of a vortex structure as the minimum energy configuration of interacting dipoles in the hexagonal lattice we shall invoke here a mechanism of spontaneous symmetry breaking. Phases of matter, such as crystals, magnets, and conventional superconductors, as well as simple phase transitions can be described by spontaneous symmetry breaking (See 58,59 and references therein). In the case of crystals as periodic arrays of dipoles, for instance, they are not invariant under all translations (only under a small subset of translations by a lattice vector). The Hamiltonian of our system (3), for in-plane dipoles, has a finite rotational symmetry defined by the symmetry group of the hexagonal lattice plus Z 2 , as well as all the corresponding states depicted in Fig. 2, a symmetry that is further reduced for the ground states, the vortex states. The vortex clearly breaks this rotational symmetry. For out-of-plane dipoles, the Hamiltonian is invariant under O(2) rotations around any perpendicular axis. Although the flag state has 120 • degrees symmetry, it is indeed invariant under O (2), as well as the other out-of-plane states. Plus, the vortex configuration possesses a Z 2 symmetry, that connects configurations with opposite chiralities. Notice that Z 2 commutes with the Hamiltonian (3). Numerically, the vortex structure is rather robust, and already occurs for finite instances. It is extremely difficult to prove that the vortex state does not depend on the boundary conditions as we go to the thermodynamic limit, as opposed to the flag state. As previously mentioned, they easily appear in all sort of simulations. This fact does not prove anything by itself. We can only rely on numerical calculations to assert the robustness of the vortex state (recall that this procedure is also the case in the quantum realm).
All the previous states, either out-of-or in-plane dipole configurations, has exactly zero chirality except the vortex structure. Thus, we can consider different routes to the spontaneous symmetry breaking that lead to the formation of vortex states. They are shown in detail in Fig. 4. Path III connects the maximum and minimum energies, by a continuously development of the chirality. The same situation happens for path I, linking two states energetically close in relative terms. Path II is rather special for it brings together the flag state and the The spontaneous symmetry breaking mechanism that we propose here use the chirality as the order parameter, which induces the appearance of a special configuration of dipoles having the minimum possible energy. It does not necessarily imply that the shape has to be that of a (macro)vortex, but at least a different one from the set of possible states possessing a clear, differentiated symmetry.
The honeycomb lattice: topology of the ground state. In the honeycomb lattice, and following our procedure, a sufficiently big sample of dipoles unveils the disposition of the dipoles for the minimum energy configuration. This structure is shown in Fig. 5, along with the lattice. The minimum energy per dipole in the thermodynamic limit is found to be −2.22691 ± 6.322 × 10 −6 , higher than the one corresponding to the hexagonal lattice.
We do not observe a single dipole structure, as in the case of the vortex in the hexagonal lattice, with ferromagnetic domains. Antiferromagnetism is not present either. Instead, the system is formed by an array of local sinks and vortices. For every sink, we have six local vortices that surrounds it. As in all classical ground states with no geometric frustration 57 , it is not highly degenerate (only two) and the total magnetization is zero.
Certainly, we cannot possibly have the topology present in skyrmion/bimeron spin textures, yet we can obtain a periodic arrangement of interesting structures such as sinks and vortices.

Three-dimensional structures
Hexagonal simple lattice. The problem of obtaining the minimum energy per particle in the bulk of a system of interacting dipoles is interesting per se, regardless of the peculiarities of the corresponding dipole topology. Also, to the best of our recollection, it has not been studied previously in the literature.
Again, in order to tackle the problem, we use our methodology. After placing N z ∈ [1, 100] layers of hexagonal lattices in parallel, each one of which with side N x = N y = 20 , we conclude that no dipole projection onto the x-y plane occurs. Of course, there are finite-size effects that account for dipoles than can be found not being perfectly perpendicular, but these effects are blurred in the thermodynamic limit. In point of fact, if the www.nature.com/scientificreports/ 2D layer retains its hexagonal shape, the final state of the entire column corresponds to layers of flag states one on top of the other. On the contrary, if the 2D layers possess some other shape, such as a parallelogram, the crystal will look like a 1/6-prism of the latter. As in the 2D case, we know that the two vertical arrangements will possess the same energy, the same zero magnetization, but the flag hexagonal lattice dipole crystal will display a particular topology, not present in the 1/6-prism. In both cases, either flag layers or a 1/6-prism, we shall have multi-ferromagnetic-domain structures of dipoles pointing up or down along the z-axis. This is the case for any d value in the rhombic cell. Thus, d shall correspond to the interlayer distance, measured in units of the hexagonal lattice constant a. Of course, there is no preferred value d. Therefore, it is mandatory to perform an exhaustive numerical exploration in order to obtain the minimum energy per dipole in the bulk and for every different d value. The results of the computations are depicted in Fig. 6. The small spheres in Fig. 6 and 7 denote the position of each dipole in the crystal, which could easily be regarded as the position of the nuclear magnetic moment corresponding to the nucleus of a certain species. The representation using spheres also helps to better understand the later comparison between different three-dimensional structures in terms of hard spheres.
The total interaction energy monotonically decreases as d decreases, which implies that the system is bounded. The first computed point is for d = √ 2/3a , which corresponds (for later comparison) to the minimum interlayer distance in the case of the hexagonal closed packed lattice, for a = b = 1 . As the interlayer separation increases, the interaction between layers becomes negligible, and we recover asymptotically the expected value for one single layer, namely, − 0.919515. Here we list the computed values for several distances d. Figure 6. Plot of the evolution of the energy per particle in the thermodynamic limit as a function of the vertical cell distance d, for both the flag-state crystal and the 1/6-prism. The horizontal line corresponds to the asymptotic d → ∞ energy per particle value of the flag state. See text for details. Figure 7. Plot of the evolution of the energy per particle in the thermodynamic limit as a function of the vertical cell distance d, for the hcp lattice. The horizontal line corresponds to the asymptotic d → ∞ energy per particle value of the vortex state. See text for details.  to a), the hcp lattice is energetically more bounded, a feature which is somehow intuitively expected.
Our results for the hcp lattice, with vortex layers, are not to be confused with helical spin configurations in hexagonal crystals 60 . On the one hand, the former occur in simple hexagonal lattices, and on the other, their explanation is due to the a positive interaction between nearest neighbors and a negative interaction with secondnearest neighbors. Therefore, rotations in planar layers occur only if spins are treated classically as dipoles in the Heisenberg exchange Hamiltonian −|J 1 | cos α + |J 2 | cos(2α) , where α is the angle of rotation of the dipoles in the layers forming an helix.

Discussion
We have described the emergence of the robust vortex state as the minimum energy configuration possible mediated by a mechanism of spontaneous symmetry breaking. In a way, since T = 0 , we are describing a classic topological transition driven by the continuous chirality ξ parameter.
The real arena where the symmetry breaking should be considered is at finite temperature, and not to employ the internal energy E or U, but the (Helmholtz) free energy F = E − T · S instead (a collection of particles will always seek to minimize its free energy). There, one should take into account that the entropy S clearly scales differently as opposed to the total energy E. Scaling here is referred to the size of the system. In other words, depending on the density of states of the system, the entropy may have a scaling with system size different than the internal energy. Therefore, the ensuing study should become more involved.
However, all the present contribution is described at zero temperature. It the temperature was finite, as it was lowered, the entropy would provide a decreasing contribution to the free energy, and the system would at some point fall into another state. This usual temperature-driven analysis has not been implemented here, for again we are treating the system at zero temperature. However, the Landau approach for finite T may not to be much different from what has been presented here. We are confident that, in a finite temperature analysis, our conclusions about the feasibility of the spontaneous symmetry breaking explanation for the existence of vortex would still hold.
Energies and configurations have been computed either by recourse to the Luttinger-Tisza method, or by employing our energy decomposition approach, which has been proved successful in other scenarios involving the computation of energies based on the dipole-dipole energy interaction of periodic systems in the thermodynamic limit.
The topology of the equilibrium configurations for out-of-plane dipoles has provided a new state, the flag state. We have noticed that it is very sensitive to the boundary conditions (as opposed to the vortex case), much like topological phases in the quantum case. Besides the usual states in the in-plane case for the hexagonal lattices, we have computed the exact minimum energies per dipole for the three-dimensional cases as a function of the interlayer distance. Remarkably, in the simple hexagonal case the layers are formed by flag states, whereas in the hexagonal closed packed configuration, the layers are the vortices themselves. Also, an energy analysis shows that the hexagonal closed packed system is more bound.
Summing up, we have presented that, within the limits of the classic dipole-dipole interaction, one can observe non-trivial topological structures that have special interest, in particular the minimum energy vortex configuration. Thus, we showed that particular configurations of dipoles are not exclusive to quantum spin textures such as skyrmions or bimerons, but also occur in the special case of triangular or hexagonal lattices.

Data availability
The data are available upon reasonable request at jbv276@uib.es.