Léo Dontot

Fernand Spiegelman

Mathias Rapacioli
email: mathias.rapacioli@irsamc.ups-tlse.fr

Structures and Energetics of Neutral and Cationic Pyrene Clusters

  L'archive ouverte pluridisciplinaire

Introduction

The investigation of Polycyclic Aromatic Hydrocarbons (PAHs) clusters properties from both experimental or theoretical sides is motivated by the understanding of processes involved in di↵erent scientific domains. In astrophysics, PAHs are suspected to be responsible for a family of infrared emission bands [START_REF] Léger | Identification of the 'unidentified' IR emission features of interstellar dust?[END_REF][START_REF] Allamandola | Polycyclic aromatic hydrocarbons and the unidentified infrared emission bands -Auto exhaust along the Milky Way[END_REF] and to contribute to UV-Visible absorption features, known as Di↵use Interstellar Bands. [START_REF] Van Der Zwet | Polycyclic aromatic hydrocarbons and the di↵use interstellar bands[END_REF][START_REF] Leger | Are polycyclic aromatic hydrocarbons the carriers of the di↵use interstellar bands in the visible?[END_REF][START_REF] Crawford | Ionized polycyclic aromatic hydrocarbons and the di↵use interstellar bands[END_REF][START_REF] Salama | Assessment of the Polycyclic Aromatic Hydrocarbon-Di↵use Interstellar Band Proposal[END_REF] In order to explain spectral bands broadening observed in di↵erent objects, [START_REF] Rapacioli | Spectroscopy of polycyclic aromatic hydrocarbons and very small grains in photodissociation regions[END_REF][START_REF] Berné | Analysis of the emission of very small dust particles from Spitzer spectro-imagery data using blind signal separation methods[END_REF][START_REF] Berné | Extended Red Emission and the evolution of carbonaceous nanograins in NGC 7023[END_REF] it has been proposed that these molecules could aggregate in regions where they are protected from stars UV photons. PAH clusters have also been proposed to be at the origin of some of the DIBs and the extended red emission. [START_REF] Rhee | Charged polycyclic aromatic hydrocarbon clusters and the galactic extended red emission[END_REF][START_REF] Lemmens | Far-IR and UV spectral signatures of controlled complexation and microhydration of the polycyclic aromatic hydrocarbon acenaphthene[END_REF] They are also suspected to play a role in the growth of astrophysical PAHs themselves thanks to the creation of chemical bonds between the di↵erent units following a photoexcitation. [START_REF] Chen | Formation of Covalently Bonded Polycyclic Aromatic Hydrocarbons in the Interstellar Medium[END_REF][START_REF] Zhen | Laboratory Photochemistry of Pyrene Clusters: An E cient Way to Form Large PAHs[END_REF] Few studies have investigated the competition between formation and evaporation of such clusters under interstellar conditions. [START_REF] Rapacioli | Formation and destruction of polycyclic aromatic hydrocarbon clusters in the interstellar medium[END_REF][START_REF] Joblin | Insights into the condensation of PAHsin the envelope of IRC +10216, PAHs and the Universe[END_REF][START_REF] Montillaud | Absolute evaporation rates of non-rotating neutral polycyclic aromatic hydrocarbon clusters[END_REF] The latter studies were restricted to neutral systems, whereas singly ionized clusters are expected to be more strongly bonded than their neutral counterparts and easier to form due to long range polarisation interactions between a cationic cluster and a neutral PAH. Many PAH clusters studies are also motivated by the understanding of their role as a precursors to soot formation in flames. [START_REF] Fiedler | The e↵ect of temperature on nanoparticle clustering[END_REF][START_REF] Eaves | The importance of reversibility in modeling soot nucleation and condensation processes[END_REF][START_REF] Violi | Soot primary particle formation from multiscale coarse-grained molecular dynamics simulation[END_REF][START_REF] D'anna | A model of particle nucleation in premixed ethylene flames[END_REF][START_REF] Totton | A quantitative study of the clustering of polycyclic aromatic hydrocarbons at high temperatures[END_REF][START_REF] Raj | A study on the coagulation of polycyclic aromatic hydrocarbon clusters to determine their collision e ciency[END_REF][START_REF] Veshkini | Application of PAHcondensation reversibility in modeling soot growth in laminar premixed and nonpremixed flames[END_REF][START_REF] Eaves | Assessing relative contributions of PAHs to soot mass by reversible heterogeneous nucleation and condensation[END_REF][START_REF] Adkins | Computed electronic structure of polynuclear aromatic hydrocarbon agglomerates[END_REF][START_REF] Giaccai | Examination of the electronic structure of oxygencontaining PAH dimers and trimers[END_REF][27][START_REF] Herdman | Intermolecular Potential Calculations for Polynuclear Aromatic Hydrocarbon Clusters[END_REF][START_REF] Saggese | Kinetic modeling of particle size distribution of soot in a premixed burnerstabilized stagnation ethylene flame[END_REF][START_REF] Aubagnac-Karkar | Soot particles inception and PAH condensation modelling applied in a soot model utilizing a sectional method[END_REF][START_REF] Totton | Modelling the internal structure of nascent soot particles[END_REF][START_REF] Wang | Revealing the molecular structure of soot precursors[END_REF][START_REF] Mosbach | Towards a detailed soot model for internal combustion engines[END_REF][START_REF] Sabbah | Exploring the Role of PAHs in the Formation of Soot: Pyrene Dimerization[END_REF] The knowledge of PAH clustering followed by oligomerization is necessary to build flames kinetic models. Let us also mention that PAH clusters are investigated as prototypes to understand organic crystal properties or to design organic solar cells devices. Electronic gaps and transport properties within clusters are investigated in that context. [START_REF] Adkins | Computed electronic structure of polynuclear aromatic hydrocarbon agglomerates[END_REF][START_REF] Heinemann | Coarse-grained electrostatic interactions of coronene: Towards the crystalline phase[END_REF][START_REF] Grančič | Evaluation of All-Atom Force Fields for Anthracene Crystal Growth[END_REF][START_REF] Sancho-García | Determining the cohesive energy of coronene by dispersion-corrected DFT methods: Periodic boundary conditions vs. molecular pairs[END_REF][START_REF] Fedorov | Properties of crystalline coronene: Dispersion forces leading to a larger van der Waals radius for carbon[END_REF][START_REF] Fioressi | E↵ects of cluster formation on spectra of benzopyrene and benzopyrene[END_REF][START_REF] Sancho-García | Theoretical study of stability and chargetransport properties of coronene molecule and some of its halogenated derivatives: A path to ambipolar organic-based materials[END_REF][START_REF] Sancho-García | Charge-transport properties of prototype molecular materials for organic electronics based on graphene nanoribbons[END_REF][START_REF] Pérez-Jiménez | Conductance Enhancement in NanographeneGold Junctions by Molecular -Stacking[END_REF][START_REF] Yoshida | Conducting Columns of Highly Symmetric Coronene, The Smallest Fragment of Graphene[END_REF][START_REF] Sanyal | E↵ect of Imide Functionalization on the Electronic, Optical, and Charge Transport Properties of Coronene: A Theoretical Study[END_REF][START_REF] Scholz | Quantifying charge transfer energies at donor-acceptor interfaces in small-molecule solar cells with constrained DFTB and spectroscopic methods[END_REF][START_REF] Darghouth | Assessment of Density-Functional Tight-Binding Ionization Potentials and Electron A nities of Molecules of Interest for Organic Solar Cells Against First-Principles GW Calculations[END_REF] From the experimental side, mostly driven by the identification of astrophysical nanograins, the spectroscopic properties of PAH clusters have been reported in the IR-FIR as well as in the UV range. [START_REF] Lemmens | Far-IR and UV spectral signatures of controlled complexation and microhydration of the polycyclic aromatic hydrocarbon acenaphthene[END_REF][START_REF] Fioressi | E↵ects of cluster formation on spectra of benzopyrene and benzopyrene[END_REF][47][START_REF] Ricca | The Infrared Spectroscopy of Neutral Polycyclic Aromatic Hydrocarbon Clusters[END_REF][START_REF] Friha | Electronic spectra of cationic PAH and PAH clusters[END_REF] Most of the experiments devoted to the characterization of PAH clusters structural properties reported stacked structures for small clusters, sometime suggesting a Parallel Displaced (PD) pattern. [START_REF] Totton | Modelling the internal structure of nascent soot particles[END_REF][START_REF] Goulart | The structure of coronene cluster ions inferred from H2 uptake in the gas phase[END_REF][START_REF] Birer | Dimer formation of perylene: An ultracold spectroscopic and computational study[END_REF] For instance, Wang et al. [START_REF] Wang | Revealing the molecular structure of soot precursors[END_REF] produced nanoclusters of hydrocarbon molecules by annealing coronene films deposited on Pt(111) and observed stacking patterns with scanning tunneling microscopy. Ion mobility spectra suggest a formation as a single stack perturbed by side molecule when the number of units increases. [START_REF] Beitz | Ion Mobility Spectrometric Investigation of Aromatic Cations in the Gas Phase[END_REF] Energetic properties were also investigated, for instance measuring ionization energies of PAH clusters. [START_REF] Joblin | Size E↵ect in the Ionization Energy of PAH Clusters[END_REF] The stability of PAH clusters at thermal equilibrium or submitted to low or high energetic collision has also been investigated. [START_REF] Schmidt | Coronene cluster experiments: Stability and thermodynamics[END_REF][START_REF] Gatchell | Failure of hydrogenation in protecting polycyclic aromatic hydrocarbons from fragmentation[END_REF][START_REF] Holm | Ionization and fragmentation of cold clusters of PAH molecules -collisions with keV ions[END_REF][START_REF] Johansson | Ionization and fragmentation of polycyclic aromatic hydrocarbon clusters in collisions with keV ions[END_REF][START_REF] Holm | Ions Colliding with Cold Polycyclic Aromatic Hydrocarbon Clusters[END_REF][START_REF] Gatchell | Ions colliding with mixed clusters of C 60 and coronene: Fragmentation and bond formation[END_REF][START_REF] Gatchell | Knockout driven reactions in complex molecules and their clusters[END_REF][START_REF] Rousseau | Low-energy ions interacting with anthracene molecules and clusters[END_REF][START_REF] Seitz | Polycyclic aromatic hydrocarbon-isomer fragmentation pathways: Case study for pyrene and fluoranthene molecules and clusters[END_REF] Oligomerization of PAHs within the cluster induced by such a collision, [START_REF] Delaunay | Molecular Growth Inside of Polycyclic Aromatic Hydrocarbon Clusters Induced by Ion Collisions[END_REF] photoabsorption [START_REF] Zhen | Laboratory Photochemistry of Pyrene Clusters: An E cient Way to Form Large PAHs[END_REF][START_REF] Gámez | Ultraviolet laser desorption/ionization mass spectrometry of single-core and multi-core polyaromatic hydrocarbons under variable conditions of collisional cooling: insights into the generation of molecular ions, fragments and oligomers[END_REF][START_REF] Bréchignac | Photoinduced products from cold coronene clusters[END_REF] or high pressure conditions [START_REF] Chanyshev | Oligomerization and carbonization of polycyclic aromatic hydrocarbons at high pressure and temperature[END_REF][START_REF] Chanyshev | Temperature-induced oligomerization of polycyclic aromatic hydrocarbons at ambient and high pressures[END_REF] has also been reported.

From the theoretical side, due to the size of the systems, most of the reported investigations were performed at the DFT level and were most of time limited to dimers. [START_REF] Fioressi | E↵ects of cluster formation on spectra of benzopyrene and benzopyrene[END_REF][START_REF] Sancho-García | Charge-transport properties of prototype molecular materials for organic electronics based on graphene nanoribbons[END_REF][START_REF] Ricca | The Infrared Spectroscopy of Neutral Polycyclic Aromatic Hydrocarbon Clusters[END_REF][START_REF] Birer | Dimer formation of perylene: An ultracold spectroscopic and computational study[END_REF][START_REF] Zhao | A Prototype for Graphene Material Simulation; Structures and Interaction Potentials of Coronene Dimers[END_REF][START_REF] Rapacioli | Correction for dispersion and Coulombic interactions in molecular clusters with density functional derived methods: Application to polycyclic aromatic hydrocarbon clusters[END_REF][START_REF] Obolensky | Interplay of electrostatic and van der Waals forces in coronene dimer[END_REF][START_REF] Bartolomei | Modeling Coronene Nanostructures: Analytical Potential, Stable Configurations and Ab Initio Energies[END_REF][START_REF] Piacenza | Van der Waals Complexes of Polar Aromatic Molecules: Unexpected Structures for Dimers of Azulene[END_REF] One should however mention a few studies at the MP2 [START_REF] Birer | Dimer formation of perylene: An ultracold spectroscopic and computational study[END_REF][START_REF] Piacenza | Van der Waals Complexes of Polar Aromatic Molecules: Unexpected Structures for Dimers of Azulene[END_REF][START_REF] Yurtsever | Stacking of triphenylene: characterization of the potential energy surface[END_REF] or SAPT [START_REF] Podeszwa | Interactions of graphene sheets deduced from properties of polycyclic aromatic hydrocarbons[END_REF] levels for dimers and the investigation at the DFT level of dimers with an hetero-atom. [START_REF] Giaccai | Examination of the electronic structure of oxygencontaining PAH dimers and trimers[END_REF][START_REF] Philpott | Bonding and magnetism in transition metal sandwich structures with the aromatic hydrocarbon coronene C24H12 outer layers[END_REF][START_REF] Philpott | Chemical bonding in metal sandwich molecules MnR2 with R=pyrene C16H10 and tetracene C18H12[END_REF][START_REF] Philpott | Transition metal sandwich molecules with large (Cn, n24) zigzag poly aromatic hydrocarbons[END_REF] Although benzene dimers present a low energy T-shape conformation, [START_REF] Lee | Understanding of Assembly Phenomena by Aromatic; Aromatic Interactions: Benzene Dimer and the Substituted Systems[END_REF][START_REF] Podeszwa | Potential Energy Surface for the Benzene Dimer and Perturbational Analysis of Pi-Pi Interactions[END_REF] it is now widely accepted that dimers of larger PAHs like coronene form stacks eventually with a twisted or parallel displaced configuration, [START_REF] Zhao | A Prototype for Graphene Material Simulation; Structures and Interaction Potentials of Coronene Dimers[END_REF][START_REF] Rapacioli | Correction for dispersion and Coulombic interactions in molecular clusters with density functional derived methods: Application to polycyclic aromatic hydrocarbon clusters[END_REF][START_REF] Bartolomei | Modeling Coronene Nanostructures: Analytical Potential, Stable Configurations and Ab Initio Energies[END_REF][START_REF] Podeszwa | Interactions of graphene sheets deduced from properties of polycyclic aromatic hydrocarbons[END_REF][START_REF] Rapacioli | Stacked Clusters of Polycyclic Aromatic Hydrocarbon Molecules[END_REF][START_REF] Hernández-Rojas | Coarse-graining the structure of polycyclic aromatic hydrocarbons clusters[END_REF] with the exception of the work by Obolenski et al. [START_REF] Obolensky | Interplay of electrostatic and van der Waals forces in coronene dimer[END_REF] Let us also mention investigations concerned with gaps and electronic transport properties of these dimers. [START_REF] Adkins | Computed electronic structure of polynuclear aromatic hydrocarbon agglomerates[END_REF][START_REF] Fioressi | E↵ects of cluster formation on spectra of benzopyrene and benzopyrene[END_REF][START_REF] Sancho-García | Theoretical study of stability and chargetransport properties of coronene molecule and some of its halogenated derivatives: A path to ambipolar organic-based materials[END_REF][START_REF] Sancho-García | Theoretical study of stability and chargetransport properties of coronene molecule and some of its halogenated derivatives: A path to ambipolar organic-based materials[END_REF][START_REF] Sancho-García | Charge-transport properties of prototype molecular materials for organic electronics based on graphene nanoribbons[END_REF][START_REF] Pérez-Jiménez | Conductance Enhancement in NanographeneGold Junctions by Molecular -Stacking[END_REF][START_REF] Yoshida | Conducting Columns of Highly Symmetric Coronene, The Smallest Fragment of Graphene[END_REF][START_REF] Sanyal | E↵ect of Imide Functionalization on the Electronic, Optical, and Charge Transport Properties of Coronene: A Theoretical Study[END_REF] Addressing larger clusters has been generally carried out using force field schemes [START_REF] Totton | A quantitative study of the clustering of polycyclic aromatic hydrocarbons at high temperatures[END_REF]27,[START_REF] Herdman | Intermolecular Potential Calculations for Polynuclear Aromatic Hydrocarbon Clusters[END_REF][START_REF] Totton | Modelling the internal structure of nascent soot particles[END_REF][START_REF] Grančič | Evaluation of All-Atom Force Fields for Anthracene Crystal Growth[END_REF][START_REF] Fioressi | E↵ects of cluster formation on spectra of benzopyrene and benzopyrene[END_REF][START_REF] Bartolomei | Modeling Coronene Nanostructures: Analytical Potential, Stable Configurations and Ab Initio Energies[END_REF][START_REF] Piuzzi | Spectroscopy, dynamics and structures of jet formed anthracene clusters[END_REF][START_REF] Totton | A First Principles Development of a General Anisotropic Potential for Polycyclic Aromatic Hydrocarbons[END_REF][START_REF] Chakrabarti | A survey of the potential energy surface for the (benzene)13 cluster[END_REF][START_REF] Totton | A transferable electrostatic model for intermolecular interactions between polycyclic aromatic hydrocarbons[END_REF][START_REF] Easter | A Monte Carlo Study of Isomers and Structural Evolution in BenzeneCyclohexane Clusters: (C6H6)(C6H12)n, n = 37, 12[END_REF][START_REF] Bowal | Partitioning of polycyclic aromatic hydrocarbons in heterogeneous clusters[END_REF][START_REF] Pascazio | Simulating the morphology of clusters of polycyclic aromatic hydrocarbons: The influence of the intermolecular potential[END_REF][START_REF] Chen | Size-dependent melting of polycyclic aromatic hydrocarbon nano-clusters: A molecular dynamics study[END_REF][START_REF] Nakamura | Stability of multiply charged clusters of polycyclic aromatic hydrocarbons[END_REF][START_REF] Takeuchi | Structures, stability, and growth sequence patterns of small homoclusters of naphthalene, anthracene, phenanthrene, phenalene, naphthacene, and pyrene[END_REF][START_REF] Calvo | Possible Formation of Metastable PAH Dimers upon Pickup by Helium Droplets[END_REF] or coarse grains models [START_REF] Violi | Soot primary particle formation from multiscale coarse-grained molecular dynamics simulation[END_REF][START_REF] Heinemann | Coarse-grained electrostatic interactions of coronene: Towards the crystalline phase[END_REF][START_REF] Hernández-Rojas | Coarse-graining the structure of polycyclic aromatic hydrocarbons clusters[END_REF][START_REF] Hernández-Rojas | Dynamics and thermodynamics of the coronene octamer described by coarse-grained potentials[END_REF] which can be combined with extensive PES exploration schemes to derive thermodynamical properties [START_REF] Heinemann | Coarse-grained electrostatic interactions of coronene: Towards the crystalline phase[END_REF][START_REF] Chen | Size-dependent melting of polycyclic aromatic hydrocarbon nano-clusters: A molecular dynamics study[END_REF][START_REF] Hernández-Rojas | Dynamics and thermodynamics of the coronene octamer described by coarse-grained potentials[END_REF] or to find lowest energy structures. [START_REF] Bartolomei | Modeling Coronene Nanostructures: Analytical Potential, Stable Configurations and Ab Initio Energies[END_REF][START_REF] Hernández-Rojas | Coarse-graining the structure of polycyclic aromatic hydrocarbons clusters[END_REF][START_REF] Chakrabarti | A survey of the potential energy surface for the (benzene)13 cluster[END_REF][START_REF] Easter | A Monte Carlo Study of Isomers and Structural Evolution in BenzeneCyclohexane Clusters: (C6H6)(C6H12)n, n = 37, 12[END_REF][START_REF] Bowal | Partitioning of polycyclic aromatic hydrocarbons in heterogeneous clusters[END_REF][START_REF] Takeuchi | Structures, stability, and growth sequence patterns of small homoclusters of naphthalene, anthracene, phenanthrene, phenalene, naphthacene, and pyrene[END_REF] Most of these studies agree on the fact that, for coronenes, clusters growth proceed by stacking for small clusters followed by assemblies of stacks for larger systems. [START_REF] Herdman | Intermolecular Potential Calculations for Polynuclear Aromatic Hydrocarbon Clusters[END_REF][START_REF] Totton | Modelling the internal structure of nascent soot particles[END_REF][START_REF] Ricca | The Infrared Spectroscopy of Neutral Polycyclic Aromatic Hydrocarbon Clusters[END_REF][START_REF] Bartolomei | Modeling Coronene Nanostructures: Analytical Potential, Stable Configurations and Ab Initio Energies[END_REF][START_REF] Rapacioli | Stacked Clusters of Polycyclic Aromatic Hydrocarbon Molecules[END_REF][START_REF] Hernández-Rojas | Coarse-graining the structure of polycyclic aromatic hydrocarbons clusters[END_REF][START_REF] Pascazio | Simulating the morphology of clusters of polycyclic aromatic hydrocarbons: The influence of the intermolecular potential[END_REF][START_REF] Hernández-Rojas | Dynamics and thermodynamics of the coronene octamer described by coarse-grained potentials[END_REF] For smaller PAHs like anthracene or pyrene, assemblies of dimers or non stacked configuration compete for being the lowest energy structures. [START_REF] Ricca | The Infrared Spectroscopy of Neutral Polycyclic Aromatic Hydrocarbon Clusters[END_REF][START_REF] Takeuchi | Structures, stability, and growth sequence patterns of small homoclusters of naphthalene, anthracene, phenanthrene, phenalene, naphthacene, and pyrene[END_REF] One should mention that the previously reported studies deal with neutral systems and that very few is known about structural and energetic properties of singly charged clusters which is the main topic of the present paper, apart from the work by Bouvier et al. [START_REF] Bouvier | A Model Potential Approach to Charge Resonance Phenomena in Aromatic Cluster Ions[END_REF] who investigated, at a semi-empirical level, clusters of anthracene and naphtalene.

The treatment of singly charged clusters requires to use a method able to describe correctly the interplay between the intermolecular interaction with a computational cost allowing for dealing with large systems. This can be conveniently done via charge resonance models as an alternative to delocalized molecular orbitals (MO) schemes. Charge resonance in valence-bond theory concerns singly ionized systems made of equivalent or similar subunits where several charge-localized forms with identical or close energies are in competition.

Hence the wavefunctions are correctly described as combinations of those valence bond forms and not a single one. [START_REF] Shaik | A Chemist's Guide to Valence Bond Theory[END_REF][START_REF] Pauling | The nature of the Chemical Bond[END_REF] These models were largely used to describe charge localisation in singly ionsed rare gas clusters. [START_REF] Milko | Ab initio excited states calculations of Kr3+, probing semi-empirical modelling[END_REF] In the past decade, we have developed a model combining the Density Functional based Tight Binding with a small Configuration Interaction scheme (CI-DFTB) [START_REF] Rapacioli | Modeling Charge Resonance in Cationic Molecular Clusters: Combining DFT-Tight Binding with Configuration Interaction[END_REF][START_REF] Rapacioli | Modelling singly ionized coronene clusters[END_REF][100][101] which allows to treat properly the charge resonance between the di↵erent units. The scope of the present work is to determine the low energy structures of cationic clusters. The global exploration of PAH clusters, even with a rather e cient computational scheme such as CI-DFTB is not easy, due to the quite complex shape of the potential energy surface resulting of multiple degrees of freedom involving stacks or multi-stack structure, parallel displacements, twists along the axis or mutual orientation of stacks or of isolated molecules around stacks. We have developed a hierarchical strategy consisting for each size in (i) achieving global optimization within a frozen molecule framework and using a simplified however computationally very e cient parameterized version of the CI-DFTB scheme and (ii) finally relaxing a sample of the lowest-energy structures thus found at the actual CI-DFTB level. The methodology is introduced in the next section. We then report and discuss the structural, energetic and electronic properties of the most stable isomers found for neutral and cationic pyrene clusters.

Method Modelling neutral clusters DFTB

The Density Functional based Tight-Binding (DFTB) 102-107 is an approximated DFT scheme with a reduced computational cost particularly well suited to study large systems. In this approach, molecular orbitals { i } are expressed in a minimal basis set of atomic orbitals

(AO) {' µ } i = X µ c µi ' µ (1)
The DFT energy is then expanded with respect to the electronic density fluctuation ⇢ around a reference density ⇢ 0 . The expansion up to second order is known as the Self Consistent Charge (SCC-)DFTB. Neglecting three centers integrals allows then to write the SCC-DFTB energy as :

E = X ↵> E rep ↵ + occ X i X µ⌫ n i c µi c ⌫i H 0 µ⌫ + 1 2 X ↵ ↵ q ↵ q (2)
with E rep ↵ is a repulsion potential between atoms ↵ and , n i is the occupation of the orbital i, H 0 µ⌫ is the matrix element of the KS operator at the reference density expressed in the AO basis. The second order term depends on atomic Mulliken charges q ↵ and a function of biatomic distances ↵ (R ↵ ).

Several modifications have to be considered for a correct treatment of PAH clusters.

The dispersion, which is poorly described by traditionally used DFT functionals and consequently in the SCC-DFTB framework can be corrected by adding an empirical correction term [START_REF] Rapacioli | Correction for dispersion and Coulombic interactions in molecular clusters with density functional derived methods: Application to polycyclic aromatic hydrocarbon clusters[END_REF]108,109 and the intermolecular Coulomb potential can be improved by going beyond Mulliken charges definition to compute the second order term, in particular computing atomic charges based on the Charge Model 3 approach. 110 In the present work, we use the combination of SCC-DFTB using the mio-set 104 parameterization with a dispersion correction term and CM3 charges [START_REF] Rapacioli | Correction for dispersion and Coulombic interactions in molecular clusters with density functional derived methods: Application to polycyclic aromatic hydrocarbon clusters[END_REF] in consistency with our previous work on ionization potentials. [START_REF] Joblin | Size E↵ect in the Ionization Energy of PAH Clusters[END_REF] To accelerate the global exploration scheme, we have first developed a simplified version of DFTB which consists in removing the self consistent process, i.e. performing a single diagonalization step. This scheme is di↵erent from zeroth-order DFTB 102,103 as we still use the second order term of Eq. 2 but we only consider the first order perturbation with respect to the density (charges) computed at the previous geometry. We call this partial SCC scheme SCC-DFTB p , the idea being to perform several SCC-DFTB p steps between two SCC-DFTB steps involving the full self-consistent process.

Fast Force Field

In addition to SCC-DFTB and SCC-DFTB p , we have developed a Fast Force Field (FFF) scheme, which still uses the DFTB parameterisation. In this model, we initially determine the SCC-DFTB atomic charges q of an isolated pyrene monomer. In the following, the molecular units are kept frozen and clusters geometries are governed by the inter-molecular energy which is decomposed as:

E = E coul + E disp + E rep (3)
E coul is the coulomb energy depending on the precomputed atomic charges in the monomer and defined by :

E coul = X A>B X ↵, ↵ q↵ q (4)
where ↵ and are respectively atoms on molecules A and B. We use the same dispersion energy expression as in the SCC-DFTB scheme. The exchange-repulsion energy between molecules can be based on the orbital overlap, 111 we define as :

E rep = K 1 X A>B X µ2A,⌫2B S 2 µ⌫ (5)
where S µ⌫ is the DFTB atomic orbital overlap matrix. K 1 is an empirical parameter (0.85) determined to reproduce the SCC-DFTB results for specific geometries. Figure 1 illustrates that the main trends of SCC-DFTB are well reproduced by the FFF model along a generic path, consisting of geometric linear interpolations between characteristic structures.

Modelling cationic clusters CI-DFTB

In the case of singly ionized molecular clusters, the problem of self-interaction has to be carefully addressed because, in such systems, the wavefunction of the system is strongly multiconfigurational which is problematic for most of traditionally used DFT functionals 112,113 and consequently also at the DFTB level. To address such situations, we have developed the CI-DFTB scheme, which is and adaptation of the CI-DFT method 114 to the DFTB framework. In this approach, the global wavefunction of the system is expressed on a basis { i } corresponding to configurations where the charge is fully localized on one specific molecule

A. = N X A C A A (6)
Each charge localized configurations A is calculated by constrained-DFTB ( [START_REF] Rapacioli | Modeling Charge Resonance in Cationic Molecular Clusters: Combining DFT-Tight Binding with Configuration Interaction[END_REF][START_REF] Rapacioli | Modelling singly ionized coronene clusters[END_REF]115 ) which is the DFTB version of Constrained-DFT. 116 A small Configuration Interaction scheme can then be constructed

0 B B @ E 1 F 12 . F 21 E 2 . . . . 1 C C A 0 B B @ C 1 C 2 . 1 C C A = E 0 B B @ 1 T 12 . T 21 1 . . . . 1 C C A 0 B B @ C 1 C 2 . 1 C C A (7)
where the diagonal elements E i are the charge-localized configuration energies. The o↵diagonal elements F are coupling terms between charge-localized configurations and are computed following the formula of derived by Wu and Van Voorhis. 114,117 The overlap between the charge-localized configurations T is computed using the MO coe cients obtained from the charge localized calculation and the DFTB atomic overlap basis. Let us mention that the dispersion and charge corrections mentioned for the neutral DFTB calculations are also used and that the "SCC-" term is omitted in the name of CI-DFTB, for sake of conciseness, although the SCC process is actually involved. More details are given in the original CI-DFTB paper [START_REF] Rapacioli | Modeling Charge Resonance in Cationic Molecular Clusters: Combining DFT-Tight Binding with Configuration Interaction[END_REF] as well as benchmarks showing the ability of this scheme in reproducing ab initio binding energies of benzene and water dimers.

Let us finally mention that, in order to gain in computational e ciency and similarly to SCC-DFTB p scheme, the self-consistent process can be removed from the constrained-DFTB calculation, leading to the CI-DFTB p scheme.

CI-FFF

Finally, we have also developed a force-field version for the CI scheme, still solving the CI Eq. 7 but with a simpler scheme to compute the matrix elements, in the spirit of the approximations of the FFF model (Eq. 3). Configuration A is built using the MOs of an isolated cation for A and those of isolated neutrals for the other units. Consistently, the energies E A of the charge constrained forms are computed taking the frozen monomer charges to compute the coulomb energy from Eq. 4, namely the Coulomb energy of a configuration

A is calculated with the charges of an isolated cation molecule for unit A and those of a frozen neutral molecule for the other units.

Moreover, we have parameterized the CI-DFTB matrix elements according to the empirical following expressions

T AB = erf ⇣ K a 2 |h A | B i| K b 2 ⌘ (8)
F AB = ✓ E A + E B 2 + K 3 ◆ T AB (9)
where K a 2 , K b 2 and K 3 have been fitted to reproduce the CI-DFTB matrices elements over a sample of geometries (K a 2 = 10.546, K b 2 = 0.8 and K 3 = 0.0572). As can be seen in Figure 1, the CI-FFF model reproduces the main trends of the CI-DFTB energies along typical geometrical paths.

Global exploration

The currently most used global optimization techniques include: basin hopping and derived techniques, 118-121 simulated annealing 122 and developments concerned with ergodicity (multi-tempering, parallel-tempering in Molecular Dynamics or Monte Carlo schemes 123,124 ), genetic/evolutionary algorithms, 125 and particle swarm algorithms. [126][127][128][129] In this work, we essentially use a parallel temperiing Monte Carlo (PTMC) scheme completed by relaxation.

The detailed scheme consists in three steps: We first perform an 'All-Exchange Parallel Tempering Monte-Carlo' 130,131 exploration. In this step, the molecules are kept frozen (rigid-body approximation) and the potential energy is computed with the FFF model for neutrals and the CI-FFF model for cations, respectively. We have performed about 10 7 Monte Carlo steps for each temperature, the distribution of the latter following a geometrical suite between 10 and 500 K. At each step, the probability of attempting the all-exchange procedure is 1 %.

From these dynamics, we select up to 2000 lowest energy structures for a given geometrical pattern. The geometrical pattern is defined by the size of sub-stacks and the moment of inertia. In a second step, these selected structures are optimized locally through a conjugated gradient algorithm still using the rigid molecules approximation and making use of the potentials SCC-DFTB p and CI-DFTB p for neutrals and cations, respectively. The last step consists in a local optimization (conjugated gradient), relaxing all degrees of freedom and performed at the highest level of theory used in this work to compute the potential, namely SCC-DFTB and CI-DFTB for neutrals and cations, respectively. Obviously, an advantage of global exploration schemes as used here is the access not only to the most stable structure (presumably) but also to higher energy isomers.

Results

In this section, we present and discuss the most stable structures found for pyrene clusters up to 7 units for cations and neutrals. We define the nomenclature that will be used to characterize these clusters based on the number of units in each stack. For instance, the hexamer of Figure 2 (bottom left) is labelled 3/2/1 as it consists of two stacks of 3 and 2 units respectively plus one isolated monomer and the heptamer (bottom right) is labelled 3/2/2. The energies for evaporating a neutral monomer from a neutral cluster, a neutral monomer from a cation cluster and a cation monomer from a cation cluster are : where E 0 N and E + N are respectively the energy of a neutral and cationic cluster of size N . We also define the cohesive energy per molecule for each stable isomer

E 0/0 evap (N ) = (E 0 N E 0 N 1 E 0 1 ) (10)
E +/0 evap (N ) = (E + N E + N 1 E 0 1 ) (11)
E +/+ evap (N ) = (E + N E 0 N 1 E + 1 ) (12)
E 0 co/mol (N ) = (E 0 N NE 0 1 )/N (13) E + co/mol (N ) = (E + N (N 1)E 0 1 E + 1 )/N (14)

Neutral pyrene clusters

In Figures S1-S6 of the supporting information, we report the global minimum structure, the two or three most stable isomers that belong to the family of the global minimum, the most stable isomer for the following most stable families 1 and the most stable stack structure.

These isomers are in the following denoted by a number (the number of units) and a letter given following the energetic ordering. The most stable structures found for pyrene clusters (labelled with a "a") from 2 to 7 units are represented on Figure 2. It can be seen that the most stable dimers and trimers form stacks whereas larger clusters are made of adjacent stacks containing 1, 2 or 3 units. Clearly, the most stable structures correspond to compact forms, which is an expected trend to maximise the van der Waals interactions. We can also notice that, when possible, hydrogen atoms (slightly positive) lie pointing towards external carbon atoms (slightly negative) of another unit, as expected to maximise the stabilisation through Coulomb interactions. Whereas the best structures for the dimer and the trimer are found to be stacks, the 4-molecule stack (4g) is only the fifth family of the tetramer and lies at 0.187 eV above the lowest energy structure 4a that belongs to the 3/1 family. The structural excitation energies from the most stable structure to the other reported isomers are reported in Table 1. The excitation gap towards another family remains small for all sizes (if one excepts the special case of the dimer), the largest gap is obtained for the pentamer with a value of 56 meV. For each families, lot of stable structures generally exist with a quasi- [START_REF] Léger | Identification of the 'unidentified' IR emission features of interstellar dust?[END_REF] Only one isomer of the most stable family is reported in the case of the hexamer as the other isomers found in this family are less stable than the most stable isomers of other families degeneracy energy (almost 10-15 meV) due to the twisted quasi-degeneracy structures in the dimer/trimer stacks. In order to discuss the geometries in more details, we have calculated the distances between the centers of mass of the constituting molecules, the mutual twist angles around the axis in the case of stack subblocks and the dihedral angle of neighbouring molecules planes within the stacks. If we consider for reference the stack structures, the typical distances (between the centers of mass) are almost the same, namely 3.43 Å in the dimer 2a, 3.42/3.42 Å in the trimer 3a and 3.41/3.40/3.41 Å in the tetramer stack 4g. Also the mutual 1-2 and 2-3 twist angles between two neighboring molecules are around 68 degrees.

The main structure di↵erence between 3a and 3b is the clockwise or counter-clockwise twist angle (i.e. in 3a, clockwise/counter-clockwise angles lead to a twist angle close to 0 degree between the two external 1-3 molecules and in 3b, clockwise/clockwise angles result into a 1-3 twist angle close to 44 degrees). This twist isomerization of course also concerns subblocks stacks of larger size clusters producing energetically close-lying isomers. Isomers 5a, 5b and 5c all belong to the 3/1/1 family. Actually in 5a, the trimer stack is symmetric (bond lengths of 3.91/3.91 Å) while in 5b and 5c the trimer stack is non-symmetric (3.49/3.88 Å in 5b and 3.63/3.82 Å in 5c). Furthermore isomers 5b and 5c di↵er via the orientational order of the molecules around the axis in the trimer stack, namely respective twists of 2/84 degrees in 5c instead of 12/4.75 degrees in 5b. Hexamers 6b and 6c also belong to the same 3/3 family.

In 6b , both pairs of trimers stacks are non symmetric (3.98/3.55 Å and 3.71/3.52 Å) while they are symmetric in 6c (3.73/3.73 Å and 3.73/3.73 Å). Heptamers 7a, 7b and 7c also all belong to the 3/2/2 family. They di↵er either via twist angles in one of the dimer stack or via mutual locations of the stacks. The observation that several families are energetically competing for being the most stable at each size can also be seen from their cohesive energy per unit as reported in Figure 3 as for a given size, several famillies compete for having the largest cohesive energy. In this plot, lines connect clusters which only di↵er by changing the number of units in the largest central stack, for instance the red N-1/1 (resp. blue N-2/1/1) line connects the most stables structures with a single (resp. two) side molecule, namely the structures 1/1, 2/2, 3/1, 4/1 (resp. 1/1/1, 2/1/1, 3/1/1, 4/1/1/1). It can be seen that the shapes of these curves look similar allways increasing with the number of units, the increase being larger for the smaller number of units. Let us consider the N-1/1-red curve. It's cohesive energy increases strongly between 2 and 3 units in the stack (clusters of size 3 and 4). This familly is the most stable at for tetramer (largest cohesive energy) and the familly N-2/1/1 is not competitive at all for this size. However, the cohesive energy increase between tetramer and pentamer is stronger for the familly N-2/1/1 as it corresponds to central stack passing from 2 to 3 units than for the N-1/1/1 familly where the number of units in the main stack increases from 3 to 4 units. To generalize, we can say that, although the trends are similar for these families, their relative positions depend of course on the number of units removed from the main stack and how the latter are positioned around the main stack. The crossing between these curves govern the stability switches between the di↵erent families as a function of size. The intermolecular distance in the neutral dimer 2a is 3.42 Å. One can notice that this distance only slightly varies when the twist angle is changed: in isomer 2b, characterized by a twist angle of 50 degrees, the intermolecular distance is 3.34 Å. The intermolecular distances in the trimer stack 3a are the same as in the dimer, which is consistent with the fact that dispersion-repulsion is the most important contribution to the bonding. In larger clusters, the main stacks are significantly deformed and bent, due to the interactions with the side molecules. Correlatively, the intermolecular distances within the stacks increase, namely 3.91/3.91 Å, 3.49/3.89 Å, 3.99/3.99 Å in the trimer stacks of the tetramer, pentamer and hexamer respectively. The same is observed in the dimer stacks of the hexamer (3.79 Å) and of the heptamer (two dimer stacks at 3.88 Å). The twists become irregular, with a juxtaposition of large and small twists angles for molecules along a same stack. Also and correlatively with the axis bending, the dihedral angles become significant.

N name family E (eV) E co/mol E (Ha)
The evaporation of a monomer is always the energetically lowest fragmentation channel.

For instance in the lowest energy hexamer 6a which is composed of a trimer stack, a dimer stack and a single molecule, the channel corresponding to the evaporation of a dimer is located 0.34 eV above that corresponding to a monomer evaporation. Likewise, the channel for evaporation of a dimer from the 3/2/2 heptamer 7a is located 0.38 eV above that corresponding to the evaporation of a monomer. As can be seen from Figure 4 and Table 2, the evaporation energy is an increasing value with size, starting at about 0.5 eV for the dimer, which tends to stabilize around 0.9 eV for the largest clusters. Within a pair-additive approximation of the dispersion and Coulomb interactions, the attachment of one unit to the clusters is expected to grow with the number of units inside the cluster. On the other hand, as these interactions decrease with the distance, especially the van der Waals interaction, the evaporation energy in large clusters is expected to reach a limit which would correspond to the evaporation of a monomer from a bulk surface.

Let us notice that in our previous work 80 using a force field scheme, the single stack structure was found to be the most stable up to 6 units instead of 3 as in the present work.

Interestingly, by comparing some typical structures (single stacks, intertwined, herringbone, ...), Ricca et al. [START_REF] Ricca | The Infrared Spectroscopy of Neutral Polycyclic Aromatic Hydrocarbon Clusters[END_REF] 1 for all computed sizes, but the evolution with size is similar. We however notice that our value for the cohesive energy of the dimer is in reasonable agreement (0.248 vs 0.215 eV/mol) with SAPT calculations. 132 From a global exploration with a force field scheme, Takeuchi et al. [START_REF] Takeuchi | Structures, stability, and growth sequence patterns of small homoclusters of naphthalene, anthracene, phenanthrene, phenalene, naphthacene, and pyrene[END_REF] reported for pyrene clusters made of subblocks arranged in an amorphous way instead of herringbone patterns. This is in line with the present results except that the subblocks contains two or three units in our case and not more that two in the former work.

E co/mol (eV/Molecule) N N N-1/1 N-2/2 N-2/1/1 N-3/3 N-3/2/1 N-3/1/1/1 N-4/2/2 N-4/3/1 N-4/2/1/1 N-5/2/
E co/mol (eV/Molecule) N N N-1/1 N-2/2 N-2/1/1 N-3/3 N-3/2/1 N-3/1/1/1 N-4/2/2 N-4/3/1 N-4/2/1/1 N-5/2/2/1

Cationic pyrene clusters

Following the same way for labelling and presenting the results as the one used for neutral clusters, the most stable isomers are presented in Figure 5 whereas all the discussed isomers are presented in Figures S7-S12 of Supporting Information. At first glance, the most stable structures found for cationic pyrene clusters (Figure 5) look similar to those obtained for the neutral clusters, at least up to the tetramer whereas di↵erences are observed for larger clusters, in particular the hexamer and the heptamer. The latter clusters (6a and 7a) are made of a central stack consisting of four molecules with isolated side molecules whereas in their neutral counterpart, the main stack is smaller (three units) and the outer molecules tend to form dimer blocks. Beyond this purely geometric aspect, can examines the distribution of the charge over the various molecular units which is reported in Table 3 for each cluster size and the various isomers. It appears that, whatever the number of units, the charge is mostly localized on less than four units of the central stack. In cluster 7a, the charge is mostly carried by two units (a charged dimer core with charges 0.43/0.40) whereas it is delocalized over three units for the trimer (charges 0.27/0.46/0.26) and the tetramer (charged trimer core with charges 0.24/0.49/0.25), the central unit of the core carrying twice more charge than each of the two edge unit. The case of the hexamer and the pentamer occurs as intermediate between a charged dimer and trimer core. All clusters beyond N=3 can therefore be seen as a charged core of two or three units with the other units organizing to form a solvation shell around the charge instead of maximizing the dispersion term that would favor, as in the case 2a 3a 4a 5a 6a 7a

Figure 5: Most stable structures of cationic pyrene clusters.

of neutral clusters, a more compact structure with external dimer blocks. Let us notice that Bouvier et al. [START_REF] Bouvier | A Model Potential Approach to Charge Resonance Phenomena in Aromatic Cluster Ions[END_REF] already concluded from an empirical valence bond model investigation that the charge in cationic clusters of benzene, naphtalene and anthracene is mostly carried by two to four units. Again, as in the case of the neutral clusters, quasi-degeneracy situations occur between families and within a same family. For instance isomers 5b and 5c geometrically di↵er via the orientational order within the tetramer stack core. The intermolecular distances are 3.12/3.14/3.94 Å in 5b and 3.25/3.07/3.58 Å in 5c, while the orientational twist angles are 60/61/3 degrees in 5b and 3/59/37 degrees in 5c.

Despite global structural similarities with neutrals mentioned above, the details of the geometries show significant quantitative di↵erences. We first discuss as reference the intermolecular bond lengths of the single stacks clusters (not the lowest structures for N > 3). In dimer 2a, the intermolecular distance is reduced from 3.43 in the neutrals species to 3.04 Å in the charged species while the twist angle is also reduced from 67 to 55 degrees. This bond shortening can also be evidenced in the stack trimer 3a, where the bond length decreases from 3.42/3.42 Å in the neutral to 3.09/3.09 Å in the charged species. This is clearly associated with a stronger bonding due to charge resonance. For the larger sizes, as discussed previously, single stack clusters are not the lowest structures and the intermolecular bond lengths within the central the stack present noticeable variations along the axis, namely The intermolecular separation tends to increase for the largest clusters, namely 3.28/3.26/3.22/3.26/3.28 Å in the stack hexamer with a vanishing variation between the edge and the center, and 3.24/3.24/3.24/3.24/3.24/3.24 Å in the stack heptamer. All single stack clusters present systematic twist angles in the range 50-56 degrees. Correlatively in these single stack and rather regular structures, the charge distribution approximately follows a square cosine dependence along the stack with a charge maximum about the center of the stack and decreasing values down to the the edges, as would be obtained in a simple topological Hückel model or even in a one-dimensional square well. [START_REF] Rapacioli | Modelling singly ionized coronene clusters[END_REF] We now discuss in more details the geometries of the lowest energy structures for N > 3. Obviously the increased longitudinal compactness is also observed in the charged dimer or trimer cores of the larger clusters, namely bond lengths of 3.10/3.10 Å in the trimer core of 4a, 3.31/3.33 Å for 5a. The intermolecular distances in the tetramer core of the cation hexamer 6a are 3.33/3.18/3.25 Å, quite similar to those of the tetramer core in the cation heptamer 7a, namely 3.45/3.11/3.30 Å. One distance at a stack edge is thus somewhat larger and this is consistent with the fact that the charge is dominantly carried by the two central molecules of the stack as already mentioned. Finally, one can examine the closest distances between the side neutral molecules to the closest ones in the core stack. These core-neutral distances range is between 5.63 Å (nearest core-neutral intermolecular distance in isomer 4a for instance) and 6.55 Å depending on the cluster sizes and the isomers. The dihedral angles between the molecules in the central stack are smaller and the stack axes for cations with size N > 4 have less bending than in the neutral equivalent geometries.

Evaporation of a neutral monomer from a cluster cation is always more favorable than evaporation of a cationic one (see In the following we discuss the neutral monomer evaporation channel only, dissociation energies for other channels can be easily computed from the absolute energy values given in tables 1 and 3. Contrary to the case of neutrals, the evaporation energy (Figure 4) is not monotonous anymore. For the dimer, it corresponds to a binding energy almost twice larger than for neutral dimer and this is due to the strong stabilisation arising from charge resonance. In the trimer, the less bonded molecules are the two external molecules, each of them being less involved in the charge transfer (they only carry one fourth of the total charge). This explains why removing an external unit is less costing in term of charge resonance stabilisation energy. For the larger systems, the less bonded (external molecules) are even less involved in the charge resonance stabilisation (carrying less that 2 % of the total charge), this should make them easier to remove but, at the same time, dispersion interactions increase with the cluster size. This explains why the evaporation energy remains almost constant for clusters between four to seven units, with values close to those observed for neutrals.

The first structural excitation energy toward another family (see table 3) can vary between 0.04 eV in the case of the pentamer or the heptamer up to 0.15 eV for the tetramer, while it was shown to be smaller than 0.06 eV in all cases for neutrals (except the dimer).

The most stable isomers of the di↵erent families are spread over a larger energy range than what is observed for neutrals. For instance, in the case of the neutral heptamer, 10 families present isomers with energies lower than 0.2 eV above the most stable structure whereas only four families could be identified within the same energy range for cations. Indeed, an energy range of 0.49 eV is necessary to observe at least 9 di↵erent families for the cationic heptamer. The cohesive energies are reported in Figure 3. For the stack family, we can observe the strongest stabilisation when the central stack grows from two to three units due to better charge resonance e↵ect. Adding more units to the stack brings only van der Waals and Coulomb stabilisation energy, resulting in a decrease of the cohesive energy per unit from the tetramer. This increase/decrease behavior is also visible for the N-1/1, N-2/1/1 and N-2/1/1 families allways with a cohesive energy maximumobtained when the central stack contains 4 units. For the other families, only the increase is visible, the deacrease occuring probably for sizes larger that the ones investigated here.

Ionization potentials

Figure 6 and table 4 presents the adiabatic ionization potentials obtained as di↵erences between the energies of the most stable structures presented in Figures 2 and5. It presents a decrease up the tetramer which comes from the fact that, in the cation, the charge resonance stabilisation energy increases with size. However, and as already discussed previously, the charge resonance mostly concerns the small dimer/trimer core and, as a consequence, the ionization potential is expected to converged fast. The limit should be the energy required to ionize the central core, which is approached for the largest clusters. Of course, the polarization of the molecules surrounding the core is expected to continue decreasing the ionization potential but it seems to be a smaller contribution here.

Conclusions

The most stable structures for neutral and cationic pyrene clusters containing up to seven units have been searched through a global exploration scheme combining Parallel Tempering Monte Carlo algorithm followed by local Conjugated Gradient quenches on the potential energy surface. In the final stage, the potential energy is computed for neutral clusters using the SCC-DFTB model and for cations using the CI-DFTB model. The latter accounts for charge resonance by performing explicit interaction between charged localized configurations obtained at the constrained SCC-DFTB level. We have also developed faster approaches which are used in the global exploration procedure : the Fast Force Field (FFF) scheme for neutrals, which consists to summing dispersion, coulomb, repulsion terms, computed from the DFTB parameters with some empirical formulae. For cations, we developed the CI-FFF scheme, which is build on the same CI scheme as in the CI-DFTB method but the energies of the charged localized configurations are computed at the FFF level and an empirical formula is derived to compute the coupling terms from the overlap of the unperturbed HOMO of the sub-units. FFF and CI-FFF are shown to behave similarly to DFTB and CI-DFTB on some selected PES pathways.

The most stable structures obtained for pyrene neutral clusters are made of assemblies of small sub-units containing less than four units. These sub-units are arranged to form compacts structures, maximizing the pair-pair Coulomb and dispersion interactions. In cations, the charge remains localized on a dimer or trimer stack core and the other molecules remain essentially neutral and arrange around this core either at the top of the stack or on the sides, initiating a first solvation shell not yet completed at N=7. We have analyzed the correlation between the intermolecular distances and the localization of the charge in the cations. We have shown that in the dimer and the trimer, despite a similar stack arrangements for neutrals and cations, intermolecular distances are reduced respectively to those in the neutrals, associated with the charge resonance stabilization. The same distance contraction a↵ects the core parts of the larger clusters investigated in the present work.

The energetic analysis shows that the cohesive energy per molecule increases regularly from 0.25 eV to 0.61 eV in the neutrals, from 0.54 eV to 0.74 ev for the cations. For all clusters, neutrals and cations, evaporation of a neutral monomer is the lowest channel. The energy required to evaporate a molecule from a neutral cluster increases with the number of units, as it essentially consists in atomic pairwise additive contributions. Although similar energies for neutrals and cations are obtained for the largest clusters investigated, the situation is di↵erent for small clusters as removing a molecule from the charged core results in weakening the charge resonance contribution to the stabilization of the cation. As a consequence, the pyrene dimer appears to be the most stable cluster with respect to evapo-ration. After reaching a minimum, the evaporation energy of cations increases again, due to increased stabilization via polarization of a larger number of molecules by the charged core.

Finally, the adiabatic ionization potential appears to decrease smoothly with size whereas the vertical ionization potential may present size e↵ects (irregularity at N=6) resulting from significant di↵erences in the relaxed geometries of neutrals clusters and those of cations.

Obviously, a strong quasi-degeneracy is observed in the competition for the lowest energy structures. The present work has been achieved in the context of DFTB, which is subject to various uncertainties in describing molecular clusters. The present version is an improvement to the basic scheme since dispersion is added and self interaction problems are circumvented in the CI-DFTB scheme. Nevertheless DFTB and CI-DFTB remain parametric methods and we hope that the present work will stimulate further investigations with higher status methods such as double hybrid DFT or even wavefunction methods, which have up to now only have applied to smaller PAH systems, essentially neutral dimers.

Finally, we also hope that the energetic data determined in the present work and in particular will be of interest for various groups working on PAH chemistry in astrochemistry or atmospheric chemistry context and also in laboratory experiments. Chemistry of PAH clusters could be a step in the stabilization of carbonaceous species either via bottom-up formation of larger systems or top-down dissociative processes. Finally clusters might provide an alternative or complementary source to large single PAH molecules with more than several tens carbon atoms in the explanation of some of the unindentified IR bands.

Figure 1 :

 1 Figure 1: Comparison of FFF and SCC-DFTB (neutral dimer, left panel) and CI-FFF and CI-DFTB (cation dimer, right panel) potential energies along specified pathways.

Figure 2 :

 2 Figure 2: Most stable structures for neutral pyrene clusters.

  also concluded from DFT B3LYP+D2 calculations that the single stack is the most stable structure up to three units, suggesting structures made of dimer block assemblies. The cohesive energies per molecule reported by Ricca et al. are systematically

Figure 3 :

 3 Figure 3: Cohesive energy per molecule for of the various structural families of neutral (left) and cationic (right) pyrene clusters. The family labelling explained in the text and is the same as used in tables 1 and 3

Figure 4 :

 4 Figure 4: Evaporation energy of a neutral monomer from neutral and cationic clusters.

Table 1 :

 1 Properties of neutral pyrene clusters (C

	1	a	1	-	-	-30.98132298
	2	a	2	-	0.248	-61.98086275
	2	b	2	0.005	0.246	-61.98069401
	2	c	1/1 (T-shape)	0.197	0.150	-61.97363821
	3	a	3	-	0.339	-92.98131455
	3	b	3	0.001	0.338	-92.98126463
	3	c	2/1	0.009	0.336	-92.98098180
	3	d	1/1/1	0.108	0.303	-92.97733323
	4	a	3/1	-	0.431 -123.98871559
	4	b	3/1	0.006	0.430 -123.98850709
	4	c	3/1	0.012	0.429 -123.98828068
	4	d	2/1/1	0.042	0.421 -123.98718280
	4	e	1/1/1/1	0.043	0.421 -123.98711967
	4	f	2/2	0.063	0.416 -123.98639958
	4	g	4	0.187	0.385 -123.98183460
	5	a	3/1/1	-	0.513 -155.00083290
	5	b	3/1/1	0.002	0.512 -155.00075753

Table 2 :

 2 Evaporation energy of a neutral monomer from neutral clusters (E

	0/0 evap ),

Table 3 :

 3 3.15/3.18/3.15 Å in the tetramer and 3.16/3.25/3.25/3.16 Å in the pentamer. Properties of pyrene cluster cations (C 16 H 10 ) + N : isomerization energy E, cohesive energy per molecule E co/mol , charge distribution over the constituting monomer units and absolute CI-DFTB energy E in Hartree.

	N label	family	E (eV) E	co/mol	charge distribution (%)	E (Ha)
	1	a	1	-		-	100	-30.70613754
	2	a	2	-	0.538	50/50	-61.72702634
	2	b	1/1 (T-shape)	0.513	0.282	88/12	-61.70818114
	2	c	1/1 (coplanar)	0.706	0.185	50/50	-61.70108551
	3	a	3	-	0.622	46/27/26	-92.73740273
	3	b	3	0.004	0.621	48/26/26	-92.73727340
	3	c	2/1	0.058	0.603	49/48/3	-92.73528438
	3	d	1/1/1	0.462	0.468	56/40/5	-92.72041806
	4	a	3/1	-	0.687	49/25/24/2	-123.75107186

Table 2

 2

	). Neutral monomer evaporation is actually the

Table 4 :

 4 Vertical (VIP) and Adiabatic (AIP) Ionization Potentials of pyrene clusters .Figure6and table 4 also report the vertical ionization potentials (i.e. di↵erence between neutral and cation energies computed both the neutral most stable structure). Although the trend is similar with that of adiabatic IPs, we can clearly see non continuous size e↵ects around the hexamer. Actually, the energetic contribution driving the geometries of the most stable neutrals clusters lead to structures which may or may not be favorable to charge resonance when the cluster is ionized or in other word, how large the relaxation in the ion is vs in the neutral. Interestingly, considering the ions at the geometries of the neutral, although the pentamer, hexamer and heptamer clusters present a main stack of three units, only the

	N VIP (eV) AIP (eV)
	1	7.55	7.49
	2	7.19	6.91
	3	6.98	6.64
	4	6.85	6.47
	5	6.79	6.44
	6	6.84	6.44
	7	6.67	6.51

whereas in the case of the hexamer, the charge is carried by only two units with a very unequal distribution (0.68/0.23). The energetic contribution in neutral clusters lead to a geometry not favorable to charge delocalization in the case of the the hexamer. Let us mention that, in agreement with our calculations, this size e↵ect has been shown to be clearly present in experimental measurements of pyrene IP's.

[START_REF] Joblin | Size E↵ect in the Ionization Energy of PAH Clusters[END_REF]

Acknowledgment

We acknowledge the financial support of the Agence Nationale de la Recherche through the GASPARIM project Gas-phase PAH research for the interstellar medium (ANR-10-BLAN-0501), the computing ressources provided by the CALMIP HPC mesocenter for allocation P0059 and GDR CNRS EMIE. L.D. thanks the ERC for support under grant ERC-2013-Syg-610256-NANOCOSMOS and EUR grant NanoX n 080876 in the framework of the "Programme des Investissements dAvenir".

Supporting Information

The supporting information presents all the structures of neutral and cationic pyrene clusters discussed in this paper.