Clar Covers of Overlapping Benzenoids: Case of Two Identically-Oriented Parallelograms

We present a complete set of closed-form formulas for the ZZ polynomials of five classes of composite Kekuléan benzenoids that can be obtained by overlapping two parallelograms: generalized ribbons Rb, parallelograms M, vertically overlapping parallelograms MvM, horizontally overlapping parallelograms MhM, and intersecting parallelograms MxM. All formulas have the form of multiple sums over binomial coefficients. Three of the formulas are given with a proof based on the interface theory of benzenoids, while the remaining two formulas are presented as conjectures verified via extensive numerical tests. Both of the conjectured formulas have the form of a 2× 2 determinant bearing close structural resemblance to analogous formulas for the number of Kekulé structures derived from the John-Sachs theory of Kekulé structures.


Introduction
Benzenoid hydrocarbons are an intensively studied class of aromatic molecules attracting attention of scientists specializing in various fields of chemistry including organic synthesis [1,2] and characterization [3][4][5], combustion chemistry [6,7], nanoscience [8,9], atmospheric chemistry [10,11], astrochemistry [12,13], computational chemistry [14,15], and chemical graph theory [16,17]. Enumeration and construction of Clar covers for benzenoid hydrocarbons constitutes one of the important branches of chemical graph theory. From the point of view of a chemist, a Clar cover is a generalized resonance structure that can be constructed for a given benzenoid B, in which the tetravalent character of each carbon atom is ensured by making it participate either in a two-center π bond with one of its closest-neighbor carbon atoms or in a delocalized aromatic π sextet with the other five carbon atoms in the same benzene ring. From the point of view of a mathematician, the benzenoid B is a finite plane graph embedded in a hexagonal lattice [16] and a Clar cover is a spanning subgraph of B whose components are either K 2 or C 6 (K 2 denotes a complete graph on two vertices and C 6 denotes a cycle graph of length six). To unify these two points of view, we establish the following correspondences to be used whenever convenient: a vertex in B ≡ a carbon atom in B, an edge in B ≡ a bond in B, the edge in K 2 ≡ a double bond, C 6 ≡ an aromatic ring, and an edge in C 6 ≡ an aromatic bond.
In general, a large number of Clar covers can be constructed for every Kekuléan benzenoid B. Finding this number-usually denoted by C and referred to as the Clar count-can be quite a complex task, particularly for large benzenoids of irregular shapes. The main goal of the current report is to show that for a large group of complex-shape benzenoids obtained by overlapping two identically Symmetry 2020, 12, 1599 2 of 20 oriented parallelograms, M(m 1 , n 1 ) and M(m 2 , n 2 ), the number C can be readily computed from the structural parameters m 1 , n 1 , m 2 , and n 2 of the two parallelograms M(m 1 , n 1 ) and M(m 2 , n 2 ) and the structural parameters µ and ν of the parallelogram-shaped overlapping region M(µ, ν). To derive these expressions, we need to introduce several auxiliary concepts.

Counting Clar Covers
The set of Clar covers can be naturally partitioned with respect to their order, i.e., the number of the aromatic rings C 6 used to construct a given Clar cover. Each such subset of Clar covers containing exactly k aromatic rings is referred to as the Clar covers of order k and the cardinality of this subset is denoted as c k . The set of Clar covers of order 0 is identical to the set of Kekulé structures of B, giving the identification between c 0 and the Kekulé count K. The finite nature of B (consisting of n carbon atoms) equips us with a natural upper bound n 6 for the maximal number of aromatic sextets C 6 -usually referred to as the Clar number Cl-that can be placed in B. Consequently, the number of subsets in the order partitioning is also finite and their cardinalities can be expressed by the sequence (c 0 , c 1 , . . . , c Cl ). A generating function of this sequence is usually referred to as the Clar covering polynomial or as the Zhang-Zhang (ZZ) polynomial, after the names of two Chinese mathematicians, Fuji Zhang and Heping Zhang, who first introduced it to the field of chemical graph theory about 25 years ago [18][19][20][21]. The ZZ polynomial of B decodes its most important topological invariants:  [22].
The most obvious way to compute the ZZ polynomial of B is to construct the full set of Clar covers of B and simply count the Clar covers of each order. This is most conveniently performed by a recursive assignment of single, double, and aromatic coverings to selected edges of B until the covering character of all edges in B is determined. An example of such a procedure is given in Figure 1, where the set of 25 Clar covers of the parallelogram M(2, 3) is generated in a recursive manner by forming a tree, in which the root corresponds to the original structure M(2, 3) and every leaf to a distinct Clar cover of M(2, 3). Since the number of Clar covers of order 0, 1 and 2 is 10, 12, and 3, respectively, the resulting ZZ polynomial is given by ZZ(M(2, 3), x) = 10 + 12x + 3x 2 . Unfortunately, this approach is feasible but for the smallest benzenoids. For M(10, 10) with 240 carbon atoms, the number of Clar covers comprises more than 8 millions resonance structures and it is completely impractical to represent it graphically in any way. Fortunately, ZZ polynomials possess attractive recurrence properties (see Theorems 3-6 of [18] or Properties 1-7 of [23]) allowing computing them in a recursive manner. These properties made it possible to use a recursive algorithm [18,23,24], which performs consecutive decompositions of B to compute the final ZZ polynomial recursively from the ZZ polynomials of the nodes of the decomposition tree. Our group designed an automatized computer program called ZZDecomposer [25][26][27], which can be used first for a graphical definition of the studied benzenoid by selecting a subset of cells in a hexagonal lattice and then for computation or derivation of its ZZ polynomial. Application of ZZDecomposer to M(10, 10) yields the desired result within a few seconds. (For details, see Figure 2) Another exemplary application of ZZDecomposer is the otherwise very tedious determination of ZZ polynomials of all the isomers of carbon fullerenes also presented in the recent volume of "Symmetry" [28]. Figure 1. 25 Clar covers of the parallelogram M(2, 3) can be constructed by a recursive assignment of single (no covering, symbol S), double (K 2 covering, symbol D), or aromatic ring (C 6 covering, symbol R) character to the edges of M(2, 3) marked with a red circle. Each Clar cover of order 0 (i.e., Kekulé structure, depicted in cyan frames) contributes a term 1 to ZZ(M(2, 3), x), each Clar cover of order 1 (purple frames) contributes a term x to ZZ(M(2, 3), x), and each Clar cover of order 2 (orange frames) contributes a term x 2 to ZZ(M(2, 3), x). The resulting ZZ polynomial of M(2, 3) is thus equal to 10 + 12x + 3x 2 . More importantly, ZZDecomposer can be used for discovering closed-form formulas for ZZ polynomials of whole families of isostructural benzenoids. First such formulas were discovered without ZZDecomposer by Zhang and Zhang [18][19][20] in their seminal papers for the simplest families of benzenoids: polyacenes L(n) [19], hexagonal chains [19], branched catacondensed benzenoids and catacondensed ladders [19], catacondensed all-benzenoids [19], and prolate rectangles [21]. Soon, other developments along this line appeared allowing finding closed-form formulas for other Symmetry 2020, 12, 1599 4 of 20 simple benzenoid structures: 2-and 3-tier benzenoid strips [29], one family of oblate rectangles [24], cyclo-polyphenacenes [30,31], and parallelograms [32]. In particular, the work performed by Gutman and Borovićanin [32,33] permits us to write down the formula for the ZZ polynomial of the parallelogram M(m, n) with arbitrarily large values of the parameters m and n as However, it was really the development of ZZDecomposer [25][26][27] which made such discoveries possible even for complicated structures [23,25,34,35] including generalized prolate rectangles [36], 3-, 4-, and 5-tier regular benzenoid strips [37,38], chevrons and generalized chevrons [33], zigzag-edge coronoids and fenestrenes [39], multiple zigzag chains [40], hexagons with corner defects [41], and ribbons [42]. At the moment, the only two remaining families of basic benzenoids for which a closed-form ZZ polynomial formula is not available are hexagonal graphene flakes O(k, m, n) and oblate rectangles Or(m, n), despite of substantial effort invested in such a development [34,35,37,38,41,[43][44][45]. We consider the problem of finding a closed formula for ZZ(O(k, m, n), x) with arbitrary values of the indices k, m, and n the Holy Grail of the ZZ polynomial theory. Other important aspects associated with the work on closed-form ZZ polynomial formulas concern new theoretical concepts, ideas, and connections emerging during such work and including among others a connection to transfer-matrix methodology [46], relations to cube polynomials [47][48][49], a relation to tiling polynomials [50], the derivation of interface theory of benzenoids [51][52][53], the development of the fragment and interface connectivity graphs [52][53][54], a connection to John-Sachs theory of Kekulé structures [44,45], and a discovery of classes of benzenoids for which the ZZ polynomials have only real zeros [55].

ZZ Polynomials of Complex Benzenoid Structures
We explained in the previous section that for most families of basic benzenoids, closed-form expressions for their ZZ polynomials are known [18][19][20]24,[29][30][31][32][33][36][37][38][39][40]53,54]. A topic that receives more and more attention in the chemical graph theory community is the problem of how the Clar covers of complex benzenoids built from basic benzenoids can be enumerated using as the counting blocks the Clar covers of the underlying basic benzenoids. Extensive numerical experimentation [36,44,45,56] suggests that in many cases the ZZ polynomial of the complex structure can be expressed as a function of the ZZ polynomials of the basic benzenoids, which form a kind of Gröbner basis here. The functional form of these expressions often takes on the form of a determinant, agreeing in spirit with the previously developed John-Sachs theory of Kekulé structures. We recently demonstrated formally that the ZZ polynomial of the parallelogram chain M k (m, n) formed by merging together k parallelograms M(m, n) can be expressed as a determinant of a tridiagonal, symmetric k × k Toeplitz matrix, whose diagonal elements are ZZ(M(m, n), x) and the non-zero off-diagonal elements are equal to 1 + x [56]. This result can be immediately generalized to a chain of k different parallelograms M(m 1 , n 1 ), . . ., M(m k , n k ) [56]. A similar discovery-but without a formal proof-was made for the hexagonal graphene flake O(2, m, n) [44]. It turns out that the geometrical observation that O(2, m, n) can be formed by overlapping two, slightly shifted identical parallelograms M(m, n), immediately translates to the result that ZZ(O(2, m, n), x) can be expressed as a determinant of a 2 × 2 matrix, whose diagonal elements are equal to ZZ(M(m, n), x), while the off-diagonal elements are highly structured polynomials in m, n, and x whose meaning remains to be understood. These two examples stimulated us to look for further instances of direct extensions of the John-Sachs theory of Kekulé structures to the world of Clar covers and, indeed, it was discovered that almost in every studied case we could construct a determinantal expression for the ZZ polynomial of a complex benzenoid in terms of ZZ polynomials of its simpler substructures [45], suggesting that it could be possible to formulate a general John-Sachs theory of Clar covers.
In the current report, we want to analyze a slightly different problem. Instead to taking a complex benzenoid and expressing it as a geometrical overlap of basic benzenoids, we rather select a collection of simple, basic benzenoids and form all complex benzenoids that can be created by overlapping the basic ones, and subsequently check if the ZZ polynomials of all of the resulting structures can be expressed as simple functions of the ZZ polynomials of the basic structures. In particular, due to the complexity of a general problem of this type, we analyze here all possible complex benzenoids that can be formed by overlapping two identically oriented (but not identical) parallelograms, M(m 1 , n 1 ) and M(m 2 , n 2 ). A graphical list of all possible Kekuléan structures of this type is presented in Figure 3. The next section reports closed-form formulas for the ZZ polynomials of all these structures together with derivations and justifications of these formulas (whenever available). Four kinds of structures of this type, presented in Figure 4, are obviously non-Kekuléan, as can be easily verified by checking that their number of peaks n ∧ and their number of valleys n ∨ are not identical [16,57,58] (for more detailed discussion, see for example p. 61 of [16]).

Ribbons Rb
, and Rb (n 1 , n 2 − n 1 , m 2 , m 1 ) shown in the bottom row of Figure 3 are usually referred to as generalized ribbons [58]. A closed-form formula of the ZZ polynomial of a ribbon was recently discovered [42] using recursive decompositions and interface theory considerations [51,52]. Adapting these results to the current case, we can give the following expression for the ZZ polynomials of the studied here generalized ribbons where the parameter τ assumes the following values: τ = m 2 , µ, and 0 for the three structures shown from left to right, respectively, in the bottom row of Figure 3.

Two Vertically
Overlapping Parallelograms MvM(m 2 , n 2 , m 1 , n 1 , µ, ν) Let us denote by MvM(m 2 , n 2 , m 1 , n 1 , µ, ν) the benzenoid obtained by overlapping vertically the parallelogram M(m 2 , n 2 ) located above the parallelogram M(m 1 , n 1 ) with the overlapping region M(µ, ν). Before we proceed to the determination of the ZZ polynomial of MvM(m 2 , n 2 , m 1 , n 1 , µ, ν), let us notice that three other structures of this type, MvM(n 2 , m 2 , n 1 , m 1 , ν, µ), MvM(m 1 , n 1 , m 2 , n 2 , µ, ν), and MvM(n 1 , m 1 , n 2 , m 2 , ν, µ) all share the same ZZ polynomial as their sets of Clar covers can be obtained by an appropriate symmetry operationŜ (respectively: vertical reflectionσ vert , horizontal reflectionσ hor , or in-plane rotationĈ 2 ) from the set of Clar covers of the original structure. In the following, we consider these four structures as a single equivalence class and we choose a representative (briefly denoted as MvM) for this class such that µ ≥ ν and m 1 ≥ m 2 . Since the geometrical requirements for this class of structures dictate that µ < m 1 , m 2 and ν < n 1 , n 2 , we always assume in the following that m 1 ≥ m 2 > µ ≥ ν and ν < min(n 1 , n 2 ). The generalization of the obtained formulas to structures not satisfying the conditions µ ≥ ν and m 1 ≥ m 2 will be given at the end of this Subsection.
The ZZ polynomial of MvM can be derived from the interface theory of benzenoids [51,52]. Before we proceed to the derivation, first we present basic tenets of this theory adapted here to MvM and necessary for completing this task. We introduce in MvM a system of parallel horizontal lines referred to as the elementary cuts I k in the way shown in Figure 5. Each such elementary cut is perpendicular to some vertical edges of MvM and dissects them into halves. The number of elementary cuts introduced in this way is For convenience, we augment this system with two additional elementary cuts, I 0 and I n 2 +m 2 +n 1 +m 1 −µ−ν in the way shown in Figure 5. • The set of vertical edges of MvM intersected by the elementary cut I k is referred to as the interface i k . Each edge belonging to the interface i k is referred to as an interface bond. Simple geometrical considerations allow establishing that the number of interfaces in MvM is It is beneficial to augment this set again with two additional empty interfaces i 0 and i n 2 +m 2 +n 1 +m 1 −µ−ν in the way shown in Figure 5. • The set of edges of MvM located at least partially between the elementary cuts I k−1 and I k is referred to as the fragment f k . MvM has n 2 + m 2 Interface bonds in each fragment f k are numbered from left to right. The leftmost interface edge in f k is referred to as e first ≡ e 0 and the the rightmost interface edge in f k , as e last . • Each fragment f k can be assigned an attribute of shape (W ≡ wider, N ≡ narrower, R ≡ to-the-right, and L ≡ to-the-left), which is defined in the following way W when e first ∈ i k and e last ∈ i k N when e first ∈ i k−1 and e last ∈ i k−1 R when e first ∈ i k−1 and e last ∈ i k L when e first ∈ i k and e last ∈ i k−1 Following this convention, it is possible to assign the attribute of shape to the whole structure MvM, simply by listing the shape of each fragment from the top to the bottom. For example, Let us now consider an arbitrary Clar cover C of MvM. For every interface edge e in MvM, we define a covering order function ord, which assumes the following values where E(K 2 ) and E(C 6 ) denote the sets of edges of K 2 and C 6 , respectively. We refer to a bond e with covering order 0 as a single bond, a bond e with covering order 1 as a double bond, and a bond e with covering order 1 2 as an aromatic bond.

•
The concept of covering order (or briefly: order) can be naturally extended to interfaces. We define the covering order of the interface i as • Since the interface i 0 is empty, we naturally have ord (i 0 ) = 0. The orders of the remaining interfaces can be recursively computed from the First rule of interface theory [51,52], which for an arbitrary Clar cover C relates the covering order of the interface i k to the covering order of the interface i k−1 and the shape of the fragment f k in the following way (a) If f k has the shape W, then ord(i k ) = ord(i k−1 ) + 1.
If f k has the shape N, then ord If f k has the shape R or L , then ord(i k ) = ord(i k−1 ).
• The interface orders obtained in this way are actually independent of the choice of the Clar cover C, as they are completely determined by the condition ord (i 0 ) = 0 and the shape sequence WWWWWRRNLWWWLLLNNNNNNN. Therefore, the interface orders are identical for every Clar cover C of MvM and can be treated as an inherent property of MvM allowing enumerating and constructing the set of Clar covers of MvM. The interface covering orders computed in this way are listed in red for the two structures MvM shown in Figure 5. • The number of interface bonds in every non-empty interface i k of MvM is larger by 1 from the order of this interface, ord(i k ), as can be easily seen from Figure 5. This property holds for a general structure of this type, as both the interface orders and the numbers of interface bonds in consecutive interfaces depend in the same manner on the shape of the fragment between the interfaces, except for the first and the last fragment. • An explicit formula for the interface order ord(i k ) as a function of the interface number k is somewhat cumbersome. It can be shown that where Symmetry 2020, 12, 1599 9 of 20 when k < min(m 2 , n 2 ) min(m 2 , n 2 ) when min(m 2 , n 2 ) ≤ k < max(m 2 , n 2 ) m 2 + n 2 − k when max(m 2 , n 2 ) ≤ k ≤ m 2 + n 2 − 1 0 when k > m 2 + n 2 − 1 (10) where p 1 (k), p 2 (k), and p 12 (k) correspond to the order of the same interface in the pristine structures M (m 1 , n 1 ), M (m 2 , n 2 ), and M(µ, ν), respectively. The deduced so far geometric properties and interface covering characters of MvM present us with an opportunity for designing an algorithm to determine a closed-form ZZ polynomial formula for this structure. To this end, let us consider explicitly the interface coverings of two particular interfaces of MvM, i n 2 +m 2 −µ and i n 2 +m 2 −ν . These two interfaces-intersected by the blue elementary cuts I 8 and I 9 (a) and I 7 and I 14 (b) in Figure 5-can be characterized as the interfaces containing, respectively, the leftmost and rightmost interface bonds of the intersection region M(µ, ν). Let us denote by o µ ≡ ord(i n 2 +m 2 −µ ) and o ν ≡ ord(i n 2 +m 2 −ν ) the interface orders of these two interfaces as determined from Equation (8). For the two structures MvM shown in Figure 5 o ν = min(µ, n 1 ).
Let us now construct all possible interface coverings for the interface i n 2 +m 2 −µ . Following the discussion above, we can construct o µ + 1 distinct coverings of the interface i n 2 +m 2 −µ corresponding to the partition 1 + . . . + 1 + 0 and o µ distinct coverings corresponding to the partition 1 + . . . + 1 + 1 2 + 1 2 . In the first case, the multiset of interface bonds covering orders consists of o µ 1s and one 0; the distinction between Clar covers belonging to this category concerns the location of the single bond. Similarly, in the second case, the multiset of interface bonds covering orders consists of o µ − 1 1s and two 1 2 s with the distinction between Clar covers belonging to this category concerning the location of the two consecutive aromatic bonds. Each of these 2 o µ + 1 partitions generates a system of fixed bonds in MvM not limited only to interface i n 2 +m 2 −ν but extending through both overlapping parallelograms. The systems of fixed bonds associated with each of these partitions for the structure (a) from Figure 5 is shown in Figure 6. This system of fixed bonds is a simple consequence of the Second rule of interface theory [51,52], which-figuratively speaking-ensures that consecutive covered interface bonds in each fragment f k alternate between the upper an lower interfaces, i k−1 and i k . The system of fixed bonds propagates to upper and lower interfaces of MvM in the following way forming two disconnected, wedge-shaped, non-fixed subgraphs of MvM, one located above the interface i n 2 +m 2 −ν and one below.
Consequently, the total set of Clar cover of MvM can be partitioned into 2 µ + 1 disjoint sets, each of them corresponding to one of the 1 + . . . + 1 + 0 or 1 + . . . + 1 + 1 2 + 1 2 partitions. The ZZ polynomial of MvM can be then constructed as a sum of ZZ polynomials of the non-fixed disconnected components of MvM over all the partition classes where the additional factor x in the second line of Equation (14) symbolizes the aromatic ring in the fixed component of MvM for each 1 + . . . + 1 + 1 2 + 1 2 partition class. Since both lines of this formula are isostructural, we are able to write Equation (14) in even more compact form as In the case when i n 2 +m 2 −µ = i n 2 +m 2 −ν (i.e., when µ > ν), to obtain a formula for the ZZ polynomial of MvM analogous to Equation (14), we need to consider all the combinations of the coverings of the interfaces i n 2 +m 2 −µ and i n 2 +m 2 −ν simultaneously. Two important aspects of this process need to be addressed: (1) Not every covering of i n 2 +m 2 −µ is compatible with every covering of i n 2 +m 2 −ν (and vice versa); the precise compatibility conditions and their origin are discussed below. (2) A selection of compatible coverings of both interfaces induces a system of fixed bonds in MvM, which separates three disconnected, parallelogram-shaped, non-fixed subgraphs of MvM. The actual dimensions of the parallelograms, depending on the positions of the single bond(s) and/or aromatic ring(s) in both interface coverings, are derived below. As an illustration we show in Figure 7 four different systems of fixed bonds (depicted in gray) and non-fixed bonds (depicted in black) associated with four selected covering types of the interfaces i 8 and i 14 in MvM formed by overlapping vertically M(11, 6) with M (11,9) by M(9, 3). Note that in the leftmost covering, the central parallelogram located between both interfaces formally corresponds to M(0, 5) and allows only a single Clar cover; therefore it can be treated either as a fixed-bond region or non-fixed-bond region. Four examples of partially covered structures MvM (11,6,11,9,9,3) show that the non-covered regions consists of three parallelograms with their shapes determined by the location of the single and/or aromatic bonds in the both covered interfaces i n 2 +m 2 −µ and i n 2 +m 2 −ν .
In this paragraph we discuss the compatibility conditions of both interface coverings. As already mentioned above, selecting a particular covering of i n 2 +m 2 −µ induces a system of fixed bonds in MvM, which extends vertically to the other interface, i n 2 +m 2 −ν , imposing double bond covering of some of its bonds. The available range of non-fixed interface bonds in i n 2 +m 2 −ν can be characterized as follows.

•
Let us assume that a 1 + . . . + 1 + 0 partition was selected for i n 2 +m 2 −µ with the single bond in position k ∈ [0, . . . , o µ ]. Somewhat involved geometric considerations show that the range of indices of non-fixed interface bonds in i n 2 +m 2 −ν associated with this choice is given by Consequently, a single bond in interface i n 2 +m 2 −µ in position k permits placing a single bond in Let us assume now that a 1 + . . . + 1 + 1 2 + 1 2 partition was selected for i n 2 +m 2 −µ with the aromatic ring in hexagon k ∈ [1, . . . , o µ ]. Again, geometric considerations show that the range of indices of non-fixed interface bonds in i n 2 +m 2 −ν associated with this choice of covering for i n 2 +m 2 −µ is given Consequently, an aromatic ring in hexagon k of interface i n 2 +m 2 −µ permits placing a single bond in interface i n 2 +m 2 −ν in position l ∈ [max(k + o ν − µ, 0), o ν − max(0, ν − k + 1)] or permits placing an aromatic ring in hexagon l ∈ [1 Let us now give the dimensions for the three parallelograms induced by the systems of fixed bonds in MvM associated with placing a single bond (i = 0) or an aromatic ring (i = 1) in position k in the interface i n 2 +m 2 −µ and a single bond (j = 0) or an aromatic ring (j = 1) in position l in the interface i n 2 +m 2 −ν . It is possible to demonstrate via elementary geometric considerations that in this situation the upper parallelogram M u ≡ M upper , the lower parallelogram M l ≡ M lower , and the middle parallelogram M m ≡ M middle are given by the following expressions Consequently, the ZZ polynomial of MvM with µ > ν is given by the following expression where explicit expressions for ZZ(M u , x), ZZ(M m , x), and ZZ(M l , x) are given by Equation (2). Substituting these expressions into Equation (19) produces a formula with a septuple summation and six binomial coefficients, which is not given explicitly here.
As mentioned earlier, the structure MvM with µ ≥ ν and m 1 ≥ m 2 , for which the ZZ polynomial formula is given by Equation (19) (when µ > ν) and by Equation (15) (when µ = ν), is a representative for an equivalence class of four structures sharing the same ZZ polynomial. We have thus two possible pathways leading to the determination of the ZZ polynomial for a general structure MvM(m 2 , n 2 , m 1 , n 1 , µ, ν), for which the conditions µ ≥ ν and m 1 ≥ m 2 are not satisfied.
The first natural way is to identify the equivalence class representative MvM for MvM(m 2 , n 2 , m 1 , n 1 , µ, ν) by producing four symmetry-related structures (Ê(MvM(m 2 , n 2 , m 1 , n 1 , µ, ν)),Ĉ 2 (MvM(m 2 , n 2 , m 1 , n 1 , µ, ν)), σ hor (MvM(m 2 , n 2 , m 1 , n 1 , µ, ν)), andσ vert (MvM(m 2 , n 2 , m 1 , n 1 , µ, ν)), where D 2 = {Ê,Ĉ 2 ,σ hor ,σ vert } is the 2D dihedral group) and selecting a structure for which µ ≥ ν and m 1 ≥ m 2 . Perhaps such a solution is aesthetically pleasing, but from a practical standpoint it is more convenient to use a symmetry-adapted version of Equations (15) and (19) given by the following expression where MvM stands already for an arbitrary structure MvM(m 2 , n 2 , m 1 , n 1 , µ, ν) and where the following symbols were used In case of formulas of high complexity, such as those given by Equations (20)- (29), it is important to make sure that no typographical errors are present in the final equations. To this end, we show in Figure 8 a MAPLE procedure MvM for computing the ZZ polynomial of an arbitrary structure MvM(m 2 , n 2 , m 1 , n 1 , µ, ν) together with several MvM structures and their ZZ polynomials computed in a twofold way: black text reports the ZZ polynomials computed in a brute force manner using ZZDecomposer and the blue text in red frames reports the ZZ polynomials computed using the provided MAPLE procedure. Note that the computational time needed to evaluate the ZZ polynomials of larger structures using recursive decompositions with ZZDecomposer exceeds a few minutes, while the MAPLE procedure produces identical results instantaneously.

Two Horizontally
Overlapping Parallelograms MhM(m 1 , n 1 , m 2 , n 2 , µ, ν) A closed-form formula for the ZZ polynomial of the composite benzenoid MhM(m 1 , n 1 , m 2 , n 2 , µ, ν) shown schematically in Figure 3 has not been reported yet. Extensive numerical experimentation shows that the appropriate formula can be written as a determinant of a certain 2 × 2 matrix ZZ(MhM(m 1 , n 1 , m 2 , n 2 , µ, ν), x) = t 1 u v t 2 (30) whose diagonal elements t 1 and t 2 correspond to the ZZ polynomials of the parallelograms M (m 1 , n 1 ) and M (m 2 , n 2 ), respectively, and the off-diagonal elements u and v are given explicitly by the following expressions We have no formal proof demonstrating the correctness of Equation (30), but at the same time the overwhelming numerical evidence clearly shows that Equation (30) is correct and applicable to structures MhM(m 1 , n 1 , m 2 , n 2 , µ, ν) with arbitrary values m 1 , m 2 > µ ≥ 0 and n 1 , n 2 > ν ≥ 0 of the structural parameters. It might be instructive to describe briefly the thought process leading to the discovery of Equation (30). In 1952 Gordon and Davison observed that there is a one-to-one correspondence between the number of Kekulé structures for a given benzenoid B and the number of monotonic path systems in B [59]. This results was rigorously proved in 1984 by Sachs [60]. Subsequently, John, Sachs, and Rempel [61,62] and later Cyvin and Gutman [57,58,63] designed a simple strategy allowing computing the Kekulé count of B by taking a determinant of the so-called John-Sachs path matrix. The John-Sachs path matrix is a square matrix of size p × p, where p is the number of peaks (or equivalently, valleys) in the Kekuléan benzenoid B. The rows of the matrix are associated with the peaks of B and the columns, with the valleys of B. The (i, j) entry of the matrix is the Kekulé count for the subgraph of B containing only those vertices of B that can be reached from the peak p i by going exclusively downward and simultaneously that can be reached from the valley v j by going exclusively upward. For details, we refer the reader to the excellent exposition with multiple examples given by Gutman and Cyvin [63] or by Bodroža and collaborators [43]. For the studied here composite benzenoid MhM(m 1 , n 1 , m 2 , n 2 , µ, ν), the John-Sachs path matrix is of size 2 × 2, because this structure contains 2 peaks and 2 valleys. The diagonal elements of this matrix correspond to the Kekulé counts of M (m 1 , n 1 ) and M (m 2 , n 2 ), and the off-diagonal elements, to the Kekulé counts of M (m 1 + m 2 − µ, ν) and M (µ, n 1 + n 2 − ν). We asked ourselves whether it is possible to extend the John-Sachs path matrix construction to the problem of determination of ZZ polynomial for this structure and to our surprise we discovered that indeed such an extension can be constructed. The diagonal entries of the original John-Sachs matrix are naturally replaced by the ZZ polynomials of the (p i , v i ) subgraphs (i.e., by ZZ (M (m 1 , n 1 ) , x) and ZZ (M (m 2 , n 2 ) , x)), which are given explicitly by t 1 and t 2 in Equations (31) and (32). Discovering the actual analytic form of the off-diagonal elements u and v requires quite an extensive experimentation, but the expressions given by Equations (33) and (34) are not surprising and show a clear connection to the (p i , v j ) subgraphs (i.e., to M (m 1 + m 2 − µ, ν) and M (µ, n 1 + n 2 − ν), even if the details of this connection are quite murky at the time and the presence of the multiplicative factor (1 + x) and the positive and negative shift of 2 in the parallelograms dimensions escapes our understanding. The presented here example of extending the John-Sachs theory of Kekulé structures to the world of Clar covers suggests that there might exist a generalization of the John-Sachs theorem valid not only for monotonic paths consisting of the K 2 graphs, but also for generalized monotonic paths consisting of the K 2 and C 6 graphs. Note that once such a general theory of general monotonic path systems is conceived, it will provide a proper justification to the conjecture expressed by Equation (30).

Two
Intersecting Parallelograms MxM(m 2 , n 2 , m 1 , n 1 , µ, ν) A union MxM(m 2 , n 2 , m 1 , n 1 , µ, ν) of two intersecting parallelograms M (m 1 , n 1 ) and M (m 2 , n 2 ), shown schematically in Figure 9, can be formally treated as a union of two overlapping ribbons, Rb (n 1 , ν, m 2 , µ) and Rb (n 1 , n 2 − n 1 − ν, m 2 , m 1 − m 2 − µ), where the first ribbon was rotated by π with respect to the standard orientation [42,58]. This second interpretation allows us to use exactly the same technique as in the last subsection to find the ZZ polynomial formula for MxM(m 2 , n 2 , m 1 , n 1 , µ, ν). We again seek for an appropriate ZZ polynomial formula in the form of 2 × 2 determinant, with the diagonal elements t 1 and t 2 corresponding to the ZZ polynomials of ribbons Rb (n 1 , ν, m 2 , µ) and Rb (n 1 , n 2 − n 1 − ν, m 2 , m 1 − m 2 − µ) and the off-diagonal elements u and v somehow related to the two intersecting parallelograms M (m 1 , n 1 ) and M (m 2 , n 2 ). Numerical experimentation shows that indeed the sought ZZ polynomial can be represented in the anticipated way as ZZ(MxM(m 2 , n 2 , m 1 , n 1 , µ, ν), with u, v, t 1 , and t 2 given explicitly by the following expressions Again, we have no formal proof demonstrating the correctness of Equation (35) but at the same time multiple numerical checks convince us that Equation (35) is correct and applicable to all the composite benzenoids MxM(m 2 , n 2 , m 1 , n 1 , µ, ν) with the structural parameters satisfying the following conditions: µ, ν > 0, n 2 > n 1 + ν, and m 1 > m 2 + µ. Following the discussion at the end of Section 4.3, we believe that Equation (35) will remain a conjecture until a theory of general monotonic path systems (i.e., an extension of John-Sachs theory applicable to Clar covers) is developed.  Figure 3 can be formally treated as a superposition of two parallelograms, M (m 1 , n 1 ) (red shading) and M (m 2 , n 2 ) (blue shading), or two ribbons, Rb (n 1 , ν, m 2 , µ) (red shading) and Rb (n 1 , n 2 − n 1 − ν, m 2 , m 1 − m 2 − µ) (blue shading).

Conclusions
We present a complete set of closed-form formulas for the ZZ polynomials of composite benzenoid structures obtained by overlapping two parallelograms, M(m 1 , n 1 ) and M(m 2 , n 2 ). Depending on the mutual location of the parallelograms and the dimension of the overlap parallelogram-shaped region, M(µ, ν), the resulting composite structure can be Kekuléan (see Figure 3) or non-Kekuléan (see Figure 4). We identify five families of Kekuléan structures: parallelograms M(m 1 + m 2 − µ, n 1 = n 2 ) and M(m 1 = m 2 , n 1 + n 2 − ν) with ZZ polynomials given by Equation (2), generalized ribbons Rb with ZZ polynomials given by Equation (4), horizontally overlapping parallelograms MhM(m 1 , n 1 , m 2 , n 2 , µ, ν) with ZZ polynomials given by Equation (30), intersecting parallelograms MxM(m 2 , n 2 , m 1 , n 1 , µ, ν) with ZZ polynomials given by Equation (35), and vertically overlapping parallelograms MvM(m 2 , n 2 , m 1 , n 1 , µ, ν) with ZZ polynomials given by Equation (20). All formulas are given in the form of multiple sums over binomial coefficients, bearing a close structural similarity to the ZZ polynomial of a pristine parallelogram M(m, n) given by Equation (2). In case of the subset of the composite benzenoids consisting of parallelograms M, generalized ribbons Rb, and vertically overlapping parallelograms MvM, we demonstrate formally the correctness of the derived formulas. For intersecting parallelograms MxM and horizontally overlapping parallelograms MhM the derived formulas are given without a proof, in the form of conjectures; pervasive numerical evidence suggests correctness of these conjectures. Since both the conjectured formulas have the form of a 2 × 2 determinant with a striking resemblance to the John-Sachs determinantal formulas for the number of Kekulé structures [57,58,[61][62][63], we expect that their formal proof could be possible only once the John-Sachs theory of Kekulé structures is generalized to the world of Clar covers. We hope that our results will stimulate the graph theoretical community to accomplish this task and simultaneously will provide an appropriate testing ground for such a generalization. The formal demonstration of correctness of the ZZ polynomial formulas for generalized ribbons and vertically overlapping parallelograms is based on the recently announced interface theory of benzenoids [51][52][53], which seems to be a promising tool for solving numerous challenging problems in chemical graph theory. All the derived and reported here ZZ polynomial formulas were tested against numerous instances of actual composite structures of a given type, for which the ZZ polynomials were computed in a brute-force decomposition manner using ZZDecomposer [25][26][27].