Gabriele D'avino
email: gabriele.davino@neel.cnrs.fr

Rainer Hegger

Dominik Brey

Praveen K Budakoti

Stéphane Méry

Irene Burghardt
email: burghardt@chemie.uni-frankfurt.de

Electron-hole separation in perylene diimide-based self-assembled nanostructures: Microelectrostatics analysis and Kinetic Monte Carlo simulations

Electron-hole separation in self-assembled mesomorphic nanostructures composed of donor-acceptor (DA) co-oligomers is investigated by a combined microelectrostatics and Kinetic Monte Carlo study. The relevant DA dyads are based on perylene diimide (PDI) acceptor moieties covalently bound to fluorene-thiophene-benzothiadiazole donor moieties, which form highly ordered, stacked structural motifs upon self-assembly. 1 These are characterized by efficient electron transport along PDI stacks, while hole transport is almost three orders of magnitude slower. Based upon an atomistic structure obtained by electron diffraction [Biniek et al., J. Mater. Chem. C 3, 3342 (2015)], the energetics of charge separation is characterized by a microelectrostatics analysis. This information is subsequently employed to compute electron-hole separation rates and dissociation yields by Kinetic Monte Carlo simulations. The latter have been calibrated against recent quantum dynamical studies for a reduced one-dimensional representation of the DA system [Brey et al., J. Phys. Chem. C 125, 25030 (2021)].

It is shown that charge separation of "cold" e-h pairs is characterized by dissociation rates around 10 9 s -1 , which are associated with two-dimensional transport features, where the predominant electron transport in the PDI stacking direction is assisted by a secondary mechanism that involves neighboring stacks.

Introduction

The splitting of tightly-bound molecular excitations into free charge carriers in organic materials is a phenomenon that keeps attracting and challenging researchers since decades, motivated by the prospect for their application in photovoltaics. A fascinating aspect of charge separation is that this is the result of phenomena spanning very different length and time scales. The dynamics right after sunlight absorption, as well as individual charge hopping events, are genuine quantum phenomena taking place at the molecular scale, yet the distance charges have to reach to escape their mutual Coulomb attraction exceeds several nanometers. Moreover, charge carriers move in an energy landscape dictated by long-range electrostatic interactions that involve the multipole moments of other molecules in the local environment.

As for the theoretical description of charge separation in organic solar-cell materials, quantum dynamical and quantum-classical simulations of coupled electrons and holes have been successfully applied to describe the first steps of charge separation on the ultra-fast timescale. [START_REF] Popp | Quantum Dynamics of Exciton Transport and Dissociation in Multichromophoric Systems[END_REF][START_REF] Yan | Understanding the free energy barrier and multiple timescale dynamics of charge separation in organic photovoltaic cells[END_REF][START_REF] Kelly | Exciton dissociation and charge separation at donor-acceptor interfaces from quantum-classical dynamics simulations[END_REF][START_REF] Kato | Non-Markovian Quantum-Classical Ratchet for Ultrafast Long-Range Electron-Hole Separation in Condensed Phases[END_REF][START_REF] Brey | Quantum Dynamics of Electron-Hole Separation in Stacked Perylene Diimide-Based Self-Assembled Nanostructures[END_REF][START_REF] Tamura | Ultrafast charge separation in organic photovoltaics enhanced by charge delocalization and vibronically hot exciton dissociation[END_REF] On the other hand, the process of splitting into free charges has been the playground for phenomenological models -specifically Onsager-Braun (OB) theory [START_REF] Onsager | Initial Recombination of Ions[END_REF][START_REF] Braun | Electric field assisted dissociation of charge transfer states as a mechanism of photocarrier production[END_REF] and its improved versions [START_REF] Wojcik | Accuracies of the empirical theories of the escape probability based on Eigen model and Braun model compared with the exact extension of Onsager theory[END_REF][START_REF] Arkhipov | Why is exciton dissociation so efficient at the interface between a conjugated polymer and an electron acceptor?[END_REF][START_REF] Nenashev | Theory of exciton dissociation at the interface between a conjugated polymer and an electron acceptor[END_REF][START_REF] Baranovskii | Calculating the Efficiency of Exciton Dissociation at the Interface between a Conjugated Polymer and an Electron Acceptor[END_REF] -or for kinetic Monte Carlo (KMC) simulations. [START_REF] Bortz | A new algorithm for Monte Carlo simulation of Ising spin systems[END_REF][START_REF] Deibel | Origin of the Efficient Polaron-Pair Dissociation in Polymer-Fullerene Blends[END_REF][START_REF] Lin | Theoretical insights into the charge transport in perylene diimides based n-type organic semiconductors[END_REF][START_REF] Andersen | A Practical Guide to Surface Kinetic Monte Carlo Simulations[END_REF][START_REF] Bässler | Charge Transport in Disordered Organic Photoconductors a Monte Carlo Simulation Study[END_REF][START_REF] Volpi | Theoretical Study of the Charge-Transfer State Separation within Marcus Theory: The C60-Anthracene Case Study[END_REF] The latter are well suited for multiscale simulations, as KMC can be fed with rates from non-adiabatic electron transfer theories, such as the popular Marcus theory, its Marcus-Levich-Jortner (MLJ) extension, along with the related Fermi's Golden Rule (FGR) approach, [START_REF] Nitzan | Chemical Dynamics in Condensed Phases[END_REF][START_REF] Landi | Hole Hopping Rates in Organic Semiconductors: A Second-Order Cumulant Approach[END_REF][START_REF] Stehr | Comparison of different rate constant expressions for the prediction of charge and energy transport in oligoacenes[END_REF] whose parameters can be obtained from accurate quantum chemistry calculations. The adoption of realistic morphologies for materials and interfaces obtained with atomistic molecular dynamics, and the description of long-range electrostatic effects on the energy landscape, have been paths explored by researchers in the field in order to improve the realism in the simulation of electronic processes in photovoltaic materials. In this context, microelectrostatic (ME) calculations [START_REF] D'avino | Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics[END_REF] have emerged as an important tool to capture electrostatic and polarization effects.

Most of the research on organic photovoltaics has focused on donor-acceptor heterojunctions, [START_REF] Yu | Polymer Photovoltaic Cells -Enhanced efficiencies via a Network of Internal Donor-Acceptor Heterojunctions[END_REF][START_REF] Günes | Conjugated Polymer-Based Organic Solar Cells[END_REF] which are the systems realizing the highest photovoltaic performances. A considerable effort has been, however, put in single-component solar cells, motivated by the promises for a simpler sample preparation, higher structural control and possibly lower photovoltage losses. [START_REF] Shi | Improved efficiency of single-component active layer photovoltaics by optimizing conjugated diblock copolymers[END_REF][START_REF] Madhu | Free Charge Carriers in Homo-Sorted π-Stacks of Donor-Acceptor Conjugates[END_REF] Except for the recent report of efficient charge generation in molecular homojunctions made of polycrystalline α-sexithiophene, [START_REF] Dong | Orientation dependent molecular electrostatics drives efficient charge generation in homojunction organic solar cells[END_REF] most research efforts have focused on co-oligomer molecular dyads made by covalently linked electron donor and acceptor units. [START_REF] Madhu | Free Charge Carriers in Homo-Sorted π-Stacks of Donor-Acceptor Conjugates[END_REF][START_REF] Lee | Challenges and Opportunities in the Development of Conjugated Block Copolymers for Photovoltaics[END_REF][START_REF] Sommer | Donor-acceptor block copolymers for photovoltaic applications[END_REF] Fullerene-based systems have been rather popular given the excellent properties of C 60 as an electron acceptor material, which have remained unrivalled until recently. Perylene diimide (PDI) derivatives have also emerged as valuable candidates in this context, [START_REF] Qu | Donor-spacer-acceptor monodisperse conjugated co-oligomers for efficient single-molecule photovoltaic cells based on non-fullerene acceptors[END_REF][START_REF] Bu | Monodisperse Co-oligomer Approach toward Nanostructured Films with Alternating Donor-Acceptor Lamellae[END_REF][START_REF] Logsdon | Charge Separation Mechanisms in Ordered Films of Self-Assembled Donor-Acceptor Dyad Ribbons[END_REF][START_REF] Schwartz | Perylenediimide-Based Donor-Acceptor Dyads and Triads: Impact of Molecular Architecture on Self-Assembling Properties[END_REF][START_REF] Hartnett | Influence of Anion Delocalization on Electron Transfer in a Covalent Porphyrin Donor-Perylenediimide Dimer Acceptor System[END_REF] in view of the remarkable electron transport properties [START_REF] Jones | High-Mobility Air-Stable n-Type Semiconductors with Processing Versatility: Dicyanoperylene-3,4:9,10-bis(dicarboximides)[END_REF][START_REF] Stoeckel | Analysis of External and Internal Disorder to Understand Band-Like Transport in n-Type Organic Semiconductors[END_REF] and their stability under ambient conditions. [START_REF] Jones | Tuning Orbital Energetics in Arylene Diimide Semiconductors. Materials Design for Ambient Stability of n-Type Charge Transport[END_REF] Power conversion efficiencies (PCE) up to 2.7% have been reached solar cells employing PDI-based dyad as the only materials in the optically active layer. [START_REF] Qu | Donor-spacer-acceptor monodisperse conjugated co-oligomers for efficient single-molecule photovoltaic cells based on non-fullerene acceptors[END_REF] More recently, a new generation of engineered PDI materials, specifically involving In the present contribution, we consider a co-oligomer dyad composed of a PDI acceptor and a fluorene-thiophene-benzothiadiazole donor unit, which forms mesomorphic selfassembled nanostructures. [START_REF] Schwartz | Perylenediimide-Based Donor-Acceptor Dyads and Triads: Impact of Molecular Architecture on Self-Assembling Properties[END_REF][START_REF] Biniek | Zipper-like molecular packing of donor-acceptor conjugated cooligomers based on perylenediimide[END_REF] Specifically, the oriented lamellar structures were shown to include a characteristic zipper-like molecular arrangement. [START_REF] Biniek | Zipper-like molecular packing of donor-acceptor conjugated cooligomers based on perylenediimide[END_REF] The dyad under investigation belongs to a family of push-pull systems, differing in the number and organization of donor and acceptor units (DA, ADA, DAD), [START_REF] Schwartz | Perylenediimide-Based Donor-Acceptor Dyads and Triads: Impact of Molecular Architecture on Self-Assembling Properties[END_REF] as well as the presence of specific bridging units introduced with the purpose of tuning the photophysics. [START_REF] Liu | Controlling charge separation and recombination by chemical design in donor-acceptor dyads[END_REF] These compounds have been extensively characterized by electron diffraction [START_REF] Biniek | Zipper-like molecular packing of donor-acceptor conjugated cooligomers based on perylenediimide[END_REF] (ED) and scanning probe microscopies, [START_REF] Grévin | High-resolution noncontact AFM and Kelvin probe force microscopy investigations of self-assembled photovoltaic donor-acceptor dyads[END_REF] which have shown the formation of ordered structures characterized by segregated channels for hole and electron transport, which makes these systems interesting for both photovoltaic and ambipolar transport applications.

In this work, we focus on a specific system within this class, depicted in Figure 1 and denoted D 0 δ + A, following the nomenclature of Ref. [40]. This system combines a fluorenethiophene donor part (D 0 ), a benzothiadiazol containing (δ + ) spacer moiety, and a PDI acceptor unit (A). For simplicity, we continue to refer in the following to the abbreviation DA, where D = (D 0 δ + ). The presence of the δ + spacer leads to direct, comparatively slow formation of an intramolecular charge transfer state (τ CT 0.1 ns). [START_REF] Liu | Controlling charge separation and recombination by chemical design in donor-acceptor dyads[END_REF] Regarding transport properties, the measured electron mobility has been found to be by far superior (µ e = 0.02 cm 2 /Vs) to hole mobility (µ h = 3.3×10 -5 cm 2 /Vs), differently from analogous systems with longer donor moieties (i.e., D n δ + A D n , species with n > 0). 42 Considering the D 0 δ + A system as a case study, the purpose of the present work is twofold: First, we present an ME analysis of the abovementioned zipper-like morphology 33,39 of the DA system, based upon an available ED structure. Second, based on the ME analysis, we carry out KMC calculations in order to predict field dependent electron-hole (e-h) dissociation rates and yields. The KMC simulations were calibrated with the help of a recent quantum dynamical analysis for a one-dimensional model adapted to preferential transport in the PDI stacking direction. [START_REF] Brey | Quantum Dynamics of Electron-Hole Separation in Stacked Perylene Diimide-Based Self-Assembled Nanostructures[END_REF] Even though the present analysis is carried out for a system exhibiting a moderate PCE in photovoltaic applications, the present approach can be transposed to a broad class of novel non-fullerene acceptor (NFA) materials which often comply with highly ordered morphologies. [START_REF] Sharma | Perylene diimide based non-fullerene acceptors: top performers and an emerging class featuring N-annulation[END_REF][START_REF] Armin | A History and Perspective of Non-Fullerene Electron Acceptors for Organic Solar Cells[END_REF] In the following, we first provide details of the methods that are employed, focusing on the ME approach (Sec. 2.1) and the KMC simulations (Sec. 2.2). Following this, we report the charge carrier energetics for the DA system under study (Sec. 3.1) which provides the basis for the KMC treatment (Sec. 3.2). Sec. 4 concludes with a discussion of the features of the e-h dissociation process characterized in the present work.

Methods

In this section, we summarize the methods employed in this paper, which combine a characterization of the energetics of the relevant e-h states, using a microelectrostatic analysis (Sec. 2.1), with a KMC propagation approach (Sec. 2.2) that leads to the prediction of e-h separation rates and yields.

Microelectrostatics analysis

Effective on-site energies of the relevant e-h states were computed within the established perturbative framework for charges localized on molecular units. [START_REF] D'avino | Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics[END_REF] Intermolecular electrostatic and induction interactions are described with a classical ME model, parameterized with density functional theory (DFT) calculations. [START_REF] D'avino | Electronic Polarization in Organic Crystals: A Comparative Study of Induced Dipoles and Intramolecular Charge Redistribution Schemes[END_REF] The present ME model assumes implicit hydrogens and sets the polarizability of charged species equal to that of the neutral molecule. Molecular parameters, i.e. atomic charges from electrostatic potential fitting (ESP scheme, [START_REF] Besler | Atomic charges derived from semiempirical methods[END_REF] heavy atom-only fit) and the molecular polarizability tensor, were computed at the ωB97XD/def2-TZVP level. ME calculations were performed for the ED structure published in Ref. 39, using the MESCal code. [START_REF] D'avino | Electronic Polarization in Organic Crystals: A Comparative Study of Induced Dipoles and Intramolecular Charge Redistribution Schemes[END_REF] The energy of a charged or neutral optical excitation localized on a given molecule reads

E x = E g x + ∆ x (1)
where E g x is the gas-phase excitation energy and ∆ x is the environmental contribution. The The environmental energy ∆ x consists of two distinct contributions,

∆ x = ∆ E x + ∆ I x ( 2
)
namely the electrostatic (∆ E x ) and induction (∆ I x ) parts. Here, ∆ E x corresponds to the energy of an excitation localized on a given molecule in the field of the permanent and induced multipoles of the surrounding molecules in the neutral crystal or aggregate, i.e.

∆ E x = σ x i δq x i φ n i ( 3
)
where the sum extends over the atoms i of the probed molecule, δq x i = (q x i -q n i ) is the differential charge, with q x i (q n i ) the atomic charges of the excited (neutral ground-state) molecule, and σ x = -1 for x = e, CT, FE or σ x = 1 for x=h. Further, φ n i is the potential at the atomic sites obtained from a self-consistent ME calculation on the neutral systems, performed with periodic boundary conditions in three dimensions.

The induction term

∆ I x = σ x 2 i δq x i φ x i (4)
measures the polarization response of the medium to a localized excitation. The sum is again over the atoms of the probed molecules, with φ x i being the potentials sourced from the dipoles induced in the environment by the excitation differential charges δq x i . The very large size of the the DA dyads considered in this work poses a computational challenge to the extrapolation of induction energies in the bulk limit. We hence resort to a coarse graining of the environment, placing the probed molecule in a cubic mesh of polarizable points (5 Å lattice spacing), which reproduces the anisotropic dielectric response of the crystal. To this aim, we first compute the dielectric tensor of the crystal with our atomistic ME model, following the approach described in Ref. 47. This yields principal components κ aa =1.74, κ bb =2.08 and κ cc =3.87, directed along the orthogonal crystal axes.

The energy of intermolecular e-h pairs at distance r reads

E eh (r) = E e -E h + E b (r) (5)
where E b (r) = V eh (r) + ∆ I eh (r) is the binding energy accounting for the Coulomb interaction between electron and hole screened by the polarizable environment. This quantity is obtained as the sum of the unscreened electrostatic e-h interaction, V eh = ij (δq h i δq e j )/r ij , and the induction term ∆ I eh . The latter is assessed upon embedding the two molecular ions in the anisotropic lattice of polarizable points described above. ∆ I eh is computed as in Eq. ( 4), with the exception that the sums extend over the atoms of the two molecules, and that φ x i = φ eh i is the electrostatic potential induced by the differential charge density of both ions.

Kinetic Monte Carlo simulations

In the following, we describe several ingredients of the KMC simulations that are employed in the present study. First, we summarize the Bortz-Kalos-Lebowitz (BKL) algorithm that is used (Sec. 2.2.1), then we turn to the computation of the transfer integrals that define the connectivities within the BKL scheme (Sec. 2.2.2), followed by a description of the KMC hopping rates (Sec. 2.2.3) and the computation of e-h dissociation times and yields (Sec.

2.2.4).

KMC algorithm

KMC simulations are carried out using the BKL algorithm [START_REF] Bortz | A new algorithm for Monte Carlo simulation of Ising spin systems[END_REF][START_REF] Deibel | Origin of the Efficient Polaron-Pair Dissociation in Polymer-Fullerene Blends[END_REF][START_REF] Lin | Theoretical insights into the charge transport in perylene diimides based n-type organic semiconductors[END_REF][START_REF] Andersen | A Practical Guide to Surface Kinetic Monte Carlo Simulations[END_REF] which represents a discretespace continuous-time Markov process. That is, a sequence of configurations, {C tn }, are defined which belong to the accessible state space {C} and are visited at times { t 0 < t 1 < . . . < t n < . . . }. In our case, the configurations are e-h configurations on a two-dimensional lattice, where each configuration consists of a single e-h pair with the electron and hole localized at respective lattice sites.

The BKL algorithm belongs to the class of so-called rejection-free KMC algorithms; it is equivalent to the N -Fold Way and Gillespie algorithms. [START_REF] Andersen | A Practical Guide to Surface Kinetic Monte Carlo Simulations[END_REF] At each step of the algorithm, rates k ij are defined for moving from configuration C i to configuration C j , where either the electron or the hole are displaced. As detailed below, these rates depends on the energy difference ∆E ij between the relevant configurations, such that the rates corresponding to individual hopping processes involving the relevant molecular pairs need to be explicitly computed.

Prior storage of the rates as a function of ∆E ij can be carried out, as implemented in the present calculations. The set of rates can be subsumed in a transition matrix {k ij }, which is sparse in our case since only nearest-neighbor lattice sites are connected by physical couplings within our model Hamiltonian, as illustrated in Figure 2.

Starting from a configuration C i , we assume that the configurations C l , l = 1, . . . , N , are accessible such that an overall rate K N i can be defined at a given time, along with partial rates K n i , n = 1, . . . , N -1,

K n i = n l=1 k il , n = 1, 2, . . . N (6)
A uniformly distributed random number is then chosen from an interval going up to the total rate, r ∈ [0, K N i ). Following this, the list of partial sums

K n i is searched until K n-1 i ≤ r < K n i
. This is the condition for selecting the transition from C i to C n . Further, the move from C i to C n is associated with a waiting time ∆t i which is generated from another random number, u ∈ (0, 1],

∆t i = - 1 K N i ln(u) (7)
A simulation run (trajectory) is carried out either until a certain threshold time is reached or a specified condition is fullfilled. In the applications reported below, typical ensembles of 5×10 3 realizations were employed.

Transfer integrals

Transfer integrals were computed by several electronic structure methods including a fragmentbased diabatization scheme [START_REF] Blancafort | Exciton delocalization, charge transfer, and electronic coupling for singlet excitation energy transfer between stacked nucleobases in DNA: An MS-CASPT2 study[END_REF] which is employed preferentially, and several orbital-based and DFT based approaches, [START_REF] Voityuk | Fragment charge difference method for estimating donor-acceptor electronic coupling: Application to DNA π-stacks[END_REF][START_REF] Schober | Critical analysis of fragment-orbital DFT schemes for the calculation of electronic coupling values[END_REF][START_REF] Wu | Constrained Density Functional Theory and Its Application in Long-Range Electron Transfer[END_REF] as detailed in the Supporting Information (Sec. S1). Transfer integrals in PDI assemblies are known to be highly sensitive to the molecular stacking pattern. [START_REF] Lin | Theoretical insights into the charge transport in perylene diimides based n-type organic semiconductors[END_REF][START_REF] Geng | Theoretical discussions on electron transport properties of perylene bisimide derivatives with different molecular packings and intermolecular interactions[END_REF] All transfer integrals are assumed to be constant; the effect of thermal fluctuations will be introduced later in the rate expressions. Figure 2 show the relevant pathways, including a dominant quasi-one-dimensional electron transport along the stacking direction (b axis) and additional electron transfer pathways towards neighboring stacks, orthogonal to the stacking direction (i.e., along the a axis), including a "close contact" (cc) pathway between neighboring PDI stacks. (The latter was determined from a slightly modified structure rather than the original ED structure, since the latter generates an unphysical contact between PDIs belonging to neighboring stacks, see the Supporting Information, Sec. S2, for details.) Meanwhile, hole transport is restricted to the b axis direction with a transfer integral which is almost two orders of magnitude smaller than the one for the dominant electron transfer pathway along the PDI columnar stacks.

Table 1 summarizes the relevant transfer integrals, only three of which -i.e., V e b , V e cc , and V h b -represent pathways that are of appreciable importance in the KMC simulations.

In practice, the negligible electron transfer integrals V e a and V e a,b , which are several orders of magnitude smaller than V e b , were discarded in the KMC set-up.

Table 1: Transfer integrals corresponding to the pathways illustrated in Fig. 2. As detailed in the Supporting Information (Sec. S1), the transfer integrals were computed using a quasidiabatization scheme [START_REF] Blancafort | Exciton delocalization, charge transfer, and electronic coupling for singlet excitation energy transfer between stacked nucleobases in DNA: An MS-CASPT2 study[END_REF] (and by the Fragment Orbital DFT 50 approach in the case of V e a and V e a,b ). The respective electron transfer integrals (labeled with superscript e) describe couplings along the b stacking direction (V e b ), a close-contact type coupling between neighboring π stacks (V e cc ), and much weaker couplings along the a axis (V e a ), and along the combined a and b axes (V e a,b ). The hole transfer integral (V h b ) is also small as compared with the dominant electron transfer integrals.

Transfer integral

[eV] V e b 9.7×10 -2 V e cc 2.8×10 -2 V e a 2.0×10 -5 V e a,b 5.6×10 -5 V h b 2.0×10 -3

KMC hopping rate determination

The hopping rates employed in the BKL algorithm are obtained within a Fermi's Golden Rule (FGR) treatment. [START_REF] Nitzan | Chemical Dynamics in Condensed Phases[END_REF][START_REF] Landi | Hole Hopping Rates in Organic Semiconductors: A Second-Order Cumulant Approach[END_REF][START_REF] Stehr | Comparison of different rate constant expressions for the prediction of charge and energy transport in oligoacenes[END_REF] Here, the general rate expression obtained within a second-order perturbation theory treatment is re-cast in terms of the overlap of spectral functions f D (E) and f A (E) for the charge donor (D) and acceptor (A) species, respectively,

k DA (∆E) = 2π |V | 2 ∞ -∞ dE f D (E + ∆E)f A (E) ( 8
)
where V is the electronic coupling (i.e., the transfer integrals of Table 1 in our case) and ∆E is the energy gap between the initial (D ± A) and final (DA ± ) configurations for a charge transfer event D ± A -→ DA ± . In our system, the energy gap is composed of three contributions,

∆E = ∆E ME + ∆E field + ∆E dis
where ∆E ME , ∆E field and ∆E dis relate to the electrostatic interaction captured by the ME treatment, the external electric field, and the static energetic disorder, respectively. The ∆E ME values are either directly accessible from the computed ME data (see Figure 4) or can be constructed based on these data, e.g., for the case of hole migration where the symmetry equivalence operations on the 2D lattice are used. The static disorder contribution ∆E dis is realized by a Gaussian disorder model (GDM) whose standard deviation is chosen as σ d = 50 meV, as a suitable estimate for an ordered mesomorphic system. Further, f D and f A are given as

f D (E) = m D n D g D,i m D (T ) χ D,i m D |χ D,f n D δ(E D,i m D -E D,f n D -E) f A (E) = m A n A g A,f n A (T ) χ A,i n A |χ A,f m A δ(E A,i n A -E A,f m A + E) (9)
where i and f refer to the initial and final states (e.g., charged or neutral) of the donor vs. acceptor species and the vibrational overlap integrals of type χ

D,i m D |χ D,f n D determine the contribution of the vibrational transition m D (i) -→ n D (f )
to the rate, together with the thermal weights g D,i m D (T ) and g A,f n A (T ).

In our implementation, the functions f D and f A are computed from a finite number of N eff = 10 effective modes, 1,53-55 which are constructed from the full set of normal modes of the donor and acceptor fragments (i.e., N = 132 modes for the PDI fragment and N = 234 modes for the fluorene-thiophene-benzothiadiazole fragment) by an orthogonal transformation. The convergence properties of this procedure are illustrated in the Supporting Information (Sec. S3.1). This approach can be considered as intermediate between the numerically expensive inclusion of all normal modes, and approximate procedures relying on a single-mode approximation. Duschinsky effects [START_REF] Nitzan | Chemical Dynamics in Condensed Phases[END_REF] are not included at this level of treatment, i.e., the normal mode set is assumed to be unchanged upon charge transfer. As further detailed in the Supporting Information (Sec. S3.1), the normal-mode frequencies and geometry displacements between the charged and neutral species were obtained from gradient information based upon DFT calculations.

Finally, the finite sums over δ-functions in Eq. ( 9) are smoothed by convolution with a Gaussian whose width is chosen as σ = 40 meV at T = 300 K; this value is adjusted such that the short-time rate matches results from our recent quantum-dynamical calculations for a 1D variant of the present system. [START_REF] Brey | Quantum Dynamics of Electron-Hole Separation in Stacked Perylene Diimide-Based Self-Assembled Nanostructures[END_REF] Details of this procedure are reported in the Supporting Information (Sec. S3.2).

The FGR treatment as described above can be seen as an alternative to traditional Marcus rates and the Marcus-Levich-Jortner (MLJ) variant which accounts for an effective high-frequency quantum mode in the charge transfer event. [START_REF] Nitzan | Chemical Dynamics in Condensed Phases[END_REF][START_REF] Landi | Hole Hopping Rates in Organic Semiconductors: A Second-Order Cumulant Approach[END_REF][START_REF] Stehr | Comparison of different rate constant expressions for the prediction of charge and energy transport in oligoacenes[END_REF] As reported in the literature, these rate approaches do not necessarily yield identical results. In the case of molecular crystals, the definition of the classical reorganisation energy appearing in the Marcus and MLJ expressions is not unambiguous, [START_REF] Mcmahon | Evaluation of the External Reorganization Energy of Polyacenes[END_REF] in contrast, e.g., to solute-solvent systems. For this reason, the FGR approach tends to be better suited for the present system.

E-h dissociation rates and yields

The e-h dissociation rate ("escape rate") is computed as the inverse of a mean first passage time (MFPT). [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF] Here, each KMC time series -or trajectory -is taken to represent a realization of a Markov process which leads to the escape from a spatially confined domain.

That is, for each realization within a set of N KMC trajectories, the time τ n , n = 1, 2, . . . , N , represents the time required to leave the bound domain for the first time. The estimate for

T escape = lim N →∞ 1 N N n=1 τ n (10)
The escape rate is defined as the inverse of the MFPT,

k escape = 1 T escape (11)
In our set-up, the bound domain is taken to correspond to e-h configurations with finite binding energies, up to the threshold where ME data were collected (see Sec. 3 for a detailed discussion). This definition coincides with the boundary employed in our recent quantumdynamical study [START_REF] Brey | Quantum Dynamics of Electron-Hole Separation in Stacked Perylene Diimide-Based Self-Assembled Nanostructures[END_REF] where a flux-over-population approach [START_REF] Farkas | Keimbildungsgeschwindigkeit in übersättigten Dämpfen[END_REF][START_REF] Kramers | Brownian motion in a field of force and the diffusion model of chemical reactions[END_REF][START_REF] Hänggi | Reaction-rate theory: fifty years after Kramers[END_REF] was used to compute the eh dissociation rate. Indeed, the computation of escape rates based upon the MFPT vs. flux-over-population approaches have been shown to be equivalent for Kramers escape type problems under rather general conditions. [START_REF] Reimann | Universal equivalence of mean first-passage time and Kramers rate[END_REF] Using the same definition of the bound domain, dissociation yields φ KMC are computed in terms of the percentage of trajectories that reach the boundary. In the context of the MFPT computation, the simulation time is chosen such that a given trajectory either escapes from the bound domain, or else recombines. In the latter case, the trajectory does not contribute to the MFPT.

For comparison, we also report e-h dissociation rates computed using the Onsager-Braun model, 7,8,12,62

k OB escape = 3 µ e 4π 0 r r 3 b exp - ∆E k B T J 1 (2 √ -2b) √ -2b (12)
where b = e 3 E/(8π 0 r (k B T ) 2 ), µ is the effective mobility at the DA interface, r is an effective dielectric constant, ∆E is the Coulomb binding energy of the initial electron-hole pair (often taken as ∆E = e 2 /(4π 0 r r b ) within the Coulomb approximation), J 1 is the Bessel function of the first order, E is the electric field and r b is the initial e-h separation.

In the context of OB theory, the dissociation yield is defined as

φ OB = k OB escape k rec + k OB escape (13)
In Sec. 3.2.1, Onsager-Braun results will be contrasted with the escape rates and dissociation yields obtained from a KMC treatment.

Results

Using the various ingredients summarized above, we now report the energetics of e-h dissociation as obtained from the ME treatment (Sec. 3.1), along with KMC simulations (Sec.

3.2) based on the ME data combined with the transfer integrals of Sec. 2.2.2 and the rate computation approaches explained above.

Charge carrier energetics: Environmental effects

Localized charge carriers

The system under investigation is a large DA dyad where the D and A moieties each present a distinct charge density and ensuing electrostatic properties. DFT calculations (ωB97XD/def2-TZVP level) reveal a negligible D to A charge transfer in the ground state, as can be concluded from the analysis of the atomic charges from Mulliken and Löwdin population analysis, which yield a charge on the A unit of -0.03 and -0.01 e, respectively. Hence, the electrostatic properties can be understood as a superposition of the two sub-units. The dyad features a dipole moment of 1.9 Debye that points perpendicular to the D-A direction, approximately parallel to the benzothiadiazole axis. This dipole is essentially sourced from the D unit and the benzodithiazole spacer, given that A is a centrosymmetric, dipole-free, molecule.

The electrostatic potential generated by the dyad shown in Figure 3a, is, however, not governed by the modest value of the dipole, but rather by the pronounced quadrupolar character of its sub-units. The D moiety, including the δ + spacer, is characterized by a negative electrostatic potential above the π-conjugation plane, which is characteristic of hole-transporting materials, e.g. pentacene or poly/oligo-thiophenes. [START_REF] D'avino | Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics[END_REF] On the other hand, the quadrupole of the PDI A moiety is determined by the strongly electron-withdrawing diimide groups, leading to a positive potential above the molecular plane that is typical for electron-transporting materials, e.g. perfluoropentacene. [START_REF] D'avino | Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics[END_REF] Figure 3b provides a visual rendering of the electrostatic potential in the crystal, as obtained from self-consistent ME calculations for the neutral system. We note that the crystal structure employed in the calculations presents inversion symmetry, with the small dipoles of different dyads compensating each other. The electrostatic field in the crystal is hence determined by the superposition of the quadrupole moments of the dyad sub-units in the unit cell. Figure 3b presents a characteristic striped character, with positive and negative potential in the regions occupied by A and D, respectively. This is the result of the molecular electrostatic properties discussed above, and of the zipper-like supramolecular arrangement, with segregated stacks of D and A units, characterizing this family of DA dyads.

The crystal potential in Figure 3b governs the energy landscape of charge and energy carriers in these systems. Specifically, it affects the electrostatic component ∆ E that can be seen as an average of the potential, weighted by the differential charge density of the corresponding charged or neutral excitation, see Eq. 3. It is therefore essential to identify the spatial region involved in the charge rearrangement upon excitation. Figure 3c shows the HOMO and LUMO amplitude isosurfaces, that can be considered to a good approximation as the boundaries of the hole and electrons density. Consistent with expectations, hole and electron densities are localized on the D and A units. Similarly, the superposition of the two isosurfaces is representative of the differential charge density of the intramolecular CT excitation.

The joint analysis of Figure 3b,c hence fully rationalizes the trends obtained with ME calculations for ∆ E , reported in Table 2. The electrostatic contribution ∆ E has indeed very similar magnitude and opposite sign for hole and electron, reflecting the difference in the potential for the regions occupied by D and A units. The net effect of ∆ E is therefore a 0.5 eV reduction of the transport gap.

We emphasize that this is a unique property of these DA dyads that has no counterpart in common single-component organic semiconductors, which have been extensively investigated with similar approaches. [START_REF] D'avino | Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics[END_REF][START_REF] D'avino | Electronic Polarization in Organic Crystals: A Comparative Study of Induced Dipoles and Intramolecular Charge Redistribution Schemes[END_REF][START_REF] Topham | Ionization in organic thin films: Electrostatic potential, electronic polarization, and dopants in pentacene films[END_REF][START_REF] Ryno | Electronic Polarization Effects upon Charge Injection in Oligoacene Molecular Crystals: Description via a Force Field[END_REF][START_REF] Poelking | Impact of mesoscale order on open-circuit voltage in organic solar cells[END_REF][START_REF] Li | Host dependence of the electron affinity of molecular dopants[END_REF] In the latter, electrostatic interactions in the solid state determine an approximately rigid shift of IP and EA, leaving the gap almost unaltered.

The rigid shift of IP and EA in single-component organic semiconductors is determined by the fact that the electron and hole density occupy the same region of the π-conjugated core, hence probing the electrostatic potential in a very similar way. Conversely, in our DA dyads, hole and electrons are spatially confined in different regions as a result of the structural segregation of D and A units, ultimately probing an electrostatic potential of opposite polarity. We expect that similar conditions could be met in other organic materials, featuring covalently-bound D and A units and suitable supramolecular packing, which represent potential candidates for non-fullerene acceptors (NFAs) for organic photovoltaics. [START_REF] Armin | A History and Perspective of Non-Fullerene Electron Acceptors for Organic Solar Cells[END_REF] The proper mastery of electrostatic interactions in DA-based systems may represent a powerful tool for the bandgap engineering of organic semiconductors, as already demonstrated in two-component systems. [START_REF] Schwarze | Band structure engineering in organic semiconductors[END_REF] The electrostatic term ∆ E for the intramolecular CT excitation and for a Frenkel-type exciton localized on the PDI A are also reported in Table 2. ∆ E CT is approximately the difference of the electron and hole analogue, in line with the above considerations on the effect of the crystal field on individual carriers and gap. The electrostatic term for the Frenkel exciton is very small, mirroring the minor charge rearrangement following this localized excitation.

While the electrostatic contribution to environmental energies reflects the intimate re-lationship between structure and molecular electrostatic properties, the induction term ∆ I follows from the dielectric response of the material. Owing to the similar dielectric constants of conjugated organic materials (ε r ∼ 3), this contribution is much less sensitive to the fine structural details and is here evaluated with a coarse-grained model for the anisotropic dielectric response of the system (see Sec. 2). The induction term has the well-known effect of reducing the IP-EA gap, due to the stabilization of the individual charge carriers by the induced dipoles in the polarizable medium. Unlike common single-component systems, we note that the gap reduction is quite asymmetric in the present case, with a contribution 0.2 eV larger in magnitude for the electron with respect to the hole. This result can be ascribed to the more localized character of the electron, determining a stronger reaction field from induced dipoles in the surroundings, as captured at a qualitative level by the Born equation. [START_REF] D'avino | Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics[END_REF] The magnitude of the induction term for the intramolecular CT exciton (henceforth referred to as CT intra ), ∆ I CT = -0.67 eV is much larger than previously reported values for intramolecular excitations, [START_REF] Londi | Comprehensive modelling study of singlet exciton diffusion in donor-acceptor dyads: when small changes in chemical structure matter[END_REF] actually comparable comparable to that of a single carrier. This is moderately surprising, considering the very large dipole moment associated with an electron transfer between D and A units. The induction term for the Frenkel exciton is very small (see Table 2), which is again due to the negligible charge reorganization upon excitation.

Energetics of charge separation

In order to discuss the energetics of charge separation in this system, we next consider intermolecular e-h pairs as function of the distance r. For this purpose, we pin a hole on a given dyad in the lattice and perform ME calculations for each position of the electron on molecules in the ab plane that is relevant to charge transport, for r up to 6 nm. Since all dyads in the unit cell are equivalent by symmetry, this approach completely characterizes the energetics of e-h separation.

The results shown in Figure 4a provide an energetic picture of the first steps of the e-h separation spatially resolved in two dimensions. We note that the e-h pair for nearest-neighbors molecules along the π-stacking direction b (denoted CT inter in the following), characterized by a binding energy of -0.47 eV, has practically the same energy as the intramolecular CT excitation (E CT intra = 0.88 eV), within a 20 meV tolerance that is considerably smaller than the accuracy of ME calculations. This is consistent with the fact that the electron addition and removal energies of the dyad are essentially determined by the A and D units, and with the very similar e-h exciton binding energies for the CT intra and CT inter states, i.e., intramolecular and intermolecular (nearest-neighbor) CT.

The e-h separation from the intramolecular CT state, either directly photoexcited or populated after an ultrafast relaxation in the singlet exciton manifold, is predicted to be characterized by a first barrierless electron transfer CT intra -→ CT inter , and by a series of states at larger e-h distance, whose energies approximately follow a screened Coulomb potential, see Figure 4b.

The effective Coulomb potential shown in Figure 4c has been fitted to the functional form

V (r) = - e 2 4π 0 r (r -r 0 ) + E 0 (14)
with the reference energy E 0 = 90.7 meV, the relative permittivity factor r = 2.548, and the reference distance r 0 = 8.66 Å (noting that r 0 differs from the e-h distance in the CT intra state, which corresponds to 21.4 Å).

In the context of the KMC calculations reported below, realizations of the e-h separation dynamics based upon either the full set of ME data or else the above fit will be shown. We will also consider a subset of the ME data in the 1D π-stacking direction (along the b axis),

shown in Figure 4b (red framed squares).

KMC results

In the following, we will describe KMC results based on the general scheme outlined in Sec. 4a and Figure S13 of the Supporting Information. E-h configurations which involve a displacement of the hole are obtained by symmetry equivalence;

2.2.1, combined with the transfer integrals of

these configurations are relevant since hole migration is included in the KMC calculations, even though the hole mobility is far inferior to the electron mobility in the present system (i.e., the electron mobility was experimentally found as 2.0×10 -2 cm 2 /Vs, and the hole mobility as 3.3×10 -5 cm 2 /Vs 42 ). The simulations describe a bulk system with open-boundary conditions and no explicit modeling of electrodes. A simplified simulation set-up for the computation of escape times is used, where KMC trajectories are terminated once they reach the boundary of the domain defined by the ME data, to avoid recrossings. (However, lattice positions beyond the ME data set are necessary to compute time-dependent mobilities as described below; for this purpose, the on-site energies were set to zero for e-h configurations beyond the maximal e-h distances defined by the ME data.) Simulations were carried out at T = 300 K, with a static disorder value of σ d = 50 meV, which is suitable to describe an ordered system. [START_REF] D'avino | Electrostatic phenomena in organic semiconductors: fundamentals and implications for photovoltaics[END_REF] All calculations were started from an initial e-h configuration corresponding to CT intra , where the electron and hole are both localized on the dyad which is at the center of the ME domain. This initial condition describes the system after the break-up of the photogenerated Frenkel exciton, similarly to the situation considered in our complementary quantum dynamical simulations. [START_REF] Brey | Quantum Dynamics of Electron-Hole Separation in Stacked Perylene Diimide-Based Self-Assembled Nanostructures[END_REF] From time-resolved spectroscopic observations, it can be concluded that the formation of the CT intra state in the D 0 δ + A system is comparatively slow, of the order of 100 ps, [START_REF] Liu | Controlling charge separation and recombination by chemical design in donor-acceptor dyads[END_REF] such that the e-h pair dissociation can be taken to correspond to a "cold" process, i.e., in the absence of vibronically hot states.

Recombination of e-h pairs on a D 0 δ + A dyad was included as a loss process in the simulations, restricted to geminate recombination (that is, for a single e-h pair with a quasiimmobile hole, the electron has to diffuse back to the hole to recombine on the same dyad).

Since the recombination time of the dyad system in the thin film material could not be reliably determined experimentally, 42 two approximate recombination times were referred to, which were experimentally obtained for the present dyad system in solution phase, [START_REF] Liu | Controlling charge separation and recombination by chemical design in donor-acceptor dyads[END_REF] i.e., τ rec = 410 ps for the D 0 δ + A dyad in chlorofom, and τ rec = 2.8 ns for the closely related D 1 δ + A system in trichloroethylene (TCE). For related systems, even shorter recombination times in the thin-film material were reported. [START_REF] Nübling | Block Junction-Functionalized All-Conjugated Donor-Acceptor Block Copolymers[END_REF]

Electron-hole dissociation rates and yields

Field-dependent escape rates and e-h dissociation yields were computed with a BKL simulation protocol as described in Sec. 2.2. As explained in Sec. 2.2.4, e-h dissociation rates ("escape rates") were obtained from the mean first passage time (MFPT), which corresponds to the averaged time needed by the KMC trajectories to reach a pre-defined boundary. This boundary is here taken to delimit the domain defined by the ME data set, corresponding to a distance of twelve monomer units from the initial, localized e-h configuration. This boundary is similar to the definition employed in our quantum dynamical simulations for the reduced 1D system 5 where the escape rate was determined from a flux-over-population approach. In the following, we also report dissociation yields, which were equally computed with reference to the boundary as discussed above; i.e., the dissociation yield φ KMC refers to the percentage of trajectories that have reached the boundary within the simulation time, see Sec. 2.2.4.

Figure 5 shows the resulting field-dependent e-h dissociation rates ("escape rates") that were computed with the BKL set-up as described above, compared to a standard computation based on Onsager-Braun (OB) rate theory (see Sec. Several situations are considered using the BKL set-up: First, the full 2D ME data set as compared with the Coulomb fit according to Eq. ( 14). More precisely, a hybrid description is chosen in the latter case, such that short e-h distances are described by the original ME data (see Fig. 4c). Second, a reduced one-dimensional KMC dynamics is considered, since electron transfer along the b axis, i.e., the PDI stacking direction, can be assumed to give the dominant contribution to the e-h separation rate. Again, a comparison is carried out between the ME data in 1D and the corresponding 1D Coulomb fit according to Eq. ( 14); in the latter case, the same hybrid approach is chosen as above, where the most strongly bound e-h pairs are included at the ME level (see Figure 4c). The aim of this comparative analysis is to examine whether (i) a 1D transport picture is appropriate, (ii) the approximation of the ME data set by the Coulomb form of Eq. ( 14) gives reasonable results.

Even though the Coulomb fit is combined with the original ME data at short e-h distances as mentioned above (see Sec. 3.2.1 and Figure 4c), the effect of the fit at larger e-h distances will be to smooth over local energetic traps. Indeed, the spread of ME energies, at a given e-h distance, can exceed the 50 meV static disorder and reach ∼150 meV (see Figure 4b), such the smoothing effect of the Coulomb fit will likely promote the transport efficiency, especially in 1D.

Figure 5a shows that the e-h dissociation rates obtained from the KMC treatment lie around 10 -9 s within the experimentally relevant field range of E = 5 × 10 6 Vm -1 to E = 2 × 10 7 Vm -1 . These rates are plausible in comparison with preliminary time-resolved experiments conducted on the present system 42 showing decay components of ∼0.2 ns -1 and ∼10 ns -1 (however, these experiments also indicate that competing processes may occur, which are not accounted for in our model). For the KMC simulations based on the 2D ME data and the 2D Coulomb fit, similar rates are obtained, showing that a Coulomb fit is valid, at least beyond a certain minimal e-h distance as discussed above. Further, it is seen that simulations based on the restricted 1D Coulomb fit also result in dissociation rates in the same range as the 2D simulations, and even yield slightly larger rates. However, this contrasts with the observation that simulations based on the restricted 1D ME data yield significantly smaller rates, almost an order of magnitude below those obtained from the 2D ME data. This discrepancy seems striking, given the good agreement between the ME vs.

Coulomb fit realizations in 2D. As already mentioned, we attribute this result to the fact that the ME data taken along the b axis exhibit local energetic traps up to about 0.15 eV, exceeding the static disorder of 50 meV and substantially slowing down the escape process in 1D. (For further reference, we note that an energy difference of 0.15-0.2 eV typically lies in the tail of the thermally activated transitions, as illustrated in Figure S11 of the Supporting Information). By comparison, the 2D process permits to circumvent these local traps, due to the non-negligible orthogonal transfer via the V e cc rate. More generally, this is a manifestation of the observation that the lower is the dimensionality, the more the transport becomes sensitive to disorder. From this, we may conclude that the e-h separation process is not strictly reducible to 1D transport. Besides, we also notice that the OB profile shown in Fig. 5a does not match the KMC rates, i.e., the KMC rates tend to become larger and less sensitive to the field, especially at lower field strengths which are here most relevant for experiment.

Complementary to the computed e-h dissociation rates, Figure 5b shows the dissociation yield as a function of field strength. As mentioned above, two values for the recombination rates are considered (τ rec = 410 ps, τ rec = 2.8 ns), which are seen to strongly influence the e-h dissociation yield. Indeed, typical dissociation rates are of the order of several nanoseconds, such that recombination strongly interferes. In general, the dissociation yield is found to be large (∼80% or higher) for all KMC simulations based on the 2D data, both for the ME based simulations and for those employing the Coulomb fit. For the 1D realizations, though, the yield is significantly smaller, especially so for the simulations based upon the restricted 1D ME data set (resulting in yields between 15% to 50%, depending on τ rec ) while the 1D Coulomb fit gives an improved yield (between 80% and 95%, again depending on τ rec ). As in the case of the dissociation rates, this is due to the fact that the Coulomb fit avoids the effects of local traps on the dynamics and renders the dissociation more efficient, also by taking advantage of a favorable entropic contribution. [START_REF] Gregg | Entropy of Charge Separation in Organic Photovoltaic Cells: The Benefit of Higher Dimensionality[END_REF] Finally, the OB results of Figure 5b show that the OB yields φ OB tend to be smaller and, importantly, follow a steeper profile as a function of the external field than the KMC results. In the Supporting Information (Sec. S3.4), similar comparisons are shown for slightly different rate regimes, showing that the main trends of the analysis remain unchanged.

To summarize the key features, we find both from the dissociation rate and the yield that the 2D simulations cannot be straightforwardly reduced to a 1D scenario, as one might expect from the predominant electron transfer along the b direction. The additional, weaker transport pathways play an important role in rendering the process more efficient. The Coulomb fit in 2D, and also in 1D, leads to a smoother transport behavior, by the elimination of local traps, resulting in larger rates and higher dissociation yields.

Time-dependent mobilities and electron-hole distributions

Even though KMC is an approximate propagation method, it does reveal aspects of the real-time dynamics, permitting to characterize the e-h dissociation process in more detail.

To this end, we now consider the time-dependent mobilities and the underlying electron-hole distributions.

Figure 6 shows time-dependent electron mobilities (panel b) that were computed from the time-derivative of the mean e-h distance (panel a),

µ e (t) = 1 E d r eh dt (15)
where we assumed that the hole is approximately static. As can be inferred from Figure 6a, obtained from an ensemble of 10 4 KMC trajectories, at least two dynamical regimes are observed: Within the first ∼5 ns, the time-evolving e-h distance is near-quadratic, leading to a near-linear increase of the mobilities shown in Figure 6b, while the longer-time evolution (> 5 ns) of the e-h distance is linear, resulting in almost constant mobilities as expected. The ratio of the field-dependent mobilities at t = 20 ns is quite well captured by the Poole-Frenkel expression, 17 µ e = µ e0 exp(aE 1/2 ), where a ∼ 2 × 10 -4 (m/V) 1/2 . For the experimentally most relevant field strength of 5×10 6 V/m, the electron mobility reaches ∼ 0.4 cm 2 /(Vs);

for reference, the computed free carrier mobility is given as 0.55 cm 2 /(Vs) as also indicated in Figure 6b. Even though the latter value is far larger than the experimentally observed value (0.02 cm 2 /Vs) obtained by field effect transistor (OFET) measurements, 42 the discrepancy can be explained by the role of disorder whose description is approximate in the present simulation set-up. Notably, we only include site energy (diagonal) energetic disorder, disregarding the fluctuations acting on transfer integrals (off-diagonal disorder).

The computed mobilities can be compared with similar quantities discussed, e.g., in Ref.

71. Interestingly, we see that µ e (t) increases as a function of time rather than the decrease observed in Ref. 71, where the observed mobilities are small and field-independent. This difference is likely related to the strictly localized, "cold" initial e-h state in our simulations and the low-dimensional character of the dynamics. By contrast, in the P3HT:PCBM junction studied in Ref. 71, the initial phase of the dynamics is determined by the injection of a hot carrier, such that the first portion of the trajectory is ruled by relaxation in the density of states (DOS). Hence, in the latter case, the instantaneous mobility drops with the thermalization of the charge.

Complementary to the time-dependent mobilities, Figure 7 shows time-dependent e-h distributions for a field of 5×10 6 V/m, with a static disorder of σ d = 50 meV. The distributions were obtained by computing sub-ensembles of 112 trajectories; overall, 4000 such sub-ensembles were created, whose time-dependent mean values constitute the ensembles that are shown. (This particular sampling was employed to show distributions that are independent of the original lattice discretization; further, we separated the random sampling with respect to the KMC steps from the static disorder sampling.) As expected, a gradual broadening of the e-h distribution occurs, along with a displacement to larger distances.

The displacement away from the initial tightly bound domain also suggests that the e-h dissociation exhibits a high yield, in line with the results shown in Figure 5b.

Discussion and Conclusions

The present DA material represents a paradigm example of a highly ordered self-assembled donor-acceptor segregated morphology where PDI π-stacking determines a dominant electron transport channel. The ME analysis that we presented for a specific, zipper-like structural motif [START_REF] Schwartz | Perylenediimide-Based Donor-Acceptor Dyads and Triads: Impact of Molecular Architecture on Self-Assembling Properties[END_REF][START_REF] Biniek | Zipper-like molecular packing of donor-acceptor conjugated cooligomers based on perylenediimide[END_REF] is ideally suited to gain a quantitative understanding of the local electrostatic effects on the e-h dissociation process. The electrostatic potential is determined by the superposition of the mostly quadrupolar character of the sub-units in the present dyad, yielding the remarkable striped pattern of the electrostatic map represented in Figure 3.

An assumption of the ME treatment is that the (electronic) polarization response of the dielectric environment to the moving charge is instantaneous, leading to an equilibrium treatment; this should be justified in the present mesomorphic structure characterized by a molecular dynamics that is expected to be similar to a molecular crystal. Complemented by KMC calculations based upon the ME-based on-site energies, a full picture of the e-h dissociation process in the present material is obtained. This treatment is complementary to our recent work [START_REF] Brey | Quantum Dynamics of Electron-Hole Separation in Stacked Perylene Diimide-Based Self-Assembled Nanostructures[END_REF] where a quantum-dynamical study of e-h separation was carried out within a reduced one-dimensional model representing an isolated PDI stack.

The intuitive picture of the e-h separation process, as suggested by the available ED structure along with experimentally determined mobilities, points towards an effectively onedimensional electron transport process along the PDI columnar stacks, while hole mobility is much inferior. Furthermore, a "cold" e-h dissociation process is expected since the initial formation of the charge-separated state within a DA dyad species is comparatively slow, excluding a vibronically hot initial state. While this interpretation essentially remains valid, a number of additional aspects emerge from the present analysis. Notably, it is shown that an effective Coulomb potential for e-h separation can be defined, but a detailed dynamical description requires individual ME data within a tightly bound domain of e-h states within a radius involving at least five PDI units. Furthermore, the role of energetic traps is prominent in the 1D transport along stacked PDI columns, as evidenced by a large difference in the e-h dissociation rates resulting from the ME on-site energies as compared with the effective Coulomb potential, if restricted to a 1D description. This difference disappears in the full 2D treatment, due to the fact that additional pathways communicating between neighboring PDI columns are activated, such that energetic traps can be circumvented. As a matter of fact, the transport is not strictly reducible to a one-dimensional description.

Overall, typical escape rates around 1 ns -1 are observed, in line with the quantum dynamical analysis 5 which was used for calibration of the present KMC simulations. The yields obtained from the KMC simulations are comparatively large, of the order of 80% or more, suggesting a quite efficient charge separation process. However, a direct comparison with experimental data for the e-h separation process is not yet available for the present system.

Preliminary time-resolved experiments 42 confirm that rates of the order of 1-10 ns -1 are plausible, but also indicate that competing processes may occur, which are not accounted for in our model. Besides geminate recombination, [START_REF] Liu | Controlling charge separation and recombination by chemical design in donor-acceptor dyads[END_REF][START_REF] Roland | Sub-100 fs charge transfer in a novel donor-acceptor-donor triad organized in a smectic film[END_REF] which could act more rapidly than in our model, excimer formation [START_REF] Son | Spectroscopic Demonstration of Exciton Dynamics and Excimer Formation in a Sterically Controlled Perylene Bisimide Dimer Aggregate[END_REF][START_REF] Singh | Excimer formation effects and trap-assisted charge recombination loss channels in organic solar cells of perylene diimide dimer acceptors[END_REF][START_REF] Würthner | Perylene Bisimide Dye Assemblies as Archetype Functional Supramolecular Materials[END_REF] have been discussed as a source of loss pathways in PDI based systems. The present analysis can, in principle, be adapted such as to include such additional pathways.

From a more general vantage point, the KMC analysis represents an approximate approach to the real-time dynamics of the DA system, but it is generally well adapted to describe slow (and cold) e-h separation processes like the one under study. Further, the method was calibrated by explicitly matching the escape rates to full quantum rates obtained by quantum-dynamical simulations for a reduced 1D realization of the same system. [START_REF] Brey | Quantum Dynamics of Electron-Hole Separation in Stacked Perylene Diimide-Based Self-Assembled Nanostructures[END_REF] Hence, quantum effects like delocalization of charge carriers are to some extent implicitly included, even though a standard KMC protocol was employed. (For reference, see, e.g., Ref.

14 for generalized KMC protocols where delocalization effects are explicitly included in the KMC approach.) As a result, the rates that were obtained via a mean first-passage time tend to be significantly larger than typical Onsager-Braun rates for the same system, as illustrated in Figure 5. Also, the e-h dissociation yields tend to be larger than the corresponding yields from Onsager-Braun theory. The present simulations could be augmented by using master equation approaches [START_REF] Landi | Hole Hopping Rates in Organic Semiconductors: A Second-Order Cumulant Approach[END_REF][START_REF] Landi | Second-Order Cumulant Approach for the Evaluation of Anisotropic Hole Mobility in Organic Semiconductors[END_REF] that are expected to yield more accurate results than KMC.

To summarize, the present protocol, which combines the ME description of a highly ordered DA lattice with KMC simulations of e-h separation, permits a realistic description of the charge separation process. Various generalizations and improvements are in reach, especially as far as the dynamical treatment is concerned. The lattice-based description is shared by a hierarchy of methods, from KMC to master equations and quantum dynamical approaches, with the ME parametrization serving as a common denominator for an accurate description of the charge separation energetics underpinning the dynamics. These approaches

  subscript labels hole charges (x=h, E h = -IP, where IP denotes the ionization potential) or electron charges (x=e, E e = -EA, where EA denotes the electron affinity), or an intramolecular D-to-A charge-transfer (i.e., CT, leading to D + A -) excitation, or else a Frenkel exciton localized on the D moiety (FE). The gas-phase energies of charge carriers, E g h and E g e , were computed from DFT total-energy differences (∆SCF approach) at the CAM-B3LYP/def2-TZVP level. Time-dependent DFT (TD-DFT) calculations with the same basis and functional were employed to compute E g CT and E g FE . Electronic structure calculations were performed with the Gaussian suite, for the molecular geometry in the crystal, as provided in Refs.39-42.

  2.2.4). The comparison between OB and KMC results should be considered with care, given the different possible choices for the effective parameters entering the OB model. OB results presented in Figure5have been obtained using an effective interfacial mobility µ = 0.01 cm 2 /(Vs), consistent with experimental data,42 the effective dielectric constant was set to r = 2.55, and the Coulomb binding energy of the initial exciton was E b = 0.35 eV (obtained from the Coulomb fit of Eq. (14) for r = 21.4 Å and r 0 = 8.66 Å). For comparison, results obtained for larger effective interfacial mobilities are shown in Sec. S3.5 of the Supporting Information; these larger mobilities are closer to the KMC results discussed in Sec. 3.2.2.

Figure 1 :

 1 Figure 1: (a) Molecular structure of the D 0 δ + A building blocks, (b) view of the zipper-type structure under study in the ac plane, (c) view of the structure in the bc plane, i.e., along the PDI stacking direction.

Figure 2 :

 2 Figure 2: Transfer integrals and associated KMC pathways. The green shade of the arrows indicates the coupling strength, with darker shade indicating larger couplings; see also Table 1 which summarizes the values of the transfer integrals. (a) Electron transfer integrals for the PDI acceptors, where V e b determines the dominant transfer along the b axis and V e cc corresponds to the close contact pathway discussed in the text. The remaining transfer integrals V e a along the a axis and V e a,b involving both the a and b directions, are discarded in the numerical treatment due to their small values. (b) Hole transfer integral V h b along the b axis.

Figure 3 :

 3 Figure 3: Electrostatic potential (a) of the D-A dyad at the van der Waals surface and (b) in the ac plane of the quasi-crystalline system -average over the b direction. The potential in the quasi-crystal presents a peculiar striped pattern, with positive and negative potential at the regions occupied by A and D sub-units, respectively. (c) HOMO (red) and LUMO (blue) amplitude isosurface, showing the localization of the hole on the D and of the electron on the A.

Figure 4 :

 4 Figure 4: Binding energies E b (r) of intermolecular e-h pairs, calculated at ME level for a fixed hole at the center and different positions of the electron. (a) Schematic representation of the ab planes of the D-A dyad crystal, with annotated E b (r) values. Blue squares mark the position of PDI acceptors. Red circles correspond to D units, their position is largely offset along the c axis (not shown). The b axis corresponds to the π-stacking direction. The gray-shaded box delimits the unit cell, measuring 1.77 × 0.74 nm 2 . (b) Distance-dependence of the e-h binding energy, approximately following a screened Coulomb potential. Blue squares mark the full ME data set, while blue squares with red frames indicate the e-h states associated with 1D transport along the b axis. The CT intra and CT inter state, which exhibit very similar binding energies, are marked explicitly. (c) The same ME data are shown as in panel (b), together with an approximate Coulomb fit. The e-h pairs at short distances, marked in green, are included separately in the KMC computations when employing the Coulomb fit for the e-h states above -0.3 eV.

Figure 5 :

 5 Figure 5: KMC escape rates and dissociation yield, as compared with Onsager-Braun (OB) theory. All quantities are shown as a function of the external electric field; the experimentally relevant range between E = 5×10 6 V/m and E = 2×10 7 V/m is indicated by a yellow shade. (a) Escape rates, where the recombination time of the e-h pair was fixed to τ r = 2.8 ns. For the curves labeled ME, the original micro-electrostatics data were used, while for the curves labeld fit, the e-h potential was generated from the Coulomb fit shown in Figure 4(c). (b) Dissociation yield computed as the fraction of trajectories that escaped from the bound domain. Results are shown for two recombination times, τ r = 2.8 ns (as in panel (a)) and τ r = 0.41 ns, see text for details.

Figure 6 :

 6 Figure 6: Time-dependent mobilities. (a) Average e-h distance r eh as a function of time for different strengths of the external electric field; the KMC average was calculated using 10 4 trajectories. (b) Resulting time dependent mobilities. The constant dashed line shows the free-carrier mobility which represents an upper threshold for the electron mobilities in the presence of the hole.

Table 2 :

 2 Energies (eV units) of charged and neutral excitations localized on a single molecule, including intramolecular and environmental components.

	Excitation	E x	E g x	∆ x	∆ E x	∆ I x
	h	-5.34 -6.04 0.70 0.24 0.45
	e	-3.97 -3.07 -0.90 -0.24 -0.66
	CT intra	0.88 2.07 -1.19 -0.52 -0.67
	FE	2.09 2.12 -0.03 -0.01 -0.02

Table 1

 1

and the ME data of Sec. 3.1.2. ME data pertaining to the ab plane were computed for 288 e-h configurations, with a fixed hole position as illustrated in Figure

Acknowledgments

We gratefully acknowledge preliminary contributions to the present analysis by Dr. K. Schwinn, Prof. M. Schmidt, Dr. C. Schlesinger, W. Popp, Dr. R. Binder, and Prof. S. Haacke. Dr. M. Brinkmann is acknowledged for providing the ED structure of the donor-acceptor system under study. Funding for this project was granted by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the project MolNanoMat (BU-1032-2).

Associated Content Supporting Information

Supporting Information is available free of charge, including the computation of transfer integrals, details on the close-contact transfer integral between PDI stacks, the computation of Fermi's Golden Rule rates, and the microelectrostatics data set.

Author Information

Corresponding authors * E-mail: gabriele.davino@neel.cnrs.fr, burghardt@chemie.uni-frankfurt.de

Notes

The authors declare no competing financial interest.


Gilles Millérioux
email: gilles.millerioux@esstin.uhp-nancy.fr

Jamal Daafouz
email: jamal.daafouz@ensem.inpl-nancy.fr

Flatness of Switched Linear Discrete-time Systems

Keywords: Switched discrete-time systems, flatness, invertibility I

This note is devoted to flatness for switched linear discrete-time systems. For this class of hybrid systems, algebraic conditions are derived to check whether a given output is flat. Then, a feedforward flatness-based control strategy for trajectory tracking is proposed.

by some so-called defining matrices (see [START_REF] Levine | Flat output characterization for linear systems using polynomial matrices[END_REF]). On the other hand, such a characterization is no longer valid for switched systems. A simple but special situation occurs when each linear mode of the switched system is active during a sufficient large time, in others words when a minimum "dwell-time" is guaranteed. In such a case, we can check for flat outputs for each linear models and one can then design a hybrid flatness-based control that takes care of each local dynamics. However, that amounts not only to discard the hybrid aspect but also this strategy may fail if the minimum "dwell-time" is not guaranteed. In this note, we relax the dwell-time constraint and thus we really consider the hybrid issue. We derive conditions which enable to check whether a given output of a switched linear discrete-time system is flat and a flatness-based control strategy for trajectory tracking is suggested.

Notation : For any integer l, 1 l refers to the l-dimensional identity matrix and O l×l stands for the l × l zero matrix. If irrelevant, the dimension of the zero matrix will be omitted and we shall merely write O. For a matrix X, X T stands for transposition. X † stands for the Moore-Penrose generalized inverse of X, namely the matrix of the same dimension as X T so that XX † X = X, X † XX † = X † , XX † and X † X are Hermitian.

II. PRELIMINARIES AND DEFINITIONS

We examine switching linear discrete-time systems of the form   

x k+1 = A σ(k) x k + B σ(k) u k y k = C σ(k) x k + D σ(k) u k (1)
where x k ∈ R n , u k ∈ R m and y k ∈ R p are the states, the inputs and the measurements, respectively. All the matrices, namely

A σ(k) ∈ R n×n , B σ(k) ∈ R n×m , C σ(k) ∈ R p×n and D σ(k) ∈ R p×m belong to the respective finite sets (A j ) 1≤j≤J , (B j ) 1≤j≤J , (C j ) 1≤j≤J and (D j ) 1≤j≤J .
At a given time k, the index j corresponds to the mode of the system and results from a

switching function σ : k ∈ N → j = σ(k) ∈ {1, . . . , J}. {σ} k 1 +T k 1
refers to the mode sequence {σ(k 1 ), . . . , σ(k 1 + T )}. For a given switching rule σ, the set of corresponding mode sequences over any interval of time of length T + 1 is denoted by Σ T . This set may contain either all possible mode sequences (also called paths) if there are not any repetitive switching patterns or can be reduced if any. We assume that the mode is known, either accessible or reconstructed (see [START_REF] Balluchi | Design of Observers for Hybrid Systems[END_REF] for this reconstruction issue). A key point lies in that no restriction on the time separation October 10, 2008 DRAFT between switches ("dwell time") must be imposed.

Let U be the space of input sequences over [0, ∞) and Y the corresponding output space.

At time k, for each initial state x k ∈ R n , when the system (1) is driven by the input sequence We introduce the subsequent vectors and matrices.

{u} k+T k = {u k , . . . , u k+T } ∈ U ,
For i < 0 M i σ(k) = 0, for i = 0 M 0 σ(k) = D σ(k)
and for i > 0,

M i σ(k) =            D σ(k) 0 p×m . . . . . . . . . C σ(k+1) B σ(k) D σ(
C σ(k+i) A σ(k+i-1) σ(k+1) B σ(k) C σ(k+i) A σ(k+i-1) σ(k+2) B σ(k+1) . . . C σ(k+i) B σ(k+i-1) D σ(k+i)            (2)
with the transition matrix

A σ(k 1 ) σ(k 0 ) = A σ(k 1 ) A σ(k 1 -1) . . . A σ(k 0 ) if k 1 ≥ k 0 = 1 n if k 1 < k 0 Īm = (1 m 0 m×(m•r) ), Īp = (0 p×(p•r) 1 p ) (3)
O i σ(k) =        C σ(k) C σ(k+1) A σ(k) . . . C σ(k+i) A σ(k+i-1) σ(k)        , C i σ(k) = A σ(k+i-1) σ(k+1) B σ(k) A σ(k+i-1) σ(k+2) B σ(k+1) . . . B σ(k+i-1) (4)
u i k =        u k u k+1 . . . u k+i        , u i k =        u k u k+1 . . . u k+i        , y i k =        y k y k+1 . . . y k+i        , y i k =        y k y k+1 . . . y k+i        (5)
October 10, 2008 DRAFT When (1) is driven by an input sequence {u} ∞ k ∈ U and a mode sequence {σ} ∞ k , one has:

y i k = O i σ(k) x k + M i σ(k) u i k (6)
Definition 1. When the system (1) is square (m = p), it is said to be flat if there exists a set of independent variables y k , referred to as flat outputs, such that all system variables can be expressed as a function of the flat outputs and a finite number of its backward and/or forward iterates. In particular, there exist two functions F and G which obey

   x k = F (y k+k F , . . . , y k+k F ) u k = G(y k+k G , . . . , y k+k G ) (7)
where k F , k F , k G and k G are Z-valued integers.

Both deriving algebraic conditions to check whether a given output of ( 1) is flat and providing a feedforward flatness-based control for trajectory tracking require a central notion namely the input invertibility.

III. INPUT INVERTIBILITY

We must distinguish input left and right invertibility. Some recent papers [START_REF] Sundaram | Designing stable inverters and state observers for switched linear systems with unknown inputs[END_REF][6] have addressed the input left invertibility for switched discrete-time systems and the results are recalled in Subsect. III-A. In Subsect. III-B, we derive results on input right invertibility.

A. Input left invertibility

Definition 2. The system ( 1) is input left invertible if there exists a nonnegative integer r < ∞ such that, for all mode sequences in Σ r , for two any inputs {u} r 0 , {u } r 0 ∈ U , the following implication applies:

{y(x 0 , σ, u)} r 0 = {y(x 0 , σ, u )} r 0 ⇒ u 0 = u 0 ∀x 0 (8)
In others words, by input left invertibility, we mean the ability to recover the input u 0 from a finite number of r + 1 measurements y i (i = 0, . . . , r), the state vector x 0 at time k = 0 and the mode sequences {σ} r 0 being known. The term input has been introduced to emphasize the fact that only the input u 0 is expected to be recovered. The least integer r for which (1) is input left invertible is called the left inherent delay or merely the delay.

October 10, 2008 DRAFT Theorem 1. The system ( 1) is input left invertible if there exists a nonnegative integer r < ∞ such that for all mode sequences in Σ r ,

rank M r σ(k) -rank M r-1 σ(k+1) = m (9) Definition 3.
A system is a left r-delay inverse for (1) if, under identical initial conditions x 0 and identical mode sequences {σ} ∞ 0 , when driven by y r k , its output fulfills ûk+r = u k for all k ≥ 0 Theorem 2. Assume that ( 1) is input left invertible with left inherent delay r. The system

   xk+r+1 = P r σ(k) xk+r + B σ(k) Īm M r † σ(k) y r k ûk+r = -Īm M r † σ(k) O r σ(k) xk+r + Īm M r † σ(k) y r k ( 10
) with P r σ(k) = A σ(k) -B σ(k) Īm M r † σ(k) O r σ(k) (11)
is a left r-delay inverse system for [START_REF] Balluchi | Design of Observers for Hybrid Systems[END_REF].

i) The system (1) is input right invertible

ii) There exists a nonnegative integer r < ∞ such that, for all mode sequences in Σ r , the

equation with unknown Q r σ(k) M r σ(k) Q r σ(k) = ĪT p ( 13
)

has a solution

iii) There exists a nonnegative integer r < ∞ such that for all mode sequences in Σ r ,

rank ĪT p M r σ(k) -rank M r σ(k) = 0 (14)
iv) There exists a nonnegative integer r < ∞ such that for all mode sequences in Σ r ,

rank M r σ(k) -rank M r -1 σ(k) = p (15) Proof: i) ⇔ ii) On the interval of time [k, k + r ],
according to (6), one has

y r k = O r σ(k) x k + M r σ(k) u r k (16)
Besides, according to the Definition 4 of right invertibility, it is assumed that x k = 0 at time k.

Thus, from (12) and ( 16) one has

y r k = M r σ(k) u r k = ĪT p ỹk+r
October 10, 2008 DRAFT Input right invertibility means that this equality must be fulfilled for any ỹk+r . As a result, the existence of u r k is strictly equivalent to the existence a matrix Q r σ(k) fulfilling

M r σ(k) Q r σ(k) = (0 p×(p•r) 1 p ) T The matrix Q r σ(k) reads explicitly Q r σ(k) = M r † σ(k) ĪT p ii) ⇔ iii)
The relation ( 14) is the immediate consequence of the algebraic result stating that for any pairs (W, Z) of matrices, the equation W X = Z with unknown X has a solution if and

only if rank(W Z) = rank W . iii) ⇔ iv) We first notice that M r σ(k) =   M r -1 σ(k) 0 (p•r )×m C σ(k+r ) C r σ(k) D σ(k+r )   Consequently, ĪT p M r σ(k) =   0 (p•r )×p M r -1 σ(k) 0 (p•r )×m 1 p C σ(k+r ) C r σ(k) D σ(k+r )   Hence rank ĪT p M r σ(k) = rank M r -1 σ(k) + rank 1 p D σ(k+r )
Finally, noticing that rank 1 p D σ(k+r ) = p, substituting the last equality into (14) gives (15).

When (1) is input right invertible, we are looking for a second dynamical system which drives the system (1) so that its output tracks a reference trajectory ỹk , ỹk+1 , . . . for k ≥ 0. Let us define the following compound vector: 1) is input right invertible with right inherent delay r . The system

ỹi k =        ỹk ỹk+1 . . . ỹk+i        Theorem 4. Assume that (
   xk+1 = P r σ(k) xk + B σ(k) Īm M r † σ(k) ỹr k u k = -Īm M r † σ(k) O r σ(k) xk + Īm M r † σ(k) ỹr k (17)
October 10, 2008 DRAFT with

P r σ(k) = A σ(k) -B σ(k) Īm M r † σ(k) O r σ(k)
is a right r -delay inverse system for [START_REF] Balluchi | Design of Observers for Hybrid Systems[END_REF].

Proof: Notice that 1) is input right invertible, (13) holds and is equivalent to

y k+r -ỹk+r = Īp y r k -ỹk+r = Īp (O r σ(k) x k + M r σ(k) u r k ) -ỹk+r = Īp O r σ(k) x k -Īp M r σ(k) M r † σ(k) O r σ(k) xk + Īp M r σ(k) M r † σ(k) ỹr k -ỹk+r Since (
M r σ(k) M r † σ(k) ĪT p = ĪT p Therefore, Īp (M r σ(k) M r † σ(k) ) T = Īp
The Moore-Penrose generalized inverse being hermitian, one has

(M r σ(k) M r † σ(k) ) T = M r σ(k) M r † σ(k)
and the following equality applies

Īp M r σ(k) M r † σ(k) = Īp (18)
As a result, taking into account (18), y k+r -ỹk+r finally reads

y k+r -ỹk+r = Īp O r σ(k) (x k -xk ) (19)
On the other hand, it is easy to see that

x k+1 -xk+1 = A σ(k) (x k -xk ) (20)
In view of the Definition 5, we must consider that ( 1) and ( 17) are initialized such that x0 = x 0 .

By virtue of (20), we infer that x k -xk = 0 for all k ≥ 0 and thus y k+r -ỹk+r = 0 for all k ≥ 0.

When the system (1) is both left and right invertible with left and right inherent delays r = r , it will be merely said that the system (1) is input invertible with inherent delay r = r = R.

October 10, 2008 DRAFT Remark 4. It is interesting to notice that when the system ( 1) is square, if Σ r = Σ r and the sets contain the all possible mode sequences over any interval of time of length r + 1 = r + 1, input left invertibility is equivalent to input right invertibility. Indeed, in such a case, ( 9) is equivalent to (15).

IV. FLATNESS

A. Flat output characterization

Let us define the inverse transition matrix as

P σ(k 1 ) σ(k 0 ) = P r σ(k 1 ) P r σ(k 1 -1) . . . P r σ(k 0 ) if k 1 ≥ k 0 = 1 n if k 1 < k 0 with P r σ(k) = A σ(k) -B σ(k) Īm M r † σ(k) O r σ(k)
Theorem 5. A componentwise independent output y k of the system (1) assumed to be square (m = p) and left input invertible with inherent delay r, is a flat output if there exists a positive integer K < ∞ such that, for all mode sequences in Σ r+K-1 , the following equality applies for all k ≥ 0: Proof: Assume that the system (1) is input invertible with inherent delay r, thus the left r-delay inverse system (10) exists. Iterating (10) yields:

P σ(k+K-1) σ(k) = 0 (
xk+r+l = P σ(k+l-1) σ(k) xk+r + l-1 i=0 P σ(k+l-1) σ(k+i+1) B σ(k+i) Īm M r † σ(k+i) y r k+i ûk+r+l = Īm M r † σ(k+l) (y r k+l -O r σ(k+l) xk+r+l ) (22)
If ( 21) is fulfilled with l = K, (22) turns into:

xk+r+K = K-1 i=0 P σ(k+K-1) σ(k+i+1) B σ(k+i) Īm M r † σ(k+i) y r k+i ûk+r+K = Īm M r † σ(k+K) (y r k+K -O r σ(k+K) xk+r+K ) (23)
revealing that xk+r+K and so ûu+r+K is independent of xk+r . In particular, (23) holds for xk 0 +r = x k 0 for all k 0 ≥ 0 which implies that xk+r+K = x k+K and ûk+r+K = u k+K for all October 10, 2008 DRAFT k ≥ 0. Therefore, after performing the change of variable k → k -K, we obtain

x k = K-1 i=0 P σ(k-1) σ(k+i+1-K) B σ(k+i-K) Īm M r † σ(k+i-K) y r k+i-K (24)
and

u k = Īm M r † σ(k) (y r k -O r σ(k) ( K-1 i=0 P σ(k-1) σ(k+i+1-K) B σ(k+i-K) Īm M r † σ(k+i-K) y r k+i-K )) (25)
Equations ( 24) and (25) reveal that x k and u k can be expressed as a function of y k and a finite number of its backward and/or forward iterates in the form [START_REF] Sira-Ramirez | Differentially Flat Systems[END_REF], that is to say y k is a flat output.

B. Flatness-based control for trajectory tracking

If we are interested in delivering the flat control (25) in a recursive way, the following theorem can be useful.

Theorem 6. Assume that ( 1) is (left and right) input invertible with inherent delay R. If there exists a positive integer K < ∞ such that for all mode sequences in Σ R+K-1 , 26) is fulfilled, the control u k reads:

P σ(k+K-1) σ(k) = 0 (26) with P R σ(k) = A σ(k) -B σ(k) Īm M R † σ(k) O R σ(k) then the control u k delivered by    xk+1 = P R σ(k) xk + B σ(k) Īm M R † σ(k) ỹR k u k = Īm u R k = -Īm M R † σ(k) O R σ(k) xk + Īm M R † σ(k) ỹR k (
xk + K-1 i=0 P σ(k+K-1) σ(k+i+1) B σ(k+i) Īm M R † σ(k+i) y R k+i u k+K = Īm M R † σ(k+K) (y R k+K -O R σ(k+K) xk+K ) Finally, if (
u k = Īm M R † σ(k) (y R k -O R σ(k) ( K-1 i=0 P σ(k-1) σ(k+i+1-K) B σ(k+i-K) Īm M R † σ(k+i-K) y R k+i-K ))
revealing that u k is a flat control since it only depends on the output.

The problem of assessing the stability of (20) turns into a usual switched stability analysis problem. We can resort to polyquadratic stability ( [START_REF] Daafouz | A poly-quadratic stability based approach for switched systems[END_REF]).

C. Illustrative example

We shall consider a MIMO switched linear discrete-time system of the form (1) and having two distinct modes (J = 2). The respective finite sets (A j ) 1≤j≤2 , (B j ) 1≤j≤2 , (C j ) 1≤j≤2 and

(D j ) 1≤j≤2 are A 1 =        0.3 1
       , B 1 = B 2 =        -1 0.5 0 -1.5 1 0 0 1        , C 1 = C 2 =   1 1 1 1 1 3 1 4   , D 1 = D 2 = 0
The system is driven by some inputs {u} ∞ 0 ∈ U , the set of input sequences ranging between two arbitrary bounds, say 0 and 1 for the example. We shall examine three distinct switching rules: a switching rule which obeys σ(k) = 1 for all k which amounts to consider a mere linear system (case 1), a switching rule which obeys σ(k) = 2 for all k which also amounts to consider October 10, 2008 DRAFT a mere linear system (case 2), a switching rule for which σ(k) alternates randomly between 1 and 2, the occurrence of two consecutive modes 1 being however not allowed.

We first consider the case 1. It turns out that (9) and (15) hold for r = r = 2. Hence, according to the Theorem 1 and Theorem 3, the system is left and right input invertible with inherent delay r = r = 2. The set Σ r is merely composed of a single mode sequence {1, 1, 1}. According to the Theorem 5, y k is a flat output because (21) is fulfilled with K = 2. For the case 2, the situation is exactly the same than the case 1 except that the set Σ r is composed of {2, 2, 2}. As a result, it can be stressed that the minimum dwell time equals 2. Then, we consider the case 3. Such a case is interesting in that the switching rule does not always guarantee the minimum dwell time. It turns out that (9) and (15) hold for r = r = 2. Hence, according to Theorem 1 and Theorem 3, the system is input invertible with inherent delay r = r = R = 2. The set Σ R does not contain all possible mode sequences but only the patterns {1, 2, 1}, {1, 2, 2}, {2, 1, 2}, {2, 2, 1}, {2, 2, 2} , the only ones which are compatible with the switching rule. Equation ( 26) is fulfilled with K = 2.

Consequently the assumptions of Theorem 6 hold and a feedforward flatness-based control can be designed. that (20) is u.a.s, more precisely, it can be checked that it is polyquadratically stable. Figure 1 represents the results of the trajectory tracking based on the flat control delivered by (27). The time plots illustrate that the desired trajectories are well tracked. The convergences toward zero are consistent with the theoretical results.

V. CONCLUSION For switched linear discrete-time systems, algebraic conditions in terms of the state description matrices to check whether a given output is flat have been derived. Then, it has been shown how a feedforward flat control for trajectory tracking can be carried out in the case when a minimum dwell time is not guaranteed. Let us finally mention that one can resort to these conditions and to such a flat control when considering periodic systems or monovariate systems with time-varying relative degree.

  for a mode sequence {σ} k+T k , {x(x k , σ, u)} k+T k refers to the solution in the interval of time [k, k + T ] of (1) emanating from x k and {y(x k , σ, u)} k+T k ∈ Y refers to the corresponding output sequence in the same interval of time [k, k + T ].

Remark 1 .Remark 2 .Definition 4 .Remark 3 .Theorem 3 .

 12433 In the Definition 2 and Definition 3, the initial condition is considered at the particular discrete time k = 0 but can be replaced by any other initial condition x k taken at the discrete time k. It is also shown that the state vector of the left r-delay inverse (10) fulfills xk+r = x k for all k ≥ 0B. Input right invertibilityLet us define a reference trajectory as a sequence of desired outputs {ỹ k , ỹk+1 , . . .} for k ≥ 0. The system (1) is input right invertible if there exists a nonnegative integer r < ∞ such that, for all mode sequences in Σ r , for any reference output ỹr at time k = r and for the initial condition x 0 = 0, there exists an input segment {u} r 0 such that{y(x 0 = 0, σ, u)} r 0 = {0} r -1 0 * ỹr(12)where {0} r -1 0 * ỹr denotes a sequence of null vectors of dimension p over the interval of time [0, . . . , r -1] and ỹr corresponds to the reference output at time k = r . The term input October 10, 2008 DRAFT has been introduced as well to emphasize the fact that we are interested in finding out the input segment {u} r 0 only. The least integer r for which (1) is input right invertible will be called the right inherent delay or merely the delay. Definition 5. A system, with input ỹr k is a right inverse for[START_REF] Balluchi | Design of Observers for Hybrid Systems[END_REF] if, under an identical initial condition x 0 and identical mode sequences {σ} ∞ 0 , it can drives[START_REF] Balluchi | Design of Observers for Hybrid Systems[END_REF] such that y k+r = ỹk+r for all k ≥ 0 Similarly to the Remark 1, in the Definition 4 and Definition 5, the initial condition is considered at the particular discrete time k = 0 but can be replaced by any other initial condition x k taken at the discrete time k. The following statements are equivalent.

  21)Σ r+K-1 stands for the set of mode sequences over the interval of time [k, . . . , k + r + K -1].

  Fig. 1. A : solid line ỹ(1) k and dashed line y

( 2 ) 2 Two distinct sinusoidal trajectories {ỹ ( 1 )

 221 k . C : |y kỹk | 2 . D : |x kxk | 0 have been planified for the two respective entries y

( 2 )

 2 k of the output y k . It can be verified through a LMI approach ([2]) October 10, 2008 DRAFT

  in the proof of Theorem 4, after substituting r by R, since x0 may differ from x 0 , one has October 10, 2008 DRAFT lim k→∞ y k+R -ỹk+R = 0 whenever (20) is u.a.s. On the other hand, iterating (27) K -1 times

	27) is a flat control and guarantees a stable trajectory tracking whenever the system (20) is uniformly asymptotically stable (u.a.s). Proof: On one hand, if (1) is input invertible with inherent delay R, by virtue of (19) and (20) yields σ(k+K-1) xk+K = P σ(k)

Acknowledgement This work is partially funded by a grant from the Agence Nationale pour la Recherche in France (Ref. ANR-05-JCJC-0112-01)