A quantitatively consistent, scale-spanning model for same-material tribocharging

We propose a quantitative, scale-spanning model for same-material tribocharging. Our key insight is to account for mesoscale spatial correlations in donor/acceptor surface properties, which dramatically affect the macroscopic charge transfer and quantitatively reconcile previous inconsistencies related to the microscale. We furthermore identify a viable mechanism by which the mesoscale features emerge, which may help constrain the list of donor/acceptor candidates. As the only free-parameters in our model involve the atomic scale, data analyzed in light of it could help resolve the detailed mechanism of tribocharging.

Tribocharging, i.e. charge transfer between materials during contact [1], plays a critical role in natural phenomena [2][3][4][5][6], industrial processes [7,8], and energy harvesting devices [9][10][11], yet resists interpretation. One fundamental roadblock has been the inability to identify the atomic-scale mechanism, and in particular the charge carriers, i.e. ions vs. electrons [1]. An equally important roadblock is the lack of quantitative agreement (or even comparison) between experiments and theory. Experimentally, issues such as the difficulty of measuring contact areas render much data qualitative (e.g. the sign or scale of charging) [12][13][14][15][16][17]. Theoretical advances are stymied by the multi-scale nature of the effect, where one must simultaneously account for probabilistic effects at the atomic scale (<1 nm), unexplained emergent features at the mesoscale (∼1 μm), and then through these explain the familiar behavior of the macroscale (>1 mm). These challenges have perhaps led to a general disinterest in quantitative reconciliation, with some authors characterizing the outlook on achieving a scale-spanning description as 'impossible' [18].
Same-material tribocharging, where charge is exchanged between identical materials, is perhaps the most puzzling manifestation of the phenomenon. It has been attributed to trapped electrons [19][20][21][22][23][24][25], induced polarization [26][27][28][29][30], or mechanochemistry [7,[31][32][33]]. Yet, experiments on same-material tribocharging have produced some valuable clues. Using soft (Young's modulus ∼1 MPa), atomically smooth (roughness < 1 nm) polymers to achieve conformal contact, Apodaca et al. found that the magnitude of charge transfer grows with the square root of the contacting area [34]. They proposed that the surface consists of equally-sized, randomlyassigned donor/acceptor sites (Fig. 1a), each capable of giving/receiving one unit charge. This allowed them to recover |∆Q| = C √ A, where the prefactor, C depends on the length scale of a single site, l 0 . Though l 0 should be on the scale of one atom, fits to their model required a value of ∼0.005Å-more than 100 times smaller than the Bohr radius of hydrogen. Thus, while their idea is compelling, this discrepancy suggested it is still somehow a 'toy model.' suggested that same-material tribocharging may arise from equally-sized, randomlydistributed donor/acceptor sites, where unit charges, e, are transferred during contact [34]. (b) Surfaces with identical microscopic size L, microscopic scale l0, but different mesoscale correlation lengths, l. On the left, l = l0, and on the right l = 5 l0 (where l0 is the site size). (c) The single-contact distribution in ∆Q is broader when l = 5 l0 (violet, ∼200 e) than when l = l0 (pink, ∼100 e). (d) For sequential contacts, charge transfer is enhanced when l = 5 l0, and fluctuations that dominate the l = l0 case are suppressed.
Nonetheless, features suggestive of their idea have been observed, but at significantly larger scales. Us-ing Kelvin Force Probe Microscopy on the same polymers, Baytekin et al. found 'mosaics'-regions of positive/negative charge after contact-that were spatially correlated over length scales up to ∼450 nm [35]. This observation raises fundamental questions: How can the impossibly small scales implied by Apodaca be connected to the mesoscopic ones seen by Baytekin? What would the introduction of an intermediate scale change in Apodaca's model? And why do the intermediate scale correlations emerge in the first place? In this work, we use analysis supported by numerical simulations to investigate the multi-scale aspects of same-material tribocharging. Our key development is to properly account for the mesoscale spatial correlations, which allows us to quantitatively resolve the discrepancies in the Apodaca framework and produce a scale-spanning model.
We start by explaining a first set of numerical simulations, where we mimic charge transfer between 'synthetic' surfaces (in contrast to physically derived surfaces later) by creating two N -element matrices involving three length scales: l 0 , l, and L. The smallest, l 0 , corresponds to the elementary donors/acceptors of the atomic scale, and is represented by a single matrix element. The largest, L, corresponds to the macroscopic system size. We assume that there is a single intermediate scale, l, that characterizes the mesoscopic correlations observed by Baytekin. (Note they measured charge, but this implies donor/acceptor correlations.) Each matrix element is assigned as donor or acceptor, with probabilities p and 1 − p, respectively. We account for correlations in assignments via thresholding a random scalar field (see Supplemental Material [36]). To perform a 'contact', we first generate a 'left' and 'right' surface from identical input length scales and probabilities (Fig. 1a). Charge transfer of one unit, e, between matrix elements [i,j ] occurs if (1) [i,j ] on the left/right is a donor, (2) [i,j ] on the right/left is an acceptor, (3) the value of an independent random uniform variable is less than the transfer probability, α, and (4) for sequential contacts, transfer at [i,j ] hasn't yet occurred. The net charge transferred is the difference between left-to-right ('right') and right-to-left ('left') transfers. Figure 1b shows two representative surfaces, one with l = l 0 and the second l = 5 l 0 . Although their only difference involves the length scale, l, we see stark changes in the charging behavior. In Fig. 1c, we plot distributions of the charge transferred in the first contact, ∆Q, for 1000 pair-instances. As we do not assume any post-contact discharge [37], both distributions are Gaussian and centered at zero, but while l = l 0 produces a width of ∼100 e, the l = 5 l 0 width is ∼200 e. Correlations also change the behavior during sequential contacts. Fig. 1d shows two examples of the accumulated charge, |Q|, vs. the number of contacts, n c . Like the initial transfer, the final charge, Q f , is typically larger when l = 5 l 0 . Additionally, the fluctuations are on the order of Q f when l = l 0 , When l ≈ l0, transfer is probabilistic (Eq. 1) and |∆Q| ∝ √ α.
but are hardly discernible when the l = 5 l 0 (consistent with experiments [34]).
We now examine the first contact behavior of these simulations in detail. Figure 2a shows the average absolute value of charge exchanged, |∆Q|, for increasing system size and several values of l. The scaling |∆Q| ∝ √ A ∝ L/l 0 is recovered in all cases, but the prefactor steadily increases with l > l 0 . Figure 2b shows that the dependence on the transfer probability, α, exhibits unexpected non-linear behavior as a function of the correlation length, l. At every value of l, we see a trend consistent with a power law, i.e. |∆Q| ∝ α k . However, the exponent, k, increases with l, starting at k = 0.5 for l = l 0 and saturating at k = 1.0 for l l 0 (Fig. 2d).
To analytically explain these observations, we first consider the case L l = l 0 , here sketching our argument (details are in the Supplemental Material [36]). We momentarily focus on right transfers, which occur with compound probability p(1 − p)α. Absent correlations, all sites [i, j] are independent, hence the total right charge transfer is Gaussian with mean eN p(1 − p)α and width e N p(1 − p)α(1 − p(1 − p)α). A similar distribution exists for left transfer, but technically only when considered independently-simultaneous left/right transfer cannot occur. Nonetheless, the probability for this is small, and we therefore approximate the left/right distributions as independent. The net transfer, ∆Q, is thus also Gaussian distributed, with zero mean and width This recovers the √ A scaling, but points to a slight mistake in the earlier work [34] in that the α-dependence is square root rather than linear-exactly as our simulations in Fig. 2d. In the Supplemental Material [36], we verify that Eq. 1 collapses our simulated data for wide ranges of p and α.
Next we consider the case L l l 0 (again with details in the Supplemental Material [36]). This fundamentally alters the argument above as the site identities exhibit spatial correlations. To handle this, we first imagine rescaling the system by l/l 0 , leading to surfaces with N = N/(l/l 0 ) 2 larger 'patches', each consisting of many sites. Identities of entire patches still occur with probabilities p and 1 − p. Next, we rescale back to deal with transfers, which still occur independently for each site. During contact, regions of donors face acceptors with the characteristic size of a patch. If the number of sites in these regions (n = (l/l 0 ) 2 ) is large, the mean transfer per patch (αn) effectively hides the fluctuations ( nα(1 − α))-hence we treat α as a rate. Thus, when l l 0 we have Here, like in the Apodaca work, the α-dependence is linear. The critical difference, however, is the dependence on the intermediate length scale, which amplifies the charge transfer by the factor, l/l 0 . We confirm this with our simulated data in Fig. 2c, where the prefactor F = αl/l 0 4p(1 − p)/π collapses |∆Q| when l l 0 . Qualitatively, the explanation for this amplification is that the variability (i.e. standard deviation) in the number of donors/acceptors on a surface increases with the scale of spatial correlations. One can quickly grasp why by considering the extreme case l = L, where each surface is either purely donor or acceptor, and consequently |∆Q| ∝ αeN ∝ A. This highlights that sameand different-material tribocharging are two manifestations of a similar underlying phenomenon, only appearing different depending on the scale at which one looks. We now turn to sequential contacts. As the simulation results of Fig. 3a show, repeated contacts with the same surfaces leads to curves in accumulated charge, |Q| vs. n c , that level off at some value, Q f . In the Supplemental Material, we show that the underlying trend is given by a saturated exponential, i.e., (3) Figure 3a also illustrates the presence of fluctuations on approach to Q f . For large l/l 0 , these aren't noticeable, but for small l/l 0 they are overwhelming, and cannot be suppressed even when we increase the system size, L. We quantify their scale by repeatedly performing first contacts between individual surface pairs ('resetting' each time) and measuring the standard deviation, δQ (Fig. 3b). In Fig. 3c we show δQ for a few pairs, which reveals that the fluctuations grow with the macroscale length, L, and depend strongly on α, but are independent of the mesoscale, l. To analyze why, we note that a particular pair has a fixed donor/acceptor arrangement, which means δQ arises solely from α. Denoting the number of donors on the left/right that face acceptors on the right/left as N , one finds that δQ = e (N + N )α(1 − α). In the Supplemental Material [36], we justify using the averages N = (L/l 0 ) 2 p(1 − p) to find the ensemble expression, This establishes that the fluctuations, like |∆Q|, grow linearly with L, but unlike |∆Q| are independent of l. Consequently, they cannot be suppressed by increasing system size, but can be suppressed by the introduction of the intermediate scale, l. In Fig. 3d, we collapse the simulated δQ data to our predicted line with the appropriate rescaling.
The last question we posed remains: Why do the intermediate scale correlations emerge in the first place? A few mechanisms have been proposed. For example, with inelastically deformed materials it has been shown that micron-sized voids form, which are somewhat suggestive of mosaics [38]. Other authors make connections to 'islands' of adsorbed surface water, but the sizes of these has not been measured nor compared to the charge mosaics [16,17,37,39,40]. In the case of void formation, the process for a growing lenghtscale is clear-the materials are ripped apart-yet the polymers used by Apodaca/Baytekin are highly elastic and not intentionally stretched. In the case of water islands, and indeed virtually all 'patch'-type models, one must assume features of a certain size exist, but how that length scale emerges is still unclear.
We propose that a nucleation process which energetically favors neighboring donors (or equivalently, neighboring acceptors) is a viable candidate. We support this by developing a second, distinct set of simulations to mimic the physics of surface formation, which are timedependent and latticed-based, and which again produce L/l 0 × L/l 0 donor/acceptor matrices. At each time step, we assume a donor site can transition into an acceptor site, and vice-versa. This could represent widely suspected mechanisms such as adsorption of an ambient donor species (i.e. H 2 O), but also novel ideas such as phase separation during polymer curing. The transition probabilities of a given site depend on its neighbors i.e., where ν is the number of neighbors that are donors (i.e. ν ∈ [0, 4]). The exponential form is motivated from an Arrhenius-like process where each neighbor modifies a local energy barrier by /kT = K. Such as process could (a) A generic nucleation process can lead to donor/acceptor regions with a characteristic size l > l0. (b) We measure l/l0 using the radial correlation function C(r) (Eq. 6), averaged over several surfaces (error bars represent the standard deviation). The inset shows the whole vertical range. (c) We vary l/l0 through parameter K in Eq. 5. The error is calculated from the error on C(r). We keep the same color scale for l/l0 in all subfigures. The inset shows that p remains constant. (d) As before, the charge transferred after a contact is amplified by the introduction of l.  Fig. 4a . To characterize these surfaces in the context of our analysis, we measure p and l/l 0 . Determining p is trivial. To get l/l 0 we start by calculating the correlation function, where s is the site identity at a position R, r is the dis-tance from the point R, and averages denoted by are over all sites separated by r. Figure 4b shows an example of C(r) for a particular surface. We use the first zero crossing (within the standard deviation of C(r)) to define the correlation length, l, as it corresponds to the typical distance before having an equal probability of switching from donor to acceptor (or vice versa). In the Supplemental Material, we show that this identically recovers the input correlation lengths in of our other set of simulations [36]. In Figure 4c, we show how sweeping through the parameter K (P 0 fixed) allows us to explore nearly two orders of magnitude in l/l 0 (with p ≈ 0.5).
After they have reached a dynamic equilibrium, we freeze these surfaces and then use them as input for contact experiments, with the same transfer rules as before. We generate 20 surfaces for various combinations of correlation length l/l 0 and system size L/l 0 , and calculate |∆Q| for every permutation (Fig. 4d). The effect of spatial correlations in this physically-derived system is the same as with the synthetic surfaces, i.e. the magnitude of transfer increases with the correlation length. Figure 4e presents the results rescaled using Eq. 2, which largely collapses them onto a single line with the predicted unity slope. In the Supplemental Material [36], we show that the deviation for intermediate K is due a broader 'spectrum' of length scales for l than in our synthetic simulations.
Based on widely held assumptions and consideration of existing experiments, we have developed a quantitative, scale-spanning model for same-material tribocharging. We have shown that the intermediate scale corresponding to donor/acceptor spatial correlations, l, plays a crucial role, amplifying the amount of charge transferred in a single contact, and suppressing otherwise overwhelming fluctuations in sequential contacts. We have furthermore introduced a candidate mechanism for how the intermediate length scale emerges, which is based on the energetics of donor/acceptor interactions during the formation process. Our results allow us to quantitatively resolve the inconsistency encountered by Apodaca [34]: using our analysis with their data and Baytekin's correlation length [35], we infer an elementary site size l 0 ≈ 4Åprecisely on the atomic scale [36]. Furthermore, our work has implications regarding the carrier and mechanism. By elevating the Apodaca framework from a 'toy model' to a 'quantitative model,' we propose that confident extractions of α from experimental data can be made, giving information directly related to underlying atomic-scale mechanism [44]. Although our model is restricted to surfaces where the macroscopic contact area is known, these considerations become even more important in systems where roughness or stiffness play a role.
This project has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie Grant Agreement No. 754411. supplemental information Galien Grosjean, * Sebastian Wald, Juan Carlos Sobarzo, and Scott Waitukaitis

SYNTHETIC SURFACES
As discussed in the main text, our numerical simulations used in Figs. 1-3 of the main paper involve making pairs of 'synthetic' surfaces, i.e. matrices, whose elements are assigned as either donor (1) or acceptor (0). We call these 'synthetic' as the process by which we create them is not physically derived, but rather comes from a numerical algorithm that allows us to efficiently produce random surfaces with a desired donor probability p and legnthscale l/l 0 . Later in this supplement, we will discuss the 'physical' algorithm based on nucleation, which applies to Fig. 4 of the main paper.
The input parameters to construct the synthetic surfaces are the numbers of rows and columns, L/l 0 , the typical size of the donor and acceptor patches, l/l 0 , and the donor probability p. The size l 0 corresponds to a single matrix element in our simulations. To fill our matrices, we first create a coarse matrix of size L/l × L/l and fill it with random values drawn from a normal distribution centered at zero and with standard deviation one (Suppl. Fig. 1a; if L/l is not an integer we round up). We use these values to fill every (l/l 0 ) th element of the final matrix, e.g. if l/l 0 = 10 then we fill points [i, j] = [0, 0], [10,0], [20,0], . . . , [10,0], [10,10], . . . etc. We assign values to the matrix elements between by performing a 5 th order 2-dimensional interpolation (Suppl. Fig. 1b). Next, we threshold this smoothly varying field to make a binary field. The value of the threshold, T , depends on the probability p, and is determined by the equation After solving this equation numerically, we assign all matrix elements with values less than T to be 1 (donor), and all with values higher to be 0 (acceptor). This results in a final binary matrix, as shown in Suppl. Fig. 1c or in Fig. 1 of the main paper. Qualitatively, these surfaces have heterogeneous features reminiscent of those seen by KPFM in Ref. [2]. We perform several tests to confirm that our procedure creates surfaces with the characteristics we desire. First, in Suppl. Fig. 2a, we show that the measured probability for being a donor (i.e. p meas = N d /N for several hundred generated matrices) is approximately equal to the input probability, p in . We therefore can simply write p meas = p in ≡ p. To ensure the length scale of the generated donor/acceptor features matches our input, l in , we measure the radial correlation function, C(r), exactly as prescribed in the main text in the context of our nucleation simulations, i.e., where s is the site identity (0 for acceptors and 1 for donors), R is the position of an element, r is the distance away from that element, and the averages are done over all positions separated by a distance r. Suppl. Fig. 2b shows C(r) for a few values of l in where, as mentioned in the main paper, we take the zero crossing to define l meas . In Suppl. Fig. 2c we plot l meas as a function of p in . This plot reveals that the measured value, l meas , is stable for intermediate values of p in , but exhibits deviations for p in > 0.9 and p in < 0.1; we therefore restrict our work to this range. In Suppl. Fig. 2d, we plot the averages over the stable region to show that l meas = (1.01 ± 0.03)l in . We therefore can simply write l meas = l in ≡ l. This validates that our 'synthetic' protocol allows us to create random surfaces with features of a desired size.
Finally, we make a brief aside to look at the behavior of the number of donors, N d , as a function of the length scale l/l 0 . Supplementary Fig. 2e shows distributions in N d for different l/l 0 . As can be seen, the mean value of N d (or equivalently the probability for a site to be a donor, p) does not change with l (Suppl. Fig. 2f). However, the widths of these distributions do change with l (inset of Suppl. Fig. 2f). As mentioned in the main paper, such broadening turns out to be the critical ingredient in how the introduction of the intermediate length scale l leads to charge-transfer enhancement. Left-to-right outcome probability  We begin by determining the distribution in the total charge transfer we expect in a first contact solely from the left surface to the right surface, P (Q ). Considering a given site pair [i, j], there are two relevant and mutually exclusive outcomes: (1) left-to-right transfer does occur, which has probability p(1 − p)α, and (2) left-to-right doesn't occur, which has probability 1 − p(1 − p)α (see Suppl. Table 1). The probability for transfer at any one site pair is independent from all others. Thus, for large L/l 0 , the probability distribution for total left-to-right transfer is Gaussian with mean e(L/l 0 ) 2 p(1 − p)α and width eL/l 0 p(1 − p)α(1 − p(1 − p)α), i.e., If we had instead begun by considering right-to-left transfer, we would have found similar probabilities (second part of Suppl. Table 1) and a corresponding distribution, P (Q ). Strictly speaking, the expressions for P (Q ) and P (Q ) are only valid when considering left-to-right and right-toleft transfer independently. This is because, at a given site pair, left-to-right and right-to-left transfer cannot occur simultaneously. The exact way to handle this is to consider the situation where three mutually exclusive outcomes are possible: (1) right-to-left transfer occurs, (2) left-to-right transfer occurs, or (3) neither right-to-left nor left-to-right transfer occurs. These outcomes and their corresponding probabilities are listed in Suppl. Table 2. Rather than solving this three-outcome problem exactly, we make progress by assuming that left-to-right transfer and right-to-left transfer are independent. This leads to four outcomes, with probabilities as in the column 'relaxed prob.' in Suppl. Table 2. In this case, the probabilities are products of those in Suppl. Table 1, and all involve terms on the order of 1, p(1 − p)α, or (p(1 − p)α) 2 . At its maximal value, which occurs when p = 0.5 and α = 1, the squared term is 1/4 th as large as the first order term. However, the experiments by Apodaca [1] suggest α is much smaller than one (more on the order of 0.1), thus realistically we expect the second order term to be 1/40 th as large as the first order term. This means that we can ignore the squared terms, which renders the relaxed case approximately equal to the strict case (column 'approx. relaxed prob.' in Suppl. Table 2).
The preceding paragraph reveals that, to a good approximation, we can treat the distributions P (Q ) and P (Q ) as independent. Since they are both Gaussian distributed, so too is P (∆Q), but with mean ∆Q = Q − Q = 0 and width σ = σ 2 + σ 2 . If we again neglect terms on the order of (p(1 − p)α) 2 , this leads to the l = l 0 probability distribution, As the mean value of ∆Q is zero, |∆Q| is necessarily proportional to the distribution width. The exact expression is This result therefore reproduces the |∆Q| ∝ L ∝ √ A scaling found previously [1], except for the square root dependency in α. In Suppl. Fig. 3a and 3c, we run multiple simulations to plot |∆Q| vs. L/l 0 for several values of α and p. In Suppl. Fig. 3e, we show that rescaling by the appropriate factor, F = 4p(1 − p)α/π, collapses all of our simulated numerical data for l = l 0 onto a single master curve.

First contact when l >> l0
In the analysis for l = l 0 , we started by considering the total left-to-right charge transfer, Q . A fundamental assumption was that the probability for transfer at any site pair was independent from all others. When patches are present, we can no longer make this assumption. This is because if a given site is a donor (acceptor), it is now more likely that its neighbor is also a donor (acceptor). We account for this by imagining a rescaled matrix consisting of larger patches, each patch containing n = (l/l 0 ) 2 elementary sites. There is therefore a total of N = N/n such patches. We randomly assign the patches of this rescaled matrix as donors/acceptors (still with probabilities p and 1 − p) to properly account for correlations. If a patch is assigned as a donor (acceptor), then all n elementary sites within it are donors (acceptors). Critically, we will assume that the transfer of charge still occurs independently for each site with probability α.
We now consider how many donor patches on the left surface face acceptors the right, N (holding off on the effect of the transfer probability, α, until later). The probability that a donor patch faces an acceptor is given by p (1 − p), thus we can expect P (N ) to be Gaussian with mean L 2 /l 2 p(1 − p) and width L/l p(1 − p)(1 − p(1 − p)) (and similarly for P (N )). Next, we repeat the strategy from the l = l 0 case and treat N and N as independent. The probability distribution for N − N is Gaussian with zero mean and width L/l 0 2p For every donor patch that faces an acceptor patch, n = (l/l 0 ) 2 elementary sites interact. If l >> l 0 , then the mean number of transfers between patches, αn, is much larger than the fluctuations, nα (1 − α). Thus, to a good approximation, we can expect an amount of charge, eα(l/l 0 ) 2 , to be transferred for every donor patch that faces an acceptor patch, which leads us to Eq. 2 of the main text, The width of this distribution, σ = αeLl/l 2 0 2p(1 − p), differs from the l = l 0 case in two ways. First, its αdependence is linear as opposed to square root in Suppl. Eq. 4. Second, it grows proportionally with the length l. This ultimately derives from the increased variability in the number of donors on each surface (consistent with the simulation results in Suppl. Fig. 2d-f). The mean absolute charge transferred is given by In Suppl. Fig. 3b and 3d, we run many simulations to plot |∆Q| vs. L/l 0 for l = 10l 0 and several values of α and p.

Sequential contacts: average behavior
Experiments show that if the same two surfaces are contacted repeatedly, charge transfer continues until one has a final total charge, Q f , and the other has −Q f [2]. In our numerical simulations, we see similar behavior. Suppl. Fig. 4 shows plots of Q/Q f for one surface vs. the number of contacts with its partner, n c . As mentioned in the main text, the relative size of fluctuations on approach to the final value are larger for smaller l. For larger l, when the fluctuations are small, the curves approach a well-defined limiting curve.
The average behavior of sequential contacts is straightforward to explain. We again start by looking at a particular surface-pair instance, and consider the number of left-to-right and right-to-left facing donor acceptor pairs (N and N , respectively). Regardless of the size of l, the average amount of charge that will be transferred left-to-right in the first contact is αe (N − N ). For the subsequent contact, this reduces N by a factor (1 − α) (on average). Accordingly, in the second contact we expect an amount of charge transfer αe(1 − α) (N − N ). Continuing like this, we expect that during the i th contact the amount of charge transferred is αe(1 − α) i−1 (N − N ). The total charge transferred up to contact number n c is therefore given by the sum To get the limiting behavior, we evaluate the sum in the limit n c → ∞, which is given by Thus rescaling by Q f in any given trial, we expect that sequential contact curves should approach the geometric series More conveniently, when α << 1 the charge accumulation can be approximated as a continuous expression. Assuming that the number of contacts n c is a continuous variable, an infinitesimal increase in Q can be written By integrating this equation, we find The accumulated charge on the surface is given by considering that Q f = eN 0 . In Suppl. Fig. 4, we plot this prediction on top of our simulated curves, which illustrates that it accurately describes the average behavior.

Sequential contacts: fluctuations
Since the arrangement of facing donors/acceptors is assumed fixed, the fluctuations of charge transferred in sequential contacts must arise solely from the transfer probability, α. To understand how, we start by considering a specific situation with two surfaces. We denote the number of donors on the left (right) surface that face acceptors on the right (left) as N (N ). We imagine contacting the two surfaces repeatedly, each time measuring the actual number of transfers left/right as N * , and each time re-initializing. The total charge transferred in one of these trials is The mean charge transferred is ∆Q = eα(N − N ), and is in general non-zero since, in a particular instance, it is not true that N = N . We are interested in the fluctuations about the mean, which arise due to the fluctuations σ * 2 = N α(1 − α) in the amounts of charge transferred right-to-left and left-to-right. The fluctuations about the mean are given by We would now like to use this result for a particular surface-pair instance to determine the average size of fluctuations, δQ, for an ensemble of pair instances. The exact way to do this would be to evaluate weighted average Conveniently, we can use a shortcut to avoid evaluating this integral. In our simulations, the number of sites is large (and even more so in a real experiment). This means that fluctuations of N (given by σ = N p(1 − p) (1 − p(1 − p))) are small compared to the mean (given by N = N p(1 − p)). Therefore the only region where the contributions to the integral in Suppl. Eq. 18 are appreciable is near to the mean; accordingly, we can simply replace N by N in Suppl. Eq. 17, which gives Eq. 3 of the main text, i.e., In Suppl. Fig. 5, we run many simulations to plot δQ vs. L/l 0 for l = l 0 and l = 10l 0 and several values of α and p. Rescaling by the factor, 2p(1 − p)α(1 − α), we confirm we collapse all of our numerically simulated data (Suppl. Fig. 5e).

NUCLEATION
In Fig. 4 of the main paper, a physically-derived nucleation process is used to generate heterogeneous surfaces, which can then be characterized in terms of donor probability p and length scale l/l 0 and used in contact charging simulations. This process is formulated from general energy considerations, and without assuming the nature of the sites. The main assumption of this model is that the probability to form a site is dependent on the presence of neighboring sites of the same type. Indeed, the length scales observed in the experiment seem to indicate that donor and acceptor sites tend to form patches of more than one site. One possible explanation is that this situation is more energetically favorable than randomly distributed sites. We will therefore assume that the probability for a site to form increases with the number of neighboring sites of the same type.
Suppose that the rates at which sites of each type form follows the general Arrhenius-like expression where the activation energy E is a function of the number of neighbors. This rate will determine the probability to form a new site at each time step. Considering only the first neighbors, we denote ν D = ν ≤ 4 the number of neighboring donors, ν A = 4 − ν the number of neighboring acceptors, and P A (ν) (P D (ν)) the probability for an acceptor (donor) to form on the surface. For example, we can consider that one species on the surface could act as a donor site, and the absence of that species would define an acceptor site. In that case, P D (0) would correspond to the probability for one such site to form on a cell surrounded by donor sites, while P A (0) would correspond to the probability for a single isolated donor site to disappear. We assume that each neighbor affects the activation energy linearly, and denote K the change in energy due to the addition or removal of a single neighboring site. In our example, K can therefore be understood as a (dimensionless) interaction energy between two adjacent donor sites. If donors are more likely to form near donors, we can expect P D to increase monotonously with ν. Conversely, P A should decrease monotonously with ν. We therefore propose the general expression In the main paper, we also have P A, 0 = P D, 0 = P 0 , so that the donor/acceptor probabilities are symmetric, reducing the parameters to just K and P 0 . All possible configurations are shown in Suppl. Fig. 6. 6. We assume that the probability for a given site (dotted cell) to change state is dependent on the configuration of the first neighbors (Suppl. Eq. (21)). If the site is an acceptor (resp. donor), its probability to become a donor (resp. acceptor) is given by PD(ν) (resp. PA(ν)), where ν is the number of donors among the first neighbors. With one more neighboring donor, the activation energy decreases (resp. increases) by K, which in turns increases (resp. decreases) the probability of transition.

Protocol
The characteristics of a surface are its total size L/l 0 , the typical size of the patches l/l 0 and the probability for a site to be a donor p. First, a matrix of size L/l 0 × L/l 0 is randomly filled with donor and acceptor sites. The output value of p will be determined by the equilibrium between the creation of new donor and acceptor sites. Assuming P A, 0 = P D, 0 = P 0 , we can expect p out ≈ 0.5. We therefore choose p in = 0.5 in order to reach the dynamic equilibrium faster. The patch size l/l 0 is influenced by K, P A, 0 and P D, 0 . While measuring p is straightforward, determining l/l 0 is not trivial. The spatial correlation function is calculated on the surface, with s = 1 for a donor and 0 for an acceptor. The first zero crossing therefore corresponds to the typical correlation length between sites. A typical run consists in the following steps. The input parameters are the change in activation energy due to a single neighbor K and the probability P 0 .
1. Either create a random surface of size L/l 0 and probability p in or start from a previously generated surface. In the main paper, we compute increasing values of K and use the previous surface as a starting point.
2. For each cell on the surface: (a) count the number of neighbors ν; (b) calculate the probability to change state (P A (ν) or P D (ν)); (c) either change state or not.
3. Repeat step 2 until an equilibrium is reached (5000 iterations in the main paper).
4. Determine p and calculate the correlation function C(r) as defined in Suppl. Eq. 22.
5. Determine l/l 0 by averaging C(r) over several surfaces and considering the first zero crossing within the standard deviation (see Fig. 4 of the main paper).
The steps of a given run are shown in Suppl. Fig. 7. The parameters are L/l 0 = 100, p in = 0.5, P 0 = 0.5 and K = 0.8. In Suppl. Fig. 7a-c, we see the emergence of a length scale on the surface. This increases the typical fluctuations of p, which is given by σ p = l/L p(1 − p). Supplementary Fig. 7d shows the evolution of p for a few example surfaces. On a long enough timescale, a given surface will explore all of the possible values of p, which can be seen in Suppl. Fig. 7e. The evolution of the length scale l/l 0 is given in Suppl. Fig. 7f. It starts close to unity, and grows to reach an equilibrium value after typically 100 to 1000 iterations. This equilibrium value of l/l 0 depends on the input parameters, particularly K.

CHARGING WITH MULTIPLE LENGTH SCALES
The KPFM data by Baytekin, et al. [2], reveals two length scales for the patches: one at ∼45 nm, and another at ∼450 nm. This brings up a natural question: if two (or more) length scales are present at the meso scale, how is the charging behavior affected? Here we will argue that the larger length scales dominate. We generate a surface with two donor/acceptor length scales by populating a matrix with normal random values on two sub-matrices, one at every (l 1 /l 0 ) th point, and another at (l 2 /l 0 ) th point (l 1 < l 2 ). Interpolating between these points at the l 0 scale and thresholding as before, we create surfaces such as shown in Suppl. Fig. 8c. This is one of many ways to introduce two length scales-looking into all such possible ways is beyond the scope of our work. Qualitatively, however, this procedure (1) preserves the feature p meas = p in and (2) creates surfaces where the smaller donor (acceptor) patches 'intrude' into the larger acceptor (donor) patches. In this latter point, it is qualitatively similar to the KPFM data presented by Ref. [2].
In Suppl. Fig. 8d, we show how |∆Q| is affected when the contacting surfaces have two length scales. Holding the smaller scale, l 1 , fixed and increasing the larger one, we see that the magnitude of charging increases for increasing l 2 . When l 2 >> l 1 , the growth in |∆Q| is approximately linear with l 2 /l 0 (inset to Suppl. Fig. 8d), which indicates that if the scales are sufficiently separated then the larger one dominates. In Suppl. Fig. 8e, we attempt to collapse our data by the previously found factor αl 2 /l 0 4p(1 − p) (but where we've replaced l by l 2 ). This collapses all data for l 2 >> l 1 , but with a prefactor that is lower than in the single length scale case. Qualitatively, this can be explained by considering that the larger donor (acceptor) patches are effectively smaller than l 2 given the intrusive acceptor (donor) regions of size ∼ l 1 that now exist inside of them. Nonetheless, under the conditions we examine here this effect is relatively small, and to within a factor of order 1 our model can still account for the charging behavior simply  Fitting each curve (dashed colored lines) to |∆Q| = kL/l0 shows that the prefactor, k, is proportional to l2/l0 for l2 >> l1 (inset, where black dashed line represents k ∝ l2/l0). (e) Rescaling all data by the factor αl2/l0 4p(1 − p) collapses data in the l2 >> l1 regime, albeit with a prefactor, s, that is less than 1. This indicates that, while charge transfer is dominated by larger donor/acceptor length scales if multiple length scales are present, the efficiency is somewhat reduced.
by assuming one length scale at l = l 2 . The presence of multiple length scales might explain the imperfect collapse from Fig. 4e in the main text. For some values of l/l 0 , the measured |∆Q| is smaller than the prediction by a factor of order 1, similarly to what we see in Suppl. Fig. 8e. To determine if it is due to the presence of multiple length scales, we remove features smaller than the dominant length scale on each surface, then perform the same single-contact measurement again. To remove the small features, we multiply the Fourier transform of every surface by a Butterworth low-pass filter of order 5 where the cutoff is set at the wavenumber that corresponds to the dominant length scale. In Suppl. Fig. 9a, we compare the collapse of |∆Q| by the prefactor F = αl/l 0 4p(1 − p)/π with and without the filter. One can see that the collapse is significantly improved by the removal of features smaller than l/l 0 , which suggests the presence of smaller length scales that are not accounted for in the model. An example of filtering is shown in Suppl. Fig. 9a, where the filtered surface (shades of gray) is superimposed with the original surface. The sites that have been affected by the filter are shown in light brown (acceptor to donor) and yellow (donor to acceptor).

RECONSIDERATION OF PREVIOUS EXPERIMENTAL RESULTS
In this section, we compare our model with the experimental data from Apodaca et al. [1]. Two types of experiments have been performed: single-contact, and multiple successive contacts, as show in Suppl. Fig. 10a and 10b, respectively. In multiple-contact experiments, the charge increases progressively until reaching a plateau, as described by Suppl. Eq. 12 and 15. In Supplementary Fig. 10b, we fit the continuous expression from Suppl. Eq. 15 to the experimental data from [1]. We find α = 0.143 and Q f = [0.972, 1.489, 2.187], where each value for Q f correspond to a different area for the surface (respectively 16.00, 35.22 and 61.62 mm 2 ). We verify that Q f ∝ √ A. The value of α obtained can  then be used to determine the length scale l 0 using the single-contact data in Suppl. Fig. 10a. Indeed, the average charge after a single contact is given by Suppl. Eq. 7 and depends on L, l 0 , l, p and α. We already know L and α, and the KPFM data from Baytekin et al. [2] gives us a value for l. Baytekin et al. identified two length scales, 450 nm and 44 nm. However, as established in the previous section, charging is dominated by the larger length scale, granted that they are sufficiently separated. We therefore consider here that l = 450 nm. Finally, we make the assumption that donors and acceptors are equally probable, i.e. p = 0.5. The fit of Eq. 7 on the data from [1] shown in Suppl. Fig. 10a gives us l 0 = 3.67Å. In other words, a square of side ∼ 4Å would donate/accept a single elementary charge. By comparison, if we do not consider the influence of donor/acceptor patches, we can fit the data using Eq. 5 which yields l 0 = 0.003Å. Similarly, the value found in [1] corresponds to l 0 = 0.005Å, differing slightly due to the different α dependency. This would imply an unreasonably high number of charges (∼ 10 5 ) on an area the size of an atom. Furthermore, one can verify that no combination of α and p would yield a plausible value for l 0 , indicating that the introduction of a length scale l > l 0 is necessary.