A simple scheme for finding magnetic aromatic hydrocarbon molecules
The identification of magnetic organic molecules is most relevant for technological applications and is of particular interest to develop rapid preliminary checks to identify likely candidates for both theoretical and experimental pursuits. Here is shown that an analysis based on a second-order perturbation treatment of electronic correlations for the Hubbard model qualitatively predicts the outcome of extensive and accurate computational methods. Based on these results it is proposed a simple protocol for screening molecules and identifying those worthy of a more sophisticated analysis on the magnetic nature of their ground states. Using this protocol two new magnetic molecules made from the combination of two naphthalene monomers and two azulene (34 carbon atoms, C 34 H 20 ) were identified. These molecules are shorter than the smallest magnetic fused-azulene oligomers so far reported. For further confirmation of this result, these molecules were also studied by means of density matrix renormalization group and density functional theory.
INTRODUCTION<!>METHODOLOGY Model Hamiltonians<!>Density matrix renormalization group<!>Density functional theory<!>RESULTS<!>New magnetic molecules<!>Building fuzed naphtha-azulene chains<!>Searching for magnetic molecules<!>(U ST < U r<!>(U<!>DFT and DMRG semi-empirical results by using Pariser-Parr-Pople Hamiltonian model<!>DISCUSSION<!>CONCLUSIONS
<p>The technological applications of carbon based materials spans from engineering of organic electronic devices [1][2][3][4] to biomedicine [5][6][7] . Purely organic magnetic compounds are a reality 8 , but the quest for functional materials with specific magnetic properties such as multiferroicity, superparamagnetism and ferromagnetism is just starting. Building blocks of such materials are molecules instead of single atoms so, depending on the number of atoms to be considered, these systems can be hard to explore computationally (for a recent review see Ref. 9 ). Quasi-unidimensional molecules, such as oligo-acenes (C 4n+2 H 2n+4 ) and fused azulene (C 8n+2 H 4n+4 ) oligomers (n stands for the number of monomers), are reliable physical systems that have attracted a great deal of attention. These planar molecules are built from laterally fused polycyclic aromatic hydrocarbons (as, for example in the configurations shown at Fig. 1 (a) and (b)). Oligo-acene (n-acene) is a one-dimensional graphene stripe with zigzag edges while fused azulene (or n-azulene) is a linear chain of conjugated rings with an odd number of carbon atoms, with opposite edges being zig-zag and armchair. In this work we present a protocol to easily screen a large amount of fused polycyclic aromatic hydrocarbons molecules in the search of an specific property, illustrated here with magnetism.  Solid state devices of n-acene display the promising feature of conductivity increasing with the number of aromatic rings [1][2][3][4] . These molecules have already been synthesized up to nine monomers 10,11 and are reported to be non-magnetic. Furthermore, n-acene beyond nine monomers with equal bond lengths have been theoretically investigated by different manybody techniques and have also been reported to not display magnetic polarization [12][13][14][15] . In marked contrast with benzene oligomers, these same numerical approaches predict fused azulene to display a singlet-triplet transition with spin polarization of the oligomer when six azulene molecules are linearly fused 9,15,16 . Recent density matrix renormalization group (DMRG) studies of fused azulenes using spin models also show many interesting features such as a reentrant nonmagnetic state when the exchange interaction strength of the common bond between rings is varied 17 . Fused azulene oligomers have also been suggested as possible organic multiferroic materials 15 , displaying as such both ferroelectric and ferromagnetic properties.</p><p>For large conjugated systems, it is well known that the σ-bands are energetically well separated from the π-bands, and their mixing can be ignored. Hence for dealing with electronic states of large conjugated systems only π-orbitals are considered. Low energy electronic excitations, understood as arising from delocalized electrons in the molecular orbital (πelectrons), have been appropriately described by the phenomenological Pariser-Parr-Pople (PPP) Hamiltonian 18,19 , extensively used for studying conjugated organic molecules. For symmetric fused-azulene molecules the ground state (GS) of the PPP model −with intersite Coulomb interaction parametrized by the Ohno 20 formula and fixed on-site Coulomb interaction U = 4.7 (in units of transfer integral between bonded sites)− is a singlet for up to 5-azulene while for 6-azulene to 11-azulene the GS is a triplet 15 . The fused azulene geometry has also been explored, in the strongly correlated regime, with the spin-1/2 antiferromagnetic Heisenberg model 17 and the magnetic ground state moment has also been observed to increase with the number of monomers.</p><p>While the PPP model correctly incorporates the long-range interactions which are only weakly screened in semiconductors, the truncated PPP model, namely, the Hubbard model further simplifies interactions by retaining only on-site electron-electron repulsions. In the strong correlation limit, the magnetic properties of the low-lying states are captured correctly by both models. Therefore, in most of our studies in this paper, we have used the simpler Hubbard model to explore the magnetism of the ground state of the quasi-unidimensional PAH molecules.</p><p>We employ DMRG and perturbation theory (PT) up to second-order to investigate the magnetic GS of oligo-acenes and fused azulene molecules as a function of the strength of electronic correlations. The first is a powerful numerical variational technique for calculating the ground state of interacting quantum systems and the second is a well established analytical approach. Different implementations of second-order perturbation theory (SPOT) has been successfully used to explore the GS of electronic correlated systems [21][22][23][24] . Therefore here we use second-order Rayleigh-Schrödinger PT to investigate magnetic transitions of fused benzene and azulene geometries modeled by using the Hubbard Hamiltonian.</p><p>Through DMRG studies we show in the RESULTS section that for n-fused azulenes, with n ≥ 7, molecules transition to a magnetic state at a low value of U/t, inside the validity of SOPT. This magnetic states survives for much larger values of U and even on the U/t = ∞ (Heisenberg) limit as shown on Ref. 15,17</p><!><p>We used the Hubbard and PPP models 18,19,25 to describe the electronic interactions in conjugated systems of the kind shown in Fig. 1. These model Hamiltonians contain a noninteracting part Ĥ0 and a term Ĥ1 that incorporates the electron-electron interactions:</p><p>The non-interacting part is a single-band tight-binding Hamiltonian,</p><p>which describes the kinetic energy with a fixed constant hopping t between nearest-neighbor carbon atoms (henceforth we will take t as the energy unit). In the framework of π-electron theories, the operator ĉ † iσ ( ĉiσ ) creates (annihilates) an electron of spin σ, localized in a p z -orbital at site i.</p><p>The interacting part of Ĥ can be written in the form Ĥ1 = ĤU + ĤV .</p><p>(</p><p>The Hubbard model is obtained if only the first term H U is considered, Ĥ1 = ĤU . This term accounts for the interaction between electrons in the same π-orbital:</p><p>where U is the magnitude of the on-site Coulomb strength and N denotes the number of sites (π-orbitals), which is equal to the number of electrons for all the molecules we have studied. The local particle number operator for electrons of spin σ at a site i is niσ = ĉ † iσ ĉiσ , and the total occupation of a π-orbital is ni = ni↑ + ni↓ . The factor 1/2 in Equation (4) fixes the chemical potential for occupation ni = 1 on each site for a half-filled band.</p><p>The PPP model is realized when the second term in Equation ( 3) is also included:</p><p>The inter-site interactions between all electrons V ij are parametrized using the Ohno 20 formula.</p><p>For the PPP model the transfer integral between bonded sites is taken to be t = First-and second-order perturbation theory for the Hubbard model</p><p>In the limit of weakly correlated electrons we can treat the on-site Coulomb interaction (including only the Hubbard term in Ĥ1 ) as a small perturbation in the total energy. The energy for the unperturbed system is simply the energy of a tight-binding model</p><p>where |n (0) is the ground state of the unperturbed system, which can easily be obtained since Ĥ0 is exactly diagonalizable.</p><p>Using the Rayleigh-Schrödinger perturbation theory 26 up to second-order, we can write the energy of the perturbed system as:</p><p>where λ = U/t, λE</p><p>with |k (0) a state of the unperturbed system. Notice that Ĥ1 is not diagonal in the {|k (0) } basis.</p><p>Details on perturbation theory for the Hubbard model in the low-correlated regime and at the singlet (S=0) and triplet (S=1) state are available in the Supplementary Material 1 .</p><p>Values of noninteracting low-lying energy E (0) 0 , first E</p><p>(1) 0</p><p>and second order E</p><p>(2) 0 correction to the lowest energy for different oligo-acenes and fused azulene oligomers can be found in Table 1 for M S = 0 and M S = 1 spin projections. We take U/t 2 as the upper limit of validity for the results obtained from perturbation theory, as for larger U/t values the condition H 1 H 0 is no longer fulfilled 27,28 .</p><p>On second-order PT the spin contamination, S 2 , is zero for the S = 0 state 29 . This results from the fact that the correcting wave function |n (1) can be written as:</p><p>(see 30 ). As H 0 and H 1 commute with S 2 , |n (1) is also an eigenstate of S 2 with the same eigenvalue as |n (0) , zero for this case.</p><p>Henceforth we will refer to first-order PT and second-order PT, for the Hubbard model in the regime of weakly correlated electrons, by the acronyms FOPT and SOPT, respectively.</p><!><p>Although the density matrix renormalization group (DMRG) algorithm was first presented to treat a one-dimensional problem 31 , different low-dimensional systems can also be accurately treated with this approach 32,33 for limited extension in the second dimension. For the fused benzene and fused azulene molecules in Fig. 1 (a-b), we can map the Hubbard Hamiltonian with first-neighbours interactions, into a one-dimensional Hamiltonian with up to secondneighbours interactions, which is easily implemented on a DMRG algorithm. Also, in the Hubbard Hamiltonian, the potential energy does not depend on the spatial coordinates of 1 See supplementary material (1) for more detail on the Rayleigh-Schrödinger perturbation theory for the Hubbard model at the weakly-correlated regime. A script for performing this calculation is also available in supplementary material (2)  the molecule. Then it is usual to simplify the computing work by taking a ladder-like representation of the system, as illustrated in Fig. 1 (c-d).</p><p>Oligomers described by the Hamiltonian in Eq. 1 have SU(2) symmetry so we can take advantage of the degeneracy of spin projections M S for a total spin S, and compute the lowlying energy states of the system in a chosen subspace. If the ground state has spin S, then the lowest energies for different spin projections subspaces M S satisfy E 0 (M S ) = E 0 (M S + 1) for M S = −S, (−S +1), ..., (S −1) and the low-lying energy E 0 (M S = S) < E 0 (M S = S +1) 34,35 .</p><p>In order to establish the accuracy of our DMRG calculations we compared DMRG results</p><p>for the energy at U = 0 with the exact solution of the tight-binding Hamiltonian. We computed the GS energy with different sizes of the Hilbert space, by varying the cutoff in the number of block states, m, between 800 and 1200. In these conditions the energy precision is in the fourth decimal digit. The relative energy error is comparable to the DMRG weight lost which is kept lower than 10 −5 in the worst cases.</p><!><p>The magnetic nature of the addressed molecules was independently obtained by means of density functional theory (DFT) calculations as implemented in the nwchem code 36 . The convergence criteria were set to the nwchem default values and the 6-311+G* basis set was chosen to expand the electronic wave function. For each molecule, spin-polarized geometry optimizations and electronic structure calculations were carried out, enforcing in turn both the spin-singlet and the spin-triplet configurations. The spin configuration with the lowest energy was identified as the actual GS. The electronic correlations enter the DFT calculation through the exchange-correlation functional. We worked with the B3LYP functional 37 and also with the PBE0 one 38 , both belonging to the state-of-the-art hybrid family. This was intended as a means for checking whether the functional choice affected the obtained magnetic nature of the molecule's GS and hence the agreement between DFT and the predictions made by SOPT and DMRG.</p><!><p>Oligo-acenes and fuzed-azulene molecules in the Hubbard Model</p><p>In this section we consider DMRG simulations and Rayleigh-Schrödinger PT up to secondorder, to compute the energy and the spin gap of fused-azulene and oligo-acene molecules with up to 74 π-orbitals (9 azulene or 18 benzene units). We explore the Hubbard model phase space with correlations ranging from the non-interacting limit U = 0 to U = 6t, and systems of different size. We calculate the unperturbed lowest energy, E</p><p>0 (M S ), and perturbed energies corrections in the first-order, E</p><p>(1) 0 (M S ), and second-order E</p><p>(2) 0 (M S ), for the system with spin projections M S = 0 (singlet state) and M S = 1 (triplet state). The results obtained are summarized in Table 1.</p><p>0 (M S ) and second E</p><p>(2) 0 (M S ) order coefficients of perturbation theory and critical value of electronic correlation obtained by</p><p>ST ), and DMRG (U ST ) calculations respectively (all in units of t) for different number of monomers (N atoms) in oligo-acene and fused-azulene molecules.</p><p>The ground state of the system will be either a triplet (S = 1) or a singlet (S = 0) depending on the U/t ratio. For a given U/t value, the energy difference between spin configurations ∆E(M S ) = E 0 (M S = 1) − E 0 (M S = 0) defines the singlet to triplet spin gap.</p><p>Upon exploring the spin gap functional dependence on U/t, we observe that the GS is a singlet for U/t → 0 but, as U/t increases, we may come across a critical value of U/t = U ST , for which the spin gap closes. A further increase in U/t would then change the magnetic nature of the GS to a triplet. This is a singlet-triplet quantum phase transition.</p><p>The dependence of the low-lying energies E 0 (M S = 0), E 0 (M S = 1) and the spin gap ∆E on the interaction U/t is illustrated in Fig. 2 (a) and (b) for a two azulene system (18 atoms, see the molecule structure in Fig. 1). Panel (a) shows the energy computed on the M S = 0 and M S = 1 spaces. We observe that the E 0 (M S = 0) curve lies slightly below the curve for E 0 (M S = 1), which indicates a singlet GS and a finite spin-gap. Second-order PT (SOPT) data is closer to DMRG results than first-order PT (FOPT). In panel (b) we plot the singlet-triplet spin gap as a function of electronic correlations. A positive finite gap indicating a singlet state is found both with DMRG and the PT calculations, in agreement with Refs. 9,15,16 . Note that the FOPT gap points to a singlet-triplet transition, though for values of U ST for which the approximation is not valid (not displayed in the plot).</p><p>In Fig. 3, we plot the singlet-triplet spin gap as a function of electronic correlations for molecules with 5 azulenes, Fig. 3 (a), and 6 azulenes, Fig. 3 (b). The gap is obtained from FOPT, SOPT and DMRG calculations. Mark that for PT we compute the spin gap as the energy difference between states with S = 1 and S = 0 and hence a triplet GS is signaled by ∆E(M S ) < 0. In contrast, for DMRG the spin gap is given by the difference between states with spin projection M S = 1 and M S = 0 and thus it is ∆E(M S ) = 0 that signals a triplet state. For U/t < 2, SOPT and DMRG give very similar and finite gaps, indicating singlet GSs in both cases. Interestingly, as depicted in Table 1, SOPT estimates U ST below 2t when the chain is larger than 6 azulene units, and the U ST agreement between the different approaches improves with the number of monomers.</p><p>In Fig. 4 we show that for oligo-acene with 10 rings (10-acene with 42 carbon atoms), FOPT incorrectly predicts a singlet-triplet transition while SOPT correctly predicts that the ground state for oligo-acenes is always a singlet. DMRG and SOPT exhibit a significant spin gap while this gap is closed under FOPT. We checked these results for oligo-acenes with up to 18 benzene rings (see Table 1) and we can reproduce the expected behavior of spin excitations of oligo-acene chains, i.e. a finite spin gap which indicates that the ground state of the system is a singlet, independently of the electronic correlation magnitude and system size 15 .</p><!><p>In the previous section, our analysis was restricted to PAH molecules consisting of either only-azulene or only-benzene monomers (blocks). However, it has been recently reported that the stability of synthesized PAH chains increased upon the introduction of azulene molecule into fused benzene oligomers 39 . Motivated by this result, in this section we will include in our considerations chains which are built from the combination of naphthalene (two fused benzene rings) and azulene molecules to a total of four fused building blocks (eight rings). The molecules share the formula C 34 H 20 , with 34 π-orbitals.</p><!><p>The molecules to be addressed in this section can be represented topologically as two-leg ladders with slanted rungs corresponding to the intra-azulene link. This scheme is graphically explained in Fig. 5, where we depict the ladder-like representations for naphthalene and azulene molecules to be used as building blocks for the new PAH chains (Fig. 5 (a)). Also in this figure, we define a code for each building block considered. Configuration "0" stands for a naphthalene block, while for azulene the two possible building blocks, combined with the two possible fused geometries result in 4 possible configurations named "1", "2" , " 1", " 2". Examples of ladder-like representation and code for two molecules having eight-rings are provided in Figs. 5 (b) and (c). Regarding the latter molecule, mark that it includes a building block named "X" made of a four-carbon ring and an eight-carbon ring. This molecule will be briefly mentioned in subsection Searching for magnetic molecules and Fig. 6 but, as it is not built from fused naphthalenes and azulenes, it is not part of the set of molecules on which this study focuses.</p><p>For a molecule built of two blocks from the set in Fig. 5 (a) (18 carbon atoms) we have 5 2 possibilities, many of which are related by symmetry. Mark however that there exist more possible fusing geometries not included in our selected block set, giving for example configuration 41 of Ref. 40 . For 4 fused blocks from the set 5 4 = 625 possible molecules exist, out of which 169 are not related by symmetry. These numbers are intended to give an idea of the prohibitive calculation times that would be required to analyze, with numerically costly methods, the magnetic nature of the GS for every possible molecule.</p><!><p>In this subsection we will discuss the dependence of the spin gap (and thus of the GS magnetic nature) on the strength of the on-site Coulomb interaction, given by the U value.</p><p>The different scenarios will be discussed in terms of the ability of SOPT to identify a possible triplet GS and the confirmation provided by DMRG calculations regarding this feature. Note that U is not a widely tunable parameter but is fixed for a molecule. Our analysis intends to provide information on the nature of the GS regardless of the specific value of U .</p><p>In order for SOPT to be able to identify a triplet GS, the spin gap must close for a value of the on-site Coulomb interaction such that U ST /t 2 (weak electronic correlations).</p><p>However, the condition U ST /t 2 does not guarantee that the GS remain magnetic for large U values. In fact, it has been shown 17 that in skewed ladder systems, the behavior of the ground state as a function of the Hamiltonian parameters can be quite surprising.</p><p>In fused azulene systems, in fact reentrant singlet phase has been found. Therefore if the SOPT ground state transitions to a triplet state at a given value of the on-site electronic correlation U = U ST , then there could be a value for U = U r T S > U ST for which we could expect that the ground state becomes a singlet again. In such a scenario, if U r T S /t 2 the spin configuration of the GS for a real molecule with a large value of the on-site interaction may well be non-magnetic.</p><p>From the previous discussion, we conclude that the second-order Rayleigh-Schrödinger PT for the Hubbard model at the weakly correlated regime can signal for a possible triplet GS as long as U ST /t 2 U r T S /t. In Fig. 6 we compare second-order perturbation theory (dashed lines) with DMRG values (symbols) of the triplet-singlet spin gap for a selection of linear PAH molecules with the same molecular formula (C 34 H 20 ) and different geometric configurations (see Figs. 5 and 7 for the labeling and geometry of these molecules). The chosen cases are illustrative of the three different possibilities for the GS of the system as a function of the electronic correlation depending whether U ST and U r T S are above or below the PT validity limit of 2t (U r T S is always greater than U ST ):</p><p>1. (2t < U ST < U r T S ) In the top panel we see that a molecule could always be in a singlet state as the singlet-triplet spin gap never goes to zero (as for the 4-azulene oligomer, 2222 configuration), or it goes to zero for value of U much larger than shown at Fig. 6.</p><!><p>T S < 2t) Still in Fig. 6 (a), we see that the spin gap could close for a small value of the on-site Coulomb interaction (U ST /t < 2) and raise again for values of U still in the range of validity of second-order perturbation theory (U r T S /t 2 as it is the case for the X100 molecule). Notice that the X100 molecule is not a combination of just naphthalene and azulene monomers. It includes a monomer, labeled X, consisting of a 4-atom ring and an 8-atom one (this configuration is being addressed for illustrative purposes).</p><!><p>In the bottom panel we see that the spin gap could close for a small value of U (U ST /t < 2) and raise again for values of U away from the range of validity of second-order perturbation theory (U r T S /t > 2 as for the 2112 molecule). Lastly, also in Fig. 6 (b), we see that the spin gap could close for a small value of U (U ST /t < 2)</p><p>and never again open in the studied range of U/t (as for the 0210 molecule).</p><p>Molecules belonging to categories 1 and 2, illustrated in Fig. 6 (a), can be expected, from a SOPT analysis, to have a singlet GS in the U/t > 2 regime. Notice that this criteria can miss some magnetic configurations like the 5-azulene molecule. Molecules in category 3, illustrated in Fig. 6 (b), are identified by SOPT as good candidates for further sophisticated analysis.</p><p>These considerations lead us to propose the following scheme to search for magnetic molecules:</p><p>1. Check whether in SOPT the system becomes magnetic for small values of U. To be explicit, check if U ST /t < 2.</p><p>2. Check whether in SOPT the system remains magnetic in the range of validity of this method. To be explicit, check if U r T S /t > 2.</p><p>3. If the molecule fulfills the above conditions, it is a good candidate to present magnetism and suited for testing with sophisticated numerical calculation by using DMRG or DFT methods.</p><p>Notice that despite not capturing all magnetic molecules (for example 5 and 6 fused azulene oligomers), the importance of this scheme is that it allows, through a very simple second-order calculation (the first two steps), to check whether it is advisable or not to perform a more sophisticated calculation like DMRG to address the magnetic nature of the GS.</p><p>Making use of this protocol, the SOPT analysis for the Hubbard Hamiltonian at weakly correlated regime, we found no triplet GS for 2 or 3 fused blocks. However when 4 molecules are fused the SOPT analysis signals seven (out of the 169 configuration tested) as possibly magnetic. DMRG calculations, for the same Hamiltonian model, were then carried out for the seven candidates and two of them were identified as magnetic. In Table 2, the values of U ST and U r T S for the seven configurations are shown. In order to confirm the magnetic nature of the ground state regardless of the U value, in the large U/t region we performed DMRG calculation for the reference value U/t = 5 15 , employing the Hubbard model instead of the PPP Hamiltonian. Indeed, configurations 0210 and 0120 present a U ST value in the region of validity of SOPT and a closed gap under DMRG. The remaining five configurations display very small, though not closed spin gaps.</p><!><p>The seven molecules identified as possibly magnetic by SOPT for the Hubbard model in the weak correlation regime, were also analyzed, within the DMRG approach with the PPP Hamiltonian. The values of the singlet-triplet spin gap are shown at the Table 2 where we can see that three of the seven configurations were found to be magnetic, including the two pointed as magnetic when using the Hubbard model. These seven molecules (from Table 2) were further examined by means of DFT calculations, with the exchange-correlation energy approximated by B3LYP and PBE0 hybrid functionals. Their geometries were confined to be planar and are shown in Fig. 7</p><p>ST , and triplet-singlet, U  and ∆E DFT P BE0 , computed with DFT respectively for the B3LYP and PBE0 exchangecorrelation functionals (all in units of Kcal/mol). Notice that DFT computes the energy difference between states belonging to the S = 1 and S = 0 subspaces. So, for DFT, ∆E < 0 signals a triplet state. Molecules identified as magnetic by DMRG (with both models) and DFT are shown in bold. spin gaps, are shown in the last two columns of Table 3, and are contrasted with those obtained with DMRG simulations.</p><!><p>In subsection Oligo-acenes and fuzed-azulene molecules in the Hubbard Model we see that first-order PT incorrectly predicts a singlet-triplet transition as illustrated for both oligo-acenes and oligo-azulene molecules. In contrast, second-order PT qualitatively reproduces the DMRG predictions: a singlet ground state for n-acene and a singlet-triplet GS transition for n-azulene (with n ≥ 7). We observe also that the agreement between different approaches improves with the oligomer length. In subsection Searching for magnetic molecules we report that for some molecular configurations the magnetic state (in the Hubbard model) becomes the ground state for low values of U/t (inside the validity of SOPT) and remains stable for the largest values of the Coulomb repulsion numerically tested. These two results lead us to propose a protocol based on SOPT to screen a large set of molecules for those configurations whose ground state change from singlet to triplet inside the validity of SOPT. With the additional requirement that within SOPT the magnetic state remains stable up to the upper validity limit of the approximation, a small set of geometries is selected.</p><p>The molecules selected by this protocol could well be no magnetic in the large value of the interaction (U/t ∼ 5) to be assigned inside the Hubbard model. Or The Hubbard model could not be a fair description of the electronic interactions in these systems. Even more, these interactions could lead to deformations in the geometry that stabilize a singlet (non magnetic) ground state. We checked for all these issues using as a test case molecules made up of the fusion of four azulene or naphthalene monomers C 34 H 10 (169 configurations):</p><p>On Table 2 we see that indeed DMRG calculations for U/t = 5 confirm that inside the Hubbard model two out of the seven configurations selected by SOPT (0120 and 0210) are in fact magnetic. On Table 3 we show that the PPP model, that does include long range interactions, confirms that these set has three magnetic molecules, where two agree with those confirmed by the Hubbard model (adding the 1212 configuration). In accord with DMRG (with both Hubbard and PPP models), molecules 0120 and 0210 are also identified as magnetic by DFT with both B3LYP and PBE0 calculations, yielding large (negative) spin gaps, when compared with the remaining five molecules. Molecules 2121 and 2122 present the largest (positive) gaps both in DMRG and DFT and are thus confirmed as non-magnetic.</p><p>Regarding the molecule 1212 (non-magnetic according to DMRG for the Hubbard model and magnetic when using PPP model) and 2112, we observe that DFT gives small gaps, the sign of which depends on the functional chosen. Finally, by using DFT, the 2212 molecule is predicted to be magnetic although with a gap much smaller than those of the 0120 and 0210 molecules. We thus reckon DFT cannot confirm or reject the DMRG prediction of the non-magnetic nature of the GS for these three latter cases. The magnetic stability of triplet states on these molecules as computed by DFT could be favoured also by the performed optimization of the geometry 12 , not done with the Hubbard and PPP model Hamiltonians.</p><p>We see that the proposed protocol allow us to search magnetic molecules on a largely reduced set of configurations where indeed a significant proportion (two out of seven) are confirmed to be magnetic by state of the art numerical methods. Notice that this protocol allow us to find at least two new magnetic molecules with non trivial geometries that are shorter (four monomers) than the ones so far reported that consist of at least six azulene monomers (C 50 H 28 ) 15 . The remaining 5 configurations proposed by the protocol could be magnetic or have a very small singlet-triplet gap which is also an interesting feature (for solar cells, as an example).</p><!><p>In this article we have introduced a protocol for analyzing the magnetic nature of the ground state of PAH molecules, particularly oligo-acenes and oligo-azulenes. The protocol is intended as a means to minimize the otherwise prohibitive numerical cost of addressing every possible molecule with computationally demanding techniques such as DMRG or DFT. It is based on a second-order Rayleigh-Schrödinger perturbation treatment of electronic correlations, on a Hubbard Hamiltonian with weak on-site Coulomb interactions. A preliminary SOPT study of the molecules under consideration provides a cheap and efficient spotting method for likely π-conjugated magnetic molecules. The fulfillment of the conditions of a singlet-triplet gap closure inside the validity rang of pertubative theory (U/t < 2) and a triplet-singlet gap opening outside of the validity of PT signals a candidate worthy of further examination with more sophisticated (and costly) machinery such as DMRG or DFT.</p><p>The DMRG calculation is then performed for a strong correlation scenario, (U/t = 5) for the Hubbard model and with semi-empirical PPP Hamiltonian. In combination with the qualitative accord which is obtained between DMRG and SOPT spin gaps in the perturbative range, a closed spin gap in the latter calculation is indicative of a magnetic GS for the molecule under consideration. This feature was also supported by DFT calculations.</p><p>We have profited from this scheme to study more than one hundred possibilities of linear molecules (with 34 carbon atoms) built from combining naphthalenes and azulenes. The SOPT analysis for the Hubbard Hamiltonian in the weak correlation limit yielded seven candidates, of which two were confirmed as magnetic by DMRG (with Hubbard and PPP Hamiltonian models) and DFT. These two new magnetic molecules, made of two naphthalenes and two azulenes, are shorter than the smallest magnetic fused azulene oligomers so far reported.</p><p>Finally, although for these molecules perturbations around the PPP model might be more appropriate than around the Hubbard model (specially away from half filling 9 ) and the fact that our protocol does miss out some magnetic molecules such as 5 and 6 fused azulenes, we have nevertheless successfully predicted new magnetic molecules as confirmed by our DMRG (PPP) and DFT calculations. The method could also be used to study other properties such as electric polarization.        2 and 3.</p><p>Dark grey and light grey stand, respectively, for carbon and hydrogen atoms. Configurations (e) and (f) are confirmed to be magnetic within all different approaches we use to scrutinize the second order PT predictions.</p>