# Randomly sampling bipartite networks with fixed degree sequences

Zachary P. Neal<EMAIL_ADDRESS>Michigan State University

###### Abstract

Statistical analysis of bipartite networks frequently requires randomly
sampling from the set of all bipartite networks with the same degree sequence
as an observed network. Trade algorithms offer an efficient way to generate
samples of bipartite networks by incrementally ‘trading’ the positions of some
of their edges. However, it is difficult to know how many such trades are
required to ensure that the sample is random. I propose a stopping rule that
focuses on the distance between sampled networks and the observed network, and
stops performing trades when this distribution stabilizes. Analyses
demonstrate that, for over 300 different degree sequences, using this stopping
rule ensures a random sample with a high probability, and that it is practical
for use in empirical applications.

###### keywords:

sampling , monte carlo , mixing time , bipartite

††journal: Submitted to arXiv

##  1 Introduction

Statistical analysis of bipartite networks frequently requires randomly
sampling bipartite networks with fixed degree sequences [e.g., Neal et al.,
2021, Jasny, 2012, Gotelli, 2000]. The sampled networks provide a comparison
(i.e. a null model) against which to compare and thereby evaluate the
statistical significance some network property of interest, such as the
nestedness of a bipartite network [Strona & Fattorini, 2014] or the strength
of an edge in a bipartite projection [Neal, 2014]. Trade algorithms offer an
efficient way to iteratively randomize a bipartite network while preserving
its degree sequences by trading the positions of edges [Strona et al., 2014,
Godard & Neal, 2022]. If enough trades are performed, these algorithms are
known to yield random samples of all such networks [Carstens, 2015]. However,
the number of trades required to ensure that the sample is random remains
unknown.

In this paper, I propose a stopping rule for these algorithms that focuses on
the distribution of sampled networks’ distances from the observed network.
Using this rule, trades are performed until this distribution has stabilized
according to a Kolmogorov–Smirnov (KS) test of distribution equality
[Kolmogorov, 1933, Smirnov, 1948]. Examining this stopping rule’s performance
in over 300 degree sequences, I demonstrate that it yields a random sample
over 90% of the time.

The paper is organized in four sections. In the background section I review
trade algorithms for randomizing and sampling bipartite networks, illustrating
their performance in a small bipartite network. I then propose a stopping rule
for deciding how many trade randomizations to perform. In the methods section,
I describe a test of its performance for bipartite networks with varying
degree sequences. In the results section, I present the results of this test,
and illustrate the stopping rule’s practicality using empirical examples drawn
from ecology, sociology, and political science. Finally, in the discussion
section, I summarize the findings, and discuss the method’s limitations and
implications.

##  2 Background

###  2.1 Sampling bipartite networks using trade randomization

A network is bipartite when its nodes can be partitioned into two sets, such
that edges exist only between nodes in different sets. A bipartite network is
also two-mode when these two sets of nodes represent distinctly different
kinds of entities. Such networks are common in the literature, and can capture
many kinds of phenomena including event attendance where people are connected
to events [Breiger, 1974], bill co-sponsorship where legislators are connected
to bills [Neal, 2014], ecosystems where species are connected to habitats
[Gotelli, 2000], and scholarship where authors are connected to papers
[Newman, 2001]. Following Latapy et al. [2008], I generically refer to one set
of nodes as ‘top nodes,’ and to the other set of nodes as ‘bottom nodes.’ A
bipartite network is characterized by two degree sequences, one for the top
nodes, and another for the bottom nodes. For example, in a bipartite network
capturing event attendance, where people are top-nodes and events are bottom-
nodes, the top-node degree sequence captures the number of events attended by
each person, while the bottom-node degree sequence captures the number of
people attending each event.

For any pair of degree sequences, there exists a (possibly empty) set
$\mathcal{B}$ of all bipartite networks that have the given degree sequence.
For example, given a top node degree sequence {1,1,2} and a bottom node degree
sequence {1,1,2}, $\mathcal{B}$ contains the five networks shown in Figure 1A,
and therefore the size or cardinality of $\mathcal{B}$ is five (i.e.,
$|\mathcal{B}|=5$). The goal of randomly sampling bipartite networks with
fixed degree sequences is to sample networks from $\mathcal{B}$ such that each
member of $\mathcal{B}$ has an equal probability of being sampled. The degree
sequences {1,1,2} and {1,1,2} define a set that is small enough that it is
possible to enumerate all its members, making it easy to directly sample from
this list. However, because longer degree sequences define sets that are too
large to enumerate, and indeed whose sizes are generally unknown [Barvinok,
2010], random sampling in such cases requires a different approach.

Many different approaches for sampling have been proposed [Verhelst, 2008,
Miller & Harrison, 2013, Admiraal & Handcock, 2008, Bezáková et al., 2007,
Chen et al., 2006, Fout et al., 2020, Wang, 2020, Strona et al., 2014, Godard
& Neal, 2022], however they vary in the evidence of their randomness, their
computational complexity, and their practicality in empirical contexts. The
most widely-used sampling approach involves starting with a network
(typically, the observed network being analyzed) that has the given degree
sequences, then repeatedly rewiring its edges in a way that preserves the
degree sequence, to obtain a new network that can be regarded as randomly
sampled from $\mathcal{B}$. Among the most efficient algorithms proven to
perform such rewiring randomly are the ‘trade’ algorithms [Strona et al.,
2014, Godard & Neal, 2022, Carstens, 2015], so-called because they trade edges
in the same way that children might trade baseball cards.

Figure 1B illustrates each step of the trade algorithm. First, the starting
network (e.g., Network B from Figure 1A) is represented as a neighbor list
that indicates, for each top-node, the bottom-nodes to which it is connected.
Next, the algorithm picks two top-nodes at random from the neighbor list
(e.g., nodes Y and Z). Then, the two chosen top-nodes trade a random number of
their unique neighbors. A single ‘trade’ can involve swapping none, one, or
many neighbors, which makes trade algorithms more efficient than checkerboard
swap algorithms that swap only one neighbor in each iteration [Ellison et al.,
2024]. The outcome of the trade is represented as a new list. The process can
be repeated, performing additional trades among pairs of randomly-selected
top-nodes. Two variants of the trade algorithm exist that differ only in how
efficiently they perform these steps: The ‘curveball’ [Strona et al., 2014]
algorithm performs trades in O(n log n) time, while the newer ‘fastball’
[Godard & Neal, 2022] algorithm performs them slightly more efficiently in
O(n) time. After performing the chosen number of trades, the resulting
neighbor list represents an ending network. In this example, performing one
trade transformed the starting network into a new network with the same degree
sequences (e.g., Network A from Figure 1A), or put another way, resulted in
sampling Network A from the set of five possible networks.

Figure 1: (A) The set $\mathcal{B}$ of all bipartite networks with top-degree
sequence {1,2,1} and bottom-degree sequence {1,2,1}; (B) Example of bipartite
network randomization using a trade algorithm.

Figure 2 illustrates the composition of a sample of 1000 networks drawn using
the fastball trade algorithm, starting from Network B. Performing only one
trade is insufficient to obtain a random sample of $\mathcal{B}$ because, as
Figure 2A shows, it yields a sample in which Network B is over-represented and
Network C is under-represented. The statistically significant $\chi^{2}$
statistic confirms that there is significant variation in the relative
frequency of each network in the sample. This variation remains even after two
(Figure 2B) or five (Figure 2C) trades have been performed. However, after ten
trades have been performed (Figure 2D), each network occurs roughly the same
number of times. The uniformity of the sample is confirmed by a non-
significant $\chi^{2}$ statistic, which indicates that the sample is now a
random sample from $\mathcal{B}$.

Carstens [2015] proved that trade algorithms will always yield a random sample
from $\mathcal{B}$, provided a sufficient number of trades are performed.
However, although trade algorithms are known to yield random samples, the
number of trades required to ensure random sampling (i.e., the mixing time) is
unknown and “intractable to compute” (p. 5). Therefore, in the next section, I
propose a stopping rule to decide when a sufficient number of trades have been
performed.

Figure 2: Distribution of networks in a sample drawn using a trade algorithm,
following the given number of trades, with a $\chi^{2}$ test of uniform
distribution.

###  2.2 A trade randomization stopping rule

As Figure 2 illustrates, when $|\mathcal{B}|$ is small, a $\chi^{2}$ test can
be used to determine when enough trades have been performed to obtain a random
sample from $\mathcal{B}$. However, typically $|\mathcal{B}|$ is too large to
enumerate all its members, and the $\chi^{2}$ test cannot be used to verify a
sample’s randomness or to know when to stop performing trades. Therefore, a
different method is required for deciding when a sufficient number of trades
have been performed to obtain a random sample.

One alternative focuses not on the frequency of each network in the sample,
but on the distance between each sampled network and the starting network. The
distance $d$ between two bipartite networks is defined as the fraction of
dyads that are different, or when the networks are represented as binary
incidence matrices, the fraction of cells that are different [Strona et al.,
2014]. For example, figure 1 illustrates that the distance between the
starting network (Network B) and the ending network (Network A) is
$d_{A,B}=0.444$.

Given a sample of bipartite networks obtained by randomizing a starting
network using a trade algorithm, $D_{t}$ is the list of distances $d$ between
the starting network and each network in the sample after performing $t$
trades. Figure 3 illustrates the distribution of $D$ in each sample whose
composition is shown in figure 2. For example, figure 3A shows that, after one
trade, the sample contains many networks at a distance of 0 from the starting
network (these are Network B), some networks at a distance of 0.444 (these are
Networks A, D, and E), and no networks at a distance of 0.666 (these would
have been Network C). Note that this pattern mirrors the frequency of networks
in the sample shown in figure 2A: many instances of Network B, some instances
of Networks A, D, and E, and no instances of Network E.

Figure 3: Distribution of distance between a starting matrix and networks in a
sample ($D$) drawn using a trade algorithm, following the given number of
trades, with a Kolmogorov–Smirnov (KS) test of equality with an earlier
distribution.

As more trades are performed, the distribution of $D$ changes. Specifically,
because trade algorithms randomly sample from $\mathcal{B}$ [Carstens, 2015],
the distribution of $D$ converges on the distribution of distances
$\mathcal{D}$ between the starting network and every member of $\mathcal{B}$.
This convergence is visible in figure 3, and in this small example can be
verified because $\mathcal{D}$ is known: {0,0.444,0.444,0.444,0.666}. When the
distribution of $D$ matches the distribution of $\mathcal{D}$, this suggests
that the members of the sample are representative of, and thus a random sample
from, $\mathcal{B}$.

Typically $\mathcal{D}$ is unknown, so it is not possible to directly test
whether the distribution of $D$ matches the distribution of $\mathcal{D}$.
However, because $D$ converges on $\mathcal{D}$, it is only necessary to test
when the distribution of $D$ has stabilized. A Kolmogorov–Smirnov (KS) test of
the equality of distributions [Kolmogorov, 1933, Smirnov, 1948] provides such
a test when it is used to compare the distributions of $D_{t}$ and $D_{t-n}$.
For example, the KS test reported in figure 3B ($KS=0.175,p<0.001$) shows that
the distribution of $D$ after two trades is statistically significantly
different from the distribution of $D$ after one trade; the distribution has
not stabilized. In contrast, the KS test reported in figure 3D
($KS=0.023,p=0.954$) shows that the distribution of $D$ after ten trades is
not statistically significantly different from the distribution of $D$ after
five trade; after ten trades the distribution has stabilized. The stability of
the distribution of $D_{1}0$ suggests that the sample is random after 10
trades, which was confirmed directly using a $\chi^{2}$ test in figure 2D.

Therefore, I hypothesize that sampling bipartite networks with fixed degree
sequences by performing trades until the distribution of $D_{t}$ is not
statistically significantly different from the distribution of $D_{t-n}$ will
ensure with high probability that the resulting sample is a random sample from
$\mathcal{B}$.

##  3 Method

###  3.1 Validation

The small example above illustrates one set of degree sequences (i.e.,
{1,1,2}, {1,1,2}) where, when the distribution of $D_{t}$ became not
statistically significantly different from the distribution of $D_{t-n}$ (see
figure 3D), the sample of networks was a random sample from $\mathcal{B}$ (see
figure 2D). I follow a similar approach to test the hypothesis that the
proposed stopping rule yields random samples, evaluating how often it yields a
random sample when used to sample bipartite networks with other degree
sequences.

For each degree sequence, I apply the stopping rule to draw a sample of 1000
bipartite networks, then evaluate whether the resulting sample is a random
sample from $\mathcal{B}$. I repeat this 1000 times, recording the percentage
of samples that are random. Applying the stopping rule requires computing $D$
and performing a KS test after every $n$ trades. In these analyses, to avoid
needing to perform these computations too often, $n$ is twice the number of
top nodes. The statistical significance of the KS test is evaluated using the
conventional $p<0.05$ threshold.

To determine whether a given sample is random sample from $\mathcal{B}$, I
require that it exhibit two properties: coverage and uniformity. First, the
sample must include at least one instance of every member of, and therefore
completely cover, $\mathcal{B}$. Second, the relative frequency of each
network in the sample must be equal, and therefore be uniformly distributed,
as indicated by a non-significant $chi^{2}$ test.

Using this method for determining whether a sample is a random sample from
$\mathbf{B}$ requires that $|\mathcal{B}|$ is small enough for its members to
be enumerated, and for a random sample of 1000 networks to plausibly contain
at least one instance of each member. Therefore, I restrict the scope of this
analysis to 335 different degree sequences where $|\mathcal{B}|<100$.

###  3.2 Application

The validation test evaluates whether the stopping rule consistently ensures a
random sample, but it is only possible for shorter degree sequences that
define smaller bipartite networks. To examine the stopping rule’s practicality
in larger empirical bipartite networks, I illustrate its use to draw samples
of bipartite networks in three different contexts: ecology, sociology, and
political science.

First, I consider the ‘Darwin’s Finches’ network, which records whether Darwin
observed each of 13 species of finches living on each of 17 islands
[Sanderson, 2000]. It is defined by the top-node degree sequence
{14,13,14,10,12,2,10,1,10,11,6,2,17}, which describes the number of islands on
which each species of bird was observed, and by the bottom-node degree
sequence {4,4,11,10,10,8,9,10,8,9,3,10,4,7,9,3,3}, which describes the number
of species that were observed on each island. In the ecological context,
randomly sampling bipartite networks with these degree sequences is necessary
to identify, for example, species’ significant co-habitation on islands that
suggests their symbiosis or lack of competition. Interestingly, Chen et al.
[2005] computed that $|\mathcal{B}|=67,149,106,137,567,626$ for these degree
sequences. Thus, even this modestly-sized bipartite network highlights the
large size of $\mathcal{B}$ from which networks are to be randomly sampled,
and thus the difficulty in determining whether the sample is random.

Second, I consider the ‘Southern Women’ network, which records whether each of
18 women attended each of 14 social events [Homans, 1950, Breiger, 1974]. It
is defined by the top-node degree sequence
{8,7,8,7,4,4,4,3,4,4,4,6,7,8,5,2,2,2}, which describes the number of events
each woman attended, and by the bottom-node degree sequence
{3,3,6,4,8,8,10,14,12,5,4,6,3,3} who attended each event. In the sociological
context, randomly sampling bipartite networks with these degree sequences is
necessary to identify, for example, women’s significant co-attendance at
events that suggests their friendship.

Finally, I consider the ‘US Senate’ network, which records whether each of 100
US Senators sponsored each of 5357 in the 117th session of the U.S. Senate
[Fowler, 2006, Neal, 2022]. It is defined by a top-node degree sequence that
describes the number of bills each Senator sponsored, and by a bottom-node
degree sequence that describes the number of Senators who sponsored each bill.
In the political science context, randomly sampling bipartite networks with
these degree sequences is necessary to identify, for example, Senators’
significant co-sponsorship of bills that suggests their political alliance.

Each of these empirical examples involve a $\mathcal{B}$ that is too large to
use a $chi^{2}$ test to determine whether a sample obtained using the stopping
rule is random. Instead, in these illustrative applications, I focus on three
outcomes that are observable. First, does the distribution of $D$ converge to
a stable distribution that is expected to match the distribution of
$\mathcal{D}$? Second, is each member of the resulting sample unique,
suggesting that the sample provides broad coverage of $\mathcal{B}$ and does
not over-represent certain members of $\mathcal{B}$. Third, can a likely-to-
be-random sample be drawn in a practical amount of time using the stopping
rule?

##  4 Results

###  4.1 Validation

Table 1 reports, for 10 selected degree sequences and for the average over all
335 degree sequences, the percent of samples that were random when the
stopping rule is used to stop performing trade randomizations. Each row
represents a pair of bipartite network degree sequences, which define a
$\mathcal{B}$ containing a specific number of networks ($|\mathcal{B}$). For
example, the first row corresponds to the small example used throughout the
background section. In this small example, using the proposed stopping rule to
decide when to stop performing trade randomizations yielded a random sample of
bipartite matrices from $\mathcal{B}$ 94.7% of the time. Using the stopping
rule for other degree sequences yielded a random sample with similar
frequency, and on average over 335 different degree sequences it yielded a
random sample 93.2% of the time. This finding provides supports the
hypothesis: the proposed stopping rule ensures with high probability that the
resulting sample is a random sample from $\mathcal{B}$.

Table 1: Performance of stopping rule for selected degree sequences Degree Sequences |  | Mean Trades | Random
---|---|---|---
Top Nodes | Bottom Nodes | $|\mathcal{B}|$ | Required | Samples
{1,1,2} | {1,1,2} | 5 | 12 | 94.70%
{1,2,3} | {1,1,2,2} | 8 | 14 | 94.20%
{2,2,3} | {1,1,1,2,2} | 31 | 12 | 95.60%
{3,3,3} | {1,1,1,2,2,2} | 93 | 12 | 95.10%
{1,4,6} | {1,1,1,2,2,2,2} | 19 | 17 | 94.30%
{1,5,6} | {1,1,1,1,2,2,2,2} | 52 | 16 | 94.80%
{1,2,2,2} | {1,1,2,3} | 27 | 22 | 94.80%
{1,1,2,3} | {1,1,1,2,2} | 68 | 23 | 94.90%
{1,1,3,5} | {1,1,1,2,2,3} | 26 | 25 | 93.30%
{2,2,4,4,4} | {2,2,4,4,4} | 72 | 31 | 93.50%
… | … | … | … | …
Average for 335 sequences | 43 | 23 | 93.20%
Results for all 335 degree sequences are provided in Supplementary
Information.

###  4.2 Application

Figure 4 shows the distribution of $D$, after a specified number of trade
randomizations, in a sample of bipartite networks with the same degree
sequences as the (A) Darwin’s Finches, (B) Southern Women, and (C) US Senate
networks. In each case, the stopping rule stopped performing trade
randomizations when this distribution stabilized, which required a different
number of trades depending on the network (e.g., 156 trades in the finches
network, 1400 trades in the Senate network). It is impossible to verify that
the resulting samples are random samples from the corresponding $\mathcal{B}$.
However, it is possible to examine whether the stopping rule performs as
expected and in practical time.

Figure 4: Distribution of sampled networks’ distance from observed (A)
Darwin’s Finches, (B) Southern Women, and (C) US Senate networks following
trades

First, as expected, in all three cases the distribution of $D$ converges to a
stable distribution. Notably, the distribution on which it converges is unique
to the degree sequence, which is consistent with the expectation that the
distribution of $D$ converges on the distribution of $\mathcal{D}$.

Second, as required of a random sample, all 10,000 members of these samples
are unique. This suggests that applying the stopping rule yields a sample with
broad coverage of $\mathcal{B}$, and does not yield a sample in which any
networks from $\mathcal{B}$ are over-represented.

Finally, using the stopping rule yielded a random sample of 10,000 bipartite
networks in a practical amount of time: less then 4 seconds for the Finches
and Women, and about 600 seconds (10 minutes) for the Senate. This is roughly
3.5 times longer than performing the same number of trades without
periodically computing $D$ or a KS test (i.e., without applying the stopping
rule). Although the additional computations required to use the stopping rule
require more time, they offer confidence in the resulting sample’s randomness.
For example, although 1400 trades can be performed on 10,000 copies of the
Senate network in just 175 seconds, the resulting sample would come with no
certainty about its randomness because more (or fewer) trades may have been
required.

##  5 Discussion

Researchers performing statistical analysis of a bipartite network need to
draw a random sample of all bipartite networks with the same degree sequences
as the observed network. Efficient ‘trade’ algorithms exist for sampling
bipartite networks through repeated randomization trades [Strona et al., 2014,
Godard & Neal, 2022]. However, it is difficult to know how many trades needed
to ensure the sample is random. In this paper, I proposed a stopping rule that
focuses on the distribution of the sampled networks’ distances from a starting
network, and stops performing trades once this distribution stabilizes. For a
wide range of degree sequences, this stopping rule ensures with high
probability (93.2% on average) that the sample is a random sample of all
bipartite networks with the given degree sequence. Additionally, in three
larger empirical examples, the stopping rule performs as expected, and returns
a sample in a practical amount of time.

The stopping rule offers a solution to a practical challenge encountered by
researchers: determining when a sample of bipartite networks is random,
regardless of how many trades this might require. These analyses do not
directly answer a related question: how many trades are required, or what is
the mixing time of a trade algorithm? While the mixing time of trade
algorithms remains unknown, these results offer some insight, while the
stopping rule offers a tool for future research on their mixing time.
Consistent with prior work [Strona et al., 2014, Carstens, 2015], the number
of trades required to achieve a random sample according to the stopping rule
is closely associated with the number of top-nodes ($r=0.979$). However,
whereas prior work has suggested five times the number of top-nodes as an
approximate mixing time [Strona et al., 2014], these results suggest this is
an underestimate. This rule-of-thumb suggests that 25 trades should be
sufficient for a network containing with five top-nodes, however the stopping
rule required an average of 31.6 trades. Similarly, the rule-of-thumb suggests
that 500 trades should be sufficient for the Senate network containing 100
top-nodes, however the stopping rule required 1400 trades. While prior
approximations of trade algorithms’ mixing time may be too low, the stopping
rule may provide a way for future research to estimate mixing time more
accurately and for larger degree sequences.

These results must be interpreted in light of three limitations. First, it was
only possible to evaluate the performance of the stopping rule for degree
sequences that define a $\mathcal{B}$ for which $|\mathcal{B}|\leq 100$. This
is primarily a practical computational limitation; with additional computing
time, it is theoretically possible to evaluate it for longer degree sequences.
However, $|\mathcal{B}|$ can grow very rapidly as a function of network size
and in most cases is unknown [Barvinok, 2010], making it difficult to
enumerate all members of $\mathcal{B}$.

Second, although the stopping rule provides confidence that the resulting
sample is random, this confidence comes at the cost of greater time and memory
requirements. In terms of time, the additional computations necessary to use
the stopping rule (i.e., periodically computing $D$ and a KS test) take
approximately 3.5 times longer than the trade algorithm alone. In terms of
memory, if the required number of trades was known in advance, then a trade
algorithm can be used to sample networks from $\mathcal{B}$ one at a time. In
contrast, using the stopping rule requires holding all the sampled networks in
memory as they are incrementally randomized and compared in terms of their
distance from the starting network. However, as noted above, future research
may use the stopping rule to investigate the mixing time of trade algorithms,
thereby providing guidance on choosing a number of trades in advance and
eliminating the need for the stopping rule.

Finally, although on average the stopping rule yields a random sample in 93.2%
of cases, this still means that 6.8% of the time it yields a non-random
sample. While the risk that the stopping rule yields a non-random sample
should not be ignored, it is worth noting that the sample is likely still
nearly-random and that this level of error is close to the error tolerance
(e.g., $\alpha=0.05$) typically accepted in conventional statistical
inference. When greater certainty of a sample’s randomness is necessary, one
can impose a higher threshold on regarding the KS test as non-significant. The
results presented in table 1 view the KS test used by the stopping rule as
non-significant when $p>0.05$. However, if a $p>0.95$ threshold is used
instead, then on average the stopping rule yields a random sample in 95.17% of
cases. However, this increased certainty comes at a computational cost of
requiring more trades (on average 33 trades vs. 23).

The analysis of bipartite networks and their projections has become common in
a range of contexts [e.g., Breiger, 1974, Neal, 2014, Gotelli, 2000, Newman,
2001]. The statistical analysis of these networks typically requires building
null models based on random samples of bipartite networks with the same degree
sequence as an observed network. Although many methods exist for generating
samples of bipartite networks with fixed degree sequences, it is difficult to
be sure that these samples are random, that is, that they uniformly sample
from all bipartite networks with the given degree sequences. The proposed
stopping rule offers a way to be reasonably sure the sample is random, while
simultaneously ensuring that the randomization algorithm runs only as long as
necessary to ensure the sample’s randomness. It may also offer a way to more
precisely estimate the mixing time of trade randomization algorithms.
Therefore, it offers a useful tool for statistically analyzing bipartite
networks, and for investigating network randomization algorithms.

## Data availability statement

All data and code necessary to replicate this study are available at
https://osf.io/rmwpz/.

## Disclosures

The author declares no conflicts of interest.

## References

  * Admiraal & Handcock [2008] Admiraal, R., & Handcock, M. S. (2008).  Networksis: a package to simulate bipartite graphs with fixed marginals through sequential importance sampling.  Journal of Statistical Software, 24. doi:10.18637/jss.v024.i08.
  * Barvinok [2010] Barvinok, A. (2010).  On the number of matrices and a random matrix with prescribed row and column sums and 0–1 entries.  Advances in Mathematics, 224, 316–339. doi:10.1016/j.aim.2009.12.001.
  * Bezáková et al. [2007] Bezáková, I., Bhatnagar, N., & Vigoda, E. (2007).  Sampling binary contingency tables with a greedy start.  Random Structures & Algorithms, 30, 168–205. doi:10.1002/rsa.20155.
  * Breiger [1974] Breiger, R. L. (1974).  The duality of persons and groups.  Social forces, 53, 181–190. doi:10.1093/sf/53.2.181.
  * Carstens [2015] Carstens, C. J. (2015).  Proof of uniform sampling of binary matrices with fixed row sums and column sums for the fast curveball algorithm.  Physical Review E, 91, 042812. doi:10.1103/PhysRevE.91.042812.
  * Chen et al. [2005] Chen, Y., Diaconis, P., Holmes, S. P., & Liu, J. S. (2005).  Sequential monte carlo methods for statistical analysis of tables.  Journal of the American Statistical Association, 100, 109–120. doi:10.1198/016214504000001303.
  * Chen et al. [2006] Chen, Y., Dinwoodie, I. H., & Sullivant, S. (2006).  Sequential importance sampling for multiway tables.  The Annals of Statistics, 34, 523–545. doi:10.1214/009053605000000822.
  * Ellison et al. [2024] Ellison, D., Jouve, B., & Stone, L. (2024).  Switching checkerboards in (0, 1)-matrices.  Linear Algebra and its Applications, 680, 274–292. doi:10.1016/j.laa.2023.10.008.
  * Fout et al. [2020] Fout, A., Fosdick, B. K., & Hitt, M. P. (2020).  Non-uniform sampling of fixed margin binary matrices.  In Proceedings of the 2020 ACM-IMS on Foundations of Data Science Conference (pp. 95–105).  doi:10.1145/3412815.3416887.
  * Fowler [2006] Fowler, J. H. (2006).  Legislative cosponsorship networks in the us house and senate.  Social networks, 28, 454–465. doi:10.1016/j.socnet.2005.11.003.
  * Godard & Neal [2022] Godard, K., & Neal, Z. P. (2022).  fastball: a fast algorithm to randomly sample bipartite graphs with fixed degree sequences.  Journal of Complex Networks, 10, cnac049. doi:10.1093/comnet/cnac049.
  * Gotelli [2000] Gotelli, N. J. (2000).  Null model analysis of species co-occurrence patterns.  Ecology, 81, 2606–2621. doi:10.1890/0012-9658(2000)081[2606:NMAOSC]2.0.CO;2.
  * Homans [1950] Homans, G. C. (1950).  The human group.  Harcourt, Brace, & World.
  * Jasny [2012] Jasny, L. (2012).  Baseline models for two-mode social network data.  Policy Studies Journal, 40, 458–491. doi:10.1111/j.1541-0072.2012.00461.x.
  * Kolmogorov [1933] Kolmogorov, A. N. (1933).  Sulla determinazione empirica di una legge didistribuzione.  Giorn Dell’inst Ital Degli Att, 4, 89–91.
  * Latapy et al. [2008] Latapy, M., Magnien, C., & Del Vecchio, N. (2008).  Basic notions for the analysis of large two-mode networks.  Social networks, 30, 31–48. doi:10.1016/j.socnet.2007.04.006.
  * Miller & Harrison [2013] Miller, J. W., & Harrison, M. T. (2013).  Exact sampling and counting for fixed-margin matrices.  Annals of Statistics, 41, 1569–1592. doi:10.1214/13-AOS1131.
  * Neal [2014] Neal, Z. P. (2014).  The backbone of bipartite projections: Inferring relationships from co-authorship, co-sponsorship, co-attendance and other co-behaviors.  Social Networks, 39, 84–97. doi:10.1016/j.socnet.2014.06.001.
  * Neal [2022] Neal, Z. P. (2022).  Constructing legislative networks in r using incidentally and backbone.  Connections, 42, 1–9. doi:10.2478/connections-2019.026.
  * Neal et al. [2021] Neal, Z. P., Domagalski, R., & Sagan, B. (2021).  Comparing alternatives to the fixed degree sequence model for extracting the backbone of bipartite projections.  Scientific Reports, . doi:10.1038/s41598-021-03238-3.
  * Newman [2001] Newman, M. E. (2001).  Scientific collaboration networks. i. network construction and fundamental results.  Physical review E, 64, 016131. doi:10.1103/PhysRevE.64.016131.
  * Sanderson [2000] Sanderson, J. G. (2000).  Testing ecological patterns.  American Scientist, 88, 332. doi:10.1511/2000.29.332.
  * Smirnov [1948] Smirnov, N. (1948).  Table for estimating the goodness of fit of empirical distributions.  The annals of mathematical statistics, 19, 279–281.
  * Strona & Fattorini [2014] Strona, G., & Fattorini, S. (2014).  On the methods to assess significance in nestedness analyses.  Theory in Biosciences, 133, 179–186. doi:10.1007/s12064-014-0203-1.
  * Strona et al. [2014] Strona, G., Nappo, D., Boccacci, F., Fattorini, S., & San-Miguel-Ayanz, J. (2014).  A fast and unbiased procedure to randomize ecological binary matrices with fixed row and column totals.  Nature Communications, 5, 1–9. doi:10.1038/ncomms5114.
  * Verhelst [2008] Verhelst, N. D. (2008).  An efficient MCMC algorithm to sample binary matrices with fixed marginals.  Psychometrika, 73, 705–728. doi:10.1007/s11336-008-9062-3.
  * Wang [2020] Wang, G. (2020).  A fast mcmc algorithm for the uniform sampling of binary matrices with fixed margins.  Electronic Journal of Statistics, 14, 1690–1706. doi:10.1214/20-EJS1702.


# The master relation for polynomiality and equivalences of integrable systems

Xavier Blot  X. B.: Korteweg-de Vriesinstituut voor Wiskunde, Universiteit van
Amsterdam, Postbus 94248, 1090GE Amsterdam, Nederland<EMAIL_ADDRESS>,
Adrien Sauvaget  A. S.: CNRS, Université de Cergy-Pontoise, Laboratoire AGM,
UMR 8088, 2 av. Adolphe Chauvin 95302, Cergy-Pontoise Cedex, France
<EMAIL_ADDRESS>and  Sergey Shadrin  S. S.: Korteweg-de
Vriesinstituut voor Wiskunde, Universiteit van Amsterdam, Postbus 94248,
1090GE Amsterdam, Nederland<EMAIL_ADDRESS>
###### Abstract.

We prove the so-called master relation in the tautological ring of the moduli
space of curves that implies polynomial properties of the Dubrovin–Zhang
hierarchies associated to different versions of cohomological field theories
as well as their equivalences to the corresponding double ramification
hierarchies.

##  1\. Introduction

The main goal of the paper is to prove the so-called master relation in
$R^{*}({\overline{\mathcal{M}}}_{g,n+m})$, $g\geq 0$, $n,m\geq 1$,
$2g-2+n+m>0$, conjectured in [3, Conjecture 3.4]. This master relation is
proven in _op. cit._ in the Gorenstein quotient and used to imply the strong
form of the DR/DZ equivalence conjecture due to Buryak-Dubrovin-Guéré-Rossi
[5, 6] (see also the foundational work of Buryak [4] for an earlier weaker
version of the DR/DZ equivalence conjecture).

The master relation in the Gorenstein quotient is sufficient for the most
important applications to integrable systems associated with ordinary
Cohomological field theories (CohFTs), partial CohFTs, (ordinary or partial),
and F-CohFTs. Indeed, the examples considered in the literature in the context
of integrable systems are constructed in terms of the tautological classes.
Yet, it was believed for a long time that polynomiality and equivalence
properties of integrable systems of topological type are governed by actual
relations in the tautological ring of the moduli space of curves. To this end,
a number of conjectures were made, see [7, Conjecture 2.5] and [8, Conjectures
1, 2, 3]. All these conjectures were shown through a sequence of works [8, 1,
3, 2] to follow from the master relation. Therefore, the present paper settles
these conjectures as well.

Despite the fact that the master relation that we prove in this paper is of
crucial importance for integrable systems, we won’t survey these applications
as it is already done in a very detailed way in [8, Section 4]. In order to
keep this paper short, we focus on the statement and proof of the master
relation itself. It is worth noticing that our proof follows from a direct
application of the virtual localization formula for moduli spaces of
(relative) stable maps to $\mathbb{P}^{1}$ [11], while the previous proof in
the Gorenstein quotient relied on a series of reduction steps to control the
intersection of the master relation with generators of the tautological rings.

###  1.1. Acknowledgments

X. B. and S. S. are supported by the Netherlands Organization for Scientific
Research.

##  2\. The master relation

###  2.1. Definition

We follow the exposition in [3].

####  2.1.1. Pre-stable star rooted trees

We fix integers $m\geq 1$, $n\geq 1$, and $g\geq 0$, such that $2g-2+n+m>0$.
We denote by $\mathrm{PSSRT}_{g,n,m}$ the set of $n+m$-pre-stable star rooted
trees of genus $g$, that is the set of pre-stable graphs

$T=(V,H,\iota\colon H\to H,g\colon V\to\mathbb{Z}_{\geq
0},H^{\iota}\simeq\\{\sigma_{1},\ldots,\sigma_{n+m}\\})$

of genus $g$ with $n+m$ legs, satisfying the following constraints:

  * •

If $E(T)$ stands for the set of edges $H^{\iota}$, then the graph $(V,E)$ is a
rooted tree, and all edges are between the root and another vertex (hence the
term “star”).

  * •

The legs $\sigma_{n+1},\dots,\sigma_{n+m}$ are attached to the root vertex,
while $\sigma_{1},\dots,\sigma_{n}$ are attached to non-root vertices. In
particular, the root vertex $v_{r}$ is uniquely determined by the rest of the
data.

  * •

For each $v\in V(T)$, let $H(v)$ be the set of half-edges incident to $v$
(including the legs). The pre-stability condition means that
$2g(v)-2+|H(v)|\geq 0$ for every $v\in V(T)$.

  * •

There is at least one leg attached to each vertex.

We will consider the ring $Q=\mathbb{Q}(a_{1},\ldots,a_{n})$ (where the
$a_{i}$’s are formal variable). If $T$ is a graph in $\mathrm{PSSRT}_{g,n,m}$,
then each edge $e$ connects the root vertex $v_{r}$ to a non-root vertex that
we denote $v_{e}$. We set $a(e)$ to be the sum of $a_{i}$’s associated with
the legs attached to $v_{e}$. If a half-edge $h$ is the leg $\sigma_{i}$, then
$a(h)$ stands for $a_{i}$.

####  2.1.2. Classes assigned to a tree

We want to assign to $T$ a class $\Xi(T)$ in
$R^{*}({\overline{\mathcal{M}}}_{g,n+m})\otimes_{\mathbb{Q}}Q[u,u^{-1}]$
(where $u$ is an extra formal variable used to control the degree). In order
to construct this class, we first assign a class to each vertex.

Let $v$ be a vertex of $T$. We will consider the moduli space of curves
${\overline{\mathcal{M}}}_{g(v),|H(v)|}$. If $v$ is the root vertex $v_{r}$ we
have a natural isomorphism $H(v_{r})\cong
E(T)\sqcup\\{\sigma_{n+1},\dots,\sigma_{n+m}\\}$ and we associate the first
$|E(T)|$ marked points in the corresponding space
${\overline{\mathcal{M}}}_{g(v_{r}),|H(v_{r})|}$ with the edges, and the
remaining $m$ points with the legs $\\{\sigma_{n+1},\dots,\sigma_{n+m}\\}$. We
set

(2.1) $\displaystyle\Psi(v_{r})$ $\displaystyle\coloneqq\prod_{e\in
E(T)}\frac{1}{1-a({e})\psi_{e}}\in
R^{*}({\overline{\mathcal{M}}}_{g(v_{r}),|E(T)|+m})\otimes_{\mathbb{Q}}Q$

where $\psi_{e}$ stands for the $\psi$-class attach the half-edge part of the
edge $e$. In the exceptional unstable case $g(v_{r})=0,m=1,|E(T)|=1$, we
formally assign to the root vertex the following class:

(2.2) $\displaystyle\Psi(v_{r})\coloneqq a(e)^{-1}\in
R^{-1}({\overline{\mathcal{M}}}_{0,2})\otimes_{\mathbb{Q}}Q,$

where $R^{-1}({\overline{\mathcal{M}}}_{0,2})$ is identified with $\mathbb{Q}$
and the negative cohomological degree is just formally assigned to allow to
treat this case non-exceptionally in what follows.

If $v$ is a non-root vertex $v_{e}\not=v_{r}$ for some edge $e$, the we
associate the first $|H(v_{e})|-1$ marked points in the space
${\overline{\mathcal{M}}}_{g(v_{e}),|H(v_{e})|}$ with the legs attached to
$v_{e}$, and the last marked point with the edge $e$ that connects $v_{e}$ to
$v_{r}$. We set

(2.3) $\displaystyle{\mathbb{D}}(v_{e})$
$\displaystyle\coloneqq\frac{\lambda_{g(v_{e})}\mathrm{DR}_{g(v_{e})}\big{(}a(h_{1}),\dots,a(h_{|H(v_{e})|-1}),-a(e)\big{)}}{1-a(e)\psi_{|H(v_{e})|}}\in
R^{*}({\overline{\mathcal{M}}}_{g(v_{e}),|H(v_{e})|})\otimes_{\mathbb{Q}}Q.$

As $a(e)=a(h_{1})+\cdots+a(h_{|H(v_{e})|-1})$, the class
$\mathrm{DR}_{g(v_{e})}\big{(}a(h_{1}),\dots,a(h_{|H(v_{e})|-1}),-a(e)\big{)}$
is the corresponding double ramification cycle. In the exceptional unstable
case $g(v_{e})=0,H(v_{e})=2$ we formally assign to this vertex the following
class:

(2.4) $\displaystyle{\mathbb{D}}(v_{e})\coloneqq a(e)^{-1}\in
R^{-1}({\overline{\mathcal{M}}}_{0,2})\otimes_{\mathbb{Q}}Q,$

where, as it was for the root vertex, the negative cohomological degree is
just formally assigned to allow the treatment of this case non-exceptionally.
With this system of notation, we set

(2.5) $\displaystyle\Xi(T)\coloneqq u^{2g-2+m}\left(\prod_{e\in
E(T)}\frac{a(e)}{u}\right)(b_{T})_{*}\left(\left(\sum_{d=-1}^{\infty}\frac{\Psi(v_{r})_{d}}{(-u)^{d}}\right)\otimes\bigotimes_{e\in
E(T)}\left(\sum_{d=-1}^{\infty}\frac{{\mathbb{D}}(v_{e})_{d}}{u^{d}}\right)\right),$

where $(b_{T})_{*}\colon\bigotimes_{v\in
V(T)}R^{*}({\overline{\mathcal{M}}}_{g(v),|H(v)|})\to
R^{*}({\overline{\mathcal{M}}}_{g,n+m})$ is the boundary push-forward
(composed with the contraction of unstable components) map extended by
linearity with the coefficients in $Q[u,u^{-1}]$, and $\alpha_{d}$ stands for
the projection to $R^{*}({\overline{\mathcal{M}}}_{g(v),|H(v)|})$ of a
tautological class.Note that the resulting formula is a polynomial in
$a_{i}$’s of degree bounded by $3g-3+n+m$.

###  2.2. Statement and proof

For any $g\geq 0$, $m,n\geq 1$, $2g-2+n+m>0$, we set

(2.6)
$\displaystyle\Xi^{m}_{g,n}\coloneqq\sum_{T\in\mathrm{PSSRT}_{g,n,m}}\Xi(T).$

######  Theorem 2.1.

We have

(2.7) $\displaystyle\Xi^{m}_{g,n}\in
R^{*}({\overline{\mathcal{M}}}_{g,n+m})\otimes_{\mathbb{Q}}\mathbb{Q}[a_{1},\dots,a_{n},u].$

(in other words, the coefficients of all negative degrees of $u$ vanish).

###### Proof.

Since $\Xi^{m}_{g,n}$ is a polynomial in $a_{1},\dots,a_{n}$ of degree at most
$3g-3+m+n$, it is sufficient to prove that $\Xi^{m}_{g,n}(a_{1},\dots,a_{n})$
is a polynomial in $u$ for each $a_{1},\dots,a_{n}\in\mathbb{Z}_{>0}$. Note
also that it is sufficient to consider the cohomology with the complex
coefficients. Then the proof directly generalizes the argument in [8, Proof of
Theorem 5.3] along the lines of computation done [10, Section 3.2].

We consider the moduli space
${\overline{\mathcal{M}}}_{g,n+m}(\mathbb{P}^{1},a_{1},\dots,a_{n})$ of
relative stable maps to $(\mathbb{P}^{1},\infty)$: the first $n$ marked points
are mapped to $\infty$ with contact orders prescribed by the $a_{i}$’s. We
consider the standard $\mathbb{C}^{*}$-action on $\mathbb{P}^{1}$ and induced
action on ${\overline{\mathcal{M}}}_{g,m}(\mathbb{P}^{1},a_{1},\dots,a_{n})$.
Let $\pi\colon
U\to{\overline{\mathcal{M}}}_{g,m}(\mathbb{P}^{1},a_{1},\dots,a_{n})$ be the
universal curve and $f\colon U\to\mathbb{P}^{1}$ the universal map. The curve
$U$ is endowed with the $\mathbb{C}^{*}$-action making $f$ equivariant. For
$i\in\\{n+1,\ldots,n+m\\}$ we denote ${\rm
ev}_{i}\colon{\overline{\mathcal{M}}}_{g,n+m}(\mathbb{P}^{1},a_{1},\dots,a_{n})\to\mathbb{P}^{1}$
the evaluation morphism at the $i$-th marking.

We will consider the line bundle $\mathcal{O}(-1)\to\mathbb{P}^{1}$. We lift
the $\mathbb{C}^{*}$-action to $\mathcal{O}(-1)\to\mathbb{P}^{1}$ with the
fiber weights at $0$ and $\infty$ equal to $-1$ and $0$ respectively. We set

(2.8) $\displaystyle
I^{m}_{g,n}\coloneqq(-1)^{g+m}\mathrm{ev}^{*}_{n+1}([0])\cup\cdots\cup\mathrm{ev}^{*}_{n+m}([0])\cup
e_{\mathbb{C}^{*}}(R^{1}\pi_{*}f^{*}\mathcal{O}(-1))$

in $H_{*}^{\mathbb{C}^{*}}({\overline{\mathcal{M}}}_{g,n+m})\cong
H_{*}({\overline{\mathcal{M}}}_{g,n+m})\otimes_{\mathbb{C}}\mathbb{C}[u]$. In
this expression, the class $[0]\in H^{2}_{\mathbb{C}^{*}}(\mathbb{P}^{1})$
stands for the $\mathbb{C}^{*}$-equivariant class dual to the point $0$, while
$e_{\mathbb{C}^{*}}$ is the $\mathbb{C}^{*}$-equivariant Euler class (here,
the derived pushforward $R^{1}\pi_{*}f^{*}\mathcal{O}(-1)$ is a vector bundle
of rank $g+a_{1}+\cdots+a_{n}-1$).

Let
$\epsilon\colon{\overline{\mathcal{M}}}_{g,m}(\mathbb{P}^{1},a_{1},\dots,a_{n})\to{\overline{\mathcal{M}}}_{g,n+m}$
and
$[{\overline{\mathcal{M}}}_{g,m}(\mathbb{P}^{1},a_{1},\dots,a_{n})]^{\mathrm{vir}}$
denote the $\mathbb{C}^{*}$-equivariant virtual fundamental class of
${\overline{\mathcal{M}}}_{g,m}(\mathbb{P}^{1},a_{1},\dots,a_{n})$. We will
prove that

(2.9)
$\displaystyle\Xi^{m}_{g,n}(a_{1},\dots,a_{n},-u)=\epsilon_{*}\left(I^{m}_{g,n}\cap[{\overline{\mathcal{M}}}_{g,m}(\mathbb{P}^{1},a_{1},\dots,a_{n})]^{\mathrm{vir}}\right).$

To this end, we need to apply the localization formula to (2.9) along the
lines of [11] (see also a survey in [8, Appendix]), and the computation of
$e_{\mathbb{C}^{*}}(R^{1}\pi_{*}f^{*}\mathcal{O}(-1))$ performed in [10,
Equation (3.13)]. As the right-hand side of (2.9) belongs to
$H_{*}^{\mathbb{C}^{*}}({\overline{\mathcal{M}}}_{g,n+m})\cong
H_{*}({\overline{\mathcal{M}}}_{g,n+m})\otimes_{\mathbb{C}}\mathbb{C}[u]$, the
theorem follows.

The localization formula expresses
$[{\overline{\mathcal{M}}}_{g,m}(\mathbb{P}^{1},a_{1},\dots,a_{n})]^{\mathrm{vir}}$
as a sum over bipartite graphs with some vertices corresponding to the
components mapped to $0$ and $\infty$. We allow these vertices to correspond
to the unstable moduli spaces as well — formally speaking, in such cases, a
separate computation is required, but the final result can be written
uniformly with the contributions of the unstable vertices introduced formally
as we did above in the definition of $\Xi^{m}_{g,n}$. Only part of these
bipartite graphs contributes when we intersect the virtual fundamental class
with $I^{m}_{g,n}$:

  * •

The first $n$ legs are incident to vertices over $\infty$ because of the
relative conditions, while the last $m$ legs are over $0$ because of the
factors ${\rm ev}_{i}^{*}([0])$

  * •

Our specific choice of lift of the $\mathbb{C}^{*}$-action to
$\mathcal{O}(-1)$ reduces the computation to an expression on star trees.
Indeed, the expression of
$e_{\mathbb{C}^{*}}(R^{1}\pi_{*}f^{*}\mathcal{O}(-1))$ has a factor
$0^{|E(v)|-1}$ for each vertex $v$ over infinity, where $E(v)$ is the set of
edges attached to $v$ (see [10, Equation (3.13)])111This observation was used
in [12] and [9] to obtain the multiple cover formula for rigid rational curves
in CY 3-folds.. Therefore, each vertex over $\infty$ should have exactly one
edge (otherwise, the class vanishes).

Hence, the bipartite graphs that we consider are, in fact, the graphs in
$\mathrm{PSSRT}_{g,n,m}$: the root is the unique vertex over $0$ while the
other vertices are over $\infty$. Note that the degree assignment at edges,
which is part of the bi-partite graph data, is uniquely determined by the
relative conditions at infinity and the fact that the graphs are of compact
type. Altogether, the virtual localization formula provides an expression of
the form

(2.10)
$\displaystyle\epsilon_{*}\left(I^{m}_{g,n}\cap[{\overline{\mathcal{M}}}_{g,m}(\mathbb{P}^{1},a_{1},\dots,a_{n})]^{\mathrm{vir}}\right)=\sum_{T\in\mathrm{PSSRT}_{g,n,m}}(b_{T})_{*}(\alpha_{1}(T)\alpha_{2}(T))$

where $\alpha_{1}(T)$ corresponds to the contribution of
$[{\overline{\mathcal{M}}}_{g,m}(\mathbb{P}^{1},a_{1},\dots,a_{n})]^{\mathrm{vir}}$
from [11] (see also a convenient reminder in [8, Appendix]), while
$\alpha_{2}(T)$ corresponds to the contribution of $I_{m,n}$ deduced from [10,
Equation (3.13)]. The explicit expression of these classes is

(2.11) $\displaystyle\alpha_{1}(T)=$
$\displaystyle\,-u^{-|E(T)|-1}\left(\sum_{i=0}^{g(v_{r})}(-1)^{i}\lambda_{i}u^{g(v_{r})-i}\prod_{e\in
E(T)}\frac{a(e)^{a(e)+1}}{a(e)!}\frac{u^{-a(e)}}{1-u^{-1}a(e)\psi_{e}}\right)$
$\displaystyle\,\bigotimes_{e\in
E(T)}\left(\frac{\mathrm{DR}_{g(v_{e})}\big{(}a(h_{1}),\dots,a(h_{|H(v_{e})|-1}),-a(e)\big{)}}{1+u^{-1}a(e)\psi_{|H(v_{e})|}}\right),$

and

(2.12) $\displaystyle\alpha_{2}(T)=$
$\displaystyle\,(-u)^{|E(T)|-1}(-1)^{g+m}u^{m}\left((-1)^{g(v_{r})}\sum_{i=0}^{g(v_{r})}\lambda_{i}u^{g(v_{r})-i}\prod_{e\in
E(T)}\frac{a(e)!}{a(e)^{a(e)}}u^{a(e)-1}\right)$
$\displaystyle\,\bigotimes_{e\in E(T)}(-1)^{g(v_{e})}\lambda_{g(v_{e})}.$

Now we use Mumford’s relation [13]

(2.13)
$\displaystyle\left(\sum_{i=0}^{g(v_{r})}(-1)^{i}\lambda_{i}u^{g(v_{r})-i}\right)\left(\sum_{i=0}^{g(v_{r})}\lambda_{i}u^{g(v_{r})-i}\right)=u^{2g(v_{r})}$

to deduce

(2.14) $\displaystyle\alpha_{1}(T)\alpha_{2}(T)=$
$\displaystyle\,(-u)^{2g-2+m}\left(\prod_{e\in
E(T)}\frac{a(e)}{-u}\frac{1}{1-u^{-1}a(e)\psi_{e}}\right)$
$\displaystyle\bigotimes_{e\in
E(T)}\left(\frac{(-u)^{2g(v_{r})}\lambda_{g(v_{r})}\mathrm{DR}_{g(v_{e})}\big{(}a(h_{1}),\dots,a(h_{|H(v_{e})|-1}),-a(e)\big{)}}{1+u^{-1}a(e)\psi_{|H(v_{e})|}}\right).$

Together with (2.10), we indeed obtain the identity (2.9). ∎

## References

  * [1] X. Blot, D. Lewański, P. Rossi, S. Shadrin. Stable tree expressions with Omega-classes and double ramification cycles. arXiv:2403.05190.
  * [2] X. Blot, D. Lewański, S. Shadrin. Rooted trees with level structures, $\Omega$-classes and double ramification cycles. arXiv:2406.06205.
  * [3] X. Blot, D. Lewański, S. Shadrin. On the DR/DZ equivalence conjecture. arXiv:2405.12334.
  * [4] A. Buryak. Double ramification cycles and integrable hierarchies. Comm. Math. Phys. 336 (2015), no. 3, 1085-1107.
  * [5] A. Buryak, B. Dubrovin, J. Guéré, P. Rossi. Tau-structure for the double ramification hierarchies. Comm. Math. Phys. 363 (2018), no. 1, 191-260.
  * [6] A. Buryak, B. Dubrovin, J. Guéré, P. Rossi. Integrable systems of double ramification type. Int. Math. Res. Not. 2020 (2020), no. 24, 10381-10446.
  * [7] A. Buryak, J. Guéré, P. Rossi. DR/DZ equivalence conjecture and tautological relations. Geom. Topol. 23 (2019), no. 7, 3537-3600.
  * [8] A. Buryak, S. Shadrin. Tautological relations and integrable systems. arXiv:2210.07552.
  * [9] C. Faber, R. Pandharipande. Hodge integrals and Gromov-Witten theory. Invent. Math. 139 (2000), no. 1, 173-199.
  * [10] A. Giacchetto, R. Kramer, D. Lewański, A. Sauvaget. The Spin Gromov-Witten/Hurwitz correspondence for $\mathbb{P}^{1}$. arXiv:2208.03259.
  * [11] T. Graber, R. Vakil. Relative virtual localization and vanishing of tautological classes on moduli spaces of curves. Duke Math. J. 130 (2005), no. 1, 1-37.
  * [12] Yu. I. Manin. Generating functions in algebraic geometry and sums over trees. The moduli space of curves (Texel Island, 1994), 401-417, Progr. Math., 129, Birkhäuser Boston, Boston, MA, 1995.
  * [13] D. Mumford. Towards an enumerative geometry of the moduli space of curves. Arithmetic and geometry, Vol. II, 271-328, Progr. Math., 36, Birkhäuser Boston, Boston, MA, 1983.


# A derivation of the electric field inside MAPS detectors from beam-test data
and limited TCAD simulations

A<EMAIL_ADDRESS>N. Tal<EMAIL_ADDRESS>Department of Particle Physics and Astrophysics,Weizmann Institute of
Science,Rehovot,7610001, Israel

###### Abstract

Solid semiconductor sensors are used as detectors in high-energy physics
experiments, in medical applications, in space missions and elsewhere. A
precise knowledge of the electric field inside the basic cells of these
sensors is highly important for their characterization and performance
understanding. The field governs the charge propagation processes and
ultimately determines the size and quality of the electronic signal of the
cell. Hence, the simulation of these sensors relies strongly on the electric
field knowledge. For a certain voltage applied to the cell, the field depends
on the specifics of the device’s growth and fabrication. The information about
these is often commercially protected or otherwise very difficult to encode in
state-of-the-art technology computer-aided-design (TCAD) software. It is
therefore practically impossible to obtain the field without some minimal
knowledge. In this work, we show that combining public beam-test data and a
very limited public TCAD-based knowledge, we are able to effectively
reconstruct the 3D electric field function in the pixel cell of one important
and widely used example, namely the ALPIDE sensor, a monolithic active pixel
sensor (MAPS). Despite its broad usage worldwide, the ALPIDE field is not
available to the community as it is under proprietary restriction. We provide
the 3D effective field function of the ALPIDE sensor and comment on the
process by which it is derived with the help of the Allpix2 software. We also
comment on how starting from the same grounds, similar work can be performed
for other devices.

###### keywords:

ALPIDE, Pixel, TCAD, $E$-Field, Charge, Allpix2

††journal: Nucl. Instr. Meth. A

##  1 Introduction

Solid semiconductor sensors are used for imaging, tracking and calorimetry in
different disciplines. These detectors may have different level of tolerance
for damage caused by different radiation types and doses, depending on the
application. When an ionizing particle traverses the elementary cell of the
sensor (pixel, strip, pad, etc.), an electric current pulse is produced at the
collection node of that cell (the electronic signal). This current is in
effect produced from electrons or holes (charge carriers), which are created
in pairs along the trajectory of the ionizing traversing particle. The
collected charges are processed by the front-end electronics and a “hit” is
formed if the collected charge is above certain threshold. The number of
charge carriers and the way these propagate to the collection node depend
therefore on the specifics of the sensor material. The electronic signal
depends in turn strongly on the material, but it also depends on a number of
other parameters such as the sensor fabrication processes, the sensor
geometry, the external voltage applied to the sensor, the damage the sensor
may have suffered due to continuous exposure to radiation, etc.

A key ingredient in the characterization and understanding charge carriers’
propagation in the cell is the electric field inside it. The field governs
this propagation process and ultimately determines the size and quality of the
electronic signal output of the cell. Hence, the simulation of these sensors
relies strongly on the electric field knowledge. An important task in the
overall design process of a sensor is the implementation of the corresponding
technology computer-aided-design (TCAD) simulation of the device. This can be
done using “Sentaurus” software from Synopsys [1] or “Atlas” software from
Silvaco [2]. One of the outputs of such TCAD simulation is the fine-grained 3D
electric field mesh inside the cell, which is later used to simulate the
charge propagation. Hence, the field is an essential piece in the
characterization of the device performance. Implementing the device’s growth
and fabrication procedures in TCAD is an extremely complicated task. Moreover,
the information about these procedures is usually protected commercially. In
some cases, the field may be one dimensional, e.g., $\vec{E}=E_{z}$ (with $z$
being the coordinate along the depth of the cell). Furthermore, this single
component may be modeled to a good approximation as a simple function of $z$
only. In these cases, the absence of a detailed TCAD mesh description of the
field is not a limiting factor for the simulation chain of the device.
However, in other cases, the field may be highly non-linear and generally, it
may have non-trivial components also in the transverse ($x,y$) coordinates,
i.e., $\vec{E}=(E_{x},E_{y},E_{z}$). The three components can themselves
depend non-trivially on all three coordinates, i.e., $E_{i}=f_{i}(x,y,z)$,
where $i=x,y,z$ and $f_{i}$ is some 3D function that is not necessarily
analytical.

This is particularly important for new monolithic active pixel sensors (MAPS)
[3], which are used more and more extensively in recent years thanks to the
development of sub-micron capabilities in different foundries worldwide.
Unlike in hybrid sensors, the readout circuitry of MAPS is integrated in the
same substrate as the sensing volume. Consequently, the material budget may be
substantially reduced, the production costs can also be reduced and the
production yield can be higher. MAPS are now becoming popular, e.g., for the
next-generation high-energy physics experiments (for tracking, vertexing and
calorimetry) and elsewhere.

In some of these cases it may happen that users of a specific sensor do not
have access to the TCAD mesh of the non-trivial field. This poses a serious
problem for these users, who effectively lack the first building block of the
device simulation and therefore cannot properly simulate its response. To
overcome this limitation, we show that combining public beam-test data and a
very limited public TCAD-based knowledge, we are able to effectively
reconstruct the 3D electric field function in the pixel cell of one important
and widely used example, namely the ALPIDE monolithic silicon sensor [4, 5, 6,
7, 8]. The ALPIDE (ALice PIxel DEtector) sensor was developed by and for the
ALICE [9] experiment at the LHC together with TowerJazz [10]. The inner
tracking system (ITS) of ALICE [11] has been upgraded successfully in 2021
with this new technology and it is now the largest pixel detector ever built.
The successful ALPIDE technology is also being used elsewhere in other
experiments/facilities like LUXE [12], sPHENIX [13, 14], DESY’s beam-test
facility [15], the CSES-2 satellite mission [16]) in space, Proton Computed
Tomography (pCT) [17, 18], and more. However, despite its wide usage, the TCAD
sampling points (referred to hereafter as “mesh”) are not available to the
community as it is under proprietary restriction.

In this work, we discuss how the 3D effective field function (denoted EFF
hereafter) in the ALPIDE pixel sensor is derived. We use the Allpix2 software
[19, 20] to compare the performance in simulation between two cases, once
starting from our 3D EFF and once starting from public results that are using
the actual TCAD mesh, which remains unknown to us. The two cases are denoted
hereafter as Allpix2+EFF and Allpix2+TCAD, respectively. For the comparison,
we take help of the Allpix2+TCAD results shown in [21]. The authors of [21] go
further to compare their Allpix2+TCAD simulation with data collected in a
beam-test [22]. The results in [21] and [22] are obtained using a prototype of
the ALPIDE sensor, called “Investigator”. The pixel design of the Investigator
sensor is close to that of the ALPIDE design, except for the pitch (the
Investigator has square pixel pitch of different sizes).

The performance of the three scenarios is shown to be very similar. This
agreement gives confidence in our EFF as the cornerstone for further
simulation campaigns, where the ALPIDE sensor is used. We comment on the
derivation process and provide the EFF as a set of ROOT [23] TFormula-
compatible strings. We provide the code that produces our EFF such that
further fine-tuning may be done, and such that similar work can be started for
other sensors. We also provide the Allpix2 configuration files used in our
Allpix2+EFF simulation. Our EFF generation code and Allpix2 configurations are
available in [24].

##  2 The Allpix2 setup

The Allpix2 software [19, 20] is a lightweight C++ framework that simulates
the processes triggered inside different semiconductor devices when ionizing
particles traverse these. Particularly, it simulates the electron-hole pair
creation, the drift or diffusion of these charge carriers given some electric
field, the charge collection by the electrodes and finally the digitization of
collected charge along with the electronic noise in the front-end electronics.
The task of simulating the energy deposition process by the ionizing particles
(an input for the charge carriers generation) is external to Allpix2. Usually
this is done by Geant4 [25, 26, 27], which is interfaced directly to Allpix2
for convenience, but this can also be provided from other sources. The Allpix2
simulation flow is implemented using several modules, accessed sequentially
using a configuration file. While for a full description the reader is
referred to the Allpix2 user manual [20], we briefly summarize the most
relevant modules related to the purpose of this work:

  * 1.

GeometryBuilderGeant4: constructs the Geant4 geometry given the description of
the sensors and their arrangement (dimension, position and orientation).

  * 2.

DepositionGeant4: implements the simulation of energy deposition by ionizing
particles in the active volume of the detector.

  * 3.

ElectricFieldReader: implements the electric field inside the sensor volume.
The field can be applied using several methods. These include TCAD mesh input
files, simple built-in functions (constant, linear, parabolic), and most
importantly for this work, any user-defined analytic functions in the form of
ROOT TFormula-compatible strings.

  * 4.

GenericPropagation: simulates the propagation of electrons and holes in the
sensor volume given the electric field.

  * 5.

DefaultDigitizer: translates the collected charge into a digitized signal
proportional to the input charge. The noise in the electronics is also
simulated here.

These modules are called sequentially in the simulation chain. The output of
the chain contains the fired pixels, i.e., pixels which have a collected
charge passing some predefined threshold. It is important to stress that an
induced charge generation in adjacent pixels (to the pixel where the ionizing
particle passes through) is also simulated by Allpix2. The list of fired
pixels can then be used in a subsequent analysis, which depends on the
specific application. This analysis may include, e.g., calibrating the charge
to the energy-loss, counting of hits, clustering of the pixels, tracking, etc.

The Allpix2 simulation configuration for this work is kept as close as
possible to that of [21]. This is necessary for the validation of the EFF by
comparing the Allpix2+EFF and Allpix2+TCAD performances. Therefore, the total
thickness of the sensor in our Allpix2+EFF runs is taken to be 100 $\mu$m,
with a depletion depth of 25 $\mu$m, like the sensor used in [21]. The pixel
pitch of the Investigator sensor used is $28\times 28~{}\mu{\rm m}^{2}$ in
$x\times y$. While we configure our code that generates the 3D EFF to have the
same parameters as for the Investigator sensor, we make sure that it is easy
to reconfigure these parameters using the production version of the ALPIDE.
For example, the pixel pitch is defined as a global parameter in the code and
hence, it can be simply re-set at one place such that the EFF would scale
naturally.

A simulated $\pi^{+}$ beam particles with an energy of 120 GeV is set to hit
the detector. This is similar to the beam condition of the beam-test setup of
[22]. The beam axis is defined as the positive $z$ direction, normal to the
sensor plane. After the $\pi^{+}$ beam particles deposit energy and generate
the electron-hole pairs along their trajectory in the sensor, the charge
carriers are propagated by Allpix2 towards the electrodes along the electric
field lines. In the upper part of the sensor (within depletion depth of 25
$\mu$m), both drift and diffusion are the preferred modes of charge carrier
motion, whereas in the region beyond the depletion depth, only diffusion
occurs due to the lack of electric field in that region. The charges reaching
the electrodes within a time window of 22.5 ns (as in [21]) are collected.
Finally, the collected charge is converted into a digital signal. Equivalent
electronic noise is added to the signal by randomly drawing from a Gaussian
distribution with a width of 10 e and mean of 0 e. Different charge collection
thresholds are checked in the range between 40 e and 700 e. A threshold
smearing option is also added to the simulation by sampling from a Gaussian
distribution with width 5 e and mean 0 e. While we mention some key simulation
parameters explicitly above, there are in fact more parameters required for
the process. All parameters are kept as in [21].

The bias voltage applied to the ALPIDE p-well in [21] and [22] is $-6$ V. A
voltage of $+0.8$ V is applied to the collection electrode. When using the
TCAD mesh, the voltage settings are dictated by that. In the user-defined
function case, however, the settings are completely set by the function’s
normalization. This normalization can be capped very roughly by $\sim-6$ V,
when integrating the EFF along the depletion region, but it strongly depends
on the shape as well.

##  3 The process of deriving the effective field function

As mentioned before, the detailed TCAD mesh of the electric field inside the
ALPIDE sensor is unavailable publicly. However, its magnitude, $\sqrt{\Sigma
E_{i}^{2}}$ (with $i=x,y,z$), in the $x-z$ plane for $y=0$ is shown in figure
4 of [21]. The field lines in this plane are also overlaid in this figure. The
magnitude of the field at the faces of the pixel in three-dimensions are shown
in figure 3 of [21]. As discussed above, due to various restrictions, these
two figures are given without axes scale-labels. The meaning of the colors in
these figures can be loosely interpreted from the field of a similar
“investigator” chip shown in figures 8.4 and 8.5 of an older work from 2018
[28]. This work used the CLIC Tracker Detector (CLICTD), a monolithic chip
[29] which is also produced in a 180 nm imaging CMOS333Complementary metal-
oxide-semiconductor. process on a high-resistivity epitaxial layer like the
ALPIDE. The extensive TCAD studies in [28] show a similar field behavior to
the one of the ALPIDE, where unlike [21], more details are provided, in
particular the scale of the color map. This information was also used to
define our EFF. As it can be seen in the figures mentioned above, the field’s
magnitude shape is highly non-linear and the EFF cannot be easily deciphered
from the available information. At this step, we therefore only have a very
rough idea on how the magnitude of the 3D EFF should look like in a few slices
of the pixel. We also know that besides mimicking the shape visually, its
components’ normalizations should result in a voltage of $\sim-6$ V when
integrated along $z$ in the upper 25 $\mu$m of the sensor.

This information is clearly not enough to derive the field, even not
effectively. However, we do have another few indirect pieces of invaluable
information from [21] detailing the performance of the sensor with the
$\pi^{+}$ particle beam discussed above. This includes the charge
distributions (figure 9), the cluster size distributions (figures 10 and 11)
and the position residual distributions in $x$ (figure 15), all given for a
specific threshold. The behavior for different thresholds can be seen in the
cluster size graphs (figure 13), the spatial resolution graphs in $x$ and $y$
(figure 16) and the efficiency graphs (figure 18). With this, we can
iteratively plug in one EFF ansatz at a time such that

  1. 1.

its 3D magnitude, $\sqrt{\Sigma E_{i}^{2}}$ (with $E_{i}=f_{i}(x,y,z)$ and
$i=x,y,z$), visually resembles the magnitude of the TCAD field shape in the
slices from figures 3-4 of [21] including the field lines,

  2. 2.

it results in an integral $-\int{\vec{E}\cdot{\rm d}\vec{\ell}}\simeq-6$ V,
where $\vec{\ell}$ runs along the negative $z$-axis direction at the top 25
$\mu$m of the sensor in its center ($x=y=0$), and

  3. 3.

it gives a good agreement with the performance figures 9-18 of [21] after
simulating the same scenarios with our Allpix2+EFF setup.

For step 1 above, we add the main features of the field, one by one and plot
the field magnitude in 3D and a few 2D slices, similar to the ones available
from TCAD. This is done initially with simple functional shapes (sphere, arcs,
stripes, Gaussian and exponential shapes, etc.). Upon adding a new feature, we
verify that it merges properly with the existing features in terms of the
relative normalization and the smoothness in the transition regions. In some
cases, adding a new prominent feature may lead to irregularities in the
transition regions. In these cases, one can add new, compensating, features
limited to the relevant regions to remove these irregularities. Regardless of
the field magnitude, in case the field lines are also available from TCAD in
some 2D slices, they must always match in direction. This is adjusted by
changing both the feature’s relative normalization and sign in all three
components.

For step 3 above, we initially simulate 1000 primary $\pi^{+}$ beam particles
to see if the cluster charge and size distributions roughly agree with those
in figures 9-10 of [21] for a nominal threshold of 120 e. If these agree, we
continue to check a wider range of thresholds (40-700 e) with larger
statistics (20,000 primary $\pi^{+}$ beam particles) and check the agreement
with figures 9-18 of [21]. For the comparison, we digitize the figures from
[21] using the ‘WebPlotDigitizer’ software [30] to get the
distributions/graphs and their uncertainties wherever applicable. A good
agreement in that respect is therefore judged according to the ratio between
our Allpix2+EFF and the Allpix2+TCAD simulation taking into account the
uncertainties.

The iterative process of adding field features to the EFF ansatz and testing
the performance (per ansatz) may be stopped when all three criteria above are
satisfied to a large extent, noting that the first two are only qualitative
tests. Specifically for the ALPIDE sensor, more than a hundred of such
iterations were needed in order to arrive at a satisfactory EFF. The final EFF
string generated automatically by our code [24], using a set of parametric
logical expressions, has $\sim 2600$ characters for the $z$ component and
$\sim 5700$ characters for each of the $x$ and $y$ components. This gives a
rough idea about how highly-complex the EFF is and particularly in the
transverse direction.

The procedure discussed above can be regarded as generic as long as there is
at least (i) a very basic knowledge on the field’s shape and (ii) reference
data for the comparison of the sensor performance with simulation.

##  4 The effective field function of the ALPIDE sensor

The results shown in this section are given after the full rigorous procedure
described in section 3. The resulting EFF magnitude of the ALPIDE sensor is
seen in figures 1, 2. These plots should be compared with figures 4 and 3 of
[21]. It can be seen that while the EFF magnitude shape is not identical to
the TCAD mesh magnitude shape, it captures all important features.

The integral $-\int{\vec{E}\cdot{\rm d}\vec{\ell}}$ with $\vec{\ell}$ running
along the negative $z$-axis direction at the top 25 $\mu$m of the sensor in
its center ($x=y=0$) results in a bias voltage of -6.08 V, as expected.

The analytical expression (in the form of a ROOT TFormula-compatible string)
is too complex and long to be included in this manuscript and hence it is
saved in [24]. In addition, the code producing the EFF as well as all the
plots below can be found in [24].

\begin{overpic}[width=201.6317pt]{fig/pixel_3D_v18.1.1.3.1.png}\end{overpic}

\begin{overpic}[width=227.64784pt]{fig/EfieldALPIDE_magnitudes_xzy0_slice_with_arrows.pdf}\end{overpic}

Figure 1: Left: the 3D EFF magnitude, $\sqrt{\Sigma E_{i}^{2}}$ (with
$E_{i}=f_{i}(x,y,z)$ and $i=x,y,z$), at the ALPIDE pixel sensor faces (sides
and top). Right: the EFF magnitude in the $z$ vs $x$ plane sliced at $y=0$
overlaid with the field lines shown as white arrows (the arrows are positioned
at the bin centers). All shapes are given for the upper 25 $\mu$m of the
sensor, where the EFF is non-zero. These plots should be compared with figures
4 and 3 of [21], respectively.
\begin{overpic}[width=429.28616pt]{fig/EfieldALPIDE_magnitudes_pixel_faces.pdf}\end{overpic}
Figure 2: The 3D EFF magnitude, $\sqrt{\Sigma E_{i}^{2}}$ (with
$E_{i}=f_{i}(x,y,z)$ and $i=x,y,z$), at the ALPIDE pixel faces overlaid with
the field lines (clockwise from top left: the bottom face of the pixel, its
top face, its $y$ face and its $x$ face). All shapes are given for the upper
25 $\mu$m of the pixel where the EFF is non-zero.

##  5 Sensor performance with the effective field function

The Allpix2 results shown in this section are using the EFF discussed in
section 4. The Allpix2+EFF results are compared with the ones of the
Allpix2+TCAD and the beam-test data from [21]. The comparison is done at the
level of a pixels-cluster in terms of the cluster’s charge, size and position.
The clustering of fired pixels (pixels with charge above threshold) follow a
simple “packman” algorithm, where all adjacent fired pixels (in both $x$ and
$y$ directions) are added together to form a cluster. In this simple
algorithm, no pixel sharing is allowed between two clusters. The same
algorithm is used also in [21]. The cluster size is simply the number of
pixels associated with it. Likewise, the cluster charge is the sum of all
pixels’ charges associated with it. The cluster position in $x-y$ is obtained
as a charge-weighted mean from the positions of all pixels of the cluster.
This can be compared with the truth position of the incoming $\pi^{+}$ beam
particles to obtain the residuals. Finally, the detection efficiency can be
defined in terms of the fraction of truth-matched clusters (i.e. clusters
matched with true incident $\pi^{+}$ beam particles) from all clusters formed.
A successful matching is achieved when the maximum distance between the
cluster position and the true incident particle position is smaller than three
times the pixel pitch in each dimension. For the comparison, we use a wide
range of thresholds, between 40 e and 700 e, where the nominal threshold is of
120 e. The number of primary $\pi^{+}$ beam particles used for the Allpix2+EFF
results below is 20,000 for the nominal as well as the other thresholds. This
gives a low enough statistical error which allows to clearly see the main
trends.

Whenever the statistical uncertainties are available in the Allpix2+TCAD
results of [21], we use those in the ratio plots shown in the bottom panels.
Otherwise, the errors in the ratio panel represent the statistical uncertainty
of our Allpix2+EFF simulation alone. Furthermore, the systematic uncertainty
bands for the Allpix2+TCAD plots shown in [21] are not overlaid here because
it was not possible to decipher these from the source. However, we do comment
qualitatively on the compatibility between the Allpix2+TCAD and Allpix2+EFF
results with respect to the systematic uncertainties from [21]. These
uncertainties should be identical between the two simulations, leading to an
overall uncertainty higher by a factor of $\sqrt{2}$ for the comparison. This
statement is true as long as one excludes systematic variations of the field
function itself, a study which we leave for future work.

The cluster charge and size distributions for the nominal threshold are shown
in figure 3. The most probable value (MPV) resulting from a fit of the
Allpix2+EFF charge distribution to a convolution of a Landau and a Gaussian
functions is $1.48$, whereas the MPV result for the Allpix2+TCAD is $1.42$. A
good agreement is also seen in the two ratio plots.

\begin{overpic}[width=212.47617pt]{fig/ClusterChargeComparisonRatio.pdf}\end{overpic}\begin{overpic}[width=212.47617pt]{fig/ClusterSizeComparisonRatio.pdf}\end{overpic}

Figure 3: The comparison of the cluster charge (left) and cluster size (right)
distributions between our Allpix2+EFF simulation and the Allpix2+TCAD
simulation (and the beam-test data) from [21]. The bottom panels show the
ratio of the Allpix2+TCAD (red) or data (grey) results to the Allpix2+EFF
results. The cluster charge distribution for the data either does not have
statistical uncertainties at the source or the errors are too small for the
digitizer software [30] to decipher. The systematic uncertainty band for the
Allpix2+TCAD charge plot shown in [21] is not overlaid here because it was not
possible to decipher it from the original plot. However, as one can see in
[21], the uncertainty on the Allpix2+TCAD histogram is as high as 1-10% in the
region where the statistics are high enough. The Allpix2+EFF uncertainties
should be identical to that (giving an overall uncertainty higher by a factor
of $\sqrt{2}$).

The distributions of the cluster size in $x$ and $y$, and the residuals in $x$
are checked as well for the nominal threshold in figures 4 and 5. A similar
conclusion about the good agreement can be drawn in these cases as well.
Particularly, using figure 5 and following the exact procedure given in [21]
to calculate the resolution, we find that it is $\sigma_{x}=3.34$ $\mu$m for
Allpix2+EFF at a threshold of 120 e, compared with $3.60\pm 0.01({\rm
stat})^{+0.24}_{-0.13}({\rm syst})$ $\mu$m for Allpix2+TCAD. As mentioned
earlier, since the procedures are identical, the overall uncertainty for our
comparison will be a factor of $\sqrt{2}$ higher.

\begin{overpic}[width=212.47617pt]{fig/ClusterSizeXComparisonRatio.pdf}\end{overpic}\begin{overpic}[width=212.47617pt]{fig/ClusterSizeYComparisonRatio.pdf}\end{overpic}

Figure 4: The comparison of the cluster size in $x$ (left) and in $y$ (right)
distributions between our Allpix2+EFF simulation and the Allpix2+TCAD
simulation (and the beam-test data) from [21]. The bottom panels show the
ratio of the Allpix2+TCAD (red) or data (grey) results to the Allpix2+EFF
results. The data histograms either do not have statistical uncertainties at
the source or the errors are too small for the digitizer software [30]. The
source also does not have systematic uncertainties for these plots.
\begin{overpic}[width=390.25534pt]{fig/ResidualXComparisonRatio.pdf}\end{overpic}
Figure 5: The comparison of the cluster position residuals distribution in $x$
between our Allpix2+EFF simulation and the Allpix2+TCAD simulation (and the
beam-test data) from [21]. The bottom panel show the ratio of the Allpix2+TCAD
(red) or data (grey) results to the Allpix2+EFF result. The three resolutions
reported are calculated in the same way discussed in [21]. The histograms from
[21] either do not have statistical uncertainties at the source or the errors
are too small for the digitizer software [30]. The systematic uncertainty band
for the Allpix2+TCAD shown in [21] is not overlaid here because it was not
possible to decipher it from the original plot.

Finally, we also see a good agreement when scanning the thresholds in the
range mentioned above and plotting the mean cluster size, the detection
efficiency and the spatial resolution in $x$ and $y$ vs the threshold. These
results are shown in figures 6 and 7. We take the systematic uncertainties of
the Allpix2+TCAD results from [21] qualitatively and note again that these
should be identical to those of our Allpix2+EFF simulation.

The largest discrepancy between the two results is seen in the spatial
resolution at low thresholds, $\lesssim 150$ e. This discrepancy is
nevertheless within the uncertainty for the resolution in $y$ (5%-15% in
40-120 e), while for $x$, the two results are compatible within the
uncertainties ($\sim 7\%$ in 40-200 e) only above $\simeq 80$ e. This lies
well below the nominal projected working range of the sensors. The second-
largest discrepancy is evident in the efficiency graph at the highest
thresholds, $\gtrsim 600$ e. This discrepancy is well within the systematic
uncertainty ($\sim 1-2\%$ above 600 e) and also here, it is manifested in a
range that is much above the nominal projected working range of the sensors.

The performance plots given here show a comprehensive set of tests in a wide
range of values each, indicating a good compatibility between our Allpix2+EFF
results and the Allpix2+TCAD results and particularly for the relevant range
of operation. These two sets of simulation results also show the same level of
compatibility with the data.

\begin{overpic}[width=212.47617pt]{fig/ClusterMeanRatio.pdf}\end{overpic}\begin{overpic}[width=212.47617pt]{fig/EfficiencyRatio.pdf}\end{overpic}

Figure 6: Left: the comparison of the mean cluster size behavior for different
thresholds between our Allpix2+EFF simulation and the Allpix2+TCAD simulation
(and the beam-test data) from [21]. Right: the same comparison for the
detection efficiency. The bottom panels show the ratio of the Allpix2+TCAD
(red) or data (black) results to the Allpix2+EFF results. The graphs from [21]
either do not have statistical uncertainties at the source or the errors are
too small for the digitizer software [30]. The systematic uncertainty band for
the Allpix2+TCAD shown in [21] is not overlaid here because it was not
possible to decipher it from the original plots. However, as one can see in
[21], the uncertainty on the Allpix2+TCAD is as high as 5-15% for the mean
cluster size plot and 1-2% for the efficiency plot. The Allpix2+EFF
uncertainties should be identical to those (giving an overall uncertainty
higher by a factor of $\sqrt{2}$).

\begin{overpic}[width=212.47617pt]{fig/ResidualXRatio.pdf}\end{overpic}\begin{overpic}[width=212.47617pt]{fig/ResidualYRatio.pdf}\end{overpic}

Figure 7:  Left: the comparison of the resolution in $x$ behavior for
different thresholds between our Allpix2+EFF simulation and the Allpix2+TCAD
simulation (and the beam-test data) from [21]. Right: the same comparison for
the resolution in $y$. The bottom panels show the ratio of the Allpix2+TCAD
(red) or data (black) results to the Allpix2+EFF results. The graphs from [21]
either do not have statistical uncertainties at the source or the errors are
too small for the digitizer software [30]. The systematic uncertainty band for
the Allpix2+TCAD shown in [21] is not overlaid here because it was not
possible to decipher it from the original plots. However, as one can see
particularly at the low end of the thresholds range, $\lesssim 150$ e, the
uncertainty on the Allpix2+TCAD is as high as 5-15% in $y$. In $x$, at the
same threshold range, the uncertainty is approximately constant at $\sim 7\%$.
The Allpix2+EFF uncertainty should be identical to that (giving an overall
uncertainty higher by a factor of $\sqrt{2}$). Hence, while the results in $y$
are fully compatible throughout the full threshold range, the results in $x$
are compatible only above 80 e.

##  6 Conclusions and outlook

We show that a carefully constructed EFF may replace a highly detailed TCAD
mesh in the simulation chain of semiconductor sensors. This is shown for the
highly non-linear and non-trivial field of the ALPIDE sensor, where the TCAD
mesh is not available publicly. Our EFF shape resembles the one from [21] and
moreover, it reproduces the same performance to a very good extent for several
quantities in a wide range of values. The procedures discussed can be regarded
as general as long as there is at least (i) a very basic knowledge on the
field’s shape and (ii) reference data for the comparison of the sensor
performance with simulation, e.g. from beam-tests, radioactive sources, cosmic
muons, etc. We point out that, while our EFF (available in [24]) can be
already used now, it may also be further optimized to reach an even better
agreement with the TCAD-based results of [21]. In that sense, the decision
where to stop the optimization process is purely driven by pragmatic
considerations, depending on the specific application.

Our work allows the many ALPIDE users worldwide, who do not have access to the
TCAD mesh, to use our EFF as an important input for their simulation. In fact,
the tracking detector simulation of the LUXE experiment [12], which is to be
built from the production version of the ALPIDE sensors, is already using our
EFF. The production version of the ALPIDE has different dimensions compared to
the “investigator” chip used in [21]. Namely the pixel size is $29.24\times
26.88~{}\mu{\rm m}^{2}$ in $x\times y$ and 50 $\mu$m deep in $z$ with a
depletion depth of 25 $\mu$m. The EFF for this version may therefore be scaled
naturally by changing the dimensions of the pixel, where effectively only the
pitch very slightly changes and where the change in thickness has no impact on
the field since it is only defined within the depletion depth. The different
thickness does impact, however, the performance that is checked with the
Allpix2 software, but this change is unrelated to the EFF. Finally, we note
that although we use the Allpix2 software extensively for the derivation and
validation of our EFF, the resulting function can be interfaced also with
other software like [31] and it is not linked specifically to Allpix2.

## Acknowledgments

We wish to thank the colleagues from the ALICE ITS project for useful
discussions and help related to the ALPIDE sensors. We would like to
especially thank Luciano Musa, Gianluca Aglieri Rinella, Antonello Di Mauro,
Magnus Mager, Corrado Gargiulo, Felix Reidt, Ivan Ravasenga, and Ruben
Shahoyan. We also wish to thank Simon Spannagel and Paul Schütze for the
detailed and dedicated help and for the useful discussions about Allpix2 and
especially the work in [21].

This work is supported by a research grant from the Estate of Dr. Moshe Glück,
the Minerva foundation with funding from the Federal German Ministry for
Education and Research, the ISRAEL SCIENCE FOUNDATION (grant No. 708/20), the
Anna and Maurice Boukstein Career Development Chair, the Yeda-Sela SABRA
(Supporting Advanced Basic Research) fund.

## References

  * [1] Synopsys, Official website, https://www.synopsys.com/company.html (2022).
  * [2] Silvaco, Official website, https://silvaco.com/tcad/ (2022).
  * [3] W. Snoeys, Monolithic pixel detectors for high energy physics, Nucl. Instrum. Meth. A 731 (2013) 125–130.  doi:10.1016/j.nima.2013.05.073.
  * [4] M. Mager, ALPIDE, the Monolithic Active Pixel Sensor for the ALICE ITS upgrade, Nucl. Instrum. Meth. A 824 (2016) 434–438.  doi:10.1016/j.nima.2015.09.057.
  * [5] ALICE Collaboration, Expression of Interest for an ALICE ITS Upgrade in LS3, ALICE-PUBLIC-2018-013 (2018).
  * [6] G. Aglieri Rinella, The ALPIDE pixel sensor chip for the upgrade of the ALICE Inner Tracking System, Nucl. Instrum. Meth. A 845 (2017) 583–587.  doi:10.1016/j.nima.2016.05.016.
  * [7] S. Kushpil, F. Krizek, A. Isakov, Recent Results From Beam Tests of the ALPIDE Pixel Chip for the Upgrade of the ALICE Inner Tracker, IEEE Trans. Nucl. Sci. 66 (11) (2019) 2319–2323.  doi:10.1109/TNS.2019.2945234.
  * [8] S. Senyukov, J. Baudot, A. Besson, G. Claus, L. Cousin, A. Dorokhov, W. Dulinski, M. Goffe, C. Hu-Guo, M. Winter, Charged particle detection performances of CMOS pixel sensors produced in a 0.18$\nu$m process with a high resistivity epitaxial layer, Nucl. Instrum. Meth. A 730 (2013) 115–118.  arXiv:1301.0515, doi:10.1016/j.nima.2013.03.017.
  * [9] ALICE Collaboration, The ALICE experiment at the CERN LHC, JINST 3 (2008) S08002.  doi:10.1088/1748-0221/3/08/S08002.
  * [10] TowerJazz, Official website, https://towersemi.com/ (2022).
  * [11] B. Abelev, et al., Technical Design Report for the Upgrade of the ALICE Inner Tracking System, J. Phys. G 41 (2014) 087002.  doi:10.1088/0954-3899/41/8/087002.
  * [12] LUXE Collaboration, Conceptual design report for the LUXE experiment, Eur. Phys. J. ST 230 (11) (2021) 2445–2560.  arXiv:2102.02032, doi:10.1140/epjs/s11734-021-00249-z.
  * [13] sPHENIX Collaboration, An Upgrade Proposal from the PHENIX Collaboration (1 2015).  arXiv:1501.06197.
  * [14] C. Dean, The sPHENIX experiment at RHIC, PoS ICHEP2020 (2021) 731.  doi:10.22323/1.390.0731.
  * [15] DESY, The ALPIDE-based new telescope at DESY’s beam-teast facility (TO APPEAR SOON) (2022).
  * [16] R. Iuppa, S. Beolé, S. Coli, G. Gebbia, E. Ricci, S. Ricciarini, E. Serra, P. Zuccon, L. de Cilladi, The innovative particle tracker for the HEPD space experiment onboard the CSES-02, PoS ICRC2021 (2021) 070.  doi:10.22323/1.395.0070.
  * [17] G. Tambave, et al., Characterization of monolithic CMOS pixel sensor chip with ion beams for application in particle computed tomography, Nucl. Instrum. Meth. A 958 (2020) 162626.  doi:10.1016/j.nima.2019.162626.
  * [18] J. Alme, et al., A High-Granularity Digital Tracking Calorimeter Optimized for Proton CT, Frontiers in Physics 8 (2020).  doi:10.3389/fphy.2020.568243.
URL https://www.frontiersin.org/articles/10.3389/fphy.2020.568243

  * [19] S. Spannagel, K. Wolters, D. Hynds, N. Alipour Tehrani, M. Benoit, D. Dannheim, N. Gauvin, A. Nürnberg, P. Schütze, M. Vicente Barreto Pinto, Allpix2: A Modular Simulation Framework for Silicon Detectors, Nucl. Instrum. Meth. A 901 (2018) 164–172.  arXiv:1806.05813, doi:10.1016/j.nima.2018.06.020.
  * [20] P. Schütze, S. Spannagel, K. Wolters, Allpix2 user manual, version v2.3.1, https://project-allpix-squared.web.cern.ch/usermanual/allpix-manual.pdf (2022).
  * [21] D. Dannheim, K. Dort, D. Hynds, M. Munker, A. Nürnberg, W. Snoeys, S. Spannagel, Combining TCAD and Monte Carlo Methods to Simulate CMOS Pixel Sensors with a Small Collection Electrode using the Allpix2 Squared Framework, Nucl. Instrum. Meth. A 964 (2020) 163784.  arXiv:2002.12602, doi:10.1016/j.nima.2020.163784.
  * [22] D. Dannheim, et al., Comparison of small collection electrode CMOS pixel sensors with partial and full lateral depletion of the high-resistivity epitaxial layer, Nucl. Instrum. Meth. A 927 (2019) 187–193.  doi:10.1016/j.nima.2019.02.049.
  * [23] R. Brun, F. Rademakers, P. Canal, A. Naumann, O. Couet, L. Moneta, V. Vassilev, S. Linev, D. Piparo, G. GANIS, B. Bellenot, E. Guiraud, G. Amadio, wverkerke, P. Mato, TimurP, M. Tadel, wlav, E. Tejedor, J. Blomer, A. Gheata, S. Hageboeck, S. Roiser, marsupial, S. Wunsch, O. Shadura, A. Bose, CristinaCristescu, X. Valls, R. Isemann, root-project/root: v6.18/02, 10.5281/zenodo.3895860 (Aug. 2019).
URL https://doi.org/10.5281/zenodo.3895860

  * [24] A. Santra, N. Tal Hod, Code for the generation of the effective field function of the alpide sensor, https://github.com/noamhod/ALPIDE.EField (2022).
  * [25] S. Agostinelli, et al., GEANT4–a simulation toolkit, Nucl. Instrum. Meth. A 506 (2003) 250–303.  doi:10.1016/S0168-9002(03)01368-8.
  * [26] J. Allison, et al., Geant4 developments and applications, IEEE Trans. Nucl. Sci. 53 (2006) 270.  doi:10.1109/TNS.2006.869826.
  * [27] J. Allison, et al., Recent developments in Geant4, Nucl. Instrum. Meth. A 835 (2016) 186–225.  doi:10.1016/j.nima.2016.06.125.
  * [28] R. M. Muenker, Test beam and simulation studies on High Resistivity CMOS pixel sensors, Ph.D. thesis, U. Bonn (main) (2018).
URL https://inspirehep.net/literature/1782549

  * [29] I. Kremastiotis, R. Ballabriga, K. Dort, N. Egidos, M. Munker, CLICTD: A monolithic HR-CMOS sensor chip for the CLIC silicon tracker, PoS TWEPP2019 (2020) 039.  arXiv:2004.02537, doi:10.22323/1.370.0039.
  * [30] A. Rohtagi, Webplotdigitizer: Version 4.5, https://automeris.io/WebPlotDigitizer (2021).
  * [31] M. Šuljić, P. Camerini, J. W. van Hoorne, Monte Carlo simulation of charge collection processes in Monolithic Active Pixel Sensors for the ALICE ITS upgrade, Nucl. Instrum. Meth. A 950 (2020) 162882.  doi:10.1016/j.nima.2019.162882.