Ljiljana Stojanović

Saadullah G Aziz

Rifaat H Hilal

Felix Plasser

Thomas A Niehaus

Mario Barbatti

Nonadiabatic dynamics of cycloparaphenylenes with TD-DFTB surface hopping

 Surface hopping with a more robust description of the excited states should be obtained via linearresponse TD-DFTB. Nonadiabatic dynamics with linear-response TD-DFTB, which has also been implemented by Mitrić et al. 4b based on an earlier version of DFTB, is the basis of our developments here.

Introduction

Nonadiabatic effects influence the underlying relaxation mechanisms of a range of processes occurring after photoexcitation in systems of all sizes, spanning from single molecules to solid-state materials. Some of these nonadiabatic processes such as charge and energy transfer, electron-hole separation, singlet fission, etc. 1 have attracted a lot of attention due to their potential to aid the development of photonic and photovoltaic materials. Therefore, a number of research groups have been aiming to develop and extend suitable theoretical algorithms for simulation of nonadiabatic processes in extended systems. [START_REF] Akimov | The PYXAID Program for Non-Adiabatic Molecular Dynamics in Condensed Matter Systems[END_REF] Since the dimensionality of extended systems is a major limiting factor in their excited-state dynamics simulations, the methods that have been used for this purpose usually rely on mixed quantum-classical nonadiabatic dynamics approaches, such as the trajectory surface hopping in its different variants. 2a, 2b, 3 These methods are often associated with fast semiempirical algorithms for electronic structure computations. 2b, 4 The density functional tight binding (DFTB) [START_REF] Elstner | Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] method bridges the gap between first-principles electronic structure methods and empirical tight-binding approaches. DFTB can be seen as an approximation of density functional theory (DFT) at the generalized-gradient approximation (GGA) level. The basic algorithms implemented in the DFTB+ code, [START_REF] Aradi | a Sparse Matrix-Based Implementation of the DFTB Method[END_REF] which we use here, are the non-self-consistent charge (non-SCC) [START_REF] Porezag | Construction of tightbinding-like potentials on the basis of density-functional theory: Application to carbon[END_REF] and SCC-DFTB (DFTB2 and DFTB3) [START_REF] Elstner | Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] for the electronic ground state. The treatment of dispersion interactions (+D) is possible in the DFTB+ code via Lennard-Jones potentials, [START_REF] Rappe | UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations[END_REF] with the Slater-Kirkwood polarizable atomic model, [START_REF] Elstner | Hydrogen bonding and stacking interactions of nucleic acid base pairs: A density-functional-theory based treatment[END_REF] and on the level of the DFT-D3 method. [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF] Electronic excitations can be calculated with linear-response time-dependent approach (LR TD-DFTB). [START_REF] Niehaus | Tightbinding approach to time-dependent density-functional response theory[END_REF] A recent extension of the method now also allows for computations using long-range corrected (LC) functionals, 12 which provide a better description of state localization and charge transfer excitations. The low computational cost of these algorithms renders levels up to LR TD-(LC)-DFTB3 (+D3) suitable for simulations of excited states in extended systems.

In this work, we implemented an interface between DFTB+ [START_REF] Aradi | a Sparse Matrix-Based Implementation of the DFTB Method[END_REF] and Newton-X, [START_REF] Barbatti | NEWTON-X: A package for Newtonian dynamics close to the crossing seam[END_REF] allowing to run decoherence-corrected fewest switches trajectory surface hopping (FSSH) with the TD-DFTB method [START_REF] Niehaus | Tightbinding approach to time-dependent density-functional response theory[END_REF] available in the recently released 1.3 version of the DFTB+ code. There are reports of few other implementations of nonadiabatic dynamics with mixed quantum-classical methods based on DFTB.

Jakowski and Morokuma have implemented Ehrenfest dynamics based on Liouville-von Neumann

Our surface hopping based on linear-response TD-DFTB implementation relies on a general interface between the newest versions of the DFTB+ and Newton-X codes, which are both publicly available free of charge. It enables diverse new features to improve the description of ground and excited state energies, such as the DFTB third-order expansion (DFTB3 model) [START_REF] Yang; Yu | Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method: Third-Order Expansion of the Density Functional Theory Total Energy and Introduction of a Modified Effective Coulomb Interaction[END_REF] and onsite corrections for TD-DFTB; [START_REF] Domínguez | Extensions of the Time-Dependent Density Functional Based Tight-Binding Approach[END_REF] to accelerate the dynamics, such as fast nonadiabatic coupling calculations based on the recently proposed orbital derivative method; [START_REF] Ryabinkin | Fast Numerical Evaluation of Time-Derivative Nonadiabatic Couplings for Mixed Quantum-Classical Methods[END_REF] and to enhance the data analysis, via a direct module for advanced transition density analysis based on the program TheoDORE. [START_REF] Plasser | TheoDORE: A package for theoretical density, orbital relaxation, and exciton analysis[END_REF] Moreover, this interface is ready to profit from analytical energy gradients for TD-DFTB based on long-range corrected functionals, 12 which are currently under development.

The FSSH/TD-DFTB dynamics (on the GGA level of functional) has a potential for the simulation of processes where the charge-transfer states are not involved in dynamics (e.g. relaxation processes of small molecules in solutions, where both solute and solvent molecules are treated with TD-DFTB, see e.g. Ref. [START_REF] Mitrić | Nonadiabatic Dynamics within Time-Dependent Density Functional Tight Binding Method[END_REF] ), or in the cases where the occurring charge-transfer states are of short-distance (like for instance in the process of hot electron injection in photovoltaic materials, see e.g. Ref. [START_REF] Hartland | What's so Hot about Electrons in Metal Nanoparticles?[END_REF] ). Even in cases where TD-DFTB is not expected to work well, mainly due to artifacts in the description of charge-transfer states, a multilayer approach may be devised, [START_REF] Grimme | A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z= 1-86)[END_REF] where dynamics is propagated with TD-DFTB and single points are recalculated at higher level. Naturally, in such a case, time-dependent properties (population lifetimes, for instance) may not be adequately described, while statistical distributions of the excited-state in equilibrium may still deliver relevant information.

As a test case for this new implementation, we choose to study the nonadiabatic dynamics of two cycloparaphenylenes, [8]CPP and [10]CPP. Cycloparaphenylenes are hoop-shaped conjugated molecules composed of repeating phenylene units bonded in para positions to form a ring. They have been synthesized by several experimental groups 24 and since then studied experimentally [START_REF] Nishihara | Exciton recombination dynamics in nanoring cycloparaphenylenes[END_REF] and theoretically [START_REF] Camacho | Origin of the size-dependent fluorescence blueshift in [n]cycloparaphenylenes[END_REF] in detail. They received a lot of attention because of their potential use as precursors in the synthesis of armchair carbon nanotubes [START_REF] Omachi | Initiation of carbon nanotube growth by well-defined carbon nanorings[END_REF] and as efficient emitters of tunable wavelengths. [START_REF] Segawa | Combined experimental and theoretical studies on the photophysical properties of cycloparaphenylenes[END_REF] CPPs show unusual optical properties, markedly distinct from the linear paraphenylenes, [START_REF] Lukes | Torsional potentials and full-dimensional simulation of electronic absorption and fluorescence spectra of para-phenylene oligomers using the semiempirical self-consistent charge density-functional tight binding approach[END_REF] such as size-independent absorption wavelengths, blue shifting of emission peaks with an increase in the ring size 25b , double-peaked emission spectra, 25b among others. CPP excitations are characterized by long relaxation lifetimes in the nanosecond domain. 25d, 30 Small CPPs are also interesting due of their large quantum yields for singlet oxygen production in the reaction with triplet oxygen. 25d, 30 As we shall discuss, nonadiabatic dynamics with LR TD-DFTB explore states with somewhat large CT character. This is not surprising given that conventional DFTB is parameterized for PBE, where exactly this same type of problem is expected to occur. [START_REF] Dreuw | Single-Reference ab Initio Methods for the Calculation of Excited States of Large Molecules[END_REF] Recent results for (LR) TD-LC-DFTB make clear that this level provides correct energy character, fixing the main deficiencies of conventional DFTB. 12b Since long-range corrected functionals are indispensible for the correct description of charge-transfer states, but the analytical gradients for long-range corrected TD-DFTB are under development at the moment, we resort to a multilayer dynamics scheme, where the conventional simulations are performed with TD-DFTB and the TD-LC-DFTB computations of the excited states populated during dynamics are performed a posteriori. This approximation precludes us of drawing out conclusions about states' lifetimes, because of the mismatch between TD-DFTB and TD-LC-DFTB gradients. Still, plenty of important insights about relaxation mechanisms could be deduced, primarily because the long-enough dynamics provides mapping of the excited states hypersurfaces, particularly the S1 hypersurface, and thus fluorescence and internal conversion processes can be studied, as well as the properties of occurring excitons.

The Newton-X/DFTB+/TheoDORE interface that we present in this paper allows also the analysis of transition density matrices of the states populated during dynamics. We utilized this feature to perform the statistical analysis of excitation energies, charge transfer, density delocalization, and multiconfigurational character for a large ensemble of points generated during the dynamics, based on different approximations of the S1 state wave functions, (LR) TD-(LC)-DFTB and SD-(LC)-DFTB. This comparison allows us to derive a comprehensive picture about what kinds of deviations in the electronic densities are caused by an exploration of the configurational space by GGA functionals and single determinant approximations. This analysis confirmed that the full LR TD-DFTB representation of the wave function is a better approximation for the excited states than the single-determinant one, often applied in the simulations of nano-sized systems. 4a This paper is organized as follows. In the section Method Description, we briefly review the fewest switches surface hopping and TD-DFTB methods, providing details on the interfaces between DFTB+, Newton-X, and TheoDORE codes. Computational Details contains the technical information about the methods used for the actual calculations. In the section Results and discussion, we first discuss the reliability of the description of the excited states of the [8]CPP and [10]CPP on the TD-DFTB level, comparing the vertical excitations with TD-DFT results calculated with three different functionals. Next, we discuss the TD-DFTB excited state dynamics of cycloparaphenylenes. Finally, we statistically analyze the quality of the TD-DFTB electronic structure information produced during the dynamics and evaluate the SD-(LC)-DFTB levels.

Method Description

Fewest switches trajectory surface hopping (FSSH)

FSSH is a mixed quantum-classical approach to quantum dynamics, which approximates the wave packet motion by a swarm of independent nuclear trajectories, each one evolving classically on a single electronic state. [START_REF] Tully | Molecular-Dynamics with Electronic-Transitions[END_REF] Nonadiabatic transitions between states are introduced by a stochastic algorithm which may instantaneously switch the surface on which the trajectory is propagated, based on evaluating of hopping probabilities. [START_REF] Tully | Molecular-Dynamics with Electronic-Transitions[END_REF] The hopping probabilities in FSSH are obtained by expanding the time-dependent electronic wave function as

(1)
is the adiabatic electronic wave function for state K, which is a function of the electron coordinates r and parametrically dependent on the classical nuclear coordinate R.

are complex time-dependent expansion coefficients. The nuclear trajectories are propagated by integrating Newton's equations on the active state M.

Simultaneously to the propagation of Newton's equations, a local-approximation of the time-dependent electronic Schrödinger equation (TDSE) is integrated to obtain the expansion coefficients, . In the adiabatic representation, it is given by (2) where is the energy of the electronic state J and the bra-ket in the second term of the right side represents the time-derivative nonadiabatic coupling (NACs) between states J and K (

.

K K K t c t t y Y = å r R r R ( ) (
)

; K t y r R ( ) K c t ( ) t R ( ) K c t , J J J JK K K dc i Ec i c dt s = -å J E 6 (3)
The coefficients are further corrected for decoherence effects with the simplified decay of mixing method (SDM) [START_REF] Granucci | Critical Appraisal of the Fewest Switches Algorithm for Surface Hopping[END_REF] according to (4) In these equations, M is the active state and the decoherence time is given by the phenomenological equation (5) where is the potential energy of state I, is the classical kinetic energy of the nuclei, and is a parameter whose recommended value is 0.1 Hartree.

The FSSH probability [START_REF] Tully | Molecular-Dynamics with Electronic-Transitions[END_REF] between the active state M and another state K at each time step is computed as (6) where Dt is the time step for the integration of the TDSE. Thus, once the initial conditions (positions and velocities) are created, to perform the FSSH simulations, one must solve the electronic TDSE at each nuclear position given by Newton's equations to obtain the excitation energies, gradients, and to evaluate time-dependent NACs. These quantities are interpolated between two classical time steps to be used in the integration of the TDSE.

Our FSSH/TD-DFTB implementation profits of the algorithmic infrastructure already implemented in Newton-X. In the standard way, time-dependent NACs are computed by finite differences method as proposed by Hammes-Schiffer and Tully. [START_REF] Hammes-Schiffer | Proton transfer in solution: Molecular dynamics with quantum transitions[END_REF] In this approach, TD-NACs are given in terms of excited-state wave function overlaps between two subsequent nuclear time steps,

. The wavefunction overlaps SJK can be evaluated on a basis of Slater determinants (determinant derivative method). [START_REF] Barbatti | Non-Adiabatic Dynamics of Pyrrole: Dependence of Deactivation Mechanisms on the Excitation Energy[END_REF] The TD-NACs can be alternatively evaluated with the recently implemented orbital derivative method [START_REF] Stojanović | New Insights into the State Trapping of UV-Excited Thymine[END_REF] developed by Izmaylov's group, [START_REF] Ryabinkin | Fast Numerical Evaluation of Time-Derivative Nonadiabatic Couplings for Mixed Quantum-Classical Methods[END_REF] which relies on time derivatives in a basis of molecular orbitals. In this case, the calculation of the overlaps is skipped, and the NACs are directly computed in terms of molecular orbital (or KS orbital) overlaps (9) Instead of the integration of the TDSE equation, having SJK computed with either method described above, the time-dependent coefficients can also be obtained via the local diabatization method [START_REF] Plasser | Surface hopping dynamics using a locally diabatic formalism: Charge transfer in the ethylene dimer cation and excited state dynamics in the 2-pyridone dimer[END_REF] proposed by Granucci and Persico. 4d In this case, the coefficients are obtained by a unitary transformation between two time steps, with the help of an adiabatic-to-diabatic transformation matrix obtained under the condition sJK = 0.

K JK J d dt y s y = I c / 1/2 2 , , 1 . KM t K K M M K K M M c c e K M c c c c t -D ¹ = " ¹ é ù = - ê ú ë û å KM t 1 1 1 , K M KM n E E K a t - - ae ö = + ç ÷ è ø I E n K a ( ) * 2 2 max 0, Re , M J MJ M J M P c c c t s ® é ù -D = ê ú ê ú ë û ( ) ( ) ( ) ( ) ( ) 1

The basics of the TD-DFTB method

The basic idea of the density functional based tight binding (DFTB) method is to perform a Volterra expansion of the Kohn-Sham (KS) DFT total energy functional around a reference electronic density

, 5 (10)
where is the sum of neutral atomic densities. Within the DFTB+ code, the second-order (SCC-DFTB or DFTB2) [START_REF] Elstner | Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] and the third-order (DFTB3) [START_REF] Gaus | DFTB3: Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method (SCC-DFTB)[END_REF] expansions of the density functional have been implemented.

The total energy within DFTB is derived employing several approximations. [START_REF] Elstner | Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] A minimal basis set of localized atomic orbitals, obtained by solving the DFT equations for atoms, is used in the expansion of KS orbitals. The matrix elements of the part of the Hamiltonian that depends only on the reference density and orbitals are evaluated based on DFT computations. Its diagonal elements are simply equal to the computed atomic orbital energies, and the nondiagonal elements are precomputed for pairs of elements for a range of

( ) ( ) ( ) . JK J K S t t t t y y = -D JK S ( ) ( ) ( ) . ij i j s t t t t c c = -D K c ( ) 0 r [ ] [ ] [ ] ( ) ( ) 2 3 0 1 2 3 0 0 0 0 1 2 3 , , , , DFTB DFTB DFTB DFTB E E E E E r r r dr r dr r dr é ù é ù = + + + + ë û ë û 0 r
internuclear distances and tabulated in so-called Slater-Koster (SK) files. [START_REF] Elstner | Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] Overlaps of atomic orbitals are likewise precomputed and stored. [START_REF] Elstner | Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] The terms of the energy expression depending on the electronic density fluctuation 𝛿𝜌 are simplified using the Mulliken approximation by means of atomic point charges. [START_REF] Elstner | Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] Remaining integrals are parameterized using Hubbard-like parameters for each atom computed with DFT. Furthermore, the nuclear repulsion term and the terms which depend only on the reference density are approximated as a sum of short-range repulsive potentials which depend only on diatomic distance. [START_REF] Elstner | Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF][START_REF] Gaus | DFTB3: Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method (SCC-DFTB)[END_REF] The LR TD-DFTB method applies a similar procedure as LR TD-DFT in the computation of excitation energies, [START_REF] Niehaus | Tightbinding approach to time-dependent density-functional response theory[END_REF] which are given as the eigenvalues of the equation (11) where the matrix elements are written in terms of orbital transitions from occupied orbitals (i, j) into virtual orbitals (a, b) as

In these equations, are two-electron integrals that involve both Coulomb and exchangecorrelation kernels and e represents the orbital energies from a previous DFTB ground state calculation. In the framework of TD-DFTB, the two-electron integrals are calculated applying the generalized Mulliken approximation, according to which the transition densities between different Kohn-Sham orbitals are represented by point charges. [START_REF] Niehaus | Tightbinding approach to time-dependent density-functional response theory[END_REF] The LR TD-DFTB method has recently been expanded beyond the generalized Mulliken approximation by introducing onsite exchange-like integrals. [START_REF] Domínguez | Extensions of the Time-Dependent Density Functional Based Tight-Binding Approach[END_REF] This correction improves the description of np* and sp* states, which were not well described in the original version of LR TD-DFTB. [START_REF] Domínguez | Extensions of the Time-Dependent Density Functional Based Tight-Binding Approach[END_REF] Moreover, long-range corrected (LC) exchange-correlation functionals have also been implemented and benchmarked in the DFTB method (LC-DFTB) 12a and in the LR TD-DFTB method for excited state simulations (TD-LC-DFTB). 12a Benchmark calculations on a set of small chromophores showed that the TD-LC-DFTB method can treat charge transfer excitations accurately. For this class of excitations, DFTB (like any DFT method) using the generalized-gradient approximation underestimates excitation energies strongly. 12b Although the lack of analytical excited-state gradients within TD-LC-DFTB precludes the use

w * * 1 0 , 0 1 w é ùé ù é ùé ù = ê úê ú ê úê ú - ë ûë û ë ûë û A B X X B A Y Y ( ) ( ) ( ) , , , . ia jb ij ab a i ia jb A i aj b B ia bj d d e e = -+ = ( ) ia jb
of this approach in surface hopping simulations at the moment; we will still use it here as benchmark results for single point calculations.

Time-dependent wave function overlaps with TD-DFTB

As discussed above, the core quantity to compute the hopping probabilities are NACs, which can be calculated in terms of the wave function overlaps SJK (Eq. ( 8)) or molecular (or KS) orbital overlaps sij (Eq. ( 9)), depending on the approach chosen. The algorithm for computation of these quantities has been already implemented in the Newton-X code for the calculation of the TD-DFT NACs. [START_REF] Barbatti | Non-Adiabatic Dynamics of Pyrrole: Dependence of Deactivation Mechanisms on the Excitation Energy[END_REF] Here we only describe them briefly to point out the approximations used in the case of TD-DFTB.

The Casida approximation-stating that the excited state wave function can be represented using a CISlike expansion on the basis of singly-excited Slater determinants-is applied. [START_REF] Casida | Recent advances in density functional methods, Part I[END_REF] This approximation was previously utilized in the implementation of the surface hopping dynamics with the TD-DFT method. [START_REF] Plasser | Surface Hopping Dynamics with Correlated Single-Reference Methods: 9H-Adenine as a Case Study[END_REF] The number of the singly-excited Slater determinants that is obtained by applying a single excitation operator to the ground state Slater determinant is 𝑁 $%& = 𝑁 ()) 𝑁 *+,& , and the form of the excited state wave function is (13) where the CIS coefficients are the eigenvectors of Eq. ( 11) and are the singly-excited Slater determinants.

The overlap of the excited state wave functions represented in the basis of Slater determinants at time steps t and is (14) The overlap of two Slater determinants is computed as the determinant of the overlap matrix of KS orbitals obtained at two subsequent time steps. [START_REF] Löwdin | Quantum Theory of Many-Particle Systems. I. Physical Interpretations by Means of Density Matrices, Natural Spin-Orbitals, and Convergence Problems in the Method of Configurational Interaction[END_REF] Furthermore, the overlaps of KS orbitals (Eq. ( 9)) are represented in terms of overlaps of atomic orbitals at two subsequent geometries from the dynamics:

(15) ( ) , , J J i a ia i a X Y y = + F å + X Y ia F t t + D ( ) ( ) ( ) ( ) ( ) ( ) * ,
.

I J I J ia jb ia jb ia jb t t t X Y X Y t t t y y + D = + + F F + D å i c µ f ( ) ( ) ( ) ( ) ( ) ( ) * . i j i j t t t c t c t t t t t µ n µ n µn c c f f ¢ ¢ + D = + D + D å 10
These atomic orbitals do not form an orthonormal set because they are computed at subsequent geometries in the dynamics, which slightly differ from each other. As already mentioned, the atomic orbital overlaps are precomputed in DFTB at different internuclear distances and stored in Slater-Koster (SK) files. Thus, during an SCC-DFTB computation, these integrals are not calculated, but they are read from the SK files.

During the dynamics, it might happen that some of the atom pairs from two subsequent geometries are displaced with respect to each other by very short distances, which are not covered within the precomputed atomic orbital overlaps because they are not relevant for simulations of single molecules. In this case, we evaluate NACs by applying a linear extrapolation using the computed values from the last two time steps. This problem does not happen with the 3 rd order DFTB and the 3ob Slater-Koster set optimized for it, because atomic orbital overlaps in the standard 3ob set are computed for interatomic distances larger than 0.02 Bohr. This problem might occur in simulations with the 2 nd order DFTB, where the mio Slater-Koster set should be used, in which the atomic orbital overlaps are defined for minimum interatomic distance of 0.4 Bohr. Therefore, we suggest using the 3 rd order DFTB and the 3ob basis set to circumvent the extrapolation of the nonadiabatic couplings.

Since DFTB is a valence-electron-only method, core occupied orbitals do not appear in the construction of Slater determinants. This fact speeds up the computation of wave function overlaps for two reasons: (1) the dimension of the Slater determinants is reduced; (2) the number of Slater determinants is decreased from to . The formal scaling of the determinant derivative method is hence reduced from to and the scaling of the orbital overlap method is reduced from to .

We tested the nonadiabatic couplings between S3 and S2 states in one sequence of the TD-DFTB surface hopping dynamics of thiophene by comparing the computed NACs between these states with the values computed from the TD-PBE wave functions imposing the geometries from the TD-DFTB dynamics. In the case of the TD-DFT wave functions, all occupied orbitals are used in the construction of ground state and excited Slater determinants. In both cases, the NACs were computed with the local diabatization method. [START_REF] Plasser | Surface hopping dynamics using a locally diabatic formalism: Charge transfer in the ethylene dimer cation and excited state dynamics in the 2-pyridone dimer[END_REF] The values of NACs obtained from TD-DFTB and TD-PBE wave functions are in excellent agreement (Figure 1), which indicates that the neglect of core molecular orbitals does not significantly affect the excited-state overlaps. Moreover, TD-DFTB excited-state wave functions seem to approximate their TD-DFT counterparts accurately also in regions of strong nonadiabatic coupling.

Excited-state transition density matrix analysis

To analyze the nature of the excited states occurring during TD-DFTB surface hopping dynamics, we interfaced the TheoDORE program with the DFTB+ program. TheoDORE is a program suite intended for the electronic wave function analysis based on state and transition density matrices between electronic states. [START_REF] Plasser | Analysis of Excitonic and Charge Transfer Interactions from Quantum Chemical Calculations[END_REF] The implemented methods are particularly suitable for obtaining insight into excitonic structure in large conjugated systems [START_REF] Mewes | Excitons in poly(para phenylene vinylene): a quantum-chemical perspective based on high-level ab initio calculations[END_REF] and are, therefore, a natural choice for analyzing CPPs. The most important part of this analysis refers to the determination of charge-transfer character of the excited states and their spatial localization properties by computing charge transfer and fragment participation ratio numbers.

Additionally, the TheoDORE program also implements the computation of the natural transition orbitals of the excited states and their participation ratio of natural transition orbitals.

Segmenting a molecule into suitable fragments, one can evaluate the extent of charge transfer between fragments upon excitation by computing the value (16) where a and b are molecular orbitals belonging to fragments A and B, is the element of the transition density matrix in the localized basis, for the transition from the ground state 0 into the excited state J. 42a

The charge transfer properties of the excited states are described with the charge transfer (CT) index, which is defined as a sum of the nondiagonal elements divided by their total sum (W), ( )

2 0 1 , 2 J J AB ab a A b B D Î Î W = å 0 J ab D J AB W 12 (17)
They have values between 0, in the case of localized excitations (Frenkel excitons), and 1, in the case of charge transfer excitations.

The participation ratio (PR) defines the number of fragments over which an exciton is delocalized. It is defined as 42a (18) where the first term quantifies the extent of the hole delocalization and the second term the extent of the electron delocalization over fragments.

A compact description of the excited states in terms of a relatively small number of transitions can be achieved in the basis of natural transition orbitals (NTOs). [START_REF] Martin | Natural transition orbitals[END_REF] The NTOs are obtained by a singular value decomposition of the transition density matrix represented in a molecular orbital basis, .

where the matrices U and V collect the sets of hole and electron orbitals, respectively, and λi are the weights of the corresponding transitions between them. 42a

The participation ratio of natural transition orbitals (PRNTO) defined as (20) quantifies the number of transitions between NTOs necessary to describe the excited state. This quantity was initially introduced as a collectivity measure in Ref. [START_REF] Luzanov | Application of transition density matrix for analysis of excited states[END_REF] and its physical meaning has recently been elucidated in the context of quantum entanglement. [START_REF] Plasser | Entanglement entropy of electronic excitations[END_REF] Its upper bound is the number of nonzero singular values ( ) of , and it is equal to this number in the case when all nonzero singular values have the same magnitude. Oppositely, the lowest value of PRNTO, reached in the case when the excited state can be represented by only one transition between NTOs, is 1. 42a 1 .

J J A B A B A CT ¹ = W W å 2 2 1 , 2 J J J A B A B A B B A PR é ù ae ö ae ö = W + W ê ú ç ÷ ç ÷ è ø è ø ê ú ë û å å å å 0,J D 0, 1 2 ( , ,...) , J T diag l l = D U V 2 2 ( ) , i i i i PRNTO l l = å å SVD N 0,J D

Computational Details

We created sets of 500 initial conditions (ICs) (geometries and momenta) using a harmonic-oscillator Wigner distribution, [START_REF] Wigner | On the Quantum Correction For Thermodynamic Equilibrium[END_REF] based on normal modes in the ground states of [8]CPP and [10]CPP. The vertical excitation energies and oscillator strengths for the transitions to the six lowest-lying states were computed with the TD-DFTB method with the third-order expansion of density [START_REF] Domínguez | Extensions of the Time-Dependent Density Functional Based Tight-Binding Approach[END_REF] and 3ob-1-1 SK set [START_REF] Gaus | Parametrization and Benchmark of DFTB3 for Organic Molecules[END_REF] (Cartesian coordinates of the ground states of [8]CPP and [10]CPP are given in the Supporting Information, SI-I).

The dispersion interactions were included via a Lennard-Jones potential [START_REF] Elstner | Hydrogen bonding and stacking interactions of nucleic acid base pairs: A density-functional-theory based treatment[END_REF] with parameters adopted from the Universal force field (UFF). [START_REF] Rappe | UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations[END_REF] The sampled geometries are used in the simulations of the photoabsorption spectra with the nuclear ensemble approach [START_REF] Crespo-Otero | Spectrum simulation and decomposition with nuclear ensemble: formal derivation and application to benzene, furan and 2-phenylfuran[END_REF] employing TD-DFTB excitation energies and oscillator strengths. The simulated spectra of both molecules are featured by broad bands peaked in both cases at ≈ 2.95 eV (Supporting Information, SI-II).

Nonadiabatic excited-state dynamics simulations were performed on TD-DFTB potential energy surfaces.

The sets of ICs for dynamics propagation were selected from the initial set of 500 ICs starting from the bright S6 and S5 states applying a very narrow energy window (2.95 ± 0.05 eV) around the maxima of the absorption spectra. This procedure produced 72 ICs in the S5 state and 35 ICs in the S6 state for [8]CPP.

We selected 35 ICs from the S5 and 17 ICs from the S6 state for dynamics propagation. In the case of [10]CPP, 64 ICs in the S5 state and 32 initial conditions in the S6 state are filtered, from which we selected 32 ICs in the S5 and 16 ICs in the S6 state for dynamics propagation. We limited the number of trajectories to ≈ 50 per molecule having in mind that the 3-ps propagation of larger number of trajectories would require substantial computational time.

The TD-DFTB dynamics (GGA functional) of [8]CPP lasted ~10 CPU-minutes per time step (on Intel(R) Xeon 3.0 GHz processor), from which the electronic structure and gradients computations elapsed ~90% of the time and ~10% of the time was spent in nonadiabatic couplings computation. For the comparison, the same simulations with the TD-PBE/6-31G method require ~110 CPU-minutes for the electronic structure and gradient computations and ~2 minutes for nonadiabatic coupling computations on the same processor. Nonadiabatic couplings were computed in both simulations with the overlap derivative method.

Thus, the construction of excited state functions in the valence-only approximation leads to ~50% time saving for computations of nonadiabatic couplings of [8]CPP with the OD method. In the case of determinant derivative methods the time saving would be even more pronounced due to the 5 th order scaling with the number of occupied orbitals used in the Slater determinant constructions.

Nonadiabatic events between excited states were taken into account by FSSH corrected for decoherence with the SDM method. [START_REF] Granucci | Critical Appraisal of the Fewest Switches Algorithm for Surface Hopping[END_REF] Since TD-DFTB as a single-reference method suffers from the same problems as TD-DFT in the description of S1-S0 NAC topology, [START_REF] Levine | Conical intersections and double excitations in time-dependent density functional theory[END_REF] we introduced a condition according to which the internal conversion between S1 and S0 states is assumed, and trajectories are stopped, whenever the energy gap between these states drops below 0.15 eV. Newton's equations of motion were integrated using the velocity Verlet algorithm [START_REF] Swope | A Computer-Simulation Method for the Calculation of Equilibrium-Constants for the Formation of Physical Clusters of Molecules -Application to Small Water Clusters[END_REF] with the time step 0.5 fs. The TDSE was integrated using the fifth order Butcher's algorithm [START_REF] Butcher | A Modified Multistep Method for the Numerical Integration of Ordinary Differential Equations[END_REF] with the time step 0.025 fs, applying interpolated electronic properties between classical steps. The nonadiabatic couplings between excited states were evaluated numerically using the orbital derivative method. [START_REF] Ryabinkin | Fast Numerical Evaluation of Time-Derivative Nonadiabatic Couplings for Mixed Quantum-Classical Methods[END_REF][START_REF] Stojanović | New Insights into the State Trapping of UV-Excited Thymine[END_REF] The dynamics was propagated for 3 ps.

The excited states occurring during dynamics were analyzed based on their transition density matrices. For this purpose, the molecules are fragmented into phenylene rings. CT, PR, and PRNTO values for the active excited state during the trajectories were computed using the TheoDORE program.

For the comparison with the TD-DFTB results, we also computed vertical excitation energies of both molecules with the long-range corrected TD-LC-DFTB method implemented in the DFTB+ code and with TD-DFT based on three different functionals (PBE, 53 B3LYP, [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF] and CAM-B3LYP [START_REF] Yanai | A New Hybrid Exchange-Correlation Functional Using the Coulomb-Attenuating Method (CAM-B3LYP)[END_REF] ), using the GAMESS-US code. [START_REF] Schmidt | General atomic and molecular electronic structure system[END_REF]

Results and discussion

Ground and excited-state geometries

The optimized geometries of the ground and the first excited states of [8]CPP and [10]CPP are shown in Their ground states are featured by alternant dihedral angles between subsequent phenylene rings with the values ±32.1˚ in [8]CPP and ±35.5˚ in [10]CPP. These values are in excellent agreement with the ones obtained in TD-CAM-B3LYP optimization (32˚ and 36˚). 26b

Vertical excitation energies

The vertical excitation energies and oscillator strengths of the first six excited states of [8]CPP and [10]CPP are compared with the corresponding TD-DFT values obtained with the PBE, B3LYP, and CAM-B3LYP functionals in Table 1.

Among the first six excited states at the ground-state geometries, both molecules have two pairs of doubly-degenerate states, S2/S3 and S4/S5. One of the pairs, which corresponds to transitions from HOMO to LUMO+1and LUMO+2, is characterized by large oscillator strengths. The other one, corresponding to transitions from HOMO-1 and HOMO-2 to LUMO, is characterized by small or vanishing oscillator strengths. TD-DFTB, TD-PBE, and TD-B3LYP predict the same state order in [8]CPP, according to which bright degenerate states are placed above dark ones. TD-CAM-B3LYP and LC-TD-DFTB invert this order, predicting the bright states to have lower energies. In the case of [10]CPP, TD-DFTB and TD-PBE predict intense transitions to degenerate S4 and S5 states, whereas the S2 and S3 states have low oscillator strengths. This order is inverted with the range-separated functionals and also with the hybrid functional, which predicts that the energies of the close-lying pairs of degenerate states differ only by 0.03 eV.

Table 1. Computed and experimental vertical excitation energies and oscillator strengths of excited states of [8]CPP and [10]CPP.

Molecule TD-DFTB TD-PBE TD-B3LYP TD-CAM-B3LYP TD-LC-DFTB Exp.

E (eV) f E (eV) f E (eV) f E (eV) f
lying double-degenerate states. The TD-DFTB energies are consistent with TD-PBE, and the TD-LC-DFTB energies are in line with the corresponding TD-CAM-B3LYP values.

According to the TD-DFTB results, the energy of S4/S5 excitation is independent of the system size, and the S0/S1 transition is blueshifted when the ring size is increased. These two features agree with the experimental results. 25b, 28

The initial conditions created using a Wigner distribution have broken symmetries compared to the ground state equilibrium geometry. For this reason, transitions to the S6 state, which are forbidden at the equilibrium geometry, become possible in the ensemble of geometries. The oscillator strength of the S6 state becomes comparable to the oscillator strengths of the degenerate S4 and S5 states. These three states are close-lying, and they can exchange order at different geometries. Thus, we selected initial conditions starting from the S5 and S6 adiabatic states including all six states in the dynamics simulations. The algorithm for the initial condition selection in the Newton-X code computes the probabilities for selecting an initial condition based on the oscillator strength for the selected state, which assures that at all selected geometries the S5 and S6 states are bright states.

Excited-state population dynamics

After the photoexcitation, the molecules undergo ultrafast internal conversion from the initially populated S6 and S5 states via the S4 state to the lower excited states. This process is shown in Figure 3, which reports the fraction of trajectories in each electronic state in each time step. For both molecules, after this initial period taking about 100 fs, the population remains mostly shared between S1 (~40%) and S2+S3 (~40%). The occupations of S4 and S5+S6 remain at about 10% each until the end of the simulations. The ground state is not reached in any of the trajectories during the 3-ps simulations. In general, the trajectories tend to relax to the S1. However, due to the lack of dissipative modes to the environment, once the molecules relax to their S1 states, a part of the population is transferred back to higher excited states, when they reach regions of strong nonadiabatic couplings between S1 and S2.

Furthermore, the S3 state can be easily repopulated from the close-lying S2 state. This recurrent back transfer explains the oscillatory behavior of the average population of these states during the dynamics (Figure 3).

Surface hopping dynamics for [n]CPPs have been reported by Adamska et al. 26b for several ring sizes n between 9 and 16. Their simulations were based on configuration interaction singles (CIS) on top of a semi-empirical Hamiltonian, propagated for 500 fs. Like our results, they also found that internal conversion to the S1 state is reached on the ultrafast scale. However, in their case, once the S1 is populated, the trajectory remains there, and the S1 state reaches 100% occupation within 200 fs. The main reason for this difference is that Adamska and co-workers propagated the trajectories using constant-temperature Langevin dynamics with solvation, 2b while we employed a microcanonical ensemble in the gas phase. The nuclear equations of motion in their case include a velocity-dependent friction term, which emulates the kinetic energy transfer to the solvent.

There is still a third set of dynamics simulations for [6]CPP, [8]CPP, and [10]CPP reported in Ref. 26c .

These simulations, based on wave-packet propagation for 200 fs, were performed involving four excited states computed on a 26-dimensional model Hamiltonian parameterized at the TD-B3LYP level.

Decay mechanisms

After the S1 state is reached during dynamics, three kinetically competitive decay processes from S1 to the ground state may occur-fluorescence, internal conversion, and intersystem crossing to the T1 state 25d followed by phosphorescence to the S0 state. Available experimental and computed quantum yields, rate constants, and lifetimes are summarized in Table 2. F =F range of applied triplet oxygen concentrations, the reaction of triplet quenching with oxygen is several orders of magnitude faster than competing phosphorescence process, meaning that the efficiency of triplet quenching is 100%. 25d The quantum yield of internal conversion was then estimated as .

From the quantum yields of these three processes, one can conclude that internal conversion (Φ ./ = 0.48, Table 2) and intersystem crossing (Φ .4/ = 0.44) are the main deactivation mechanisms in [8]CPP, whereas fluorescence has negligible quantum yield (Φ 56 = 0.08). In [10]CPP, fluorescence becomes the dominant relaxation mechanism with a quantum yield of 0.46, whereas the quantum yield of intersystem crossing is significantly decreased (Φ .4/ = 0.17) compared to its value in [8]CPP. Internal conversion remains an important deactivation pathway in [10]CPP, with a quantum yield of 0.37. A recent wave packet dynamics simulation [START_REF] Reddy | Indirect Intersystem Crossing (S1 → T3/T2 → T1) Promoted by the Jahn-Teller Effect in Cycloparaphenylenes[END_REF] of relaxation process in cycloparaphenylenes within triplet manifold, which included T3-T1 states, revealed the ultrafast internal conversion to the T1 state. Strong Jahn-Teller interactions occurring between doubly-degenerate T3 and T2 states, which induce this ultrafast relaxation, promote indirect intersystem crossing from S1 to T1 state (S1→T3/T2→T1).

Even though the intersystem crossing and phosphorescence were proven to be important, we do not consider this reaction channel, because our dynamics did not involve triplet states, and we limit the discussion to fluorescence and internal conversion. 2).

1 IC ISC fl F = -F -F ( ) ( ) 1 / 1 IC fl ISC IC k k k - = + F - 1/ . fl fl k t =
This shortening with the molecular size can be rationalized based on the dynamics results. The distributions of the excitation energies and oscillator strengths in the S1 states of both molecules are shown in Figure 4Error! Reference source not found.. The S1 excitation energy (𝐸 * ) distribution is centered at 1.5 eV in [8]CPP and it increases to 1.8 eV in [10]CPP, indicating the local minima on the S1 surfaces.

The S1→S0 oscillator strength (𝑓) in [10]CPP is distributed towards higher values compared to the one in [8]CPP. Their mean values also increase with the molecular size, from f[8]CPP = 0.04 to f[10]CPP = 0.10.

These two trends contribute to decreasing the fluorescence lifetime (𝜏 56 ) in [10]CPP in comparison to the [8]CPP, because 𝜏 56 is proportional to 𝑓 <= 𝐸 * <> . [START_REF] Crespo-Otero | Spectrum simulation and decomposition with nuclear ensemble: formal derivation and application to benzene, furan and 2-phenylfuran[END_REF] Using the mean values for the excitation energies and oscillator strengths, we may estimate that the ratio between fluorescence lifetimes is 3.6. This value is underestimated compared to the experimental value (~15), but correctly indicates the qualitative trend.

The disagreement between these two values could be partially attributed to environmental effects, since the experimental rate constant measurements were done in tetrahydrofurane, whereas the simulations were performed in gas phase. The fluorescence lifetime shortening from [8]CPP to [10]CPP should directly lead to an increase of the fluorescence quantum yield, as it hinders competing processes. The experimental ratio between the fluorescence quantum yields is 5.5 (Table 2).

The increase in the oscillator strength of the S1 state in [10]CPP in comparison to [8]CPP has been explained by decreasing the energy gap between S1 and S2 states, which increases the vibronic coupling between them 26c and by localization of excitations in the S1 states in larger cycloparaphenylenes, which is accompanied by asymmetrical geometrical distortions leading to an increase in the oscillator strength 26b . We have analyzed the vibrational relaxation using normal mode analysis. [START_REF] Plasser | Excited-State Diproton Transfer in [2,2 '-Bipyridyl]-3,3 '-diol: the Mechanism Is Sequential, Not Concerted[END_REF] The geometries at the S2/S1 hopping point for all trajectories of [8]CPP were collected, and each of them was then aligned with the [8]CPP S0 minimum. The alignment of structures is achieved by removing the translational degrees of freedom by coinciding their centers of masses. The structures are then superimposed by rotation within their center of masses systems. The optimal rotation matrix that maximizes the overlap between two geometries is computed applying the Kabsch's algorithm. [START_REF] Kabsch | A solution for the best rotation to relate two sets of vectors[END_REF] After being aligned, the Cartesian difference between the S2/S1 hopping geometry and the reference geometry was computed and projected on the normal mode basis calculated for the ground state. A map of the active normal modes at the S2/S1 hopping point was built by averaging the computed activations of normal modes over all S2/S1 hopping points.

The relative activation of normal modes at the hopping points with respect to the ground-state geometry are shown in Figure 5Error! Reference source not found.. They identify the main vibrational motions driving the excited state relaxation. The most active are low-wavenumber torsional modes (𝜐 < 120 𝑐𝑚 <= ) and radial breathing modes (RBM, 150 < 𝜐 < 220 𝑐𝑚 <= ). Several normal modes in the region between 230 𝑐𝑚 <= and 400 𝑐𝑚 <= feature small activation levels. They involve bending of Cipso-Cipso bonds with respect to the ring planes. The radial flexural modes (RFM, 𝜐 ≈ 500 𝑐𝑚 <= ) also show small activation levels. They correspond to ring deformations induced by Cipso-Cipso vibrations outward the ring and Cortho-Cortho vibrations inward the ring. Modes with the wavenumbers between 800 𝑐𝑚 <= and 900 𝑐𝑚 <= have medium activation levels. They are mostly ring deformation modes with out-of-plane C-H bending. The modes with wavenumbers about 1200 𝑐𝑚 <= , also with medium activation level, are inplane C-H bending modes. In the region between 1430 and 1790 cm -1 , a series of C-C stretching modes show small activation. They are characterized in more details below. Cycloparaphenylenes, as well as carbon nanotubes, have three different types of C-C stretching modes: G1, G2, and Gio vibrations. [START_REF] Pena-Alvarez | From linear to cyclic oligoparaphenylenes: electronic and molecular changes traced in the vibrational Raman spectra and reformulation of the bond length alternation pattern[END_REF] The G1 and G2 modes correspond to the circumferential stretching of the ortho-ortho and ipso-ipso C-C bonds. In the G1 modes all phenylene rings experience identical in-phase motions with the same magnitudes, whereas in the G2 modes, the adjacent phenylene rings oscillate out of phase. The G2 vibrations have higher frequencies compared to the G1. The Gio modes correspond to ipsoortho C-C bond vibrations and have lower frequencies compared to G1 and G2. The Gio vibrations also have two components corresponding to in-phase and out-of-phase (we denote them as Gio ' ) motions in adjacent rings. Gio (1430-1480 cm -1 ) and G2 (~1790 cm -1 ) modes show small activation levels at the S2/S1 crossing. As expected from our analysis of the vertical excitation energies, the distributions of the excitation energies show significant differences between S1 energies computed with the two methods. The TD-DFTB energies are red shifted in comparison to the TD-LC-DFTB energies by ~1 eV. Nevertheless, not all disagreement between the results is caused by TD-DFTB inaccuracy. We have already discussed that range-separated functionals tend to blue shift the excitation energies of CPPs with respect to the experimental absorption peak for approximately 0.5 eV, whereas the TD-DFTB tends to red shift these energies for similar magnitude (see discussion of Table 1). Moreover, TD-DFTB overestimates CT numbers of the S1 state of both molecules in comparison with the CT numbers obtained from TD-LC-DFTB computations. According to the TD-LC-DFTB, the S1 states are obtained by mixing of Frenkel and charge transfer excitons, with the CT distributions peaked at ~0.65. On the other hand, TD-DFTB predicts almost pure CT character of the S1 states with CT distributions peaked at ~0.9 in both molecules. This indicates that TD-DFTB, similarly as TD-DFT at the level of GGA functionals, overestimates charge-transfer character of excited states. It has been confirmed, for instance, on a set of conjugated systems 62 that TD-DFT with the GGA-level functionals predicts spurious charge transfer character for the lowest excited state whereas predominantly local character is obtained at the ADC(2) level and with TD-DFT with hybrid and long-range corrected functionals.

The PR distributions calculated with TD-DFTB indicate that the S1 states of both molecules are delocalized over 8 fragments in [8]CPP and ~9 fragments in [10]CPP (Figure 7). The inclusion of longrange correction to the functional leads to some localization of excitons in comparison to the TD-DFTB prediction, with the S1 state of [8]CPP mostly delocalized over 6 to 8 fragments and that of [10]CPP showing a broad distribution between 4 to 9 fragments and peaking at 6. (It has been already observed that the localization of excitons can be adequately described only with range-corrected hybrid DFT models. [START_REF] Körzdörfer | Long-range corrected hybrid functionals for π-conjugated systems: Dependence of the range-separation parameter on conjugation length[END_REF] ) According to the TD-LC-DFTB, the broad PR distribution for [10]CPP implies that in this molecule both localized and delocalized excitons occur in the first excited state. The maximum of the PR distribution peaked at 6 fragments indicates a more pronounced exciton localization than in the case of the smaller molecule. This result is in agreement with the dynamics results of Adamska et al., 26b who noticed that localization of excitons, absent in [9]CPP occurs in larger CPPs (n ≥ 12), where they are localized on 5 to 7 rings. This exciton localization is attributed to its self-trapping induced by geometrical distortions. In addition, Nishihara et al. 25a observed, based on the temperature dependence of the fluorescence lifetime, that the excitons of larger CPPs (n = 12) are localized at temperatures smaller than the room temperature and delocalized at higher temperatures.

It is worth mentioning here that, even though, based on energy criterium, the B3LYP functional provides the most reliable description on the excited states of cycloparaphenylenes, it does not correctly describe the delocalization properties of the exciton in the S1 state, which is an important feature of dynamics. We have optimized the geometries corresponding to the localized exciton in [10]CPP starting from several different geometries taken from dynamics at which the system was in the S1 state and the exciton was localized. The optimizations were performed with TD-B3LYP and TD-CAMB3LYP using 6-31G basis set and with a small optimization radius. TD-CAMB3LYP yielded non-symmetric structure and localized S1 state (PR~7), whereas TD-B3LYP always converged to the completely symmetrical structure and completely delocalized S1 state (PR~10). This finding indicates that TD-B3LYP predicts that localization of exciton in the S1 state is not the most favorable energetically, and probably it would be shown in dynamics performed based on TD-B3LYP energies and gradients. Most probably in the TD-B3LYP dynamics the system would converge to the symmetrical structures and delocalized excitons in the S1 state.

We wanted to verify how the single-determinant (SD) approximation of the excited-state wave functions 4a affects the description of their localization and charge transfer properties. We represented the S1 states in the form of single determinants selecting the one with the largest CI-coefficients obtained either with TD-DFTB (rendering SD-DFTB) or TD-LC-DFTB (rendering SD-LC-DFTB). The excitation energies in the SD approximation correspond to simple KS orbital energy gaps. CT and PR were computed with the renormalized SD wave functions. The three quantities for [8]CPP and [10] Some insights into the differences between compositions of the TD-DFTB and TD-LC-DFTB wave functions that induce distinct charge-transfer and exciton localization properties can be gained by analyzing their PRNTO distributions (Figure 8Error! Reference source not found.). The PRNTO distributions show that the TD-DFTB wave functions, built via Eq. ( 13), are represented mainly by a single transition, whereas the TD-LC-DFTB wave functions are represented by more than one important transition between electron/hole NTO pairs. Thus, it is clear that the LC correction not only changes the orbitals (via the action of the HF exchange) but also changes the distribution of determinants contributing to the state description. Based on this constatation, we are in a position to understand the relations between DE, CT, and PR with the different levels.

TD-DFTB and SD-DFTB do not show a significant difference in the excitation energies because the TD-DFTB wave functions are dominated by one transition. In addition, the two-electron integrals in Eq. 10 (that are included in TD-DFTB but not SD-DFTB) are usually small compared to the HOMO-LUMO gap for local/semilocal functionals. On the contrary, when LC corrections are considered, the SD approximation is an inadequate representation of the TD-LC-DFTB result, composed of many determinants. For functionals involving a fraction of Hartree-Fock exchange the mentioned two-electron integrals are typically large. These two effects lead to the strong deviation in the excitation energy between the two levels. The situation of PR is more subtle than that of DE and CT. Although TD-DFTB state is dominated by a single determinant, other determinants still contributes to the final result. And these secondary contributions are responsible for donating the delocalized character to the state, reflected on the large PR.

In SD-DFTB, however, only the main determinant is considered, resulting in more localized densities and smaller PR. The inclusion of LC correction in TD-LC-DFTB leads to the same effect-PR reduction-but due to a different reason, as now the reduced PR is obtained as result of a multideterminantal composition.

As far as we can see, the good agreement of PR computed with TD-LC-DFTB and SD-DFTB was just

Figure 1 .

 1 Figure 1. Comparison of nonadiabatic couplings (a.u.) between S3 and S2 states in the sequence of a trajectory of thiophene computed based on the TD-DFTB and TD-PBE wave functions.

Figure 2 .

 2 Figure 2. The molecules are highly symmetric in their ground states, having the symmetries of the D4d ([8]CPP) and D5d ([10]CPP) point groups. Their characteristic geometrical features are influenced by three different effects: backbone strain, steric effects, and competition between π-conjugation and quinoidalization. 25c

Figure 2 .

 2 Figure 2. Structures of [8]CPP and [10]CPP in the ground and the first excited states optimized with (TD-)DFTB.

Figure 3 .

 3 Figure 3. Fractional occupations of the excited states of [8]CPP (left) and [10]CPP (right) during dynamics.

Fujitsuka

  et al. estimated fluorescence quantum yields, reaction rates, and lifetimes from fluorescence intensity decay experiments 25b and determined triplet lifetimes from phosphorescence decay experiments. 25d Moreover, they estimated quantum yields and rate constants of S1/T1 intersystem crossing and S1/S0 internal conversion from experiments of triplet quenching with triplet oxygen. Based on the phosphorescence spectra, they confirmed that this reaction produces singlet oxygen, and from the triplet oxygen phosphorescence decay profile they determined the quantum yield of the singlet oxygen generation. They assumed that the quantum yield for singlet oxygen generation is equivalent to the quantum yield for intersystem crossing ( ). The rationale behind this hypothesis was that, in the 2 ISC O

Figure 4 .

 4 Figure 4. Distributions of the S1 excitation energies (left) and oscillator strengths (right) of [8]CPP and [10]CPP during dynamics in single trajectories.

Figure 5 .

 5 Figure 5. Relative activation of different normal modes at the geometry near the S2/S1 conical intersection of [8]CPP.

Figure 6 .

 6 Figure 6. Excitation energies of S1 -S6 states of [8]CPP at the geometries obtained by distorting the molecules along three C-C stretching modes, ν170 (Gio, 1432.5 cm -1 ), ν174 (Gio, 1448.7 cm -1 ), and ν201 (G2, 1786.0 cm -1 ). The x-axis represents a fraction of a normalized normal mode for which the geometry is distorted with respect to the equilibrium geometry.

Figure 7 .

 7 Figure 7. Distributions of the excitation energies (left), charge transfer numbers (middle), and participation ratio numbers (right) in the S1 state of [8]CPP (top of the figure) and [10]CPP (bottom of the figure), during dynamics of a single trajectory considering the TD-(LC)-DFTB and SD-(LC)-DFTB methods.

  CPP are shown in Figure 7Error! Reference source not found.. In both cases, the energies of the S1 state at the SD-DFTB approximation are very close to the TD-DFTB values, whereas the SD-LC-DFTB excitation energies are overestimated by approximately 3 eV in comparison to the LC-TD-DFTB values. CT and PR are both insensitive to range separation, and SD-DFTB and SD-LC-DFTB provide very similar distributions. In the SD approximation, CT has a sharp distribution peaked at ~0.85. (A result like that obtained with TD-DFTB). The results for PR, however, are somewhat surprising. For both molecules, PR distributions obtained with both SD approximations are in reasonable agreement with that from TD-LC-DFTB, our highest level.

Figure 8 .

 8 Figure 8. Distributions of the PRNTO of the S1 states reached during dynamics in single trajectories of [8]CPP (left) and [10]CPP (right).

  The CT distributions computed based on TD-DFTB and SD-DFTB both predict almost pure charge transfer character of the S1 states, implying that the highest-weight transitions between electron/hole NTO pairs are the charge-transfer transitions. The inclusion of the LC correction in TD-LC-DFTB reduces the CT value through the multideterminantal contribution to the density. Nevertheless, the LC correction in SD-DFTB does not change the CT values in comparison to SD-DFTB because, by construction, only one determinant is used.

Table 2 Experimental quantum yields (F), decay rates (k), and lifetimes (t) of fluorescence (fl), internal conversion (IC), intersystem crossing (ISC), and phosphorescence (ph) of [8]CPP and [10]CPP.

 2

	Molecule	Φfl	ΦIC	ΦISC	kfl (10 7 s -1 )	kIC (10 7 s -1 )	kISC (10 7 s -1 )	τfl (ns)	τph (µs)
	[8]CPP	0.084 a 0.475 b 0.440 b	0.48 a, b	2.7 c	2.5 b	208.3 d	60 b
	[10]CPP	0.46 a	0.368 b 0.171 b	7.00 a, b	5.5 c	2.6 b	14.3 d	58 b

a Ref.

25b

(The absorption spectra of CPPs were recorded in THF.) b Ref.

25d

. FISC is assumed to be the quantum yield of singlet oxygen (see text). (The experiments are performed in THF.)

c

Computed as applying the experimental rate constants for fluorescence and intersystem crossing and experimental quantum yield for internal conversion.

d Internal fluorescence lifetimes computed as

Acknowledgments

This Project was funded by the Deanship of Scientific Research (DSR) King Abdulaziz University, Jeddah, under grant no. 43-130-35-RG. The authors, therefore, acknowledge DSR support for Scientific Research. Also, the authors appreciate the kind cooperation provided by the Deanship of Scientific Research (DSR), King Abdulaziz University. MB and LS work was supported by Excellence Initiative of Aix-Marseille University (A*MIDEX) and the project Equip@Meso (ANR-10-EQPX-29-01), both funded by the French Government "Investissements d'Avenir" program. MB also acknowledges funding from HPC resources from GENCI-CINES (Grant 2017-A0010810012). TN would like to thank the Laboratoire d'Excellence iMUST for financial support.

imply that these modes play a central role for the relaxation of the excited states down to S1, and for the eventual repopulation of the higher excited states.

Charge transfer and exciton (de)localization in the S1 state

The excitation energies, charge transfer (CT) numbers, and participation ratio (PR) computed with TD-DFTB along a single trajectory of [8]CPP and

Conclusions

In the present study, we implemented the decoherence-corrected fewest switches surface hopping based on LR TD-DFTB, enhanced by transition density analysis through interfaces between Newton-X, DFTB+, and TheoDORE programs. This implementation was applied to study the excited states and relaxation dynamics of two cycloparaphenylenes, [8]CPP and [10]CPP.

The vertical excitation energies of [8]CPP and [10]CPP computed with the TD-DFTB and LC-TD-DFTB were compared with the TD-DFT energies obtained with the PBE, B3LYP, and CAM-B3LYP functionals.

TD-DFTB predicts the same ordering of the states as TD-B3LYP for [8]CPP, which gives the closest excitation energies to the experimental values. For [10]CPP, the order of the close-lying pairs of degenerate states obtained with TD-DFTB does not correspond to the one predicted with TD-B3LYP.

From the comparison of energies, we can conclude that TD-DFTB systematically underestimates the excitation energies by approximately 0.5 eV, which is a typical feature of DFT methods at the GGA functional level. However, it describes two important features of the absorption spectra well: the most intense peaks are independent of system size, and the S0-S1 energy gap is blue-shifted with the increase of system size.

Based on the results of the TD-DFTB surface hopping dynamics and normal mode analysis, we can draw the following conclusions:

(1) Both molecules relax to S1 within 100 fs. Up to 3 ps, there is no indication of further relaxation to S0. In the gas phase, repopulation of higher excited states, especially S2 and S3, does not allow for full population transfer to S1;

(2) The experimentally observed trend of increasing the fluorescence rate constant in CPP [10] in comparison to CPP [8] is reproduced by the TD-DFTB dynamics. It originates from a simultaneous increase of the S0-S1 energy gap and oscillator strength in [10]CPP in comparison to its values in [8]CPP;

(3) The vibrational dynamics driving the relaxation in the manifold of excited states is enabled by low-amplitude activation of C-C stretching modes, which induces a series of state crossings.

These same modes may be responsible for the back transfer in the gas phase.

Associated content

Supporting information

The Supporting Information is available free of charge on the ACS Publications website. Cartesian coordinates, simulated spectrum, potential energy curves.

Author information

Corresponding authors * E-mail: ljiljana.stojanovic@univ-amu.fr; * E-mail: mario.barbatti@univ-amu.fr; website: www.barbatti.org

Notes

The authors declare no competing financial interest.


N Stalin-Muller

K Dang Van

Numerical simulation of the sliding wear test in relation to material properties

Keywords: Numerical simulation, Sliding wear test, Material properties

Although many wear criteria have been proposed, none is intrinsic to the materials and to the structures in contact. We propose new computational methods for modelling the wear phenomenon and to improve the meaning of these criteria. To interpret sliding wear experiments, numerical calculations have been performed to estimate the mechanical state in the vicinity of the wear region and to infer the dissipation mechanism activated. Calculations in two and three dimensions are proposed to estimate the stress, and the total and plastic strain cycles near the contact under repeated loadings. A result is the stress and strain amplitudes which are derived for different hertzian loadings based on measured friction coefficients. To evaluate wear and fatigue properties, plastic strain amplitude and dissipation rate are also estimated. The modification induced by wear of the original contact geometry is accounted for following the proposals of Kapoor and Johnson (J. Mech. Phys., Solids, 34 (3) (1992) 223-239). The initial contact geometry (circular) is modified by wear so that the contact surface increases and the contact stress decreases. To account for these effects, computations are made with different geometries and hertzian pressure varying from initial conditions to the stabilized state corresponding to shakedown conditions. The plastic strain amplitude and the rachetting rate evolution are then obtained.

Introduction

A diversity of wear models have been proposed in the past but few can really be ll<:ed by engineers and designers. In a recli'nt paper Meng and Ludema [ l] analyse the main reasons and conclude that "It is unlikeiy that many applicable equa tions will emerge using the approach of the past'' in spite of the great number of research works already done. One of the main difficulties in modelling comes from the complexity of evaluating the thermo-mechanical state (or even the chem ico-thermo-mechanical state) in the vicinity of regions sub mitted to wear. From a fundamental point of view, several mechanisms can be at the origin of wear, such as, for example, cutting or plastic burnishing in a ploughing process between asperities in abrasive wear. Various authors have proposed that wear takes place by a fatigue process and that the K coefficient of Archard's wear law is a measure of the numbef of the load cycle (or repeated asperity encounters) to produce rupture. However, deriving wear laws from fatigue theories is not an easy task for two principal reasons. First, it is a multiaxial fatigue phenomenon for which very few criteria can be applied successfully. Second, one has to evaluate the multiaxial stress and strain tensors in the regions where wear occurs, which was not really done until now.

Recently Kim and Ludema [START_REF] Kim | A correlation between low cycle faligue properties and scuffing properties of 4340 steel[END_REF] have given an interesting correlation between scuffing and low cycle fatigue properties of 4340 steel. However they were not able to evaluate the mechanical vari.:bles such as stress or strain amplitude which enter t.'1e fatigue laws. Consequently they could not continue their quantitative analysis. Recently, Johnson [START_REF] Kl | Contact mechanics and the wear of metals[END_REF] published an important contribution entitled "Contact mechanics and the wear of metals" in which he recalled the main analysis necessary to estimate the mechanical state in the wear region. His reasoning is based on an interpretation of the wear mech anism as a fatigue or ductile fracture process. This is sug gested by experimental observations using scanning electron microscope and ferrog.-aphy. Johnson concludes that "I now take the view that delamination wear cracks are indeed shear cracks, but of a ductile nature, driven by plastic strain rather than elastic stress intensity". In this process rachetting is an important phenomenon because it induces large accumula tion of plastic deformation which can lead to the detachment of particles.

Various approaches have been followed up to now to esti mate the mechanical parameters in the wear region. Three of these deserve some discussion.

• The slip line field method for plane strain deformation of rigid, perfectly !'lastic materials of Challen and Oxley [ 41 provides the rachetting strain resulting from each cycle.

However, this method overestimates the rachetting effect because it disregards the elastic strains.

• Shake-down theorems have been applied to provide the elastic and plastic shake-down limits (shake-down map) for a hard cylinder rolling and/ or sliding on an elastic plastic half-space.

• The approximate methods of Bower and Johnson [START_REF] Bower | The influence of stnin hardening on cumulative plastic defonnalion in rolling and sliding line contact[END_REF] pro vide an estimate of the ratchet rate for specific values of the friction coefficient.

Another procedure to estimate the mechanical siate near the contact is to consider an elasto-plastic model and a finite element method. This was done, for instance, by Barghava et al. [ 6) with an incremental translation of the loading. Such an approach is rather involved and the authors are limited to a small number of cycles, which is insufficient to provide a precise idea of the stabilized ratchetting rate.

More recently, Dang Van and Maitoumam developed a systematic and reliable finite element algorithm for thennoe lasto-plastic solids which is well adapted for repeated rorn.1g and sliding contact problems [START_REF] Van | Sleldy-Slale now in class ical elastoplaslicity: Applicalions to iepeared rolling and sliding contact[END_REF][START_REF] Dangvanandh | Maitoumam, Thennomec:hanicalstatenearrolling contact area[END_REF][START_REF] Van | Ekuro-pwtic Calculations o/IM Mechanical State in Reciprocating Moving COlllaCis: App lication /0 Frtning Fatigw[END_REF]. Their method also avoids lengthy and repeated calculation resulting from an incremental analysis of the loading history. This method is based on a steady state assumption for an observer moving with the contact load. Two calculations procedures were pro posed. The first is the pass by pass stationary method ( PPSM) and the second the direct stationary method (DSM) for a repeated moving load. This second provides directly the sta �!!!zed state. Recently these methods have been extended to three-dimensional cases by Maitournam et al. [ 10].

The purpose of this paper is to present applications of these new computational methods to wear problems. Our objective is to derive stress, plastic strain tensors and their variations, in particular the ratchetting rate, in contact regions where wear may happen in relation with the loading conditions and the local friction coefficient p.. This work constitutes a first step for a quantitative understanding of wear phenomena.

The numerical method

The steady-state assumption is proposed for an observer moving with the loading. This means that, in this reference frame, all physical and mechanical quantities are time inde pendent and only position dependent along the moving direc tion. Furthermore, the hypothesis of small strains, pemtits derivation of the evolution eciuations governing our problem.

Disregarding the thennomechanical couplings, these equa tions are, the equation of motion or of equilibrium, the bound ary and initial conditions and the constitutive equations. In practice, the main difference between our method and a classical integration scheme lies in the determination of the plastic strain rate and the related internal hardening parame ters. For instance, for a pin pressed against a rotating disk The plastic strain tensor at any point (defined by its coor dinates r, 9 and z) can be evaluated by integration along the loading path:

2 .. [�(r,8.z)]�"= -llt. *� �(r, 9. z) = J(l*�)<r. E. z) dE+�(r. 21T. z)

Description of the studied system

Our numerical simulation concerns the classical experi ment set-up used by, for example, Papaphilippou [START_REF] Papaphilippou | Elude tribologique d'une Conte l gnphiie sp6roidal avlDICI� trempeL[END_REF] on the The friction coefficient is found experimentally to increase during the wear test as shown in Fig. 2. This variation is explained by the fonnation of small oxide islands ( l in Fig. 2), a generalization of oxidation [START_REF] Kim | A correlation between low cycle faligue properties and scuffing properties of 4340 steel[END_REF], and the removing of the oxide film {3) prior to the stabilization ( 4). We dis regard the transition phase during the numerical simulation.

The local friction coefficient is thus taken as constant. The mesh used for the simulation concerns the section of the wear path which belongs to the (r, :z) plane. The surface of contact being very small, the mesh situated near the loading is very fi ne. In particular, the mesh size near the contact is approximately one tenth of a. Fig. 4 gives an idea of the discretization. The two zones A and B are fi nely meshed (up to 300 elements).

The cast iron is assumed to obey a linear kinematic hard ening law. Its mechanical characteristics are listed in Table 1. l, so that the pressure distribution should differ from that predicted by Hertz's theory. Nevertheless, this proposed distribution is an approximation for our initial computations.

For repeated loadings the distribution should tend towards the hertzian distribution even if the imposed load exceeds the static elastic limit providing that P 0/ k is below the elastic shakedown limit.

As an example, let us consider the following case corr e sponding to a 10 N load and a friction coefficientµ equal to ;

' - 'f'e& (SLIDINGDIRl!C110N)

••• • ••• ••4 .. . � ' � •• • •• • • • 4 . • • • ' •;
I • • .-�---->•' . •• � (TRANSVERSE) . . . • .
�st-�-t-�-+��1--�+-�-+-�-+��1--�...J-�-l-�......J §4t-�-t-�-+��1--�+-�-+-�-+��1--�+-�-l-�......J �a+-�-+-�--+��1--�+-�-l...�-L-�-+-�-i-�-.i�--J !:; i2 t--�+-�-+-�-+-�- +�-+-�-1 -�-1-�--i.�--1�-.J �1i-�+-�f- ---il---4�-+�-4-�-1-�-l-�-L-___J IO.j.....--l----i.-_..£1--��-l---+-�-l-•�l---l---I 0 10
Fig. 13. Variation dissipation rate per cycle for P=20N and µ;=0.8.

Calculations have been perfonned for !he lhree nonnal loadings discussed previously and for the following values of the friction coefficient: 0.2, 0.5, 0.6, 0.7, 0.8, 0.9 and J.

The plastic sheal' strain ampfaude t:. e&, whirh see ms to be the most significant parameter for fatigue, is presented in To account for the evolution of the contact area during a wear test, it is necessary to perform other computations with varying contact area and pressure �istribution. As suggested by Kapoor and Johnson [START_REF] Kapoor | Effect of changes in contact geometry on shakedown of surfaces in rolling/sliding conw:t[END_REF], the wear phenomena will induce variations of confonnity and consequently the area of contact and the pressure disuibution. Following these authors the stabilized shape corresponds to the elastic shakedown state. They also give the shakedown loads with the aspect To estimate the contact parameters which will be used we proceed as follows. We measure the width of the sliding track A and assume that A corresponds to one of the transverse lengths of the contact. Taking P0/k• 1.S, we compute the 3 where a represents the initial contact radius.

' • , le-0.S • -.... ... ,•�•-• -•-• -•-• -•-•-•-•-•-• -•-•-•-•-• -• • 1.c-
Note that the contact dimension in the sliding direction remains rather constant. It is for this reason that we suppose B to be nearly constant from the beginning to the end of the wear experiment. In the following calculations, only the transverse dimension of the contact area A changes.

For instance, for P = 20 N and µ. = 0.8 die following cases were studied and presented in Fig. 11 and summarized in Table4.

The variation in resulting shear strain amplitude �E., is presented in Fig. 12. The results obtained with the modified Hertz's distribution, which are much lower than the initial one, seem to be more realistic because the evolution of the geometry of the contact was accounted for. Fig. 13 presents the dissipation rate under the loading.

Finally, another important quantity, the ratchetting rate, can also be estimated by our computation method. The ratch etting rate depends explicitly on the constitutive equation [START_REF] Bower | The influence of stnin hardening on cumulative plastic defonnalion in rolling and sliding line contact[END_REF]. Using the same non-linear kinematic hardening law they considered, we derive numerically the ratchetting rate. To avoid lengthy computations, only 20 cases were examined for P 0/ k. For instance, the evolution of the rate betting rate for plastic strain components is given in Fig. 14 for µ.=O.S.

The ratchetting rate decreases rapidly and tends toward a stabilized value of lo-s per cycle, which confirms the order of magnitude given in Ref. [START_REF] Challen | An explanalion of different regimes of friction and wear using asperity defonnalion models[END_REF].

S. Conclusions

The first difficulties to overcome for modelling wear phe nomena are to evaluate the physical quantities which are responsible for the degradation of material. Until now, only approximate computational methods for particular constitu tive equations, e.g. perfect plasticity, are available. We try to develop systematic and reliable computational tools to respond to this deficiency. These mell:ods should be appli cable to general elasto-plastic constjtutive equations includ ing non-linear kinematic hardening and isotropic hardening, improving our material properties description. This is achieved by a steady-state elasto-plasticity algorithm. In this work, we try to apply this method to wear problems in sim ulating a pin-on-disk experiment. The flexibility of the method for evaluating the mechanical quantities near the con tact region is demonstrated. This constitutes a first step in interpreting wear phenomena in relation to other material properties such as fatigue or ductility.

  = 40 mm. b =I mm. Sphere: d=6 mm. 0isangularpositionoflbeal_,.... llidill&oadle-pllb. Fig. I. Description of the experimental set-up modelled numerically in the space (r, B, z).

(

  Fig. I), one has d d -=-w-<lr d8 in which "' is !he angular velocity, so that the usual plastic strain r3te �=lt.� au should be replaced by �=-lt. *� a9 au where lt. is the plastic milltiplier andfis the yield function.

Fourier's series

  are convenient for these three -dimensional problems in view of the contact geometry. This method, pro posed by Zienkiewicz, was adapted for the cyclic elasto plastic problem by Maitoumam et al. [ 10].

Fig. 2 .Fig. 4 .

 24 Fig. 2. Evolution in time of the friction coefficient during wear tests on a cast iron disk [II].

  .. � .. .. 1 ...... ......... ..... ... 'fable i Material propenies of cast iron Young's modulus, E Poisson's ratio . P Shear yield limi� k Kinematic hardening coefficient, During the test, the hard alumina sphere remains elastic but the cast-iron disk undergoes plastic deformation. This deformation increases the conformity between the two sur faces, thereby reducing the contact stress and promoting shakedown. The evolution of the contact characteristics is schematized in Fig. 3(a), 3(b) and 3(c) which represents the evolution in three steps of the contact area pressure dis tribution during the wear test.

4 .

 4 Numerical simulations Three levels of loading are considered: 10 N, 15 N and 20 N. Contact parameters (the radius a of the initial contact zone, P0/k which is the maximum normalized pressure) are first evaluated with Hertz's elasticity theory. The results obtained are presented in Table 2. .. .. lo... . . . • t ... � ANGUl..AR POsmoN TETA (nil Fig. 5. Distribution ofthetangentialdisplacemcnt U1versus II (rad), I'-ION, p.-0.S. P0/ k exceeds the elastic limit which is approximately equal to 3 as shown by Johnson [ 12). Consequently the real radius should be greater and the maximum P0 lower than the values in Table

:

  ANGULAR POSmON TETA (rd) Fig. 6. Distribution of the nonnal displacement U, vc111 us 0 (rad), F= ION, ii-O.S. � ... .... 1... ..... ..... ... 'A(al} Fig. 7. Distributionoflhe plastic shear strain •,lvCllius II (rad), F= JON, ii=O.S {plasticshalceJown Slate). TErA (ni) Fig. 8. Themaximalplasticshearsttain �,_versus 6(radj, F-ION, ii=O. S (plasticshakedownS181e).

•

  Fig. 9. '!be plastic shear strain amplitude .4•,& for lhe three levels of loading.

Fig. 10 . 5 .

 105 Fig. 10. Dissipation per cycle rate, E •• for three levels of loadings.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Contact area relaled to loading. ment u,, the plastic shear strain E: and the equivalent plastic shear strain � ••• ,. This case corresponds to plastic shake down state as can be seen in Figs. 7 and 8.

Fig. 9

 9 Fig. 9 versus the load parameter P 0/ k and the various values of the friction parameters are considered. The dissipation per cycle near the surface given by:

  05 2 12 22 32 42 ,2 62 72 12 92 aumbetofqcle Fig. 14. Evolution percycleoflhcratcheaillgrateP.lk-3.S, p.-0.S. other dimension using the classical Hertz equation: P= (2/3) XP01TABwithP0= l.Sk. Theresultsobtainedaresum marized in Table