Frédéric Marazzato
email: marazzato@lsu.edu

Alexandre Ern
email: alexandre.ern@enpc.fr

Laurent Monasse
email: laurent.monasse@inria.fr

Quasi-static crack propagation with a Griffith criterion using a variational discrete element method

A variational discrete element method is applied to simulate quasi-static crack propagation. Cracks are considered to propagate between the mesh cells through the mesh facets. The elastic behaviour is parametrized by the continuous mechanical parameters (Young modulus and Poisson ratio). A discrete energetic cracking criterion coupled to a discrete kinking criterion guide the cracking process. Two-dimensional numerical examples are presented to illustrate the robustness and versatility of the method.

Introduction

Discrete element methods (DEM) are popular in the modeling of granular materials, soil and rock mechanics. DEM generally use sphere packing to discretize the domain as small spheres interacting through forces and torques [START_REF] Labra | High-density sphere packing for discrete element method simulations[END_REF], but the main difficulty is to derive a suitable set of parameter values for those interactions so as to reproduce a given Young modulus E and Poisson ratio ν at the macroscopic level [START_REF] Jebahi | Discrete element method to model 3D continuous materials[END_REF][START_REF] Celigueta | Accurate modelling of the elastic behavior of a continuum with the discrete element method[END_REF]. Advantages of DEM are their ability to deal with discontinuous materials, such as fractured or porous materials, as well as the possibility to take advantage of GPU computations [START_REF] Spellings | GPU accelerated discrete element method (DEM) molecular dynamics for conservative, faceted particle simulations[END_REF]. A first DEM parametrized only by E and ν has been proposed in [START_REF] Monasse | An energy-preserving discrete element method for elastodynamics[END_REF] for elastic computations on Voronoi meshes. In a consecutive work [START_REF] Marazzato | A variational discrete element method for quasistatic and dynamic elastoplasticity[END_REF], a variational DEM has been proposed for elasto-plasticity computations on polyhedral meshes using cell-wise reconstructions of the strains. The numerical results reported in [START_REF] Marazzato | A variational discrete element method for quasistatic and dynamic elastoplasticity[END_REF] confirmed that the macroscopic behaviour of elastic continua is indeed correctly reproduced by the variational DEM. The method developed in [START_REF] Marazzato | A variational discrete element method for quasistatic and dynamic elastoplasticity[END_REF] takes its roots in [START_REF] Eymard | Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilization and hybrid interfaces[END_REF] which is indeed a hybrid finite volume method. It is called variational DEM since it is possible to reinterpret the method as a consistent discretization of elasto-plasticity with discrete elements. In particular, a force-displacement interpretation of the method is derived from the usual stress-strain approach. Also, the mass matrix is diagonal and the stencil for the gradient reconstruction is compact as in usual DEM.

DEM for cracking have been developed in [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF] and [START_REF] André | A novel DEM approach for modeling brittle elastic media based on distinct lattice spring model[END_REF] with cracks propagating through the facets of the (Voronoi) mesh and using a critical stress criterion (initiation criterion). Coupled FEM-DEM techniques for crack computations, as [START_REF] Zárate | A simple FEM-DEM technique for fracture prediction in materials and structures[END_REF] (2d) and [START_REF] Zárate | A three-dimensional FEM-DEM technique for predicting the evolution of fracture in geomaterials and concrete[END_REF] (3d), have been introduced to take advantage of the FEM ability in computing elasticity and of the ability of DEM to handle cracked media. A similar approach, but using a different reconstruction of strains based on moving least-squares interpolations, can be traced back to [START_REF] Belytschko | Element-free Galerkin methods[END_REF] (2d) and [START_REF] Sukumar | An element-free Galerkin method for three-dimensional fracture mechanics[END_REF] (3d). Crack propagation can be based instead on the Griffith criterion which relies on the computation of the stress intensity factors (SIF) at the crack tip when coupled with the Irwin formula. Virtual element methods (VEM) have been recently applied to crack propagation [START_REF] Hussein | A virtual element method for crack propagation[END_REF]. Cracks were allowed to cut through the polyhedral mesh cells as in the extended finite element method (XFEM) which is based on an extended space of basis functions [START_REF] Chahine | Crack tip enrichment in the XFEM using a cutoff function[END_REF] and a level-set description of the crack [START_REF] Moës | X-FEM, de nouvelles frontières pour les éléments finis[END_REF]. Phase-field methods instead smooth the crack and have been developed among others in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] and subsequent work. Phase-field methods are not based on SIF computations but rather on a variational formulation of cracking [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. Furthermore, DEM using cohesive laws have been developed for fragmentation computations [START_REF] Mariotti | Modeling of the fragmentation by discrete element method[END_REF] with a view towards uniting initiation and propagation. These methods allow one to devise an initiation criterion and also to control the energy dissipation as with a Griffith criterion. The cracks still go through the mesh facets. This is also the case for similar methods of higher-order such as discontinuous Galerkin methods [START_REF] Hansbo | A discontinuous Galerkin method for cohesive zone modelling[END_REF].

The main goal of the present work is to develop a variational DEM using a Griffith criterion to compute crack propagation through the mesh facets. The method supports in principle polyhedral meshes, but the present numerical experiments are restricted to triangular meshes. The proposed method is close to [START_REF] Marazzato | A variational discrete element method for quasistatic and dynamic elastoplasticity[END_REF] (where there is no cracking) but the degrees of freedom (dofs) are different. Only cell dofs are used in the present work. The cracking algorithm hinges on two main ingredients. The first ingredient is an approximation of the energy release rate at every vertex along the crack. The second ingredient is a kinking criterion used to determine the next breaking facet and thus the crack path. The kinking criterion, in the spirit of [START_REF] Sih | Strain-energy-density factor applied to mixed mode crack problems[END_REF], consists in selecting for the crack path the inner facet of the mesh that maximizes a quantity representing the local density of elastic energy.

The present work is organized as follows. Section 2 briefly recalls the equations of elasticity and cracking in a Cauchy continuum. Section 3 introduces the proposed variational DEM and presents the space discretization of the governing equations. Moreover, a numerical test is reported to assess the convergence of the space discretization in the presence of a singularity. Section 4 addresses the full discretization of the quasi-static cracking problem. Section 5 contains numerical results on quasi-static crack propagation problems in two space dimensions. Finally, Section 6 draws some conclusions.

Governing equations for quasi-static cracking

We consider an elastic fragile material occupying the domain Ω ⊂ R 2 in the reference configuration and evolving over the finite pseudo-time interval [0, T ], T > 0, under the action of a volumetric force f and boundary conditions. The pseudo-time interval [0, T ] is discretized by means of (K + 1) discrete pseudo-time nodes (t k ) k∈{0,...,K} with t 0 := 0 and t K := T . The strain regime is restricted to small strains so that we use the linearized strain tensor ε(u) := 1 2 (∇u + (∇u) T ) ∈ R 2×2 , where u is the R 2 -valued displacement field. The material is supposed to be homogeneous and isotropic. The stress tensor σ(u) ∈ R 2×2 is such that

σ(u) := C : ε(u), ( 1
)
where C is the fourth-order stiffness tensor. The elastic material is characterized by the Young modulus E and the Poisson ratio ν or equivalently by the Lamé coefficients λ and µ. The boundary of Ω is partitioned as ∂Ω = ∂Ω D ∪ ∂Ω N , a Dirichlet condition is prescribed on ∂Ω D , and a Neumann condition on ∂Ω N , so that we enforce for all k = 0,

• • • , K, u = u D (t k ) on ∂Ω D , σ(u) • n = g N (t k ) on ∂Ω N . ( 2
)
Since cracking can occur, we denote Γ(t k ) the crack at the pseudo-time node t k and the actual domain at the pseudo-time node t k is

Ω(t k ) := Ω \ Γ(t k ). (3)
This implies that

∂Ω(t k ) = ∂Ω D ∪ ∂Ω N ∪ Γ(t k ). We enforce a homogeneous Neumann condition on Γ(t k ) for all k = 0, • • • , K, i.e., σ(u) • n = 0 on Γ(t k ). ( 4
)
Since we are interested in crack propagation, we assume that Ω(0) already contains a crack, i.e., Γ(0) = ∅. The crack Γ(t k ) is supposed to be a countably rectifiable 1-manifold for all k = 0, • • • , K (see [START_REF] Maso | Generalised functions of bounded deformation[END_REF]). This hypothesis ensures the almost everywhere (a.e.) existence of a normal vector n and a tangent vector τ to Γ(t k ) at any point y ∈ Γ(t k ) [START_REF] Simon | Lectures on geometric measure theory[END_REF]. Figure 1 illustrates these quantities.

• y n τ y • n τ e 3 Crack Γ(t k ) Ω(t k )
Figure 1: Sketch of a crack in the two-dimensional domain Ω(t k ).

The stress intensity factors (SIF) at any point y ∈ Γ(t k ) are usually defined for a purely elastic material as

       K 1 (y) := lim y →y σ nn (y ) 2πd(y, y ), K 2 (y) := lim y →y σ nτ (y ) 2πd(y, y ), (5)
where d(•, •) is the Euclidean distance in R 2 . If the stresses remain bounded in the vicinity of y ∈ Γ(t k ), then the SIF are null. Using the Irwin formula, one can define the energy release rate G(y) in the plane strain hypothesis as

G(y) := 1 -ν 2 E K 1 (y) 2 + K 2 (y) 2 . ( 6
)
Admissible states are characterized by the inequality

G(y) ≤ G c , ∀y ∈ Γ(t k ), (7)
where G c is a material property associated with the capacity of the material to sustain loads without locally failing and thus opening cracks. The material remains healthy at the point y ∈ Γ(t k ) if G(y) < G c and breaks if G(y) = G c . The material parameter G c is assumed to be homogeneous for simplicity.

To formulate the governing equations for quasi-static cracking, we consider the following functional spaces depending on the pseudo-time node t k :

V D (t k ) := v ∈ H 1 (Ω(t k ); R d ) | v |∂Ω D = u D (t k ) , V 0 (t k ) := v ∈ H 1 (Ω(t k ); R d ) | v |∂Ω D = 0 , ( 8
) where standard notation is used for the Hilbert Sobolev spaces. The weak solution is searched as a pair (u, Γ) such that for all k = 0,

• • • , K, u(t k ) ∈ V D (t k ), Γ(t k )
⊂ Ω is a 1-manifold satisfying the above assumptions, and

a(t k ; u(t k ), ṽ) = l(t k ; ṽ), ∀ṽ ∈ V 0 (t k ), G(y) ≤ G c , ∀y ∈ Γ(t k ). ( 9
)
Here we introduced the stiffness bilinear form such that for all (v, ṽ)

∈ V D (t k ) × V 0 (t k ), a(t k ; v, ṽ) := Ω(t k ) ε(v) : C : ε(ṽ), ( 10
)
and the linear form acting on V 0 (t k ) as follows:

l(t k ; ṽ) := Ω(t k ) f (t k ) • ṽ + ∂Ω N g N (t k ) • ṽ. ( 11
)
Note that the Dirichlet condition on ∂Ω D is enforced strongly, whereas the Neumann condition on ∂Ω N ∪ Γ(t k ) is enforced weakly.

Space semi-discretization

In this section, we present the space semi-discretization of (9) using a variational DEM.

Discrete sets and degrees of freedom

The domain Ω is discretized with a mesh T h of size h made of polygons with straight edges. We assume that Ω is itself a polygon so that the mesh covers Ω exactly. We also assume that the mesh is compatible with the initial crack position Γ(0) and with the partition of the boundary into the Dirichlet and Neumann parts. Recall that the space dimension is d = 2. Let C denote the set composed of the mesh cells and, for all k = 0, . . . , K, let F(t k ) denote the set composed of the mesh facets. This set depends on the pseudo-time node t k since a facet F ∈ F(t k ) is replaced, after cracking, by two boundary facets F -, F + ∈ F(t k ) (F -, F + are the same geometric object, but are different objects regarding the data structure since each one belongs to the boundary of a different mesh cell). The barycentre of a mesh cell c ∈ C is denoted by x c and the barycentre of a mesh facet F ∈ F(t k ) is denoted by x F .

Let t k be a pseudo-time node with k = 0, • • • , K. We partition the set of mesh facets as

F(t k ) = F i (t k ) ∪ F b (t k )
, where F i (t k ) is composed of the internal facets shared by two mesh cells and F b (t k ) is the collection of the boundary facets sitting on the boundary ∂Ω(t k ) = ∂Ω D ∪ ∂Ω N ∪ Γ h (t k ), where Γ h (t k ) denotes the discrete crack at t k . Notice that every boundary facet belongs to the boundary of only one mesh cell. The subsets F i (t k ) and F b (t k ) depend on the pseudo-time node t k since, as the facet F ∈ F i (t k ) cracks, it is replaced by the facets C) are placed at the barycentre of every mesh cell c ∈ C. We use the compact notation v h := (v c ) c∈C for the collection of all the cell dofs and we write 2 illustrates the position of the displacement dofs.

F + , F -∈ F b (t k ). The discrete crack Γ h (t k ) is composed of facets belonging to a subset of F b (t k ). This subset is denoted by F Γ (t k ) ⊂ F b (t k ).
(t k ) = F b N (t k ) ∪ F b D . Vector-valued volumetric degrees of freedom (dofs) for a generic displacement field (v c ) c∈C ∈ R d#(
v h ∈ V h := R d#(C) . Figure

Discrete bilinear and linear forms

The discrete stiffness bilinear form hinges on a reconstruction operator that provides a displacement value at every mesh facet by an interpolation formula from neighbouring cell dofs. Specifically, using the cell dofs of v h ∈ V h and the Dirichlet boundary conditions, we reconstruct a collection of displacements v

F := (v F ) F ∈F (t k ) ∈ R d#(F (t k )
) on all the mesh facets. The reconstruction operator is denoted R(t k ; •) and we write

v F := R(t k ; v h ) ∈ R d#(F (t k )) . ( 12
)
The reconstruction operator depends on t k because of the connectivity modifications due to the crack propagation. Let us first describe the reconstruction operator on boundary facets. Let F ∈ F b D be a Dirichlet boundary facet. Then the reconstruction is simply defined by evaluating the Dirichlet

Ω(t k ) ∂Ω(t k ) (u c ) c∈C Γ h (t k )
Figure 2: Domain Ω(t k ) covered by a polygonal mesh and vector-valued degrees of freedom for the displacement.

boundary condition at x F . Let F ∈ F b N (t k ) be a Neumann boundary facet. The main idea to define v F is to use a barycentric combination of the cell dofs close to F . A similar idea has been considered for finite volume methods in [13, Sec. 2.2] and for cell-centered Galerkin methods in [START_REF] Di | Cell centered Galerkin methods for diffusive problems[END_REF]. We thus select a subset of neighboring cell dofs of F , say I F ⊂ C, and set

v F := i∈I F α i (x F )v i , ( 13
)
where the v i 's are the dofs of v h and the coefficients α i (x F ) are the barycentric coordinates of the facet barycenter x F in terms of the selected positions of the dofs. For this construction to be meaningful, all the points associated with the selected dofs must not lie on the same line, so that, in particular, the cardinality of I F is at least (d + 1) = 3. Let us then describe the reconstruction for an inner facet F ∈ F i (t k ). We use a reconstruction similar to the one presented above except that the two cells sharing the inner facet F play symmetric roles. We refer to this construction as symmetric reconstruction. Specifically, let c + and c -be the two cells sharing the inner facet F ∈ F i (t k ). Then, we select I -(resp. I + ) as being composed of the cell c + (resp. c -) and of all the other cells sharing an inner facet with c -(resp. c + ). Notice that these two sets are disjoint. We then set

v F := 1 2 i∈I -∪I + α i (x F )v i , ( 14
)
so that, in the case of a simplicial mesh, 2(d + 1) dofs are used for the reconstruction (always including c -and c + ). Note that i∈I -α i (x F ) = i∈I + α i (x F ) = 1 here. Figure 3 presents an example where c -= c i , c + = c j , I -= {j, j 2 , j 3 } and

I + = {i, i 2 , i 3 }.
Having defined the reconstructed facet displacements, it is now possible to devise a discrete R d×d -valued piecewise-constant gradient field for the displacement that we write

G C (v F ) := (G c (v F )) c∈C ∈ R d 2 #(C) . Specifically, we set in every mesh cell c ∈ C, G c (v F ) := F ∈∂c |F | |c| v F ⊗ n F,c , (15)
where the summation is over the facets F of c and n F,c is the outward normal to c on F . Note that ( 15) is motivated by a Stokes formula and that for all v h ∈ V h , we have since F ∈∂c |F |n F,c = 0. We define a constant linearized strain tensor in every mesh cell c ∈ C such that

G c (R(t k ; v h )) = F ∈∂c |F | |c| (R(t k ; v h ) F -v c ) ⊗ n F,c , ( 16
)
x i x j x i 2 x i 3 x j 2 x j 3 F x F
ε c (v F ) := 1 2 (G c (v F ) + G c (v F ) T ) ∈ R d×d , ( 17
)
and a constant stress tensor in every mesh cell c ∈ C such that

Σ c (v F ) := C : ε c (v F ) ∈ R d×d . ( 18
)
Finally, we define an additional reconstruction that is used to formulate the stabilization bilinear form in the discrete problem (see below). This operator is a cellwise nonconforming P 1 reconstruction R c defined for all c ∈ C by

R c (t k ; v h )(x) := v c + G c (R(t k ; v h )) • (x -x c ), ∀x ∈ c. ( 19
)

Discrete problem

We set

V hD (t k ) := {v h ∈ V h | R(t k ; v h ) F = u D (t k ; x F ), ∀F ⊂ ∂Ω D }, ∀k = 0, • • • , K, V h0 (t k ) := {v h ∈ V h | R(t k ; v h ) F = 0, ∀F ⊂ ∂Ω D }, ∀k = 0, • • • , K. ( 20
)
The discrete stiffness bilinear form is such that for all (v h , ṽh )

∈ V hD (t k ) × V h0 (t k ) (compare with (10)) a h (t k ; v h , ṽh ) := c∈C |c|ε c (R(t k ; v h )) : C : ε c (R(t k ; ṽh )) + s h (t k ; v h , ṽh ), (21)
where the stabilization bilinear form s h is intended to render a h coercive and is defined as

s h (t k ; v h , ṽh ) = F ∈F i (t k ) 2µ h F |F |[R(t k ; v h )] F • [R(t k ; ṽh )] F + F ∈F b D 2µ h F |F |[R(t k ; v h )] F • [R(t k ; ṽh )] F , ( 22
)
where h F is the diameter of the facet F ∈ F(t k ). For an interior facet F ∈ F i (t k ), writing c - and c + the two mesh cells sharing F , i.e., F = ∂c -∩ ∂c + , and orienting F by the unit normal vector n F pointing from c -to c + , the jump of

R(t k ; v h ) across F is defined as [R(t k ; v h )] F := R c -(t k ; v h )(x F ) -R c + (t k ; v h )(x F ). ( 23
)
The sign of the jump is irrelevant in what follows. The role of the summation over the interior facets in ( 22) is to penalize the jumps of the cell reconstruction R across the interior facets. For a Dirichlet boundary facet F ∈ F b D , we denote c -the unique mesh cell containing F , we orient F by the unit normal vector n F := n c -which points outward Ω, and we define

[R(t k ; v h )] F := R(t k ; v h ) F -R c -(t k ; v h )(x F ). ( 24
)
Let us recall that for

u h ∈ V hD (t k ), R(t k ; u h ) F = u D (t k ; x F ) and for v h ∈ V h0 (t k ), R(t k ; v h ) F = 0.
The role of the summation over the Dirichlet boundary facets in ( 22) is to penalize the jumps between the cell reconstruction R and the value interpolated in the Dirichlet boundary facets. The bilinear form s h is classical in the context of discontinuous Galerkin methods (see [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF][START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF] for instance, see also [START_REF] Di | Cell centered Galerkin methods for diffusive problems[END_REF] for cell-centred Galerkin methods). It is possible to replace the coefficient 2µ in ( 22) by βµ with a user-dependent dimensionless parameter β of order unity. The numerical experiments reported in [START_REF] Marazzato | A variational discrete element method for quasistatic and dynamic elastoplasticity[END_REF] indicate that this choice has a marginal influence on the results.

Verification test case

This section presents a verification test case related to the convergence rate with a singularity at the crack tip. The crack does not propagate, i.e., we consider a steady setting using the above discrete stiffness bilinear form and load linear form. The convergence rate of the method in the presence of a singularity is tested in the case of an infinite plate under mode 3 loading at infinity as presented in Figure 4. A convergence rate of O(h

1 2
), similar to that obtained with Lagrange

a a σ • n = 0 θ r ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ ⊗ -τ e z τ e z σ • n = 0 σ • n = 0 e x e y e z
u(r, θ) = 2τ µ ar 2 sin θ 2 e z , ( 25
)
where τ is the modulus of the antiplane shear stress imposed at infinity. The displacement defined in [START_REF] Moës | X-FEM, de nouvelles frontières pour les éléments finis[END_REF] verifies the statics equation in a strong form since div(u) = 0. The stresses are

σ(r, θ) = τ a 2r sin θ 2 e r -cos θ 2 e θ ⊗ e z . ( 26
)
The domain shown in Figure 4 being symmetric with respect to the red dashed line, only its right part is considered. As the analytical solution ( 25) is only valid close to the crack tip, a small ball around the crack tip, which corresponds to the green dashed circle in Figure 4, is meshed. The setting is presented in Figure 5. The convergence towards the analytical solution is checked on the meshed ball with the reference solution imposed as Dirichlet boundary condition over the whole boundary including the crack lips.

, ( 27
)
where e 1 , e 2 denote the errors on the computations with mesh sizes h 1 , h 2 and the number of dofs n 1 , n 2 .

Quasi-static crack propagation

In this section, we formulate the discrete problem for quasi-static crack propagation. The space discretization is achieved by means of the variational DEM scheme presented in the previous section. At every pseudo-time node t k , the problem is solved iteratively with inner iterations enumerated by m ∈ {0, . . . , M }. Since the crack can change at each inner iteration, we use the notation Γ h (t k,m ) for the crack and the notation F i (t k,m ) and F b (t k,m ) for the partition of the mesh facets at the inner iteration m, with the facets located in the crack collected in the subset F Γ (t k,m ). Each inner iteration consists in two steps. First, freezing the position of the crack, we find the discrete displacement

u h (t k,m ) ∈ V hD (t k ) solving the quasi-static problem a h (t k,m ; u h (t k,m ), ṽh ) =
l h (t k ; ṽh ) for all ṽh ∈ V h0 (t k ) (the bilinear form a h depends on t k,m since the reconstruction operator changes as the crack propagates). Then we use the newly computed displacement field u h (t k,m ) to determine whether crack propagation occurs and update accordingly the subsets F i (t k,m+1 ), F b (t k,m+1 ), and F Γ (t k,m+1 ). We iterate this procedure until there is no more crack propagation in the second step. The inner iteration in the discrete quasi-static crack propagation scheme can thus be summarized as follows: For all m ∈ {0, . . . , M },

(i) u h (t k,m ) ∈ V hD (t k ) s.t. a h (t k,m ; u h (t k,m ), ṽh ) = l h (t k ; ṽh ), ∀ṽ h ∈ V h0 (t k ), (ii) (F Γ (t k,m+1 ), F b (t k,m+1 ), F i (t k,m+1 )) = CRACK QS(F Γ (t k,m ), F b (t k,m ), F i (t k,m ), u h (t k,m )).
(28) The rest of this section is devoted to the description of the procedure CRACK QS. This procedure consists in the three consecutive steps outlined in Figure 6. The first step involves the procedure ESTIMATE which considers all the vertices of F Γ (t k,m ) and computes for each of these vertices an approximate energy release rate. The second step involves the procedure MARK which flags among all the inner facets sharing a vertex with an energy release rate larger than the maximum value G c the facet that will indeed break. The selection is made by using a discrete kinking criterion. The last step uses the procedure UPDATE and simply consists in updating the data structure according to the crack propagation. The procedure is repeated from the recomputation of the solution of the first line of Equation ( 28) until no facet is marked in the procedure MARK.

Procedure ESTIMATE

Let V Γ (t k,m ) be the set of all vertices in the crack Γ(t k,m ). The procedure ESTIMATE computes an approximate energy release rate

G h (v) for all v ∈ V Γ (t k,m ). Let F Γ v (t k,m ) be the set of cracked facets sharing a vertex v ∈ V Γ (t k,m ). (The set F Γ v (t k,m ) reduces to a single facet if v is the crack tip.) Let F i v (t k,m ) be the set of inner facets sharing a vertex v ∈ V Γ (t k,m
). An approximate energy release rate for the vertex v ∈ V Γ (t k,m ) is evaluated as

G h (v) := max F ∈F Γ v (t k,m ) max F ∈F i v (t k,m ) πn F • {Σ h (t k,m )} F • [u h (t k,m )] F , ( 29
)
where

[u h ] F := u c --u c + , {Σ h } F := 1 2 (Σ c -+ Σ c + )
, and n F is the normal vector to F pointing from c -to c + . This expression is rooted in the fact that the elastic energy contained in a facet [START_REF] Marazzato | A variational discrete element method for quasistatic and dynamic elastoplasticity[END_REF]. The factor π comes from the fact that the density of elastic energy per facet must be multiplied by 2π to take into account the surface created by cracking (see [19, p. 48]). This is linked to the concept of the crack closure integral. The output of the procedure ESTIMATE is the collection of approximate energy release rates {G h (v)} v∈V Γ (t k,m ) .

F writes 1 2 n F • {Σ h (t k,m )} F • [u h (t k,m )] F |F | as motivated in

Procedure MARK

The goal of the procedure MARK is to identify the unique inner facet F ∈ F i (t k,m ) through which the crack will propagate. The criterion is based on an adaptation of the maximisation of the strain energy density which was introduced in [START_REF] Sih | Strain-energy-density factor applied to mixed mode crack problems[END_REF]. The vertices of V Γ (t k,m ) are ordered as they break during a computation and we select the last N vertices in V Γ (t k,m ) to define the subset V Γ N (t k,m ). The integer parameter N is set to N = 6 in our computations; this choice gives satisfactory results while avoiding excessive branching of the crack path. Finally, we select the vertices in V Γ N (t k,m ) whose approximate energy release rate is larger than the material parameter

G c : V Γ * N (t k,m ) := {v ∈ V Γ N (t k,m ), G h (v) ≥ G c }. (30)
Among all v ∈ V Γ * N (t k,m ), we select the single vertex through which the crack will propagate at t k,m as z := Argmax

v∈V Γ * N (t k,m ) G h (v). ( 31
)
If there is more than one maximizer, one is picked randomly. Note that in most situations, the vertex z is located at the crack tip.

Having selected the vertex z, we now mark one facet F ∈ F i z (t k,m ) for cracking. We impose only one restriction on the selection process of the facet to be broken: we limit the number of facets broken per cell to one. This limit is justified by the fact that when a facet breaks, the resulting geometric singularity creates very high stresses that lead to breaking the other facets of the cells containing the facet thus creating many fragments. The limitation we impose is to avoid this situation. The setting is illustrated in Figure 7. The output of the procedure MARK is the facet F, through which the crack will propagate, defined as

z F Γ (t k,m ) F ∈ F i z (t k,m ) \ {F} V Γ * (t k,m ) V Γ (t k,m ) F
F := Argmax F ∈F i z (t k,m )\F i C (t k,m ) 1 2 {Σ h (t k,m )} F • {ε h (t k,m )} F , ( 32
)
where

F i C (t k,m
) denotes the set of inner facets contained in a cell with one facet already broken.

Procedure UPDATE

The subsets F Γ (t k,m+1 ), F i (t k,m+1 ), and F b (t k,m+1 ) can now be updated as follows:

       F Γ (t k,m+1 ) := F Γ (t k,m ) ∪ {F}, F i (t k,m+1 ) := F i (t k,m ) \ {F}, F b (t k,m+1 ) := F b (t k,m ) ∪ {F -, F + }, (33)
where we recall that F -and F + are the same geometric object as the inner facet F, but are now each one on the boundary of a single mesh cell.

Remark 1 (Update of a h ). The updates in [START_REF] Zárate | A three-dimensional FEM-DEM technique for predicting the evolution of fracture in geomaterials and concrete[END_REF] affect the reconstruction operator used to evaluate the discrete stiffness bilinear form. Figure 8 presents a sketch of an inner facet whose reconstruction has to be recomputed after a neighbouring inner facet breaks. The purpose of recomputing the reconstruction on certain inner facets is to avoid using dof values on both sides of the crack in the same reconstruction.

x j

x i F x k x l x m x n F x F x j x i x k F x F F - F + x l x m x n

Numerical experiments

Several numerical experiments are presented to show the versatility of the proposed numerical method. The python scripts1 for these numerical experiments use the finite element library FEniCS [START_REF] Logg | Automated Solution of Differential Equations by the Finite Element Method[END_REF] and scipy2 . Although the proposed method is able to handle polyhedral meshes, our computations only use triangular meshes. This is a consequence of the current restriction of FEniCS to simplicial meshes.

Crack speed with prescribed crack path

We consider a test case taken from [START_REF] Li | Numerical investigation of dynamic brittle fracture via gradient damage models[END_REF]. The test case consists of an already cracked plate under antiplane shear loading. The crack is forced to propagate along a straight line represented by the dashed line in Figure 9. The goal of this test case is to study the crack propagation velocity. The dimensions of the plate are L = 5m and H = 1m and the initial length of the crack is l 0 = 1m. The constant increment in boundary loading is written ∆u D . The material parameters are µ = 0.2Pa and G c = 0.01kN/mm. We are interested in the length of the crack with respect to the cumulated boundary loading displacement u D , where the final displacement load is u D = 1m. The reference solution for the crack speed S with respect to the loading speed, taken from [START_REF] Li | Numerical investigation of dynamic brittle fracture via gradient damage models[END_REF], is µH Gc ≈ 4.47. As this solution is only valid when L → ∞, we checked that see that the results with the two meshes are very similar. The results with ∆u D = 10 -3 m and ∆u D = 10 -2 m are very similar and are in agreement with the analytical solution. For these two values, ∆u D h is less than 0.5, so that the increment in the imposed Dirichlet condition is smaller than the mesh size. This is not the case for ∆u D = 0.1m. The different aspect of the curves for ∆u D = 0.1m is explained by the fact that as ∆u D is large in that case, a large number of facets can break at some of the displacement increments, thus leading to this staircase shape. However, one can notice that at the end of every other displacement increment, the curve for ∆u D = 0.1m reaches the same value as the curves computed with the other ∆u D values. Table 2 contains the errors of the crack speeds (computed with a least-squares fit on the two numerical computations) with respect to the analytical solution. The agreement of the computed crack ∆u D / h 0.1 0.05 0.1 3.6% 5.1% 0.01 2.0% 0.68% 0.001 1.9% 0.70% speeds with the analytical solution is very satisfactory for all ∆u D .

L σ • n = 0 H l 0 σ • n = 0 σ • n = 0 u = 0 u = -u D (t) u = u D (t)

Opening mode with unknown crack path

The setting for this test case is presented in Figure 11. The dimensions of the plate are L = 32mm and H = 16mm and the initial length of the crack is l 0 = 4mm. The material parameters are E = 3.09GPa, ν = 0.35 and G c = 300kN/mm. First, we use a structured mesh of size h = 0.4mm leading to 25, 920 dofs. The increment in boundary conditions is defined as ∆u D = h. Figure 12 presents the obtained crack path. We notice an unstable crack propagation, as expected, in the sense that when the propagation starts, it breaks the entire sample at a given t k .

L σ • n = 0 H l 0 u • n = u D (t) u • n = u D (t) σ • n = 0 σ • n = 0
Figure 12: Opening mode: u y in colors and crack path in white on a structured mesh, u D = 2.5mm.

We also perform computations on two unstructured meshes of sizes h = 1.4mm and h = 0.74mm corresponding respectively to 2, 792 dofs and 11, 044 dofs. Both meshes do not contain facets with a direction that could lead to a totally straight propagation of the crack. The finer mesh is not a refinement of the coarser one. Figure 13 shows the crack paths obtained on the two meshes. The crack paths obtained are satisfactory as the propagation is rather straight and the results on the two meshes are quite similar.

Single-edge notched shear test

The setting of this test case comes from [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF]. It consists in a square with an already initiated crack loaded in shear on its top surface. The lower surface is recessed while the upper surface is loaded in shear. The two lateral parts are free of stress as well as the crack. Figure 14 illustrates the setting. The crack is of initial length l 0 = 0.5mm and the dimension of the sample is H = 1mm. Three computations are performed on unstructured meshes of size h = 2.8 • 10 -2 mm (coarse mesh), h = 1.3 • 10 -2 mm (fine mesh), and h = 7.7 • 10 -3 mm (finest mesh), leading respectively to 13, 396, 65, 956, and 210, 328 dofs. Figure 15 shows the computed crack paths. Our results can be compared with [START_REF] Muixí | A hybridizable discontinuous galerkin phase-field model for brittle fracture with adaptive refinement[END_REF] which uses a phase-field model discretized by a hybridizable discontinuous Galerkin formulation. The computations are in satisfactory agreement with those of [27] regarding the general orientation of the crack and the number of branches. We observe in Figure 15 that the crack propagates downwards along a somewhat curved path (with rather close predictions between the two finer meshes). The trajectory is sightly different from the one predicted in [START_REF] Muixí | A hybridizable discontinuous galerkin phase-field model for brittle fracture with adaptive refinement[END_REF] where the crack propagates along a rather straight line which forms a sharp angle with respect to the initial crack. Experimental results would be needed to assess the correctness of these numerical results.

H H l 0 u = (u D (t k ), 0) u = 0
The load-displacement curves are displayed in Figure 16 along with the values of the imposed displacement and the resulting force when the crack starts propagating. The force is computed through an integration of the tangential component of the reconstructed normal stress Σ h • n on the upper and lower surfaces of the sample. The force has also been computed through a residual method and the difference has been found to be negligible. One can first notice that up to an imposed displacement of 9.5µm, all the curves are superimposed and exactly reproduce the elastic response of the sample with the fixed initial crack. As the imposed displacement increases beyond the above value, jumps in the load-displacement curves appear progressively. These jumps are a consequence of facets cracking, and the slope of the elastic response is reduced after each jump owing to the propagation of the crack. This explains the observed zigzag behavior of the response curves. Altogether, crack propagation thus induces a softening of the sample as expected. The load-displacement curve obtained on the coarse mesh stops at the value u D = 0.2mm for which the crack reaches the rightmost boundary of the sample. Instead, the computations on the two finer meshes support larger values for u D and lead to rather similar predictions. Furthermore, one can see that the crack starts propagating around an imposed displacement of 10µm, which is similar to the value reported in [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF]. The value of the force, however, is different. We believe that this difference can be attributed to the sharp interface representation of the crack in the present method. To substantiate this claim, we performed some additional computations on the finest mesh using a fixed interface position, P 1 -Lagrange finite elements, and an imposed displacement u D = 5µm. With a sharp interface, the load is 0.13kN (consistently with the DEM prediction on the same mesh), whereas it is 0.32kN if there is no crack (the sample is fully sound). If instead the initial crack is represented as a damage field [?, 27] with a smoothing length = 5h, the load is close to the value reported in [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF][START_REF] Muixí | A hybridizable discontinuous galerkin phase-field model for brittle fracture with adaptive refinement[END_REF], namely 0.20kN (notice that this value is as expected in the interval (0.13, 0.32)kN).

Notched plate with a hole

This test case comes from [START_REF] Muixí | A hybridizable discontinuous galerkin phase-field model for brittle fracture with adaptive refinement[END_REF]. The material parameters are E = 6GPa, ν = 0.22, and G c = 2.28•10 -3 kN/mm. We use fixed displacement increments of ∆u D = 10 -2 mm. Figure 17 presents a sketch of the sample. The dimensions of the plate are L = 65mm and H = 120mm. The two holes on the left of the sample have a diameter of 10mm and the hole on the right of the sample has a dimeter of 20mm. The initial length of the crack is l 0 = 10mm. One also has a = 20mm, b = 55mm, d = 69mm and e = 36.5mm. The right hole is free of stress, the lower hole is recessed and the upper hole has an imposed displacement u = (0, u D (t k )). We use three unstructured meshes with h = 2.8mm, h = 1.5mm and h = 0.78mm having respectively 9, 926, 39, 380 and 157, 340 dofs.

Figure 18 shows the computed crack paths. We compare our results with [START_REF] Muixí | A hybridizable discontinuous galerkin phase-field model for brittle fracture with adaptive refinement[END_REF] without taking into account the secondary crack starting from the largest hole as we restrict ourselves to crack propagation and not crack initiation. We notice that for the three computations, the crack goes towards the largest hole in a similar fashion which also seems consistent with [START_REF] Muixí | A hybridizable discontinuous galerkin phase-field model for brittle fracture with adaptive refinement[END_REF]. The loaddisplacement curves are given in Figures 19 and20, together with the values of the imposed displacement and the resulting force when the crack starts propagating and when the crack mesh are still rather inaccurate for higher values of u D , whereas Figure 20 indicates that the predictions on the two finer meshes are in satisfactory agreement as far as the load-displacement curves are concerned. The predictions of the path of crack propagation are also similar on both meshes, but the value of the imposed displacement when the crack reaches the hole is different, as reflected in the caption of Figure 18.

Conclusions

We have presented a variational Discrete Element Method (DEM) to compute Griffith crack propagation. The crack propagates through the facets of the mesh and thus between discrete elements. The variational DEM is a consistent discretization of a Cauchy continuum and only requires three continuum macroscopic parameters for its implementation: the Young modulus, the Poisson ratio, and the critical energy release rate. The displacement degrees of freedom are attached to the barycentre of the mesh cells. A discrete Stokes formula is used to devise a piecewise constant gradient and linearized strain reconstructions. An approximation of the energy release rate is computed in the procedure ESTIMATE. The procedure MARK then determines the breaking facet at each pseudo-time node t k . Finally, the procedure UPDATE updates the necessary discrete quantities after the facet that has been marked has been broken. A convergence test in antiplane shear has confirmed the efficiency of the variational DEM discretization as well as the O(h 1 2 ) convergence rate in energy norm. The robustness of the method regarding the computation of the crack speed has been verified. Also, several numerical experiments have shown that the method can provide reasonable crack paths.

This work can be pursued in several directions. A first idea would be to adapt the present methodology to three-dimensional problems with two-dimensional cracks. A second direction concerns the regularity of the crack surface. Indeed, in the spirit of [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], a crack should be a surface that minimizes energy. To achieve this goal, the variational DEM could be coupled to gradient flows used for surface lifting, as in [START_REF] Romon | Introduction à la géométrie différentielle discrète[END_REF], with the goal of moving the crack surface vertices. One would then have to verify the convergence of the discrete crack area with tools similar to [START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF]. A third direction for further study is to approximate cohesive cracking laws instead of a Griffith cracking law so as to enable the simulation of crack initiation as well as crack propagation. Inspiration can be found in [START_REF] Mariotti | Modeling of the fragmentation by discrete element method[END_REF] which uses a DEM with a linear cohesive law. Finally, a last direction can be to consider an enrichment similar to [START_REF] Chahine | Crack tip enrichment in the XFEM using a cutoff function[END_REF] close to the crack tip so as to obtain a convergence with order O(h).

  We also introduce the partition between boundary facets with Neumann boundary conditions F b N (t k ) (recall that homogeneous Neumann boundary conditions are imposed on newly created crack lips) and with Dirichlet boundary conditions F b D which does not depend on t k . One thus has F b

Figure 3 :

 3 Figure 3: Dofs associated with the interior facet F used in the reconstruction.

Figure 4 :

 4 Figure 4: Sketch of the antiplane shear experiment in an infinite plate.

Figure 6 :

 6 Figure 6: Details of the procedure CRACK QS.

Figure 7 :

 7 Figure 7: Sketch of the discrete sets considered in the procedure MARK.

Figure 8 :

 8 Figure 8: Recomputation of the reconstruction stencil associated with the inner facet F after the breaking of the neighbouring inner facet F . Left: reconstruction before cracking. Right: reconstruction after cracking. (The two cells separated by the crack are drawn slightly apart.)

Figure 9 :

 9 Figure 9: Crack speed: problem setup.

Figure 10 :

 10 Figure 10: Crack speed: crack length as a function of loading displacement u D . Left: coarse mesh. Right: fine mesh.

Figure 11 :

 11 Figure 11: Opening mode: setup.

Figure 13 :

 13 Figure 13: Opening mode: u y in colors and crack path in white on a unstructured mesh. Left: coarse mesh, u D = 2.5mm. Right: fine mesh, u D = 3.4mm.

Figure 14 :

 14 Figure 14: Single-edge notched shear test: setup.

Figure 15 :

 15 Figure 15: Single-edge notched shear test: u x in colors and crack path in white, u D = 0.2mm. Left: coarse mesh. Middle: fine mesh. Right: finest mesh.

Figure 17 :

 17 Figure 17: Notched plate with a hole: setup.

Figure 18 :Figure 19 :Figure 20 :

 181920 Figure 18: Notched plate with a hole: u y in colors and crack path in white (zoom around the main hole). Left: coarse mesh, u D = 0.65mm. Middle: fine mesh, u D = 1.4mm. Right: finest mesh, u D = 3.0mm.

Table 1 :

 1 The results of the computation, which are reported Number of dofs, L 2 -error and convergence rate, L 2 -error on the gradient and convergence rate.

		r
	e	θ
	e x

y Figure 5: Sketch of the meshed ball around the crack tip.

in Table

1

, corroborate an O(h 1 2 ) convergence rate in the energy-norm, as expected. We also

nb dofs u -R(u h ) L 2 Convergence rate ∇u -G h (u h ) L 2 Convergence rate

Table 2 :

 2 Crack speed: error with respect to analytical solution depending on the choice of h and ∆u D .

https://github.com/marazzaf/DEM_cracking.git

https://scipy.org/

Acknowledgements

Partial support by CEA is gratefully acknowledged.


Maria Saladina

Pablo Simón Marqués

Anastasia Markina

Safakath Karuthedath

Christopher Wöpke

Clemens Göhler

Yue Chen

Magali Allain

Philippe Blanchard

Clément Cabanetos

Denis Andrienko

Frédéric Laquai

Julien Gorenflot

Carsten Deibel

Charge Photogeneration in Non-Fullerene Organic Solar Cells: Influence of Excess Energy and Electrostatic Interactions

Keywords: organic solar cells, non-fullerene acceptors, dissociation, excess energy, molecular quadrupole moment

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The field of organic solar cells has shifted its focus to non-fullerene acceptors (NFA) since 2015. 1,2 Due to the fast development of new molecules and the variety of their chemical structures, understanding of the underlying photophysics in non-fullerene donor-acceptor (D-A) blends is far from complete. Knowledge of the factors that influence photogeneration in these systems is essential to move the field forward. Charge generation in organic solar cells is a two-step process, owing to their low dielectric constant. [3][4][5] Following the charge transfer from donor to acceptor (or vice versa), a nonrelaxed chargetransfer (CT) state is formed. It undergoes thermalisation in the CT manifold and results in formation of a bound electron-hole pair with initial separation distance r 0 . The electron-hole pair can subsequently relax to the lowest CT state, and recombine geminately to the ground state or dissociate into free charge carriers. 6,7 The probability of CT dissociation η diss depends on the energy barrier between the CT and the charge-separated (CS) state. [8][9][10] An activation energy for free charge formation is found to be consistently smaller than predicted for a Coulombically bound electron-hole pair, [11][12][13][14] although the reasons for this discrepancy are still under debate. This barrier lowering has been explained by energetic disorder effects, interfacial dipoles, entropy etc. [15][16][17][18][19][20] In fullerene-based organic solar cells, the energetics of CT and CS states are influenced by the quadrupole moment of the donor. 21 The electrostatic potential at the D-A interface was shown to reduce the dissociation energy required for charge separation. 22 The magnitude of the quadrupole moment depends largely on the molecular structure. In efficient non-fullerene organic solar cells, small molecule acceptors with the A-D-A structure usually possess a large quadrupole moment. [23][24][25] Surprisingly, its influence on charge separation has received little experimental attention up to now.

The binding energy of an electron-hole pair depends on its separation distance. 26,27 Greater delocalisation 28 and higher local mobility 29 of CT states were shown to be important for efficient photogeneration in polymerfullerene based solar cells. It is plausible that during thermalisation in the CT manifold, the electron and hole with excess kinetic energy are driven apart to a larger thermalisation distance r 0 . Thus excess energy may increase the probability for a CT state to overcome the potential barrier to CS states. 6,30 The influence of excess energy on photogeneration yield up to this date was investigated in donor-fullerene systems. Several authors have reported a positive effect of excess photon energy on the yield of free charge formation in these systems. Ohkita et al. demonstrated that the yield of free charge formation in polythiophene:fullerene blends increased with the free energy difference between the singlet and the CS state. 31 By em- ploying pump-push photocurrent measurements to the polymer BTT-DPP blended with the fullerene acceptor PC 61 BM ( [6,6]-phenyl-C 61 -butyric acid methyl ester), Dimitrov et al. could observe an increase in the photocurrent yield after applying an infrared push pulse providing bound CT states with excess energy. 32 Transient absorption spectroscopy (TAS) experiments have shown that the free charge carrier generation yield is dependent on the excitation energy above the optical gap, hence on the excess energy of singlet excitons. 32,33 The results were corroborated by the energy dependent internal quantum efficiency (IQE) above the optical gap. It was proposed that higher energy singlet states transformed into more delocalised hot interfacial charge transfer states that could split more easily.

In contrast, several studies of photogeneration pathways in fullerene systems employing steady-state measurements have reported negligible influence of excess energy of both singlet excitons and charge transfer states on charge separation. 12,34,35 By using a combination of time delayed collection field (TDCF), electroluminescence, photothermal deflection spectroscopy and sensitive EQE measurements, Vandewal et al. found no difference in the IQE or the field dependence of charge generation upon above and below the optical gap excitation in small-molecule and polymer-based organic solar cells. 36 In non-fullerene organic solar cells the influence of excess energy on the yield of photogeneration has been investigated in only one material system. Recently, Perdigón-Toro et al. studied a highly efficient PM6:Y6 non-fullerene system, in which free charge formation was probed by TDCF and found to proceed through relaxed CT states virtually independent of excitation energy. 25 However, at this point little is known about the generality of this result.

In this work, we aim to determine the influence of photogeneration pathways and electrostatic potential at the D-A interface on the yield of free charge carrier formation in non-fullerene organic solar cells. For this purpose, we employ TDCF measurements in two NFA systems at excitation energies above and below the optical gap, allowing us to assess the efficiency of electron-hole pair separation at the interface between donor and acceptor. To support our analysis, we examine photogeneration and recombination in these systems with sensitive photocurrent and TAS measurements.

Through careful assessment of the CT dissociation efficiency at various excitation energies, we demonstrate that non-fullerene organic solar cells benefit from photogeneration via nonrelaxed CT states exhibiting higher dissociation yield. Using the quantum chemical calculations and the kinetic model of photogeneration we relate the quadrupole moments of acceptor molecules to the CT binding energy. We find that the solar cells comprising of the small molecule acceptor with a larger molecular quadrupole moment show higher CT separation yield due to the reduced activation energy for free charge carrier formation.

II. RESULTS

We focus on two non-fullerene acceptor systems, PTB7-Th:ITIC and PTB7-Th:h-ITIC. h-ITIC was designed as the dipolar analogue of the reference acceptor ITIC (Figure 1(a)). To that end, the new acceptor presents a push-pull π-conjugated structure, with an indenothienothiophene as electron-donating group connected to only one electron-withdrawing moiety, namely dicyanovynilindanone (DCI). In contrast, ITIC exhibits an indenodithienothiophene central core, which is functionalised with two DCI groups, showing lateral symmetry. 37 Therefore, both compounds open an excellent scenario for comparing two analogues, one with a quadrupole moment (ITIC) and its counterpart with a dipole moment (h-ITIC). The polymer donor PTB7-Th was chosen to complement the absorption spectrum of h-ITIC molecule (Figure S10), while PTB7-Th and ITIC absorb in the same region.

We perform gas-phase quantum chemical calculations to obtain the dipole and quadrupole moments of the ITIC and h-ITIC small molecule acceptors. The mutual arrangement of molecules, dipole (D) and quadrupole moments (Q total 20 , Q mol 20 ) for single and dimer configurations are shown in Figure 1(c). ITIC molecules have zero charge and dipole moment, while their quadrupole moments are non-zero due to the A-D-A architecture. Both single and dimer states yield approximately the same quadrupole moment per molecule (Q mol 20 ≈ 49.2 ea 2 0 ). h-ITIC acceptor is an asymmetric molecule with an A-D architecture. Electro-negative and electro-positive areas are distributed in such a way that they lead to the formation of a non-zero dipole moment. However, this dipole moment is cancelled out at a certain molecular orientation in a dimer, leading to the rise in a quadrupole moment value.

Indeed, in Figure 1(c) we see that in the dimer configuration, the dipole moments are significantly reduced, while the quadrupole moment is non-zero (Q total 20 ≈ 38.8 ea 2 0 ). Overall dipole moment of the unit cell is nearly zero D = 0.04 ea 0 (cf. Figure S12). However, the quadrupole moment of h-ITIC dimer is lower than the quadrupole moment of a single ITIC molecule while the effective van der Waals volume is higher. Using a lattice model approximation of point quadrupoles placed on a grid one can estimate electrostatic contribution in both cases of h-ITIC and ITIC. Lower quadrupole moments and higher grid space lead to the lower electrostatic potential in the case of h-ITIC (cf. Figure S13). When blended with a donor polymer, this difference in electrostatic potential of the acceptors influences the device performance through charge separation at the D-A interface and charge transport. 38 Current-voltage characteristics of the devices in cells points to substantial geminate recombination and collection efficiency losses in this material system.

We employ TAS to determine the relative importance of geminate and nongeminate recombination in ITIC and h-ITIC blend films in the absence of an external electric field. In Figure 2 we plot the normalised decay dynamics of the charge carrier density n at the excitation energy of 2.33 eV for different pump fluences with pump delays from 1 ns to 70 µs. The data is fitted globally assuming that the charge carrier population consists of bound electron-hole pairs and free charge carriers. 39

n(t) = (1 -f ) n 0 exp - t τ + λγt + (f n 0 ) -λ -1/λ , (1)
where n 0 is the initial charge carrier density, τ is the CT lifetime, γ is the nongeminate recombination rate constant and λ + 1 the recombination order. 40 The fractions of CS and CT states are denoted f and 1 -f , Table 1. Parameters extracted from the fits with Equation (1) to the charge carrier density kinetics from TAS.

Parameter PTB7-Th:ITIC PTB7-Th

:h-ITIC f 0.82 0.71 τ [ns] 2.80 2.83 λ + 1 2.40 2.64 γ [(cm 3 ) λ s -1 ] 6.6 • 10 -17 4.2 • 10 -21 γ eff [cm 3 s -1 ] 1.7 • 10 -10 6.3 • 10 -11
respectively. The effective Langevin recombination rate constant was calculated according to γ eff = γn λ-1 with n = 1 • 10 16 cm -3 . n 0 was estimated from TDCF measurements at similar conditions (cf. Figure S15).

From the fits we obtain recombination parameters for the PTB7-Th:ITIC and PTB7-Th:h-ITIC blend films, listed in Table 1. The fraction of geminate recombination is 18% and 29%, respectively, with a typical CT lifetime of 2.8 ns in both cases, thus the dissociation of electron-hole pairs at the D-A interface is less efficient in PTB7-Th:h-ITIC. The faster nongeminate recombination in PTB7-Th:ITIC blend is overcompensated by better collection efficiency due to higher and more balanced charge carrier mobilities, 41,42 as shown in Figures S16 andS17. In PTB7-Th:h-ITIC the collection efficiency limits the fill factor, in addition to the primary limitation, the field dependent dissociation of CT states.

Field and excess energy dependence of CT state dissociation efficiency. In the following, we examine the influence of an external electric field on the CT state dissociation efficiency of our NFA systems. In addition, we aim to determine whether charge separation depends on the excitation energy of the CT state. If photogeneration proceeds through low energy relaxed CT states, we expect similar field dependence of CT dissociation yield for all the photon energies, as the hot CT states will thermally relax to the lowest lying CT state before separation. If on the other hand, dissociation of CT states is assisted by excess kinetic energy, then an applied external field would have a weaker effect on the hot than on the relaxed CT states.

We use the TDCF technique to investigate the field dependence of charge carrier photogeneration in PTB7-Th:ITIC and PTB7-Th:h-ITIC bulk heterojunction solar cells upon photoexcitation above and below the optical gap, down to the relaxed CT states. TDCF is a transient method with an optical pump and an electrical probe which is sensitive at pump fluences as low as several nJ cm -2 . This allows us to work with low photogenerated charge carrier densities during the experiment to minimise the influence of nongeminate recombination. The photogenerated charge carriers are collected at high reverse bias that further decreases the chance of their recombination. Consequently, when the system is kept at a certain prebias voltage during photogeneration, the field dependence of TDCF signal is governed by the changes in CT dissociation.

TDCF measurements on PTB7-Th:ITIC were performed at six excitation photon energies E ph ranging from 1.21 eV to 1.91 eV. The lower limit of E CT in this NFA system is 1.19 eV, as determined from the temperature dependent V oc measurements in Figure S18. PTB7-Th:h-ITIC devices were excited with E ph of 1.23 eV (lowest energy with detectable TDCF signal), 1.65 eV and 2.33 eV, last two corresponding to donor and acceptor absorption. The optimum experimental conditions for charge photogeneration studies were determined by applying different laser pulse fluences and delay times dur-

!" "!# "!$ %& % " '()' *+, )-'. !" "!/ "!# "!0 1 2 + 3 ! ) 4 * 5 2 6 2 ( ) 7 8 9 : :
) -. ing TDCF measurements (see Figures S19 andS20). The excitation density was adjusted to remain in a charge carrier density range of 1.5 • 10 14 -3.5 • 10 15 cm -3 for all excitations (cf. Figure S21). We measured the total extracted charge as a function of applied prebias V pre during photogeneration in PTB7-Th:ITIC and PTB7-Th:h-ITIC. Using the evaluation method employed by Zusan et al., 35 we quantify the efficiency of CT state dissociation in Figure 3 as η diss = EGE(V pre )/EGE(-2.0 V), where EGE denotes the external generation efficiency. The dissociation is assumed to be equal to unity for PTB7-Th:ITIC at V pre = -2.0 V. In the PTB7-Th:h-ITIC blend the photocurrent does not saturate at -2.0 V and the normalised EGE gives a higher limit of η diss in this system. EGE is calculated as the charge carrier density normalised to the incident photon density, which is estimated from the fluence at the sample at given wavelength.

;<=0%<5>?<?@ !& ),')) !AA),' !&/),')) !BC),' !C$),')) !/ ),' ;<=0%<5>5%?<?@ ) !&C),' ) !$B),' )&!CC),' (a) (b)
The dissociation efficiency η diss in PTB7-Th:ITIC in Figure 3(b) follows the photocurrent at higher photon energies (1.44 eV, 1.53 eV and 1.91 eV). At low photon energies (1.21 eV, 1.29 eV and 1.36 eV) the charge generation yield is lower and has a stronger field dependence, indicating that relaxed CT states in this system dissociate less efficiently. There is no gradual change from one regime to the other, the data rather forms two clusters. These results suggest that there is a limit to photogeneration in PTB7-Th:ITIC below the absorption edge, but this limit is at higher energy than the relaxed CT energy determined by linear extrapolation of V oc (T ) to 0 K. In contrast, the electric field dependence of η diss in PTB7-Th:h-ITIC in Figure 3(a) does not change with excitation energy, indicating that charge separation in this system does not benefit from the excess photon energy.

Stronger bias dependence of the photocurrent in PTB7-Th:h-ITIC as compared to PTB7-Th:ITIC is a sign of a stronger field dependence of CT dissociation at the D-A interface. In the presence of an electric field, geminate recombination is suppressed and the number of free charge carriers in this system increases. Comparing the values at a forward bias of 0.5 V to the ones at a reverse bias of -2.0 V, around 20% of charge carriers are lost to geminate recombination in PTB7-Th:ITIC, while in PTB7-Th:h-ITIC geminate losses are much larger-more than 40%, a sign of less efficient splitting of electronhole pairs at the D-A interface. For both systems the lower photocurrent yield around the maximum power point is most likely originating from the transport resistance caused by imbalanced electron and hole mobilities (see Figures S16 andS17), which is more pronounced in the PTB7-Th:h-ITIC device. All in all, the results from TDCF confirm the findings from TAS and indicate that the difference in molecular structure of ITIC and h-ITIC has a decisive influence on the process of charge separation.

Internal quantum efficiency. The differences in the charge separation yield should be reflected in the IQE of the devices. Recently, Felekidis et al. have used EQE measured under short circuit and reverse bias to approximate IQE spectra of polymer-fullerene solar cells. 43 We apply the same method and obtain excitation energy dependent IQE of PTB7-Th:ITIC and PTB7-Th:h-ITIC by normalising the photocurrent yield to the value at -2 V and -4 V, respectively. Due to low illumination intensities, the influence of nongeminate recombination on the change in photocurrent is minimised and the wavelength dependence of IQE at short circuit is predominantly governed by CT dissociation, as corroborated by the comparison to TDCF at delay times of only 5 ns in Figure S22.

The IQE along with the sensitive EQE data is presented in Figure 4. Both PTB7-Th:ITIC and PTB7-Th:h-ITIC blends share the same optical gap, as shown in Figure S10. For photon energies above this gap, the IQE remains constant for both blends with the mean value of 92% for ITIC and 62% for h-ITIC. In the subgap region, the field dependence of the IQE becomes stronger as the dissociation efficiency drops, verifying the findings from TDCF measurements. In the h-ITIC blend, the IQE becomes photon energy dependent roughly below the optical gap (grey dashed line), while in PTB7-Th:ITIC the IQE decline occurs below 1.44 eV (cyan dashed line) as determined from the normalised IQE spectra at different bias voltages. This transition occurs 0.25 eV above the relaxed CT energy of 1.19 eV.

The EQE of the h-ITIC blend in the subgap region resembles that of neat PTB7-Th, the former being slightly redshifted, probably due to larger energetic disorder in the blend. In the EQE spectra of the ITIC blend, there is a broadening in the subgap region due to additional

! "# ! "$ ! "% ! "& ! "' ! "( ! " ! ! ! ! ( ) * ) + , - . / 0 1 0 2 (0% (0! 0% 0! ) 34 +,562 0! !0% !0! "!0% " 0! 7 * ) + , 2
+89:#"94;4"797<+ +89:#"94;797<+ +=5->+89:#"94+ Figure 4. Sensitive EQE and IQE of PTB7-Th:ITIC and PTB7-Th:h-ITIC. EQE of neat PTB7-Th is plotted for comparison. IQE is calculated from the photocurrent normalised to the value at high reverse bias as in ref. 43 . IQE of h-ITIC blend is photon energy dependent below the optical gap of the donor 1.68 eV (grey dashed line), IQE of ITIC blend is energy dependent below 1.44 eV (cyan dashed line).

states, whether from the acceptor or the charge transfer states. The difference in wavelength dependence of the IQE of both systems below the optical gap suggests that a contribution from these states is the reason for better charge photogeneration in ITIC compared to h-ITIC. At the same time, these lower energy states could be the cause of lower V oc in PTB7-Th:ITIC devices.

III. DISCUSSION

We now consider the charge separation properties of the ITIC and h-ITIC small molecule acceptors by examining their molecular structure. In the absence of a dipole moment, quadrupole moments are the main contribution to the electrostatics of ITIC and h-ITIC molecules. The electrostatic potential induces a shift of energy levels of the highest occupied (HOMO) and the lowest unoccupied molecular orbitals (LUMO) of both acceptors. At the D-A interface the density of acceptor molecules is smaller than in the bulk, which gives rise to the bending of the electrostatic potential. Therefore, the shift of the HOMO and LUMO energy levels is smaller at the interface with a donor, and the energy gap between the CT and the CS states is reduced. 22 The A-D-A architecture of the ITIC molecules provide higher quadrupole moment than h-ITIC molecules. As was shown before for small molecule:C 60 bulk heterojunction solar cells, 21 a larger quadrupole moment leads to a lower energy barrier for free charge formation. The lower electric field dependence of η diss in PTB7-Th:ITIC blend is likely the consequence of such a quadrupole-induced barrier lowering, leading to higher F F and j sc of the solar cell devices. Additionally, charge splitting is more favourable in the A-D-A architecture, due to charge delocalisation along the backbone reducing the Coulomb interaction. Therefore, ITIC molecules with their higher quadrupole moment per unit volume in comparison to h-ITIC are better candidates for charge photogeneration, as verified by our experimental results.

We now discuss the influence of singlet exciton energy in excess of the optical gap on the photogeneration yield. Above the optical gap, the IQE and the field dependence of η diss are independent of the photon energy in both PTB7-Th:ITIC and PTB7-Th:h-ITIC. These findings indicate that photogeneration is not dependent on the singlet exciton excess energy, suggesting that free charge carriers are rather originating from a common intermediate precursor. Consistent with our reasoning, previous investigations of photogeneration pathways in fullerene 34,36 and non-fullerene 25 systems inferred negligible effect of excitation energy above the optical gap on photogeneration. The absence of η diss dependence on the singlet exciton excess energy supports the notion of two-step photogeneration and makes the concept of direct photogeneration from singlet excitons to free charge carriers rather improbable.

Next, we turn to the question whether charge photogeneration proceeds through the energetically lowest CT states or not. We argue that the distinct field dependence of η diss in PTB7-Th:ITIC upon excitation above and below 1.44 eV in Figure 3(b) corresponds to two pools of relatively hot and relatively relaxed CT states. We expect the existence of a threshold energy level (in this case 1.44 eV), above which dissociation of CT states into free charges is more efficient. The electric field dependence of η diss in PTB7-Th:ITIC is comparable to j photo only when exciting the system with photon energies above the threshold. This implies that the main contribution to photogeneration in PTB7-Th:ITIC is from the pool of CT states with energy in excess of the relaxed CT. At excitation below the threshold of 1.44 eV, the pool of relaxed CT states dissociates into free charges less efficiently. Interestingly, the IQE becomes energy dependent also below 1.44 eV. We note that the energy dependence of IQE is only obvious when plotted on the linear scale.

In contrast to PTB7-Th:ITIC, the electric field dependence of η diss in PTB7-Th:h-ITIC is independent of the excitation energy. In view of the low IQE and strongly field dependent charge separation, the photogeneration in the h-ITIC blend is likely originating from the pool of CT states below the threshold. Taken together, these findings indicate that the barrier for free charge formation is too high in this system, even for the states with energy in excess of the relaxed CT. Thus, we propose that the primary limitation of CT separation in NFA organic solar cells stems from the energy barrier between the CT and CS states, which is influenced by the electrostatic potential. Only when this barrier is low enough, as in the case of PTB7-Th:ITIC, separation of CT states can benefit from the excitation energy in excess of the threshold.

To further develop our understanding of the difference in separation of relaxed and hot CT states, we evaluate the parameters influencing η diss directly from TDCF data. Considering a kinetic competition between charge separation and recombination, the yield of CT dissociation is given by

η diss (F ) = k d (F ) k f + k d (F ) (2)
where the dissociation rate constant k d depends strongly on the electric field F , while the rate constant of recombination to the ground state k f is assumed to be field independent. Due to the distribution of separation distances between electrons and holes at the interface, η diss represents the effective dissociation probability. 44 We find that the field dependence of

k d (F ) = η diss (F ) 1 -η diss (F ) • k f (3)
can be described by an exponential function in the electric field range relevant to working conditions of organic solar cells, as indicated by a linear slope in a semi-log plot of k d vs F in Figure 5(a). At higher fields, above ≈ 2.5 • 10 5 V cm -1 , k d (F ) exhibits a second exponential behaviour as η diss (F ) approaches unity. Interestingly, the slope of k d (F ) changes with excitation energy E ph . The commonly used Onsager function for the field dependent increase of k d depends on the Coulomb radius r c , 8,26 i.e. the distance at which electronhole binding energy is equal to the thermal energy. Thus, as the slope of k d changes, this function does not fit the data at all the photon energies with a fixed value of r c .

To describe the electric field dependent increase of the CT dissociation yield we use a model that was previously implemented by Popovic et al. for fitting the field-induced fluorescence quenching in organic semiconductors. 45 According to the model, the probability of dissociation is increased due to an applied electric field F as

k d (F ) = k d (0) exp eF r 0 kT . ( 4
)
where k d (0) denotes the dissociation rate constant at zero field (F = 0), e the elementary charge, k the Boltzmann constant and T the temperature. The parameter r 0 represents a mean initial separation distance between an electron and a hole after thermalisation in the CT manifold and before dissociation into free charges. We plot η diss against field F in Figure 5(b) and fit the data with Equations ( 2) and ( 4). From the fit we obtain the zero field dissociation rate constant k d (0) and thermalisation length r 0 . The field is approximated as F ≈ (V -V oc )/d with d being the active layer thickness. For k f we use the relation k f = τ -1 -k d (0), with τ determined from TAS (see Table 1) and k d (0) from fitting η diss (F ) at lowest excitation pulse energies. The geminate recombination rate constant k f obtained this way is 1.1 • 10 8 s -1 in the ITIC blend and 2.0 • 10 8 s -1 in the h-ITIC blend, respectively.

Given that the effective CT lifetimes τ were similar for both systems, dissociation is expected to be slower 2) and ( 4), plotted against E ph . Binding energies for several r0 are calculated according to E0 = q 2 /4π 0r0 with = 3.5, where 0 is the permittivity of the organic semiconductor material. The vertical dashed line indicates transition between sub-and supergap characteristic field dependence in PTB7-Th:ITIC. in h-ITIC. Indeed, the zero-field dissociation rate constant k d (0) in PTB7-Th:h-ITIC is on average lower than in PTB7-Th:ITIC (Figure 6(a)) and comparable to k f . As seen from the slope in Figure 5(a), the electric field dependence of k d in h-ITIC and low energy subgap ITIC is similar, which implies comparable thermalisation length r 0 . At the lowest excitation energies (1.23 eV and 1.21 eV, respectively) dissociation in both donoracceptor systems is equal as k d (0) = 2.5 • 10 8 s -1 . The evidently much stronger electric field dependence of η diss in h-ITIC at this photon energy in Figure 5(b) is thus primarily a consequence of faster geminate recombination. Consistent with the kinetic Monte Carlo simulations performed by Felekidis et al., 43 the simultaneous blueshift and downshift of the IQE spectrum in h-ITIC relative to ITIC in Figure 4 can also be explained by the higher k f . In the PTB7-Th:ITIC solar cell kinetic energy in excess of the relaxed CT states provides a stronger delocalisation of electron-hole pairs at the D-A interface. As shown in Figure 6(b), the thermalisation distance r 0 between electron-hole pairs in this system, extracted from the fits to η diss (F ), increases from 2.9 nm at 1.21 eV excitation energy to 5.5 nm at excitation of 1.53 eV. The latter agrees well with a recent report of electron-hole separation distance of 5.1 nm in another NFA blend PM6:Y6. 46 Due to the larger separation of electronhole pairs in PTB7-Th:ITIC at higher photon energies, their Coulomb binding energy decreases from 141 meV to 75 meV. Thus, the activation energy for free charge carrier formation is reduced for the CT states with excess kinetic energy. This finding confirms previous reports of CT excess energy assisting charge separation in organic solar cells. 31,32 The drastic rise of k d (0) in PTB7-Th:ITIC at 1.44 eV excitation, where the thermalisation length r 0 is 4.0 nm, might signify the transition between relaxed and hot CT dissociation in the PTB7-Th:ITIC system, supporting our previous argument about the existence of a threshold energy level. Excess energy is driving an electron-hole pair apart during its thermalisation in the CT manifold. At the threshold excitation energy, corresponding to a certain separation distance after the thermalisation process, the CT exciton binding is small enough to generate free charges barrierless.

IV. CONCLUSION

In summary, we have investigated the influence of photon energy and molecular quadrupole moments on the yield of CT dissociation for two non-fullerene acceptor systems, PTB7-Th:ITIC and PTB7-Th:h-ITIC. Quantum chemical calculations reveal that in the more stable dimer configuration the molecular quadrupole moment of the h-ITIC acceptor is 38.8 ea 2 0 , lower than the molecular quadrupole moment of a single ITIC molecule (51.3 ea 2 0 ). In combination with the larger effective van der Waals volume of the h-ITIC dimer, this lower molecular quadrupole moment leads to a weaker electrostatic potential at the donor-acceptor interface. As a consequence, the observed yield of free charge carrier generation in the PTB7-Th:h-ITIC blend is inferior to its counterpart with a large molecular quadrupole moment PTB7-Th:ITIC, as consistently shown by TAS, TDCF and sensitive EQE measurements. In PTB7-Th:ITIC the electric field dependence of the CT dissociation yield decreases with increasing photon energy when exciting the system below the optical gap, while in PTB7-Th:h-ITIC the field dependence of charge separation efficiency is independent of excitation energy. From empirical fits to the TDCF data the separation length of electron-hole pairs in PTB7-Th:h-ITIC is under 3 nm, roughly independent of photon energy. In PTB7-Th:ITIC, the CT separation increases from 2.9 nm to 5.5 nm when exciting the system with photon energies of 1.21 eV and 1.53 eV, respectively. Our results establish that the free charge carrier generation in PTB7-Th:ITIC proceeds through nonrelaxed more delocalised CT states with reduced Coulomb binding of 75 meV. Whereas in PTB7-Th:h-ITIC, due to its lower electrostatic potential and more localised CT states, the barrier for free charge formation is ca. 160 meV even for the CT states with excess energy. [2,3d]thiophen-2-yl)methylene)-3-oxo-2,3-dihydro-1H-inden-1-ylidene)malononitrile (h-ITIC) was synthetised as described in the Supporting Information. All reagents and chemicals from commercial sources were used without further purification unless specified. Solvents were dried and purified using standard techniques. Flash chromatography was performed with analytical-grade solvents using Sigma-Aldrich silica gel (technical grade, pore size 60 Å, 230-400 mesh particle size). Flexible plates Alugram Xtra SIL G UV254 from Macherey-Nagel were used for TLC. Compounds were detected by UV irradiation (Bioblock Scientific). NMR spectra were recorded with a Bruker Avance III 300 ( 1 H, 300 MHz and 13 C, 75MHz) or a Bruker Avance DRX500 ( 1 H, 500 MHz; 13 C, 125 MHz). Chemical shifts are given in ppm relative to TMS and coupling constants J in Hz. IR spectra were recorded on a Bruker spectrometer Vertex 70 and UV-Vis spectra with a Perkin Elmer 950 spectrometer. High-resolution mass spectrometry (HRMS) was performed with a JEOL JMS-700 B/E. Single-crystals of h-ITIC were obtained from solvent diffusion of hexanes in a solution of the named compound in chloroform. X-ray diffraction data were collected and is gathered in the Supporting Information.

Device fabrication and characterisation. The active layer blend solutions were prepared with 20 mg ml -1 concentration of PTB7-Th:ITIC (1:1.2) and PTB7-Th:h-ITIC (1:2.4) in chlorobenzene. The solutions were stirred overnight at 50 • C in the glovebox. Bathocuproine (Ossila) solution (0.5 mg ml -1 ) in methanol was left under stirring at 50 • C in the glovebox overnight, then filtered with 0.2 µm PTFE membrane filter and used directly afterwards. Solar cells were fabricated in regular architecture on glass substrates coated with prepatterned indium tin oxide (ITO). The substrates were cleaned in ultrasonic bath with detergent, acetone, isopropanol and deionised water. The substrates were then dried with nitrogen and exposed to low-pressure oxygen plasma for 5 min. After that, a 35 nm layer of poly(3,4-ethylenedioxythiophene) polystyrene sulfonate (PEDOT:PSS, Clevios AI 4083) was spin-coated on ITO and annealed at 140 • C for 10 min. The substrates were then transferred to a nitrogen-filled glovebox, where they were heated again for 10 min. The blend solutions were subsequently spin-coated to yield an active layer thickness of ca. 100 nm, followed by the filtered bathocuproine solution giving a layer thickness of 8 nm. The devices were completed by 100 nm of thermally evaporated Ag with a base pressure below 10 -6 mbar through shadow masks with active areas of 4.0 mm 2 for steady-state measurements and 0.5 mm 2 for TDCF. Current-voltage characteristics were measured using a Keithley 236 source measure unit. White light LED intensity was matched to the short circuit current calculated from the EQE spectra. The EQE measurements were performed using Bentham TM300 monochromator and a Si reference photodiode. The output current was measured using Stanford Research SR830 DSP lock-in amplifier.

Gas-phase quantum chemical calculations. DFT calculations were done using B3LYP functional and 6-311g(d,p) basis set as implemented in the GAUSSIAN package. Singlet and dimer states were extracted from the available crystal packing of ITIC and h-ITIC unit cells (see ref. 47 and Figure S9). Obtained molecule configurations were used to calculate dipole and quadrupole moments. See Supporting Information for more details.

Transient absorption spectroscopy. TAS was carried out using a previously described custom pump-probe setup. 48 The output of a titanium:sapphire amplifier (Coherent LEGEND DUO, 4.5 mJ, 3 kHz, 100 fs) reduced to 2 mJ/pulse by splitting, pumped an optical parametric amplifiers (OPA) (Light Conversion TOPAS Prime) to generate 1300 nm pulses. This beam was in turn used to seed a calcium fluoride (CaF 2 ) crystal mounted on a continuously moving stage, thereby generating a white-light super continuum from 350 to 1100 nm. The excitation light (pump pulse) was provided by an actively Q-switched Nd:YVO 4 laser (InnoLas picolo AOT) frequency-doubled to provide ca. 800 ps wide pulses at 532 nm. The pump laser was triggered by an electronic delay generator (Stanford Research Systems DG535) itself triggered by the transistor-transistor logic (TTL) sync from the Legend DUO, allowing control of the delay from pump and probe with a jitter of roughly 100 ps. Pulse length and repetition rate allowed for the delay between slightly less than 1 ps to 300 µs. The sample was kept under a dynamic vacuum of < 10 -5 mbar in a cryostat (Optistat CFV, OXFORD Instruments). The transmitted fraction of the white light was guided to a custom-made prism spectrograph (Entwicklungsbüro Stresing) where it was dispersed by a prism onto a 512 pixel NMOS linear image sensor (Hamamatsu S8380-512DA). The probe pulse repetition rate was 3 kHz, while the excitation pulses were directly generated at 1.5 kHz frequency, while the detector array was read out at 3 kHz. Adjacent diode readings corresponding to the transmission of the sample after excitation and in the absence of an excitation pulse were used to calculate ∆T /T . Measurements were averaged over several thousand shots to obtain a good signal-to-noise ratio. The delay at which pump and probe arrive simultaneously on the sample (i.e. zero time) was determined from the point of maximum positive slope of the TA signal rise, which is expected to correspond to the maximum of the pump pulse.

Time delayed collection field. For TDCF, the sample was excited monochromatically using a Light Conver-sion PHAROS femtosecond laser (400 µJ, 5 kHz, 290 fs, 50/50 beam splitter) together with an Light Conversion ORPHEUS optical parametric amplifier enabling an output beam of variable wavelength 315 nm to 2600 nm. The repetition rate of a laser was reduced to 5 kHz by a pulse picker. To ensure monochromatic irradiation, longpass filters were used to block the residual radiation emitted by the laser. The fluence arriving at the sample was set by a combination of two Thorlabs FW102C/FW112C filter wheels. Fluence and laser stability were monitored with a Newport 818-BB-40 photodetector. During the experiment, the sample was mounted in a self-built sample holder. Prebias and collection bias voltages were produced by a Keysight 81160A pattern generator, its output voltage was applied to the sample through a unity gain amplifier. The current from the solar cell was measured as a voltage drop over a 10 Ω resistor connected in series to the device. This voltage drop was amplified using a differential amplifier and subsequently measured with a GaGe CS121G2 digitizer. During photogeneration experiment the time delay of electrical probe was 5 ns. Collection bias was at least -4.0 V for efficient extraction of all the remaining mobile charge carriers. Total extracted charge for all the excitation wavelengths was of the same order of magnitude (Figure S21). Fluences were varying from 3 nJ cm -2 to 13 µJ cm -2 , being higher when exciting sub-gap states.

Sensitive EQE. Sensitive EQE was measured with a home-built monochromatic excitation setup. A MSH-300D double monochromator (LOT Quantum Design) with an 100 W quartz-tungsten halogen lamp was used as monochromatic light source, with additional optical bandpass filters to reduce stray light from higher harmonic wavelengths. The photocurrent from the solar cells was measured with a Zurich HF2LI lock-in-amplifier in combination with a variable Zurich HF2TA transimpedance amplifier. External bias to the solar cell was provided by an output of the latter, where applicable. A small fraction of the monochromatic light was recorded using a Hamamatsu K1718-B two-color photodiode involving a Silicon and an Indium-Gallium-Arsenide detector, and measured with a SR-DSP830 lock-in amplifier and the second input of the HF2LI, respectively. During the measurements, the solar cells were held in a Helium atmosphere in a custom built closed-cycle liquid helium cryostat (Cryostat). Sensitive EQE measurements were performed at bias voltages ranging from -3.0 V to +0.5 V for PTB7-Th:ITIC and from -4.0 V to +0.5 V for PTB7-Th:h-ITIC. Photocurrent of PTB7-Th:ITIC saturated at -2.0 V. In the case of PTB7-Th:h-ITIC there was no saturation, and EQE at -4.0 V gives a higher limit of the calculated IQE.

Figure 1 .

 1 Figure 1. (a) Chemical structures of acceptor molecules. (b) j(V ) curves of PTB7-Th:ITIC and PTB7-Th:h-ITIC solar cells in the dark and under one sun illumination intensity. (c) Mutual arrangement of acceptor molecules in a single and dimer configurations and corresponding values of dipole (D) and quadrupole moments (Q total 20 -total quadrupole moment of a system, Q mol 20 -quadrupole moment per molecule). Some of the side chains were removed to simplify the representation.

  Figure 1(b) reveal that h-ITIC blend devices suffer from lower short-circuit current (j sc ) and fill factor (F F ), but have higher open circuit voltage (V oc ) than ITIC devices. The extremely low fill factor in PTB7-Th:h-ITIC solar

Figure 2 .

 2 Figure 2. Kinetics of charge carriers in (a) PTB7-Th:h-ITIC and (b) PTB7-Th:ITIC (symbols) at different fluences and global fits (solid lines) according to Equation (1).

Figure 3 .

 3 Figure 3. Dissociation yield η diss from TDCF (symbols) as a function of prebias voltage Vpre for different excitation photon energies for (a) PTB7-Th:h-ITIC and (b) PTB7-Th:ITIC, compared to j photo at one sun illumination intensity, normalised to the value at -2.0 V (black solid lines).

Figure 5 .

 5 Figure 5. (a) k d (F ) of PTB7-Th:ITIC at E ph of 1.29 eV, 1.44 eV and 1.91 eV, and PTB7-Th:h-ITIC at E ph of 1.65 eV (the rest omitted for clarity). The fits with Equation (4) (solid lines) reveal an exponential dependence of k d (F ). The fit to k d (F ) of ITIC at 1.29 eV is scaled to illustrate the difference in the slopes (cyan dashed line). The Onsager-Braun model has fixed slope with rc = const (grey dashed line). (b) η diss (F ) of PTB7-Th:ITIC and PTB7-Th:h-ITIC, and fits to the data with Equations (2) and (4).

Figure 6 .

 6 Figure 6. (a) k d (0) and (b) r0 of CT states in PTB7-Th:h-ITIC and PTB7-Th:ITIC obtained from the fits to η diss (F ) with Equations (2) and (4), plotted against E ph . Binding energies for several r0 are calculated according to E0 = q 2 /4π 0r0 with = 3.5, where 0 is the permittivity of the organic semiconductor material. The vertical dashed line indicates transition between sub-and supergap characteristic field dependence in PTB7-Th:ITIC.

V

  . EXPERIMENTAL SECTION Materials. Poly[4,8-bis(5-(2-ethylhexyl)thiophen-2yl)benzo[1,2-b;4,5-b ]dithiophene-2,6-diyl-alt-(4-(2-ethyl hexyl)-3-fluorothieno[3,4-b]thiophene-)-2-carboxylate-2-6-diyl)] (PTB7-Th) and 3,9-bis(2-methylene-(3-(1,1dicyanomethylene)-indanone))-5,5,11,11-tetrakis(4-hexyl phenyl)-dithieno[2,3-d:2 ,3 -d ]-s-indaceno[1,2-b:5,6-b ]di thiophene (ITIC) were purchased from 1-Material. (Z)-2-(2-((9,9-bis(4-hexylphenyl)-9H-indeno[1,2-b]thieno

ACKNOWLEDGMENTS

The work of M.S., P.S.M. and C.D. has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie Grant Agreement No. 722651 (SEPOMO). The work of C.W. and C.G. was funded by the DFG (project DE 830/19-1). The work of S.K., F.L. and J.G. was supported by funding from King Abdullah University of Science and Technology (KAUST), S.K. specifically acknowledges the Office of Sponsored Research (OSR) under Award No: OSR-2018-CARF/CCF-3079 and Award No. OSR-CRG2018-3746. D.A. acknowledges KAUST for funding his sabbatical stay. A.M. acknowledges the financial support received from Marie Curie Individual Fellowship (SMOLAC). The authors thank the MATRIX SFR of the University of Angers.

Appendix A: Corresponding Author Email: deibel@physik.tu-chemnitz.de

Appendix B: Notes

The authors declare no competing financial interest.