Cationic vacancies as defects in honeycomb lattices with modular symmetries

Layered materials tend to exhibit intriguing crystalline symmetries and topological characteristics based on their two dimensional (2D) geometries and defects. We consider the diffusion dynamics of positively charged ions (cations) localized in honeycomb lattices within layered materials when an external electric field, non-trivial topologies, curvatures and cationic vacancies are present. The unit (primitive) cell of the honeycomb lattice is characterized by two generators, J1,J2∈SL2(Z)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_1, J_2 \in \mathrm{SL}_2({\mathbb {Z}})$$\end{document} of modular symmetries in the special linear group with integer entries, corresponding to discrete re-scaling and rotations respectively. Moreover, applying a 2D conformal metric in an idealized model, we can consistently treat cationic vacancies as topological defects in an emergent manifold. The framework can be utilized to elucidate the molecular dynamics of the cations in exemplar honeycomb layered frameworks and the role of quantum geometry and topological defects not only in the diffusion process such as prediction of conductance peaks during cationic (de-)intercalation process, but also pseudo-spin and pseudo-magnetic field degrees of freedom on the cationic honeycomb lattice responsible for bilayers.

entailing M atoms coordinated with oxygen around D atoms, arranged in a honeycomb fashion (as shown in Fig. 1a). Thus, of specific interest to our work is the diffusion dynamics of the de-localized cations when an external electric field and non-trivial topology, curvatures and specific defects are present. 41,42 Experimental investigations reveal the diffusion of cations to be largely restricted along honeycomb pathways in honeycomb layered tellurates such as Na 2 Ni 2 TeO 6 (as shown in Fig. 1b). 15 While computational studies are consistent with this observation, they further suggest diffusion of cations in honeycomb pathways is restricted in honeycomb layered oxides for other materials as well, such as NaKNi 2 TeO 6 and A 2 Ni 2 TeO 6 where A = Na, K, Rb and Cs are cations with a large ionic radius, exhibiting a prismatic coordination with oxygen atoms of the Ni and Te octahedra forming the inter-layers. 16,43,44 In particular, Van der Waals forces and Coulomb repulsive forces tend to localize the cations in honeycomb lattices, creating a loosely-bound 2D non-Bravais hexagonal lattice with a two-cation basis known as the honeycomb lattice, which favors de-localization leaving vacancies in the hexagonal vertices only when sufficient activation energy from thermal fluctuations or the electric field, � E = (E x , E y , 0) is present. 45,46 Consequently, provided the ground state of the system devoid of activation energies was initially vacancy-free, the number of vacancies, h in the honeycomb lattice is expected to closely correlate with the number of de-localized (mobile) cations, ν ≃ h ∈ N . Whence, the number of vacancies directly impacts the performance of the material as an effective cathode. 47 Meanwhile, in thin layers of superfluids, superconductors and liquid crystals deposited on curved 2D surfaces, topological defects are known to couple to 2D curvature degrees of freedom, leading to the identification of the number of topological defects as the Euler characteristic of the surface. [48][49][50][51][52][53][54] This lends credence to analogous treatments for cationic vacancies in layered materials. 2,41 Moreover, layered materials demonstrating a bilayer arrangement of metal atoms (with each layer arranged in a triangular lattice) have been found, a vast majority being Ag-based layered oxides and halides such as Ag 2 MO 2 ( M = Co, Cr, Ni, Cu, Fe, Mn, Rh ), Ag 2 F , Ag 6 O 2 (or equivalently as Ag 3 O ), Ag 3 Ni 2 O 4 , and more recently Ag 2 M 2 TeO 6 (where M = Ni, Mg, Co, Cu, Zn). [55][56][57][58][59][60][61][62][63][64][65][66] Preliminary experimental and computational studies reveal that the bilayers represent a monolayer-bilayer phase transition of the honeycomb lattice, with the bifurcation mechanism not clearly understood. 63,67 Herein, we consider how the honeycomb lattice of cations and emergent geometry constrains the model of cationic diffusion in such honeycomb layered oxides, leading to a rich topological description. 41 Consistently treating the number of vacancies, h as the number of holes and handles (genus) of an emergent 2D manifold with a conformal metric, we conclude that the primitive basis and the corresponding Euler characteristic must obey the modular transformation 68 , invariant under the special linear group with integer entries, J ∈ SL 2 (Z ) with α, β, γ , δ ∈ Z and det(J) = αδ − βγ = 1 , where N = 2k is the total number of cation sites enclosed within the parallelogram with the primitive basis labeled by dx and dy as shown in Fig. 2.
In particular, the thermodynamics of the diffusive system of cations can be described by the partition function, Z of the form, which is invariant under the discrete re-scaling ( k → k + 1 ≃ k , for large k) and discrete rotations ( k → −1/k ) of the primitive cell for small k, generated respectively by, when the Euler characteristic can be written in a Fourier series, is the p-th Betti number 69 associated with the n-th primitive cell and topology labeled by h, k =βM/2 ∈ N is an integer, M is the average potential energy of the cations (Strictly speaking, it can be taken as the inverse of the Green's function of the system dominated by its potential energy, and hence can be complex-valued when the system has a finite life-time. 70 ), β = 1/2πk B T is the 'reduced' inverse temperature with k B Boltzmann's constant and f h is a constant independent of k but otherwise is dependent on the topology, h and temperature, T = 1/k Bβ .
The conformal geometry and resultant Betti numbers can be taken to emerge from the underlying matrix group field theory of the self-interactions and quantum correlations of the cations in each primitive cell. In particular, a matrix large N group field theory can be considered, where N = exp(2πk) is the size of the group. 17 Consequently, this construction effectively treats cationic vacancies in honeycomb layered materials as topological defects in the underlying field theory, related to modular symmetries of the honeycomb lattice in the context of emergent 2D quantum geometries. 19 For instance, we show that Eq. (2a) follows from the partition function of k pairs of cations forming the primitive cells of the honeycomb lattice, whereby each pair interacts via the Ising Hamiltonian due to emergent pseudo-spin and pseudo-magnetic degrees of freedom associated with the modular symmetry and broken conformal symmetry respectively. Analogous pseudo-degrees of freedom have also been considered in graphene-based systems. 13,14 Thus, the framework can be utilized to elucidate the molecular dynamics of the cations in exemplar honeycomb layered frameworks and the role of quantum geometry and topological defects not only in the diffusion process such as prediction of conductance peaks during the cation (de-)intercalation process, but also pseudo-spin and pseudo-magnetic field degrees of freedom on the cationic honeycomb lattice whose interactions predict cationic bilayered frameworks. 67 Hereafter, we shall set Planck's constant, the speed of electromagnetic waves in the material, c , Boltzmann's constant, k B and the elementary charge of the cations, q e to unity, =c = k B = q e = 1 , and employ Einstein summation convention unless explicitly stated otherwise.

The model
To build an intuitive geometric and topological picture of cationic diffusion, we begin by summarizing crucial results and clues from an idealized model, previously considered by the present authors in a separate publication. 41 Whilst the cations are positively charged, charge conservation requires the cationic vacancies created www.nature.com/scientificreports/ through de-localization by the electric field, E present in the 2D plane to be considered electrically neutral. Nonetheless, the vacancies can be treated as possessing a fictitious 'magnetic moment' given by, where n is the unit normal vector to the 2D layer comprised of the honeycomb lattice of cations. Thus, the cations diffusing through and around the vacancies created by the electric field, E within the honeycomb lattice will introduce the ' Aharonov-Casher' phase 71 , where ∂A is the boundary of a 2D patch, A of the manifold associated with the 2D layer spanning the primitive cells, where locally the Cartesian coordinates are given by � x = (x, y, z) . Aligning the manifold with the x − y plane, the unit normal vector becomes � n = (0, 0, 1). Applying Stokes' theorem and substituting respectively Maxwell's equations (Gauss' law) and the normalization, where ρ(x, y) = ρ 2D (x, y)/β , ρ 2D (x, y) is the 2D cation number density, β is set as the integration cut-off scale in the z-direction and ν is the total number of mobile cations, yields the condition, where ν ≃ h ∈ N ≥ 0 is also the number of vacancies.
Moreover, assuming the diffusion paths trace arc lengths defined by the 2D conformal metric, where g ab is the 2D metric tensor, �(x, y) must satisfy Liouville's equation 41 , with K(x, y) the Gaussian curvature of the manifold. In fact, one can relate the two scalar functions, and AC by requiring that, Whence, by Stokes' theorem, the Euler characteristic of the manifold is given by the Gauss-Bonnet/Poincaré-Hopf theorem, where we have used Eqs. (5), (6) and (7) to arrive at our result. For instance, for a compact orientable 2D manifold homeomorphic to h number of simply-connected 2-tori, the Euler characteristic is given by, χ(h) = 2 − 2h where h is the genus of the surface given by, which satisfies ν ≃ h for a large number of diffusing cations, ν → ∞ . Thus, this avails the avenue to treat the number of cationic vacancies as the genus, h which uniquely defines the emergent topology of the manifold. Moreover, using F 0i = � E = (E x , E y , 0) and 1 2 ε ijk F jk = � B = 0 with ε ijk the 3D Levi-Civita symbol normalized as ε 123 = 1 , Eqs. (3) and (7) follow from the phase equations of motion 72,73 , on the Minkowski metric, are Maxwell's equations, ξ µ = (1, � 0) and n µ = ( � 0, 1) are time-like and space-like unit normal four-vectors respectively, F µν = ∂ µ A ν − ∂ ν A µ is the electromagnetic field strength, A µ is the electromagnetic (U(1)) gauge field, * F µν = 1 2 ε µνσρ F σρ is the dual field strength with ε µνσρ the 4D Levi-Civita symbol normalized as ε 1234 = 1 and J µ = (ρ, � J) is the current density of the cations. Thus, by Eq. (10), the 2D charge density is related to the Gaussian curvature by 41 , To incorporate diffusion in the formalism, we introduce the diffusion current given by, corresponding to Fick's first law with D the diffusion coefficient and p the center of mass momentum of the cations. Taking the cation number density to satisfy Boltzmann distribution at equilibrium, with 1 2 M� AC (x, y) a 'gravitational' potential energy governing the diffusion dynamics and M = 2/D a peculiarly defined center of mass effective mass using the diffusion coefficient, and applying Eq. (10) yields, Thus, Eq. (14) correspond to the Hamilton-Jacobi equations for the cations with AC corresponding to Hamilton's principal function, the second equation to the 2D Langevin equation 74,75 and β to the mean-free time/ path between collisions (friction term). Thus, this serves as the motivation for β appearing as the cut-off time and length scale in Eqs. (4) and (10). Moreover, the peculiar relation, M = 2/D can be better understood by applying the Virial theorem 76 , where the averages are evaluated at equilibrium using, In this study, we shall consider the particular Hamiltonian for the cations, with momenta, p j , displacement vectors, r j , m = 1/β a mass per cation parameter defined as the inverse of the mean time/path between collisions, β and V (r j ) ≃ 1 2m µ −1 N k=1 � r k · � r j the leading interaction term in the potential energy definedproportional to µ , the mobility of the cations. Typically, other terms such as the Vashishta-Rahman potential 77 , which capture interactions of the cations with the slabs atoms especially oxygen, contribute higher order terms neglected herein. This requires that the diffusion coefficient, including cation-cation correlation terms 78  Observe that, when cation-cation correlation ( j = k ) terms vanish, the diffusion coefficient becomes the self-diffusion coefficient, whereas √ µ takes the role of frequency of the harmonic oscillator. Moreover, since the mobility is a constant, we can re-define it as 16πG ≡ µ , where G ∼ a 2 and a is taken to be the lattice constant with dimensions of length. We thus have, β = 8πGM and N = 2k =βM = 4GM 2 , where G is a gravitational constant, in obvious comparison with Schwarzschild black hole thermodynamics. 18

Results
Conductance spikes. We note that, Eqs. (10) and (13) require that the diffusion current takes the Chern-Simons form 79 , where σ = 2πβDρ = µρ is the conductivity of a single primitive cell, µ = 2πβD is the mobility (Einstein-Smoluchowski equation) and, is the Chern-Simons level. However, the current (density) we are interested in is not necessarily the Hall current, J AC but the spatial part of J µ , which couples to the electromagnetic tensor, F µν in Maxwell's equations given in Eq. (11). Consequently, this predicts conductance spikes whenever k → k + 1 primitive cells containing N = 2k → N + 2 = 2k + 2 cation sites are activated in the (de-)intercalation process. However, due to typically low measurement sensitivity in existing experimental data, the integer nature of the conductance cannot be ascertained. Nonetheless, distinct conductance spikes can be observed at low resolution (i.e k → k + w with w ∈ N ≫ 1 ) at specific voltage values as shown Fig. 3. The low resolution is an artifact of the multi-layered nature of the materials, with each honeycomb lattice not only contributing active primitive cells with cationic sites during the extraction but also multiple sites getting activated at once at spread out external voltage values. Moreover, experimental data suggests that a large activation energy, E a ≫ E K a ≃ 121 meV, where E K a is the activation energy of potassium, K cations in the honeycomb lattice 46 and the presence of cationic vacancies before the extraction process would tend to disfavor the cation extraction process from occurring at evenly spread-out (low to high) voltage values during cycling, often leading to a solitary broad current peak centered at the high voltage regime, for instance, as can be seen in the I-V cycling characteristics for A 2 Ni 2 TeO 6 with A = Li, Na, K , since the activation energy for Li and Na is vastly greater than that of K, i.e. E Li a > E Na a > E K a . 2 Consequently, we have plotted only the Current, I -Voltage, V characteristics of the extraction process for K 2 Ni 2 TeO 6 in Fig. 3, which exhibits several distinct current peaks across varied voltage values.
Moreover, we can employ Eq. (8), which suggests the partition function, is given by the sum over different geometries of the manifold with distinct topology, www.nature.com/scientificreports/ where R = R ab g ab is the 2D Ricci scalar, R ab = R acbd g cd is the 2D Ricci tensor and R acbd = K(g ab g cd − g ad g bc ) is the general form of the 2D Riemann tensor in Riemannian geometry with the equivalence, assumed to be valid. Evidently, Eq. (19) is the partition function of 2D quantum geometry in Euclidean signature 19 , where the coupling constant corresponds to κ = 1/k . It is worth noting that, considering emergent geometries within crystals to describe defects is not entirely a novel idea, since it has been considered in great detail for disclinations and dislocations within the context of classical geometries with torsion. 20 (20)) obeys the necessary modular symmetries defined in Eq. (2b), which are imposed by the primitive cell of the cations in honeycomb layered materials. Nonetheless, Eq. (2) approaches Eqs. (19) and (20) in the limit κ → 0 , as required.

Modular symmetries.
The honeycomb lattice is spanned by the primitive basis, dx and dy defining a parallelogram enclosing k = 1 pair of cation sites. Since the unit cell in Fig. 2 is a rhombus, we have dx/dy = 1 . Moreover, transforming the basis by the matrix, where, we find that J 1 and J 2 , given in Eq. (2b), correspond to the re-scaling, dx ′ /dy ′ = dx/dy + 1 = 2 and the discrete rotation, dx/dy → dx ′ /dy ′ = −dy/dx = −1 , as illustrated in Figs. 4 and 5 respectively. Moreover, we shall consider the modular form 68 defined as, g(dx, dy) = f (k)dk to completely characterize the honeycomb lattice in an invariant manner, under J ∈ SL 2 (Z) with dx/dy = k . By definition, g(dx ′ , dy ′ ) = g(dx, dy) is invariant under the modular transformations, i.e. k → J · k = (αk + β)/(γ k + δ) . Consequently, f(k) transforms as a modular form of weight 2, Proceeding, we must take the large limit, k → ∞ , which spans the entire honeycomb lattice. Moreover, assuming df (k)/dk = 0 , we obtain, g(dx, dy) = kf (k). Now, consider the diffusion dynamics of the cations given in Eqs. (6) and (14). Defining the velocity, � u =β� p = exp(�)d� x/ds , where ℓ is the arc length interval along the proper length, ds, we obtain, where we have used β M = N = 2k from Eqs. (17b) and (8). Thus, setting f (k) = −2πχ(k) defines the modular form as the action of a particle of mass, M/2 in 2D Riemannian geometry, g(dx, dy) = (M/2) ds or equivalently the exponent of Eq. where g ab = exp(−2ωφ)g ab is the 2D metric tensor with g ab given in Eq. (6), K is the Gaussian curvature associated with g ab and R is the Ricci scalar associated with g ab , Q(ω) is a parameter dependent on k and genus h. Setting � = ωφ and x a = ωX a , the metric in Eq. (24) reduces to the 2D identity matrix, g ab = exp(−2�)g ab = δ ab requiring that the Ricci scalar vanishes, R = 0 getting rid of the last term even when Q(ω) = 0.
The Liouville action reduces to, is the divergent vacuum energy of with g ab and approximated as separate non-interacting fields, C is a constant that will be set to vanish by regularization, A = 1 2π d 2 x is the area element, A (2π) 2 d 2 p → p , p are the allowed momenta/energies of the bosonic field, , ω j = (ω,ω) , ω = ik and ω = −ik . Thus, Riemann zeta function regularization requires 87 , To elucidate the nature of Eq. (26), we consider the Virasoro algebra 88 ,  Proceeding, it is known that the field V (α) = exp(2αφ) is primary when the conformal dimension is given by 89,90 , while the marginal condition for the primary field that guarantees conformal invariance of the theory is L =L = 1 with α = ω = −ω which yields, with ω(k) = ik,ω(k) = −ik and Q (ω) =ω + 1/ω . The central charge is given by, where T = 1/β is the temperature, E = M/2 is a defined energy in order for k = E/T =βM/2 as required. The conformal dimension and spin are given by, respectively. However, the theory is known to be unitary only for c = 1 and c = ∞ corresponding to k = 1 and k → i∞ respectively. Since k is considered real, we must have k = 1 . Thus, k → ∞ not only breaks discrete rotation symmetry but also breaks unitarity. Indeed, unlike Eq. (30), Q ∼ ik is not invariant under k → −1/k , which breaks the discrete rotation symmetry generated by J 2 , where k = N/2 is the number of primitive cells in the honeycomb lattice. Now, we consider the partition function, where q(ω j ) = exp(i2πω j ) , q(ω j ) = exp(−i2πω j ) and ω j = (ω,ω) , ω j = (ω, ω) , arriving at the energy, with k =βM/2 from Eq. (17b). Consequently, where we have used Eq. (32) and defined the momentum of the primary field, V (α) = exp(2αφ) as, Recall that the marginal condition is guaranteed by L =L = 1 , which now translates to, for which χ(h) = 2 − 2h yields, c =c = 1 + 24h or equivalently Q 2 (k) = −(k − 1/k) 2 = 4h.
In this case, unitarity is only achieved for h = 0 , which corresponds to the two-sphere satisfying k = 1 . Moreover, to reconcile with Eq. (26) the extra factor of −1/12 must correspond to the vacuum energy, θ = θ(s = 0) after regularization. This requires the momenta p ∈ N be positive integers, which yields the expression θ = 1 + 2 + 3 + · · · = −1/12 by regularization. Note that, for k ≥ 1 , all other positive values of h cannot satisfy the marginal condition and hence must break conformal invariance.

Discussion
If the Euler characteristic of the n-th primitive cell can be associated with a manifold defined by the Poincaré polynomial, where b (p=0,1,2) n is the p-th Betti number, the Euler characteristic of emergent manifold can be calculated as the Euler characteristic of the connected sum of emergent manifolds corresponding to the primitive cells, n with topology, h, which can be decomposed into a Fourier series given by Eq. (2c), where q(k) = exp(2πik) = 1 with k ∈ N and we have used the fact that only the p = 1-st Betti numbers are additive in a connected sum. Consequently, the simplest real-valued partition function that respects the modular symmetries, J 1 and J 2 in their respective limits of k and is proportional to the partition function in Eq. (19) (and hence Eq. (20)) corresponds to Eq. (2), where the limit h, k → ∞ breaks the discrete rotation symmetry, J 2 of the partition function, Z whilst promoting scale invariance, J 1 . Mathematically, this is required since there exists no holomorphic modular forms of weight 2, invariant under both J 1 , J 2 ∈ SL 2 (Z) transformations. 68 Nonetheless, this can be remedied by considering the Euler characteristic proportional to the almost holomorphic modular form of weight 2 in the large genus limit ( h → ∞) 68 , where for such h primitive cells in a connected sum, we set, Thus, unlike in Liouville CFT, modular invariance and hence conformal invariance is guaranteed for all integer values of h and k. However, it is broken for the primitive cell when χ n = 0 , corresponding to a phase transition. 67 To discuss the effects of such a phase transition, we note that, J 2 2 = −I 2 , where, which implies the unit basis acquires a minus sign under J 2 2 . Since J 2 exchanges the basis dx with dy and viceversa, it corresponds to a discrete rotation when acting on the primitive cell. There are 4 such discrete rotations such that J 4n 2 correspond to complete 2πn rotations of the primitive cell, where n ∈ N is a real number. In addition, J 2(2n+2) 2 exchanges one cationic site in the primitive cell with the other. Thus, under the exchange of two cations belonging to the same primitive cell, the transformation picks up a minus sign. This can be understood as the origin of the pseudo-degree of freedom we shall refer to as pseudo-spin, which distinguishes the two sublattices of the honeycomb lattice. 13,67 Under specific conditions, the pseudo-spin of the graphene lattice can be linked to the spin of the electrons localized on the carbon atoms in the sub-lattice. 13 However, for cations in honeycomb layered oxides, no such identification can be affirmed. Nonetheless, the pseudo-spin degree of freedom, coupled with the SL 2 (Z) group imply the partition function given in Eq. (2a) that the underlying theory of cations is a conformal field theory whose ground state must avoid pseudo-spin frustration by pseudo-spin anti-ferromagnetic behavior, but nonetheless prevents the cations from forming a stable honeycomb lattice due to a repulsive exchange interaction which can be offset by pairing of opposite pseudo-spin degrees of freedom. 67 The effective theory for two pseudo-spin cations ( j = 1, 2 ) in k = βM/2 = N/2 non-interacting honeycomb primitive cells corresponds to the 1D Ising Hamiltonian 91 , where, B(h) = 2πMχ(h)/2 is the pseudo-magnetic field 14,67 in the z-direction interacting with the pseudo-spins, σ z which is taken to be proportional to the Euler characteristic, χ(h) , while A(h) is the Heisenberg term representing the exchange interaction, assumed to depend on the genus, h = ν + 1 with ν the cationic vacancy number.
This Ising model is exactly solvable, where standard calculation for the partition function yields 92 , where Tr h,s is the trace over the genus h and spins σ , and ± are the eigenvalues of the transfer matrix, given by, with, Thus, the non-interacting system of k primitive cells occupied by pairs of N = 2k cations interacting via their pseudo-spins and the pseudo-magnetic field yields the partition function in Eq. (2a) where f h = 2 exp(βA(h)) , which takes on varied values for different topology configurations, h. Moreover, the exponents of components of (38) E 2 (h, τ ) = −24 h n=0 σ 1 (n)q n (τ ), .
(40)  www.nature.com/scientificreports/ the transfer matrix correspond to the pseudo-spin energy states, where E ↑ − E ↓ = 2B corresponds to a gapped phase where the honeycomb lattice bifurcates into bilayers with energies E ↑ and E ↑ due to a finite pseudomagnetic field, B = 0 , whereas E ↑↓ = E ↓↑ = A correspond to the ferromagnetic ( A > 0 ) and anti-ferromagnetic ( A < 0 ) alignment of the pseudo-spins.
In conclusion, we have constructed a consistent framework to treat cationic vacancies in honeycomb layered materials as topological defects, h, by relating the Euler characteristic of the manifold to modular symmetries and 2D quantum geometries. 19 The framework predicts integer conductance spikes during (de-)intercalation process, proportional to the number of active cation sites, k → k + w participating in the diffusion process at high resolution, w ∼ 1 ∈ N , which remain unobserved. Nonetheless, the framework greatly elucidates the geometric nature of the diffusion process which occur in these novel materials, and the crucial role played by cationic vacancies as topological defects, and hence should find great utility in finding avenues for performance optimization of such cathode materials for energy storage. 42