Quantum walks in polycyclic aromatic hydrocarbons

Aromaticity is a well-known phenomenon in both physics and chemistry, and is responsible for many unique chemical and physical properties of aromatic molecules. The primary feature contributing to the stability of polycyclic aromatic hydrocarbons is the delocalised π-electron clouds in the 2p z orbitals of each of the N carbon atoms. While it is known that electrons delocalize among the hybridized sp 2 orbitals, this paper proposes quantum walk as the mechanism by which the delocalization occurs, and also obtains how the functional chemical structures of these molecules arise naturally out of such a construction. We present results of computations performed for some benzoid polycyclic aromatic hydrocarbons in this regard, and show that the quantum walk-based approach does correctly predict the reactive sites and stability order of the molecules considered.


I. INTRODUCTION
The structure and properties of Benzene and other arenes has been a subject of special interest in the field of quantum chemistry for a long time.It is now known that each of the six carbon atoms in Benzene have their orbitals hybridized into the sp 2 state, and the 2p z orbitals host a delocalized π-electron cloud of the molecule.This type of cyclic conjugated system has been a subject of significant interest in the chemical sciences, and the description of such systems is studied under the concept(s) of (anti)aromaticity [1][2][3][4][5].These molecules are characterized by their electron-rich clouds, which stabilize the structure through cyclic resonant structures.While it is known that the delocalization does occur in the electron cloud, the underlying physical process by which the process takes place is an open question as of yet.
In order to understand the chemical characteristics of the cyclic conjugated system, bond order of the bonds in the system is studied.Bond order is a widely used concept in chemical sciences, and works as a tool to help predict chemical behaviour.Certain reactions can take place only for bonds of certain order (e.g.electrophilic and nucleophilic substitution reactions can only occur for bond order ≥ 1).It can also help in explaining the reactive tendencies of some molecules by consideration of the sum of bond order [6] of the atoms preferred by the molecule.The strength of a chemical bond depends on two thingsthe degree of overlap between the interacting orbitals, and the difference in energies of the atomic orbitals involved in bonding, which is reflected by bond polarity.The bond polarity is used to qualify the ionic nature of the bond, while the degree of overlap quantifies its covalent nature [7,8].
A comprehensive method to calculate bond order from first principles was only available in 2017 [9].As we investigate polycyclic aromatic rings in this work, we can use a slightly simpler method that relies on the system's geometry.It requires a basis relative to which bond order can be calculated.We make the traditional choice of C − C single bond in Ethane (H 3 C − CH 3 ) to have bond order 1, and the C = C double bond of Ethene (H 2 C = CH 2 ) to have a bond order of 2. This gives us a quantity known as Relational Bond Strength Order (RBSO) [10], and provides a robust description of bond strength as long as the Born-Oppenheimer approximation is valid.This provides a useful tool to study the bonds in an aromatic molecule, as a bond with bond order closer to 1 would have characteristics that are closer to a single bond, but in a bond with order closer to 2, the double-bond characteristics will dominate.
With advancements in the field of quantum information and computation, a quantum walk has emerged as one of the most efficient ways to model the controlled dynamics of quantum states and quantum particles [11][12][13][14][15].The quantum walk can be regarded as a quantum counterpart of random walks which have served as an efficient way to model dynamics defined by classical physics.Therefore, it is natural to explore the potential of quantum walk to model quantum dynamics in the range of physical systems where quantum physics plays a defining role.It is described in two main formalisms -the discrete-time quantum walk (DTQW) and the continuous-time quantum walk (CTQW).The dynamics of DTQW require a coin and a position Hilbert space in order to be defined, however, CTQW dynamics can be defined with only a position Hilbert space.A quantum walk spreads quadratically faster than a classical walk on its position space and its dynamics can be engineered using a set of evolution parameters, and therefore it has been used as a basis for design and implementation of quantum algorithms and quantum simulations [16][17][18][19][20][21][22][23][24][25][26][27][28][29], to study problems such as graph isomorphism [30], quantum percolation [31][32][33], and to develop schemes for implementation of universal quantum computation [34,35], among others.Quantum walks are thus very versatile tools indeed, and their practical significance has been demonstrated by way of implementation in many quantum systems, such as NMR [36], integrated photonics [37][38][39], ion traps [40,41], and cold atoms [42].
Modelling the dynamics of the FMO complex using quantum walks is one of the important works reported [22,43] in the direction of using quantum walks for modelling quantum dynamics in physical systems.However, not much progress has been reported in modelling of dynamics in chemical systems beyond FMO complex despite various algorithms having been proposed for using quantum walks for several computational tasks.Quantum simulation of dynamics in chemical complexes is one of the promising applications envisioned using quantum computers, and quantum algorithms using quantum walks could play an important role in that direction.
In this work, we will make use of both the discrete-time and continuous-time quantum walk formalisms in order to study the structure, stability, and the relative chemical reactivity of different sites when exposed to an electrophile for each of the molecules considered.We also qualitatively establish that the order of stability of the molecules arises naturally from a quantum walk-based framework.We thus propose the hypothesis that the electrons in the 2p z orbitals delocalize in the π-electron cloud via a quantum walk.
The paper is divided into six sections.In Sec.II, we discuss the discrete-and continuous-time quantum walks, and elucidate briefly the specific variants of the quantum walks used in this work.A quick recap of the concepts of relative bond strength order is given in Sec.III, and the methods we use to model the dynamics in aromatic hydrocarbons are detailed in Sec.IV.We present our results in Sec.V, and conclude in Sec.VI.

II. QUANTUM WALK A. Continuous-time quantum walk
A quantum walker performing CTQW has its position space defined by a graph Γ(V, E), where V is the set of its vertices and E the set of edges.Given the cardinality of

2
. Let A be the matrix defined as, A, defined as such on Γ, is known as the adjacency matrix of the graph, and offers a representation of Γ itself.It is a real-valued matrix, and the vertices are labeled by the computational basis states {|1 , |2 , ..., |N }.The quantum state of the entire graph is defined at an arbitrary point of continuous time t by the wavefunction |ψ(t) , which is defined as, The CTQW is a quantum process, and thus obeys the Schrödinger equation, while its classical counterparts obey the Markovian master equation.The Hamiltonian H Γ used in the expression of the Schrödinger equation is given by, Here D is a diagonal matrix, and each diagonal element d jj corresponds to the degree of vertex j. γ is a constant, A is the adjacency matrix (defined as in Eq. ( 1)), L is known as the Laplacian of Γ, and γ is the rate of transition for the graph.It may be observed from the expression in Eq. ( 4), that H Γ is independent of time, and thus the solution to Eq. ( 3) is given by the time evolution operator U , and may be expressed as, where the expression of ψ is in the units of .Since the adjacency matrix A is always real and symmetric, due to it representing an undirected graph, so is the Hamiltonian -which ensures that U is necessarily a unitary operation.The coefficients γ in the Hamiltonian in Eq. ( 4) are rates of transition between different nodes.Traditionally, in unweighted graphs, the rates are all normalized to 1.In the case of arenes, however, we use the values of RBSO as mentioned in the beginning of this paper in Sec.I, and formally defined in Sec.III as the values of the weights of the graph.This corresponds to how easily the electrons are able to move across the bonds, and physically represents a measure of the bond length, as bonds of higher order are shorter.

B. Discrete-time quantum walk
The quantum evolution of a walker executing a DTQW on a one-dimensional lattice may be described on a Hilbert space H = H C ⊗ H P , where H C and H P are coin and position Hilbert spaces, respectively.The coin space in this case is a 2-dimensional space, and thus its basis set can be assumed to consist of two vectors |↑ , |↓ , such that they are mutually orthogonal.The coin space represents a Hilbert space which is internal to each walker.The position space is an infinitedimensional Hilbert space with the basis chosen to be the set { | x | x ∈ Z }.Each vector |n represents the site x = n in the position space.
The initial state of the walker is therefore represented as a tensor product of its states in the two spaces, and may be written in the form shown in Eq. (6).
Here, α, β ∈ C are amplitudes of the walker's internal coin states.The evolution is defined as a unitary operation in the coin space, followed by a coin-dependent shift operator in the position space.Both the operators are unitary operations, and a typical DTQW evolution with a single-parameter coin is described as, where, , and Where (a, b) ∈ R represent the amount of shifting in the position space experienced by the component in the eigenspaces corresponding to |↑ and |↓ respectively.As can be seen from Eq. ( 7) the shift operation effects the traversal of the components of the probability amplitude in different directions.In order to obtain the order of reactivity of sites, the algorithm in [28] requires a variant of DTQW known as a directed DTQW (D-DTQW) [45].In case of the D-DTQW, the shift operator only allows a single component of the probability amplitude to traverse the graph.Thus depending on which component is allowed to traverse the graph, the D-DTQW shift operation may be defined in one of the two ways, The coin operation for the node-ranking algorithm has the same form as C θ , however, it uses a position dependent form which may be expressed as, This operation is also a special unitary matrix, similar to the C θ defined for the DTQW in Eq. (7).In case of the procedure to arrange the nodes in order of reactivities, we use a specialized coin operation of the form defined in Ref. [28], were α x is the proportion of the incoming weight with respect to the total incoming and outgoing weights at the node represented by |x .In this case, the incoming and outgoing weights are equal as the graph is undirected, hence the expression for α x may be simplified to α x = dx 2 , where d x represents the degree of node |x .The shift operator used to implement the algorithm is defined as where k ∈ Z, and U kx is a unitary matrix that restricts the walker to Markovian jumps between specific nodes in the position space.Since the coin operator rotates between the 'stored (ψ ↑ c )' and 'travelling (ψ ↓ c )' components of the probability distribution , it results in different amounts of storage at different nodes, depending on the amount of information that passes through them, generating a ranking of the nodes.The position space in this case is the graph Γ(V, E) with the set of nodes V and edges E. The adjacency matrix A of this graph is defined as a ij = RBSO(i, j), where RBSO(i, j) is the relative bond strength order of the bond between the i th and j th carbon atoms in the molecule under consideration, as defined in Sec.III.

III. RELATIVE BOND STRENGTH ORDER
Many computational approaches to calculate bond order exist, however, most have severe fundamental limitations.Approaches requiring categorization of electrons to be spin-up or spin-down fail to achieve universality as they cannot describe noncollinear spin magnetism.Some methods consider bond order to be an explicit functional of the total electron density and spin magnetization density functions.However, in the limit of a complete basis set, the density matrix is an overcomplete representation of the distribution of electron density.Thus, a functional of the density matrix may not necessarily be a functional of the electron density.This causes the bond order results to be inconsistent across different quantum chemistry methods, even if they yield the same electron density, spin magnetization density, and energy [44,46].This dependence renders this particular kind of formulation unphysical.
Other approaches that can be used in its stead include Wheatley-Gopal and Laplacian correlations between overlap and bond order [47,48], Mayer Bond Index [49] and First Order Delocalization Index [50] applied to density-based charge partitions, Natural bond orbital [51], Adaptive natural density partitioning [52,53], among others.
In this work, we shall consider electron delocalization in benzoid polycyclic aromatic systems, specifically the C−C conjugated double bonds in Benzene, Naphthalene, Anthracene and Phenanthrene.Since all the bonds under consideration are C − C bonds, we can use a relative bond strength order [54] (hence called bond order) here.We determine the bond order by using experimentally determined vibrational frequencies of the bonds [55][56][57][58] to calculate the force constant matrix.We then find the local modes of vibration that correspond to specific bonds, and thus obtain local stretching modes [59,60], which can be used to design a sensitive measure of bond strength.
Briefly, the method proceeds as follows.We calculate the force constant matrix corresponding to the complete set of vibrational frequencies of the molecule under consideration by Wilson's GF method [61].The internal degrees of freedom ({ q 1 , ..., q 3N −6 }) describing the potential energy surface (PES) in an optimal manner are often nonlinear, so it is assumed that the displacements with respect to the internal coordinates are small.This enables linearization of the internal coordinates as the set { Q i }, where i = 1, .., 3N − 6.A PES V can be expanded in a Taylor series around its minima in terms of { Q i }, and the force derivative matrix F is then given by the Hessian of V .The first term in the Taylor expansion is adjusted with the zero point energy and the second term vanishes due to the evaluation at minima.Thus, we obtain the relation, which can then be compared with the classical vibrational kinetic energy of the form, where g ij is an element of the metric tensor of the internal curvilinear coordinates, and Ȧ ≡ ∂A ∂t .Evaluating the metric tensor g in the minimum q of the PES V gives the Wilson's G-matrix, This leads us to Wilson's equation, given by where F is the force constant matrix in terms of the internal coordinates ({ q i } 3N −6 i=1 ), D is the matrix of vibrational eigenvectors d µ , each of which forms one of its columns, G is the Wilson matrix from Eq. ( 13), and Λ is a diagonal matrix consisting of the eigenfrequencies (normal modes) corresponding to d µ .
Diagonalizing F with D, we obtain the force constants corresponding to each local mode as, In this work, we use the vibrational frequencies as listed in [55] to calculate our force constants.According to the extended version of Badger rule [62], the bond order can be calculated from a power relationship between it and local mode corresponding to each bond.We use the CC bonds in Ethane and Ethene as having bond orders 1 and 2, respectively, which fixes the relationship between bond order and the local vibrational mode force constants as, [54] It may be noted that the bond order of a bond that does not exist in the molecule is by definition considered to be zero.In Table I, we show the bond order values corresponding to the different C − C bonds as per Fig. 1.

IV. METHODS
In this section, we consider the behavior of the π-electrons in the aromatic system.The behavior of the delocalised electrons in the conjugated system is modeled in the form of quantum walk.We will use both, the CTQW and DTQW formalisms to model the electrons' behavior.We will use the CTQW to model the delocalization process in order to obtain results relating to delocalization modes and stability of the molecule, and the DTQW is employed in an algorithmic form defined in [28] to rank sites of the molecule in order of their reactivity towards an electrophile.
In both the quantum walk formalisms, the molecules we have considered are modeled as a graph Γ(V, E), where V is the set of its nodes and E the set of edges.Only the bonds and carbon atoms that form the conjugated system are included in the  The atom numbers in the 'Carbon atoms' column correspond to the respective carbon atoms of the molecules, marked as shown in Fig. 1.
sets E and V , respectively.The graph is represented in the form of its adjacency matrix, and the weight of edge is defined to be its bond order, as it represents the ease by which an electron may traverse the graph, thus qualitatively describing the ease by which an electronic wavefunction may delocalize over the bond.The graphs corresponding to each molecule are illustrated in Fig. 1.The π-electrons are modeled as independent, noninteracting quantum walkers that are free to delocalize over this graph.The interelectronic interactions between delocalized π electrons in a polycyclic aromatic system are primarily nearest-neighbor pairwise interactions.They are taken into consideration when the bond order is calculated, as the electrons find it easier to delocalize over a shorter bond, i.e., a bond with a higher bond order in general.Thus, the quantum walk dynamics of the electrons can be assumed independent of each other, subject to the constraint that edge weights must take into account the effect of appropriate pairwise interactions along the network.
In the CTQW study, we consider the state of the system after some time t, and look at the probabilities of finding some electron at a particular position and its variation with time.This gives us a general idea of the dynamics of the delocalizing electronic wavefunctions.We also look at how the maximum probability of any electron to exist at a particular position, hereafter called MAXP, varies over time.It is defined for each position basis vector as the maximum probability of a single electron to exist at that particular basis vector in the position space.A high mean value here indicates that one of the π-electrons is regularly found here, and therefore is evidence of localization, indicating that at least one of the π-bonds connecting the considered vertex has a higher bond order, i.e., it is not a dominant part of the delocalization mode of the molecule, and thus has a higher double-bond character.A lower value, on the other hand, corresponds to the fact that all the bonds at the considered vertex are involved in delocalization.This study helps to ascertain the bond delocalization mode of each molecule out of several possible choices, as illustrated in Fig. 2. We also consider the truncated mean of this data as a variable called TRP, in which we discard the highest and lowest values of the probability, and calculate the mean of the remaining probabilities at a particular position over time.This measure discards the outlier(s) created by electrons that may localize, and provides a reasonable qualitative description of where the wavefunction is likely to exist.This provides a way to qualitatively characterize the reactivities of the various molecules, and also a way to arrange them in order of their reactivities.A higher value of TRP implies the electronic probability distribution often has a peak at certain places, which provides an estimate of the likelihood of an electron being present at a certain point in the chain -thus qualitatively estimating the availability of that site to form chemical bonds.
The DTQW formalism characterizes the symmetries of the considered graph with respect to a diffusing quantum particle.The DTQW-based algorithm used [28] requires a coin Hilbert space mapped to the connections in structure in addition to the position Hilbert space for its implementation, and estimates the ability of each site to accumulate information as a quantum particle diffuses across it.This provides a qualitative estimate of site activities, and thus enables a qualitative overview of the behavior of the conjugated system in the presence of an electrophile.

A. Evolution of probability with time
In this subsection, we take a look at the evolution of the probability of each electron to exist at different points for each molecule as a function of time.This shows a oscillatory behavior symmetric about the initial state in the case of Benzene, one period of which is shown in Fig. 3.It is clear from Fig. 3 that over time, the walkers diffuse on the network representing Benzene, as illustrated in Fig. 1(a) in an oscillatory manner, with a period of 12.60.The total probability of finding an electron at a particular position is always unity, however, as can be seen from Fig. 3, it is always uncertain which electron is detected, in accordance with the concept of indistinguishability of electrons.The probabilities corresponding to each electronic wavefunction are also verified to sum to unity over all positions.This may be derived from the fact that the results for this calculation may be expressed as a bistochastic N × N matrix B. The j th row of B represents the probability of the j th walker to be found at each of the N sites.
Unlike the behavior for Benzene, the walkers on the network representing Naphthalene do not display any oscillatory behavior, and only partially localize in time.There is significantly more delocalization as compared to Benzene, however, the total probability of finding an electron at any position point is always unity.This is clearly illustrated in Fig. 4.
Similar to the case of Napthalene, the electrons in the case of Anthracene and Phenanthrene also do not have any oscillatory behavior.A few snapshots of the state of the two systems are shown in Figs. 5 and 6.It may be seen qualitatively by inspection that the walkers in the case of Phenanthrene (Fig. 6) display significantly more mixing in the position basis than the case of Anthracene (Fig. 5), thus alluding to the fact that Phenanthrene is more stable than Anthracene.This idea is further developed and illustrated in Fig. 10.

B. Reactivity of sites
In this subsection, we look at the results obtained upon application of the DTQW-based algorithm.This analysis takes into consideration a single quantum walker diffusing over the entire network through a directed variant of the DTQW, and generating an ordered arrangement, i.e. ranking, of the nodes.This ordered arrangement is based on the amount of information passing through each node.The algorithm used is invariant of the starting position of the walker in case of a finite network over sufficient run-times, and thus it is only necessary to run it for a single walker.
The sites in Fig. 7 are numbered according to the scheme illustrated in Fig. 1.From Fig. 7(a), it is apparent that each of the Benzene sites are equivalent, and hence an electrophilic substitution is equally likely to occur at any of the carbon atoms.This is the same result as obtained by symmetry -the Benzene molecule has D 6h symmetry, and thus each of its carbon atoms and the bonds in the network are equivalent.This is not the same case, however, for Naphthalene (three unique sites, D 2h symmetry group), Anthracene (four unique sites, symmmetry group D 2h ) and Phenanthrene (six unique sites, C 2v symmetry group).In each case, a lower rank for a particular implies a higher probability of an electrophilic substitution occurring at that site.Physically, this corresponds to the fact that the position most likely to be substituted by an electrophile has the least information passing through it.The reaction, therefore, tends to occur via the pathway causing the least loss of information in the molecule, as is expected by the second law of thermodynamics.
These observations are also in line with the fact that an electrophile would preferentially create a substituted aromatic compound at site C3 for Naphthalene (or equivalently, at sites C5, C8 or C10).It also explains naturally the tendency of Anthracene to form substituted compounds at equivalent locations C5 and C12, as well as the high tendency of electrophilic substitutions at C11 and C12 in the case of Phenanthrene.

C. Delocalization modes
In this section, we take a look at the the variation of maximum probability of an electron to exist at a particular site (MAXP), for each unique site in the molecule.For the Benzene molecule, only the vertex 1 is unique, and all others are equivalent sites.In case of the Napthalene, vertices 4, 6 and 10 are unique.Anthracene and Phenanthrene, despite having the same number of sites, have different numbers of unique sites due to their structures.In case of Anthracene, sites 2, 5, 11 and 14 are unique, while in case of Phenanthrene, the unique sites are the ones labeled 1, 2, 3, 9, 10 and 13.In Fig. 8, we show a plot showing how the quantity MAXP varies with time.We have only plotted this for positions that are unique.Every equivalent position has the same curve, and the same mean.The values are sampled at a time interval of every 0.01 units.
From Fig. 8, we observe that the plot for Benzene (Fig. 8(a)) shows a periodic behavior, while the others do not show any periodicity as such.The mean values for each position are plotted in Fig. 9. Through a representation of the mean value of MAXP measured over time, Fig. 9 essentially presents a way of looking at the bond delocalization mode of the different bonds in the molecule.A vertex with a high mean value implies there is at least one higher-order bond at the corresponding site on the molecule, while a vertex with a lower value implies all the bonds at the corresponding site have delocalized.This enables us to verify the delocalization mode for each of the molecules.Benzene has only one possible pattern shown in Fig. 2(a), but Naphthalene prefers to exist in the pattern depicted in Fig. 2(c).Similarly, Anthracene prefers to exist in a mixture of the patterns shown in Fig. 2(d) and Fig. 2(f), however, the pattern in Fig. 2(d) is slightly more dominant.In the case of Phenanthrene, the dominant arrangement is that of a biphenyl unit connected by a bridge, illustrated in Fig. 2(i), mixed with the slightly less preferred peripherally delocalized pattern shown in Fig. 2(g).

D. Order of stability
The final metric that we have considered is the truncated mean of the probabilities of the electrons to exist at various positions, also known as TRP.The plots of TRP for unique (i.e.inequivalent) positions for each molecule are shown in Fig. 10.Just as the case of MAXP, a lot of dynamical variation can be observed, however, a higher mean value here implies that electrons do not fully localize at that position, even for a very small time, and the position is a part of the delocalized cloud.A higher TRP over all the positions also implies the species is more stable overall as all bonds tend to have a higher degree of delocalization, and a consequently a higher resonance energy.It thus naturally gives rise to the known stability order of the four aromatic molecules considered, i.e.Benzene > Naphthalene ∼ Phenanthrene > Anthracene.

VI. CONCLUSION
We have studied the structure and properties of four benzoid polycyclic aromatic hydrocarbons using quantum walks.We have characterized the evolution of the probability of finding electrons at different points, and developed metrics with the help of which we are able to qualitatively understand the structure of these molecules.We also use some of the developed metrics to characterize the chemical properties of these molecules, namely, the characterization of the electron-rich sites in the structure, which are preferentially targeted by electrophiles in solution.We are also able to perform a relative stability analysis of these species, and our results agree with the previously established results from previous chemical studies.
This work also proposes a formalism in the analysis of aromatic compounds, wherein the quantum walk is the fundamental physical process by which electrons diffuse in the delocalized π-electron cloud.This has applications in simulations of chemical reactions via quantum simulators and/or quantum computers capable of realizing quantum walks on networks.For smaller molecules, the calculations may be done on a classical computer as well, making this formalism accessible to both classical and quantum computing paradigms.We aim to extend this analysis to non-planar molecules, as well as substituted aromatic molecules.This approach also prompts the use of machine learning techniques, such as deep neural networks, in order to calculate the bond orders by heuristic data such as the one generated by the node-ranking algorithm.

FIG. 1 .
FIG. 1.A description of the benzoid polycyclic aromatic molecules considered in this study.The labels refer to the respective carbon atoms.The figures show the following molecules -(a) Benzene, (b) Naphthalene, (c) Anthracene, and (d) Phenanthrene.

FIG. 2 .
FIG. 2. An illustration of the different possible delocalization modes for each of the benzoid polycyclic aromatic hydrocarbons considered in this work.Fig.(a) illustrates the only possible delocalization mode for Benzene, (b) and (c) depict the two possible ways by which naphthalene can delocalize.The possibilities for Anthracene are depicted in Figs.(c), (d), and (e), and those for Phenanthrene are represented in (g), (h), and (i).

FIG. 3 .
FIG. 3. The evolution of the probability distributions of the different walkers in the network representing Benzene, with the position basis superimposed on a plot of the Benzene ring for easier visualization.Each walker is represented by a different color.The evolution in case of Benzene exhibits an oscillation about the initial state.

FIG. 4 .
FIG. 4. The evolution of the probability distributions of the different walkers in the network representing Naphthalene, with the position basis superimposed on a representation of the network representing Naphthalene.Each walker is represented by a different color.

FIG. 5 .
FIG. 5.The evolution of the probability distributions of the different walkers in the network representing Anthracene, with the position basis superimposed on a representation of the network representing Anthracene.Each walker is represented by a different color.

FIG. 6 .
FIG. 6.The evolution of the probability distributions of the different walkers in the network representing Phenanthrene, with the position basis superimposed on a representation of the network representing the Phenanthrene molecule.Each walker is represented by a different color.

FIG. 7 .
FIG. 7. Results of applying the node-ranking algorithm from [28] on each of the networks representing the different molecules.The sites which are equivalent have the same ranks, and therefore the same reactivity.The plots in (a), (b), (c), (d) show the results of the algorithm when applied to networks representing Benzene, Naphthalene, Anthracene and Phenanthrene molecules, respectively.

FIG. 8 .
FIG. 8. Variation of MAXP with time at every node on each network considered in this work.The plots in (a), (b), (c), (d) depict the results obtained when applied to networks representing Benzene, Naphthalene, Anthracene and Phenanthrene molecules, respectively.A large amount of dynamical variation is observed, however, the mean values are plotted with dotted lines.The plot for Benzene shows a periodic variation, as expected.

FIG. 9 .
FIG. 9. Mean of MAXP for every node of each network considered in this work, averaged over 200 units of time, sampled at every 0.01 units.The plots in (a), (b), (c), (d) show the results obtained when applied to networks representing Benzene, Naphthalene, Anthracene and Phenanthrene molecules, respectively.As expected, the plot for Benzene shows that every site is equivalent.

FIG. 10 .
FIG. 10.TRP of the probability of an electron to exist at a particular site on each molecule, plotted for 200 units of time, sampled at every 0.01 units.The plots in (a), (b), (c), (d) depict the results obtained when applied to networks representing Benzene, Naphthalene, Anthracene and Phenanthrene molecules, respectively.The dotted line shows the mean value of the plot, averaged over time.

TABLE I .
Table showing the values of bond order obtained for various bonds in polycyclic aromatic molecules considered in this work.