# Grammars and reinforcement learning for molecule optimization

Egor Kraev

###### Abstract

We seek to automate the design of molecules based on specific chemical
properties. Our primary contributions are a simpler method for generating
SMILES strings guaranteed to be chemically valid, using a combination of a new
context-free grammar for SMILES and additional masking logic; and casting the
molecular property optimization as a reinforcement learning problem,
specifically best-of-batch policy gradient applied to a Transformer model
architecture.

This approach uses substantially fewer model steps per atom than earlier
approaches, thus enabling generation of larger molecules, and beats previous
state-of-the art baselines by a significant margin.

Applying reinforcement learning to a combination of a custom context-free
grammar with additional masking to enforce non-local constraints is applicable
to any optimization of a graph structure under a mixture of local and nonlocal
constraints.111The source code to produce these results can be found at
_https://github.com/ZmeiGorynych/generative_playground_

##  1 Introduction

An important challenge in drug discovery is to find molecules with desired
chemical properties. While ultimate usefulness as a drug can only be
determined in a laboratory or clinical context, that process is expensive, and
it is thus advantageous to pre-select likely candidates in software.

While deep learning has been extensively investigated for molecular graph
encoding ([Duvenaud et al., 2015], [Kearnes et al., 2016], [Gilmer et al.,
2017]), molecule generation is still subject of active research. The simplest
natural approach to candidate molecule generation is to generate some sort of
a linear representation, such as a string of characters in the SMILES format
[Weininger, 1988], using an encoder-decoder network architecture similar to
that used in machine translation, as done in [Gómez-Bombarelli et al., 2016].

This approach’s performance was comparatively poor because a molecule’s
structure is not linear, but rather a graph which typically includes cycles,
so it falls to the model to learn how to generate SMILES strings that
correspond to chemically valid molecules - a nontrivial task that leaves the
model with little spare capacity to additionally optimize a given chemical
metric of the molecules produced.

A way to partially remedy that involves generating not the actual SMILES
strings, but a sequence of production rules of a context-free grammar (CFG)
for SMILES, as done by [Kusner et al., 2017]. That guarantees that the SMILES
strings produced are grammatically valid, putting less burden on the model to
ensure validity and thereby achieving better metrics.

However, [Kusner et al., 2017] give two reasons why this is still not
guaranteed to produce _chemically_ valid molecules: firstly, a grammatically
valid SMILES string is not guaranteed to be chemically possible (because of
atom valences being wrong, for example), and secondly, because a model
typically assumes a maximum number of steps per molecule, and the rule
expansion is not guranteed to terminate by that number of steps.

Rather than generating SMILES strings, one could directly generate a molecular
graph out of a vocabulary of valid components, an approach taken by [Jin et
al., 2018]. An important insight of that approach is that one can regard a
molecule as an acyclic graph if one considers cycles to be graph nodes, along
with atoms not part of a cycle. Thus, by first constructing such a graph and
then resolving that into an actual molecule, the molecules produced by this
approach are guaranteed to be chemically valid, so model optimization can
focus on optimizing the properties of the actual molecule, resulting in state-
of-the-art scores.

Unfortunately the approach of [Jin et al., 2018] is also quite complex, using
a two-part internal latent representation, a graph message-passing network, a
tree message-passing network, and scoring for subgraphs as they are being
assembled. Information transfer between nodes was done via message passing
along the graph edges using an RNN, with the information distance between
nodes thus being dependent on the graph traversal order. Furthermore, it seems
to require every possible cycle we generate to be contained in a predefined
vocabulary, leading to potentially huge vocabularies, especially if we wish to
take into account features such as chirality and stereoisomers.

Our major contribution is an approach that also guarantees the chemical
validity of any molecule produced, but is simpler, more expressive, and more
memory-efficient, and results in scores substantially exceeding those of [Jin
et al., 2018].

We consider that a context free grammar (CFG) is a very natural way of
representing an acyclic graph. We formulate a new grammar for SMILES strings
that guarantees correct atom valences by construction, and allows for
arbitrarily long, branching chains of aliphatic cycles, and the more common
kinds of aromatic cycle structures. Due to the recursive nature of a CFG, we
can represent a combinatorially large number of possible cycle structures with
a small number of expansion rules (for example, an aliphatic 5-cycle can
consist of any 5 sub-structures with 2 free valences each, some of which can
also form part of other cycles).

To address the issues identified by [Kusner et al., 2017] we introduce
additional masking into the process of selecting the next production rule. One
kind of masking tracks the ‘terminal distance’ of the current sequence, making
sure the rule expansion terminates before the maximum number of steps allowed
by the model.

Another kind of masking makes sure all cycles are closed, and enforces minimum
and (optionally) maximum cycle length (cycles of length 2 are illegal in the
SMILES format). We achieve that by attaching additional information, such as a
unique cycle ID, to any grammar token that indicates creation of a new cycle,
and propagating it when expanding that token.

All the molecules produced using that additional masking are still valid
according to our CFG, so we can use standard CFG parsers (for example, `nltk`
[Bird et al., 2009]) to decompose the molecules from the ZINC database
[Weininger, 1988] into the production rules of our grammar.

As discussed in Section 2.7, our grammar is only expressive enough to
represent a little over a third of the molecules in the ZINC database, not due
to any inherent limitation of the approach, but merely because that is
sufficient to represent the optimized molecules from the literature, eg those
in [Kusner et al., 2017] and [Jin et al., 2018], and to generate new molecules
with scores exceeding state of the art.

The second innovation we make is that instead of using an RNN, be it over the
sequence of rules, or via graph message-passing, we base our neural model on
the Transformer architecture [Vaswani et al., 2017]. The Transformer’s
information distance between any two inputs is always one, so we have the full
information about the sequence so far at our immediate disposal when expanding
the next rule. As the memory requirements of the Transformer grow
approximately as a square of the number of generation steps, it helps that the
new grammar gives a concise representation of a molecule (with approximately
half as many rules needed per atom as in the grammar used in [Kusner et al.,
2017], see Section 2.7)

We could have used this architecture in an autoencoder, followed by Bayesian
optimization, as has been done in the literature so far. However, as the
properties of that approach are well known by now, we decide to try a
different tack, namely consider the molecule generation problem as a
reinforcement learning problem, and optimize using best-of-batch policy
gradient.

We first train a baseline model off-policy to teach it the distribution of
molecules in the ZINC database, and then conduct a series of optimization
experiments. Using best-of-batch policy gradient, with varying degrees of
anchoring to the baseline model, demonstrates a smooth tradeoff between
objective optimization and similarity to the molecules in the ZINC dataset,
and yields scores exceeding state of the art even with the more conservative
parameter settings.

##  2 Generating guaranteed valid SMILES strings

###  2.1 SMILES

SMILES [Weininger, 1988] is a common format for describing molecules with a
character string.

Briefly, atoms following each other in the string are assumed to be connected
with a single bond, double bonds are denoted with `=`, and triple bonds with
`#`, thus for example $CO_{2}$ is represented as `O=C=O` and $N_{2}$ as `N#N`.
Side branches are represented by brackets, thus $SO_{3}$ is `O=S(=O)=O`,
aromatic atoms (those contributing one valence to an aromatic ring) are
represented as lowercase letters, and cycles are represented by numerals
following an atom, with atoms being followed by the same numeral understood to
be connected; thus for example benzene is `c1ccccc1`. Valences not used in the
explicitly specified bonds are assumed to have hydrogens attached, for example
$H_{2}O$ is `O`.

Existing molecule databases, such as the ZINC database [Sterling and Irwin,
2015] we use, are in SMILES format, and can be processed by software such as
`http://www.rdkit.org` [Landrum,].

###  2.2 Context-free grammar

Non-binary trees are a natural description of molecular structure, be it using
the SMILES approach (cutting the cycles), or the approach of [Jin et al.,
2018] (treating cycles as graph nodes to be resolved into individual atoms at
a later stage).

A context-free grammar (CFG) is, in turn, a natural way of describing non-
binary trees. A CFG is a set of tokens, and a set of production rules, each
rule mapping a token (the rule’s ‘left hand side’) onto a sequence of tokens
(the rule’s ‘right hand side’). Tokens that occur on the left hand side of at
least one production rule are called nonterminals, the other tokens terminals.

A sequence that’s valid according to a CFG is generated by starting with the
sequence consisting of the root token, and then repeatedly taking a (typically
the leftmost) nonterminal token in the sequence, selecting a rule whose left
hand side equals that token, and replacing the token with the right hand side
of that rule. That process terminates once there are no nonterminals left in
the sequence. If we consider that process, not as replacing the token in
question, but as inserting the new tokens as its children, the process of
generating the sequence is equivalent to generating a non-binary tree whose
nodes are nonterminals and whose leaves are terminals. That tree can be
naturally mapped to the tree describing molecular structure, as we demonstrate
by providing a CFG that achieves that.

A context-free grammar thus allows us to efficiently represent molecules as
non-binary trees with some constraints on tree structure. That efficiency
matters for model training: RNN-style models have problems passing information
over very long distances. The Transformer architecture we use doesn’t have
that problem, but its RAM requirements grow approximately as a square of the
number of generation steps (in our case, rules needed to generate a given
molecule). Finally, for any architecture, length of the sequence used to
represent a molecule directly affects training time.

Thus the average number of tokens needed to represent a molecule constrains
the maximum size of molecules it’s practical to produce.

Using a CFG allows us to combine a granular representation (yielding a
combinatorial variety of structures), with explicit representation of
particularly frequent or complex combinations (eg double aromatic cycles).

To generate guaranteed valid molecules, we introduce a new grammar (full
listing in Appendix A) that respects atomic valences and explicitly introduces
cycles.

###  2.3 Respecting valences

The initial rules are

    smiles -> nonH_bond
    smiles -> initial_valence_1 bond
    smiles -> initial_valence_2 double_bond
    smiles -> initial_valence_3 triple_bond

Any token ending in `bond` represents an open bond to which some atom must be
attached; `nonH_bond` requires that that atom not be hydrogen (see Section
2.3.1 on treatment of hydrogen). The molecule can then be grown via rules such
as

    nonH_bond -> valence_2 bond
    double_bond -> ’=’ valence_2
    double_bond -> ’=’ valence_3 bond
    double_bond -> ’=’ valence_4 double_bond
    triple_bond -> ’#’ valence_3
    triple_bond -> ’#’ valence_4 bond

So far it seems as if we could only produce linear molecules. However, the
`valence` tokens can expand to individual atoms as well as to structures that
contain branches, for example

    valence_4 -> ’C’
    valence_3 -> ’N’
    valence_3 -> valence_4 branch
    valence_2 -> ’O’
    valence_2 -> valence_3 branch
    valence_2 -> valence_4 ’(’ double_bond ’)’

This allows us, so far, to generate arbitrary tree-shaped molecules, ie those
without cycles.

####  2.3.1 Implicit hydrogens

A complication in the above approach is that the SMILES format allows for
implicit hydrogens, ie any valence of any atom that’s not specified is assumed
to be occupied by hydrogen. Thus in effect, these hydrogen atoms are
represented by empty strings, which is a problem for us as a CFG does not
allow to expand a token into nothing (otherwise parsing would be much harder).

If the implicit hydrogen atoms were represented by for example lowercase `h`,
we could have grammar rules such as

    bond -> ’h’
    bond -> nonH_bond
    branch -> ’h’
    branch -> ’(’ nonH_bond ’)’

If all we cared about was generating new SMILES strings, we could have used
the above 4 rules and then inserted a post-processing step that replaces each
`h` with an empty string. However, we also want to be able to parse existing
SMILES strings into our grammar, so we go with a two-step approach: we start
with a grammar as described above, and then eliminate all occurrences of
`bond` by replacing them with either nothing or `nonH_bond`, and similarly for
`branch`.

Thus for example, the rule

    valence_3 -> valence_4 branch

is replaced by the two rules

    valence_3 -> valence_4
    valence_3 -> valence_4 ’(’ nonH_bond ’)’

###  2.4 Cycles

The different kinds of cycles are introduced by rules such as

    nonH_bond -> aliphatic_ring
    aliphatic_ring -> valence_3_num cycle_bond
    aliphatic_ring -> valence_4_num cycle_double_bond
    valence_2 -> vertex_attached_ring
    vertex_attached_ring -> valence_4_num ’(’ cycle_bond ’)’

The challenging bit here is the numbering - the two ends of a cycle must be
marked with the same numeral, to know we should connect them; furthermore,
SMILES reuses these numerals, so if the numeral `1` was used in a cycle that’s
already closed, we can use it again for the next cycle we start. Thus cycle
numeral management is a task that’s quite hard to solve using just a context-
free grammar.

We solve this by attaching additional properties to some nonterminal tokens.
When we expand a nonterminal which has `ring` as a substring, indicating a new
cycle, we give it a unique cycle identifier, which is attached as a property
to any tokens resulting in that expansion, which have substrings `num` or
`cycle`. When these are in turn expanded, the cycle id is again propagated.
The grammar rules for cycle propagation are formulated in a way that
guarantees there are ever at most two tokens tagged by any one cycle ID, in
fact exactly two until we expand them into numerals.

####  2.4.1 Numeral assignment

When a nonterminal is about to be expanded into a numeral, we know by
construction that it has a cycle id. We have two cases to consider: firstly,
there exists another token with this cycle id. In that case, we scan the token
sequence so far to determine the lowest available numeral, and generate a mask
that makes sure that specific numeral will be chosen; and store the pair
(cycle id, numeral) in a cache. In the second case, when there are no tokens
with this cycle id, we look up the numeral to use in the cache, by cycle id.

####  2.4.2 Cycle size

A cycle with only two atoms (same numeral attached to atoms that are already
connected by a regular bond) is considered illegal in SMILES. To prevent
these, we attach and propagate another property to nonterminals that are part
of a cycle, namely ring size, incrementing it with each rule expansion. Cycle
size masking then forbids the rules that would have led to premature cycle
closure. As cycles of length more than 8 don’t occur in the database (and our
scoring function penalizes cycles with length more than 6), we also forbid
cycle expansion beyond length 8.

####  2.4.3 Cycle chaining

Aliphatic cycles can form arbitrarily long, potentially branching chains, due
to a pair of rules that allow starting a new cycle while constructing another:

    cycle_bond -> aliphatic_ring_segment cycle_bond
    aliphatic_ring_segment -> valence_3 ’(’ cycle_bond ’)’ valence_3_num

In the first of these rules, the `cycle_bond` token on the right hand side
will get the same cycle id as that on the left hand side; while the
`aliphatic_ring_segment` token will trigger generation of a new cycle id, that
will be propagated to the appropriate nonterminals on the right hand side of
the second rule, in this case `cycle_bond` and `valence_3_num`.

This is just an example, a similar pattern is used to attach an aliphatic
cycle to aromatic cycles described below.

####  2.4.4 Aromatic cycles

Aromatic cycles and chains thereof must fulfil additional relationships
between the number of atoms and the number of double bonds in each cycle.
Solving this in a general way is beyond the scope of current work, we limit
ourselves to rules that allow us to generate any single aromatic cycle of
length 5 or 6, and common patterns of linked pairs of aromatic cycles,
possibly with chains of aliphatic cycles attached. For example, here are some
of the rules used in generating a coupled pair of aromatic rings of size 6 and
5:

    double_aromatic_ring -> ’c’ num1 aromatic_atom aromatic_atom
                            aromatic_atom ’c’ num ’n’ num1 aromatic_atom aromatic_atom
                            aromatic_atom_num
    aromatic_atom -> ’n’
    aromatic_atom -> ’c’
    aromatic_atom -> ’c’ ’(’ nonH_bond ’)’
    aromatic_os -> ’o’
    aromatic_os -> ’s’
    aromatic_atom_num -> ’c’ num

###  2.5 Making sure the expansions terminate

The final challenge we address is making sure the rule expansion terminates
before the maximum allowed number of expansion rules. To do this, we define
the concept of _terminal distance_ of a token, defined as the length of the
shortest sequence of rules needed to transform that token into a sequence
consisting only of terminals.

We calculate that distance for all tokens reachable from the root token by
means of the following algorithm:

  1. 1.

Define a set $T$ of all tokens observed so far, and seed it with the root
token `smiles`.

  2. 2.

Iterate over the elements $t$ of $T$.

    1. (a)

For each $t$ we firstly apply all applicable rules and add any new tokens
generated thereby to $T$, seeding their terminal distance with 0 for terminals
and $\infty$ for nonterminals.

    2. (b)

We then calculate the terminal distance of a token $t$ as one plus the minimum
over all applicable rules of the sum of terminal distances of the tokens on
the right hand side of the respective rule:

$TD(t)=1+\min_{r\in G:r.lhs=t}\left(\sum_{t^{\prime}\in
r.rhs}TD(t^{\prime})\right)$ (1)

We repeat step 2. this until convergence, that is, until no new tokens are
observed and the terminal distance for known tokens no longer changes after an
iteration.

We then define the change in terminal distance made by a production rule as

$\Delta TD(r)=\sum_{t^{\prime}\in r.rhs}TD(t^{\prime})-TD(r.lhs)$ (2)

By definition of the terminal distance, for every nonterminal $t$ there exists
a production rule $r$ such that $r.lhs=t$ and $\Delta TD(r)=-1$.

Finally, we define the terminal distance of a token sequence as the sum of
terminal distances of the individual tokens.

We use the terminal distance concept to make sure the rule expansion
terminates before the maximum number of steps, in the following manner: at
each step, we consider the number $s$ of steps left and the terminal distance
$td$ of the sequence generated so far.

At each step, we make sure $td\leq s$, using induction. First, we choose the
maximum rule sequence length to be larger than the terminal distance of the
root symbol. Second, at each rule selection step we consider all rules $r$
whose left hand side is the next nonterminal to expand, and only allow those
where $\Delta TD(r)+td\leq s-1$. That is a nonempty set because $\Delta
TD(r^{\prime})=-1$ for at least one applicable rule $r^{\prime}$. Thus, by the
time we run out of steps, that is, $s=0$, we know $td=0$, that is our token
sequence consists only of nonterminals.

Note that in our case, we have to apply this algorithm not to the tokens of
our original grammar, but to the extended tokens (including cycle size
information); and likewise use the production rules consisting of the original
rules plus the propagation of cycle size information222This is an additional
limit maximal cycle size, namely to make sure the number of extended tokens is
finite. This is because cycle masking forbids some rules (eg those that would
create a cycle of two atoms) that would be allowed by our original grammar.
Effectively, the rules for propagating additional information, along with the
original CFG rules, induce a CFG on the space of extended tokens, and it is
the terminal distance within that CFG that we must consider.

Also for computational efficiency reasons (to not multiply the number of
tokens unnecessarily) we disregard the cycle id propety of the tokens when
calculating terminal distance, as it only affects choice of numeral but not
terminal distance.

###  2.6 Limits of context-free grammars

Of the above modifications to a pure context-free grammar, cycle size masking
could have been implemented by including the size value into the token string,
and making a copy of each rule whose left hand side is that token, for every
size value that occurs. This would increase the number of rules by a couple
dozen, but stay within the limits of a CFG.

On the other hand, numeral choice and terminal masking do not appear possible
to achieve using a CFG alone. However, these merely restrict the list of
possible production rules at any step - but the resulting SMILES string is
still valid according to the original CFG without these extra attributes. This
is imporant because it allows us to use our original CFG to parse known
molecules, eg from the ZINC database, to train our model on.

###  2.7 Grammar conciseness and expressiveness

Our approach leads to a more concise way of representing molecules using
production rules, with an average 2.85 production rules per atom and 62.8
rules per molecule in the ZINC dataset, compared to average 5.46 rules per
atom and 120.8 rules per molecule in the grammar used by [Kusner et al.,
2017]. This in turn allowed us to generate larger molecules (about 100 atoms
for a typical molecule that uses the maximum number of steps, and over 400
atoms in special cases) on commodity hardware.

Our grammar can represent arbitrarily long, branching chains of aliphatic
cycles, as well as single aromatic cycles with five or six atoms, pairs of
aromatic cycles, and arbitrarily long branching chains of aliphatic cycles
attached to aromatic ones.

Because our grammar is more restrictive than a fully generic SMILES grammar
(to make it easier for us to generate guaranteed valid molecules), it can’t
represent all the molecules in the ZINC database, but rather a little over a
third of the molecules, 92K out of 250K. Our grammar is also expressive enough
to represent the recent state-of-the art constructed molecules, such as those
in [Kusner et al., 2017] and [Jin et al., 2018].

This approach can be extended to cover the whole ZINC database by introducing
additional production rules - certainly by brute force, by adding to the
grammar every unknown pattern that the parser encounters; and most likely also
in a more elegant way. Extending our grammar to represent most or all
molecules from that dataset is subject of ongoing work.

There is something of a tradeoff between expressiveness and guaranteed
correctness - as one adds more features to the grammar, it becomes harder to
maintain the guaranteed-correctness property. Fortunately, that property is
not actually necessary for our approach (best-of-batch policy gradient) to
work - it’s sufficient for a large fraction, say 75%, of each batch, to be
valid for the optimization to succeed.

Future attempts to expand the grammar’s expressiveness might make use of that
by relaxing the guaranteed-correctness property if needed.

##  3 Model choice

###  3.1 Reinforcement learning

Most of the literature so far ([Gómez-Bombarelli et al., 2016],[Kusner et al.,
2017],[Jin et al., 2018]) on molecule generation has focussed on training an
autoencoder for the chosen molecule representation, using a database of known
molecules such as ZINC, then doing Bayesian optimization over the latent
space.

We choose a different approach, treating molecule creation as a reinforcement
learning problem. In our case, at each step the action taken by the network is
choosing the next production rule; the state is the full sequence of the
production rules so far; and the reward is nonzero for the step following
molecule completion (ie when the sequence defined by applying the chosen
production rules contains no more nonterminals), and zero for all other steps.
Details of the reward specification are in Sections 3.3 and 3.4.

The first step of the autoencoder/Bayesian optimization approach, training the
autoencoder on a database of known molecules, typically served two objectives:
firstly, ensuring that the model learns to produce chemically valid molecules
(for the majority of models where that is not guaranteed by construction);
secondly, ensuring that the molecules the model produced are ‘similar’ to the
molecules the model has been trained on.

While this approach has many desirable traits, it also has disadvantages:
firstly, once the autoencoder is trained, its notion of ‘similarity to known
good molecules’ is cast in stone, and no subsequent optimization over the
latent space would allow it to produce, for example, models with similar
structure but substantially larger than the ones it was trained on. In
contrast, casting molecule generation as a reinforcement learning problem
allows us to treat ‘similarity to known good molecules’ as the gradual concept
that it is, for example (as we do) by including a variably-weighted loss term
penalizing the coefficient distance from a version of the model trained off-
policy on known-good molecules. This will allow us to search over a greater
space of possible molecules, while maintaining control over the degree of
similarity to the known-good ones.

Secondly, although the molecules in the latent space neighborhood of a certain
molecule are guaranteed to be ‘similar’ to it in some sense, that neighborhood
is not guaranteed to contain _all_ the ‘similar’ molecules we might be
interested in. In comparison, in a reinforcement learning problem we can
explicitly specify any proximity metric we desire, and simply include it in
the reward.

Thirdly, learning a latent space representation of the space of all the
molecules in the database is a reasonably hard task, which is not necessary if
all we want is to generate molecules that maximize some metric and are similar
to existing ones. We conjecture that simply teaching a model to have a high
probability of producing those molecules is a computationally cheaper task,
which can be used to achieve the same aim via anchoring (see below).

Finally, given the extremely non-Euclidean structure of the space of all valid
molecules, it seems promising to conduct the search for optimal molecules
directly on a representation that seems close to it, in our case the sequence
of grammar production rules tailored to molecular structure, rather than on
the Euclidean latent space.

###  3.2 Architecture

As the molecular properties we seek to optimize are essentially nonlocal, and
because a grammar-rule representation of a branching tree means that even
atoms that are nearby on the tree can end up generated by rules that are far
apart in the rules sequence, we use the Transformer architecture [Vaswani et
al., 2017], chosen because any two items in the sequence are directly linked
by its attention mechanism. Our implementation of the Transformer decoder
outputs one vector of logits per call (for each molecule in the batch), to be
used in choosing the next production rule for each sequence in the batch; and
omits the calculation steps that use the sequence generated by the encoder as
it doesn’t exist in our caase.

We use 6 layers, 6 heads, $d_{k}=d_{v}=16$, $d_{model}=128$,
$d_{inner\\_hid}=256$.

Prior to calculating the log-softmax to be used for sampling the next grammar
rule, we calculate a mask as described in Section 2 and subtract 1e6 from the
logits at the indices forbidden by the mask.

###  3.3 Training

We use same maximal sequence length as [Kusner et al., 2017], 277. That allows
us a maximum batch size of 40 on an AWS p2.xlarge instance.

In the first stage, we train the model off-policy on 92K molecules from the
ZINC database [Weininger, 1988] that can be represented by our grammar. In
order to reward the model for producing molecules similar to the ones it
observes, we use simple policy gradient loss (3). Here $s$ goes over the model
steps, and $\pi_{s}\left(r_{s}\right)$ is the model-produced probability
(after applying all masking) of choosing rule number $r_{s}$ at step $s$.

$loss=-\sum_{s=0}^{S}\log\left(\pi_{s}\left(r_{s}\right)\right)$ (3)

We optimize in batches of 40 molecules, using the average loss for the batch,
and Adam optimizer with a learning rate of 1e-4. After each batch, we do an
on-policy simulation of a batch of 40 molecules with the same loss function,
to judge convergence - but without updating the model coefficients.

To judge convergence, we compare the distribution of log probability, logP, SA
score, and the number of aromatic cycles in the molecules produced by the
model to those in the database, and find that they have converged after 15
epochs of training (30K batches, about one week hours on a p2.xlarge). We take
the model thus trained as our base model.

We then proceed to on-policy training using best-of-batch policy gradient with
base model anchoring. That is, our loss is the log-likelihood _for the best-
scoring molecule in the batch_ plus the L2-distance between the coefficients
of the current model and the base model, multiplied by a weight $w_{a}$(4).
The purpose of the second term (‘anchoring’) is to prevent catastrophic
forgetting, and to control the degree of similarity of optimized model to the
base model (and thus to the molecules in the database).

$loss=-\sum_{s=0}^{S}\log\left(\pi_{s}\left(r_{s}\right)\right)+w_{a}\left|\vec{p}-\vec{p}_{base}\right|_{2}^{2}$
(4)

The reward function is used to choose the best molecule of the batch, but its
value for the best molecule does not affect the loss function, which remains a
simple log likelihood for the best molecule.

We found that using the best molecule in the batch, rather than the batch
average, during the second phase was crucial for successful optimization -
without it, as soon as the search found a somewhat successful molecule, the
optimization tended to converge to always producing that molecule; while
taking best-of-batch enabled continued exploration.

###  3.4 Optimization and the reward function

We optimize the metric used in prior literature ([Gómez-Bombarelli et al.,
2016], [Kusner et al., 2017], [Jin et al., 2018]), namely logP minus the
synthethic accessibility score, penalized for cycles of more than 6 atoms.
Each of those three terms is evaluated for each of the molecules in the ZINC
dataset using code shared by [Kusner et al., 2017], and their mean and
standard deviation are calculated. These are used to normalize each of the 3
score constituents, before adding them up to calculate a molecule’s score.

The overall reward function is thus represented by (5).

$R=logP_{norm}+SA_{norm}+C_{norm}$ (5)

For reasons explained in Section 4, we use an extended version of (5), namely
(6).

$R=logP_{norm}+SA_{norm}+C_{norm}+w_{SA}\min\left(SA_{norm},0\right)-w_{ac}\min\left(num\\_aromatic\\_cycles-5,0\right)$
(6)

Here the first three terms on the right are the standard scoring function used
in the literature, and the final two terms allow us to penalize low SA scores
and too large a number of aromatic cycles in a molecule.

##  4 Results

Depending on how harshly we penalize the distance to the base model (choice of
$w_{a}$), our model produces a range of molecules with scores substantially
exceeding state-of-the-art, with varying sizes, as shown in Table 1 and
discussed in detail below.

Source | Score | SA | Norm.SA | Aromatic rings
---|---|---|---|---
Average over ZINC database | 0.0 | -2.82 | 0.0 | 1.85
Standard deviation over ZINC database | 2.07 | 0.77 | 1.0 | 0.97
Kusner et al. | 2.93 | -1.93 | 1.35 | 1
Jin et al. | 5.30 | -2.33 | 0.87 | 5
Strong Anchor, SA penalty 1st | 5.68 | -1.96 | 1.31 | 5
Strong Anchor, SA penalty 2nd | 5.60 | -2.62 | 0.52 | 4
Strong Anchor, SA penalty 3rd | 5.33 | -2.87 | 0.21 | 5
Weak Anchor, SA and aromatic cycle penalty 1st | 9.48 | -2.83 | 0.27 | 7
Weak Anchor, SA and aromatic cycle penalty 2nd | 8.73 | -2.98 | 0.08 | 5
Weak Anchor, SA and aromatic cycle penalty 3rd | 8.44 | -2.67 | 0.46 | 6
Weak Anchor, SA penalty | 12.25 | -2.98 | 0.08 | 12
Unconstrained | 46.45 | -1.00 | 2.46 | 13
Table 1: Statistics of optimized molecules (higher scores are better)

When trying to optimize (5) without constraning similarity to known molecules
($w_{a}=0$) and using the original scoring function from the literature (that
is, $w_{SA}=0$ and $w_{ac}=0$), we’ve found that the optimization tried to get
as many aromatic cycles as possible, leading to a large logP but very low SA
score, and then attach a large number of halogen atoms to these to bring the
SA score back up (see Figure 1, right). While this allowed for an order of
magnitude increase in the molecule score (‘Unconstrained’ in Table 1), such
molecules don’t appear to be realistic candidates for having useful
properties.

We next added weak anchoring ($w_{a}=1e8$) and a penalty on negative
normalized SA score ($w_{SA}=20$). The result is shown in Table 1 as ‘Weak
Anchor, SA penalty’ and in Fig 1, left. We see that the introduced constraints
temper but don’t entirely remove the tendency of the optimization to blindly
maximize the score by adding a lot of aromatic cycles.

To counteract that, we add a penalty on the number of aromatic cycles
exceeding 5 ($w_{ac}=5$). The three top-scoring molecules are shown in Table 1
as ‘Weak Anchor, SA penalty and aromatic cycle penalty’, and in Figure 2. We
see that even after constraining the number of aromatic cycles, we can achieve
scores well in excess of 8, and a larger variety of atoms and structures than
in the previous cases. However, we still see many more halogens than in a
typical molecule from the ZINC database.

For a final set of simulations, we increase the anchoring weight to
$w_{a}=2e8$, keep the SA penalty as before, and remove the aromatic cycle
penalty. The three top-scoring molecules are shown in Table 1 as ‘Strong
Anchor, SA penalty’, and in Figure 3. We see that these molecules appear much
more similar to the ZINC database - in addition to the aromatic cycles of 6
carbon atoms we see an aliphatic cycle, aromatic 5-cycles containing sulphur
and nitrogen, and just one chlorine atom. All of these molecules’ scores
exceed those of the top molecule of [Jin et al., 2018].

Figure 1: Top molecule with weak anchoring and SA penalty (left) and no
anchoring (right) Figure 2: Top 3 molecules with weak anchoring, and SA and
aromatic cycle penalty Figure 3: Top 3 molecules with strong anchoring and SA
penalty

##  5 Conclusion

A major contribution of [Jin et al., 2018] was regarding the molecule as a
junction tree, with any cycles represented as nodes on par with non-cycle
atoms. This made it possible to create a model guaranteed to generate valid
molecules. However, their implementation was comparatively complex and also
required all possible cycles to be enumerated in a dictionary upfront.

Here we have shown that the same idea can be implemented in a simpler fashion,
using a custom context-free grammar, combined with some additional masking for
choosing the next production rule.

A properly constructed context-free grammar gives us the best of both worlds:
the ability, with a small number of production rules, to represent all
possible valid combinations of atoms to build a given structure, as in [Kusner
et al., 2017]; and at the same time, the ability to explicitly specify complex
structures (eg linked aromatic rings) upfront using a single rule per
structure, similarly to [Jin et al., 2018], while letting further rule
expansions supply the detail. All that is done as part of a single grammar,
with the decomposition of a molecule into a production rule sequence being
done using a standard CFG parser, and the generation of new molecules being
done by any model able to recursively produce a sequence of tokens (we chose
the Transformer because of its nonlocal information propagation, but an RNN
stack could have been used as well).

To solve some of the shortcomings of a purely CFG-based approach, we
propagated additional information with the tokens during rule expansion, and
used it to limit the set of allowed production rules at each step.

This approach, combining a custom CFG with additional rule masking that uses
non-local information, is applicable beyond molecules, to any domain where we
need to optimize complex graph structures under a mixture of local and
nonlocal constraints: the local constraints can be taken care of in the CFG
design, and the nonlocal ones enforced via additional masking.

Finally, we show that molecule optimization can be naturally cast as a
reinforcement learning problem, and the state-of-the-art results we get even
with the very basic policy gradient method suggest further scope for
improvement using more advanced approaches, for example methods based on tree
search.

## References

  * [Bird et al., 2009] Bird, S., Klein, E., and Loper, E. (2009).  Natural Language Processing with Python.  O’Reilly Media.
  * [Duvenaud et al., 2015] Duvenaud, D. K., Maclaurin, D., Iparraguirre, J., Bombarell, R., Hirzel, T., Aspuru-Guzik, A., and Adams, R. P. (2015).  Convolutional networks on graphs for learning molecular fingerprints.  In Cortes, C., Lawrence, N. D., Lee, D. D., Sugiyama, M., and Garnett, R., editors, Advances in Neural Information Processing Systems 28, pages 2224–2232. Curran Associates, Inc.
  * [Gilmer et al., 2017] Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E. (2017).  Neural message passing for quantum chemistry.  CoRR, abs/1704.01212.
  * [Gómez-Bombarelli et al., 2016] Gómez-Bombarelli, R., Duvenaud, D. K., Hernández-Lobato, J. M., Aguilera-Iparraguirre, J., Hirzel, T. D., Adams, R. P., and Aspuru-Guzik, A. (2016).  Automatic chemical design using a data-driven continuous representation of molecules.  CoRR, abs/1610.02415.
  * [Jin et al., 2018] Jin, W., Barzilay, R., and Jaakkola, T. S. (2018).  Junction tree variational autoencoder for molecular graph generation.  CoRR, abs/1802.04364.
  * [Kearnes et al., 2016] Kearnes, S., McCloskey, K., Berndl, M., Pande, V., and Riley, P. (2016).  Molecular graph convolutions: moving beyond fingerprints.  Journal of Computer-Aided Molecular Design, 30:595–608.
  * [Kusner et al., 2017] Kusner, M. J., Paige, B., and Hernández-Lobato, J. M. (2017).  Grammar Variational Autoencoder.  ArXiv e-prints.
  * [Landrum, ] Landrum, G.  Rdkit: Open-source cheminformatics.
  * [Sterling and Irwin, 2015] Sterling, T. and Irwin, J. J. (2015).  Zinc 15 – ligand discovery for everyone.  Journal of Chemical Information and Modeling, 55(11):2324–2337.  PMID: 26479676.
  * [Vaswani et al., 2017] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, L., and Polosukhin, I. (2017).  Attention is all you need.  CoRR, abs/1706.03762.
  * [Weininger, 1988] Weininger, D. (1988).  Smiles, a chemical language and information system. 1. introduction to methodology and encoding rules.  J. Chem. Inf. Comput. Sci., 28(1):31–36.

##  Appendix A Full Listing of the Grammar

    smiles -> nonH_bond
    smiles -> initial_valence_1
    smiles -> initial_valence_1 nonH_bond
    smiles -> initial_valence_2 double_bond
    smiles -> initial_valence_3 triple_bond
    initial_valence_1 -> ’F’
    initial_valence_1 -> ’Cl’
    initial_valence_1 -> ’Br’
    initial_valence_1 -> ’I’
    initial_valence_1 -> ’[’ ’O’ ’-’ ’]’
    initial_valence_1 -> ’[’ ’N’ ’H’ ’3’ ’+’ ’]’
    initial_valence_2 -> ’O’
    initial_valence_2 -> ’S’
    initial_valence_3 -> ’[’ ’C’ ’@’ ’H’ ’]’
    initial_valence_3 -> ’[’ ’C’ ’@’ ’@’ ’H’ ’]’
    initial_valence_3 -> ’N’
    initial_valence_3 -> ’[’ ’N’ ’H’ ’+’ ’]’
    nonH_bond -> valence_1
    nonH_bond -> valence_2
    nonH_bond -> valence_2 nonH_bond
    nonH_bond -> valence_3 double_bond
    nonH_bond -> valence_4 triple_bond
    double_bond -> ’=’ valence_2
    double_bond -> ’=’ valence_3
    double_bond -> ’=’ valence_3 nonH_bond
    double_bond -> ’=’ valence_4 double_bond
    triple_bond -> ’#’ valence_3
    triple_bond -> ’#’ valence_4
    triple_bond -> ’#’ valence_4 nonH_bond
    valence_4 -> ’C’
    valence_4 -> ’[’ ’C’ ’@’ ’]’
    valence_4 -> ’[’ ’C’ ’@’ ’@’ ’]’
    valence_4 -> ’[’ ’N’ ’+’ ’]’
    valence_3 -> ’[’ ’C’ ’@’ ’H’ ’]’
    valence_3 -> ’[’ ’C’ ’@’ ’@’ ’H’ ’]’
    valence_3 -> ’N’
    valence_3 -> ’[’ ’N’ ’H’ ’+’ ’]’
    valence_3 -> valence_4
    valence_3 -> valence_4 ’(’ nonH_bond ’)’
    valence_2 -> ’O’
    valence_2 -> ’S’
    valence_2 -> ’S’ ’(’ ’=’ ’O’ ’)’ ’(’ ’=’ ’O’ ’)’
    valence_2 -> valence_3
    valence_2 -> valence_3 ’(’ nonH_bond ’)’
    valence_2 -> valence_4 ’(’ double_bond ’)’
    valence_1 -> ’F’
    valence_1 -> ’Cl’
    valence_1 -> ’Br’
    valence_1 -> ’I’
    valence_1 -> ’[’ ’O’ ’-’ ’]’
    valence_1 -> ’[’ ’N’ ’H’ ’3’ ’+’ ’]’
    valence_1 -> valence_2
    valence_1 -> valence_2 ’(’ nonH_bond ’)’
    valence_1 -> valence_3 ’(’ double_bond ’)’
    valence_1 -> valence_4 ’(’ triple_bond ’)’
    nonH_bond -> valence_2 slash valence_3 ’=’ valence_3 slash valence_2
    slash -> ’/’
    slash -> ’\’
    nonH_bond -> aliphatic_ring
    aliphatic_ring -> valence_3_num cycle_bond
    aliphatic_ring -> valence_4_num cycle_double_bond
    valence_2 -> vertex_attached_ring
    vertex_attached_ring -> valence_4_num ’(’ cycle_bond ’)’
    cycle_bond -> valence_2 cycle_bond
    cycle_bond -> valence_3 cycle_double_bond
    cycle_double_bond -> ’=’ valence_3 cycle_bond
    cycle_bond -> valence_2_num
    cycle_double_bond -> ’=’ valence_3_num
    nonH_bond -> aliphatic_ring_segment
    nonH_bond -> aliphatic_ring_segment nonH_bond
    cycle_bond -> aliphatic_ring_segment cycle_bond
    aliphatic_ring_segment -> valence_3 ’(’ cycle_bond ’)’ valence_3_num
    aliphatic_ring_segment -> valence_4 ’(’ cycle_bond ’)’ ’=’ valence_4_num
    aliphatic_ring_segment -> valence_4 ’(’ cycle_double_bond ’)’ valence_3_num
    nonH_bond -> aromatic_ring_5
    nonH_bond -> aromatic_ring_6
    starting_aromatic_c_num -> ’c’ num
    aromatic_atom -> ’n’
    aromatic_atom -> ’c’
    aromatic_atom -> ’c’ ’(’ nonH_bond ’)’
    aromatic_os -> ’o’
    aromatic_os -> ’s’
    aromatic_os -> ’n’ ’(’ nonH_bond ’)’
    aromatic_os -> ’[’ ’n’ ’H’ ’]’
    aromatic_atom_num -> ’n’ num
    aromatic_atom_num -> ’c’ num
    aromatic_atom_num -> ’c’ num nonH_bond
    aromatic_os_num -> ’o’ num
    aromatic_os_num -> ’s’ num
    aromatic_os_num -> ’n’ num nonH_bond
    nonH_bond -> double_aromatic_ring
    double_aromatic_ring -> ’c’ num1 aromatic_atom aromatic_atom aromatic_atom ’c’ num ’c’ num1 aromatic_atom aromatic_atom aromatic_atom aromatic_atom_num
    double_aromatic_ring -> ’c’ num1 aromatic_atom aromatic_atom aromatic_atom ’c’ num ’n’ num1 aromatic_atom aromatic_atom aromatic_atom_num
    double_aromatic_ring -> ’c’ num1 aromatic_atom aromatic_atom aromatic_atom ’n’ num ’c’ num1 aromatic_atom aromatic_atom aromatic_atom_num
    aromatic_ring_6 -> starting_aromatic_c_num aromatic_atom full_aromatic_segment aromatic_atom aromatic_atom_num
    aromatic_ring_6 -> starting_aromatic_c_num full_aromatic_segment full_aromatic_segment aromatic_atom_num
    aromatic_ring_5 -> starting_aromatic_c_num aromatic_os full_aromatic_segment aromatic_atom_num
    aromatic_ring_5 -> starting_aromatic_c_num aromatic_atom aromatic_os aromatic_atom aromatic_atom_num
    aromatic_ring_5 -> starting_aromatic_c_num full_aromatic_segment aromatic_os aromatic_atom_num
    aromatic_ring_5 -> starting_aromatic_c_num full_aromatic_segment aromatic_atom aromatic_os_num
    aromatic_ring_5 -> starting_aromatic_c_num aromatic_atom full_aromatic_segment aromatic_os_num
    full_aromatic_segment -> aromatic_atom aromatic_atom
    aromatic_os -> side_aliphatic_ring
    side_aliphatic_ring -> ’c’ num ’(’ cycle_bond ’)’
    full_aromatic_segment -> side_aliphatic_ring_segment
    side_aliphatic_ring_segment -> ’c’ num ’c’ ’(’ cycle_bond ’)’
    side_aliphatic_ring_segment -> ’c’ ’(’ cycle_bond ’)’ ’c’ num
    valence_4_num -> ’C’ num
    valence_4_num -> ’[’ ’C’ ’@’ ’]’ num
    valence_4_num -> ’[’ ’C’ ’@’ ’@’ ’]’ num
    valence_4_num -> ’[’ ’N’ ’+’ ’]’ num
    valence_3_num -> ’[’ ’C’ ’@’ ’H’ ’]’ num
    valence_3_num -> ’[’ ’C’ ’@’ ’@’ ’H’ ’]’ num
    valence_3_num -> ’N’ num
    valence_3_num -> ’[’ ’N’ ’H’ ’+’ ’]’ num
    valence_3_num -> valence_4_num
    valence_3_num -> valence_4_num ’(’ nonH_bond ’)’
    valence_2_num -> ’O’ num
    valence_2_num -> ’S’ num
    valence_2_num -> ’S’ num ’(’ ’=’ ’O’ ’)’ ’(’ ’=’ ’O’ ’)’
    valence_2_num -> valence_3_num
    valence_2_num -> valence_3_num ’(’ nonH_bond ’)’
    valence_2_num -> valence_4_num ’(’ double_bond ’)’
    num -> ’1’
    num1 -> ’1’
    num -> ’2’
    num1 -> ’2’
    num -> ’3’
    num1 -> ’3’
    num -> ’4’
    num1 -> ’4’
    num -> ’5’
    num1 -> ’5’
    num -> ’6’
    num1 -> ’6’
    num -> ’7’
    num1 -> ’7’
    num -> ’8’
    num1 -> ’8’
    num -> ’9’
    num1 -> ’9’
    num -> ’%10’
    num1 -> ’%10’
    num -> ’%11’
    num1 -> ’%11’
    ...
    num -> ’%49’
    num1 -> ’%49’


Prospects for chargino pair production at CEPC
Jia-Rong Yuan 12, Hua-Jie Cheng 13,Xu-Ai Zhuang *1
1 Institute of High Energy Physics, Chinese Academy of Science, Yuquan Road
19B, Shijingshan District, Beijing 100049, China
2 University of Chinese Academy of Sciences, Yuquan Road 19A, Shijingshan
District, Beijing 100049, China
3 Department of Physics, National Taiwan University, Roosevelt Road 1 Sec. 4,
Taipei 10617, Taiwan
*<EMAIL_ADDRESS>(corresponding author)

## ABSTRACT

The proposed Circular Electron Positron Collider (CEPC), with a center-of-mass
energy $\sqrt{s}=240$ GeV, will primarily serve as a Higgs factory. At the
same time, it can offer good opportunities to search for new physics phenomena
at low energy, which can be challenging to discover at hadron colliders, but
well motivated by some theoretical models developed to explain, e.g., the
relic abundance of dark matter. This paper presents sensitivity studies of
chargino pair production, considering scenarios for both a Bino-like and a
Higgsino-like neutralino as the lightest supersymmetric particle, using full
Monte Carlo (MC) simulation. With the assumption of systematic uncertainties
at the level of 5%, the CEPC has the ability to discover chargino pair
production up to the kinematic limit of $\sqrt{s}/2$ for both scenarios. The
results have a minor dependence on the reconstruction model and detector
geometry. These results can also be considered as a reference and benchmark
for similar searches at other proposed electron-positron colliders, such as
the Future Circular Collider ee (FCC-ee) or the International Linear Collider
(ILC), particularly given the similar nature of the facilities, detectors,
center-of-mass energies, and target luminosities.

Keywords: CEPC; chargino; Bino; Higgsino

††footnotetext: * This study was supported by the National Key Programme
(Grant NO.: 2018YFA0404000).

###### Contents

  1. 1 Introduction
  2. 2 Detector, Software and Samples
  3. 3 Search for chargino pair production at CEPC
    1. 3.1 Search for chargino pair production with Bino LSP
    2. 3.2 Search for chargino pair production with Higgsino LSP
    3. 3.3 Chargino pair production search summary
  4. 4 Conclusion
  5. 5 Acknowledgements

## 1\. Introduction

Supersymmetry (SUSY) [1, 2, 3, 4, 5, 6, 7] predicts new particles, each of
whose spin differs by a half unit from their corresponding standard model (SM)
particles. In SUSY models with conserved $R$-parity [8], SUSY particles are
always produced in pairs, and the lightest supersymmetric particle (LSP) is
stable and can be considered as a potential dark matter candidate [9, 10].

The charginos $\tilde{\chi}^{\pm}_{i}$ $(i=1,2)$ and neutralinos
$\tilde{\chi}^{0}_{j}$ $(j=1,2,3,4)$ are referred to as electroweakinos. They
are the mass eigenstates formed from linear superpositions of the Bino, Wino
and Higgsino particles, the superpartners of the charged and neutral Higgs
bosons and electroweak gauge bosons. The subscripts i and j indicate states of
increasing mass. In the Minimal Supersymmetric Standard Model (MSSM), the
magnitude of the Bino, Wino, and Higgsino mass parameters are referred to as
$M_{1}$, $M_{2}$ and $\mu$ respectively. The mass splitting between the
electroweakinos mainly depends on the absolute values of $M_{1}$, $M_{2}$ and
$\mu$. Two SUSY scenarios are considered in this paper. In the first scenario,
the absolute values of $M_{1}$ and $M_{2}$ parameters are considered to be
near the weak scale and similar in magnitude, while the magnitude of $\mu$ is
significantly larger, such that $|M_{1}|<|M_{2}|\ll|\mu|$. In this case, the
LSP $\textstyle\tilde{\chi}^{0}_{1}$ is Bino-like
$\textstyle\tilde{\chi}^{0}_{1}$ and the next to lightest supersymmetric
particle (NLSP) is part of a Wino-like doublet forming
$\textstyle\tilde{\chi}^{0}_{2}$ and $\textstyle\tilde{\chi}^{\pm}_{1}$. The
second scenario considers the case that the absolute value of $\mu$ is near
the weak scale, while the magnitude of $M_{1}$ and $M_{2}$ can be
significantly larger, i.e. $|\mu|\ll|M_{1}|,|M_{2}|$, so that
$\textstyle\tilde{\chi}^{0}_{1}$, $\textstyle\tilde{\chi}^{\pm}_{1}$ and
$\textstyle\tilde{\chi}^{0}_{2}$ are Higgsino-like and almost mass degenerate.
The first scenario is favored by dark matter arguments and the second scenario
is motivated by naturalness considerations [11, 12, 13], but is a challenging
experimental signature to search for due to the small mass splittings between
NLSP and LSP.

Direct searches for chargino pair production with Bino-like or Higgsino-like
LSPs were performed previously at the Large Electron-Positron Collider (LEP)
and the Large Hadron Collider (LHC). This process was excluded by LEP for
$\textstyle\tilde{\chi}^{\pm}_{1}$ masses below 92.4 (91.9) GeV independently
of the $\textstyle\tilde{\chi}^{0}_{1}$ mass for the Higgsino (Bino) LSP case
[14, 15, 16, 17, 18, 19]. At ATLAS (CMS), chargino pair production assuming a
Bino-like LSP was excluded for $\textstyle\tilde{\chi}^{\pm}_{1}$ masses up to
420 (200) GeV assuming a W-boson-mediated chargino decay [20, 21]. For the
scenario assuming a Higgsino-like LSP, $\textstyle\tilde{\chi}^{\pm}_{1}$
masses below 193 GeV were excluded for mass splittings down to 4.6 GeV, and
mass splittings from 1.2 GeV to 30 GeV were excluded at the LEP bounds on
$\textstyle\tilde{\chi}^{\pm}_{1}$ mass by ATLAS [22]. Light higgsinos with
mass splittings below 1 GeV, preferred by naturalness conditions from the low-
energy fine-tuning measures, are still not excluded by the LHC and will be
very challenge for the HL-LHC too.

For the Bino-like (Higgsino-like) LSP case, the increased center-of-mass
energy of CEPC will extend the search sensitivity for
$\textstyle\tilde{\chi}^{\pm}_{1}$ ($\textstyle\tilde{\chi}^{\pm}_{1}$ and
$\textstyle\tilde{\chi}^{0}_{1}$) masses by more than 100 GeV compared to the
LEP experiment. With a cleaner collision environment and better reconstruction
and identification efficiency for low-energy particles [23], CEPC will have
excellent sensitivities to the very compressed electroweakino and sleptons
(superpartners of the leptons) search scenarios which are very difficult to
search for at the LHC and even HL-LHC.

The FCC-ee and ILC are proposed electron-positron colliders [24, 25]. The CEPC
and FCC-ee are both circular colliders designed to be built and operated in
several stages with center-of-mass energies from 90 GeV to 350 GeV. The ILC is
a linear collider with designed center-of-mass energies from 250 GeV to 1 TeV.
These colliders will primarily operate with center-of-mass energies at the ZH
production threshold, so an energy of 240 GeV is assumed for these studies. A
conservative systematic uncertainty of 5% is assumed in this paper and is
consistent with the LEP results [26, 27, 28]. The results presented here are
expected to be largely independent of the specific detector, trigger and data
acquisition choices made, and can be easily considered as a reference for
studies of the other two facilities with proper scaling of the luminosity.

This paper presents studies of the sensitivity to chargino pair production,
followed by the chargino decay to a W boson and a Bino-like or Higgsino-like
neutralino LSP, as illustrated in Figure 1. In order to simplify the analysis,
only leptonic W decays are considered. Events with two opposite-sign (OS)
muons and significant recoil mass ($M_{recoil}$, the invariant mass of the
recoil system against the two muons) are selected for all scenarios.

Figure 1: Representative diagram illustrating the pair production of charginos
and subsequent decay into a two-lepton final state via W bosons.

## 2\. Detector, Software and Samples

The CEPC Conceptual Design Report (CDR) supplies detailed information on the
detector and software [29]. Two kinds of CEPC detector concepts are proposed.
The baseline concept is a particle-flow oriented detector equipped with an
ultra-high granularity calorimeter, a low-material tracker and a 3 Tesla
solenoid. The alternative concept uses a dual readout calorimeter and a 2
Tesla solenoid. In this analysis, the CEPC baseline concept is used as the
detector model.

The following software was used for the study: Whizard 1.95 [30] is used to
generate the official SM MC simulated event samples. MadGraph 2.7.3 [31] and
Pythia 8.244 [32] are used for the generation of SUSY MC simulated event
samples. MokkaC 0.1.0 [33] is used to simulate the interactions between
particles and the detector. Clupatra 00-10 [34] is used to reconstruct tracks
from hits in the detector. The particle flow algorithm Arbor 3.4.2 [35] is
used to reconstruct physics object. LICH v4 [36] based on Multivariate Data
Analysis (TMVA) [37] is used for lepton identification.

The $\textstyle\tilde{\chi}^{\pm}_{1}$ are pair produced from electron-
positron collisions and each $\textstyle\tilde{\chi}^{\pm}_{1}$ decays into a
W boson and a $\textstyle\tilde{\chi}^{0}_{1}$ with a 100% branching ratio. No
other sparticle is considered in the production or decay. The signal samples
of chargino pair production in the Bino-like LSP case are parametrized as
function of the mass of the LSP and $\textstyle\tilde{\chi}^{\pm}_{1}$: the
lower bound on the chargino mass is set by the LEP limit, while the LSP mass
is bound by the mass difference between the chargino and the W boson. The
$\textstyle\tilde{\chi}^{\pm}_{1}$ mass is therefore varied in the range 90 -
119 GeV. The signal samples of chargino pair production in Higgsino-like LSP
case are parametrized as function of two SUSY parameters $\mu$ and tan$\beta$.
The $\mu$ varied in the range 90 - 118 GeV, while tan$\beta$ is varied in the
range 10 - 60. Samples containing $5\times 10^{5}$ ($1\times 10^{5}$) events
are simulated for each signal point with Bino-like (Higgsino-like) LSP.

Reference points with a $\textstyle\tilde{\chi}^{\pm}_{1}$ masses of 110 GeV
and with $\textstyle\tilde{\chi}^{0}_{1}$ masses of 1 GeV, 10 GeV or 25 GeV
are used for the Bino-like LSP scenario in this paper to illustrate typical
features of the SUSY models to which this analysis is sensitive. The
corresponding leading order (LO) cross section as computed by MadGraph is 2789
fb. For the Higgsino-like LSP scenario, three reference points with tan$\beta$
of 30 and $\mu$ of 90 GeV, 106 GeV or 118 GeV are used in this paper, with
theoretical cross sections at LO as computed by MadGraph of 1966 fb, 1522.9 fb
and 681.2 fb respectively.

This analysis only considers the dominant SM background processes with final
states with two leptons (electrons, muons or taus) and significant recoil
mass. The backgrounds are categorized into three types: Higgs processes, two
fermions backgrounds and four fermions backgrounds. The Higgs processes
considered are $\nu\nu H,H\to\tau\tau$. The two fermions backgrounds
considered are $\mu\mu$ and $\tau\tau$ processes. The four fermions
backgrounds include ZZ, WW, single Z, single W and Z or W mixing processes.
The samples are normalised to a luminosity of 5.05 ab-1. The cross sections of
the dominant background processes are shown in Table 1.

processes | cross section [fb]
---|---
$\tau\tau$ | 4374.94
$\nu\nu H,H\to\tau\tau$ | 3.07
$ZZ$ or $WW\to\tau\tau\nu\nu$ | 205.84
$ZZ\to\tau\tau\nu\nu$ | 9.2
$\nu Z,Z\to\tau\tau$ | 14.57
$ZZ$ or $WW\to\mu\mu\nu\nu$ | 214.81
$ZZ\to\mu\mu\nu\nu$ | 18.17
$WW\to\ell\ell$ | 392.96
$\nu Z,Z\to\mu\mu$ | 43.33
$\mu\mu$ | 4967.58
$e\nu W,W\to\mu\nu$ | 429.2
$e\nu W,W\to\tau\nu$ | 429.42
$eeZ,Z\to\nu\nu$ | 29.62
$eeZ,Z\to\nu\nu$ or $e\nu W,W\to e\nu$ | 249.34
Table 1:  Cross-sections of the dominant background processes in CEPC.

## 3\. Search for chargino pair production at CEPC

The event reconstruction consists of track reconstruction, particle flow
reconstruction and compound physics objects reconstruction. Tracks are
reconstructed from the hits in the detector by Clupatra [34]. Particle flow
reconstruction uses the tracks and the calorimeter hits to reconstruct single
particle physics objects. The output of particle flow reconstruction can be
used to reconstruct compound physics objects such as converted photons, taus
and jets. The identification efficiency of muons is 99.9% for energies above 2
GeV [36]. For muons with energies below 1.3 GeV and at the edge of the barrel
region or the overlap region of barrel and endcap, the identification
efficiency is lower than 90%. The recoil system consists of all the particles
except the two OS charged muons. Without considering the beam energy spread,
the resolution of the reconstructed recoil mass is between 300 and 400 MeV
[29].

(a) Bino LSP case (2 OS $\mu$, both energy $>$ 10 GeV) (b) Bino LSP case (2 OS
$\mu$, both energy $>$ 10 GeV) (c) Higgsino LSP case (2 OS $\mu$, both energy
$>$ 1 GeV) (d) Higgsino LSP case (2 OS $\mu$, both energy $>$ 1 GeV)

Figure 2: The kinematic distributions for (top) Bino-like and (bottom)
Higgsino-like LSP cases. The stacked histograms show the expected SM
backgrounds. The hatched bands represent the sum in quadrature of statistical
uncertainties of the total SM background. For illustration, the distributions
of three SUSY reference points (defined in Sec. 2) are also shown as dashed
lines. The first bin clearly includes the underflow events and the last bin
includes the overflow events.

The following variables are efficient to discriminate the signal events from
SM backgrounds:

  * •

$\Delta R(\mu^{+},\mu^{-})$, the angular distance between two muons.

  * •

$\Delta R(\mu^{\pm},recoil)$, the angular distance between the muons and the
recoil system.

  * •

$|\Delta\phi(\mu^{\pm},recoil)|$, the azimuthal distance between the muons and
the recoil system.

  * •

$|\Delta\phi(\mu^{+},\mu^{-})|$, the azimuthal distance between two muons.

  * •

$E_{\mu^{\pm}}$, the energy of the muons.

  * •

$P_{T}^{\mu^{\pm}}$, the transverse momentum of the muons.

  * •

$M_{recoil}$, the invariant mass of the two neutrinos and two LSPs.

The signal regions were optimized and defined based on the above kinematic
variables. Zn [38] was used as a sensitivity reference in the signal region
optimization, the defnition of which is shown in formula (1). The statistical
uncertainty and a 5% global systematic uncertainty were taken into account in
the Zn calculation.

$Zn=\left[2\left((s+b)\ln\left[\frac{(s+b)(b+\sigma_{b}^{2})}{b^{2}+(s+b)\sigma_{b}^{2}}\right]-\frac{b^{2}}{\sigma_{b}^{2}}\ln\left[1+\frac{\sigma_{b}^{2}s}{b(b+\sigma_{b}^{2})}\right]\right)\right]^{1/2}$
(1)

Several kinematic distributions after the two OS muon selection with $E_{\mu}$
larger than 10 GeV (1 GeV) are shown in Figure 2 (a)(b) (Figure 2 (c)(d)), and
show the good discrimination power between signal and SM processes.

### 3.1 Search for chargino pair production with Bino LSP

For the Bino-like LSP scenario, events containing exactly two OS charged muons
with energy larger than 10 GeV are selected. A selection on $\Delta
R(\mu^{+},\mu^{-})$ is used to reject background events from $\tau\tau$
process and $\mu\mu$ process. Events are required to have $P_{T}^{\mu^{\pm}}>$
30 GeV to suppress WW and Z$\nu$ processes. Most of the signal events have
large recoil mass according to the signal topology as shown in Figure 2 (a),
so $M_{recoil}>$ 130 GeV is applied to reject $\mu\mu$ and ZZ background
events and to increase the signal sensitivity. The signal region definition is
summarized in Table 2.

Signal Region
---
== 2 muons (OS)
$E_{\mu^{\pm}}>10$ GeV
$0.4<\Delta R(\mu^{+},\mu^{-})<1.6$
$P_{T}^{\mu^{\pm}}>$ 30 GeV
$M_{recoil}>$ 130 GeV
Table 2:  Summary of signal region selection requirements for the chargino
pair production with Bino LSP.

The kinematic distributions of $M_{recoil}$ and $P_{T}^{\mu^{-}}$, after
applying all signal region requirements except on the shown variable, are
given in Figure 3. The ”Others” component includes $\nu\nu H$, $H\to\tau\tau$,
$ZZ$ or $WW\to\tau\tau\nu\nu$, $ZZ\to\tau\tau\nu\nu$, $\nu Z,Z\to\tau\tau$,
$e\nu W,W\to\mu\nu$, $e\nu W,W\to\tau\nu$, $eeZ,Z\to\nu\nu$, $eeZ,Z\to\nu\nu$
or $e\nu W,W\to e\nu$ processes. The black line with arrow indicates the
signal region selection of this kinematic variable. The sensitivities shown in
the lower pad of the figures are obtained from the cumulative signal and
background events calculated in the direction of the cut arrow. The event
yields from background processes and signal reference points after signal
region requirements are shown in Table 3. The dominant background
contributions are from $ZZ$ or $WW\to\mu\mu\nu\nu$, $\mu\mu$ and
$WW\to\ell\ell$ processes. The expected sensitivities as function of the
$\textstyle\tilde{\chi}^{\pm}_{1}$ and $\textstyle\tilde{\chi}^{0}_{1}$ masses
with systematic uncertaintyies of 0% or 5% for the Bino-like LSP case are
shown in Figure 4. The discovery potential can reach up to the kinematic limit
of $\sqrt{s}/2$, and is not sensitive to systematic uncertainties.

(a) $M_{recoil}$ ($M_{recoil}$ >130 GeV) (b) $P_{T}^{\mu^{-}}$
($P_{T}^{\mu^{-}}$ >30 GeV)

Figure 3: “N-1” distributions after signal region requirements for the Bino-like LSP case. All signal region requirements are applied except on the variable shown. The stacked histograms show the expected SM backgrounds. The hatched bands represent the sum in quadrature of statistical uncertainties of the total SM background. For illustration, the distributions from the SUSY reference points (defined in Sec. 2) are also shown as dashed lines. The lower pad is the sensitivity Zn calculated with statistical uncertainty and 5% flat systematic uncertainty. The first bin clearly includes the underflow events and the last bin includes the overflow events. Process | Yields
---|---
$ZZorWW\to\mu\mu\nu\nu$ | 1632$\pm$42
$\mu\mu$ | 609$\pm$61
$WW\to\ell\ell$ | 163$\pm$13
$\tau\tau$ | 88$\pm$14
$\nu Z,Z\to\mu\mu$ | 47.9$\pm$7.3
$ZZ\to\mu\mu\nu\nu$ | 27.7$\pm$6.2
Others | 0.74$\pm$0.74
Total background | 2568 $\pm$ 77
m ( $\textstyle\tilde{\chi}^{\pm}_{1}$, $\textstyle\tilde{\chi}^{0}_{1}$) = ( 110 , 1 ) GeV | 5940 $\pm$ 130
m ( $\textstyle\tilde{\chi}^{\pm}_{1}$, $\textstyle\tilde{\chi}^{0}_{1}$) = ( 110 , 10 ) GeV | 6470 $\pm$ 140
m ( $\textstyle\tilde{\chi}^{\pm}_{1}$, $\textstyle\tilde{\chi}^{0}_{1}$) = ( 110 , 25 ) GeV | 8470 $\pm$ 160
Table 3:  The number of events with statistical uncertainties in the signal
region for signals and SM backgrounds for chargino pair production with Bino
LSP case.

(a) systematic uncertainty = 0% (b) systematic uncertainty = 5%

Figure 4: The expected sensitivities as a function of the
$\textstyle\tilde{\chi}^{\pm}_{1}$ and $\textstyle\tilde{\chi}^{0}_{1}$ masses
for chargino pair production in the Bino-like LSP scenario. Results are made
assuming a systematic uncertainty of (left) 0% or (right) 5%. Black numbers
represent expected Zn values.

### 3.2 Search for chargino pair production with Higgsino LSP

For the Higgsino-like LSP scenario, events containing exactly two OS muons
with energy above 1.0 GeV are selected. The cuts on $E_{\mu^{\pm}}$ and
$\Delta R(\mu^{\pm},recoil)$ are used to reject background events from
$\tau\tau$ processes and Z or W mixing processes. A requirement of
$|\Delta\phi(\mu^{\pm},recoil)|<2.9$ is used to suppress $\mu\mu$ process, and
$|\Delta\phi(\mu^{+},\mu^{-})|<1.4$ is used to suppress the background events
with two back-to-back muons, such as from the $\gamma\gamma\to\mu\mu$ process
where the photons are produced from the incoming electrons/positrons. Most of
the signal events have large recoil mass according to the signal topology as
shown in Figure 2 (c), so $M_{recoil}>237.5$ GeV has been applied to reject
$\mu\mu$, Z or W mixing and WW background events and to increase signal
sensitivity. The signal region definition for Higgsino-like LSP case is
summarized in Table 4.

Signal Region
---
== 2 muons (OS)
$E_{\mu^{\pm}}>1.0$ GeV
3.2 $<\Delta R(\mu^{\pm},recoil)<$ 4.6
$|\Delta\phi(\mu^{\pm},recoil)|<2.9$
$|\Delta\phi(\mu^{+},\mu^{-})|<1.4$
$M_{recoil}>237.5$ GeV
Table 4:  Summary of signal region selection requirements for the chargino
pair production with Higgsino LSP.

The kinematic distributions of $E_{\mu^{-}}$ and $M_{recoil}$, after applying
all signal region requirements except on the shown variable, are given in
Figure 5. The event yields from background processes and signal reference
points after signal region selections are in Table 5. The dominant background
contributions are from $\nu Z,Z\to\mu\mu$ and $\tau\tau$ processes. The
expected sensitivities as a function of the $\textstyle\tilde{\chi}^{\pm}_{1}$
mass and $\textstyle\tilde{\chi}^{0}_{1}$ mass with systematic uncertainty
assumptions of 0% or 5% for the Higgsino-like LSP case are shown in Figure 6.
The discovery sensitivity can again reach up to the kinematic limit of
$\sqrt{s}/2$, and is not sensitive to systematic uncertainties and mass
splitting between NLSP and LSP. Sensitivity is also expected at lower mass
splittings not explored here.

(a) $E_{\mu^{-}}$ ($E_{\mu^{-}}>$ 1.0 GeV) (b) $M_{recoil}$ ($M_{recoil}>$
237.5 GeV)

Figure 5: “N-1” distributions of used variables after signal region requirements for the Higgsino-like LSP case. All signal region requirements are applied except on the variable shown. The stacked histograms show the expected SM backgrounds. The hatched bands represent the sum in quadrature of statistical uncertainties of the total SM background. For illustration, the distributions from the SUSY reference points (defined in Sec. 2) are also shown as dashed lines. The lower pad is the Zn calculated with statistical uncertainty and 5% flat systematic uncertainty. The first bin includes the underflow events. The first bin clearly includes the underflow events and the last bin includes the overflow events. Processes | Yields
---|---
$\tau\tau$ | 107$\pm$16
$\mu\mu$ | 37$\pm$15
$\nu Z,Z\to\mu\mu$ | 27.8$\pm$5.6
$ZZ\to\mu\mu\nu\nu$ | 5.5$\pm$2.8
$ZZ$ or $WW\to\mu\mu\nu\nu$ | 3.2$\pm$1.8
$WW\to\ell\ell$ | 1.0$\pm$1.0
Others | 2.6$\pm$1.6
Total background | 183$\pm$23
( $\mu$ , tan$\beta$ ) = ( 90 GeV , 30 ) | 396$\pm$39
( $\mu$ , tan$\beta$ ) = ( 106 GeV , 30 ) | 267$\pm$28
( $\mu$ , tan$\beta$ ) = ( 118 GeV , 30 ) | 296$\pm$20
Table 5:  The number of expected events and statistical uncertainties in the
signal region for signals and SM backgrounds for chargino pair production in
the Higgsino-like LSP case.

(a) systematic uncertainty =0% (b) systematic uncertainty =5% (c) systematic
uncertainty =0% (d) systematic uncertainty =5%

Figure 6: The expected sensitivities as a function of the
$\textstyle\tilde{\chi}^{\pm}_{1}$ and $\textstyle\tilde{\chi}^{0}_{1}$ masses
for chargino pair production with the Higgsino-like LSP scenario. Results are
made assuming a systematic uncertainty of (left) 0% or (right) 5%. Black
numbers represent expected Zn values.

### 3.3 Chargino pair production search summary

The studies presented in this paper demonstrate that the discovery sensitivity
for electroweakinos with both a Bino-like LSP and a Higgsino-like LSP can
reach masses up to $\sqrt{s}/2$ at CEPC. The sensitivity can be further
improved if considering additional final states (i.e. electron and tau),
alternate signal scenarios such as disappearing track or chargino pair
production with photon initial-state radiation (ISR), or more complicate
analysis techniques (i.e. machine learning etc).

The prospective limits at CEPC compared with LEP and LHC results are shown in
Figure 7. In general, the sensitivity at CEPC is mainly constrained by the
center-of-mass energy and is not strongly dependent on the mass splitting
between NLSP and LSP, allowing CEPC to cover very compressed electroweakinos
that are very difficult to search for at the LHC due to the reduced acceptance
and reconstruction efficiency of the soft leptons.

Figure 7: Observed and expected exclusion limits on SUSY simplified models for
chargino-pair production with (a) Bino-like LSP and (b) Higgsino-like LSP
obtained by ATLAS. The observed limits obtained by LEP are also shown for
Higgsino LSP case in Figure 7 (b), while not shown in Figure 7 (a) since the
limit of chargino mass is below 100 GeV. The prospected limits at CEPC are
also shown in the dotted purple line for rough comparison.

## 4\. Conclusion

Prospective searches for chargino pair production via W boson decay in
scenarios with a Bino-like LSP or Higgsino-like LSP are performed at CEPC
using simulated samples. With the assumption of a 5% systematic uncertainty,
the discovery sensitivity in electroweakinos can reach up to $\sqrt{s}/2$. The
choice of systematic uncertainty has no large impact on the discovery
sensitivity. The sensitivity can be further improved by considering more final
states, additional signal scenarios or more complex analysis techniques. CEPC
can therefore extend the sensitivity to new Supersymmetry particles, and
possible dark matter candiates, of masses up to 120 GeV [14, 15, 16, 17, 18,
19]. For the Higgsino-like LSP scenario, CEPC has a very good sensitivity to
the compressed regions with a small mass splitting between
$\textstyle\tilde{\chi}^{\pm}_{1}$ and $\textstyle\tilde{\chi}^{0}_{1}$, which
for the ATLAS and CMS experiments are very difficult to reach [22]. The result
of this research is based on the CEPC expectations but can also be considered
as a reference for other lepton collider experiments at a center-of-mass
energy close to 240 GeV, such as ILC [24] and FCC-ee [25].

## 5\. Acknowledgements

The authors thank Cheng-dong Fu, Gang Li, and Xiang-hu Zhao for the simulation
code they provided.

## References

  * [1] Y.A. Golfand, E.P. Likhtman, JETP Lett. 13, 323 (1971)
  * [2] D.V. Volkov, V.P. Akulov, Phys. Lett. B 46, 109 (1973)
  * [3] J. Wess, B. Zumino, Nucl. Phys. B 70, 39 (1974)
  * [4] J. Wess, B. Zumino, Nucl. Phys. B 78, 1 (1974)
  * [5] S. Ferrara, B. Zumino, Nucl. Phys. B 79, 413 (1974)
  * [6] A. Salam, J.A. Strathdee, Phys. Lett. B 51, 353 (1974)
  * [7] S.P. Martin, Adv. Ser. Dir. High Energy Phys. 18, 1 (1998)
  * [8] G.R. Farrar, P. Fayet, Phys. Lett. B 76, 575 (1978)
  * [9] H. Goldberg, Phys. Rev. Lett. 50, 1419 (1983)
  * [10] J.R. Ellis, J.S. Hagelin, D.V. Nanopoulos, et al., Nucl. Phys. B 238, 453 (1984)
  * [11] R. Kitano, Y. Nomura, Phys. Lett. B 631, 58 (2005)
  * [12] R. Kitano, Y. Nomura, Phys. Rev. D 73, 095004 (2006)
  * [13] H. Baer, V. Barger, P. Huang, et al., Phys. Rev. D 87, 115028 (2013)
  * [14] LEPSUSYWG, ALEPH, DELPHI, L3 and OPAL experiments, note LEPSUSYWG/02-04.1.  URL http://lepsusy.web.cern.ch/lepsusy/Welcome.html
  * [15] ALEPH Collaboration, Phys. Lett. B 533, 223 (2002)
  * [16] The DELPHI Collaboration, Eur. Phys. J. C 31, 421 (2003)
  * [17] M. Chemarin, H. El-Mamouni, J. Fay, et al., Search For Gauge Mediated Susy Breaking Signatures At $\sqrt{s}$ = 192-208 Gev.  Tech. Rep. L3-Note-2707, CERN, Geneva (2001).  URL https://cds.cern.ch/record/571780
  * [18] OPAL Collaboration, Searches for Charginos with Small DeltaM in e+e- Collisions at sqrts = 189 - 202 GeV.  Tech. Rep. CERN-OPAL-PN-464, CERN, Geneva (2000).  URL https://cds.cern.ch/record/837835
  * [19] OPAL Collaboration, New Particle Searches in $e^{+}e^{-}$ Collisions at $\sqrt{s}$ = 200-209 GeV.  Tech. Rep. CERN-OPAL-PN-470, CERN, Geneva (2001).  URL https://cds.cern.ch/record/837791
  * [20] ATLAS Collaboration, Eur. Phys. J. C 80, 123 (2020)
  * [21] CMS Collaboration, JHEP 11, 079 (2018)
  * [22] ATLAS Collaboration, Phys. Rev. D 101, 052005 (2020)
  * [23] M. Ruan, H. Zhao, G. Li, et al., Eur. Phys. J. C 78(5), 426 (2018)
  * [24] T. Behnke, J.E. Brau, P.N. Burrows, et al., arXiv e-prints arXiv:1306.6329 (2013)
  * [25] M. Bicer, H. Duran Yildiz, I. Yildiz, et al., J. High Energy Phys. 01, 164 (2014)
  * [26] ALEPH Collaboration, Phys. Rev. B 526, 206 (2002)
  * [27] DELPH Collaboration, Eur. Phys. J. C 31, 421 (2003)
  * [28] OPAL Collaboration, Eur. Phys. J. C 32, 453 (2004)
  * [29] The CEPC Study Group, arXiv e-prints arXiv:1811.10545 (2018)
  * [30] W. Kilian, T. Ohl, J. Reuter, Eur. Phys. J. C 71, 1742 (2011)
  * [31] J. Alwall, R. Frederix, S. Frixione, et al., J. High Energy Phys. 07, 079 (2014)
  * [32] T. Sjöstrand, S. Ask, J.R. Christiansen, et al., Comput. Phys. Commun. 191, 159 (2015)
  * [33] P. Mora de Freitas, H. Videau.  Detector simulation with MOKKA / GEANT4: Present and future.  In Proceedings of International Workshop on Linear Colliders (LCWS 2002), (Korea: Jeju Island, 2002), p. 623-627
  * [34] F. Gaede, S. Aplin, R. Glattauer, et al., J Phys Conf Ser 513(2), 022011 (2014)
  * [35] M. Ruan, H. Zhao, G. Li, et al., Eur. Phys. J. C 78(5), 426 (2018)
  * [36] D. Yu, M. Ruan, V. Boudry, et al., Eur. Phys. J. C 77(9), 591 (2017)
  * [37] A. Hoecker, P. Speckmayer, J. Stelzer, et al., arXiv e-prints arXiv:physics/0703039 (2007)
  * [38] G. Cowan.  Discovery sensitivity for a counting experiment with background uncertainty.  http://www.pp.rhul.ac.uk/~cowan/stat/medsig/medsigNote.pdf, retrieved 9th April 2021