Yu L Dokshitzer

G Marchesini

Monte Carlo and large angle gluon radiation

We discuss the problem of incorporating recoil effects into the probabilistic QCD evolution scheme based on the picture of colour dipoles as done in recent Monte Carlo programs. Such a scheme correctly describes subleading soft contributions to multiplicity distributions. However we find that a simple receipt for incorporating recoil effects into the dipole multiplication picture based on the ordering of gluon energies conflicts the collinear factorization and does not lead to the correct DGLAP equation.

Preface

We share opinion of our late friend Bo Anderson who once bitterly complained about two aspects of scientific publishing that became standard in recent decades: 1) shying away from expressing emotions and 2) not making public an outcome of a study that produced a negative, undesired results. The present paper is moderately emotional; at the same time it is perfectly in line with the second Bo's demand: it reports the finding that the authors did not plan, and would rather not want, to arrive at.

The answer to the question that is posed in this paper suggests that a formulation of Monte Carlo event generation based on the multiplication picture of energy ordered colour dipoles is plugged by a serious problem.

Introduction

Generation of events using Monte Carlo methods is an indispensable tool for planning, running and analysing the results of modern high energy experiments. A possibility to generate multi-particle production as a Markov chain of successive independent parton splittings is based on the general property of factorization of collinear singularities which is characteristic for logarithmic field theories (quantum field theories with dimensionless coupling).

Event generation is necessarily an approximate procedure. Strictly speaking, in order to predict the cross section dσ n for production of n partons in a hard interaction with the hardness scale Q 2 , one needs to plug in the QCD matrix element of the (at least) n th order of the perturbation theory. Instead, one repeats n times the basic 1 → 2 splitting process of the first order in α s . By so doing one correctly samples the major part of the n-parton phase space, though not all of it. It is this part that gives the dominant, maximally collinear enhanced contribution to, e.g., the inclusive particle distribution, D (n) = O (α n s log n Q 2 ). This approximation known as leading logarithmic (LLA) can be systematically improved by including higher order parton splitting processes. Thus, the next-to-leading logarithmic accuracy corresponding to the series α n+1 s log n Q 2 is achieved by incorporating 1 → 3 parton splittings, etc. For inclusive parton distributions this leads to the DGLAP evolution equations [1] whose generalisation to multi-particle distributions can be achieved in the spirit of jet-calculus [2], with additional account of soft gluon coherence.

Beyond the one loop "evolution Hamiltonian" it starts to matter how one organises the parton chain, that is, what one chooses as "parton evolution time". In particular, due to the presence of soft gluons one encounters corrections of the type (α s log 2 x) k which are formally subleading (non-collinear) but become explosively large in the x ≪ 1 limit (here x is the gluon energy fraction). This eventuality can be cured by an appropriate tuning of the evolution time parameter. Namely, for time-like parton multiplication (jets)one has to choose the angular ordering [3]. This choice corrects the "naive" fluctuation time ordering (dictated by an examination of Feynman denominators in a cascading process) by taking into clever account collective effects leading to destructive interference in soft gluon radiation (for more details see, e.g., [4][5][START_REF] Yu | Basics of Perturbative QCD[END_REF]).

An alternative way of dealing with soft gluon interference effects for time-like parton multiplication (jets) is provided by the "dipole scheme" [4,[START_REF] Gustafson | [END_REF] in which an independently radiating parton, 1 → 2, is replaced by a colourless dipole formed by two partons neighbouring in the colour space, 2 → 3. Gluon radiation off a dipole is automatically suppressed at angles exceeding the dipole opening angle thus reproducing the angular ordering. Colour suppressed contributions O (1/N 2 c ) that lie beyond the scope of the dipole approach are numerically small and difficult to trigger experimentally [8].

Within the standard coherent parton cascade picture [3][4][5][START_REF] Yu | Basics of Perturbative QCD[END_REF] it is the small-angle multiplication processes populating jets that enjoy full all-order treatment (get "exponentiated"). The dipole formulation offers a possibility to improve this treatment by taking into consideration logarithmically enhanced effects due to multiple emission of soft gluons at large angles with respect to jets. Non-collinear soft gluons dominate inter-jet particles flows in various hard processes. They also complicate the analysis of the so-called nonglobal QCD observables [9][10][11], i.e. in observables in which recorded radiation is confined in geometrically definite phase space regions. It is then interesting to involve these corrections into a Monte Carlo, code based on dipole emission, including soft radiation awy from jets.

In this paper we discuss a dipole based scheme (in the large-N c approximation) which is well suited for deriving improved analytic predictions for observables like mean multiplicities, inclusive soft particle spectra and energy flows, correlations, etc., that incorporate large-angle soft gluon radiation effects as in [9][10][11]. At the same time, we observe that once one aims at beyond the no-recoil (soft) approximation, treating colour dipoles as independently evolving entities is likely to conflict the collinear factorization. It does not lead to the correct DGLAP equation.

The paper is organised as follows. We start by constructing in Section 2 an evolution equation for the generating functional that describes soft gluon content of e + e -annihilation events. No recoil is included. This equation is based on the factorization of the multi-soft gluon distribution, in the planar limit, and uses the centre-of-mass energy as an "evolution time". Contrary to what is done normally, no angular variable enters the evolution time here. The evolution equation in the integral form we interpret as a Markov process of successive dipole branchings, which interpretation leads to a Monte Carlo chain process for probabilistic event generation. As a first example, we solve the evolution equation for mean gluon multiplicity to analyse subleading corrections due to large-angle soft gluon radiation. Here we work in the soft approximation meaning: 1) small energies of emitted gluons and 2) neglecting the energy-momentum loss (recoil) by the primaries.

These two aspects of the "soft gluon approximation" do not necessarily go together. In Section 3 we consider an inclusive energy distribution of the leading quark in the quasi-elastic limit, (1x) ≪ 1, where emitted gluons ought to be soft but the quark recoil is essential. At the beginning, in the construction of the evolution equation no recoil effects were included. They are incorporated later by using the recoil strategy that seems natural for a dipole multiplication scheme, see [12]. We find, however, that the quark fragmentation function so obtained deviates essentially from the collinear resummation result given by the DGLAP evolution equation. This is in contrast to what is happening in Monte Carlo schemes based on angular ordered time-like cascades where collinear singularities are correctly resummed [5,13].

Discussion of the problem encountered is relegated to the Conclusions section.

Multiple soft gluons and Monte Carlo

The method of generating functionals makes it straightforward to generate exclusive events using Monte Carlo techniques. We start from the known n-gluon distribution in the soft limit and construct the evolution equation for the corresponding generating functional.

The generating functional for a quark-antiquark pair q q plus an ensemble of n secondary partons, γ * → p a p b q 1 , q 2 , . . . q n , is

G(Q, u) = ∞ n=0 dΦ (n) ab (q 1 , q 2 , . . . , q n ) u a u b n i=1 u(q i ) M (n) ab (q 1 , . . . , q n ) 2 .
(2.1)

Here p a , p b are the final on-mass-shell momenta of the quark and antiquark produced by the virtual photon γ * and q i the momenta of the secondaries. The matrix element M

(n) ab describes production of n + 2 real final state partons and contains any number of virtual ones. dΦ (n) is the full phase space factor of the ensemble of n + 2 massless particles:

dΦ (n) ab (q 1 , q 2 , . . . , q n ) = (dp a )(dp b ) n i=1 (dq i ) (2π) 4 δ (4) (p a + p b + n i=1 q i -Q), (2.2)
where (dq) = d 3 q 2ω (2π) 3 , ω = q 0 = |q| .

The "source functions" u are attached to each parton in (2.1). They help to extract an arbitrary final state observable once the functional is known. The fully inclusive measurement, that is when one allows for production of any number of particles with arbitrary momenta, corresponds to setting all u = 1. This gives G(Q, u = 1) = 1 corresponding to normalization to the total cross section. In this section we work in the no-recoil approximation.

Generating functional in the soft limit

In the large-N c ("planar") approximation the squared matrix element can be approximated at tree level as

M (n) ab (q 1 , . . . , q n ) 2 ≃ (4π 2 ᾱs ) n n! perm W ab (q i 1 , . . . , q in ) • M (0) ab 2 ᾱs = N c α s π ,
(2.3a) where W ab is the "multiple antenna" function [4] describing the production probability of n soft real gluons q i :

W ab (q 1 , q 2 , . . . , q n ) = (p a p b ) (p a q 1 )(q 1 q 2 ) • • • (q n p b ) . (2.3b)
This result holds in the soft gluon approximation q i ≪ p a , p b . In invariant terms, (q i p a ), (q i p b ) ≪ (p a p b ), (2.4) or, equivalently, (q i , p

a + p b ) ≪ 2(p a p b ) = Q 2 giving ω i ≪ E a ≃ E b ≃ Q/2, (2.5)
with E a , E b and ω i the energies in the centre-of-mass of the "dipole" (ab). In this approximation the phase space can be simplified as

dΦ (n) ab (q 1 , . . . , q n ) ≃ (dp a )(dp b )(2π) 4 δ 4 (p a + p b -Q) • n i=1 (dq i ) ≡ dΦ (0) ab • n i=1 (dq i ) . (2.6)
The expression (2.3) was derived implying strong ordering of gluon energies:

ω in ≪ ω i n-1 ≪ . . . ≪ ω i 2 ≪ ω i 1 ≪ Q/2.
The answer, however, is symmetric with respect to momenta of participating gluons. Therefore, we may ignore the ordering between gluons and set a common upper bound E = 1 2 Q for gluon energies. Extracting the quarks whose momenta in the soft gluon approximation are not affected by radiation, we may write

G(u; Q) ≃ dΦ (0) ab u(p a )u(p b ) |M (0) ab | 2 • G(p a , p b ; 1 2 Q). (2.7)
The functional G(p a , p b ; E) describes production of any number of gluons with energies limited from above by some value E, off the dipole formed by the quarks with momenta p a , p b . By setting E = 1 2 Q we obtain the multi-gluon generating functional describing e + e -annihilation process with s = Q 2 .

The generating functional contains proper virtual factors for each n-gluon contribution. Omitting for the moment the virtual corrections, we have

G (real) (p a , p b ; E) ≡ n W ab (q 1 , . . . , q n ) n i=1 ω i dω i dΩ 4π ᾱs • u(q i ) • ϑ(E -ω i ) . (2.8)
Due to the symmetric structure of the antenna function (2.3) one may choose W ab (q 1 . . . q n ) to represent any energy ordered ensemble of gluons and drop the symmetry factor 1/n!.

Evolution equation

To obtain an evolution equation for the generating functional we need to exploit the structure of the multi-gluon antenna distribution (2.3) and to deduce a recurrence relation. To construct a recurrence relation one should select an "evolution variable" to order emitted gluons. In the logic of the collinear approximation this is done by ordering the angles of successive gluon emissions, which ordering takes full care of the destructive interference contributions in the soft region [3][4][5][START_REF] Yu | Basics of Perturbative QCD[END_REF] and preserves the probabilistic parton multiplication picture. In particular, this procedure was applied to the expression (2.8) to construct the HERWIG event generator [5,13]. However, as was already stated, the distribution (2.3) is valid in the soft limit for arbitrary angles, and one is tempted to use this property in order to lift off the collinear approximation. The derivation of the evolution equation avoiding the small angle approximation was done in [9][10][11] where the upper bound E on gluon energies in the event centre-of-mass was treated as the "evolution time" parameter. Let us recall the corresponding construction employing the energy as an evolution variable.

Energy as an evolution variable

We take the ordering in the energy of emitted gluons in the center-of-mass of the event (2.5) and deduce an evolution equation for the soft generating functional. The evolution equation follows from the exact recurrence relation -the factorization property of the multi-gluon antenna function (2.3b) -which reads

W ab (q 1 , q 2 , . . . , q n ) = W ab (q ℓ ) • W aℓ (q 1 , . . . q ℓ-1 ) W ℓb (q ℓ+1 , . . . , q n ), (2.9)
where ℓ is any gluon. To construct the evolution equation one takes ℓ to be the most energetic among the gluons.

Introducing the scaled antenna function

W ab (q) ≡ ω 2 • W ab (q) = ξ ab ξ aq ξ bq , ξ ik = 1 -cos Θ ik ,
which depends only on the angles between partons, we obtain

E∂ E G(p a , p b ; E) = dΩ 4π ᾱs W ab (q) u(q) G(p a , q; ω)G(q, p b ; ω)-G(p a , p b , ω) ω=E . (2.10)
Here the subtraction term takes care of virtual corrections, so that taking all u(q) = 1 one derive the desired normalization

G(p a , p b ; E)| u=1 = 1 . (2.11)
No recoil was included up to now; it will be considered later.

Integral evolution equation and Monte Carlo

By exponentiating the total one-gluon emission probability we obtain the Sudakov form factor which describes the probability that the dipole does not radiate gluons with energies up to a given value E:

S(p a , p b ; E) = exp - E dω ω dΩ 4π ᾱs (q ab t ) W ab (q) ϑ(q ab t -Q 0 ) . (2.12)
Here q ab t is the invariant transverse momentum of the gluon with respect to the pair of quarks:

(q ab t ) 2 = 2(p a q)(p b q) (p a p b ) = 2 W ab (q) , (2.13)
and an arbitrary parameter Q 0 has been introduced as a collinear cutoff. By constructing the logarithmic derivative of the Sudakov form factor over the maximal gluon energy,

E∂ E ln S(p a , p b ; E) = - dΩ 4π ᾱs (q ab t ) W ab (q) ϑ(q ab t -Q 0 ) ω=E , (2.14)
and plugging this expression into (2.10) we can trade the virtual subtraction term for the Sudakov factor. By so doing we arrive at the equivalent equation for the soft generating functional in the integral form:

G(p a , p b ; E) = S(p a , p b ; E) + dω ω dΩ 4π u(q) • S(p a , p b ; E) S(p a , p b ; ω) • ᾱs (q ab t ) W ab (q) ϑ(q ab t -Q 0 ) • G(p a , q; ω) G(q, p b ; ω).
(2.15)

Iteration of this equation can be interpreted as parton branching which can be realised as a Monte Carlo (Markov) process. To see this, observe that the iteration kernel

dP(p a , p b , q; E) = dω ω dΩ 4π S(p a , p b ; E) S(p a , p b ; ω) • ᾱs (q ab t ) W ab (q) ϑ(q ab t -Q 0 ) (2.16)
can be written as the product of two probability distributions (here the boundary

q ab t > Q 0 is implicit): dP(p a , p b , q; E) = dP (1) ab (ω, E) • dP (2)
ab (Ω q ) , P

ab (ω, E) = S(p a , p b ; E) S(p a , p b ; ω) , dP (1)
ab (Ω q ) = ᾱs

dΩ q 4π ξ ab ξ aq ξ qb ᾱs dΩ 4π ξ ab ξ aq ξ qb -1 .
(2.17)

In the Monte Carlo process one tries to generate emission of a gluon with momentum q off a given dipole. The first distribution, dP (1) , provides the gluon energy ω, if the transverse momentum bound is satisfied, and the second one, dP (2) , its direction Ω q . If the boundary is not satisfied, the dipole does not emit. One has

dP(p a , p b , q, E) = 1 -S(p a , p b ; E), (2.18)
which shows that the Sudakov factor gives the probability of not emitting a soft gluon within the resolution q ab t > Q 0 .

Mean multiplicity

By construction the generating functional (2.15) embeds only soft gluon radiation. Incorporating double logarithmic (simultaneously soft-and collinear-enhanced) contributions, it also contains single logarithmic effects due to emission of soft gluons at large angles.

Here we demonstrate an application of the dipole-based evolution to the calculation of subleading large-angle soft gluon corrections. In the collinear approach, such corrections are treated as due to "multi-jet" configurations (contributing to the "coefficient function") rather than jet evolution ("anomalous dimension"), see [START_REF] Yu | Basics of Perturbative QCD[END_REF].

Consider mean gluon multiplicity. In order to obtain an equation for the multiplicity of secondary partons, one applies to (2.10) the variational derivative over the probing function u(k) and integrates over k, while setting u ≡ 1 for all remaining probing functions (one-particle inclusive measurement):

E∂ E N(ξ ab , E) = dΩ 4π ᾱs (k ab t ) W ab (q) 1 + N(ξ qa , E) + N(ξ qb , E) -N(ξ ab , E) ω=E . (2.19) Introducing N = 1 + N we have [14] N ′ (ξ ab ) = E∂ E N (ξ ab ) = dΩ 4π ᾱs (k ab t ) W ab (q) N (ξ qa ) + N (ξ qb ) -N (ξ ab ) . (2.20)
In order to illustrate the difference with the standard approach based on the angular ordering we will analyse (2.20) in the double logarithmic approximation which neglects hard parton splittings and recoil effects. For the time being we will ignore the running of the coupling; the corresponding corrections will be addressed later.

Representing the integrand on the r.h.s. of (2.20) as

ξ ξ 1 ξ 2 N (ξ 1 )+N (ξ 2 )-N (ξ) = 1 ξ 2 ξ ξ 1 N (ξ 1 ) -N (ξ) + 1 ξ 1 ξ ξ 2 N (ξ 2 ) -N (ξ) - ξ -ξ 1 -ξ 2 ξ 1 ξ 2 N (ξ),
and using the 1 ↔ 2 symmetry, we obtain

N ′ (ξ) = dΩ 4π ᾱs 2 ξ 2 ξ ξ 1 N (ξ 1 ) -N (ξ) - ξ -ξ 1 -ξ 2 ξ 1 ξ 2 N (ξ) . (2.21)
We can perform the integration over the azimuth of q around the direction of p a ,

dΩ 4π = 1 2 2 0 dξ 1 • dφ 1 2π .
Using the relation 1

ξ 2 ≡ dφ 1 2π 1 ξ 2 = 1 |ξ -ξ 1 | , we get ξ -ξ 1 -ξ 2 ξ 1 ξ 2 = 1 ξ 1 ξ -ξ 1 |ξ -ξ 1 | -1 = - 2 ξ 1 ϑ(ξ 1 -ξ). (2.22)
As a result,

N ′ (ξ) = ᾱs 2 0 dξ 1 |ξ -ξ 1 | ξ ξ 1 N (ξ 1 ) -N (ξ) + ᾱs 2 ξ dξ 1 ξ 1 N (ξ), (2.23)
with the E dependence implicit. The first contribution we split into integrals over small and large angles, 0 ≤ ξ 1 ≤ ξ and ξ < ξ 1 ≤ 2. Introducing the integration variable η = ξ 1 /ξ ≤ 1 in the first region and η = ξ/ξ 1 in the second, 1 2 ξ ≤ η ≤ 1, the equation (2.23) takes the form [14]

N ′ (ξ) = 1 0 dη ᾱs 1 -η 1 η N (ηξ) -N (ξ) + 1 1 2 ξ dη ᾱs 1 -η N η -1 ξ -N (ξ) . (2.24)
This result we represent as a sum of two terms,

E∂ E N (ξ) = 1 0 dη ᾱs η N (η ξ) + ∆(ξ), (2.25a)
where

∆(ξ) = 1 0 dη ᾱs 1 -η N (η ξ) -N (ξ) + 1 1 2 ξ dη ᾱs 1 -η N η -1 ξ -N (ξ) . (2.25b)
The first integral term on the r.h.s. of (2.25a) generates the standard DL anomalous dimension for the mean multiplicity; the second one constitutes a subleading correction.

The main term

We have to solve the equation

E∂ E N DL (ξ; E) = 1 0 dη ᾱs η N DL (η ξ; E) (2.26) with the initial condition N (ξ; E)| E √ 2ξ=Q 0 = 1. The solution is a function of a single variable Q ξ ≡ E √ 2ξ
(the maximal transverse momentum of partons radiated by the dipole).

The argument of the running coupling in (2.26) can be approximated as

ᾱs (k ab t ) = ᾱs E 2ξ 1 ξ 2 ξ ≃ ᾱs E 2ξ 1 .
The correction to this approximation turns out to be negligible since the logarithmic factor ln(ξ 2 /ξ) vanishes in the collinear limit ξ 1 → 0:

∝ dξ 1 ξ 1 β 0 ᾱ2 s • ln ξ 2 ξ = O ᾱ2 s .
Thus we obtain the double differential equation

d 2 d ln Q ξ 2 N DL (Q ξ ) = 2 ᾱs (Q ξ ) N DL (Q ξ ), Q ξ = E 2ξ, (2.27)
with the initial conditions N (Q 0 ) = 1 and N ′ (Q 0 ) = 0. Its solution for the fixed coupling case reads

N DL (Q ξ ) = cosh γ ln Q ξ Q 0 ≃ 1 2 Q ξ Q 0 γ , (2.28)
where the anomalous dimension γ is given the correct expression

γ = ᾱs 1 0 dη η • η 1 2 γ = ᾱs • 2 γ =⇒ γ = √ 2 ᾱs .

The correction term

Let us analyse the correction term (2.25b) keeping the contributions of the order of ᾱ3/2

s while neglecting contributions O (ᾱ 2 s
). The answer depends on the value of the opening angle ξ.

Small opening angles (ξ ≪ 1). It is the region of small opening angles that gives a dominant contribution to the integral determining the DL anomalous dimension in (2.26). In this entire region the correction (2.25b) does not contribute at the ᾱ3/2 s level. Indeed, substituting the solution (2.28) into the expression (2.25b) and setting ξ → 0 in the lower limit of the second integral we get

∆(ξ; E) N (ξ; E) = ᾱs 1 0 dη 1 -η η γ/2 + η -γ/2 -2 = ᾱs • 2ψ(1) -ψ(1 + 1 2 γ) -ψ(1 -1 2 γ) = ᾱs • 1 4 γ 2 ψ ′′ (1) + . . . ≃ 1 2 ᾱ2 s ψ ′′ (1). ( 2

.29)

Large opening angle (ξ ≃ 2). At the same time, ∆ contributes to the multiplicity of the "fully open" dipole, ξ = 2. The expression (2.25b) for ∆(2; E) contains one integral:

∆(2; E) = ᾱs 1 0 dη 1 -η N (η • 2; E) -N (2; E) .
Expanding the difference of the multiplicity factors,

N (η • 2; E) -N (2; E) = E∂ E N (2; E) • 1 2 ln η + . . . , we derive ∆(2; E) ≃ E∂ E N (2; E) • ᾱs 1 0 dη 1 -η • 1 2 ln η = - π 2 12 ᾱs • E∂ E N (2; E).
Integrating this correction over energy we finally obtain an order α s correction to the dipole multiplicity due to emission of two soft energy ordered gluons:

N (2; E) = N DL (2; E) • 1 - π 2 12 ᾱs + . . . . (2.30)
In the standard approach, this correction belongs to the coefficient function and originates from a "four-jet" configuration with two large-angle energy ordered soft gluons as additional "jets" [START_REF] Yu | Basics of Perturbative QCD[END_REF].

Other observables sensitive to large-angle gluon effects

There is a variety of observables that are sensitive to soft gluon radiation at large angles. Among them one can mention production of heavy quark pairs in jets [14], non-global jet observables [9] including various particle distributions and correlations related with particle flows in the inter-jet regions [10,11,[15][16][17][18][19]. The dipole-based evolution equation (2.15) is well suited for taking into account single logarithmic correction effects due to multiple soft gluons in these observables and elsewhere. At the same time, it remains insufficient for building a realistic Monte Carlo event generator which is impossible without incorporating full parton decay probabilities, including hard parton splittings, and of the recoil effects to ensure energy-momentum conservation.

Attempt to include recoil in dipole multiplication

Generally speaking, a gluon emitted by a dipole cannot be ascribed as offspring to either of the two partons that form the dipole. However, in the collinear limit when the gluon momentum q is quasi-parallel to one of the two hard partons, say, p a , it is the parton a that can be said to "independently split" into two, sharing its momentum, P a → p a + q. A small change in the "spectator" momentum P b → p b which is necessary to compensate for the virtuality of the (a, q) pair, vanishes in the collinear limit. A general recoil strategy has to be formulated in such a way that in the collinear limit the answer reduces to the standard DGLAP parton splitting function P a → p a + q that includes hard momentum configurations, p a ∼ q. This is a source of a significant single-logarithmic correction, this time not from soft gluon radiation at large angles but from hard collinear-enhanced emission.

Thus, in order to properly formulate a recoil strategy, one must split the soft dipole radiation function W ab (q) into two pieces, W (a) ab (q) and W (b) ab (q), each of which incorporates the collinear singularity when q || p a or q || p b , respectively:

W ab (q) = (p a p b ) (p a q)(qp b ) = W (a) ab (q) + W (b) ab (q), (3.1a)
or, graphically,

q p a p b = + (3.1b)
One can do this differently, for example by introducing [12] W (a)

ab (q) = (p a p b ) (p a q)(p a + p b , q) , W (b)
ab (q) = (p a p b ) (p b q)(p a + p b , q) . (3.2)
Another possibility to split the dipole is given by the "conditional probabilities" [5,[START_REF] Yu | Asymptotic Freedom and Local Parton-Hadron Duality[END_REF]]

W (a) ab (q) = W ab (q) • 1 2 1 + (qp a )(p b Q) -(qp b )(p a Q) (p a p b )(qQ) = 1 2ω 2 ξ qa 1 + ξ ab -ξ qa ξ qb . (3.
3) The distribution (3.3) produces the exact angular ordering upon averaging over the azimuthal angle φ of the gluon momentum q around the singular direction p a :

W (a) ab φq,p a = 1 ω 2 ξ qa • ϑ ξ ab -ξ aq . (3.4)
Note that in all cases the distributions W (c) ab (q) (c = a, b) fall fast and essentially become irrelevant when the emission angle exceeds the opening angle of the parent dipole, ξ qc > ξ ab , that is away from the angular ordered kinematics.

Recoil strategy

Consider the elementary process

P a + P b =⇒ p a + p b + q .
One ascribes definite recoil pattern separately to the two terms in (3.1). Thus, for the first term W (a) ab which is collinear singular in the direction q || p a , we choose, following the Catani-Seymour prescription [12] p It defines the two light-like vectors P a , P b that represent the momenta of the parent partons prior to the gluon emission:

p (a) a = zP a + (1 -z)yP b -k t , (3.5a) q = (1 -z)P a + zyP b + k t , (3.5b)
where k t is ortogonal to P a and P b . The on-mass-shell condition p 2 a = 0 gives (P a + yP bq) 2 = 0 =⇒ y = (qP a ) (P aq, P b ) .

(3.6)

The light-cone fraction z of the parent momentum P a carried by the final quark p a is given by the expression 1z = (qP a ) (P a P b ) .

In the collinear limit y → 0. In the soft limit both y → 0 and z → 1.

We now recast the phase space element in terms of the "parent momenta":

dΦ ab (q) = (dp a )(dp b )(dq)(2π) 4 δ (4) (p a + p b + q -Q) = (dP a )(dP b )(2π) 4 δ (4) (P a + P b -Q) • (dq) × J (a) ab (q) ≡ dΦ (0) ab • (dq) J (a) ab (q). (3.7)
The Jacobian of this transformation reads

J (a) ab (q) = 1 - (qP a ) (P a -q, P b ) 1 - (qP b ) (P a P b ) -1 = 1 -y z . (3.8)
The gluon radiation probability we represent as a sum of two contributions:

dΦ ab (q) • W ab (q) ⇒ dΦ (0) ab • (dq) W (a) ab (q) J (a) ab (q) + W (b) ab (q) J (b) ab (q) . (3.9)
Finally, the evolution equation becomes

E∂ E G(P a , P b ; E) = dΩ 4π ᾱs (q ab t ) ϑ(q ab t -Q 0 ) × c=a,b J (c) ab (q) W (c) ab (q) u(q) G(p (c) a , q; ω) G(q, p (c) b , q; ω) -G(P a , P b ; ω) , (3.10a)
which again satisfies the normalization G(P A , P b ; E) = 1 for all u(q) = 1 (cf. (2.11)).

In the graphic form it can be represented as

+ p b p a = b c E d dE - a - (3.10b)
This evolution equation is the modification of (2.10) that includes recoil as described above. The gluon radiation function W ab (q) = ω2 W ab (q) is split according to (3.9) into the two pieces which are collinear singular when q is parallel to p a or to p b , correspondingly, see [12]. By introducing the Sudakov form factors this equation can be recast in an integral form as it has been done in the soft case in (2.15). An iterative solution of this integral equation defines a Monte Carlo process for generating multi-parton ensembles with account of collinear non-enhanced single-logarithmic corrections due to large-angle soft gluon emission. Now we are going to check if the dipole recoil scheme described above is consistent with known analytical results concerning the collinear resummation. The simplest observable of this type is an inclusive energy distribution of the final-state quark.

The non-singlet quark energy distribution

The non-singlet quark fragmentation function is obtained by taking the derivative of the generating functional G ab with respect to the quark source u(p a ) and setting all the remaining source functions u = 1. This way we obtain the distribution in the momentum of the quark a accompanied by any number of gluons. If the quark energy is taken large, (1x) ≪ 1, all radiated gluons are soft and the analysis simplifies significantly. The problem becomes essentially Abelian1 and is described by multiple independent radiation of soft gluons by the quark.

Collinear approximation

We shall restrict ourselves to configurations in which all radiated gluons have small emission angles with respect to the quark direction. This -quasi-collinear -approximation is sufficient for the analysis of the anomalous dimension which accumulates collinear ("mass") singularities of the fragmentation function in all orders and describes the scaling violation.

In this kinematics the measured variable x is given by

x = 1 - i y i , ω i = y i E , p a ≃ x P a . (3.11)
What matters is the total energy 1-x carried, on average, by an ensemble of radiated soft gluons.

Looking for collinear singular contributions -the terms of the type ᾱs ln(1x) • ln Q 0 -we may omit the two last terms on the r.h.s. of the equation (3.10b) that contain the factor W (b) ab (q) which is non-singular when the gluon momentum q becomes collinear to p a .

One is then left with contributions involving only the singular antenna piece W Since our observable is inclusive, we may use unitarity arguments to simplify the equation. Indeed, by the nature of the adopted recoil strategy, production of gluons by the "upper" dipole (bc) in (3.12) does not affect the momentum of the quark a.

Therefore, since such emissions do not affect the measurement, they are subject to real-virtual cancellation and can be neglected.

One is left to consider only gluons emitted the lower blob involving the parton a:

= - c c a a b a E d dE (3.13)
This equation generates multiple gluon emissions off the quark. By construction, successive emissions are ordered in gluon energies. And here comes the crucial observation: as we have discussed above, the emission angle of the gluon c is essentially limited from above by the aperture of the parent dipole, ξ ab , see (3.4). Therefore, soft gluons generated by the evolution equation (3.13) turn out to be ordered, simultaneously, in energies and angles with respect to the radiating quark. Instead we know that the DGLA equation is obtaine from ordering only in collinear variables, disregarding the relative energies of emitted partons. The corresponding equation in the present dipole-formulation looks as follows.

Evolution equation with dipole recoil strategy

Since the radiated gluons are soft, we use the distribution W (a) ab (q) defined in (3.3) which describes the soft part of the gluon radiation probability (splitting function). We may also approximately set J (a) ab (q) = 1. In going from the generating functional G(P a , P b ; E) in (3.10) to the single particle distribution, we express the latter as a function of the angle between the partons ξ ab , the energy scale E ≃ E a and the energy fraction of the triggered quark: D(ξ ab , E, x).

Performing integration over the azimuthal angle φ aq (see (3.4)) one deduces from (3.10)

D( ξ; E; x) = δ(1 -x)+ 1 dy y ξ dξ 2ξ ᾱs ϑ yE ξ -Q 0 × D ξ, yE, x 1 -y -D ξ, yE, x . (3.14)
Here ξ = ξ ab is the angular aperture of the original dipole (ab) while ξ = ξ ac is the aperture of the secondary (ac)-dipole. Equivalently one may write2

D( ξ; E; x) =δ(1 -x) S(E ξ) + 1 dy y ξ dξ 2ξ ᾱs ϑ yE ξ -Q 0 S(E ξ) S(yE ξ) D ξ, yE, x 1 -y . (3.15)
The distribution D in the integrand is evaluated at ξ since a gluon produced by the dipole (ac) has an angle with respect to a bounded above by ξ ac , see (3.4). One has

ln S(E ξ) = - 1 dy y ξ dξ 2ξ ᾱs ϑ yE ξ -Q 0 . (3.16)
This equation correctly satisfies the sum rule for the first moment: integrating over x one obtains 1 for the number of quarks a (to see this one uses (2.18)). However, it is in conflict with the Poisson nature of the multiple soft gluon radiation and contradicts the collinear resummation results leading to the DGLAP equation. The origin of this failure can be understood already at the level of two emitted gluons.

Two gluon emission

Consider the emission of two soft gluons p 1 , p 2 off the parton q of the q q dipole. The antenna functions that potentially contribute in the collinear limit (recall that we keep all gluon angles with respect to the quark a = q to be small) are displayed here:

q q q q q q (c) 2 (b) 1 2 (a) 2 1 1 (3.17)
The first two graphs correspond to the splitting of the dipole (a1):

W (a) a1 (2) → ϑ(ξ a1 -ξ a2 ) and W (1)
a1 (2) → ϑ(ξ a1 -ξ 12 ),
while the third one, W

1q (2) → ϑ(ξ 1qξ 12 ), is the relevant part of the large-aperture dipole (1q). Due to the local recoil prescription used, only the contribution (a) affects the momentum of the quark q. In the two remaining ones, (b) and (c), the gluon 2 borrows its energy-momentum from the gluon 1 and does not produce any quark recoil. Therefore, these contributions cancel against corresponding virtual corrections in the inclusive quark measurement. In conclusion, within the adopted recoil strategy, only the graph (a) should be kept, and we obtain the following phase space for the two-gluon emission:

ξ a2 ≪ ξ a1 ≪ ξ aq and ω 2 < ω 1 ≪ E. (3.18)
The first condition comes from the angular ordering in the graph (3.17)(a), and the second condition -from the energy ordering of successive emissions. We know, however, that in order to obtain the DGLAP equation that properly resums collinear singular contributions, one needs to assemble angular (or, transverse momentum) ordered emissions (the first ordering), regardless to the order of gluon energies (the second one). At the same time, the dipole logic is leading us to the double-ordered gluon ensemble, according to (3.18). As a result, instead of a simple Poisson distribution of soft radiation (1/n!) we obtain something like a modified Bessel function distribution (1/(n!) 2 ).

What is missing here is actually the coherence of QCD radiation. As well known, a soft gluon ω 2 , with ω 2 ≪ ω 1 ≪ E a , could be emitted at large angles (ξ a2 ≫ ξ a1 ) directly by the original parton q + p 1 ≃ q. In the language of Feynman amplitudes, such radiation occurs as a coherent sum of the graphs (b) and (c).

Another way to discuss this point is the following. The dipole radiation pattern of soft gluon 2 off the composite antenna (a1q) can be represented as follows: (2). (3.19a) Beyond the large-N c approximation, the latter term (corresponding to emission of the gluon 2 by the quark-antiquark dipole) acquires a N 2 c -suppressed correction,

W a1 (2) + W 1q (2) = W a1 (2) + W 1q (2) -W aq (2) + W aq
N c → N c -1/N c : N c 2 W a1 (2) + W 1q (2) =⇒ N c 2 • W a1 (2) + W 1q (2) -W aq (2) + C F • W aq (2). (3.19b)
This is in line with the collinear factorization theorem. Indeed, the factor in square brackets is non-singular when p 2 becomes parallel to the direction of the quark or antiquark momentum so that these collinear enhancements are contained by the quark dipole term W aq (2) and must be proportional to the "colour charge" (Casimir operator) C F of the quark. The only collinear singularity of the first contribution on the r.h.s. is that of p 2 parallel to p 1 which suggest an interpretation of this term as independent emission of the gluon 2 by the previous generation (harder) gluon 1. The full two soft gluon emission probability takes the form

dw(1, 2) ∝ C F N c 2 • W aq (1) W a1 (2) + W 1q (2) -W aq (2) + C 2 F • W aq (1)W aq (2). (3.20)
In the large-N c limit, this is identical to the original dipole expression. At the same time, the representation (3.20) suggests a different recoil strategy when it comes to incorporating the energy-momentum conservation.

The first term in (3.20) describes cascade gluon multiplication: q q → 1 followed by 1 → 2. In the inclusive quark distribution this is subject to the real-virtual cancellation. Energy ordered configurations, ω 2 ≪ ω 1 cancel leaving no trace; "hard" gluon splittings with ω 2 ≃ ω 1 modify one gluon emission by giving rise to α s (k 2 1⊥ ) for the radiation intensity of the primary gluon 1 by the quark.

The second term describes independent radiation of two gluons off the quark -the second order term of the Poisson distribution of primary gluons -which is necessary to correctly recover the DGLAP quark fragmentation function.

Equation (3.20) provides a probabilistic representation for the two-gluon production in terms of a combination of independent and cascade gluon emissions. With account of hard parton splittings and of the natural energy recoil prescription, it served as the base for constructing the HERWIG parton event generator [5,13]. The corresponding probabilistic representation for production of three and four energy ordered gluons was developed in [START_REF] Yu | Asymptotic Freedom and Local Parton-Hadron Duality[END_REF] and described in [START_REF] Yu | Basics of Perturbative QCD[END_REF], which resulted in a formulation of the parton multiplication scheme based on the exact angular ordering of successive parton splittings.

The difference between the two strategies becomes transparent. In the dipole scheme, successive gluons borrowing the energy from the quark were ordered in emission angles, W

aq (1)

W (a) a1 (2) , ξ 2a < ξ 1a < 1,
while in (3.20) the angle of the softest gluon is not bounded, and the two emissions are totally independent,

W (a) aq (1) W (a) aq (2) , ξ 2a < 1, ξ 1a < 1.
Let us see how this has happened, from the point of view of the gluon-gluon splitting.

In the quasi-collinear configuration, ξ a1 ≪ 1, the radiation pattern of the softest gluon 2 consists of a "narrow" and a "wide-angle" dipoles, (a1) and (1q). When, within the dipole recoil strategy, the "wide-angle" antenna W 1q (2) is split into the pieces collinear singular with respect to the directions of the first generation gluon (1) and of the antiquark,

W 1q (2) = W (1) 1q (2) + W (q) 1q (2), the contribution W (1)
1q (2) is looked upon as radiation of the gluon 2 by the gluon 1, whichever the kinematical configuration of the gluon 1 with respect to its parent quark (dipole factorization). At the same time, within the HERWIG logic (collinear factorization), the gluon-gluon multiplication occurs only in a restricted angular aperture, ξ 21 < ξ 1a . Indeed, in the kinematical region ξ 1a ≪ ξ 21 ≪ 1 the two last terms in the combination

W a1 (2) + W 1q (2) -W aq (2) on the r.h.s. of (3.19) cancel: lim ξ 1a →0 ξ 2 1a • W 1q (2) -W aq (2) = O ξ 1a ξ 12 → 0
(while the first antenna term, W a1 (2), is restricted in angle on its own). Physically, the gluon 2 radiated off the gluon 1 in the "wide" dipole (3.17)(c) at relatively large angles, ξ 21 > ξ 1a , is produced by the grandparent parton a (quark). Correspondingly, it must borrow energy from the quark rather than locally from the dipole (1q) to which it belongs from the point of view of the colour connection topology.

Reformulating this way the local dipole recoil prescription modifies the energy ordered evolution equation. The modification reduces to substituting the full dipole opening angle ξ for the running angle parameter ξ in the argument of the distribution D on the r.h.s. of previous figure. As demonstrated in the Appendix, this way one recovers the approximate DGLAP equation for the quark fragmentation function.

Conclusions

Monte Carlo generation of QCD events is a quarter century old business, based on the structure of resummation of collinear enhanced Feynman diagram contributions. The probabilistic parton cascade picture was first established for one-particle inclusive quantities (DGLAP evolution equations for DIS structure functions and e + e - fragmentation functions, [1]), generalised via "jet-calculus" [2] and, with account of the soft gluon coherence [3][4][5][START_REF] Yu | Basics of Perturbative QCD[END_REF], has laid the base for probabilistic description of the internal structure of parton jets and their ensembles.

Collinear-non-enhanced ("large angle") soft gluon radiation provides significant NNLL corrections to global event characteristics (e.g., mean particle multiplicity); it is also responsible for inter-jet multiplicity and energy flows and determines the structure of various non-global observables [9][10][11]. Effects of multiple soft gluon radiation at large angles lie beyond the scope of the standard (collinear) approach and must be treated order by order in perturbation theory (while collinear enhanced contributions are resummed in all orders). An elegant expression [4] for the multiple soft gluon production probability (2.3) is valid for arbitrary angles and offers a possibility of improving the parton picture. The structure of multi-gluon distribution (2.3) naturally suggests an interpretation in terms of a chain of colour connected dipoles. By choosing an "evolution time" variable this chain may be generated via a Markov process of successive dipole splittings. An evolution time can in principle be chosen differently. We have used energy ordering of gluons for this purpose. The generating functional that we have constructed with the help of the corresponding evolution equation allows one to calculate specific effects due to multiple emission of soft gluons at large angles in the large-N c approximation.

In order to construct a realistic Monte Carlo generator for multi-parton ensembles it is imperative, however, to formulate an adequate recoil prescription which would ensure energy-momentum conservation at every successive step of the parton (dipole) multiplication.

In the present paper we addressed the question, whether the "dipole factorization" extends beyond the no-recoil approximation. In other words, whether the splitting of a colour dipole into two can be treated independently of the prehistory of the system, that is locally in the evolution time (which is a necessary condition for constructing a Markov process).

Having taken an inclusive energy distribution of the final quark as the simplest example of a collinear sensitive observable, we have shown that a naive implementation of the dipole recoil strategy results in violating the collinear factorization.

However, the existing Monte Carlo implementations of the dipole picture (for recent examples see for instance [21][22][23]) chose to order successive gluon emissions in (dipole centre of mass) transverse momenta rather than in (laboratory frame) gluon energies. Such a construction apparently prevents the problem under consideration from appearing in the leading order.

Indeed, as we have seen above, the problem with the quark recoil appears when a softer gluon ( 2) is emitted at a large angle than its predecessor (1). The inclusive quark energy distribution is sensitive to the region where the two gluons have comparable energies. Then, a softer one that carries a slightly smaller energy but is radiated at a large angle, acquires a larger transverse momentum, and therefore, by construction, is looked upon as "the first" by the k ⊥ -ordered MC scheme. Therefore, "the second" gluon in such a picture, that is the one with smaller transverse momentum, is automatically having a smaller emission angle; in the leading logarithmic approximation the coherence is respected.

The problem of an interplay between energy-momentum recoil and soft gluon coherence discussed in this paper gets postponed but should be kept in mind in view of attempts at constructing the next-to-leading order QCD MC generators. Now the gluons are no longer ordered in energies (the gluon q i does not need to be softer then q i-1 ) but in the emission angle with respect to the quark. In the angular ordered form, a soft gluon radiated at a relatively large angle can be said to be emitted before any other (harder or softer) gluons that move in a collinear bunch around the quark. This is in accord with the colour coherence according to which a large-angle soft gluon "sees" only the total colour charge of a collinear group of partons.

To derive the DGLAP equation we replace

dω i ω i = dy i y i = dz i 1 -z i , 1 -
y i = z 1 • • • z n (A.5a)
where now 1z i is the local fraction of energy taken away by the gluon q i from the parent parton (quark), while y i is the gluon energy fraction with respect to E, the energy of the primary parton P a :

y i = (1 -z i )z i-1 • • • z 1 . (A.5b)
We finally obtain the evolution equation

D(Q, x) = δ(1 -x)S(Q) + Q 2 Q 2 0 dq 2 t q 2 t 1- Q 0 q t 0 dz ᾱs (q t ) 2(1 -z) • S(Q) S(q t ) D q t , x z = δ(1 -x) + Q 2 Q 2 0 dq 2 t q 2 t 1-
Q 0 q t 0 dz ᾱs (q t ) 2(1z)

• D q t , x z -D(q t , x)

(A.6)
with Q = E √ ξ and q t = yE √ ξ. The upper boundary of the z-integration, 1 -Q 0 /q t , can be replaced by 1 in the logarithmic collinear approximation, which results in the integral equation that correctly gives (the soft piece of) the DGLAP anomalous dimension.

Finally, let us mention another obstacle that one faces when trying to construct the full realistic Monte Carlo scheme starting from the energy ordered gluon ensembles.

In order to match the full DGLAP anomalous dimension, one has to include the "hard part" of the splitting function by adding to the distributions W (c) ab (q) the pieces vanishing at ω = 0. In the equation (A.6) this calls for the replacement

dy i y i = dz i 1 -z i =⇒ dz i 1 1 -z i -1 2 (1 + z i ) . (A.7)
However, while we could cast the leading soft distribution in (A.1) in terms of the local energy fraction, dy i /y i = dz i /(1z i ), the finite term (1 + z i ) dz i cannot be expressed locally, via the single variable y i , see (A.5).

  (a) b = (1y)P b .

Non-Abelian effects due to final state cascading of gluons reduce to making the effective coupling in the radiation probability of a primary gluon k run with k

⊥ .

The y integration in (3.15) is actually limited by y = 1x since the distribution D(ξ, Q, z) ∝ Θ(1z)

Acknowledgements

We thank the organisers and participants of the CERN TH jet phenomenology seminar for the opportunity to extensively discuss the issues raised in the paper. We are especially grateful to Walter Giele, Gösta Gustafson, David Kosower, Leif Lönnblad, Zoltan Nagy, Gavin Salam, Jim Samuelsson, Mike Seymour, Torbjorn Sjöstrand, Peter Skands, Davidson Soper and Bryan Webber for illuminating discussions.

A. Global recoil and correspondence with DGLAP

We have no universal recipe under sleeve for extending the dipole gluon multiplication picture beyond the no-recoil approximation. However, for a simple example of the inclusive quark distribution discussed in Section 3 the situation is straightforward to cure.

As we saw from the previous discussion, the origin of the failure lay in the fact that the softest gluon emitted at angles larger than the angle of the previous branchings did not contribute to the quark recoil. So, it suffices, by brute force, to permit the gluons of all generations to contribute to the quark recoil, irrespectively to the value of their emission angles.

This amounts to replacing in the integrand of (3.15) the distribution D evaluated at ξ = ξ ac with the distribution evaluated at the full opening angle ξ = ξ ab . One then has the different equation,

which actually corresponds to the DGLAP evolution. To see this, we iterate equation (A.1) to obtain the series

is the product of with theta-functions which ensures the energy ordering of gluons. Now, due to the symmetry of the multiple integral (A.2) with respect to energy and angular variables, we can replace the energy ordering with the angular ordering:

Θ en-ord =⇒ Θ ang-ord = Θ(ξ n < • • • < ξ 1 < ξ). (A.4)


Antonio De Nicola

Andrea Correa

Andrea Giunchi

Luca Muccioli

Gabriele D ' Avino

Ju- Nji Kido

Giuseppe Milano

Bidimensional H-Bond Network Promotes Structural Order and Electron Transport in BPyMPs Molecular Semiconductors: A Computational Study

The presence of a hydrogen bond (H-Bond) network has been proved to impact significantly the OLED devices efficiency by promoting molecular orientation and structural anisotropy in thin films. The design of specific compounds to control H-Bond network formation in an amorphous material, and hence to improve OLED performances, is needed. A successful example is given by the bi-pyridyl based family n-type of organic semiconductors named BPyMPM. The experimental evidences demonstrated a surprisingly higher electron mobility in thin film composed of B4PyMPM (B4), which is almost two order of magnitude higher than mobility measured for very similar member of the family, B2PyMPM (B2). We present a comprehensive computational study, where classical and ab-initio methods were combined, to investigate the 2D H-Bond network in B4 and B2 thin films. The results indicate that B4 forms a larger number of intermolecular C-H•••N hydrogen bonds which promote a higher orientational and positional order in B4 films and superior electron transport properties.

INTRODUCTION

The Organic Light Emitting Diode research field has constantly grown in the last 35 years: Scopus web archive counts roughly 12000 publications from 1994 to 2020, while the number of publications per year is still increasing. 1 Nowadays OLEDs dominate the display market and show promise for a similar success in lighting applications in the near future, owing to the high brightness and the low power consumption they can offer. New generation OLEDs are produced by vapor-deposition of organic materials which form layers, in which the molecular orientation is documented to be a key property capable to improve outcoupling and charge transport efficiency. [START_REF] Marcato | Molecular Orientation Effects in Organic Light-Emitting Diodes[END_REF][START_REF] Boehm | The Interplay of Interfaces, Supramolecular Assembly, and Electronics in Organic Semiconductors[END_REF][START_REF] Ohisa | Molecular Interdiffusion between Stacked Layers by Solution and Thermal Annealing Processes in Organic Light Emitting Devices[END_REF][START_REF] Watanabe | Simultaneous Manipulation of Intramolecular and Intermolecular Hydrogen Bonds in N-Type Organic Semiconductor Layers: Realization of Horizontal Orientation in OLEDs[END_REF][START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF][START_REF] Watanabe | Control of Molecular Orientation in Organic Semiconductor Films Using Weak Hydrogen Bonds[END_REF][START_REF] Duhm | Orientation-Dependent Ionization Energies and Interface Dipoles in Ordered Molecular Assemblies[END_REF][START_REF] Frischeisen | Increased Light Outcoupling Efficiency in Dye-Doped Small Molecule Organic Light-Emitting Diodes with Horizontally Oriented Emitters[END_REF][START_REF] Liehm | Comparing the Emissive Dipole Orientation of Two Similar Phosphorescent Green Emitter Molecules in Highly Efficient Organic Light-Emitting Diodes[END_REF] The vapor deposition process allows to produce smooth films, with a control of film thickness on nanometer length scale, by condensing molecules from the gas phase onto a temperature-controlled substrate, [START_REF] Ediger | Anisotropic Vapor-Deposited Glasses: Hybrid Organic Solids[END_REF] but offers much less control on molecular orientation and crystallinity, which are governed by a complex interplay between kinetic and thermodynamic effects, often substrate-dependent. [START_REF] Ediger | Perspective: Highly Stable Vapor-Deposited Glasses[END_REF][START_REF] Roscioni | Pentacene Crystal Growth on Silica and Layer-Dependent Step-Edge Barrier from Atomistic Simulations[END_REF] A viable strategy to maximize the thermodynamic thrust is to design molecules with specific shapes and functional groups that promote highly directional supramolecular interactions.

In this context, Prof. Yokoyama and co-workers demonstrated the general relationship between the structural anisotropy of a molecule and its orientational order and preferred alignment direction in a thin solid film. [START_REF] Yokoyama | Horizontal Molecular Orientation in Vacuum-Deposited Organic Amorphous Films of Hole and Electron Transport Materials[END_REF][START_REF] Yokoyama | Horizontal Orientation of Linear-Shaped Organic Molecules Having Bulky Substituents in Neat and Doped Vacuum-Deposited Amorphous Films[END_REF] In another series of studies, intermolecular C-H•••N hydrogen bonds were exploited for driving the supramolecular assembling, by designing molecules containing multiple pyridine rings properly arranged in order to maximize intermolecular interactions that promote long-range structural order. References missing?

This strategy was showed to be successful and several molecules were demonstrated to form flat surfaces and regular stacking through self-assembly. [START_REF] Watanabe | Simultaneous Manipulation of Intramolecular and Intermolecular Hydrogen Bonds in N-Type Organic Semiconductor Layers: Realization of Horizontal Orientation in OLEDs[END_REF][START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF][START_REF] Watanabe | Control of Molecular Orientation in Organic Semiconductor Films Using Weak Hydrogen Bonds[END_REF][START_REF] Sasabe | Influence of Substituted Pyridine Rings on Physical Properties and Electron Mobilities of 2-Methylpyrimidine Skeleton-Based Electron Transporters[END_REF] Moreover, the presence of a bidimensional H-bonds network promotes anisotropic optical properties, like the refractive index and extinction coefficient k, and a non-negligible electron mobility. [START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] Remarkably, a difference of about two orders of magnitude in the electron mobility, from time-of-flight measurements, was found between vacuum deposited film composed of B4 (~10 $% 𝑐𝑚 ( 𝑉 $* 𝑠 $* ) and B2 (~10 $, 𝑐𝑚 ( 𝑉 $* 𝑠 $* ), respectively, [START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] whose chemical structure is reported in Figure 1.

The understanding of the driving forces regulating molecular orientation is mandatory to optimize material properties, such as charge transport, and to design processes and/or compounds to improve device efficiency. In this matter, the molecular dynamics (MD) at atomic resolution is a powerful computational method able to rationalize the microscopic origin of the molecular preferred orientation. [START_REF] Yoo | A Molecular Dynamics Study on the Interface Morphology of Vapor-Deposited Amorphous Organic Thin Films[END_REF][START_REF] Youn | All-Atom Simulation of Molecular Orientation in Vapor-Deposited Organic Light-Emitting Diodes[END_REF][START_REF] Olivier | Computational Design of Thermally Activated Delayed Fluorescence Materials: The Challenges Ahead[END_REF][START_REF] Hu | Origin of High Efficiencies for Thermally Activated Delayed Fluorescence Organic Light-Emitting Diodes: Atomistic Insight into Molecular Orientation and Torsional Disorder[END_REF][START_REF] Muccioli | Simulation of Vapor-Phase Deposition and Growth of a Pentacene Thin Film on C60(001)[END_REF][START_REF] Walters | Influence of Molecular Shape on the Thermal Stability and Molecular Orientation of Vapor-Deposited Organic Semiconductors[END_REF][START_REF] Antony | Influence of Vapor Deposition on Structural and Charge Transport Properties of Ethylbenzene Films[END_REF][START_REF] Lee | The Molecular Origin of Anisotropic Emission in an Organic Light-Emitting Diode[END_REF][START_REF] Tonnelé | Elucidating the Spatial Arrangement of Emitter Molecules in Organic Light-Emitting Diode Films[END_REF][START_REF] Han | Deposition Growth and Morphologies of C60 on DTDCTB Surfaces: An Atomistic Insight into the Integrated Impact of Surface Stability, Landscape, and Molecular Orientation[END_REF] Several studies have been reported to successfully study interfacial properties and molecular arrangement of numerous compounds used in OLED materials. [START_REF] Yoo | A Molecular Dynamics Study on the Interface Morphology of Vapor-Deposited Amorphous Organic Thin Films[END_REF][START_REF] Youn | All-Atom Simulation of Molecular Orientation in Vapor-Deposited Organic Light-Emitting Diodes[END_REF][START_REF] Hu | Origin of High Efficiencies for Thermally Activated Delayed Fluorescence Organic Light-Emitting Diodes: Atomistic Insight into Molecular Orientation and Torsional Disorder[END_REF][START_REF] Yoneya | Self-Assembly of Donor-Acceptor Semiconducting Polymers in Solution Thin Films: A Molecular Dynamics Simulation Study[END_REF] For a general description of the challenges and the strategies currently used in the computational modelling of OLED materials, we refer the interested reader to a recent perspective [19] and the references therein.

In the present study, density functional theory (DFT) calculations and all-atom MD simulations are combined to investigate the molecular orientation of B4 and B2 in thin films as final aim. The role of the hydrogen bond network and its consequences on the molecular geometries and electrical properties are investigated. The paper is organized as follows: in the Methods and Models section the atomistic model description, forcefield optimization, simulation protocol and computational details are reported. The structural and charge properties calculated from MD simulations are discussed and compared with experiments in the Results and Discussion section, while the main results are summarized in the Conclusion section.

METHODS AND MODELS

2.1 Molecular Dynamics Simulations. All the atomistic simulations reported in the paper were performed by using the GROMACS 5.0 simulation package. [START_REF] Berendsen | A Message-Passing Parallel Molecular Dynamics Implementation PROGRAM SUMMARY Title of Program: GROMACS Version 1.0[END_REF][START_REF] Hess | GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation[END_REF] A time step of 2 fs has been employed for all simulations. The temperature, in both NVT and NPT ensembles, has been held constant at 298 K by using a stochastic velocity rescaling algorithm [START_REF] Bussi | Canonical Sampling Throught Velocity Rescaling[END_REF] with a characteristic coupling time 𝜏 1 = 0.02 ps; moreover, in the NPT ensemble the pressure has been kept constant at 1 atm by using the Berendsen barostat [START_REF] Berendsen | Molecular Dynamics with Coupling to an External Bath[END_REF] (pressure coupling time 𝜏 2 = 0.2 ps, compressibility of 1•10 -5 bar -1 , anisotropic scaling of box sides with the same 𝜏 2 value for x, y, z directions). All the bonds involving hydrogen atoms have been constrained using the LINCS algorithm. [START_REF] Hess | A Linear Constraint Solver for Molecular Simulations[END_REF] Lennard-Jones interactions were calculated by using a cut-off distance of 1.4 nm, while electrostatic interactions have been treated with the Particle Mesh Ewald algorithm [START_REF] Darden | Particle Mesh Ewald: An Nlog(N) Method for Ewald Sums in Large Systems[END_REF] with a grid spacing of 0.1 nm. System compositions and total simulation times are listed in Table 1. 2 Atomistic Model. We adopted the atomistic OPLS force field (FF), developed by Jorgensen in 1996, [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF] except for soft aryl-aryl torsions. These torsions and their energy barriers are fundamental in determining the molecular shape and must be described accurately for a correct modeling. [START_REF] Muccioli | Supramolecular Organization of Functional Organic Materials in the Bulk and at Organic/Organic Interfaces: A Modeling and Computer Simulation Approach BT -Multiscale Modelling of Organic and Hybrid Photovoltaics[END_REF] For instance, a recent extension of the OPLS-AA has been optimized to improve the description of inter-ring torsions [START_REF] Dahlgren | Characterization of Biaryl Torsional Energetics and Its Treatment In[END_REF] for molecules containing biaryl groups. Following these guidelines, we made a reparameterization of the Fourier coefficients used in the Fourier series sum describing the torsional energy profile (see eq. S4 of Supporting Information). The torsional energy profiles for φ .* and φ .( dihedrals were calculated (see Figure 1 for their definition). The geometries of B4 and B2 molecules have been optimized at different fixed aryl-aryl angles (with steps of 10º), at PBE0/6-311G** level, 37-39 using Gaussian09. [START_REF] Frisch | Gaussian09 Revision D[END_REF] By fitting the torsional energy profile from DFT calculations of φ .* and φ .( dihedral angles, Fourier coefficients are obtained. The comparison between torsional energy profiles at DFT and FF level is reported in the Supporting Information together with the full list of FF parameters and the functional forms of bonded and nonbonded interactions.

Sample Preparation.

The thin films of B4 or B2 compounds, prepared via vacuum deposition, possess to some extent an ordered structure, presumably due to a columnar arrangement favored by the presence of a bidimensional H-bond network. [START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] Since the structural information in the solid state of B4 and B2 is not available, we built thin films initial configurations starting from the experimental crystal structure of a molecule, B4PyPTZ, having optical properties and chemical structure similar to B4 and B2. The crystallographic unit cell [START_REF] Watanabe | Simultaneous Manipulation of Intramolecular and Intermolecular Hydrogen Bonds in N-Type Organic Semiconductor Layers: Realization of Horizontal Orientation in OLEDs[END_REF] of B4PyPTZ, containing 16 molecules, has been replicated 2 × 2 × 1 times along in a, b, c directions to build a supercell composed of 128 molecules.

The crystalline structure of B4PyPTZ is characterized by molecules stacked in columns. Each molecule belonging to a column is oriented in a way leaving the nitrogen atoms of terminal pyridine rings pointing toward the exterior of the column. The peculiar position of nitrogen atoms maximizes the probability to form of intermolecular H-bonds between adjacent columns, forming a bidimensional H-bond network. [START_REF] Watanabe | Simultaneous Manipulation of Intramolecular and Intermolecular Hydrogen Bonds in N-Type Organic Semiconductor Layers: Realization of Horizontal Orientation in OLEDs[END_REF] The stability of B4PyPTZ sample has been checked performing a MD simulation in the NPT ensemble employing the OPLS-AA force field with the tuned torsional potential for 𝜑 1* , 𝜑 1( angles. The equilibrium mass density obtained from the MD simulation (1278 kg/m 3 ± 5 kg m -3 ) slightly underestimates the experimental density 5 of about 3%. However, the deviation from the density experimental value is similar to those observed for the majority of organic molecules (2 -4%) modelled by using OPLS force field. [START_REF] Dahlgren | Characterization of Biaryl Torsional Energetics and Its Treatment In[END_REF][START_REF] Chen | Thermodynamic Properties of the Williams, OPLS-AA, and MMFF94 All-Atom Force Fields for Normal Alkanes[END_REF][START_REF] Jorgensen | Development of an All-Atom Force Field for Heterocycles[END_REF][START_REF] Bernardes | Evaluation of the OPLS-AA Force Field for the Study of Structural and Energetic Aspects of Molecular Organic Crystals[END_REF][START_REF] Moral | Cost-Effective Force Field Tailored for Solid-Phase Simulations of OLED Materials[END_REF] Besides, the lattice constants (a, b, c) of crystalline structure are reproduced within an error of +1.9%. The time behavior of the mass density and the X-Ray pattern calculated on the equilibrium structure are reported in the Supporting Information.

The equilibrium structure of the B4PyPTZ supercell was used as initial coordinate set to build B4 and B2 systems. The description of the procedure employed to build systems and further details about it are reported in the Supporting Information.

NPT simulations revealed a difference in the equilibrium mass density between the two bulk systems, even if the only difference between the two molecules lies in the position of nitrogen atoms in the pyridine rings (Figure 1). Starting from initial configurations, at the same density for both molecules, the bulk of B4 reaches a slightly higher equilibrium density (1197 ± 5 kg m-3) than B2 (1176 ± 5 kg m-3) revealing a more efficient packing (time behavior of mass density is reported in Figure S8 of Supporting Information).

2.4 Electronic Structure Calculations. Most organic semiconductors with mobilities below ~1 cm 2 /Vs are characterized by thermally activated hopping-like transport. The most popular, semiclassical rate equation for charge hopping is due to Marcus [START_REF] Marcus | Electron Transfer Past and Future[END_REF] (Eq. 1), stating a proportionality between the rate between two molecules 𝑖 and 𝑗 and the square of the intermolecular electronic coupling 𝐽 ?@ , and an Arrhenius-like exponential dependence with activation barrier of A∆𝐺 ?@ + 𝜆F ( /(4𝜆):

𝑘 ?@ = (M ℏ O PQ R S%MTU V 1 exp Y- A∆[ PQ \TF R %TU V 1 ] (1)
Where 𝜆 = 𝜆 ? + 𝜆 ^ is the total reorganization energy, quantifying the dynamic energetic fluctuations and including both the internal (𝜆 ? ) and environmental contribution (𝜆 _ ); ∆𝐺 ?@ is the difference between the molecular site energies (averaged over dynamic fluctuations) and eventually including the effect of an applied electric field. For electron transport, the site energies correspond to the electron affinity with reversed sign, i.e.

∆𝐺 ?@ = -𝐸𝐴 @ + 𝐸𝐴 ? .

The site energy of a given molecule includes an intramolecular and an environmental component [START_REF] D'avino | Electrostatic Phenomena in Organic Semiconductors: Fundamentals and Implications for Photovoltaics[END_REF] :

𝐸𝐴 = 𝐸𝐴 bc^-Δ (2)
The term 𝐸𝐴 bc^ corresponds to the electron affinity of the isolated molecule, and we approximate it with the LUMO en-ergy, 𝐸 efgh . The environmental component Δ is instead obtained with classical microelectrostatic (ME) calculations using MESCal code, [START_REF] D'avino | Electronic Polarization in Organic Crystals: A Comparative Study of Induced Dipoles and Intramolecular Charge Redistribution Schemes[END_REF] in which molecules interact through forces between atomic permanent charges and atomic induced dipoles, the latter requiring a self-consistent treatment of mutual interactions. The environmental contribution to site energies can be further partitioned into an electrostatic (Δ i ) and an introduction (Δ j ) term, Δ = Δ i + Δ k . In particular, Δ i corresponds to the energy to charge a molecule in the electrostatic potential of the permanent and induced multipoles of the surrounding molecules in the neutral sample, while Δ j accounts for the dielectric screening of the added charge by the polarizable medium. [START_REF] D'avino | Electrostatic Interactions Shape Molecular Organization and Electronic Structure of Organic Semiconductor Blends[END_REF] The induction term Δ j affects the absolute values of the energy levels, typically increasing the 𝐸𝐴 by ~ 1 eV, but it is essentially isotropic and then contributes only marginally to the spread of the transport levels. [START_REF] D'avino | Electrostatic Phenomena in Organic Semiconductors: Fundamentals and Implications for Photovoltaics[END_REF] Then, this term is here disregarded, since our primary interest is the estimation of energetic disorder.

Our samples are characterized by the presence of structural disorder which may strongly affect the energetics of localized charge carriers (see e.g. ref. [START_REF] Gali | Energetic Fluctuations in Amorphous Semiconducting Polymers: Impact on Charge-Carrier Mobility[END_REF] for a discussion of different effects in play). In fact, each molecule experiences a different electrostatic environment, and this leads to a broad distribution of charge transport levels. These distributions often feature an approximately Gaussian shape, whose standard deviation 𝜎 is taken as a measure of the energetic disorder. The broadening of transport levels has been estimated upon sampling the 𝐸𝐴 over all the molecules, and for 11 thin film MD configurations. This allows us also to compute the static (or positional) and dynamic nature in of the energetic disorder. To such an aim, we define the average 𝑥̅ o and the variance 𝜎 o ( of the energy of molecule 𝑚 over different configurations. Static and dynamic disorder are then computed as 𝜎 ^pcp?qr S〈(𝑥̅ o -〈𝑥〉) ( 〉 and 𝜎 uvwco?qr S〈𝜎 o ( 〉, respectively, where angle brackets denote the average over the molecules in the sample. It can be shown that 𝜎 pxpcy ( = 𝜎 ^pcp?q ( + 𝜎 uvwco?q ( . For the parameterization of the ME method, molecular polarizability and ESP charges for the neutral molecule and the anion have been obtained at the equilibrium geometry by PBE0/cc-pVTZ DFT calculations using Gaussian16. 50 The simulated samples were treated as thin films, applying periodic boundary conditions in two dimensions (x, y) and constructing supercells of the original samples. We have verified that a cutoff of 44 nm (about 8 times the average thickness of the samples) ensures a proper convergence of electrostatic interactions, i.e. within 50 meV for absolute energies and 5 meV for energy differences between sites. Electronic couplings were calculated for every pair of molecules ("dimers", having at least two atoms at distance lower than 7 Å) in a single configuration of the B2 and B4 samples. Calculations were carried out at DFT level B3LYP/6-31G(d), ORCA code [START_REF] Neese | The ORCA Program System[END_REF] with the dimer projection method. [START_REF] Senthilkumar | Charge Transport in Columnar Stacked Triphenylenes: Effects of Conformational Fluctuations on Charge Transfer Integrals and Site Energies[END_REF] The intermolecular reorganization energies 𝜆 ? for B2 and B4 were calculated using the four-point method [START_REF] López-Estrada | Reassessment of the Four-Point Approach to the Electron-Transfer Marcus-Hush Theory[END_REF] via PBE0/cc-pVTZ calculations , and amount to 0.3 eV for B2 and 0.32 for B4, respectively. These values are considerably larger than the one of the prototypical n-type organic semiconductor C60 fullerene. [START_REF] Tummala | Static and Dynamic Energetic Disorders in the C60, PC61BM, C70, and PC71BM Fullerenes[END_REF][START_REF] Faber | Electron-Phonon Coupling in the C60 Fullerene within the Many-Body GW Approach[END_REF]

RESULTS

3.1 B4 and B2 Films. Since the experimental fabrication of B4 and B2 films is made by vacuum deposition, to mimic as much as possible experimental conditions MD simulations of the bulk/vacuum interface have been performed. The intent is to include the effect of the vacuum/film interface affecting molecular geometries and structural properties. To this aim, four independent configurations of B2 and B4 have been extracted from MD trajectories of bulk phase (a .configuration every 20 ns of the last 80 ns of each trajectory) and used as a starting coordinate set to simulate vacuum/bulk interface systems. The vacuum interface has been built by increasing the box side (z) by 20 nm. Then, MD simulations in NVT ensemble have been performed for each independent configuration. Representative snapshots of equilibrium configurations of B2 and B4 systems are reported in Figure 2. The visual inspection of equilibrium configurations reveals a well-defined columnar structure for the B4 system (left side of Figure 2). Instead, B2 molecules (right side of Figure 2) are much more disorderly stacked on each other, and, many of them appear tilted with the respect to the normal direction (z) of the bulk/vacuum interface. A more quantitative indication is gained by computing the X-Ray powder pattern and the radial distribution function (g(r)). As can be seen from panel A in Figure 3, both X-Ray traces show a peak at q~ 1.53 Å -1 , corresponding to molecular 𝜋-stacking distance. The B4 profile shows a sharp peak at this position, which is comparable with the experimental spacing of ~ 4 Å measured on a thin film of B4 material. [START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] The comparison with the broader peak found for B2 indicates a decrease of positional order along the stacking direction for the latter system. The difference in the positional order between B2 and B4 systems is confirmed by the g(r) displayed in Figure 3B. [START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] Using the last 80 ns of trajectories, structure factor patterns for B4 (black) and B2 (red) were calculated by using the ISAACS package. [START_REF] Le Roux | ISAACS Interactive Structure Analysis of Amorphous and Crystalline Systems[END_REF] (B) Radial distribution functions g(r) calculated for the core hydrogens. Explicit assignment of the main peaks is reported in the plot. A snapshot of B4 system, counting eight molecules belonging to adjacent columns is reported to outline the hexagonal packing.

These functions have has been computed by considering only the hydrogen atom belonging to the central pyrimidine ring; those atoms are shown as red spheres in the snapshot of Figure 3B. Due to its peculiar position, the pyrimidine hydrogen is a convenient choice to characterize structural properties at molecular (π-stacking) and bulk scales (hexagonal packing in case of B4). As can be seen from the comparison between B2 and B4 curves, it is clear that B4 has much higher positional order, in particular in the region of π-stacking (~ 3-5 Å) which is fundamental for promoting electron transport, with peaks about three times more intense than for B2. At larger distances, the B4 g(r) shows two main peaks, region I at ~ 1.5 nm and region II at ~ 2.8 nm, which represent the characteristic interplanar distances of the nearly hexagonal packing. Instead, B2 shows a featureless broad halo in the same region, which indicates lower structuration. Additional snapshots highlighting the hexagonal packing are reported in Figure S10 of Supporting Information. Since structural order and molecular orientation are correlated, a high positional order of B4 structure should correspond to a high orientational order, too. The regularity in molecular orientation can be computed using the orientational order parameter (P2), which is defined as:

〈𝑃 ( 〉 = €〈qx^R•〉$* ( ( 4
)
where 𝜃 is the angle between an arbitrary molecular axis m (see definition in Figure 4A) and a reference direction z. The angle bracket indicates the ensemble average. If we choose the molecular axis in the pyrimidine plane as in reference [6], when every molecular axis is oriented parallel to the z vector (the normal to the surface), and hence in a perfect edge-on configuration, P2 is equal to 1, for an isotropic distribution of orientations P2 equals to 0 and for a perfect face-on distribution, P2 is equal to -0.5.

In Figure 4B, computed and experimental P2 are compared: while for B2 the orientational order is slightly overestimated with respect to the measured value, the calculated values show the same trend of the experimental data, with B4 being more ordered that B2. More in detail, we observe that B4 molecules have a more pronounced "face-on" orientation than B2. The tendency to form well-ordered columnar stacks is beneficial for key properties such as charge transport and efficiency in OLEDs devices, as reported by Kido and co-workers. [START_REF] Ohisa | Molecular Interdiffusion between Stacked Layers by Solution and Thermal Annealing Processes in Organic Light Emitting Devices[END_REF][START_REF] Watanabe | Simultaneous Manipulation of Intramolecular and Intermolecular Hydrogen Bonds in N-Type Organic Semiconductor Layers: Realization of Horizontal Orientation in OLEDs[END_REF][START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF][START_REF] Watanabe | Control of Molecular Orientation in Organic Semiconductor Films Using Weak Hydrogen Bonds[END_REF][START_REF] Yokoyama | Molecular Orientation in Small-Molecule Organic Light-Emitting Diodes[END_REF] 3.2 Hydrogen Bonds. The higher solubility of B4 in chloroform with a small amount of methanol is an evident experimental proof of the strongest H-bond network formed by B4 with respect B2. [START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] In addition, the presence of a strong bidimensional H-Bond network was also demonstrated by scanning tunneling microscopy (STM) experiments carried out at the interface between a solution of B4 and a graphite surface. [START_REF] Ziener | Supramolecular Assemblies of a Bis(Terpyridine) Ligand and of Its [2×2] Grid-Type ZnII and CoII Complexes on Highly Ordered Pyrolytic Graphite[END_REF][START_REF] Meier | Weak Hydrogen Bonds as a Structural Motif for Two-Dimensional Assemblies of Oligopyridines on Highly Oriented Pyrolytic Graphite: An STM Investigation[END_REF][START_REF] Ziener | Self-Assembled Nanostructures of Oligopyridine Molecules[END_REF][START_REF] Meier | Adsorbate-Substrate-Mediated Growth of Oligopyridine Monolayers at the Solid/Liquid Interface[END_REF][START_REF] Meier | Concentration and Coverage Dependent Adlayer Structures: From Two-Dimensional Networks to Rotation in a Bearing[END_REF] To detect H-bonds features, simulation trajectories have been analyzed in detail. In particular, B2 and B4 molecular structures display several possible (weak) H-bond acceptor (A) and donor (D) sites (definitions in Figure 5A). During the simulation, a hydrogen bond is counted when two geometrical conditions are simultaneously satisfied: [START_REF] Van Der Spoel | Thermodynamics of Hydrogen Bonding in Hydrophilic and Hydrophobic Media[END_REF] First, the distance between A and D must be shorter or equal to rC cut-off distance of 0.3 nm; second, the angle 𝜃 (ƒ$"$…) must be larger than 𝜃 q †p$x ‡ ‡ = 120º. We have tested different cut-off distances and cut-off angles in the H-Bond counting to verify the correct applicability of the used values proposed in reference [63]. Additional analyses are reported in the Supporting Information.

In agreement with the experimental evidences, 6,58-62 we have found a higher average number of H-bonds/molecule (~ 2.5) for B4 with respect to B2 (1.5), as shown in the histogram in Figure 5B. According to a geometrical criterion based on the 𝜋-stacking distance (~ 4 Å), hydrogen bonds were also labelled as inplane or out-of-plane. If the D-A distance along the normal direction of bulk/vacuum interface z) between D and A is longer than 4 Å, the H-Bond is counted as out-of-plane; for longer distances the H-Bond is counted as out-plane, or in other words the bond is roughly coplanar with the molecular planes. The fraction of in-plane bonds is reported in Figure 5, panel C, alongside the intra molecular H-bond fraction in panel D. In case of B4, the majority of H-bonds is formed between molecules laying a same plane (~ 0.87), instead, a lower fraction (~ 0.58) is found for B2. Furthermore, B4 molecules do not form, in practice, intramolecular H-bonds, while B2 has a noticeable amount of them (nearly 0.12, Figure 5D). These results, in addition to previous structural results (including order parameter, X-Ray pattern, g(r) and mass density) confirm the interpretation of the assembly given by Kido and co-workers [START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] and also discussed in the interpretation of B4 two dimensional structures. [START_REF] Ziener | Supramolecular Assemblies of a Bis(Terpyridine) Ligand and of Its [2×2] Grid-Type ZnII and CoII Complexes on Highly Ordered Pyrolytic Graphite[END_REF][START_REF] Meier | Weak Hydrogen Bonds as a Structural Motif for Two-Dimensional Assemblies of Oligopyridines on Highly Oriented Pyrolytic Graphite: An STM Investigation[END_REF][START_REF] Ziener | Self-Assembled Nanostructures of Oligopyridine Molecules[END_REF][START_REF] Meier | Adsorbate-Substrate-Mediated Growth of Oligopyridine Monolayers at the Solid/Liquid Interface[END_REF][START_REF] Meier | Concentration and Coverage Dependent Adlayer Structures: From Two-Dimensional Networks to Rotation in a Bearing[END_REF] In particular, B4 pyridine nitrogen atoms are positioned at the very edge of the molecule, and then favoring in-plane connections between pairs of coplanar molecules via intermolecular H-bonds, and to form a bidimensional network. The ortho position of nitrogen atoms in pyridine rings in B2 instead implies a rotation of those rings to maximize the number of intermolecular H-Bonds. However, even if a 2D H-bond is formed, such a rotation of reduces the planarity and, as a consequence, the packing ability of B2. This is clearly reflected by the lower mass density reached by B2 bulk, reduced structural order and lower number of H-Bonds per molecule with respect to B4. In addition to those evidences, the distributions of φ .* and φ .( dihedral angles suggest a slightly higher planarity for B4 (Figure S11 of Supporting Information).

Electron Transfer Parameters.

As already mentioned, Kido and coauthors 6 measured the electron mobility of B2 and B4 with the time of flight technique, obtaining at room temperature values of the order of 10 $, and 10 $% cm 2 V -1 s -1 , respectively. In this section, exploiting the atomic resolution of the system morphologies provided by MD simulations, we aim at establishing a link between structural and electronic properties, and analyze the factors that contribute to the large offset between B2 and B4 mobility.

We begin the comparison of the charge transport properties of the two systems by inspecting the distributions of electron affinities, i.e. the density of states available to charge carriers. The histograms in Figure 6 clearly show two main differences: first, B4 has an electron affinity larger than B2, which is a positive characteristics because it improves the air stability of the semiconductor [START_REF] Bin | Efficient N-Dopants and Their Roles in Organic Electronics[END_REF] and facilitates a possible doping by electron donating impurities. Moreover, the distribution of electron affinity for B4 is narrower than the one for B2. Again, this goes in favor of superior electron transport in the B4 system, since its standard deviation 𝜎 is usually taken as a measure of the energetic disorder, which in turn is one of the main limiting factors for charge mobility. [START_REF] Gali | Energetic Fluctuations in Amorphous Semiconducting Polymers: Impact on Charge-Carrier Mobility[END_REF] For a quantitative analysis of the energetic disorder and its nature, we then analyzed the standard deviations of 𝐸𝐴, 𝐸 efgh , and ∆ Š , reported in Table 2. It is also useful to partition the disorder in a static (or positional) term, that does not change in time and affects the ∆𝐺 ?@ values in eq. 1, and a dynamic one which instead has to be included in the reorganization energy. [START_REF] Tummala | Static and Dynamic Energetic Disorders in the C60, PC61BM, C70, and PC71BM Fullerenes[END_REF] Both static and dynamic terms negatively impact charge transfer rates, ultimately reducing the mobility, although each with a distinct physical mechanism. We note that the calculated static disorder of 𝐸𝐴 is larger for B2, and qualitatively in line with the experimental values, which were obtained by fitting mobility dependence on electric field and temperature with an empirical but well established equation for amorphous semiconductors [START_REF] Bässler | Charge Transport in Disordered Organic Photoconductors a Monte Carlo Simulation Study[END_REF] (0.091 and 0.076 eV for B2 and B4, respectively [START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] ). This observation is in line with the higher structural order present in B4 samples. The amount of dynamic disorder is instead rather similar for the two systems, yet again slightly larger for B2. These room-temperature energetic fluctuations correspond to a reorganization energy

𝜆 = ‹ OE•Ž••P' R U V 1
about 0.12 eV higher for B2. From Table 2 it can be further seen that the difference in energetic disorder between the two systems, both for the static and dynamic component, is sourced from the environmental electrostatic contribution ∆ Š , while the fluctuations of the LUMO energies are quite similar, albeit again a bit larger for B2.

To explain the higher disorder of the energy levels in B2 and its chiefly dynamic nature, let us consider the four terminal pyridine units: any thermally activated rotation about the pyridinephenyl bond marginally affects the electrostatic layout of B4, where the local dipole associated with the pyridine nitrogen is coaxial with the bond. On the other hand, pyridine rotations can produce a rather large variation of charge distribution for B2, where the dipole is off-axis. This makes that, for B2, a given amount of dynamic conformational disorder translates in more important variations of the electrostatic landscape felt by the charge carriers. The argument is strengthened by the DFT calculation of the gas phase molecular dipoles along the trajectories, that yields larger values and fluctuations for B2: 2.9 ± 1.3 D versus 0.8 ± 0.3 D for B4. We now focus our attention on the magnitude and fluctuations of intermolecular electron-transfer couplings 𝐽, representing the pre-exponential factor in equation 1, and on their relationship with the structural organization in the simulated samples. Figure 7A, compares the magnitude and directionality of intermolecular couplings between the two materials, reporting the absolute value of 𝐽, computed for molecular pairs in our MD samples, against 𝑟 -/𝑟 (gray dots), being 𝑟 the module of the intermolecular distance and 𝑟 -its z-axis component. Charge transfer couplings are highly anisotropic in both systems (see also Figure S16), attaining appreciable values only when the intermolecular vector is parallel to the stacking direction (𝑟 -/𝑟 > 0.6). As better appreciated by the average values of |𝐽| over 𝑟 -/𝑟 intervals (red squares), the magnitude of the couplings is significantly larger in B4 than in B2, and this difference determines a superior electronic connectivity with respect to B2. This quantity, shown in Figure 7B, corresponds to number of neighbors that, on average, is connected to each molecule by a coupling exceeding a given value of |𝐽|. The two systems approach in a very different way the connectivity of 2 characteristic of an ideal 1D semiconductor, with e.g. B4 attaining the value of 1.5 with couplings of 41 meV, and B2 only for |𝐽|=12 meV.

The superior charge transfer properties of B4 as compared to B2 are completely unexpected at the molecular level, given the similar structure and the nearly identical LUMO orbitals (see inset Figure 7A) of the two molecules. [START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] We expect that the difference between the two materials may be due to the better supramolecular packing of B4 discussed above. To prove this hypothesis with an "alchemical" experiment, impossible to perform with real systems, we have computed the intermolecular couplings for a fictitious B4 phase, obtained from the atomic coordinates of the B2 sample by switching the positions of the N atom and one CH group in the peripheral pyridine rings. Figure 7C, compares the couplings of this hypothetical sample, labelled "B4 < B2", with the original B2 ones. This plot reveals a nearly perfect correlation between the two data sets (Pearson coefficient R=0.994), demonstrating that the differences between the charge transport properties of the two materials does not arise from the chemical structure. It has instead to be ascribed to the better molecular organization in B4, in turn resulting from the presence of a two-dimensional H-bond network.

CONCLUSIONS

Combining DFT calculations and MD all-atom simulations we investigated the role of intermolecular hydrogen bonds in the molecular orientation, stacking and electron carrier mobility of thin films composed of pyridine-based molecules B2 and B4.

The only difference between the B4 and B2 chemical structure lies in the position of nitrogen atoms in the four outer pyridine rings, which is surprisingly sufficient for causing a huge difference in the optical and electron transport properties, as experimentally proved. We found that B4 has a higher capability to connect with adjacent molecules, by intermolecular hydrogen bonds, compared with B2 (2.5 vs 1.5 bonds per molecule). The stronger 2D H-bond network formed for B4 is sufficient to stabilize a supramolecular organization in which molecules form ordered layers and columns.

The calculation of essential electron transport parameters, such as reorganization energies, site energies, and electronic couplings revealed that the microscopic origin of the higher electron mobility of B4 with respect to B2 has its roots in lower reorganization energy and electrostatic disorder and, more importantly, in the superior morphology promoted by the hydrogen bond network.

On a general vein, our findings disclose how intermolecular hydrogen bond interactions promote structural order in molecular thin films, ultimately improving their electron transport properties. The engineering of these weak interactions offers viable opportunities for the rational design of new materials with improved performances in OLED applications.

These findings provide a detailed molecular view of the Hbond interactions governing spatial organization and charge transport in thin films. The knowledge gained from the study is of primary importance in developing new materials and/or to improve performances of OLED devices.

ASSOCIATED CONTENT Supporting Information

Model description, Force-Field parameter tables, parameterization procedure for soft aryl-aryl torsions. Additional analysis of crystalline system of B4PyPTZ, procedure description for the thin film preparation of B4 and B2. Additional figures and analysis for Hbond and charge transport. The Supporting Information is available free of charge on the ACS Publications website.

Notes

Any additional relevant notes should be placed here.

Figure 1 .

 1 Figure 1. Chemical structure of B4 and B2 molecules. The different position of nitrogen atoms in the 4-and 2-pyridyl rings are indicated. In the following sections of the paper the dihedral angles φ .* and φ .( are also indicated.

Figure 2 .

 2 Figure 2. Snapshots representative of the equilibrium configuration of B4 and B2 systems. Nitrogen atoms are represented in blue, carbons in cyan and hydrogens in white. The vacuum/bulk interface is perpendicular to z direction.

Figure 3 .

 3 Figure 3. (A) X-Ray powder pattern is calculated using 𝐶𝑢𝐾 ~ radiation (𝜆 = 1.5418 Å) as in the experiments.[START_REF] Yokoyama | Molecular Stacking Induced by Intermolecular C-H•••N Hydrogen Bonds Leading to High Carrier Mobility in Vacuum-Deposited Organic Films[END_REF] Using the last 80 ns of trajectories, structure factor patterns for B4 (black) and B2 (red) were calculated by using the ISAACS package.[START_REF] Le Roux | ISAACS Interactive Structure Analysis of Amorphous and Crystalline Systems[END_REF] (B) Radial distribution functions g(r) calculated for the core hydrogens. Explicit assignment of the main peaks is reported in the plot. A snapshot of B4 system, counting eight molecules belonging to adjacent columns is reported to outline the hexagonal packing.

Figure 4 .

 4 Figure 4. (A) scheme showing the definition of molecular axis m (red arrow) and face-on, edge-on orientation on a substrate. (B) P2 calculated from MD simulations (solid triangles) and from experiments 6 (open squares).

Figure 5 .

 5 Figure 5. (A) definition of hydrogen bond acceptor and donor sites for B4. (B) average number of H-bonds/molecule. (C) In-plane H-Bond fraction, and (D) intra molecular H-bonds fraction. Data has been computed on the last 35 ns of MD trajectories.

Figure 6 .

 6 Figure 6. Distribution of electron affinities for B4 and B2, averaged on eleven MD configurations. Individual values were obtained as 𝐸𝐴 = -(𝐸 efgh +∆ Š ) by DFT and microelectrostatic calculations.

Figure 7 .

 7 Figure 7. (A) Electron-transfer couplings as a function of the normalized projection of the intermolecular vector distance 𝒓 along the stacking axis z. Gray dots show the couplings computed for molecular pairs extracted from the MD sample, red squared (error bars) show the mean value (standard deviation) over intervals of 𝑟 𝑧 /𝑟. LUMO orbitals of B4 and B2 are shown as insets. (B) Electronic connectivity, i.e. average number of neighbors per molecule having an electronic coupling larger than the abscissa value. (C) Correlation between the couplings of B2 and those of a fictitious B4 sample obtained modifying the atomistic morphologies of B2 (see text). The quadrant bisector (black line) is shown as a guide for the eye. The same plot in a log-log scale in shown as an inset.

ABBREVIATIONS

  ABBREVIATIONSB4 = (4,6-bis(3,5-di(pyridine-4-yl)phenyl)-2-methylpyrimidine B2 = (4,6-bis(3,5-di(pyridine-2-yl)phenyl)-2-methylpyrimidine

Table 1 . Characteristics of the simulated systems.

 1

					Aver-
	System	No. Molecules	Molecule Type	En-semble	age Box Length [x,y,z]	Time (ns)
					(nm)
					5.382
					± 0.1,
	Bulk	128	B4PyPTZ	NPT	3.514 ± 0.08,	100
					5.443
					± 0.1
					4.698
					± 0.1,
	Bulk	128	B4	NPT	3.596 ± 0.1,	100
					5.701
					± 0.12
					5.102
					± 0.1,
	Bulk	128	B2	NPT	3.587 ± 0.1,	100
					5.567
					± 0.17
					4.751,
	Film	128	B4	NVT	3.627,	350
					25.706
					5.069,
	Film	128	B2	NVT	3.561,	350
					25.566
	2.

Table 2 :

 2 Averages of intramolecular (𝑬 𝐋𝐔𝐌𝐎 ), intermolecular (∆ 𝑬 ), and total electron affinities (𝑬𝑨 = -𝑬 𝐋𝐔𝐌𝐎 -∆ 𝑬 ), and corresponding standard deviations separated in static and dynamic contributions 𝝈 𝒔𝒕𝒂𝒕𝒊𝒄 and 𝝈 𝒅𝒚𝒏𝒂𝒎𝒊𝒄 . Average energies are in eV, standard deviations in meV.

	B2	average	𝝈 𝒕𝒐𝒕𝒂𝒍	𝝈 𝒔𝒕𝒂𝒕𝒊𝒄	𝝈 𝒅𝒚𝒏𝒂𝒎𝒊𝒄
	𝐸 ¨©ª«	-1.681	97	50	83
	∆ Š	0.000	181	124	132
	𝐸𝐴	1.681	183	111	145
	B4	average	𝝈 𝒕𝒐𝒕𝒂𝒍	𝝈 𝒔𝒕𝒂𝒕𝒊𝒄	𝝈 𝒅𝒚𝒏𝒂𝒎𝒊𝒄
	𝐸 ¨©ª«	-2.201	81	30	75
	∆ Š	0.156	129	67	110
	𝐸𝐴	2.046	156	81	134

ACKNOWLEDGMENT

GM and ADN acknowledge the computing resources and the related technical support used for this work have been provided by CRESCO/ ENEAGRID High Performance Computing infrastructure and its staff. CRESCO/ENEAGRID High Performance Computing infrastructure is funded by ENEA, the Italian National Agency for New Technologies, Energy and Sustainable Economic Development and by Italian and European research programs, see ww.cresco.enea.it/English/ for information. LM acknowledges funding through PRIN project "From natural to artificial light-harvesting systems: unveiling fundamental processes towards a bio-inspired materials design" (HARVEST) protocol 201795SBA3.

Funding Sources

AUTHOR INFORMATION Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.