Efficient Symmetrized DMRG Algorithm for Low-lying Excited States of Conjugated Carbon Systems: Application to Coronene, Ovalene and 1,12-Benzoperylene

0.23 MB·

by  Suryoday Prodhan

Efficient Symmetrized DMRG Algorithm for Low-lying Excited States of
Conjugated Carbon Systems: Application to Coronene, Ovalene and
1,12-Benzoperylene
Suryoday Prodhan,1,a) Sumit Mazumdar,2,3,b) and S. Ramasesha1,c)
1)Solid State and Structural Chemistry Unit, Indian Institute of Science, Bangalore-560012,
India.
2)Department of Physics, University of Arizona, Tucson, Arizona 85721, USA.
3)College of Optical Sciences, University of Arizona, Tucson, Arizona 85721,
USA.
(Dated: 5 January 2016)
Symmetry adapted density matrix renormalization group (SDMRG) technique has been an efficient method
for studying low-lying eigenstates in one- and quasi-one dimensional electronic systems. However, SDMRG
6
method had bottlenecks involving construction of linearly independent symmetry adapted basis states as
1
the symmetry matrices in the DMRG basis were not sparse. We have developed a modified algorithm to
0
overcome this bottleneck. The new method incorporates end-to-end interchange symmetry (C ), electron-
2 2
hole symmetry (J) and parity symmetry (P) in these calculations. The one-to-one correspondence between
n
direct-product basis states in the DMRG Hilbert space for these symmetry operations renders the symmetry
a
J matricesinthenewbasiswithmaximumsparseness,justonenon-zeromatrixelementperrow. Usingmethods
similartothoseemployedinexactdiagonalizationtechnique forPariser-Parr-Pople(PPP)models,developed
1
in the eighties, it is possible to construct orthogonal SDMRG basis states while bypassing the slow step of

] Gram-Schmidtorthonormalizationprocedure. The methodtogether withthe PPPmodelwhichincorporates
h
long-range electronic correlations is employed to study the correlated excited state spectrum of coronene,
p
substituted coronene, ovalene and 1,12-benzoperylene; all these systems can be viewed as graphene dots. It
-
m is knownthatin linearpolyenes the correlationparametersofthe modelHamiltoniandetermine the ordering
of the lowestone-photon, two-photonand triplet states, which are experimentally important. In the systems
e
h we have studied, the correlationparameters for all the orbitals of the model are same, however,the ordering
c of the important states is very different in these systems. This underlines the importance of the topology of
s. thehoppingterminthePPPmodelformolecules. Wealsodiscussindetailthe propertiesoftheseimportant
c states in the molecules we have studied. The properties include bond orders, spin density distribution in the
i
s tripletstateandtransitiondipolemomentforthelow-lyingone-photonstatesfromthegroundstate. Besides,
y the exciton binding energies and charge gaps are also reported for these systems.
h
p
[
I. INTRODUCTION the Brillouin zone leading to formation of Dirac cones.

1 FractionalquantumHalleffectatroomtemperatureand
v Successful exfoliation of a single atomic layer of tunneling throughclassicallyforbidden regionexhibiting
3 graphene from graphite crystal has led to an explosion Klein paradox are consequences of the dispersion1.
5
of research, both experimental and theoretical, in the Being a π-conjugated system, electron-electron corre-
0
area of graphene and related systems. A sheet of car- lation should play an important role in the physics of
0
0 bon atoms, arranged in honeycomb structure can be as- thistwo-dimensionalcarbonarray,althoughthestrength
. sumed to be the parent structure from which carbon al- of correlation may get diminished from well-known π-
1
lotropes like carbon nanotubes and graphite can be con- conjugated systems like polyenes, linear polyacenes or
0
6 structed. Each of the carbon atoms in graphene being carbon nanotubes due to higher dimensionality. Tight-
1 sp2-hybridized are σ-bonded with three nearest neigh- binding calculations of graphene without explicitly con-
: borswhilethe2p orbital,perpendiculartothemolecular sidering the electron-electron repulsion provides suffi-
v plain provides thze conjugation network. The completely ciently accurate ground state description2, yet, ignoring
i
X filledσ-bandliesdeepdowninenergyandprimarilyren- explicitelectroncorrelationshasbeen criticizedrecently.
r ders stability to the structure. On the other hand, the Onthecontrary,organicpolyeneandpolyacenesystems,
a half filled π-band, due to the overlapof the 2p orbitals, which are quasi-one dimensional, have been well studied
z
ismostlyresponsibleforfascinatingpropertieslikelinear in the long-ranged Pariser-Parr-Pople(PPP) model and
dispersion of electronic energy band near the vertices of the effect of electron-electron correlation is found to be
qualitativelysignificant. Theorderingofexcitedstatesin
polyenesdifferfromsingleCIcalculationsanditprovides
anexplanationofthe opticalspectrainthesesystems3,4.
a)Electronicmail: [email protected] Correlation also plays a significant role in dynamical
b)Electronicmail: [email protected] properties of polyenes; observation of electrolumines-
c)Electronicmail: [email protected] cenceefficiencywhichis farinexcessoftheupper bound
2
placedbythespinstatistics,getsresolvedonlyaftertak- and highly accurate DMRG algorithm for utilizing end-
ing explicit electron correlations in account5. Crossover to-endinterchangesymmetry(C ),electron-holesymme-
2
of the dipole allowed or one-photon state (1B ) and the try (J) and spin-flip or parity symmetry (P). The ad-
u
dark two-photonstate (2A ) in polyacene series with in- vantage of this algorithm overthe previous symmetrized
g
creasing chain length, is observedin the PPP model cal- DMRG method is that the symmetry operators can be
culations and agrees well with experiments6. Correlated expressedashighlysparsematriceswithonlyonenonzero
electronic calculations even support the higher optical element per row. This allows representing the symme-
gaps observed in carbon nanotubes with different diam- try operators by vectors of the same dimension as of the
eters, contrary to tight binding calculations7. The ef- DMRG Hilbert space. Hence it is possible to employ
fect of an additional dimension on the correlation effect, very large DMRG cut-offs in the new symmetry adap-
therefore, is worthwhile investigating starting from sys- tation algorithm. We have used this new algorithm to
tems suchaspolycyclic peri-condensedacenes(graphene calculate the ground state, lowest-lying two-photon and
nano-discs), which resemble closely to two-dimensional one-photonstatesinM =0sectoralongwithafewlow-
S
graphene system. lying triplet states in M = 1 sector for coronene, ova-
S
lene, 1,12-benzoperylene and substituted coronene. No-
While a number of optical studies have been car-
riedoutonpolynucleararomatichydrocarbons8,mostof tably,coronenehasD6h symmetry,while the otherthree
are of lowersymmetry groups;nonetheless, in our study,
the theoretical studies are either in the non-interacting
we have considered only one C symmetry axis in all
limitorapproximatingthe electronrepulsionviaa mean 2
the molecules. We have also characterized these states
field. Moffitt and Coulson calculated the bond order in
by calculating bond orders along with the spin density
coronene exploiting both Hu¨ckel’s tight binding method
at different sites (in M = 1 sector). Charge gaps and
and valence bond method without considering high en- S
ergy ionic structures9 while Pritchard and Sumner cal- exciton binding energies for coronene, ovalene and 1,12-
benzoperylene are also mentioned here.
culated the ground state wavefunction and bond order
The paper is organizedas follows. In the next section,
for coronene and ovalene using self-consistent molecu-
lar orbital theory10. Hummel and Ruedenberg also used we have given an account of the modified algorithm and
compared it with exact calculations for 16-site polyene
similar techniques to study other large polyaromatic hy-
drocarbons. A number of electronic structure calcula- chain and tetracene with 18 pz-orbitals in conjugation.
InsectionIII,wehavediscussedourresultsforcoronene,
tions of these polycyclic compounds and their ions, em-
ovalene, 1,12-benzoperylene and substituted coronene.
ploying density functional theory (DFT) and its time-
dependent variant (TDDFT) have been reported11. The In the last section we have summarized our study and
concluded.
accuracy of these calculations depend upon the reliabil-
ity of the functional form of the exchange correlation
potential which still is not beyond doubt for strongly
correlated systems. Furthermore, implementation of the II. METHODOLOGY
TDDFT method considers only single particle-hole exci-
tations which lacks the desired accuracy. In recent past, A. Modified Algorithm For Symmetry Adaptation
excited state ordering in some of these molecules were
studiedusingmultiplereferenceconfigurationinteraction The DMRG procedureis a truncatedHilbertspace al-
(MRCI) approach,which considers upto quadruple exci- gorithmwhichadoptsiterativeblock-buildingschemefor
tations of each targeted state12. The theoretical predic- calculatinglow-lyingeigenstatesoflargeoneorquasi-one
tions are in agreement with the optical and two-photon dimensional Hamiltonians14,16. The system block or left
absorption data. Another study has speculated the ef- block(L)alongwiththeenvironmentblock,alsodenoted
fect of different molecular geometries over the ordering asrightblock (R)arelinearlygrownthroughadditionof
of excited states in polynuclear hydrocarbons13. While one new site to each block per iteration and these two
theHartree-Focksemiempericalcalculationsareoversim- together form the superblock which corresponds to the
plified, the accuracy of the other method used is still physical system to be studied. In the infinite DMRG al-
dubious. This problem can be overcome by exact di- gorithm, superblock of size 2l consists of system and en-
agonalizationtechnique in full Hilbert space, though the vironment blocks of same size (l). The reduced density
sizeoftheHilbertspaceincreasesexponentiallywithsys- matrix of the system (environment) is generated from
temsizeandmakesitcomputationallyexpensiveevenfor low-lying eigenstates of the superblock by tracing over
moderate-sized systems. the environment (system) states. Matrix elements of
The density matrix renormalization group method the reduced density matrix of the system (environment)
(DMRG), put forward by White14, is a well-accepted are represented by ρm,m′ =Pnψm,nψm∗′,n, where m,m′
methodfor calculatinglow-lyingstatesofone andquasi- aresystem(environment)stateswhilenareenvironment
one dimensional systems. Ramasesha et al. have suc- (system) states. Reduced density matrices for both the
cessfullyimplemented spatialanddiscretesymmetriesin system and environment are diagonalized and Ml den-
DMRG to target low-lying states in specific symmetry sity matrix eigenstates (DME) with highest eigenvalues
subspaces6,15. In this study, we have employed a new ({µ }for system and {µ }for environment)arestoredas
L R
3
columnvectorsofaMl−1d ×Ml matrix,whereMl−1 is The Hamiltonian of a half-filled bipartite system con-
σ
the number of truncated reduced density matrix eigen- serves electron-hole symmetry (J), where interchanging
states retained at (l −1)-th iteration and d is the di- the creation and annihilation operators at each lattice
σ
mension of the Fock space of the newly added site in sites with an associated phase factor keeps the Hamilto-
l-th iteration. The Hamiltonian and site operators of nianinvariant. Onapplicationofthissymmetryelement,
the newly generated blocks are constructed in the basis Fock states of a single site transforms in the following
µl−1σandthenrenormalizedusingtheMl−1d ×Ml ma- way:
σ
trix. Further addition of two new sites gives a system of
2l+2sites,whosebasis|µlLσiσi′µlRiisgeneratedbydirect Jˆi|0i=|↑↓i Jˆi|↑i=(−1)ηi|↑i
productof |µlLi ofthe systemblock, |µlRi ofthe environ- Jˆi|↑↓i=−|0i Jˆi|↓i=(−1)ηi|↓i
mentblockwith|σiiand|σi′iofthenewlyaddedsitesto
the left and right blocks respectively. The Hamiltonian where η is 1 if ‘i’ belongs to A-sublattice and 0 if it
i
for the 2l+2 sites system is formed in this direct prod- belongs to B-sublattice.
uctspaceanddiagonalizedtogetthe energyeigenstates, The result of parity or spin-flip symmetry (P) which
which are utilized for the next iteration. In case of sys- flips the spin at each site is given by the following:
temswithfinitesizeN,finiteDMRGalgorithmenhances
the accuracy in eigenstate calculation; infinite algorithm Pˆ|0i=|0i Pˆ|↑i=|↓i
i i
is employed starting from a small system until the de-
Pˆ|↑↓i=−|↑↓i Pˆ|↓i=|↑i
sired superblock size is reached and then sweeping over i i
thesuperblockisappliedbychangingthesystemanden-
Non-relativistic Hamiltonians remain invariant under
vironment sizes, keeping the superblock size fixed. The
parity symmetry. This symmetry can only be employed
basis at a step of the finite DMRG sweeping can be rep-
in M = 0 space, as application of this symmetry maps
resented by |µlLσiσi′µNR−l+2i, where l can vary between fromSone non-zero MS space to a −MS space.
1 and N −3. Once the size of the system (environment)
BothJ andP operatorsaresiteoperatorsandhence,
i i
block reaches its maximum, the environment (system)
the operators for the full system can be defined as:
blockstartsto growattheexpense ofthe otherblock. A
complete sweep consists of expansion of both blocks to Jˆ=YJˆi Pˆ =YPˆi
their maximum and returning back to the block size of
i i
the final infinite DMRG step.
The symmetryelementsconsideredinourmodifiedal- Thematrixelementofthesetwosymmetryoperatorsfor
gorithm are same as in the work by Ramasesha et al.15. a superblock of size 2l+2 can be represented by:
The C symmetry about an axis perpendicular to the
2 hµ′σ′σµ|Rˆ |νττ′ν′i
plane of the molecule corresponds to the end-to-end in- 2l+2
(1)
terchange of the superblock with a phase factor: =hµ|Rˆl|νihσ|Rˆi|τihσ′|Rˆi′|τ′ihµ′|Rˆl|ν′i
Cˆ2|µlLσiσi′µlRi=(−1)γ|µlRσi′σiµlLi i.e. itisdirectproductofmatrixelementsofcorrespond-
or in shorter notation (we will use onwards) ing symmetry operators in the left block DME space,
rightblockDMEspaceandintheFockspacesofthetwo
Cˆ2|µσσ′µ′i=(−1)γ|µ′σ′σµi newly added sites.
γ =nµ(nσ+nσ′ +nµ′)+nσ(nσ′ +nµ′)+nσ′nµ′ The symmetry operators of the superblock follow the
relations
i.e. the system and environment blocks inter-
change their density matrix eigenstates while the two (Cˆ )2 =1 Jˆ2 =(−1)(Ns+Ne) Pˆ2 =1 (2)
2
newly added sites also interchange their Fock states.
nµ,nσ,nσ′ and nµ′ are the occupancies in the system where Ns and Ne are number of sites and number of
block, environment block and in the two newly added electrons in a basis state, onwhich the symmetry opera-
sites. Astheinterchangeshouldnotaffectthesuperblock tions are applied. In case of half-filled fermionic system
Hilbert space, system and environment blocks should N = N ; subsequently, each of the operator becomes
s e
have same density matrix eigenspaces. In other words, their own inverse, the symmetry operators commute
employing C symmetry approximately halves the effec- among themselves and form an abelian group. However,
2
tive size of the superblock Hilbert space beside reduc- the symmetry operators may not commute in individual
ing the computational expense by avoiding environment block spaces and in Fock spaces of newly added sites.
block calculations, at every infinite DMRG step and at Asmentionedearlier,the C symmetryoperatormaps
2
the last step of the finite DMRG iteration. In this sym- onedirectproductstateexclusivelytoanotherone. How-
metry operation, each DMRG basis state is carried over ever, J and P operators do not produce a single state,
to only one unique DMRG basis state, with a phase of due to the absence of one-to-one symmetry correspon-
±1. Thus the C symmetry operationcanbe storedas a dence between different density matrix eigenstates of a
2
column vector with entry ±j , resulting from C opera- particular block, although the mapping of Fock states
k 2
tion on the state k. of a single site is exclusive. This absence of one-to-one
4
mappingofthe statesinthe sub-blocksandhence aone- where ηs are the corresponding phase factors.
to-one mapping of the states of the superblock leads to (ii){µi(ns >ne,MS =0)}: ApplicationofJL symme-
some serious problems in symmetry adaptation. Firstly, try operator on these density matrix eigenvectors gener-
the matrix representation of the operators in the DME ate eigenvectors with particle number 2n −n ; yet, P
s e L
basis is not very sparse. Secondly, it is difficult to keep symmetryshouldmaptheeigenspacetoitself. Therefore,
track of and thereby include all the symmetry partners to begin with, this group of eigenvectors is further sub-
of the J and P symmetry while introducing the DMRG divided into smaller sub-groups {µ (n = n′,M = 0)},
i e e S
cut-off. Thirdly, weeding out linearly dependent states where n′ is some particular value of n and n′ < n .
e e e s
after symmetry adaptation is a slow step as it involves Considering these smaller blocks, one at a time, we set
Gram-Schmidt orthonormalization of the symmetrized up the matrix of Pˆ in this sub-space. This is diagonal-
superblock states. Thus, a successful symmetry adap- ized to obtain the eigenvectors of Pˆ so that each of the
tation scheme must have a one-to-one unique correspon- new DMEs map into themselves with a phase factor of
dence between states in the DME basis space. ±1. These transformed µ (n =n′,M =0) provide a
k e e S
We invoke the one-to-one correspondence among dif- one-to-one mapping under the symmetry operations Jˆ,
ferent DMEs which results in only one non-zero ele- Pˆ and PˆJˆ.
ment in each row of the matrices of symmetry opera- (iii) {µi(ns =ne,MS <0)}: This set of eigenstates
tors. This is achieved by considering truncated reduced are similar to the one discussed above, differing only in
density matrix eigenstates restricted to a specific region the fact that J will map the space to itself instead of
L
of the Hilbert space, as by symmetry we can obtain the P . The space is sub-divided into smaller blocks on the
L
DMEs in the complementary spaces. For example, the basis of different negative M values and the matrices
S
DMEs with total number of electrons ne <ns, the num- of the operator Jˆ are setup in these smaller DME ba-
berofsitesinthesystemblock,arerelatedtoDMEswith sis spaces. The eigenstates of this operator serve as the
ne >ns by – new DMEs for obtaining one-to-one correspondence of
the DMEs under the symmetry operations.
Jˆ µ (n >n )=µ (n <n ) (3)
L i s e j s e (iv) {µi(ns =ne,MS =0)}: This particular set of
and similarly the DMEs with M > 0 are related to eigenvectors span a vector space which on application
S
DMEs with MS <0 by, by both JL and PL should map into itself. Hence, lin-
ear combinations are found out which are eigenvectors
PˆLµi(MS <0)=µj(MS >0). (4) of JL with eigenvalues (−1)ζJ, ζJ = 0 or 1. Using J-
eigenvectorswithaparticulareigenvalue(−1)ζJ,another
Hence,westartwiththefollowingsetofDMEofthesys-
setoflinearcombinationsareformedwhichareeigenvec-
tem block which satisfy the criteria {µ (n ≥ n ,M ≤
i s e S torsofP , witheigenvalues(−1)ζP, ζ =0or1. Tothis
0)}. L P
end, we first construct the matrix representation of Jˆin
We consider here only the left block, as the algorithm
thebasisofDMEswithn =n andM =0. We obtain
will be same for the right block as well. This particular s e S
the eigenvectors of Jˆin this space by diagonalizing the
set of eigenstates can be segregated into four groups –
(i){µi(ns >ne,MS <0)}: ApplicationofJL, PL and mwhaitcrhixt.hWe emtahternixuosfePˆthiesesiegtenuvpe.ctTohresorfesJˆulatsintgheeibgeansivseicn-
theirproductontheseeigenstatesgivethecorresponding
tors of Pˆ replace the DMEs in this space and they will
symmetry counterparts. As J and P may not com-
L L
be employed in symmetry adaptation.
mute, we only considered the product P J ; the order
L L
While imposing the cut-off in the DMEs, we only con-
of operationsof the two operatorswill become irrelevant
sider the subspace of DMEs with n ≤ n and M ≤ 0.
in the direct product basis space of a half-filled system, e s S
Since the DMEs related to these states by the operation
where J and P of the superblock commute.
of Jˆ, Pˆ and PˆJˆwill also have the same eigenvalues, it is
Jˆ µ (n > n ,M < 0) = µ (n < n ,M < 0) ensuredthatallDMEsaboveacut-offinthedensityma-
L i s e S iJ s e S
trixeigenvalueareincluded. Usuallythe cut-offimposed
Pˆ µ (n > n ,M < 0) = µ (n > n ,M > 0)
L i s e S iP s e S in the specified sector is about 0.4 times of the cut-off
Pˆ Jˆ µ (n >n ,M <0)=µ (n <n ,M >0) intended for the full density matrix. Since we have a
L L i s e S iPJ s e S
one-to-one correspondence of the DMEs within a block
The DMEs obtained from µi by symmetry operations for the operations Jˆ, Pˆ and PˆJˆ, it is a simple matter to
are indexed and stored as eigenvectors {µi}. Symmetry extendtheone-to-onecorrespondencetothedirectprod-
correspondence in every pair of these four eigenstates is uctbasisstatesformedbytheDMEsoftheleftandright
alsodeducedusingtherelationPˆLJˆL =(−1)(ns+ne)JˆLPˆL blocksandtheFockspacestatesofthenewlyaddedsites.
along with Eqn. 2 and stored separately. The matrix Symmetry adapted direct product bases in an irre-
elements of the symmetry operators are given by, ducible representation Γ can be generated by employing
projection operator, which is given by –
(Jˆ ) =(−1)ηJδ (Pˆ ) =(−1)ηPδ
L i,j i,iJ L i,j i,iP
1
(Pˆ Jˆ ) =(−1)ηPJδ P(Γ)= hXχ(Γ,Ri)Rˆi (5)
L L i,j i,iPJ i
5
whereχ(Γ,R )isthecharacterofthesymmetryoperator action:
i
Rˆ intheirreduciblerepresentationΓandhbeingtheor-
deirofthegroup. OperationbyP(Γ)oneachofthebasis H = X t0(cˆ†i,σcˆj,σ+ H.C.)+Xǫinˆi
statesinthedirectproductspaceproducesthesymmetry <i,j>,σ i
(7)
adapted linear combinations. However, if P(Γ) acts on U
+X nˆi(nˆi−1)+XVij(nˆi−zi)(nˆj −zj)
eachofthestatesinthedirectproductspace,wewillend 2
i i>j
up with a linearly dependent symmetry adapted basis.
In the earlier symmetry adaptation procedure15, weed- t is the nearest-neighbor hopping integral between
0
ing out linear dependence by Gram-Schmidt orthonor- bonded sites ‘i’ and ‘j’, ǫ is the site energy of ‘i’-th car-
i
malization technique was an extremely expensive step, bonatomwhileU istheon-siteCoulombrepulsionterm.
computationally. Withpresenttechnique,wherewehave Theintersiteinteractionenergies,V areinterpolatedby
ij
a one-to-one correspondence between symmetry related Ohno scheme18 assuming C-C bond length of 1.4 ˚A. cˆ†
i,σ
states,thisturnsouttobeextremelysimple andcompu- and cˆ are the creation and annihilation operators of
i,σ
tationally veryefficient. Firstly we note that, one-to-one an electron with spin σ in the 2p orbital of the carbon
z
correspondence between the DMRG basis states under atomatsite ‘i’andnˆ = cˆ† cˆ isthe corresponding
the symmetry operationsimplies thata givenbasis state i Pσ i,σ i,σ
number operator. For carbon atoms, the standard PPP
can appear only once (or not at all) in a given symme-
parametersaret =−2.4eVandU =11.26eVwhilethe
0
try space. Therefore, in the symmetry adaptation step,
site energy(ǫ)for unsubstitutedcarbonatomistakento
we sequentially go through the list of basis states, apply
be zero. z is the local chemical potential, expressed by
i
the projection operator and construct the symmetrized
the the number of electrons at site ‘i’ which leaves the
combinationalongwitheliminationofthepartnerswhich
site neutral; for carbon atoms in a π-conjugated system,
haveappearedinthesymmetrycombinationfromthelist
z = 1. Substitution effect by donor or acceptor groups
i
of unsymmetrized basis states. This ensures that a ba-
canbe modeledbyconsideringthesiteenergyǫ positive
i
sis state occurs at most once in the symmetrized states.
or negative respectively.
Taken together with the orthonormality of the unsym-
As DMRG calculations are done in constant M -basis
S
metrized basis states, this ensures orthogonality of the
space,contrarytothespin-adaptedbasis(Diagrammatic
basis states generated by the projection operator. Us-
valencebondbasis)calculations,labellingstatesbytheir
ing these symmetry combinations, we can generate the
spin value is ambiguous. The parity symmetry partially
transformationmatrix S, which is D×D matrix where
Γ resolves this as the singlet state lies in the even-parity
D is the dimensionality of the unsymmetrized space and
subspace while the triplet state lies in the odd-parity
D is the dimensionality of the symmetrized space.
Γ subspace. Though, higher spin states also lie in the
TheHamiltonianofthesuperblockcanbetransformed
even (odd) parity sub-spaces, they remain high in en-
into symmetrized basis space by operating with the ma-
ergy compared to singlet and triplet states in moderate
trix S as: sized carbon-based systems and so can be neglected6.
H˜ =S†H S (6)
2l+2 2l+2
C. Comparison with exact calculation results
ThetransformedHamiltonianisthendiagonalizedtoget
lowest-lying eigenvectors in different irreducible repre-
The accuracy of our modified algorithm can be estab-
sentations labeled as eA+, eB+, eA−, eB−, oA+, oB+,
lished by comparing the results obtained for small sys-
oA− andoB−,wheretheleft-handsuperscriptrepresents
temswithexactPPPresults,likein16-sitepolyenechain
characterunderP-symmetry,the rightsuperscriptisthe
andtetracenemoleculewhichhas18sites. Thedifference
character under J-symmetry and the letter A/B is the
betweengroundstateenergyof16-sitepolyenechain,ob-
characterunderC symmetryintheirreduciblerepresen-
2 tained by symmetrized DMRG technique, and the exact
tationΓ. Thesingletliesinthe‘e’spacewhilethetriplet
calculation result is plotted for four different cutoff val-
lies in the ‘o’ space. The space ‘+’ refers to ‘covalent’
ues M′l in the sector {µL(n ≥ n ,M ≤ 0)} of the
space while ‘-’ to the ‘ionic’ space. A(B) represents the i s e S
sub-block in Fig. 1. In all the calculations, the energy
space even (odd) under C symmetry. The ground state
2 values are accurate at least up to fifth decimal places
is usually in eA+ while the optically connected space is
when compared with exact calculation results. We have
eB−. The lowest triplet lies in the oB+ space.
also plotted the fluctuations in the DMRG energy about
the exact energy with number of finite DMRG sweeps
for a fixed M′l value, which depicts that the fluctuation
B. Model Hamiltonian
is negligible with number of sweeps and therefore, 2−3
sweeps will suffice for acceptable convergence of energy.
ThemodelHamiltonianemployedforthestudyofperi- In Table I, PPP energies of the ground state (11A+),
condensed polyacenes is the well known Pariser-Parr- lowesttwo-photonstate(21A+), lowestone-photonstate
Pople (PPP) Hamiltonian17, which considers long-range (11B−)andofthelowesttripletstate(13B+)ofboth16-
Coulombic interactionalongwith on-site Hubbardinter- site polyene chain and tetracene molecule are presented.
6
1.5e-06 {+ε}1 17{−ε} 3 2 x 21 20
∆E1.29ee--0067 ∆Esweep33..311.241eee---000777 541 2236311462157118051691134 112892 2201 x 46 54675 1832111321151110114 123212 181179 y
6e-07 1 2nswe3ep 4 {7−ε}8 9 24{+ε2}3 y 7 8 99 1014 15 16
3 1 23
3e-07 y 4 4 5 2 6 22 24
3 10
6 25 7 14915 1116 21 20 25 26
0 150 225 M’l 300 375 71 813 1217 1819 27 x
8 12 28
9 31 29
FIG. 1. (Color online) Deviation of ground state energy cal- 11
10 32 30
culated by symmetrized DMRGtechniquefrom exact energy
valuefor differentcutoffM′l ofthedensitymatrix eigenstate FIG.2. (Coloronline)Schematicdiagramsofcoronene,1,12-
space. The numberof finite DMRG sweeps (2) is fixed in all benzoperylene and ovalene with site (black and bold) and
the cases. Inset: Fluctuation of DMRG ground state energy bond (red and italics) indices. The sites of substitution in
with respecttoexactenergy valuewithincreasing numberof substitutedcoroneneareindicatedincurlybrackets;+ǫ(sites
finite DMRG sweeps. The number of density matrix eigen- 1 and 24) represents a donor site while −ǫ (sites 8 and 17)
states retained, M′l, is 225. represents an acceptor site.
III. RESULTS AND DISCUSSION
M′l is kept at 225 for 16-site polyene chain and 425 for
tetracene molecule, thus total number of DMEs retained We have studied a few low-lying states of coronene,
for the sub-block is ∼ 550 and ∼ 1000 respectively at ovalene and 1,12−benzoperylene in different symmetry
each step of DMRG calculation. Full CI calculations, spaces considering ∼1000 DMEs in PPP model and the
employingdiagrammaticvalencebondmethod(DVB)19, energy gaps are tabulated in Table. II and III. We have
shows excellent agreementwith the results obtained and alsocalculatedenergiesintheabovesymmetryspacesfor
the upper bound of errorin energies ofdifferent states is substituted coronene,with donorandacceptorgroupsof
of the order of 10−4. equal strength, set to 1.0 eV (|ǫ| = 1.0 eV). We have
employed only a few of the existing spatial symmetries
of the molecules, namely C symmetry axis perpendic-
2
ular to the molecular plane (in plane C axis in case of
2
1,12−benzoperylene)forourcalculations,asourDMRG
calculations cannot incorporate more than one symme-
TABLEI.Comparison oftheenergiesofgroundstate,lowest
two-photonstate,lowestopticalstateandlowest tripletstate try axis. Being treated as planar molecules, with Wan-
of 16-site polyene chain and tetracene, obtained by modified nier orbitals, the reflection in the molecular plane is a
symmetrized DMRG algorithm and exact calculation. The trivial symmetry. Yet, we can label the eigenstates ac-
exactcalculationsarecarriedoutusingdiagrammaticvalence cordingtodifferentirreduciblerepresentationsofthe full
bond method. group based on transition dipole moment direction for
excitations from the ground state. In case of coronene,
System State Exact calculation DMRGcalculation
wehaveconsidereditssubgroupD ratherthanthe full
2h
Polyeneof 11A+ −33.55575 −33.55575 group D for labelling the eigenstates. The lowest op-
6h
16 sites 21A+ −32.03385 −32.03385 tical states obtained may not be the states with high-
est transition dipole moment, which will be prominent
11B− −30.43291 −30.43291
inUV-visiblespectroscopybutcorrespondstothelowest
13B+ −32.87993 −32.87993 energystateoftheappropriatesymmetrysubspace;how-
Tetracene 11A+ −43.69083 −43.69076 ever,theresultsobtainedareingoodagreementwithgas
phase molecular spectroscopic data available for above
molecule 21A+ −40.55749 −40.55725 molecules8, suggesting that the states obtained are also
11B− −40.51439 −40.51419 the ones with considerable transition dipoles in the en-
ergy range of the experiment.
13B+ −42.47046 −42.47038
The relative ordering of the energy levels are closely
associated with the effective correlation strength in lin-
7
TABLE II. Lowest energy states in coronene and ovalene with the energy gaps in units of eV. Although coronene has D
6h
symmetry, here the energy states are labelled according to Abelian D symmetry representations. The transition dipole
2h
moment from ground state totheexcited states (µ ) along specific axes are also specified in thelast two columns in units
tr,x/y
of Debye.
Coronene
Natureof the state State level Energy Gap µ (D) µ (D)
tr,x tr,y
Two-photon (21A+/11B+)a 3.97 0.00 0.00
g 1g
(31A+/21B+) 4.09 0.00 0.00
g 1g
(41A+/31B+) 5.08 0.00 0.00
g 1g
Optical 11B− 4.83 0.81b 8.43
2u
11B− 4.87 7.51 0.49
3u
Triplet (13B+/13B+) 2.35 0.00 0.00
2u 3u
(23B+/23B+) 3.00 0.00 0.00
2u 3u
(33B+/33B+) 3.02 0.00 0.00
2u 3u
(13A+/13B+) 3.35 0.00 0.00
g 1g
Dark states (11B+/11B+) 2.82 0.00 0.00
2u 3u
(x1B−/ x1B−)c 3.88 0.28 1.74
2u 3u
(11A−/11B−) 4.73 0.00 0.00
g 1g
Ovalene
Two-photon (21A+/11B+) 3.10 0.00 0.00
g 1g
(31A+/21B+) 3.75 0.00 0.00
g 1g
(41A+/31B+) 5.48 0.00 0.00
g 1g
Optical 11B− 3.32 0.29 6.58
2u
21B− 4.53 0.33 0.58
2u
11B− 4.64 7.28 0.22
3u
Triplet (13B+/13B+) 1.49 0.00 0.00
2u 3u
(13A+/13B+) 2.77 0.00 0.00
g 1g
(23A+/23B+) 2.80 0.00 0.00
g 1g
(23B+/23B+) 2.82 0.00 0.00
2u 3u
Dark states (11B+/11B+) 2.65 0.00 0.00
2u 3u
(11A−/11B−) 4.06 0.00 0.00
g 1g
(21A−/21B−) 4.69 0.00 0.00
g 1g
a Probablelabelsforstates whichcannot beidentifieduniquelyaregivenwithinparentheses.
b Non-zerovalueoftransitiondipolemomentalongpolarizationdirectionforbiddenbysymmetryisanartifactofthecalculations as
averagedensitymatrices,calculated fromeigenstates ofdifferentsymmetrysubspaces, areemployedtodeterminethetransitiondipole
moment. However,theerrorsarenegligibleasintensitiesdependonthesquareofthetransitiondipolemoment.
c ximpliesthatthelevelnumberingisunknown.
ear π−conjugated systems. In extended systems effec- based on one-particle bandwidth which are similar in all
tive correlation strength can be thought of as the ratio themoleculeswehavestudied. InPPPorextendedHub-
of the on-site PPP correlation to the bandwidth of the bardmodel,theon-sitecorrelationstrengthtoafirstap-
one-particle spectrum. However, in molecules with dis- proximation is expressed by U ∼ (U −V ), where
PPP 12
crete excitation spectrum, it is difficult to use the width V is the nearest-neighbor PPP term. Our calculations
12
oftheenergyspectruminplaceofthebandwidth. Thus, yield the lowest two-photon gap in ovalene molecule to
it is difficult to a prioripredictthe orderingofthe states be almost twice that of the lowest spin gap, which is a
8
TABLE III. Lowest energy states in 1,12-benzoperylene and substituted coronene with the energy gaps in units of eV. The
transition dipole moment from ground state to theexcited states (µ ) along specific axesare also specified in thelast two
tr,x/y
columns in unitsof Debye.
1,12-Benzoperylene
Natureof the state Statelevel Energy Gap µ (D) µ (D)
tr,x tr,y
Two-photon 21A+ 2.99 0.00 0.00
1
31A+ 4.26 0.00 0.00
1
41A+ 4.58 0.00 0.00
1
Optical 11B− 3.72 0.00 3.17
1
21B− 4.71 0.00 3.83
1
31B− 4.83 0.00 5.20
1
Triplet 13B+ 2.84 0.00 0.00
1
23B+ 3.65 0.00 0.00
1
33B+ 4.26 0.00 0.00
1
Substitutedcoronene
Two-photon 21A 4.01 0.00 0.00
g
31A 4.81 0.00 0.00
g
41A 4.90 0.00 0.00
g
Optical 11B 2.90 0.00 0.23
u
21B 3.87 0.12 1.46
u
31B 5.18 0.17 7.52
u
41B 5.97 4.10 0.59
u
Triplet 13B 2.29 0.00 0.00
u
23B 2.92 0.00 0.00
u
33B 3.48 0.00 0.00
u
characteristic feature of regular polyene systems where al.13. The two-photon gaps reported earlier are 3.96 eV
two-photon state is predominantly two triplet excita- and 3.03 eV respectively for these systems, against 3.97
tions. Theothertwomoleculesandsubstitutedcoronene, eV and 3.10 eV, found in our calculations. The triplet
on the other hand, do not show this signature. There- gapsearlierreportedare2.38eVand1.57eVforcoronene
fore, it can be suggested that the effective electronic and ovalene respectively, which are close to 2.35 eV and
correlation strength in ovalene molecule is comparable 1.49eV obtainedby us. However,the lowestone-photon
to regular polyene systems while the π− electrons in gaps and relative positions of the lowest optical states
other molecules are less correlated. The ordering of low- withrespecttothelowesttwo-photonstatesdonotagree
lying singlet excited states in coronene and ovalene, al- wellwiththe previousstudy. Incoronene,the lowestop-
though, are in accordance with those in strongly cor- tical gap reported earlier is 4.11 eV while we have found
related systems; the lowest optical states lie above the two nearly degenerate excitations at higher energies of
lowest two-photon states in both the systems. However, values, 4.83 and 4.87 eV. A low-lying state at an energy
1,12-benzoperyleneshowssomepeculiartrendasthelow- of3.88eV,whichhasasmalltransitiondipoleappearsin
est two-photon state remains close to the lowest triplet thesamesubspaceasoftheopticalstates. Thisisproba-
state and quite low in energy compared to 11B−, show- blybecausewehaveincorporatedonlyoneC symmetry
1 2
ing that in this system, we cannot view the 21A+ state axisinourcalculations. InDMRGcalculations,forcom-
1
as made up oftwo triplets generatedindifferent partsof puting transition dipole moments we have used average
the molecule. density matrices, obtained from eigenstates of different
The lowest two-photon gaps along with the lowest symmetry subspaces and hence the spatial symmetry is
tripletgapscalculatedincoroneneandovaleneareinex- notretained. Thiscouldleadtoasmalltransitiondipole
cellentagreementwiththeearlierreportbyAryanpouret moment as an artifact for this state (3.88 eV) although
9
0.24
TseActBoLrEanIVd.inLoMweSst=tri0plseetcsttoarteweitnherχgiPes=obt−a1inefodrincoMroSne=ne1, n n>αβii00..1281 Coronene
ovalene and 1,12-benzoperylene in units of eV. <
0.15
0 1 2 3 4 5 6 7 8 9 10 11 12
M =1 M =0,χ =−1 0.24
S S P >
β 0.21
Coronene −58.67378 −58.56182 nαi0.18 GTwroou-pnhdo sttoante state
Ovalene −80.50885 −80.37843 <ni0.15 Ovalene OTrpitpicleatl ssttaattee
0.12
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1,12-Benzoperylene −52.54237 −52.32442 0.24
>
β 0.21
ni
α 0.18
<ni0.15 1,12-Benzoperylene
0.12
0 1 2 3 4 5 6 7 8 9 10 11
Site index (i)
this state correspond e−h allowed subspace and if spa-
tialsymmetryisnotretainedstrictlywewillhaveaweak
FIG.3. (Coloronline)Doubleoccupancyplotwithsiteindex
transition dipole. On the other hand, the 2.82 eV exci- for coronene (top panel), ovalene (middle panel) and 1,12-
tation in coronene is optically forbidden as it belongs to benzoperylene(bottom panel). Asall molecules aresymmet-
the same e-h symmetry subspace as the ground state. ricaboutC axis,onlyhalfofthemolecularsitesareplotted,
2
This argument is supported by the fact that this state the other half being symmetry equivalent to it. The color
acquires some intensity on breaking the e− h symme- schemefordoubleoccupancyindifferenteigenstatesaresame
try by introducing substituents, as can be seen from Ta- inallthreepanelsandareasfollows: groundstate(blackfilled
ble. III. Thus it appears that as in benzene, in coronene circles),lowesttwo-photonstate(redfilledsquare),lowestop-
ticalstate(greenfilledtriangles)andlowesttripletstate(blue
the lowest degenerate excitations are pushed up in en-
open squares). Doubleoccupancy in theother partnerof the
ergy while the higher energy non-degenerate states cross
degenerate optical states in coronene is shown by green bro-
these excitations on incorporating electron-electron in-
kenlinewithopentriangles. Linesconnectingdatapointsare
teractions. In case of ovalene, we have found the lowest
given only as a guide to theeye.
opticalgaptobegreaterthanthelowesttwo-photongap
while reverse ordering of states has been reported ear-
lier. The earlier reported energy gap between the lowest
one-andtwo-photonstatesaswellasthatobtainedbyus Hilbert space becomes half of that in M =1 sector cal-
S
are relatively small. The earlier study used interaction culation. As DMRG is a variational technique, accuracy
parameters optimized to yield experimental gaps in con- increases with the dimension of the Hilbert space and
densed phase while our study employsparameterswhich hence,theM =1sectorcalculationsaremoreaccurate.
S
are applicable to isolated molecules. In condensed phase Spin gap is also a good measure of electronic correla-
the one-photon state will shift more towards the longer tionstrength. ForaPPPchain,thespingapcanbeesti-
wavelengthregimethantwo-photonstate,theformerbe- matedfromtheHeisenbergmodel,withanantiferromag-
ing ionic in nature. Thus we expect that the one-photon netic J ∼ 4t2/U . This implies that with increasing
PPP
state will be below the two-photon state for ovalene in effectivecorrelationstrength,thegapbetweenthelowest
solution. The polarization in ovalene in the lowest one- singletstateandthelowesttripletstatereduces. Accord-
photon state is along the short axis while the higher en- ing to this line of argument, ovalene has the highest ef-
ergyopticalexcitationis longaxispolarized. This is due fectivecorrelationstrengthwhile1,12-Benzoperylenehas
to the fact that excitons with charge separation along thelowestamongthemoleculessitedinTable. IIandIII.
the short axis are more stable than those with separa- Substitution by donor-acceptor groups does not seem
tion along the long axis. Indeed, similar polarization for to have appreciable effect on correlation strength, al-
dipole excitations is also found in polyacenes6. though, the lifting of electron-hole symmetry allows op-
Thelowesttripletstateenergyineachofthemolecules tical transition to states which are labelled as dipole-
can be calculated in two ways while employing both C forbidden in unsubstituted molecules. The extent of
2
and J symmetries – either by calculating the lowest en- mixing ofdifferent symmetry states ofthe unsubstituted
ergy state in M =1 sector or by calculating the lowest system due to substitution depends on the strength of
S
energy state in the M = 0 space with χ = −1. The the donor-acceptor groups. In the present case, the
S P
energies obtained by both calculations are given in Ta- dipole-forbidden 11B+/11B+ state in unsubstituted
2u 3u
ble. IV. It can be observed that the energies obtained coronenebecomes symmetry alloweddue to substitution
in M = 1 sector calculations are lower. This can be andcorrespondstothelowestopticalstateinsubstituted
S
attributed to the larger Hilbert space of the M = 1 coronene.
S
state for a given DMRG cut-off compared to the Hilbert Itisevidentfromabovethattheeffectiveelectroncor-
space dimension of the M =0 space with χ =−1. In relation strength is not controlled by the symmetry of
S P
fact,theparitysymmetryalmostbisectsthecompletedi- themolecule;whilecoroneneisthemoleculewithhighest
rectproductspaceandtherefore,the effective size ofthe symmetry (D ), substitution by donor-acceptor groups
6h
10
0.1 0.1 of the substituents considered in this study.
0.08 Coronene 0.08
>i0.06 0.06 We have tabulated the bond order values of different
<sz0.04 0.04 bonds,asindexedinFig. 2inlowestenergystatesofeach
0.02 0.02
symmetry subspaces, for all the four molecules studied
0 0
0 2 4 6 8 10 12 14 16 18 20 22 24
0.3 0.3 here. The bond order values are related to bond lengths
>i0.2 Ovalene 0.2 asindicatedbyCoulson20;bondswithhigherbondorder
sz0.1 0.1 valuehasmoredoublebondcharacterwhilelowvaluede-
<
0 0 pictsthatthebondisprimarilyasinglebond. Thereper-
0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 cussionofsimilarbondordervaluesfordifferentbondsin
0.2 0.2
>i0.01.51 1,12-Benzoperylene 00..115 coronene establishes D6h symmetry of the ground state.
<sz0.05 0.05 This particular symmetry is not retained in the excited
0 0 statesindicatingdistortedequilibriumstructuresinthese
-0.050 2 4 6 8 10 12 14 16 18 20 22 -0.05 states. In excited states the bonds are symmetric about
Site index (i)
twoin-planeC axes,whichbisectsthemolecule. Intwo-
2
photon state, a quinonoid structure emerges at the cen-
FIG. 4. Spin density plot for coronene (top panel), ova-
tralpartof the molecule, while in the doubly degenerate
lene (middle) and 1,12-benzoperylene (bottom) in thelowest
opticalstates,whichhavealmostsimilarbondordersand
tripletstate. Thedatapointsareconnectedbyalineforsake
in the lowest triplet state the central benzene ring gets
of smooth vision.
squeezed compared to ground state and its connectivity
to the peripheralbenzene rings becomes weaker. Similar
structures appear even in substituted coronene. Small
variationinbondordersforbondsconnectedtodonoror
lowersthesymmetrytoC ,yet,theeffectivecorrelation
2h
acceptorsitesarequite expectedasthe effective hopping
strengthinboththesystemsasindicatedbythespingaps
amplitudes for bonds involving these sites will be differ-
and relative positions of the energy levels, are compara-
ent from the rest of the molecule. In ovalene, the edge
ble. Ovalene is also a lower symmetry molecule (D )
2h
bonds switch their character from double bond to sin-
than coronene, but appears to have large effective cor-
gle bond and vice versa in two-photon state and triplet
relation strength. On the contrary, 1,12-benzoperylene
state, relative to those in the ground state. In optical
(C )beingofsimilarsymmetrytosubstitutedcoronene,
2v
state, the edge bonds become more or less of same char-
seems to have rather low effective correlation strength
acterasallthebondordersbecomenearlyequal. Yet,the
compared with the latter.
mostnotablebondvariationinovaleneisforthebonds11
Spindensityatdifferentsitesinthelowesttripletstate,
and 12 and their symmetry partners, which also switch
which is the expectation value of hSˆzi in this state is
i their character like the edge bonds, but in all excited
plotted for different molecules in Fig. 4. In the coronene
states; other bonds in the inner part of the molecule re-
moleculewenotedthatthespindensityishighestatsites
mainalmostinvariant. However,the biggestvariationin
1, 8, 17 and 24 which are all related by symmetry of the
bondorderisobservedinthe lowesttripletstate of1,12-
molecule. Thefractionaldoubleoccupanciesattheabove
benzoperylene when comparedto the groundstate bond
mentioned sites are small (Fig. 3), which is consistent
orders. The Pauliblockingatthe sites 2and21converts
with the higher spin densities at these sites. However,
the bond 1 from a single bond to a double bond. Weak-
the bond order pattern along the periphery has period
ening of neighboringbonds can also be attributed to the
four andcorrespondto ( ==——== ) whichis un-
n Pauli blocking. Bond orders of most of the peripheral
usual as one would expect a quinonoidal structure from
andringbondsbecomenearlyequalwhilebondsparallel
thespindensitydata. Onthecontrary,thespindensities
to bond 1 remain approximately same in all states. In
are well-localized at sites 1 and 32 in ovalene and at 2
the two-photon state, overall expansion of the structure
and 21 in 1,12-benzoperylene (Fig. 2 and 4). In ovalene
is observed while in the optical state, most significant
and 1,12-benzoperylene,there are many sites with nega-
variation in bond order occurs in the lowermost ring as
tive spin densities unlike in the case of coronene. Again,
shown in Fig. 2. The bond order of the weakest bond of
we noted that the fractional double occupancy is least
the molecule increases both in the two-photon state and
at sites with high positive spin densities (Fig. 3). The
the opticalstate also,but rathersmall comparedto that
peripheral bond order in ovalene has period four with
in the triplet state.
kinks at the high spin density sites at the middle of the
The charge gap of a system is defined as the energy
molecule. Bond 11 and its symmetry partners (Fig. 2)
required to remove an electron from a neutral molecule
intheinteriortwobenzeneringsareratherweakascom-
andplacingitonanotherneutralmolecule. This process
pared to that in the ground state and we can expect
gives rise to two chargedspecies of different charges and
puckering of the structure in the triplet state. In 1,12-
can be represented by the reaction6:
benzoperylene also we find bond order pattern with pe-
riod four on the periphery if we exclude bonds involving 2M →P++P−
the high spin density sites. In substituted coronene, the
spin density distribution is not sensitive to the strength The charge gap for a system can be calculated by ob-