Michele Turelli

Domenico Alberga

Gianluca Lattanzi
email: gianluca.lattanzi@unitn.it

Ilaria Ciofini

Carlo Adamo
email: carlo.adamo@chimieparistech.psl.eu

Theoretical insights on acceptor-donor dyads for organic photovoltaics

Introduction

The field of organic photovoltaics has been steadily growing for the past two decades. A wide and concerted research effort has been successful in improving the efficiency of organic solar cells so as to reach commercially viable values, even if their efficiency still lags behind the conventional inorganic technologies [START_REF] Adamo | Best Research-Cell Efficiency Chart[END_REF] . This research area remains nonetheless extremely appealing thanks to the versatile chemistry of organic compounds that promises virtually unlimited possibility for fine-tuning and purpose-driven design.

However, these benefits come hand in hand with intrinsic difficulties. The physics of organic solar cells is more complex and less understood compared to that of their inorganic counterpart [2][3][4] . The persistence of gaps in our understanding is partly due to the vast chemical space that these compounds occupy which results in a variety of possible behaviours that need to be traced back to more general principles. The coming maturity of the field is signalled by the increasing necessity to rationalise the wealth of experimental data already collected into a theoretical framework capable of directing the investigation and synthesis of novel materials through a list of design principles [5][6][7][8][9] . Indeed, theoretical approaches are becoming increasingly widespread as tools to guide the search for new molecules [10][11][12][13][14][15] , the majority of which, however, are still being discovered and tested experimentally 16 .

The need for a deep theoretical understanding aimed at establishing a series of sound design principles is even more pressing for the class of materials employed in the so-called singlematerial organic solar cells (SMOSCs). The domain of SMOSCs has recently seen a boost due to the ability of these materials to circumvent the disadvantages of the more established donor-acceptor bulk heterojunction architecture [START_REF] Peumans | Materials for Sustainable Energy: A Collection of Peer-Reviewed Research and Review Articles from Nature Publishing Group[END_REF][START_REF] Yu | [END_REF][19] , such as thermodynamic instability 20,21 and the recent identification of molecules with higher conversion efficiencies [22][23][24][25][26][27][28][29][30][31][32] . The treatment of these systems brings an additional layer of complexity both in terms of the physical process involved and their molecular size. Indeed, their size is often too large to allow the application of the quantum chemical methods needed to accurately predict the relevant properties 33 and forces treatments based on less accurate semi-empirical methods 34,35 .

In this regard we present here a computational study of one of these materials employed in a SMOSC active layers that has been shown experimentally to have a remarkable performance 22 . For this system, a full characterisation is difficult and leaves the interpretation of the mechanisms underlying its efficiency unsatisfactory. This makes it an ideal test-case to prove the value and efficacy of theoretical methods aimed at complementing and directing the experimental work.

The molecule belongs to a particular class of materials employed in SMOSCs called dyads 36,37 . A dyad is a system combining an electron donor with an acceptor group by covalently linking them in a single molecular architecture, thus realising a molecular heterojunction 38 . In the majority of cases, dyads employed in organic photovoltaics are obtained via the covalent bonding of a fullerene (C 60 ), playing the role of acceptor, with a πconjugated monodisperse oligomer, acting as donor 39 . The dyad we focus on consists in a dithiafulvalene-functionalized diketopyrrolopyrrole (DPP) acting as donor and a [60]fulleropyrrolidine (Ful) as acceptor (see Figure 1). Experimentally 22 , the two moieties were selected on the basis of reported performances of the separated molecules in the literature rather than through the application of more rigorous purpose-oriented design rules.

In dyads the first step of light-to-energy conversion common to organic active materials, i.e. the dissociation of the exciton into separate charge carriers, is achieved on the molecule itself 16,32 . However, differently from inorganic semiconductors, in dyads the mechanism is more complex and nuanced, involving a series of intermediate steps and transitions that influence the overall performance 4 . Following charge separation, charges move in the active layer towards the electrodes: the efficiency of this transport is fundamental in determining the performance of the solar cell 40,41 . In this work, we will focus on transport dynamics with the aim of pointing out the electronic and structural features of the dyads that are most relevant to enhance solar cell efficiencies. It is worth to notice that, at the best of our knowledge, this is the first time that this kind of theoretical study has been carried out on dyads for SMOSC applications.

Computational details

The modelling of DPP-Ful dyad poses specific challenges due to its nature of ambipolar conductor. Carriers of both signs can be transported independently over morphologies whose microscopic details are not entirely known 22 . From a theoretical point of view this requires the prediction of both ordered and disordered morphologies and to simulate transport in both cases.

The first mandatory steps consist in the gas-phase calculations on the isolated molecule.

These calculations were performed at density functional theory (DFT) level applying the B3LYP functional in conjunction with the triple- Pople's basis set (6-311G(d,p) ) 48,49 .

Dispersion interactions were accounted for through the Grimme's empirical potential (D3) 50 .

All optimisations, including those carried out to compute reorganisation energy λ for both anionic and cationic states, were performed at the same level of theory via the Gaussian09 package 51 .

Ordered morphologies were obtained following the approach adopted in reference 52. We sampled the thermodynamically most relevant crystal structures, i.e. crystal morphologies with a non-trivial probability of being occupied at experimental conditions, via the polymorph predictor module of Materials Studio 53 . The prediction of crystal structures is carried out by arranging the packing motifs within a given space group and then by simultaneously optimizing the molecular packing inside the cell and the unit cell parameters. These steps were carried out using the Dreiding force field 54 , largely used for molecular crystal prediction 10,52 . We limit our search to the most probable space groups reported for similar organic compounds i.e. P2 1 2 1 2 1 , P2 1 /c, P , P2 1 , C2/c and Pna2 1 55 . A number of independent searches have been performed until no new low-energy structures, significantly different from those already obtained, are found. The resulting structures are then clustered and ranked according to their energy and density and those with the highest rank are selected. This ranking is then refined according to the results of full DFT geometry optimizations of the unit cells and the atomic positions performed including periodic boundary conditions (PBC) using the CRYSTAL17 code 56 . In this case, the same exchange correlation functional (B3LYP-D3)

is used together with an optimised basis set, tailored to organic crystals, available in the SI.

Amorphous configurations are extracted following the same computational procedure proposed in reference 52. We obtain the disordered configurations for a system of 500 molecules by taking 150 uncorrelated snapshots from a molecular dynamics (MD) trajectory.

The initial conditions are chosen as to accurately reproduce the disordered phase. DPP-Ful molecules are initially placed at random positions in a cubic cell 270 Å wide. This configuration is evolved through a simulation in the NPT ensemble at high temperature T=800 K and P=1 atm until the volume of the cell reaches convergence. Another NPT simulation is then performed at lower temperature, T=300 K. The system is further evolved through a third NPT simulation that produces the trajectory from which configurations are extracted. All simulations are terminated at convergence of the volume. The final trajectory is 30 ns long. Simulations were carried out via the NAMD2.12 package 57 with CGenFF parameters 58,59 . Van der Waals cut-off is set at 12 Å and electrostatic interactions are treated according to the particle mesh Ewald method.

Charge dynamics is simulated, for both crystalline and amorphous phases, using the nonadiabatic high-temperature limit of semiclassical Marcus charge transfer theory 42 . The approximation is justified by the fact that charge carriers in organic materials are mostly localised on single molecules that weakly interact with their neighbours 40,41 . The formula for the evaluation of the rate constants reads:

The two main parameters in equation 1, encoding microscopic information, are the reorganisation energy λ, a molecular level property evaluating the energy cost of exchanging a charge carrier, and the transfer integral J ab , that, computed for each pair of exchanging molecules, quantifies the propensity of the carrier to be transferred mediated by the interaction Hamiltonian of the pair.

We note that reorganisation energy in our calculations neglects the contribution due to the change in polarisation of the surrounding environment before and after the exchange and only accounts for intramolecular geometry reorganisation, according to:

where the subscripts indicate molecular state (c or n, charged or neutral) and molecular geometry (C or N, charged optimised or neutral optimised). For sake of coherence, the same DFT approach, B3LYP-D3, has been used for the calculations of reorganization energies. Our previous experience 52,60,61 showed that mobilities computed the B3LYP level are in (very) good agreement with the experiments, even if, generally speaking, theirs values could be affected by the choice of the DFT approach 62 .

Transfer integrals are computed treating the molecular pair as an isolated dimer under the assumption that its orbitals arise from the interaction of monomer frontier orbitals -HOMO for hole transport and LUMO for electron transport -φ i H/L and φ j H/L 63 :

This approximation is justified by the small coupling, between monomers in the pair 34 . All the orbitals needed to estimate the electronic coupling are obtained through the semi-empirical ZINDO method 35 , which has been proven to provide accurate enough estimates 31,32,51 .

Rates are evaluated for all pairs specified in this list. The final step to estimate mobility is to perform kinetic Monte Carlo simulations using the computed rate constants. The mobility is given by the ratio of〈v〉, the average velocity in the direction of the electric field E and the magnitude of the latter.

Owing to the fact that the molecule carries both a DPP electron donor moiety and a fullerene electron acceptor group and that the two act as separate transport channels, the two distinct conducting moieties were defined in order to simulate the charge-separated (CS) state. The separation of the two moieties is performed accordingly to the localisation of the HOMO and LUMO orbitals (see Figure 2). The separation is carried out by cutting the bond between the pyrrolidine in the fullerene moiety and the proximal thiophene ring, as sketched in Figure 1.

The bonds at the cut are saturated with hydrogen atoms and the computational procedure to simulate transport is applied to the morphologies of the separated moieties.

Results and discussion

The most widely applied theoretical framework to describe dynamics in such systems is Marcus theory [40][41][42] . An approach based on the non-adiabatic high-temperature limit of this theory is well justified since the working conditions and components, i.e. high temperature (300 K) and highly localised charge carriers, largely fit into the range of validity of the theory.

The convenience of Marcus theory lies in its ability to describe dynamics in molecular solids as a set of charge transfer reactions with associated rates. Rates depend critically on two parameters (see Eq.1) associated with molecular electronic properties and, more importantly, to bulk structural properties (i.e. crystal structure packing). Their prediction is essential in giving an informed estimate of charge transport performance.

Bulk properties are linked to solid phase morphologies that are the result of the selforganisation of large numbers of molecules. The prediction of solid phase packings, either crystalline or amorphous, remains to this day a challenging task 47 . This is especially true for relatively large organic compounds like dyads that might exhibit polymorphism 43,44 . The reason lies in the heavy computational cost of an accurate treatment of intermolecular interactions that require a high-level description of the electronic structure that includes dispersion forces, a key driver for aggregation in molecular solids 43,[START_REF] Price | Non-covalent Interactions in Quantum Chemistry and Physics: Theory and Applications[END_REF] . In order to overcome this issue, crystal packings are guessed via a less-demanding computational procedure, that employs a classical force field coupled with a Monte Carlo sampling [START_REF] Akkermans | [END_REF] . The output is a set of the most thermodynamic viable crystal polymorphs whose ranking is refined through higher level calculations. The procedure still overlooks other contributions to the ranking that might change it significantly 47 however, the comparison between properties evaluated for these packings and the available experimental data may lead to the reliable identification of the molecular packing responsible for the reported performance. This is of course extremely relevant since it provides direct knowledge of microscopic details that either boost or inhibit transport efficiency of charge carriers. Moreover, the availability of crystal morphologies not directly accessible in the experiments due to their unfavourable thermodynamic rank allows a direct comparison between them. Starting from this, a more general picture can be drawn, thus paving the road towards a set of synthetic principles to be followed in the design of dyads with enhanced performance.

Isolated molecules and crystalline polymorph

A first step has been the characterisation of molecular properties relevant to charge transfer in gas-phase. Previous works 22 have largely stressed how, in dyads, properties of the isolated donor and acceptor moieties should be retained in order to ensure an optimal transport of the charge carriers, meaning that they should remain sufficiently decoupled in terms of their respective electronic structure. A qualitative indication of this decoupling is given by the localisation of the frontier orbitals of the molecule reported in Figure 2, whose energies are reported in Table S1. The strong localisation of frontier orbitals in different parts of the molecule is clear evidence that the moieties are kept largely independent and can efficiently produce a charge-separated state 48 . Indeed, the HOMO is highly localised on the donor moiety, DPP, while the probability amplitude of the LUMO is markedly centred on the Ful moiety. This effect can be traced back to structural features of the dyad. In particular, the possibility of weak interactions of the fullerene with the closest aliphatic chain reduces the planarity of the thiophene rings linking the fullerene moiety to the diketopyrrolopyrrole group.

As a consequence, the π-conjugation is constrained to the remaining part of the donor moiety, producing a certain degree of insulation between the two parts, as can be seen in Figure 2.

This effect is further modulated by the interactions proper to each polymorph packing where it is generally amplified.

Starting from the gas-phase data, reorganisation energies were also estimated according to equation 2. The ambipolar character of the molecule implies that both anionic and cationic states should be investigated. Despite the moieties being virtually independent in terms of transport after charge separation, we consider the response of the whole molecular structure to the presence of a charge carrier. Computed values for the hole (λ(h) = 361 meV) are compatible with those reported for other diketopyrrolopyrrole derivatives, λ(h)~ 300-400 meV 64,65 . At the same time, the energy cost for the exchange of an electron (λ(e) = 212 meV) is close to the reported reorganisation energies for similar fullerene derivatives, λ(e)~130-150 meV for PCBM 66 . This agreement is another strong indication of the efficient decoupling of the two moieties in the dyad.

As already mentioned in the introduction, polymorphs can exhibit dramatically different transport properties. Thin-films employed as active layers usually have a pronounced polycrystalline character, so that it could be expected that different crystal domains, with distinct orientations and possibly packing, will coexist with less ordered regions 67 . In order to provide a reasonable description of the structures present in DPP-Ful thin-films we predicted and analysed the 5 lowest energy crystal packings obtained following the procedure described in the computational details.

In Table 1 we report the unit cell data for the crystal structures identified after full geometry optimization at the DFT level following a preliminary optimization using an empirical force field. Structural parameters for these latter are reported in Table S2, while their structures are reported in Figure S1-S5. The structures belong to the P 2 1 2 1 2 1 symmetry group (structure 0 -i.e. the most stable-, 1, 3, 4) and P n a 2 1 (structure 2). At this level, all the extracted polymorphs are within a range of about 20 kcal/mol from the most stable structure taken as reference. It is interesting to notice that when moving from empirical to DFT calculations a significant change of the ranking is observed, with variations affecting both the energies values and the relative stabilities. It is nonetheless reassuring, that the two most stable polymorphs are still the same according to both approaches, even if their structural parameters are significantly changed (see Table 1 and Table S3). Since we expect the DFT structures to be more accurate than those obtained via the molecular mechanics approach, only the former will be discussed in detail.

The DFT energy differences between the polymorphs considered are rather pronounced.

Based on these values in fact, the probability to populate any of the structures beyond polymorph 1 is absolutely negligible. Nonetheless, one has to consider that in actual thinfilms several polymorphs could co-exist due to the experimental conditions of preparation which could favour the nucleation of a particular structure over another and achieve stability thanks to large kinetic barriers. At the same time, an incertitude of about 3 kcal/mol in the stabilities should be expected 68,69 . As a consequence, all these polymorphs will be discussed in some details, starting from their structural features.

Given the character of ambipolar conductor of the molecule, the most desirable packing should arise from a balance between donor and acceptor competing arrangements, both considerably impacting the magnitude of carrier mobilities. The main factor affecting these latter is the value of transfer integrals, J ij , which depends on the orbital overlap between the two moieties involved in the charge transfer. The relative distance between exchanging moieties, d ij , is the most natural index in term of structure, but other, more complex, geometrical features emerge, among these the angle between conjugation planes in donor moieties and their relative shift.

A detailed analysis of the polymorphs obtained can give more precise examples of these geometric features in the particular case of DPP-Ful, but results can easily be generalised since many characteristics are common to other dyads and, more in general, to organic photovoltaic materials.

The optimized structure of the most stable polymorph (0), belonging to the orthorhombic P2 1 2 1 2 1 space group, features four dyad molecules in the asymmetric unit. In the supercell these molecules are arranged according to a double herringbone in the ac-plane of the crystal, as shown in Figure 3. These two herringbones are identical but have opposite orientation, i.e.

along the positive and negative directions of crystal axis a. The minimum approach distance for donor moiety dimers is d ij = 3.4 Å and it is achieved for consecutive molecules belonging to the same herringbone (Figure 3 lower-left). Layers are then stacked along b, where fullerenes form columns acting as conducting channels for electrons. The distance between fullerenes pairs in consecutive stacking planes is of 3.1 Å (10.3 Å centre-to-centre). The closest contacts are achieved in the stacking direction for fullerenes, as well as in the orthogonal plane (bc) for DPP molecules. The optimized structure of this polymorph is also the one showing, among those investigated, the best agreement with the limited structural data available. In particular, a comparison with the X-ray diffractogram reported in ref. 22 shows that many of the peaks are in reasonable agreement with the simulated XRD pattern for crystal structure 0, as shown in Figure S6. We also remark that the most stable polymorph is not the one with the highest density (see Table 1) suggesting dispersion-driven interactions between DPP-moieties are able to confer a higher stability and that these moieties play a key role in stabilising the packing from a thermodynamic point of view.

The optimized motif of the next most-stable packing, polymorph 1, features a structure in which fullerenes are distributed in a regular lattice in the bc-plane, with the space in between occupied by a close-knitted donor moiety matrix (Figure 4, right). Similar planes are then stacked in the a direction where the closest approach between fullerenes belonging to stacked planes is achieved with a distance of d ij = 2.9 Å allowing rather large values of J ij . The donor structure across the planes is composed of a network of DPP moieties with their backbones arranged almost perpendicularly so that no π-π stacking is observed for DPP moieties in close proximity (see figure 4, left). The packing is rather efficient and confers to the polymorph a density of 1.443 g/cm 3 . The closest approach distance between DPP pairs is rather large, d ij = 4.9 Å, and does not allow efficient transport of holes within the bc-plane nor along a.

Crystal motif of polymorph 2, presents one of the most efficient packings in the set. In the a direction fullerene moieties are arranged again in separated columns that constitute the backbone of two repeated herringbone structures in the ab-plane (Figure 5, right). The minimum distance between a pair of DPP moieties is realised between molecules belonging to neighbouring herringbone structures along c, with a minimum distance d ij of 3.45 Å. The approximate conjugation planes of the pair at the point of closest approach are almost parallel and consequently the orbital overlap is large. DPP pairs in the other directions have reduced contact due to the hindrance of lateral aliphatic chains that keeps them distanced. Fullerene moieties in the backbone of the herringbone structures and stacked along a also have a minimum approach distance d ij of 3.2 Å.

In polymorph 3 we have again a motif composed by two herringbone structures, evident in the bc-plane that are then stacked along crystal axis a. In this case the two rows of molecules making up each distinct herringbone structures are more intertwined (see Figure 6) with respect to other similar structures in the other polymorphs. In this configuration fullerenes are packed close to another fullerene belonging to a neighboring herringbone structure. These pairs are then stacked along a where we have again column-like stacking. All elements in the column reach a minimum distance of 3.4 Å, ensuring rather large values for transfer integrals.

The fact that rows belonging to the same herringbone motif are very intertwined in the bcplane means that donor moieties remain more distant to accommodate for the fullerenes so that DPP moieties cannot reach efficient packing. Minimum distance in the bc-plane is around 4.3 Å and the approach configuration is edge-on, detrimental to transfer efficiency. In the stacking direction a, the hindrance of the aliphatic chains brings the minimum approach distance at over 9 Å, making transport along this direction completely inefficient.

The structure of polymorph 4, is characterised by two herringbone motifs wedged into each other in the bc-plane and then stacked along a, as shown in Figure 7. In this herringbone structure the closest approach for DPP moieties is realised between elements of different herringbone motifs (Figure 7). The minimum distance is d ij = 3.5 Å but the approach is edgeon, i.e. the approximate conjugation planes of the DPP backbones are almost perpendicular to each other. This arrangement lowers considerably the orbital overlap. The stacking along a of the motif in the bc-plane results in the formation of fullerene columns or channels similar to what already observed for polymorph 1, see Figure 4. The minimum approach distance for acceptor moieties in this direction is 3.0 Å, a value small enough to allow large rates of electron transfer between fullerenes.

Polymorphs 2 and 3 boast the largest densities in the set with 1.454 and 1.463 g/cm 3 respectively.

As mentioned, the analysis of the localized structures provides valuable insight about the main features boosting or hampering transport at supramolecular level and relating packing to transport efficiency through the assessment of carrier mobilities will help identify and possibly generalise the main mechanisms at play. Results obtained for the theoretical mobilities for all DFT structures are reported in Tables 2 and3. Data computed using the optimized molecular mechanics structures are reported in Tables S4 andS5 for sake of completeness, but they will be not discussed, since they are less reliable as supported by the worse agreement with the experimental mobilities (see below).

The most stable polymorph (0) is also the one with the best overall performances in terms of both electron and hole transport. Transport takes place mainly along the a-direction for holes while transport across different planes, i.e. along c, is very limited. In the ac-plane the configuration of the closest exchanging pair allows a rather large coupling element, J ij =4.11⋅10 -3 eV, yielding a Marcus rate of 2.04⋅10 11 s -1 . The rather large orbital overlap is ensured by the proximity of the conjugated backbone of one exchanging moiety in to the diketopyrrolopyrrole unit of the second. The overall mobility along a reaches μ a (h)=4.86⋅10 -4 cm 2 /V⋅s, with the herringbone structures acting as main transport channels. Along crystal axis c on the same plane transport is significantly less efficient due to a smaller coupling between DPP moieties belonging to different herringbone structures, J ij = 4.11⋅10 -4 eV giving μ c (h)= 2.84 ⋅10 -7 cm 2 /V⋅s. In the stacking direction, corresponding to crystal axis b, the hole mobility is very low μ b (h)= 4.8 ⋅10 -11 cm 2 /V⋅s. Conversely, along this direction we see the largest value for electron mobility, that takes place along the fullerene channels, μ b (e)=4.99 ⋅10 -4 cm 2 /V⋅s, while being negligible in the other directions.

As remarked above, this packing is not the most dense, yet its structure is able to maximize dispersion interactions by allowing a pronounced interdigitation of DPP moieties that ultimately gives it the highest stability. Comparison with the experimental data, (μ(e)=1.1⋅10 -4 cm 2 /V⋅s, μ(h)=1.3⋅10 -4 cm 2 /V⋅s), reveals that our estimates (μ avg (e)=2.88⋅10 -4 cm 2 /V⋅s, μ avg (h)=2.81⋅10 -4 cm 2 /V⋅s) are able to reproduce the behaviour of the material experimentally observed. Our estimates for both charge carriers are roughly 2.8 times larger than experimental mobilities but the order of magnitude and the relative values between holes and electrons are compatible with the available data 22 . The difference in the absolute values is to be expected given the many effects normally present in a real-life device that can lower the value of mobility with respect to an ideal periodic structure. Real-life devices operate in the centimetre-scale, this means that charges are required to be transported over large distances where other effects (e.g. electrodes, interfaces, higher-scale morphology of the active layer) that lie outside the range of our treatment may affect the overall transport efficiency. cm 2 /V⋅s and μ c (h) = 8.18⋅10 -8 cm 2 /V⋅s, several orders of magnitude lower than the electronic one.

The crystal motif of polymorph 2 reaches the best compromise between the packings of the two moieties and shows very large transport efficiencies for both carriers (μ a (e)=2.16 ⋅ 10 -1 cm 2 /V⋅s and μ c (h)=2.84 ⋅10 -1 cm 2 /V⋅s). Optimal transport occurs in orthogonal directions for holes and electrons, respectively along c and a. This means that electrons move preferentially along the fullerene channels that are formed along a, while holes move preferentially within the plane perpendicularly to the direction of the herringbone structures (see Figure 5). Polymorph 3 offers an example of another limiting case. Fullerenes are optimally packed and can transport electrons efficiently reaching very large mobilities compared to the standard of organic photovoltaics materials, reaching a value of 1.34 cm 2 /V⋅s in the a direction, while in the other two directions transport is absolutely negligible . Donor moieties on the other hand have reduced orbital overlap, with the largest value of the transfer integral being J ij = 4.64⋅10 -3 eV and largest hole mobilities in the order of 10 -6 cm 2 /V⋅s.

Polymorph 4 structure shows low transport efficiencies for both carriers. While also in this case the stacking along a gives rise to channels of fullerenes for electron transport, the larger distance compared to polymorph 1 and the configuration of the pairs with the closest approach translate in a smaller electron mobility, μ a (e) = 1.58⋅10 -5 cm 2 /V⋅s along the direction of the fullerene stacking that is negligible in the other orthogonal directions. DPP pairs in closest proximity have a shortest distance of 3.5 Å which is comparable to the one found in polymorph 0 but the particular configuration of the exchanging pair makes the orbital overlap sensibly smaller (J ij =3.56⋅10 -4 eV) and results in very small mobilities along all directions, the largest being μ c (h) = 6.88⋅10 -11 cm 2 /V⋅s

The comparison we carried out offers an overview of the strongly different behaviours that polymorphs of a same molecule can have and more in general of the fundamental role played by supramolecular effects. Nonetheless, the main goal remains to prove the effectiveness of computational approaches in complementing and interpreting experimental data. From our analysis DPP-Ful reveals a dominant crystalline form that is thermodynamically favoured over the others. This is consistent with data reported in reference 22, suggesting the presence of "long-range ordering" in active layers made with this material.

Our results suggest that the computational approach we implemented is able to reproduce the main microscopic features of the crystal form observed experimentally providing a direct microscopic picture of the molecular packing responsible for the large reported efficiency of DPP-Ful-based devices.

Amorphous Phase

Having discussed crystal polymorphs, let us consider the amorphous phase, which remains interesting for two main reasons. Firstly, amorphous regions may coexist with more ordered ones and contribute to the overall performance of the device. Secondly, completely amorphous films are generally more convenient to produce compared to polycrystalline or crystalline ones 66,70 . The opposite is true in simulations, where the amorphous phase is more challenging to model compared to the ordered one. In fact, the absence of long-range order, makes it formally inaccurate to apply PBC. This requires the production of a significant number of configurations with relatively large amounts of molecules in order to have reliable estimates of the structural effects affecting the mobility (no border effects). Following the approach detailed in the computational details we produced configurations of 500 molecules in cells whose linear size is 106 Å. For all sufficiently close molecular pairs we computed the transfer integrals and the associated rates. This step of the calculation is more expensive in this case. Generally, in ordered phases periodicity limits the number of possible pair configurations and only a small set of rates repeated over the lattice is sufficient. Mobility estimates are computed for every configuration extracted. Disordered systems are to a good extent isotropic, so different directions of the applied electric field can be considered in order to collect some statistics on mobility, results are reported in Table 2. Amorphous mobilities show a marked difference between the two charge carriers.

The reason can be deduced from the analysis of the pair distribution functions g(r) for acceptor and donor moieties. The pair distribution function is defined according to the following:

Where N(r) is the number of molecules with distance r from a given molecule and the denominator is the linear density of an ideal gas with the same total density ρ I of the cell. The features of the distribution signal deviations with respect to the ideal gas behaviour and reveal how the geometry and chemical character of the material may introduce some degree of order at supramolecular scale. We calculated g(r) for both moieties separately.

For each moiety we chose a point from which to compute distances according to their respective characteristics. Fullerene moieties are symmetric enough for the centroid to be a sensible choice while diketopyrrolopyrrole moieties have a more complex geometry. For this reason, we considered as relevant point the centre of charge of the donor moiety in its positively charged state. The pair distribution functions for both moieties are reported in the plot in Figure 8.

By complementing the discussion of Figure 8 with data for the amorphous phase reported in The first peak at around 7.5 Å in Figure 8 is an indication that fullerenes tend to attract each other and form small clusters with a rather efficient packing and high exchange rates reaching 10 11 s -1 . The depletion zone following this first peak is a sign that the donor moieties attached to the fullerenes have more hindrance and require a certain volume around the cluster to be accommodated. This more or less regular alternation of clusters and regions occupied by donor moieties causes the amorphous material to have a poor performance in terms of electron transport as electrons are very unlikely to be able to move across neighbouring clusters and hence remain trapped.

The difficulty of donor moieties in assembling due to the hindrance of aliphatic chains and their linear geometry results instead in a more favourable morphology for hole transport. The structure is sufficiently intertwined and homogeneous to achieve non-negligible transport performances. Electronic couplings between neighbouring molecules have an average magnitude of 10 -3 eV (see Table 2) and the homogeneous structures allows an average transport rate of 10 11 s -1 . These structural features are reflected in the values of mobilities.

Thanks to conformational disorder, holes can achieve mobilities of μ(h)=2.29⋅10 -6 cm 2 /V⋅s on DPP moieties, a value significantly larger than that estimated mobilities for polymorphs 1, 3 and 4 (see Table 2). Conversely, due to the relative ease with which fullerene moieties are able to assemble, clusters severely affect the electron transport efficiency of the amorphous medium with a very low value of μ(e)=2.32⋅10 -9 cm 2 /V⋅s, not comparable with any of the crystal polymorphs studied exhibiting, on average, far larger mobilities for electrons.

Conclusions

Our computational study has allowed us to appropriately characterise, at the molecular level, the packing that is responsible for the transport properties of DPP-Ful active layers reported in ref. 22. The comparison of transport efficiencies for both carriers in the different polymorphs extracted enabled us to point out some general features that might be well applicable to other dyads and hence prove useful in directing future experimental work.

Linking aliphatic chains to donor moieties has been shown to improve solubility and processability of the material 71 but it also makes more difficult for DPP moieties to achieve efficient packing. Despite the possibility for these chains to promote order at ambient temperature due to steric hindrance, our analysis of the predicted structures suggests that they introduce stringent constraints to the packing of donor moieties, considerably hindering the ππ stacking and hence transfer efficiencies between molecules. At the same time, their presence on the diketopyrrolopyrrole unit seems to be promoting the interaction with the fullerene moiety so as to enhance the relative insulation of donor and acceptor moieties by diminishing the π conjugation of the linker. Clearly, these effects might be overestimated by the method used to predict crystal structures. However, even assuming a substantial geometric change of aliphatic chains in the experimental crystal structure their effect could remain important. Indeed, a modulation of their length may very well lead to dyads with better stacking properties of donor moieties while maintaining good stability of the crystal structure at ambient temperature.

In addition, we found that in general fullerene moieties and donor moieties do compete for optimal packing. This means that not all stable polymorphs have the same efficiency in transporting holes or electrons unless an equilibrium is reached against these mechanisms.

This fact results crucially in a highly directional carrier transport that is not, in general, a desirable feature for polycrystalline thin-films, with many randomly oriented ordered regions coexisting between which carriers should be able to cross as swiftly as possible. Moreover, holes and electrons in packings exhibiting high transport efficiencies tend to be optimally conducted in directions orthogonal to each other. This is possibly due to the fact that fullerenes are rather large molecules and, when stacked in columns determine orthogonal planes with considerable spacing between them.

Competing clustering of fullerene and donor moieties emerges also in the amorphous phase.

An indication that the typical hurdles linked to fullerene-derivatives in bulk heterojunction blends are also transferred to molecular heterojunctions.

A number of these features can be extrapolated to a substantial degree to the larger class of organic photovoltaics dyads. This demonstrates the validity of our computational approach that is able to estimate molecular level properties and crucially to complement it with information about packing at the supramolecular level. The latter task remains one of the most delicate steps, in an ideal implementation an appropriate description including disorder and effects of morphologies at higher scales, would allow the simulation of more realistic structures fundamental in providing more accurate predictions of photovoltaic performance.

Polymorph 1

 1 exhibits a much more unbalanced transport behaviour. Electron transport is more efficient than polymorph 0 given the closer packing in the a direction, and again only takes place along this axis. The minimum distance between two fullerenes belonging to different bc-planes is only 2.6 Å. This allows a large transfer integral, J ij =1.68⋅10 -2 eV and a large exchange rate giving a mobility of μ a (e) = 1.58⋅10 -2 cm 2 /V⋅s. Despite the proximity of DPP moieties in the close-knitted donor structure the backbones in the closest exchanging pair cross almost orthogonally, resulting in a limited overlap of orbitals and low transport efficiency across the bc-plane with mobilities μ b (h) = 1.28⋅10-7

Figure captions Figure 1 .

 captions1 Figure captions Figure 1. Chemical structure of DPP-Ful with [60]-fulleropyrrolidine (Ful) acceptor moiety in blue and the diathiafulvalene-functionalised diketopyrrolopyrrole (DPP) donor moiety in red.

Figure 2 .

 2 Figure 2. Isosurface (values +/-0.005) representation of frontier orbitals. (left: HOMO; right : LUMO).

Figure 3 .

 3 Figure 3. Structural features of the most stable polymorph (0) obtained after a full DFT optimization: the shortest distance between fullerene pairs and the distribution of acceptor moieties (top), the shortest distance between DPP pair and a view of the herringbone structures in the ac-plane (bottom).

Figure 4 .

 4 Figure 4. Structural features of the polymorph 1 obtained after a full DFT optimization: the shortest distance between DPP pairs moieties with the largest orbital overlap for the structure (left) and a view of the crystal motif along the bc-plane (right).

Figure 5 .

 5 Figure 5. Structural features of the polymorph 2 obtained after a full DFT optimization: the shortest distance between for a DPP pair (left) and a view of the bc-plane (right).

Figure 6 .

 6 Figure 6. Structural features of the polymorph 3 obtained after a full DFT optimization: the shortest distance between for a DPP pair (left) and a view of the bc-plane (right).

Figure 7 .

 7 Figure 7. Structural features of the polymorph 4 obtained after a full DFT optimization: the shortest distance between for a DPP pair (left) and a view of the bc-plane (right).

Figure 8 .Figure 8

 88 Figure 8. Radial Pair distribution functions for the donor (red) and acceptor (blue) moieties of DPP

Table 2

 2 can give us useful insight in the behaviour of DPP-Ful in disordered morphologies.

Table 1 .

 1 Computed cell parameters (a,b, c, Å ), shortest distance between donor (d h ab , Å) and acceptor groups (d e ab , Å), densities (D, g/cm 3 ) and relative stabilities (E, kcal/mol) for the five different crystalline polymorphs. All reported values refer to polymorphs fully optimized at the B3LYP-D3 level.

	Polymorph	a	b	c	d h ij	d e	ij	Space group	D	E
	0	16.75	9.63	53.95	3.4	3.1	P2 1 2 1 2 1	1.381	0.0
	1	19.26 18.97 22.80	4.9	2.6	P2 1 2 1 2 1	1.441	1.5
	2	18.97 31.92 13.65 3.45 3.2	Pna2 1	1.454	7.2
	3	12.85 40.54 15.77	4.3	3.4	P2 1 2 1 2 1	1.463	7.7
	4	9.68	31.63 32.71	3.5	3.0	P2 1 2 1 2 1	1.200 20.7

Table 2 .

 2 Data on hole transport and electron transport d ij in Å, J ij in eV, k ij in s -1 , μ in cm All reported values are obtained using B3LYP-D3 optimized structures.

	2	/V⋅s.

Table 3 .

 3 Electron and hole mobilities along all crystal axes a,b and c. All reported values are in cm 2 /V⋅s and are obtained using the B3LYP-D3 optimized structures.

	Polymorph	μ a (e)	μ b (e)	μ c (e)	μ a (h)	μ b (h)	μ c (h)
	0	1.00⋅10 -6	4.99⋅10 -4	1.17⋅10 -11	4.86⋅10 -4	4.80⋅10 -11	2.84⋅10 -7
	1	1.58⋅10 -2	7.56⋅10 -11	1.46⋅10 -10	1.98⋅10 -7	1.28⋅10 -7	8.18⋅10 -8
	2	2.16⋅10 -1	2.15⋅10 -11	1.11⋅10 -4	3.80⋅10 -7	7.04⋅10 -10	2.84⋅10 -1
	3	1.34⋅10 0	8.90⋅10 -13	1.64⋅10 -9	8.47⋅10 -10	9.09⋅10 -10	2.28⋅10 -6
	4	1.58⋅10 -5	3.08⋅10 -11	4.51⋅10 -11	5.68⋅10 -11	1.72⋅10 -10	6.88⋅10 -11

Acknowledgements

We thank Dr Frédéric Labat for his help with the CRYSTAL code. We The authors acknowledge the use of computational resources at Tr s rand entre de alcul (TGCC) provided by GENCI through project A0050810135 and the CINECA computing facility for the availability of high-performance computing resources and support. This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme ( o STRI ES; grant agreement No 648558)


Iman Hmedoush
email: iman.hmedoush@inria.fr

Cédric Adjih

Paul Mühlethaler

Lou Salaün
email: lou.salaun@nokia-bell-labs.com

Multi-Power Irregular Repetition Slotted ALOHA in Heterogeneous IoT networks

Irregular Repetition Slotted Aloha (IRSA) is one candidate member of a family of random access protocols to provide solutions for massive parallel connections in the Internet of Things (IoT) networks. The key features of this protocol are repeating the transmitted packets several times and using Successive Interference Cancellation (SIC) at the decoder to resolve the collisions, which dramatically increases the performance of Slotted ALOHA. Motivated by multiple previous studies of IRSA performance in different settings, we focus on the scenario of an IoT network where the packets of different nodes are received with different powers at the base station, either per design due to different transmission power, or induced by the fact that the nodes are at different distances from the base station. In such a scenario, the capture effect emerges at the receiver, which in turn enhances the protocol performance. We analyze the protocol behavior using a new density evolution which is based on dividing nodes into classes with different powers. By computing the probability to decode a packet in the presence of the interference, we explore the achievable throughput and its associated gain and show the excellent performance of Multi-Power IRSA.

I. Introduction

A. Modern Random Access Protocols

A scenario of a typical IoT network is to connect embedded devices, such as environmental sensors, to a central collection node. In many IoT scenarios, there will be a huge inflow of data from devices sharing the same radio spectrum. One of the main challenges of IoT networks is to assign the network resources such that we ensure the connectivity for all devices with a certain rate depending on the required application.

Non-Orthogonal Multiple Access (NOMA) is one of the proposed solutions to solve the issue of massive connectivity in IoT. NOMA often operates at the symbol level, however, an interesting alternative is a family of protocols [START_REF] Paolini | Coded Random Access: Applying Codes on Graphs to Design Random Access Protocols[END_REF], frequently known as the Coded Slotted Aloha (CSA) family, which has introduced NOMA features to classical ALOHA. They can be considered as modern random access protocols. The common idea of these protocols is to allow each active terminal to send multiple copies of its transmitted packet with the same data and preamble but also to include in each copy an indicator to the positions of the other copies in the frame. On the other side, the receiver will use Successive Interference Cancellation (SIC) to perform an iterative decoding process which removes the successfully decoded packet from all the positions of This work was partly partly supported by the Inria-Nokia Bell Labs common labs under research action "Network Information Theory". its copies. Contention Resolution Diversity Slotted Aloha (CRDSA) [START_REF] Casini | Contention resolution diversity slotted aloha (crdsa): An enhanced random access schemefor satellite access packet networks[END_REF] is one variant of CSA family where each packet is repeated only twice, while in Irregular Repetition Slotted Aloha (IRSA) [START_REF] Liva | Graph-based analysis and optimization of contention resolution diversity slotted aloha[END_REF], the packet can be repeated a different number of times by each active terminal (i.e. based on a probability distribution function). For the throughput metric, e.g. the number of received packets per slot, it has been shown in [START_REF] Liva | Graph-based analysis and optimization of contention resolution diversity slotted aloha[END_REF], that IRSA can achieve a throughput of T = 0.97 for very large frames (infinite length) and around T = 0.8 in practical implementations, and [START_REF] Krishna | Iterative collision resolution for slotted aloha: An optimal uncoordinated transmission policy[END_REF] proved that T → 1 with soliton distributions. Previous studies have explored the question of: How can we enhance the performance of IRSA?

One way to increase the achieved throughput of this protocol is through capture: captures occur when two (or more) packets are received with sufficiently different power. The packet with the highest power is recovered first, and then the second after SIC. This paves the way for potentially achieving T > 1.

B. Capture, related work and our approach Multi-Power IRSA

In a heterogeneous setting, where IoT nodes are placed in different positions from the base station, the different path loss factors experienced from different channel conditions naturally give rise to the capture effect at the receiver. The effect of such path loss in wireless communication has been extensively studied. By developing stochastic geometry methods, in [START_REF] Franc | Stochastic Geometry and Wireless Networks, Volume I -Theory[END_REF], the received interference and capture probability from a set of nodes in a uniform wireless network is computed. Other studies [START_REF] Zhang | The performance of successive interference cancellation in random wireless networks[END_REF] have considered the case of capture effect with the property of SIC in the receiver. In [START_REF] Khaleghi | Near-far effect on coded slotted aloha[END_REF], a study of Irregular Repetition Slotted Aloha (IRSA) in the presence of capture effect and SIC was presented: there, the transmission power is identical for all users and the difference of received power is entirely caused by path loss from distance.

Several schemes where the transmitter power is changed at each transmission have been introduced. In [START_REF] Shao | Noma-based irregular repetition slotted aloha for satellite networks[END_REF], a NOMA-Based IRSA scheme with different transmission powers is proposed and typical Density Evolution is used for studying the system while the Differential Evolution method is used for optimizing the parameters. In [START_REF] Alvi | Enhancing crdsa with transmit power diversity for machine-type communication[END_REF], CRDSA with transmit power diversity has been considered. The transmission power distribution was optimized by differential evolution with respect to a total power budget to enhance the throughput.

Our approach called Multi-Power IRSA (MP-IRSA) differs from previous work by introducing the important differences Fig. 1: Example of IRSA, with Λ 2 = 2 3 and Λ 3 = 13 -and of the decoding process that the replicas of the same packet have the same power, and that the users are grouped into classes corresponding to identical received power. MP-IRSA is well-suited to heterogeneous IoT networks, e.g. the signals of farthest terminals arrive at the base station with lower strength, and they might be grouped in classes with lower power: their adjusted transmission power stays low, providing energy savings. MP-IRSA allows the phenomenon of cascading decoding, can increase throughput.

We analyze MP-IRSA through a new variant of multi-class density evolution that relies on the grouping of users with the same received signal power into one class. With that, we are able to study experimentally the performance of the protocol: the influence of each user class on the decoding process and the performance of the other classes are also studied. The impact of the density and the power of the class on the achievable throughput and the maximum load is analyzed.

The paper is organized as follows: Section II introduces the concepts of IRSA with capture and multiple classes. Section III explains the new density evolution of IRSA with multiple classes. Later in section IV, the interference model inside the class is derived. The section V presents our numerical results and finally, section VI concludes our work.

II. Principle of IRSA with capture and multiple classes

We consider K classes of N c users in a class c ∈ K. Each class contains users whose signals are received with the same power by the base station. Each user has one single packet to transmit (all packets have identical sizes and fit a slot). One considers a frame of size M slots. Each user transmits identical copies of its packet on randomly chosen slots. The number of replicas is picked by each user at the beginning of the frame from a random discrete distribution Λ c which is common for all users in our assumption. Precisely, we define: Λ c,i Pr(a user of class c repeats his packet i times).

At the end of the frame, the receiver receives the superposition of the physical signals sent by collided users on each slot. We abstract the collision model as: consider a given slot, where n 1 + n 2 + . . . + n k replicas that had been sent from k classes, with n i being the number of users of class i transmitting on that slot. Let us denote the number of replicas of class i that can be recovered in the presence of the interference of other classes as D i (n 1 , n 2 , . . . n k ). Here, we also adopt an intra-class collision model, i.e. whenever n i ≥ 2, no packet from class i can be recovered on the slot. Packets from other classes, however, might sometimes still be recovered when their received power is higher.

The decoding process is iterative [START_REF] Liva | Graph-based analysis and optimization of contention resolution diversity slotted aloha[END_REF], and we consider the following rules while performing the decoding of class i:

• We consider all slots containing exactly one packet from class i. In each of these slots, this packet can be recovered if the interference from other classes is sufficiently small (lower than a given threshold). • Once one packet has been recovered on a certain slot, the copies of this packet can be removed from other slots at the physical signal level by Successive Interference Cancellation (SIC). A simple example is given in Fig. 1a (with Tanner graph in Fig. 1b: bipartite graph with users and slots, with one edge for each transmission), and the detailed explanation of the rules of the decoding process is shown in Fig. 1. This figure shows two classes with 3 nodes in each class where the received power of class 2 is lower than the received power of class 1. The transmissions from the two classes are colliding on the slots: 2, 3 and 4. The decoding process starts with class 1 (Fig. 1c) as follows: the intra-class decoding starts first by searching for the packets without collisions (singletons). The packet of user A of class 1 is a singleton received on slot 1, so it is decoded and removed from its replicas positions on slots 2 and 3 as shown in Fig. 1d and Fig. 1e. The removal of the packet from user A on slot 2 makes the packet of user B on slot 2 colliding only with a replica from user D from class 2. Since class 2 received power is considerably lower than the received power of class 1, the packet of user B could be considered as a singleton in this case with an "underthreshold" interference from class 2, and it is in turn decoded and removed from slots 2 and 3 as well as shown in Fig. 1f and Fig. 1g. Removing the packet of user B from slot 3 makes the packet of user C decodable even in the presence of a replica of user E from class 2. Note that the packets of class 1 can be removed in the presence of under-threshold interference from class 2 but the opposite is not true. The same process of iterative decoding of class 2 continues (Fig. 1k) by removing the packets of users D and E on slots 4 and 5, as shown in Fig. 1l and Fig. 1m and recovering user F in Fig. 1n. Note that in this example, performing the decoding of class 2 packets is impossible before the decoding of class 1 packets (Fig. 1j).

In classical IRSA [START_REF] Liva | Graph-based analysis and optimization of contention resolution diversity slotted aloha[END_REF], there is only one class. In our work, various decoding order strategies can be adopted depending on when we switch on decoding packets from one class, to decoding packets from another class. When all packets have been retrieved, or no packet from any class can be decoded anymore, the decoding process is considered to be ended. Note that all the strategies will lead to the same final state.

III. Density Evolution

The use of the density evolution method for IRSA was pioneered by [START_REF] Liva | Graph-based analysis and optimization of contention resolution diversity slotted aloha[END_REF] (and then developed in most of the subliterature to this family of protocols). It models stochastically the decoding process (as in Fig. 1), by computing the expected amount of decoded packets at each iteration. It allows us to find the throughput obtained for a given load. We provide a new density evolution formulation for MP-IRSA and use it later for our performance analysis. As for LDPC decoding, the following function is introduced to represent the repetition degree distribution for class c: Λ c (x) = l i=1 Λ c,i x i . It is also important to define the edge degree distributions:

λ c (x) = i=1 λ c,i x i-1 , ρ(x) = i=1 ρ i x i-1
where: λ i,c : defines the probability that a randomly selected edge from the tanner graph in Fig. 1b is connected to a user node of class c which has a degree i. ρ i : defines the probability that an edge is connected to a slot node of degree i.

This method is a variant of the message passing algorithm [START_REF] Mansour | A turbo-decoding message-passing algorithm for sparse parity-check matrix codes[END_REF], where the idea is to perform iterative decoding:

• knowing the average probability that each user's packet has been recovered at the iteration , this allows computing the average probability that each slot contains a packet that can be recovered at iteration ; and conversely, • knowing the probability that a slot contains a packet that can be recovered at iteration , assuming actual decoding followed by SIC, this leads to a new probability that each user's packet is recovered at the end of decoding iteration and the beginning of the iteration + 1.

Because it is a message passing algorithm, it is more convenient to define, use, and update average probabilities on edges, as done in the LDPC literature. We consider a random edge corresponding to one given user of class c transmitting on a random slot:

• q c,i is the probability that the edge which is connected to one user of class c is not known. • p c,i is the probability the edge that connects a packet of class c to a random slot is not revealed. where: i is the iteration number. The last packet from class c on a slot can be recovered if all the packets from the same class and the classes with a higher power have been already decoded in the previous iterations:

1 -p c,i = δ c 1 -q c,i l-1 (1)
Where: δ c is the probability that a packet from class c can be decoded in the presence of interference from other classes, and:

δ c = f c q 1,i , q 2,i , ....., q c-1,i (2)
Where f c is the function of the probability of the unknown edges of all classes which will be clearly defined in eq. ( 7) in the next section. A packet from class c can be recovered if there is at least one copy of this packet that has been decoded on another slot:

q c,i = p l-1 c,i-1 (3)
Using the polynomial representations of user degree and edge degree distributions, we can average the edge probabilities for class c in ( 1) and (3) as follows:

p c,i = 1 -δ c • ρ 1 -λ c p c,i-1 l-1 (4)
By letting M → ∞ (Asymptotic case), we can write:

p c,i = 1 -δ c e -Gc Rc λ c (p c,i-1 ) (5)
Where: G c is the load of class c which can be defined as the average number of class c users per slot, i.e, G c = N c M . R c is the rate of class c and is defined as R c 1 Λ c (1) . The necessary condition to decode more packets in each class at each decoding iteration is : p c,i < p c,i-1 which can be written more precisely using (5) as follows:

1 -δ c e -Gc Rc λ c (p c ) < p c (6)
Given the parameters (G c , Λ c , δ c ) for each class c ∈ K, one can simply follow the evolution of the decoding process to ultimately estimate the number of decoded packets at the end of the decoding process and also to define the best degree distribution to be given to each class in order to achieve the highest load.

IV. Inter-class Interference Model

In a heterogeneous IRSA setting, where capture effect is considered, the analysis of the decoding process of a packet from a given class c has to take into account two types of packets: the packets from other classes ∀n ∈ K and n c (considered as interference) that collide on the same slot and also collisions from the same-power packets (same class c). Thus, the density evolution of a class-based IRSA couples two interference models: intra-class interference that can be studied and processed as in classical IRSA in [START_REF] Liva | Graph-based analysis and optimization of contention resolution diversity slotted aloha[END_REF] and interclass interference that appears in the factor δ c for each class. In case of capture at the receiver, when the difference between the strong and the weaker signals is sufficiently large, the weaker signals can be canceled before decoding the strong signal. This gives a condition on the inter-class interference on each slot which affects the decoding process and it can be expressed as:

δ c = Pr                        K n = 1 n c X n Π n ≤ T h                        (7)
where: Π n is the received power of class n, X n is a random variable representing the number of undecoded packets from class n on the slot, and T h is the interference threshold (SINR) beyond which a packet can be decoded. Indeed, (7) integrates the interference with the effect of path loss and fading on the received power which has been widely covered in wireless networks. Notice that ( 7) is linked to stochastic geometry [START_REF] Franc | Stochastic Geometry and Wireless Networks, Volume I -Theory[END_REF], in particular when the received power directly depends on distance but also depends on the stage of the decoding process. We can write that X n follows a Poisson distribution with a mean B n which represents the average number of undecoded users from class n on the slot:

B n = λ n (q n )N n l Λ n,l l M (8)
Where: λ n (q n ) represents the probability that an edge on class n has not been revealed and N n is the total number of users in class n. The sum l Λ n,l l, represents the average number of replicas of a packet in class n. The left side of the inequality in ( 7) is a sum of scaled Poisson random variables which generally has not a closed formula for its distribution. However, a good approximation can be computed numerically.

V. Numerical Results

We mainly focus on the asymptotic performance of MP-IRSA using our density evolution which is based on classes. The primary goal is to show that isolating the users with the same received power in a class and taking into account the interference from other classes has a strong and notable impact on the decoding process, since the classes with high received power have more probability to be decoded first, and followed by decoding lower-power classes and so on. Indeed, results evidence a cascading effect on the decoding process starting from the highest-power class to the lowest-power class. Our second objective is to compare the throughput of MP-IRSA with the optimized multi-power CRDSA which was proposed in [START_REF] Alvi | Enhancing crdsa with transmit power diversity for machine-type communication[END_REF] and also to confirm that it has higher than the throughput of classical IRSA in [START_REF] Liva | Graph-based analysis and optimization of contention resolution diversity slotted aloha[END_REF]. This emphasizes the importance of favoring captures into the system design. Our third objective is to explore the impact of the various parameters.

In our numerical experiences, we study the influence of the number of classes, the density of each class (defined as the proportion of the nodes in that class), and the power of each class on the achievable throughput and its associated gain, in different scenarios. The users' degree distribution for IRSA in our analysis is a soliton distribution from [START_REF] Krishna | Iterative collision resolution for slotted aloha: An optimal uncoordinated transmission policy[END_REF]. We gradually increase the load (x-axis of the curves) and plot the corresponding throughput for each class. The total throughput of the network will be the sum of the throughput of all classes. Fig. 2 shows the throughput achieved by IRSA for 2 equal density classes and different powers. We tested three different cases with a power difference between both classes corresponding to 10 dB, 20 dB, and 30 dB respectively. Typical of IRSA, the curve starts to grow linearly (with near 0% loss) until a threshold, after what the throughput falls quickly.

For a sufficiently large power difference, i.e. with P = 100 (+20 dB) or with P = 1000 (+30 dB), the maximum achieved throughput is the near-maximal value: T max = 0.894 for the first class and T max = 0.904 for the second class. This means that the total throughput is 1.798 = 0.894 + 0.904, a near-doubling of the 1 packet/slot bound without capture. This illustrates from the cascading effect of MP-IRSA with high power difference: the decoder can finishes decoding the packets of the highest power class (intra-class decoding) almost entirely as if the other class did not exist, and then moves to the next class (inter-class decoding).

This phenomenon cannot occur with smaller power differ- Fig. 4: The throughput of IRSA with 6 classes and geometrical difference in power and equal densities and SIR = 1 ences, and this is illustrated for P = 10: where T max = 0.673 for both classes: the decoding process of class 2 can be blocked by class 1, and the receiver probably will have to decode by cycling back and forth between the classes. The maximum total throughput is 1.346, significantly less than before. Fig. 3 shows the performance of 2-classes with CRDSA with 2, 3, and 4 repetitions respectively. The goal of this comparison is to confirm that IRSA can perform better than CRDSA, as in the classical case, as was shown in [START_REF] Liva | Graph-based analysis and optimization of contention resolution diversity slotted aloha[END_REF][START_REF] Krishna | Iterative collision resolution for slotted aloha: An optimal uncoordinated transmission policy[END_REF], and confirm this even in the case of classes. With the power of classes set to (1, 100), the best degree for CRDSA is 3 as in IRSA without capture. For that value, IRSA indeed still performs better with than CRDSA which achieves a maximum throughput equal to 0.814. The cascading effect for CRDSA degree 3 and 4 as for IRSA can be guessed.

To understand more the impact of the power difference between the classes on the system performance, we study the case of 6 classes with a geometrical difference in power and for different SIR levels, as in Fig. 4. One of the key factors that affect system throughput is users differentiation. As before, when the classes have more difference in power, the obtained throughput can be higher, and this is related to the fact that the created interference from one class to another is more tolerable in the SIC receiver as seen in Fig. 4a, Fig. 4b, Fig. 4c and Fig. 4d (for clarity the result of only 5 classes out of 6 are plotted). Interestingly, for a geometric factor of 6, the total throughput with 6 classes is almost 5 (> 0.8 × 6). This implies that it is almost if there were 6 cascading decodings, with more than 80% successfully used slots for each class. For smaller geometric factors, the total throughput is less falls dramatically, and adding more intermediate classes seems to sometimes decrease performances.

Another important key factor is receiver sensitivity (i.e, the SIR). Fig. 4d shows a throughput of T max = 0.966 in case of S IR = 1, while the throughput degrades almost to the half T max = 0.469 when the SIR is 4 times greater. This can be seen easily by looking at a packet from the class 3 and one of its interferers in both cases. In case of S IR = 1, the collision between 8 packets from class 2 and one packet from class 3 can be still decoded since S IR = 512 8×64 = 1, While in case of S IR = 4, the packet of class 3 can be decoded in case of interference with one or two packets at maximum from class 2 which is a degradation by a factor of 4.

Fig. 5 shows the importance of the number of classes on the achieved throughput and the associated gain. The performance degrades notably when there are more classes in the network as in Fig. 5c due to the increase in the interference between the classes. The total throughput is around 2.7, 2.0, and 1.6 respectively. It appears as if the intermediate classes were blocking the cascading effect. Thus, the number of classes should be chosen after taking into account different parameters including the received power range and the needed SIR.

Another critical factor is the users' densities in each class. Fig. 6 shows different scenarios with different users' densities. Recalling that the decoding process starts from the highest power class, Fig. 6a illustrates what occurs when the highest power class has also the highest density. The effect of the class density is not clear in this case: even though the power ratio between classes is 10, lower power classes are often blocking high power one (as in Fig. 2), this occurs less often because of their lower density. The impact of class density is clearer in Fig. 6b and Fig. 6c, where the highest power class (in yellow) is more affected by the interference from the other classes (87.5% and 95% of the users are interferes with class 1). As Fig. 6: The effect of different densities of the classes shown in both figures, the decoding of the highest power class continues correctly (linearly) up to a certain load, and then it degrades slowly affected by the high interference from other high-density classes. Two observations can be made: in those cases, for a higher value of the load, the lower power classes unnecessarily "jam" the channel, because their packet success rate goes to 0; but on the other hand, when the degradation of the throughput of higher classes is not reached, the system is interesting as it introduces a form of priority between users.

VI. Conclusion

In this paper, we formally introduced multi-power IRSA (MP-IRSA), as a random access method, when the replicas of different users are transmitted/received with different powers. We introduced a new density evolution variant based on grouping users into classes: it allows analyzing the performance of MP-IRSA. Multi-power IRSA proves to be a better choice than CRDSA and classical IRSA and confirms the benefits of effects of the capture effect at the receiver. The impact of different system factors on the achieved throughput and the associated load was extensively studied. First, the power difference between the classes and the number of classes plays a huge role in the decoding process, and we observe the best results when the power difference is large: decoding is cascading (one class after the other), instead of cycling (switching between different classes). The receiver sensitivity and the density of each class are other important factors that we studied. Another key factor that we did not explore is the repetition degree distribution: we used the same soliton distribution for all classes. Finding methods to optimize a common degree distribution for all classes or multiple degree distributions for multiple classes, and studying its impact, is a possible future work of this paper.

Fig. 2 :

 2 Fig. 2: The throughput of 2 classes with different powers and identical density 50%. Each color is for one scenario with 2 classes.

Fig. 3 :

 3 Fig. 3: The throughput of 2 classes of CRDSA with 2, 3 and 4 repetitions with different powers and similar densities

  6 classes with a geometrical difference of 8 in power and SIR = 1 6 classes with a geometrical difference of 8 in power and SIR = 2 6 classes with a geometrical difference of 8 in power and SIR = 4

Fig. 5 :

 5 Fig. 5: The effect of increasing the number of classes in the same power range


domain: phys spi

text:
Anne Talneau
email: anne.talneau@c2n.upsaclay.fr

Grégoire Beaudoin

Gilles Patriarche

Frédérique Ducroquet

Efficient Electrical Transport Through Oxide-Mediated InP-on-Si Hybrid Interfaces Bonded at 300 °C

Keywords: 3D hybrid integration on Silicon, Electrical transport, Semiconductors hybrid interfaces

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In the domain of Silicon photonics, the ability to electrically operate the hybrid III-V on Si (or SOI) interface will open the full 3D field bringing new design opportunities and immense gains for thermal-budget, chip power consumption, footprint and integration of driving electronics. Electrical transport has been demonstrated in the case of GaAs/Si [1][2] ; it has also been demonstrated in the pioneer work of D.Pasquariello [3] in the case of both hydrophobic and hydrophilic InP/Si bonding , and recently for air-bonded InP/Si [4] . Hybrid heterostructures obtained by bonding have the advantage versus heteroepitaxial growth to offer a high crystalline quality, even though the lattice mismatch is important, being ~8% in the case of InP/Si. We have here considered bonding processes being CMOS processing compatible. Bonding has then to be operated at a temperature below or at 300°C maximum in order to preserve the already processed contacts. Thus an intermediate layer is mandatory, otherwise wafer fusion without intermediate layer is reported to be operated in the 600°C temperature range, in the case of InP/GaAs [5] or hydrophobic bonding is reported to be operated at 400°C in the case of InP/Si [3] .

The intermediate layer commonly introduced is a SiO 2 layer, due to its large compatibility with Silicon. It can be an intentionally deposited SiO 2 layer with a thickness in the [5-50] nm range; this approach has been widely implemented for hybrid devices [6][7][8][9] but such a large thickness prevents electrical transport through the interface. When bonding is operated without any intentionally deposited layer but under ambient air, this processing leads to a non-controlled thickness of the oxide layer at the interface, which is in the [2-5] nm-thick range in the case of [1] , or which thickness is not given in the case of [4] . For such a thin layer, electrical transport has been measured thanks to tunneling through the oxide layer. In order to generate an ultra-thin SiO 2 interfacial layer, we have developed a process operating under vacuum with de-oxidized surfaces activated by ozone, reaching a 1.2nm-thick oxide interface layer [10] .

Hybrid interface

We consider here bonding operated at 300°C, thus including an interfacial layer. We have investigated two interfacial layers: (i) the SiO 2 interfacial layer generated by ozone-activated bonded surfaces and (ii) an intentionally deposited thin ZnO interfacial layer, taking advantage of the conductivity of this metallic oxide [11][12] .

In order to clearly evidence the degradation of the electronic transport that could be brought by the intermediate oxide layer, we have considered a hybrid heterojunction with doping polarity and concentration leading to an ohmic behaviour without any oxide. The n-Si 10 19 cm -3 /n-InP 10 19 cm -3 doping polarity and concentration heterostructure has already experimentally demonstrated an ohmic behavior [13] ; this choice is also supported by the band diagram calculated with the CHARGE solver of the Lumerical commercial tool. Figure 1 indicates that an ohmic behavior can be expected due to the reduced potential barrier at the conduction band edge. The energy reference Fermi level is taken at 0 eV. Additionally, the high doping of both materials ensure a negligible contribution of both material sections to the measured resistance.

Oxide-mediated bonding

Both 1cm square InP and Si samples are de-oxidized, terminated by HF last. Both surfaces are then activated by ozone -O 3 -and annealed under vacuum at 300°C during 1 h 30 min, with a slight pressure. Such a process has produced a 1.2nm-thick SiO 2 interfacial layer [10] . Figure 2-left shows the STEM image of the interface, evidencing no structural defect in both lattices. We have also considered a conductive oxide intermediate layer, since even if its thickness will be several nanometers, its conductivity would not hamper the electrical transport. ZnO is well known to present a conductivity related to oxygen vacancies [11][12] . Both InP and Si surfaces are prepared oxide-free as described in [10] , a ZnO layer 1.6 nm-thick has been deposited by sputtering at room temperature on the InP surface only. Both surfaces are then immediately brought into contact and annealed at 300°C under vacuum during 1 h 30 min, the same bonding procedure used for ozoneactivated surfaces [10] . Figure 3-left shows the STEM image of the interface, evidencing no structural defect in both lattices. Figure 3-right shows in-situ EDX profiles for all the elements. The Zn element distribution has a 15% maximum, showing a large inter diffusion both in Si and in InP. This inter diffusion is enabled since the sputtering deposition at room temperature produces an amorphous layer of low density which easily diffuses during the 1 h 30 min duration annealing at 300°C for bonding.

The oxygen element is present on a distance of ~ 4nm, larger than the 1.6nm-thick ZnO deposited layer, and on the Silicon element side mainly, indicating that silicon has been oxidized, generating the stable oxide SiO 2 . The 4nm-thick oxide interface layer is then a sum of SiO 2 + Zn(Si)O x +ZnO.

Electrical processing and measurements

For electrical transport measurement, a n-doped (S) 1x10 19 cm -3 500nm-thick InP membrane has been bonded on a n-doped (P) 1x10 19 cm -3 Si substrate. Both oxide bonding, ozone-activated surfaces and ZnO intermediate layer, as well as an oxide-free bonding, have been processed and measured, the oxide-free interface being a benchmark since it has already evidenced efficient electrical transport [13][14] . After substrate removal, the InP membrane evidences so surface defect. 150µmdiameter mesa were etched by Cl 2 -based Inductive Coupled Plasma etching through the InP membrane down to the Si substrate. Electrical contacts were performed, a Ti-Au contact on top of the InP mesa, and a Ni/Au/Ge contact on the Si substrate. These two contacts have been qualified to be ohmic. The Si contact geometry has a circular pattern 100µm wide and 110µm far from the InP mesa as sketched in Figure 4 left. Electrical measurements were conducted using a four-probe set-up allowing the exclusion of probe resistance, the current being forced from the mesa contact to the outer contact. V(I) transport characteristics of the heterojunctions (figure 4 right) show ohmic behavior in the case of the oxide-free bonded interface, as already reported, and also for the oxide-mediated bonded interfaces.

In the case of the SiO 2 1.2nm-thick interface layer obtained with the ozone activation of de-oxidized surfaces, this ohmic behavior confirms the tunneling transport already proposed as the transport mechanism in the GaAs/Si interface [1] , and corroborates the measurement performed on InP/Si airbonded interfaces [4] since, even if no TEM picture is provided to clearly show what exists at this interface, an oxide layer is most probably present due to the bonding being operated in ambient air.

Related to this tunneling mechanism, the thinnest oxide layer is sought since the current transport is a decaying exponential function of the oxide layer thickness.

The ZnO bonded interface also shows an ohmic behavior comparable to the O3-activation bonded interface leading to the SiO 2 intermediate layer. Even being thicker (~4nm-thick), the multi composition of this oxide layer (SiO 2 + Zn(Si)O x +ZnO) evidenced by the in-situ EDX profiles of the elements, allows tunneling through its SiO2 part, and conduction through its Zn-based part which is Si and In-doped due the inter diffusion. As a consequence, the ZnO bonded interface has an electrical behavior very similar to the one including the SiO 2 interfacial layer obtained by O 3 -activated surfaces.

The major contribution to the resistance being the interfacial resistance, when compared to the oxidefree bonded interface, both oxide-bonded heterostructures show no significant increase of the resistance related to the presence of the oxide layer at the interface.

The resistivity calculation is beyond the scope of this paper, since its correct calculation would require a current flow simulation to quantify the electric field lines according to the structure and contacts geometry and a quantitative measurement of any possible surface currents along the etched mesa walls.

Discussion

Even though these experimental measurements demonstrate that a very thin oxide intermediate layer does not increase the interfacial resistance when compared to an oxide-free interface, thus allowing efficient electric transport, the proposed involved mechanisms remain macroscopic.

Commercial simulation tools implementing energy gaps and electronic affinities for the band structure and transport calculation does not allow a detailed description of the interface. They also do not take into account that Si and InP do not have the same crystallographic system and the same polarity. A quantum transport approach simulation should be undertaken in order to correctly describe the complex mechanisms occurring at the interfaces (e.g. bond breaking and charge transfers) and include a detailed description of the oxide interfacial layer, especially when the composition of the oxide layer becomes complex due to diffusion during bonding and its related annealing. Advanced simulation implementing full-band quantum simulation with empirical pseudopotential and Non Equilibrium Green Function formalism [15] is under development to support an optimized design of the interfacial layer for efficient electrical transport.

Conclusion

In this study, we have demonstrated that an oxide layer, being thin enough, does not hamper the electrical transport through a n-InP/n-Si hybrid heterostructure. Bonded hybrid interfaces obtained under vacuum operation at 300°C including or a SiO 2 1.2nm-thick interfacial layer or a 4nm-thick complex SiO 2 + Zn(Si)O x +ZnO interfacial layer have been measured, evidencing an ohmic behavior and no significant increase of the interfacial resistance when compared to an oxide-free bonded interface. Oxide-mediated bonding has thus two advantages, a technological one related to the efficient and easily implementing bonding processing, and a functional one allowing electrical transport through the interface, provided the right choice is made for the thin oxide interfacial layer.

In-depth electrical characterization, involving low-T measurement, is under progress, and a quantum transport approach simulation is currently undertaken in order to correctly explain the conduction mechanisms involved.

Figure 1 :

 1 Figure 1: Band diagram (Conduction Band and Valence Band) calculated by Lumerical CHARGE solver . Doping are n-Si/ n-InP both 10 19 cm -3

Figure 2 -

 2 right shows insitu EDX profiles for all the elements, Oxygen is present at the interface, without any diffusion in either lattice. The residual roughness of both surfaces being in the [0.3-0.4] nm range, such a thin layer could be very close to the thinnest possible layer for oxide-mediated bonding. Annealing experiments operated up to 500°C have evidenced the layer stability, without any migration or reordering.

Figure 2 :

 2 Figure 2 : O 3 activation bonded interface left: STEM image of the interface right : in-Situ EDX all elements distribution[9]

Figure 3 :

 3 Figure 3: ZnO-mediated bonded interface left: STEM image of the interface right: in-Situ EDX all elements distribution

Figure 4

 4 Figure 4 Left: Schematic of the electrical contacts (not at scale) Right: V(I) Transport characteristic of the hybrid interfaces