Zhang–Zhang Polynomials of Multiple Zigzag Chains Revisited: A Connection with the John–Sachs Theorem

Multiple zigzag chains Zm,n of length n and width m constitute an important class of regular graphene flakes of rectangular shape. The physical and chemical properties of these basic pericondensed benzenoids can be related to their various topological invariants, conveniently encoded as the coefficients of a combinatorial polynomial, usually referred to as the ZZ polynomial of multiple zigzag chains Zm,n. The current study reports a novel method for determination of these ZZ polynomials based on a hypothesized extension to John–Sachs theorem, used previously to enumerate Kekulé structures of various benzenoid hydrocarbons. We show that the ZZ polynomial of the Zm,n multiple zigzag chain can be conveniently expressed as a determinant of a Toeplitz (or almost Toeplitz) matrix of size m2×m2 consisting of simple hypergeometric polynomials. The presented analysis can be extended to generalized multiple zigzag chains Zkm,n, i.e., derivatives of Zm,n with a single attached polyacene chain of length k. All presented formulas are accompanied by formal proofs. The developed theoretical machinery is applied for predicting aromaticity distribution patterns in large and infinite multiple zigzag chains Zm,n and for computing the distribution of spin densities in biradical states of finite multiple zigzag chains Zm,n.


Introduction
The construction, enumeration, and characterization of Kekulé structures and Clar covers of benzenoid hydrocarbons constitute one of the most important topics of graphtheoretical characterization of chemical molecules, a scientific discipline dwelling in the intersection of theoretical chemistry and discrete mathematics [1][2][3][4][5]. For a chemist, a Kekulé structure is a resonance structure in which each carbon atoms participates in exactly one double bond [6]. Similarly, a Clar cover is a resonance structure in which each carbon atoms participates in exactly one double bond or in exactly one aromatic sextet with some other five carbon atoms [7]. These concepts are explained graphically in Figure 1 using pyrene (i.e., the multiple zigzag chain structure Z(2, 2)) as an example, together with an intuitive definition of the ZZ polynomial of Z(2, 2). For a mathematician, a Kekulé structure K is a spanning subgraph of the molecular graph G corresponding to a given benzenoid B, such that every component of K is isomorphic to a complete graph on two vertices, K 2 . Similarly, a Clar cover C is a spanning subgraph of G, such that every component of C is isomorphic to K 2 or isomorphic to a cyclic graph on six vertices, C 6 [8]. The most important questions pertinent to graph-theoretical characterization of a given benzenoid B are as follows. (i) What is the number of Kekulé structures, K ≡ K(B), that can be constructed for B? (ii) What is the number of Clar covers, C ≡ C(B), that can be constructed for B? (iii) What is the maximal number, Cl ≡ Cl(B), of aromatic sextets C 6 that can be accommodated in B? (iv) How many distinct Clar covers exist, c Cl ≡ c Cl (B), with exactly Cl aromatic sextets? (v) Assuming that k is a non-negative integer, how many distinct Clar covers exist, c k ≡ c k (B), with exactly k aromatic sextets? Note that the answer to the last question, i.e., the sequence of numbers [c 0 , c 1 , c 2 , . . .], is the most general here, being capable of also providing a solution to the preceding questions. Indeed, we have K = c 0 , C = ∑ k c k , Cl = max{k|c k = 0}, and c Cl = c max{k|c k =0} . Note that for a finite benzenoid B with N atoms, the maximal number of aromatic sextets that can be accommodated in B is naturally bounded from above by N/6, and thus for all k > N/6, we have c k = 0. The most convenient way of representing the subsequence [c 0 , c 1 , c 2 , . . . , c Cl ] is given in the form of its generating function which is most often referred to as the Clar covering polynomial or, from the names of its inventors, as the Zhang-Zhang polynomial or the ZZ polynomial of B [9][10][11][12][13][14][15]. Substantial research effort has been invested in the determination of ZZ(B, x) for elementary families of benzenoids [8,. The rapid development of Clar theory stimulated by these discoveries in recent years has led to many new interesting applications and connections to other branches of chemistry, graph theory, and combinatorics [8,[17][18][19]21,28,34,[42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59][60][61]. In one of our recent contributions to the field of characterization of Clar covers of basic pericondensed benzenoids, we have derived [33] a closed-form generating function for the ZZ polynomials of multiple zigzag chains Z(m, n) and generalized multiple zigzag chains Z k (m, n) [2,62]. For a graphical definition of both families of benzenoids, see Figure 2. For Z(m, n), the derived formula (Equation (20) of the work in [33]) had an aesthetically pleasing and compact structure of a finite, regular continued fraction ∞ ∑ m=0 ZZ(Z(m, n), x)t m = 0; −t, (−1) 2 z t, (−1) 3 z t, . . . , (−1) n z t, 1 + (−1) n+1 z t (2) where, as usual in the theory of ZZ polynomials, z = 1 + x, and where the continued fraction notation [a 0 ; a 1 , a 2 , a 3 , . . . , a n ] should be interpreted as in the following example [a 0 ; a 1 , a 2 , a 3 , a 4 ] = a 0 + 1 multiple zigzag chains Z(m, n) (≡ Z n (m − 1, n))  The formula in Equation (2) is quite regular, but probably its high symmetry is most elegantly reflected in a slightly different continued fraction representation: Note that the length of the continued fraction in Equation (4) is n + 1, i.e., the number of times the product z t appears in it is n. Similar expressions, in terms of products of generating functions of the form (4), were obtained for the generalized multiple zigzag chains Z k (m, n) shown in Figure 2.
Despite the internal beauty of the ZZ(Z(m, n), x) and ZZ(Z k (m, n), x) generating functions, it turned out that extracting the original ZZ polynomials from the generating functions in Equations (2) and (4) is a rather formidable task. A binomial expansion of the continued fraction of type (4), starting from the most shallow level, together with an appropriate multiple sum rotation produces a rather lengthy and cumbersome expression for the ZZ polynomial of Z(m, n) given by with the coefficient c j expressed by Similar expansions for ZZ(Z k (m, n), x) were not rigorously attempted, in conviction that their mathematical structure is much too complicated and of little practical importance. Equations (5) and (6) provide formally a closed-form formula for ZZ(Z(m, n), x), but at the same time, they do not bring full intellectual satisfaction of possessing such a formula, as Equation (5) provides little, if any, conceptual insight into the problem of enumeration of Clar covers of Z(m, n). We concluded our previous work [33] by expressing a hope that the evaluation of multiple sums in Equation (5) could lead to a more compact and transparent expression, but to date no such a result has been reported.
In the current work, we report the discovery and demonstrate formal validity of new ZZ polynomial formulas for the multiple zigzag chains Z(m, n) and the generalized multiple zigzag chains Z k (m, n). The new formulas possess both a transparent algebraic form and deep internal structure, which reveal a direct connection to the structural parameters m, n, and k, and offer a possibility for fast and robust determination of the ZZ polynomials in a form of a determinant of (n + 2)-diagonal Toeplitz (or almost Toeplitz) matrices of size m 2 × m 2 or m+1 The new formulas are numerically consistent with the previously obtained Equations (5) and (6). In a sense, the reported formulas are a direct conceptual extension of the formulas in Equations (14.34)-(14.36) and Equations (14.38)-(14.39) of the seminal work of Cyvin and Gutman [2], which reported the determinantal John-Sachs expressions for the number of Kekulé structures for Z(m, n) and Z k (m, n). The actual form of the reported Toeplitz matrices reported here originates from an anticipated extension of the John-Sachs theorem to the world of Clar covers. This hypothesized extension encouraged us to recently conduct an intensive search of examples of elementary benzenoids B, for which it is possible to construct a generalization of the John-Sachs matrix, whose determinant would be equal to the ZZ polynomial of B and which, upon the evaluation x = 0, would reduce to the original John-Sachs matrix of B.
Among other examples given in [63], we have been able to achieve this task for multiple zigzag chains with n = 1, 2, and 3. Here, we have been able to achieve this goal for general values of the structural parameters m and n for Z(m, n), as well as for an arbitrary-sized structures Z k (m, n). In addition, the somewhat ad hoc derivations of these formulas are accompanied here with a formal mathematical proof of their validity using the usual recurrence relations obeyed simultaneously by the ZZ polynomials of Z(m, n) and Z k (m, n) and the determinants of the discovered Toeplitz matrices.
In the last section of the current manuscript, we present an application of the developed theoretical framework for computing the aromaticity distribution patterns and the biradical spin population patterns in finite and infinite multiple zigzag chains Z(m, n). These structures, usually referred by physicists as graphene nanoribbons, received considerable attention from the physics community as promising candidates for future nanoelectronics and spintronics [64][65][66][67][68].

Preliminaries
To explain our reasoning leading to the discovery of determinantal formulas for the ZZ polynomials of the multiple zigzag chains Z(m, n) and generalized multiple zigzag chains Z k (m, n), we need to start with some formal preliminaries. Assume, like always, that a benzenoid B is a planar graph embedded in a hexagonal lattice. We define a Kekulé structure K of B as a spanning subgraph of B whose components are K 2 (i.e., complete graphs on 2 vertices). Similarly, we define a Clar cover C of B as a spanning subgraph of B whose components are either K 2 or C 6 (i.e., cycle graphs of length 6) [7]. The number of the C 6 components in a given Clar cover C will be referred to as the order of this Clar cover. Denoting by c k the number of distinct Clar covers of order k, the generating function for the sequence [c 0 , c 1 , c 2 , . . . , c Cl ] is usually referred to as the Clar covering polynomial (or Zhang-Zhang polynomial or ZZ polynomial) of B and is formally defined as where the maximal number Cl of the C 6 components that can be accommodated within B is referred to as the Clar number of B [9][10][11][12][13][14][15]. In practice, the ZZ polynomial of an arbitrary benzenoid B can be robustly computed using recursive decomposition algorithms [4,20,26] or can be conveniently determined using interface theory of benzenoids [31,39,[69][70][71][72]. A useful theoretical tool for finding ZZ polynomials for an arbitrary benzenoid is ZZDecomposer [25,26]. With this program, one can conveniently define a benzenoid using a mouse drawing pad and subsequently use the underlying graph representation of the benzenoid to find its ZZ polynomial, manipulate its Clar covers, and determine its structural similarity to other related benzenoids. ZZDecomposer has been successfully applied to find close-form formulas of ZZ polynomials for numerous families of basic elementary benzenoids [20,[22][23][24][25]27,29,30,[32][33][34]36,37,63,73,74]. Assume further that the benzenoids Z(m, n) and Z k (m, n) are drawn in a plane in the way shown in Figure 3 with some of their edges oriented vertically. We say that a vertex p i is a peak, if all its neighbors are located below p i . Similarly, we say that a vertex v j is a valley, if all its neighbors are located above v j . It is a well-known fact that if the structure B is Kekuléan, and both families of B = Z(m, n) and B = Z k (m, n) possess this quality for any admissible values of the structural parameters m, n, and k, then the number of peaks in B is equal to the number of valleys in B; let us denote this number as l. It is clear from Figure 3 that l = m 2 for Z(m, n) with even m, l = m+1 2 for Z(m, n) and Z k (m, n) with odd m, and l = m+2 2 for Z k (m, n) with even m. The next two important concepts are the wetting region of a peak p i , which is defined as a subgraph of B spanned by all the vertices of B that are accessible from p i by going exclusively downward, and the funnel region of a valley v j , which is defined as a subgraph of B spanned by all the vertices of B that are accessible from v j by going exclusively upward. Finally, the p i → v j path region is defined as the intersection of the wetting region of p i and the funnel region of v j . All these concepts [75] are illustrated schematically in Figure 4.   Let us finally denote by K p i → v j the number of Kekulé structures for the p i → v j path region. Then, for an arbitrary Kekuléan benzenoid B with l peaks p 1 , . . . , p l and l valleys v 1 , . . . , v l , the John-Sachs matrix P(B) is defined as a square, l × l matrix with elements P ij = K p i → v j . The celebrated John-Sachs theorem (see, for example, Theorem 1 in [76] or Equation (5) in [77]) states that the number of Kekulé structures for B is equal to the determinant of P(B), In our recent work [63], we have asked whether the concept of the John-Sachs matrix P(B) can be generalized to encompass the theory of Clar covers. The natural extension meant to make the matrix P(B) x-dependent, i.e., replacing P ij = K p i → v j by P ij = ZZ p i → v j , x , does not work, as the determinant of the resulting matrix does not evaluate to the ZZ polynomial of B. Notwithstanding, we have discovered that for numer-ous families of basic benzenoids, it is actually possible to find a generalized John-Sachs matrix P ZZ (B) that possesses all the qualities expected for such a generalization: 1.
The determinant of the John-Sachs matrix P(B) is equal to the number of Kekulé structures as stipulated by Equation (8). Similarly, we request that the determinant of the generalized John-Sachs matrix P ZZ (B) is equal to the ZZ polynomial of B 2.
As the Kekulé structures of B are simply the Clar covers of B of order 0, the number K(B) of Kekulé structures can be obtained from the ZZ polynomial of B by evaluating it at which effectively corresponds to removing from the set of Clar covers those which contain at least one component C 6 . It is then only natural to request that in the same limiting process the John-Sachs path matrix P(B) is obtained from the generalized John-Sachs path matrix P ZZ (B) upon the evaluation at Numerous examples of generalized John-Sachs path matrices P ZZ (B) have been given in our recent work [63]. Here, we augment this collection with two further examples corresponding to multiple zigzag chains Z(m, n) and generalized multiple zigzag chains Z k (m, n). Note that most of the previous results were given without a proof, while here for both Z(m, n) and Z k (m, n) we are able to furnish formal proof of all of the presented results. Interestingly, the off-diagonal elements of the generalized John-Sachs path matrices P ZZ (B) have a rather unexpected form, for which no plausible explanation nor interpretation has been yet discovered. We hope that a formal extension of the John-Sachs theory to the world of Clar covers will explain this conundrum. We also hope that the results given in the current work will stimulate the community to work on such the generalization of John-Sachs theorem to the world of Clar covers. We finally truly hope that the designed here strategy will be sufficient to discover one of the last missing gems in the theory of ZZ polynomials of basic benzenoids, the ZZ polynomial formula for an arbitrary hexagonal flake O(k, m, n), similarly like the original John-Sachs theory was used to prove analogous formula for K(O(k, m, n)).

Discovery of the Determinantal Formulas
In this section, we briefly discuss the process that led us to the discovery of the determinantal formulas of the ZZ polynomials of multiple zigzag chains Z(m, n) and generalized multiple zigzag chains Z k (m, n). Our motivation for presenting this somewhat obsolete reasoning is purely pedagogical, as it is hoped that similar heuristic reasoning processes can lead to a discovery of other missing formulas in the theory of ZZ polynomials. We limit our discussion to the case of Z(m, n) with an even value of m. In our previous work [63], we have discovered that the (multiple) zigzag chains Z(m, 1), Z(m, 2), and Z(m, 3) are characterized by the following generalized John-Sachs path matrices: where, as before, z = 1 + x. A number of regularities can be immediately noticed: • The matrices P ZZ (Z(m, n)) are m 2 × m 2 Toeplitz matrices with n + 2 diagonals: 1 subdiagonal (z), 1 diagonal (w 0 ), and n consecutive superdiagonals (w 1 , . . . , w n ). • The value on the subdiagonal, z = 1 + x, does not depend on n.

•
The value of the diagonal element w 0 is equal to ZZ(M(2, n), x), given explicitly by This value is consistent with the natural extension of the John-Sachs matrix P(B) diagonal elements, i.e., with replacing The value w l on the l-th superdiagonal is a product of a multiplicative factor z l , a numerical factor c nl , and a polynomial p nl (z) of degree 2, 1, or 0. Further analysis requires slightly larger amount of data. Similar techniques to those used in our previous work [63] allow us to find the values of w l for P ZZ (Z(m, 4)) and P ZZ (Z(m, 5)); thus, we have Equipped with this knowledge it is relatively easy to identify the last two missing pieces of the puzzles: the numerical factor c nl and the polynomial p nl (z). 1.
The polynomials p nl (z) always start from 1 and contain coefficients, which factorize into small primes. This suggests that they are hypergeometric polynomials, i.e., hypergeometric functions p F q a 1 ,...,a p b 1 ,...,b q ; z with at least one of the upper indices a j being a negative integer. As the formula for the diagonal elements in Equation (12) contains a hypergeometric function 2 F 1 , it is natural to seek p nl (z) also in this form. The lower parameter b 1 is suggested by the denominator in the linear term of the polynomials p nl (z), and it is equal to 3 for w 1 , 5 for w 2 , 7 for w 3 , etc. Therefore, in a general case, we can expect for w l a value of b 1 = 1 + 2l. Another observation is that the coefficient in the linear term of the polynomial, equal to a 1 a 2 b 1 , is always positive, which, taking into account that at least one of the upper indices should be negative, shows that actually both a 1 and a 2 are negative. One of these numbers, say a 1 , is immediately recognizable as −2, because of the constant degree of the polynomials p nl (z) for l ≤ n − 2. The other index, a 2 , is a function a 2 = a 2 (n, l), which for l = n should be 0 (because of the degree 0 of the polynomial p nl (z) for l = n) and for l = n − 1 should be 1 (because of the degree 0 of the polynomial p nl (z) for l = n). All these facts suggest that a 2 = l − n. A straightforward verification with Maple [78] shows that indeed the hypergeometric function 2 ; z reproduces all the reported polynomials p nl (z) in the matrix elements w l given above. Note that this polynomial reproduces somewhat fortuitously also the diagonal entry w 0 . 2.
The last remaining task is the identification of the two-dimensional sequence of numbers c nl , which numerical values for small n and l are given by the following triangle: This task can be readily performed by typing, for example, the last rows of this triangle in the The On-Line Encyclopedia of Integer Sequences [79], which recognizes it as the sequence A085478 generated by These two identifications conclude our heuristic deduction and allow to express the ZZ polynomial of multiple zigzag chains Z(m, n) with even m as a determinant of the generalized John-Sachs path matrix P ZZ (Z(m, n)) where the matrix elements w l are given explicitly by The hypergeometric function in Equation (17) can be expanded as a power series in (1 + x), which gives a binomial-like definition convenient for evaluation: While the definition of w l given by Equation (18) is probably more transparent and practical, we want to stress that the hypergeometric form in Equation (17) has been essential in the process of the identification of w l , furnishing a convenient unified framework for the search process.
It is further possible to make the form of the matrix P ZZ (Z(m, n)) slightly more uniform by noticing that the presence of the binomial coefficient ( n+l n−l ) in Equation (17) (or the binomial ( n+l 2l+j ) in Equation (18)) imposes that which means we could omit the triangle of zeros in Equation (16) and simply fill all the superdiagonals with the entries w l . We also note in passing that Equation (18) naturally extends also to l = −1, giving which allows to rewrite the entry z = 1 + x on the first subdiagonal of Equation (16) simply as w −1 . This gives the most transparent form of the ZZ polynomial of a multiple zigzag chain Z(m, n) with even m as where the matrix elements w l are computed using Equation (18). The correctness of Equation (16) has been extensively tested by comparing it with the ZZ polynomials of multiple zigzag chains Z(m, n) computed by brute-force recursive calculations using ZZDecomposer for various values of the structural parameters m and n. A formal proof of correctness of Equation (16) is presented in Section 4.
So far, we have discussed the determinantal formula for the ZZ polynomial of multiple zigzag chains Z(m, n) with even m. Now, we explain how analogous formulas for Z k (m, n) and for Z(m, n) with odd values of m have been discovered. It is probably simplest to start with Z k (m, n) with even m. Figure 3 suggests that for even m, the generalized John-Sachs matrix P ZZ (Z k (m, n)) corresponds to a matrix of size m 2 + 1 × m 2 + 1 . As, for even values of m, the first m 2 peaks and m 2 valleys of Z k (m, n) coincide with those of Z(m, n), it is natural to expect that the first, diagonal m 2 × m 2 block of P ZZ (Z k (m, n)) is identical to P ZZ (Z(m, n)). The remaining diagonal entry corresponds to the p m+2 2 → v m+2 2 path and according to the hypothesized generalization of the John-Sachs theorem it can be expressed as ZZ(M(1, k), x) = 1 + k(1 + x). It is also natural to assume the value on the first subdiagonal still is z = 1 + x. Discovering the values on the superdiagonals requires intensive numerical experimentation similar in character to the efforts needed to discover w l . Fortunately, this process is relatively straightforward, as independent variations of two structural parameters-k and m-supply a sufficient amount of numerical data to complete the identification in an almost linear fashion. It turns out that the value v k l located on the l th superdiagonal in the last column of P ZZ (Z k (m, n)) with even m can be expressed by the following hypergeometric formula: or its binomial-like equivalent where s = k. Both formulas bear close resemblance to Equations (17) and (18). It is obvious from the presence of the binomial coefficient ( s+l k−l ) in Equation (22) that v s l = 0 for l > s, which means that the last column contains only k + 1 non-zero entries, including the diagonal element. Somewhat fortuitously, the diagonal element, v k 0 = 1 + kz, can be also expressed by Equations (22) and (23) with l = 0 and s = k.
No effort is needed to discover the form of the generalized John-Sachs matrix P ZZ (Z(m, n)) for Z(m, n) with odd m; it is immediately identified as P ZZ (Z n (m − 1, n)), because of the structural identity Z(m, n) = Z n (m − 1, n).
The last remaining structure, for which we need to find the generalized John-Sachs matrix, is Z k (m, n), with an odd value of m. This is a considerable harder task. Here, P ZZ (Z k (m, n)) corresponds to a matrix of size m+1 (for details, see Equation (10) in [35]). Extensive numerical experimentation allows to establish that the value u k l located on the l th superdiagonal in the last column of P ZZ (Z k (m, n)) with odd m can be expressed by where t = k, with v n l and v n−j l given by Equations (22) or (23). Direct evaluation of Equation (25) gives the following hypergeometric representation of this term: or its binomial equivalent Explicit determinantal formulas for ZZ polynomials of multiple zigzag chains Z(m, n) and generalized multiple zigzag chains Z k (m, n) discovered in the previous paragraphs have been presented in Figure 2 together with binomial representation of all the matrix elements required for their determination and graphical representation of the analyzed structures. Figure 2 can be thought as the graphical abstract of the current work summarizing the most important discoveries reported here. For those who prefer hypergeometric representation of the elements of the relevant path matrices, the following compact formulas should suffice:

Formal Proof of Determinantal Formulas for ZZ(Z(m, n), x) and ZZ(Z k (m, n), x)
We start the proof by demonstrating that the formulas presented in Figure 2 produce the ZZ polynomials of Z(m, n) and Z k (m, n) for low values of the index m. We argue in the next paragraph that ZZ polynomials of multiple zigzag chains Z(m, n) with m ≥ 3 and ZZ polynomials of generalized multiple zigzag chains Z k (m, n) with m ≥ 1 can be computed in a recursive fashion from analogous ZZ polynomials of Z(m, n) and Z k (m, n) with lower values of the index m. Such a recursive algorithm works neither for Z(m, n) with m = 1 or m = 2 nor for Z k (m, n) with m = 0, so we treat these cases separately here. Simple geometric considerations show that Z(2, n) = M(2, n), Z(1, n) = M(1, n), and Z k (0, n) = M(1, k). As the ZZ polynomial of a parallelogram M(m, n) is well known (ZZ(M(m, n), x)=∑ m j=0 ( m j )( n j )(1 + x) j , see for example Equation (4) of [24]), we obtain the following relations: which confirm that the ZZ polynomials of Z(1, n), Z(2, n), and Z k (0, n) can be computed from the corresponding generalized John-Sachs matrices presented in Figure 2. For Z(m, n) with m ≥ 3 and for Z k (m, n) with m ≥ 1, the formal proof is a consequence of the general recursive relations derived by Zhang and Zhang for ZZ polynomials. For the multiple zigzag chains Z(m, n) and generalized multiple zigzag chains Z k (m, n), the relevant recursive relations are derived in Figure 5 by choosing a particular edge in Z(m, n) and Z k (m, n) and assigning it to K 2 , C 6 , or none. The three possibilities are labeled in Figure 5 as D, R, and S, respectively, where the designated letter abbreviates the corresponding chemical terms: a double bond, an aromatic ring, and a single bond. The resulting recurrence relations interconnecting the ZZ polynomials of various structures Z(m, n) and Z k (m, n) should be consistently satisfied for all the permissible values of the structural parameters m, n, and k. Noting that Z(m, n) ≡ Z 0 (m, n) ≡ Z n (m − 1, n), we see that Equation (34) is a specialized version of Equation (35) with k = n and m replaced by m − 1. It is therefore sufficient to consider only Equation (35), which is valid for the following values of the structural parameters: 1 ≤ m and 1 ≤ k ≤ n together with the initial conditions given by Equations (31)-( 33) and corresponding to the edge cases Z(1, n), Z(2, n), and Z k (0, n). Equipped with this input, we are able to demonstrate below that the deduced earlier determinantal formulas for ZZ(Z(m, n), x) and ZZ(Z k (m, n), x) (see Figure 2) satisfy the recurrence relation in Equation (35).

Z(m, n) with Even m
Let us compute both sides of Equation (34) assuming the correctness of the determinantal formulas given in Figure 2. LHS of Equation (34) evaluates to |P ZZ (Z(m, n))| given by Equation (21). RHS of Equation (34) where the value of u n−1 l defined in Equation (27) evaluates to the following explicit expression Substituting the values of u l given by Equation (37) to the first determinant in Equation (36) and splitting it with respect to the last column into a sum of two determinants, (see, for example, in [80]) allows us to rewrite Equation (36) as a sum of |P ZZ (Z(m, n))| and Laplace expansion of the first determinant in Equation (38) with respect to the last column shows that the both addends in Equation (38) cancel each other, retaining only |P ZZ (Z(m, n))| as the evaluation of RHS of Equation (34). This is identical to the LHS of Equation (34), proving the consistency of the determinantal formulas reported here with the recurrence relation in Equation (34) for even m.

Z(m, n) with Odd m
The analogous consistency demonstration for the ZZ(Z(m, n), x) determinantal formulas with odd m is slightly more complicated. We start with LHS of Equation (34) evaluated to the determinantal form with odd m, and perform the following consecutive operations on the following determinant: • Use v s l given by Equation (23) • Identify ( ) the second determinant in Equation (40) (40) by subtracting its last column from its last-but-one column obtaining • Direct calculations using the definitions in Equations (18) and (23) allow us to establish that Using these relations simplifies the determinant in Equation (41) to the following form: • Laplace expansion of the determinant in Equation (42) with respect to the last row gives allowing us to identify ( ) the expression in Equation (43)

Z k (m, n) with Even m
The consistency check process in this case is quite similar to that discussed in the previous section, so we minimize the discussion, giving only its most important intermediate steps. We start with the LHS of Equation (35) evaluated for even m to the determinantal form with the help of Figure 2 obtaining We perform the following consecutive operations on this determinant: • We rewrite the last column of the determinant in Equation (44)  In the first determinant, we subtract the last column from the last-but-one column obtaining • Direct calculations using the definitions in Equations (18), (23) and (27) allow us to establish that using these relations simplifies the determinant in Equation (45) to the following form: • Laplace expansion of the determinant in Equation (46) with respect to the last row gives allowing to identify ( ) this expression as (1 + x)ZZ(Z n−k (m − 1, n), x). The sum of both performed identifications, ( ) and ( ), representing LHS of Equation (35) by construction, is identical to RHS of Equation (35) evaluated to the determinantal form for even m using the formulas presented in Figure 2.

Z k (m, n) with Odd m
The consistency check process in this case relies on the identity u k l = u k−1 l + (1 + x)v n−k l easily derivable directly from the definition of u k l given by Equation (25). We start with LHS of Equation (35) evaluated to the determinantal form for odd m with the help of Figure 2, obtaining Substituting the aforementioned identity for for every u k l in the last column of the determinant in Equation (48) and representing it as a sum of two determinants (decomposition with respect to the last column) gives The first determinant in Equation (49) is identified using Figure 2 as ZZ(Z k−1 (m, n), x), and the second determinant as ZZ(Z n−k (m − 1, n), x). Altogether, the expression in Equation (49)

Chemical Applications
In the current section, we use the derived equations to determine various chemical properties of multiple zigzag chains Z(m, n). In particular, we focus on the aromaticity pattern distribution within a multiple zigzag chain Z(m, n) and on the distribution of unpaired spins in biradical excited states of multiple zigzag chains Z(m, n). First, we define the ZZ aromaticity indicator and the ZZ spin populations and compare them with analogous quantities obtained from rigorous quantum chemical calculations using selected systems. Subsequently, the ZZ indices are computed for families of isostructural multiple zigzag chains and their variation with growing size of the graphene flake is analyzed. The developed ZZ methodology allows one to treat substantially larger flakes than the conventional quantum chemical calculations and in some cases permits computing the asymptotics in transition to infinite graphene nanoribbons.

Local ZZ Aromaticity Indicators for Multiple Zigzag Chains Z(m, n)
The ZZ aromaticity index for some hexagon h i of a multiple zigzag chain Z(m, n) is defined as a ratio of the number of Clar covers of Z(m, n) in which the hexagon h i has aromatic character to the total number of Clar covers of Z(m, n). The denominator in this ratio is easily computed as ZZ(Z(m, n), 1), while the hexagon-dependent numerators must be computed independently for each hexagon of Z(m, n). The maximal hexagon aromaticity computed in this way for an arbitrary graphene flake cannot exceed the theoretical value characterizing the infinite graphene sheet, i.e., 33%; the minimal value has a natural lower bound of 0%. The absolute aromaticity scale obtained in this way is very convenient, as it not only allows for finding most and least aromatic hexagons within a given graphene flake, but it also allows for aromaticity comparisons between flakes with various shapes.
The ZZ aromaticity indices computed in this way are shown schematically in Figure 6 for two multiple zigzag chains: Z(5, 3) and Z(4, 6). The maximal ZZ hexagon aromaticity in Z(5, 3) is approximately 67% of the value expected for the infinite graphene sheet, while for Z (4,6), the maximal value is smaller (45%). A few general observations can be made. (i) Various quantum chemical aromaticity indicators show rather large discrepancies in the predicted aromaticity patterns. In particular, the Bird aromaticity indicator differs substantially from all the other indicators. (ii) The ZZ aromaticity pattern for Z(5, 3) agrees quite well with analogous patterns predicted by HOMA, NICS, PDI, and FLU, i.e., with the majority of selected quantum chemical aromaticity indicators. (iii) The ZZ aromaticity pattern for Z(4, 6) is similar to the patterns predicted by majority of other indicators (NICS, PDI, FLU, SCI, and PLR). (iv) The hexagons with maximal aromaticity in Z(5, 3) and Z (4,6) are correctly predicted by the ZZ approach. (v) The hexagons with minimal aromaticity (i.e., with maximal antiaromaticity) in Z(5, 3) and Z(4, 6) are again correctly predicted by the ZZ approach. In particular, the ZZ approach is the only one to indicate very clear antiaromatic character of the hexagons characterized by positive chemical shift (blue dots in Figure 6) in the NICS approach. The two benchmark comparisons presented in Figure 6 suggest that the ZZ polynomial approach to graphene flakes aromaticity can be very useful, allowing to treat systems much larger than those accessible with quantum chemical calculations, as the computational cost of computing ZZ polynomials is only fractional in comparison to quantum chemical calculations. (left) and Z(4, 6) (right) are compared to analogous quantities obtained from rigorous quantum chemical calculations. We use the following quantum chemical aromaticity indicators: NICS [81], HOMA [82,83], Bird [84], PDI [85], SCI [86], PLR [87][88][89][90], and FLU [91]; for more details, see in [92]. The size of the red circle corresponds to the aromaticity of each hexagon, with the largest circle corresponding to the maximal aromaticity obtained with each method. The blue circles in NICS aromaticity patterns correspond to anti-aromatic hexagons with positive chemical shift.
The ZZ approach is now used to study the convergence of the aromaticity patterns for the family of multiple zigzag chains Z(5, n) with the growing value of n. In Figure 7, we present the computed ZZ aromaticity indicators for all the hexagons of Z(5, n) with n = 3, 6, 9, 12, 15, 21, and 45. A few observations are in place. (i) The aromaticity pattern in Z(5, n) seems to converge with n to a uniformly non-aromatic infinite graphene strip. (ii) The location of the most aromatic and the most anti-aromatic hexagons is not altered by the elongation of the flake. (iii) In each of the five polyacene chains in Z(5, n), the aromaticity decreases monotonically from the convex end to the concave end. By studying the recursive decomposition of Z(5, n), it is possible to derive an analytical formula for the aromaticity of the most aromatic hexagon in each polyacene row. Denoting by A k (n) the aromaticity of the most aromatic hexagon in the row k of Z(5, n), the results of the analysis can be given as follows: Each of these formulas converges like c/n to 0 when n → ∞ with the coefficient c equal to 25/32 for A 1 (n), 40/32 for A 2 (n), and 30/32 for A 3 (n), showing that an infinitely long graphene strip of width 5 is completely non-aromatic, i.e., that all the bonds have either single or double character. At the same time, the formulas show that in each finite nanoribbon Z(5, n) the most aromatic hexagon is located at the convex end of the second polyacene row. Note that obtaining these conclusions would not be possible directly from quantum chemical calculations without the developed here ZZ polynomial methodology.

Spin Densities in Biradical Multiple Zigzag Chains Z(5, n)
Another application of the developed here methodology is the computation of spin densities in biradical multiple zigzag chains Z(5, n). The probability of finding an unpaired electron on the p z orbital of a selected carbon atom is computed from the associated ZZ polynomial in a similar manner as described above for aromaticities. Namely, we consider all possible distributions of two unpaired electrons in a biradical Z(5, n) structure by localizing them at the p z orbitals of two selected carbon atoms-C i and C j , and for each such distribution, we compute the number of Clar covers corresponding to it by evaluating the associated ZZ polynomial at x = 1 (or equivalently, at z = 2). This number, divided by the total number of Clar covers in the non-radical Z(5, n) structure, gives the probability of finding the two unpaired electrons at positions C i and C j simultaneously. Summing such contributions over all the possible distributions of two unpaired electrons gives the total atomic probabilities of finding an unpaired electron over each carbon atom.
The probabilities obtained in this way are then normalized to 1 and the resulting ZZ spin density probability distributions are plotted as red (spin up) and blue (spin down) circles centered at carbon atoms; for more details, see in [92]. The results computed in this manner are plotted in Figure 8 for two families of multiple zigzag chains: Z(5, n) with n = 4, 5, 6, 7, 8, 9, 10, 12, 15 and 21 (upper panel) and Z(m, 10) with m = 4, 5, 6 and 7 (lower panel). All the presented spin density distributions look alike and consist of a single ↑ spin wave at one zigzag edge of Z(m, n) and a single ↓ spin wave at the other zigzag edge of Z(m, n). The spin densities at the carbon atoms in the interior of Z(m, n) are negligible, except for the smallest studied here structures. The spin waves at each zigzag edge have no nodal points with the maximum of the distribution located approximately at the middle of the edge. By extrapolating this picture to infinite graphene strips, one can expect that the zigzag edges of the strip will have a constant spin density for each carbon atom located at each zigzag edge, with one side corresponding to spin up, and the other side, to spin down. The resulting spin density pattern is consistent with the predicted previously edge-state magnetism of infinite graphene nanoribbons with zigzag edges, for which it was suggested that the antiferromagnetic state may at some conditions be lower in energy than the nonmagnetic state, leading to spontaneous magnetization of such nanoribbons [64,[93][94][95][96][97]. Similar considerations for finite-length nanoribbons did not confirm this effect, suggesting that the ferromagnetic correlation of unpaired electrons at the edge of finite nanoribbons is rather an excited state property [98,99]. Nevertheless, it is remarkable that simple topological considerations performed with Clar covers and ZZ polynomials can lead to similar conclusions as those drawn from sophisticated ab initio quantum mechanical calculations. The computed ZZ spin population pattern is composed of two spin waves of opposite signs (depicted in blue and red) located at the zigzag edges of each multiple zigzag chain. The value of the spin population is visualized as a circle centered on each carbon atom. The predicted patterns very closely resemble the antiferromagnetic states of graphene nanoribbons discovered previously in quantum mechanical calculations; for more details, see text.

Conclusions
Inspired by the recently reported collection of generalized John-Sachs matrices P ZZ (B) for various classes of elementary aromatic hydrocarbons B [63], which allowed us to express the ZZ polynomial of B simply as ZZ(B, x) = |P ZZ (B)|, we have extended this approach in the current work to two important families of benzenoids: multiple zigzag chains Z(m, n) and generalized multiple zigzag chains Z k (m, n), defined graphically in Figure 2. The extension allowed us to discover compact and simple determinantal formulas for ZZ(Z(m, n), x) and ZZ(Z k (m, n), x), offering a possibility of considerable simplifications in the process of computation of ZZ polynomials for these two types of structures. The previously known ZZ(Z(m, n), x) formula, reported originally as Equation (58) in [33] and reproduced here as Equation (5), has the form of multiple sums of multiple products of complicated binomials and is too complex to shed light on the structure and classification of Clar covers of Z(m, n). In contrast, the determinantal formulas reported here hide all the algebraic complexity in the form of a determinant of a very simple object: a highlystructured Toeplitz or almost-Toeplitz matrix. In particular, we have discovered that for multiple zigzag chains Z(m, n) with even values of the structural parameter m, their ZZ polynomials have particularly transparent form given by ZZ(Z(m, n), x) = w 0 w 1 w 2 · · · w m−2 where w l is a simple polynomial in (1 + x) given explicitly by Note that w l ≡ 0 for l > n, so maximally only the lowest n + 2 diagonals do not vanish identically. An analogous formula for the ZZ polynomial of multiple zigzag chains Z(m, n) with odd values of the structural parameter m is only slightly more complicated ZZ(Z(m, n), x) = where v s l is again a simple polynomial in (1 + x) given explicitly by The corresponding formulas for ZZ(Z k (m, n), x) closely follow the same pattern; we have for odd values of m, where u t l is again a low-order polynomial in (1 + x) explicitly given by The reported determinantal ZZ polynomial formulas for Z(m, n) and Z k (m, n) also implicitly define the corresponding generalized John-Sachs matrices P ZZ (Z(m, n)) and P ZZ (Z k (m, n)) for these two families of elementary benzenoids, which upon the evaluation x = 0 reduce to the regular John-Sachs matrices P(Z(m, n)) and P(Z k (m, n)) given previously by Equations (14.33)-(14.41) in [2] (see also [62]). We want to stress that in contrast to our earlier work on generalized John-Sachs matrices [63] where most of the discovered facts were only conjectures; here, all the presented results are furnished with appropriate formal proofs. We hope that the presented here results will stimulate the graphtheoretical community to discover further examples of generalized John-Sachs matrices P ZZ also for other benzenoids, which will pave the road to conception and formulation of the generalization of John-Sachs theorem [75][76][77][100][101][102][103][104][105][106] to the world of Clar covers. We also hope that the presented here techniques will suggest an appropriate line of attack on the most difficult unsolved problem in the theory of ZZ polynomials: discovering the closed-form formula for the ZZ polynomial of hexagonal flake O(k, m, n) with arbitrary set of parameters [37,41,107].
The developed mathematical machinery of ZZ polynomials for Z(m, n) and Z k (m, n) has been applied to two practical chemical problems: determination of aromaticity patterns in finite and infinite graphene nanoribbons and in the determination of spin populations in biradical states of graphene nanoribbons. These results show good agreement with previously reported quantum chemical data, reproducing similar aromaticity patterns for Z(5, 3) and Z(4, 6) as other, well-established aromaticity indicators and similar antiferromagnetic spin patterns as those obtained from solid state calculations. In contrast to the quantum chemical calculations, the developed here methodology allows for studying the transition from the molecular to crystal regime and establishing that the infinite zigzag ribbons are completely antiaromatic. Note that similar conclusions could not be reached from a standard quantum chemical view point.

Data Availability Statement:
The data presented in this study are not available on request from the corresponding author.