# Quantifying the performance of high-throughput directed evolution protocols

Adèle Dramé-Maigné  Anton Zadorin  Iaroslava Golovkova  Yannick Rondelez

###

Laboratoire Gulliver, CNRS, ESPCI Paris, PSL Research University, 10 rue
Vauquelin, Paris, France. * equally contributed to this work.

## Abstract

Most protocols for the high-throughput directed evolution of enzymes rely on
random encapsulation to link phenotype and genotype. In order to optimize
these approaches, or compare one to another, one needs a measure of their
performance at extracting the best variants. We introduce here a new metric
named the Selection Quality Index (SQI), which can be computed from a simple
mock experiment with a known initial fraction of active variants. As opposed
to previous approaches, our index integrates the random co-encapsulation of
entities in compartments and comes with a straightforward experimental
interpretation. We further show how this new metric can be used to extract
general trends of protocol efficiency, or reveal hidden mechanisms such as a
counterintuitive form of beneficial poisoning in the Compartmentalized Self-
Replication protocol.

## Introduction

Molecular Directed Evolution (DE) is a technique to obtain biomolecules with
desirable or improved functions, by iteratively generating pools of randomized
variants and screening or selecting them to extract the best performers from
these pools. This technique, initially developed for the search of nucleic
acid ligands[22, 12], has been extended to proteins and used to explore other
chemistries and functions [6, 2]. Most importantly, it is now also used to
select biopolymers with tailored catalytic properties, i.e. artificial
enzymes[3].

The first high-throughput enzymes DE experiments used living cells as a
selection medium for the mutant phenotypes. In such _complementation_
approaches, the cell expression machinery is employed to translate an
exogenous randomized genetic element into the corresponding mutated
polypeptide. Additionally, survival and replication of the host cells is made
dependent on the catalytic activity provided by the exogenous gene. Selection
between desirable and undesirable phenotypes becomes possible for two reasons:
first, the cell boundaries ensure that the genotype-phenotype linkage is
maintained, and that selection, acting on phenotypes, yields improved
transferable genotypes; second, cell transformation strategies provide a
straightforward way to obtain populations of cells where each member contain
one and only one of the possible mutant genotypes. Alternatively, when no
strong coupling between the exogenous phenotype and host cell survival can be
achieved, screening methods, where the cells can be sorted according to an
observable marker of phenotype, have been designed.

While these _in vivo_ implementations are very efficient, they still suffer
from a number of technical issues. First, complementation approaches are only
possible for catalytic activities that are related to the host cell’s
essential processes. In the screening protocols, permeability of the cell
walls to metabolite may interfere with the diffusion of substrate to the
catalyst. Finally, in both cases, the strategy has to be compatible with the
metabolic function of the living host; activities that are toxic for the host
gets naturally counter-selected, in some case disrupting the intended DE
process.

To bypass these problems, a number of recent developments have introduced the
idea of performing high-throughput directed enzyme evolution _in vitro_[14,
11, 15, 4, 18, 7, 23, 20, 1, 5, 13]. In these protocols, mutant genotype-
phenotypes are randomly distributed in artificial micro-compartments, with an
internal volume ranging between the femtoliter and the nanoliter. As is the
case for DE using living organisms carriers, two strategies can be used to
extract the best genotypes: screening and selection. In screenings, the
phenotypes induce a detectable modification of the physico-chemistry of the
compartment, which is used to sort the compartments one-by-one into desirable
and non-desirable bins[4]. In selection approaches, the artificial compartment
additionally enables a biochemical genetic amplification process, directed at
the gene of interest and whose activation is conditional on the desired
activity.

A common important feature of these high-throughput _in vitro_ protocols is
that they use a _random_ distribution of the mutant library in the artificial
micro-compartments (usually droplets, but also liposomes, micro-patterned
arrays[21], etc). By contrast with the case of _in vivo_ approaches, _in
vitro_ protocols provide no internal mechanism to insure that one and only one
genotype will be present in each compartment. As a consequence, one can only
control the _average_ number of mutants per compartment (a value referred as
$\lambda$ in the following discussion) but not their actual repartition in the
compartment population. If the partitioning process is truly random, the
distribution of occupancy will follow the Poisson law. For example, if
$\lambda$ is set to 1 (i.e. there are as many compartments as mutants
genotypes), the proportion of compartment containing 0, 1 and more than one
genotypes will be roughly one third each. Lowering the values of $\lambda$
yields more wasted (unoccupied) compartments, while higher values increases
the number of compartments containing multiple genotypes.

High-throughput protocols are therefore fundamentally stochastic, and their
purpose is to increase the proportion of highly functional variants after each
round. Their performance at enriching libraries will depend on a number of
experimental choices or constraints. These factors include the value of
$\lambda$, but also, for example, the contamination by parental genotype, the
homogeneity of the compartment’s physicochemical parameters, or the robustness
of the genotype-phenotype linkage, etc. Various experimental approaches have
been proposed, with different throughput and performance. Assessing the
relative efficiency of the methods and the particular contribution of these
factors on a given DE protocol can be challenging.

### Enrichment factor

In an attempt to characterize the quality of a given high-throughput directed
evolution protocols, many papers report a so-called enrichment factor
$\varepsilon=\frac{p^{\prime}}{1-p^{\prime}}\cdot\frac{1-p}{p}$, where $p$ and
$p^{\prime}$ are the frequencies of the allele of interest before and after
one selection cycle, respectively (in some cases the enrichment factor is
defined using the simple ratio of frequencies, $\varepsilon=p^{\prime}/p$). To
evaluate the value of this enrichment factor for a given protocol, an
experiment is set up where the initial library contains a small proportion of
active variants among many inactive ones, and this proportion is measured
again after a single selection round. Generally, to perform these experiments,
a mock library containing only two versions of the target protein, an active
and an inactive one, is used. A compilation of reported values of
$\varepsilon$ from various reports is given in Fig 1, while the full tables
are included in SI Material & Methods (Table S2, S3, S4).

Figure 1: Reported enrichment factors from a selection of references. The
enrichment factor is plotted in linear (right) or log (left) scale against the
initial fraction of active mutants used in the experiment. Lines connect data
points from the same manuscript, and obtained with a single value of $\lambda$
(indicated next to the line). Disks are used for protocol using screening and
polydisperse emulsions; diamonds, for protocols using screening and
monodisperse (microfluidic) emulsions; Stars, for protocols using selections
in polydisperse emulsion. Some values were recomputed from the reported date.
See SI for details.

However, the usage of $\varepsilon$ as a measure of quality comes with a
number of issues:

  * •

First, $\varepsilon$ depends strongly both on the initial fraction of
functional mutants $p$ and on the value of $\lambda$. As can be seen in Fig 1,
measurements of enrichment factors made with different values of $p$ generally
yield values of $\varepsilon$ that vary over many orders of magnitude.
Similarly, changes in $\lambda$ yield different enrichment factors, with no
clear interpretation.

  * •

Second, the upper theoretical limit of $\varepsilon$ is infinity when $p$ or
$\lambda$ become small. This reflects the fact that in principle, a flawless
selection or screening assay should get rid of the poor phenotype in a single
round. It is therefore difficult to estimate what should be a correct, or
acceptable, value of an experimental $\varepsilon$. In practice, smaller
values of $p$ or $\lambda$ tend to increase the observed value of
$\varepsilon$, suggesting, against common knowledge, that the protocol
performs better on less rich libraries. Additionally, measurements at very
small $\lambda$ are not really useful in assessing a protocol meant to provide
the highest possible throughput, as the throughput is actually proportional to
$\lambda$. Note that the theoretical lower limit of $\varepsilon$ is 0 at high
$p$ or high $\lambda$. Every value below 1 would indicate a process that
counter-selects for the good genotype.

Altogether, interpreting the enrichment factor is difficult and it is clear
that one cannot use it directly to compare different protocols in terms of
their ability to enrich libraries. We introduce here an improved performance
metrics, the Selection Quality Index (SQI). The SQI is expressed as a fraction
of the theoretical (achievable) value for the selection of a functional
variant from a pool of lethal variants and normalizes for the effect of $p$
and $\lambda$. An SQI of 1 denotes a perfect protocol, whereas an SQI of 0
indicates a complete failure, i.e. the absence of any enrichment effect
(negative values indicate counter-selection). Accordingly, this index -in
combination with the throughput of the protocol- can be used to assess the
true absolute potential of a given experimental approach to concentrate rare
catalytic variants. It can also be used to compare different protocols, or
different variants of the same protocol.

### Model

The factors that control the change in frequency of a mixture of active and
inactive mutants in a high-throughput protocol can be separated in two groups.

The first group contains controllable or known parameters of the experimental
design: $\lambda$, the sharing status, and the replication/selection function.
These effects are detailed below:

  * •

As mentioned above, when $\lambda$ increases, a fair fraction of the
compartments actually contain more than one phenotype, leading to the
possibility of “hitch-hiking”. Inactive mutants, which should not have
survived the screening or selection process, are carried over by sharing a
compartment with a better phenotype. In principle, knowledge of the
compartmentalization process allows us to describe the distribution of
occupancies. For example, in the case of random compartmentalization, this
distribution will follow a Poisson law.

  * •

The second effect, also linked to the random co-encapsulation of variants, is
fitness sharing. This occurs when the presence of the poor phenotype is
actually detrimental to the good one. This is typically the case in selection
experiments: an equal fraction of the replication potential associated with
one compartment will be used to amplify each local variant, whether or not
they actually contributed to the replication activity. This sharing effect
does not happen in screening protocols, which works by discarding the low
activity compartments, rather than replicating the genotypes in the good ones.
In other words, in screenings, a variant present in a supra-threshold
compartment contributes one genotype to the next generation irrespective of
the presence of co-encapsulated variants in his compartment. These two
situations are refereed as sharing and non-sharing, respectively, in the rest
of the manuscript.

  * •

Third, the replication/selection function[25], which describes how many
genotypic copies are pushed to the next generation, according to the
phenotypic activity observed in a compartment. In screenings for example, the
value of the function is one if the phenotypic activity is above a given
threshold, and 0 otherwise. In selections, it is set by the internal
replication chemistry.

The second group consists of experimental idiosyncrasies that are not directly
controllable. Together, they explain the observed deviation from a perfect
assay. These factors can be based on a variety of causes, some known or
guessable, other altogether unknown. For example, in screening test, the
sorting machinery may not be 100% efficient or generally makes a number of
errors (typically increasing when sorting is attempted at higher frequencies
[4]). Noise in expression can also play a role when compartments containing
the same genotype may end up with various level of phenotypic activity
(typical when bacterial expression systems is used, i.e. lysate assays [18]).
Additionally, leaks between compartments, poor genotype-to-phenotype linkage
(e.g. in mRNA display strategies), hidden selection biases (e.g. linked to
toxicity effect), non-homogeneity or polydispersity of the encapsulating
compartments (see below), carry-over of the parent genomes in selections, and
a myriad of other causes may affect the result. One typically does not have
enough information to precisely describe or model these effects, and is rather
interested in quantifying their collective negative impact on the selection
protocol, in order to empirically try to minimize it.

### The Selection Quality Index (SQI)

We have recently shown that all effects belonging to the first group can be
modeled using the single equation (1). This equation describes the evolution
of a population of variants submitted to random-compartmentalized selection
rounds, characterized by an additive replication/selection function $f$, a
sharing rule $\varphi(n)$ (the fraction of the value given by $f$ allocated as
copies to each individual in the droplet of size $n$), and a distribution of
phenotypes $\rho$ ([24]).

$\rho^{\prime}=\frac{{\displaystyle\sum_{n=0}^{\infty}\frac{\lambda^{n}}{n!}\varphi(n+1)\langle\delta_{x}*\rho^{*n},f\rangle}}{{\displaystyle\sum_{n=0}^{\infty}\frac{\lambda^{n}}{n!}\varphi(n+1)\langle\rho^{*n+1},f\rangle}}\,\rho$
(1)

where$*$ means convolution of distributions, $\delta_{x}$ is the delta-
function centered at $x$, and $\langle\rho,f\rangle$ means application of the
distribution $\rho$ to the function $f$.

Here, we are interested in quantifying the quality of a typically mock
selection protocol where only two variants are used. One of them is active and
the second completely dysfunctional (i.e., in isolation, it should not survive
the selection process). Therefore $\rho$ has only two discrete values.

Figure 2:  Graphical representation of the expected new fraction $p^{\prime}$
as a function of initial fraction $p$ for various replication/selection
function $f$ and for two values of $\lambda$. The exact analytical results
derived from (1) are shown by solid lines. The linear approximation is shown
by a dotted line for the cases both with and without sharing. $\Delta$
corresponds to the intercept of the curves with the vertical axis and is
independent of _f_.

Equation 1 allows us to compute the theoretical (expected) value of
$p^{\prime}$ as a function of $p$ for different $\lambda$, which we will call
$p^{\prime}_{theo}(p)$ (Fig. 2). Importantly, this function is discontinuous
at $p=0$, reflecting the effortless invasion of the population of
dysfunctional mutants by a rare functional one. In an ideal selection with no
co-encapsulation, the jump observed at low $p$, corresponding to the selection
of one active variant among an infinity of inactive ones would be equal to 1\.
Therefore, the value of this jump, $\Delta$, provides a good absolute
characterization of the selection efficiency at low $p$, and depends only on
$\lambda$. Indeed, one can show that the general equation(1) yields an
expression for $\Delta$, which does not depend on the replication function
$f$:

$\Delta=p^{\prime}_{theo}(p\rightarrow
0)=\frac{{\displaystyle\sum_{n=0}^{\infty}\frac{\lambda^{n}}{n!}\varphi(n+1)}}{{\displaystyle\sum_{n=0}^{\infty}\frac{\lambda^{n}}{n!}(n+1)\varphi(n+1)}}.$

However, it is generally not practical to perform experiments at very small
$p$, so we will also need the behavior of $p^{\prime}$ at $p\rightarrow 0$. As
this behavior does depend on $f$, we propose to correct using a linear
approximation of the following form.

$p^{\prime}_{theo}(p)\simeq\Delta+(1-\Delta)p.$ (2)

This expression is exact for linear replication functions and performs well
for nonlinear ones (especially, for not very large $\lambda$, see Fig. 2).
Altogether, we find that all possible scenario collapses into only two cases:

In the non-sharing situation, which is typical of screening protocols, we have
$\varphi(n)=1$ and the theoretical value of $\Delta$ as a function of
$\lambda$ is given by:

$\Delta=\frac{1}{1+\lambda}.$ (3)

In the sharing situation, which is more typical of selections assay, there is
a cost associated with the presence of co-encapsulated variants and
$\varphi(n)=1/n$. We obtain:

$\Delta=\frac{1-e^{-\lambda}}{\lambda}.$ (4)

We conclude that, knowing $\lambda$ and $p$, and the sharing behavior of a
given experiment setup, one can obtain a good approximation of the
theoretically expected gene frequency $p^{\prime}_{theo}(p)$ after one round
of selection of a very rare functional mutant in a population of dysfunctional
ones. This $p^{\prime}_{theo}$ can be used as a reference to assess the
quality of the experimental protocol, independently of the controllable
parameters $\lambda$ and $p$. In principle, one can thus use any convenient
value of $\lambda$ and $p$ and simply measure the experimental frequency
$p^{\prime}_{exp}$ after the selection cycle. The SQI of the protocol can then
be expressed as the ratio of the experimental to theoretical frequency jumps,
i.e:

$SQI=(p^{\prime}_{exp}-p)/(p^{\prime}_{theo}-p).$ (5)

Finally,

$SQI_{screen}=\frac{(1+\lambda)(p^{\prime}-p)}{1-p}\quad and\quad
SQI{}_{selection}=\frac{\lambda
e^{\lambda}(p^{\prime}-p)}{(e^{\lambda}-1)(1-p)}.$

These formula provide an absolute measurement of protocol quality and
selection efficiency, that corrects for $p$, $\lambda$ and the sharing
behavior. It is independent of the replication/selection function $f$, and
reflects only the experimental contingencies belonging to group 2.

### General trends for SQI values

In Fig. 3 we use the reported data shown in Fig. 1 to compute the SQI of a
variety of experimental protocols and conditions. This immediately reveals
some interesting trends.

First, some of the data series, taken at various but relatively high initial
fractions (e.g. for $p\geq 1/100$) collapse to a single SQI value. For example
the points from [Fallah-Araghi 2012] and [Beneyton 2014], now give roughly the
same SQI, independent of the initial fraction. Additionally, for these
relatively high initial fractions, we note a cluster of data points that
achieve an SQI close to one, indicating that the best achievable enrichment
performance has been attained.

Second, for the lower $p$ values, the selection efficiency always decrease
with decreasing initial fractions of active mutants. This stands in striking
contrast with the trend observed in Fig. 1, where enrichment factors seemed to
indicate that all assays performed much better at lower initial fraction. This
apparent trend was just a mechanical consequence of the fact that it is easier
to invade a population of dysfunctional mutants, rather than a population
already containing a significant fraction of functional mutants. Our analysis
corrects for that effect. It thus highlights the fact that other experimental
contingencies actually dominate at low $p$, making experimental protocols less
and less ideal as the initial fraction decreases. This is due to a decreasing
signal-to-background ratio, where the background noise can result from a
variety of effects, such as contamination by parent DNA in selections, sorting
errors in screenings, or the difficulty to recover very small amounts of
genetic material in the general case.

Third, protocols can be separated according to their type: only screening
using microfluidic compartmentalization appear to yield efficiencies close to
1. In contrast, protocols using screening of bulk emulsions typically have a
lower performance. This indicates that the size-dispersity of the compartments
is an important factor contributing to the efficiency of screenings. In
contrast, some selection processes are able to obtain high index without
resorting to monodisperse emulsion, which could be related to the superior
robustness of selections versus screenings, with respect to co-encapsulation
(already noted in [25]). Indeed, reports concerning selections tend to use
higher $\lambda$ values.

Besides these general cross-studies trends, one can also use the SQI to
analyze individual experimental setups. For example, In 2009, Baret and
colleagues tested the Fluorescent Activated Droplet Sorting method [4] with
different conditions, and concluded that poisoning co-encapsulation explained
the lower efficiency at higher $\lambda$. Re-analysis of these results shows
that the SQI, although correcting for $\lambda$, decreases clearly as well.
This suggests that another cause linked to the increased concentration of
cells, e.g. leaks, could explain the behavior.

Figure 3: SQI calculated for the data in Fig 1. Lines connect data point
originating from the same manuscript, and done with a single value of
$\lambda$ (indicated next to the line). Disks are used for protocol using
screening and bulk (non-monodisperse) emulsions; diamonds, for protocols using
screening and microfluidic (monodisperse) emulsion; Stars, for protocols using
selections in bulk emulsion. The dotted line at $=1$ represents the
theoretical maximum performance, once initial fraction and random partitioning
have been taken into account. Note that the CSR experiment[14] is off scale
here (see Fig4).

### The SQI reveals deviation to Poisson statistics

The SQI seems to be able to stratify experimental platforms in terms of both
the quality of the experimental design, and their reaction to increasingly
challenging libraries (that is, containing a smaller and smaller fraction of
active variants). One approach, however, seems to stand out in terms of its
SQI, which is much above 1 when measured for a high value of $\lambda=1.7$.
This approach is an _in vitro_ selection process, termed CSR, for
Compartmentalized Self-Replication[14]. It targets a bacterial library of
variants of the Taq polymerase gene, which is encapsulated in droplet with
Taq-specific primers. After lysis at high temperature and PCR thermal cycling,
mutants polymerases with higher activity replicate their own genes better than
less active ones. When the droplets are broken, the retrieved genetic
population is thus enriched in active Taq variants. According to its SQI, CSR
apparently selects active mutants better than the theoretical limit due to
random coencapsulation (Fig. 4). This could indicate that this particular
selection process (and possibly some other) has a build in mechanism that
limits the carryover of inactive variants. To investigate this point, and
since, the SQI in this case was evaluated from relatively sparse reported
information, we decided to re-implement the experiment.

After a number of failed attempts, we were able to reproduce the CSR with some
modifications. First, we used Klentaq[19], a lysate-robust polymerase variant,
instead of Taq. We found that his enzyme was expressed at sufficient levels in
KRX cells (while the initial report used TG1 cells). And finally we used
monodisperse droplets in fluorinated oils, instead of the polydisperse
emulsion in mineral oil initially reported (see SI experimental section).
Given these adjustments, we were indeed able to observe the self-replication
of the polymerase gene, both in bulk solution, and in 24 $\mu$m droplets. We
thus measured SQI at different values of $\lambda$, as shown in Fig 4. These
measurement showed that the CSR process is indeed able to perform optimally at
low $\lambda$ values, were the SQI is roughly 1. It also confirmed that higher
$\lambda$ result in SQI clearly above 1, and therefore that the CSR reaction
is somehow immune to random co-encapsulation.

Figure 4: Left: SQI calculated for CSR protocols at various $\lambda$. Red:
original CSR report[14], in polydisperse emulsion. Blue and green: this work,
two independent experimental replicates in monodisperse emulsions. The dotted
line at SQI $=1$ represents the theoretical maximum efficiency. Middle: Change
in expected frequency jump assuming that droplets containing more than n
bacteria are poisoned and do not participate in the reaction. Right: The SQI
is recomputed with a cutoff set to n = 2, corresponding to the experimental
observation that PCR is quickly poisoned by excess lysate.

In the course of the experiments, we also noted that the self-PCR reaction
performed less efficiently when higher concentration of bacteria were used in
the master mix. This is most likely associated with the strong toxicity of the
bacterial lysate toward the PCR reaction[17]. We thus reasoned that this
phenomenon could explain the observed behavior: if droplets containing more
than one bacteria are unable to support PCR, whatever the mixture of active
and inactive variants they contain, then they simply do not contribute to the
new generation. Therefore, it is as if co-encapsulation would not occur and
the jump in frequency should be less affected (Fig. 4 Middle). We therefore
performed the reaction, in test tubes, with a mixture of variants that
reproduce the content of droplet containing one active, plus increasing amount
of inactive variants (see SI experimental section). Indeed, we find that the
self-PCR reaction yield decreases a lot as soon as a single inactive variant
is encapsulated with an active wild-type. The reaction does not happen with
two or more encapsulated inactive bacteria. Fig. 4 right shows the recomputed
SQI taking this observation into account. In the case, the corrected SQI stays
close to 1 for all experiments, irrespective of the value of $\lambda$.

### Discussion

In this manuscript, we derive in a rigorous way a metric that can be used to
evaluate and compare experimental high-throughput directed evolution
protocols. Like previous approaches it requires a mock experiment, where a
know mixture of functional and dysfunctional variants is submitted to one
selection round, and the change in frequency is evaluated. However, contrary
to its predecessors, our metric separates experimental design in two
categories, depending on whether they are controllable design parameters or
not. We incorporate the value of initial fraction $p$ and the average
occupancy $\lambda$ in the calculation of the score, to provide a metric that
is independent on these factors and only reflects other, less controlled
causes that may affect the efficiency of a protocol.

When applied to reported data, we observe that the SQI indeed provide an
informed evaluation of protocol quality. The fact that it corrects for $p$ and
$\lambda$ allows us to extract trends by comparing different experimental
approach and make hypothesis about the cause for the observed differences in
efficiency. It also clearly reveals the effect of decreasing signal-to-
background values, when the initial fraction of active mutants becomes very
small. This approach allowed to detect abnormal selection efficiency, and link
this behavior to a deviation from the model of additive fitness, most probably
due to the toxicity of bacterial lysate. Counterintuitively, these toxic
effects have a positive effect of the selection efficiency.

The equations presented above are valid when individual mutants are _randomly_
distributed in the available compartments. In this case, one indeed obtains a
poissonian distribution of the number of mutant per compartment, on which we
have based our analysis. However, other distributions can be more relevant in
some cases. For example, some microfluidic devices use physical effects to
encapsulate objects in compartments with distribution sharper than the poisson
distribution[16, 8, 9, 10]. Alternatively, mutants candidates can have a
tendency to stick or aggregate to each other (for instance in pairs, triples,
tetrads, etc.), which will also distort the distribution away from the poisson
law. Polyploidy represent a related case, although closer in context to
population dynamics than directed evolution protocols. Fortunately, our
mathematical approach is general enough to handle many of these cases, and
yield theoretical values of $\Delta$ in closed analytical form, when
available.

The corresponding mathematical derivations are provided in SI. For example,
Fig 5 shows the effect of aggregated partitioning in pairs or triplets. When
$\lambda$ is small, we recover the intuitive result that the best possible
outcome in terms of frequency jump is the inverse of the aggregation cluster
size (i.e., in the case of $k$-clusters, the good variant can not invade more
than $1/k$ of the population in one round). However, the negative effect of
aggregation gets offset when $\lambda$ increases, as all curves gather on the
same asymptote, which depends only on the sharing behavior. Interestingly,
clustered selection tends to be more resilient to reasonable increases in
$\lambda$ than unclustered selection (Fig 5, left inset).

Figure 5: Right, the effect of aggregation of the mutants in pairs or triplets
on the theoretical value of $\Delta$ (the achievable frequency jump in one
round, valid for low $p$). Full lines correspond to sharing while dotted lines
indicate no sharing. Red; encapsulation of discrete individual; Blue,
encapsulation in pairs; green, encapsulation in triplets. Inset, the rescaled
curves $\Delta/\Delta_{\lambda=0}$. Left, the effect of polydispersity of the
compartment in the case of screenings. Here we assume that the volumes are
Gamma-distributed with mean one, and various shape parameter $\alpha$. The
corresponding distributions are shown in inset.

Many high throughput _in vitro_ DE protocols use emulsions as a convenient way
to provide compartments. Our bibliographic analysis has also highlighted the
role of the mono/poly-dispersity of these emulsions on the selection process.
At equal average compartment volume, polydisperse emulsion will have more
compartments containing multiple mutants. This is expected to decrease their
efficiency at selecting one particular phenotype from a mixture. Indeed, in
the present case of selection from a active/inactive mixture of variants, we
can show that the polydisperse emulsion always perform worse than the
monodisperse setting. This effect is shown for Gamma-distributed compartment
volumes, in the case of screening, in Fig. 5 right.

In conclusion, our analysis suggests that experimentalists using high
throughput methods to search for catalysts should put their efforts on
improving the efficiency of their protocols at very low initial fractions.
Indeed, many protocols target a throughput somewhere between $T=10^{6}$ and
$T=10^{9}$ and it seems unlikely that a very rare active mutant can be
purified in a single round. However, at low $p$, according to the equation 5,
the SQI can be written as

$SQI=p^{\prime}/p^{\prime}_{theo}(p)\sim p^{\prime}/\Delta.$

This means that knowing $\lambda$, and hence $\Delta$, the SQI directly
provides an estimation of the final fraction of active mutant (see 2). Many
(mostly microfluidic) approaches already provide a good SQI at $p>10^{-3}$,
implying immediate fixation in these cases. It is therefore the ability to
bring rare variants (with an initial fraction that can be estimated as the
inverse of the throughput of the protocols) to $p^{\prime}=10^{-3}$ that is
going to be critical for real applications. For example, if an SQI $>10^{-3}$
is obtained for $p$ as low as $1/T$, then one needs only two rounds to harness
the full potential of the experimental protocol, that is, fix a functional
variant initially present as a single copy in the library.

### Acknowledgements

This work was supported by the ERC Grant “ProFF” (number 647275).

## References

  * [1] Jeremy J. Agresti, Eugene Antipov, Adam R. Abate, Keunho Ahn, Amy C. Rowat, Jean-Christophe Baret, Manuel Marquez, Alexander M. Klibanov, Andrew D. Griffiths, and David A. Weitz.  Ultrahigh-throughput screening in drop-based microfluidics for directed evolution.  Proceedings of the National Academy of Sciences of the United States of America, 107(9):4004–4009, March 2010.  WOS:000275131100014.
  * [2] Frances H. Arnold.  Directed Evolution: Bringing New Chemistry to Life.  Angewandte Chemie International Edition, 57(16):4143–4148, April 2018.
  * [3] Frances H. Arnold and George Georgiou, editors.  Directed Enzyme Evolution: Screening and Selection Methods.  Humana Press, Totowa, N.J, 2003 edition edition, May 2003.
  * [4] Jean-Christophe Baret, Oliver J. Miller, Valerie Taly, MichaÃ«l Ryckelynck, Abdeslam El-Harrak, Lucas Frenz, Christian Rick, Michael L. Samuels, J. Brian Hutchison, Jeremy J. Agresti, Darren R. Link, David A. Weitz, and Andrew D. Griffiths.  Fluorescence-activated droplet sorting (FADS): efficient microfluidic cell sorting based on enzymatic activity.  Lab on a Chip, 9(13):1850–1858, July 2009.
  * [5] Thomas Beneyton, Faith Coldren, Jean-Christophe Baret, Andrew D. Griffiths, and Valerie Taly.  CotA laccase: high-throughput manipulation and analysis of recombinant enzyme libraries expressed in E. coli using droplet-based microfluidics.  Analyst, 139(13):3314–3323, 2014.  WOS:000337125600014.
  * [6] K Chen and F H Arnold.  Tuning the activity of an enzyme for unusual environments: sequential random mutagenesis of subtilisin E for catalysis in dimethylformamide.  Proceedings of the National Academy of Sciences of the United States of America, 90(12):5618–5622, June 1993.
  * [7] Pierre-Yves Colin, Anastasia Zinchenko, and Florian Hollfelder.  Enzyme engineering in biomimetic compartments.  Current Opinion in Structural Biology, 33:42–51, August 2015.  WOS:000365362400007.
  * [8] David J. Collins, Adrian Neild, Andrew deMello, Ai-Qun Liu, and Ye Ai.  The Poisson distribution and beyond: methods for microfluidic droplet production and single cell encapsulation.  Lab on a Chip, 15(17):3439–3459, August 2015.
  * [9] Dino Di Carlo, Daniel Irimia, Ronald G. Tompkins, and Mehmet Toner.  Continuous inertial focusing, ordering, and separation of particles in microchannels.  Proceedings of the National Academy of Sciences of the United States of America, 104(48):18892–18897, November 2007.
  * [10] Jon F. Edd, Dino Di Carlo, Katherine J. Humphry, Sarah KÃ¶ster, Daniel Irimia, David A. Weitz, and Mehmet Toner.  Controlled encapsulation of single-cells into monodisperse picolitre drops.  Lab on a Chip, 8(8):1262–1264, July 2008.
  * [11] Jared W Ellefson, Adam J Meyer, Randall A Hughes, Joe R Cannon, Jennifer S Brodbelt, and Andrew D Ellington.  Directed evolution of genetic parts and circuits by compartmentalized partnered replication.  Nature biotechnology, 32(1):97–101, January 2014.
  * [12] Andrew D. Ellington and Jack W. Szostak.  In vitro selection of RNA molecules that bind specific ligands.  Nature, 346(6287):818–822, August 1990.
  * [13] Ali Fallah-Araghi, Jean-Christophe Baret, Michael Ryckelynck, and Andrew D. Griffiths.  A completely in vitro ultrahigh-throughput droplet-based microfluidic screening system for protein engineering and directed evolution.  Lab on a Chip, 12(5):882–891, 2012.  WOS:000300047800006.
  * [14] F J Ghadessy, J L Ong, and P Holliger.  Directed evolution of polymerase function by compartmentalized self-replication.  Proceedings of the National Academy of Sciences of the United States of America, 98(8):4552–4557, April 2001.
  * [15] Andrew D Griffiths and Dan S Tawfik.  Directed evolution of an extremely fast phosphotriesterase by in vitro compartmentalization.  The EMBO journal, 22(1):24–35, January 2003.
  * [16] Evelien W. M. Kemna, Rogier M. Schoeman, Floor Wolbers, Istvan Vermes, David A. Weitz, and Albert van den Berg.  High-yield cell ordering and deterministic cell-in-droplet encapsulation using Dean flow in a curved microchannel.  Lab on a Chip, 12(16):2881–2887, July 2012.
  * [17] Milko B. Kermekchiev, Lyubka I. Kirilova, Erika E. Vail, and Wayne M. Barnes.  Mutants of Taq DNA polymerase resistant to PCR inhibitors allow DNA amplification from whole blood and crude soil samples.  Nucleic Acids Research, 37(5):e40, April 2009.
  * [18] Balint Kintses, Christopher Hein, Mark F. Mohamed, Martin Fischlechner, Fabienne Courtois, Celine Leine, and Florian Hollfelder.  Picoliter Cell Lysate Assays in Microfluidic Droplet Compartments for Directed Enzyme Evolution.  Chemistry & Biology, 19(8):1001–1009, August 2012.  WOS:000309152500012.
  * [19] H. Klenow and I. Henningsen.  Selective Elimination of the Exonuclease Activity of the Deoxyribonucleic Acid Polymerase from Escherichia coli B by Limited Proteolysis*.  Proceedings of the National Academy of Sciences of the United States of America, 65(1):168–175, January 1970.
  * [20] Tadas Povilaitis, Gediminas Alzbutas, Rasa Sukackaite, Juozas Siurkus, and Remigijus Skirgaila.  In vitro evolution of phi29 DNA polymerase using isothermal compartmentalized self replication technique.  Protein Engineering, Design and Selection, 29(12):617–628, December 2016.
  * [21] Yannick Rondelez, Guillaume Tresset, Kazuhito V. Tabata, Hideyuki Arata, Hiroyuki Fujita, Shoji Takeuchi, and Hiroyuki Noji.  Microfabricated arrays of femtoliter chambers allow single molecule enzymology.  Nature Biotechnology, 23(3):361–365, March 2005.
  * [22] C. Tuerk and L. Gold.  Systematic evolution of ligands by exponential enrichment: RNA ligands to bacteriophage T4 DNA polymerase.  Science, 249(4968):505–510, August 1990.
  * [23] Asao Yamauchi, Toshihiro Nakashima, Nobuhiko Tokuriki, Masato Hosokawa, Hideki Nogami, Shingo Arioka, Itaru Urabe, and Tetsuya Yomo.  Evolvability of random polypeptides through functional selection within a small library.  Protein Engineering, Design and Selection, 15(7):619–626, July 2002\.
  * [24] Anton S. Zadorin and Yannick Rondelez.  Natural selection in compartmentalized environment with reshuffling.  arXiv:1707.07461 [q-bio], July 2017.  arXiv: 1707.07461.
  * [25] Anton S. Zadorin and Yannick Rondelez.  Selection strategies for randomly distributed replicators.  arXiv:1711.04350 [q-bio], November 2017.  arXiv: 1711.04350.

See pages - of SI/Appendix.pdf See pages - of SI/Experimental_SI.pdf


#  Playing with Words at the National Library of Sweden -
Making a Swedish BERT

###### Abstract

This paper introduces the Swedish BERT (“KB-BERT”) developed by the KBLab for
data-driven research at the National Library of Sweden (KB). Building on
recent efforts to create transformer-based BERT models for languages other
than English, we explain how we used KB’s collections to create and train a
new language-specific BERT model for Swedish. We also present the results of
our model in comparison with existing models—chiefly that produced by the
Swedish Public Employment Service, Arbetsförmedlingen, and Google’s
multilingual M-BERT—where we demonstrate that KB-BERT outperforms these in a
range of NLP tasks from named entity recognition (NER) to part-of-speech
tagging (POS). Our discussion highlights the difficulties that continue to
exist given the lack of training data and testbeds for smaller languages like
Swedish. We release our model for further exploration and research here:
https://github.com/Kungbib/swedish-bert-models.

Keywords: BERT, Language Understanding, Swedish Language Models

languageresourceLanguage Resources

Playing with Words at the National Library of Sweden -
Making a Swedish BERT

Martin Malmsten, Love Börjeson, Chris Haffenden
---
KBLab, National Library of Sweden
Humlegården, Stockholm
www.kb.se/kb-labb
{martin.malmsten, love.borjeson<EMAIL_ADDRESS>

Abstract content

##  1\. Introduction

Recent advances in neural network architectures have enabled significant new
developments within natural language processing. With the appearance of the
first BERT model—Bidirectional Encoder Representations from
Transformers—Google’s researchers could point towards state-of-the-art
performance across a variety of NLP tasks in English [Devlin et al., 2019].
This was followed by the development of M-BERT, a multilingual model capable
of generalizing across languages [Pires et al., 2019, Wu and Dredze, 2019],
which in turn prompted various monolingual models beyond English—i.e.
CamemBERT and FinBERT—that have been shown to outperform M-BERT for deep
transfer learning of these languages [Martin et al., 2020, Virtanen et al.,
2019, Rönnqvist et al., 2019]. While demonstrating impressive performance on a
number of tasks, these latter models have needed to overcome the challenges of
a lack of data and testbeds for languages other than English.

Here we present the new Swedish model we have developed at the National
Library of Sweden (KB): KB-BERT. As with the likes of FinBERT, this has meant
finding solutions to the lack of data for a lower-resourced language. We
elaborate upon how our work at the recently-established KBLab for data-driven
research has allowed us to address this problem — at least in part (see
discussion below). Technical issues of access to data converge with wider
questions about the representativeness of language models: we have
deliberately used parts of the library’s collections with significant elements
of colloquial language, including the newspaper archive and government
reports—both containing large amounts of direct reported speech from a wide
spectrum of social actors—to create a BERT trained on a broad range of
Swedish. We have also included historical materials from the newspaper archive
dating back to 1945 to capture something approximate to the living language of
the national community [Anderson, 2016, Ahearn, 2017]. There has thus been a
significant democratizing of data in the training of KB-BERT.

This paper makes several contributions. Firstly, we explain how we used KB’s
collections to create the training data necessary for the making of our
Swedish BERT. Secondly, we describe and show the results from our testing
process, where we compare the performance of KB-BERT for various NLP tasks in
Swedish with that of existing models. Thirdly, we evaluate this performance
and highlight the particular steps for future research that we think this
raises.

##  2\.  Corpus disposition

The main challenge that we faced in developing KB-BERT—apart from the amount
of computational resources needed—was in compiling and cleaning a sufficiently
large and diverse corpus of Swedish text, cf. [Linder et al., 2020].
Fortunately, we have been able to draw upon KB’s vast textual resources, which
result from both large-scale digitization projects and the library’s role as
the national archive for electronic and traditional legal deposits.

To limit the corpus solely to modern Swedish, we have selected only resources
created from the 1940’s until late 2019. More precisely, we have included
material from the following categories:

  * •

Digitized newspapers—text extracted from OCR’d newspapers from KB’s archives.
This newspaper corpus contains a wide variety of text such as quotes, speech,
articles, short stories and adverts.

  * •

Official Reports of the Swedish Government—the collection of government
publications (Statens offentliga utredningar or SOU) contain around 5,000
volumes, ranging from 1940 to today. These texts are highly vetted before
publication and represent a large amount of correctly written Swedish. The
sentences tend to be longer than in other texts.

  * •

Legal e-deposits—resources deposited at KB under the e-deposit law (2012:
492). These include reports from Swedish authorities, e-books from publishers
and e-magazines.

  * •

Social media—comments collected from internet forums. This is by far the most
challenging corpus in terms of sentence structure and generally creative use
of the Swedish language. Punctuation is, for example, apparently optional in
such texts.

  * •

Swedish Wikipedia—all Wikipedia articles in Swedish.

Corpus | Words | Sentences | Size
---|---|---|---
allnews.txt | 2 997M | 226M | 16 783MB
sou.txt | 117M | 7M | 834MB
e-deposit.txt | 62M | 3.5M | 400MB
social.txt | 31M | 2.2M | 163MB
sw-wiki.txt | 29M | 2.1M | 161MB
Total | 3 497M | 260M | 18 341MB
Table 1: Size of cleaned corpora

As can be seen in the table above, there is an obvious skew towards newspaper
text. On the other hand, the allnews.txt corpus is the most diverse of the
corpora we have included. It is also our understanding that the amount of a
certain type of text is not proportional to the model’s subsequent
understanding of that type of text. In other words, incorporating smaller
corpora of different types of text, rather than simply adding more text of an
existing type, is still beneficial not only to the model’s representativity
but also its capacity and robustness as a whole.

The one notable omission is the comparatively low amount of fiction included
in the total corpus. This results from the fact that there is only a small
number of e-books in the combined corpus. The newspaper corpus has a certain
level of fictional texts, but the extent is unknown and it is mostly short
stories rather than full books. The legal e-deposit corpus does, however,
include a few hundred works of fiction.

##  3\. Text wrangling

The text in its raw form needed extensive cleaning and filtering before it
could be used to create pre-training data. Text preparation was carried out
using a combination of available and custom-made tools to deal with
discrepancies between and within corpora, as we describe below.

###  3.1. OCR quality fix

One major hurdle has been the presence of a few very specific OCR errors,
which most likely stem from an overly aggressive use of word lists, or a
faulty word list, in combination with the OCR engine confusing the letter m
with the letter r followed by n—and vice versa. This is especially problematic
for the word om, meaning about in Swedish, being OCRed as orn (which is not a
Swedish word) and then replaced by örn (meaning eagle in Swedish). Since it is
a very common word showing up in ca 82% of all newspaper pages, where more
than half of all occurrences were incorrectly OCRed on this point, the
resulting model would have ended up with a faulty understanding of Swedish.

This also became problematic with common endings in Swedish such as
“arna”/”erna”, which would be OCRed as “ama” / “ema”. For example, “tanterna”
[ “the old ladies” ] becomes “tantema” (not a Swedish word). To correct this
we developed a fix to identify candidates where m had been replaced with rn
and vice versa, by checking both versions against The Language Bank in
Gothenburg’s (Språkbanken) morphological thesaurus: Saldo [Borin et al.,
2013]. If the original word did not exist in Saldo, but the transformed one
did, the word was added as a candidate in a mapping file.

In the end, over 30,000 such candidates were identified—the majority appearing
only a few times each— and corrected. In the case of om – örn the word was
simply replaced, essentially deleting mentions of actual eagles from the
Swedish newspaper corpus. In total, 13 million occurrences were corrected.

###  3.2. Sentence- and paragraph splitting

The pre-training data creation process ultimately demands “paragraphs”. This
can mean quite different things, depending on the type of text at hand. In
newspaper text it is usually a block of text spanning part of a page, while in
monographs it is not uncommon for a paragraph to cover the whole page. In
newspapers it is also not uncommon for single words or sentences to show up as
paragraphs in the OCR. These are not valuable since they have little or no
context. To reduce noise, we therefore set a somewhat arbitrary lower limit
for a paragraph at ten sentences. This meant that we discarded quite a large
amount of paragraphs that were deemed too short. Sentences at the beginning
and end of paragraphs that did not seem to be whole sentences were also
discarded.

However, the major obstacle to creating paragraphs that we faced was in
splitting text into sentences. In the end, we used a custom script together
with spaCy[Honnibal and Montani, 2017] to try and recreate sentences from a
stream of text. This consisted of a pre-processing step to deal with
hyphenation—which is very common since Swedish has fairly long composite words
and columns of text can be narrow—followed by Spacy’s sentence splitter, and
lastly a post-process to deal with special cases not covered by spaCy. It is
worth noting that SpaCy does not yet have a statistical model for Swedish,
which is partly why so much custom scripting was required.

###  3.3. Emojis

A special case is that of the emojis that are occasionally used as punctuation
in texts on social media. For example, the phrase “Hej ☺Vad heter du?” [“Hi
☺What’s your name?”] is two sentences since the U+1F600 GRINNING FACE works as
end-of-sentence punctuation in this context. On the other hand, “Menvafan
☹☹☹varför gjorde du så!?” [“What the hell ☹☹☹why did you do that?”] is one
sentence even though it contains three emojis. We felt it important to retain
emojis both in the tokenization and in their use as sentence delimiters since
they are prevalent in some types of discourse, and if the model was to have
any chance of being used to analyze and/or understand conversations in social
media, this was necessary.

###  3.4. Tokenization

Once the text was cleaned and split into paragraphs of sentences, we used the
SentencePiece [Kudo and Richardson, 2018] library to create a tokenizer file
and then manually added emojis and skin color modifiers that had been filtered
out for some reason (Kudo and Richardson, 2018). Following the lead of both
the German [Branden et al., 2019] and Finnish BERT [Virtanen et al., 2019]
projects, we chose a dictionary size of approximately 50,000 tokens, with the
rationale that Swedish contains many compound words.

##  4\. Pre-training

We used the code and instructions provided in Google’s BERT repo to create
pre-training data and subsequently pre-train the model [Devlin et al., 2019].
We trained for one million steps with a max sequence length of 128 and a batch
size of 512. Then we trained for 100,000 steps with a msl of 512 and batch
size of 128. We continued training for approximately another million steps
with msl of 128 and batch size of 512 to see if this would yield better
results downstream.

Model pretraining was made partly in-house at the KBLab and partly, for
unproblematic material without active copyright, with generous support of
Cloud TPUs from Google’s TensorFlow Research Cloud (TFRC).

##  5\. Results and performance comparison

To evaluate the performance of the model and compare it with others we used
NER- and POS-tagging as a downstream task. These are both token classification
task that require some language understanding to perform well.

###  5.1. NER tagging

For this task, we used the Stockholm-Umeå Corpus 3.0 suc3 dataset hosted by
The Language Bank in Gothenburg (Språkbanken) converted to be used by a
modified version of the run_ner.py program included in the Huggingface[Wolf et
al., 2019] framework. The file was split into a training, test and evaluation
set, with 70% used for training, and 20% and 10% for test and evaluation
respectively.

The results show that our model, KB-BERT, outperforms existing BERTs either
trained for multilingual understanding by Google or specifically for Swedish
by Arbetsförmedlingen. Results for HFST-SweNER [Kokkinakis et al., 2014] are
also included in the table below, but a difference in evaluation method means
that the numbers are not necessarily comparable.

Tag | AF-AI | M-BERT | KB-BERT | HFST*
---|---|---|---|---
PER | 0.913 | 0.945 | 0.961 | 0.913
ORG | 0.780 | 0.834 | 0.884 | 0.534
LOC | 0.913 | 0.942 | 0.958 | 0.780
TME | 0.655 | 0.888 | 0.906 | –
MSR | 0.828 | 0.853 | 0.890 | –
WRK | 0.596 | 0.631 | 0.720 | 0.275
EVN | 0.716 | 0.792 | 0.834 | 0.513
OBJ | 0.710 | 0.761 | 0.770 | 0.437
AVG | 0.876 | 0.906 | 0.927 | –
Table 2: FB1 for NER tags

The model from Arbetsförmedlingen is trained using only a fraction of the data
used for KB-BERT, which shows in its limited performance.

###  5.2. Part-of-spech tagging

Using the SUC dataset we evaluated Part-of-speech tagging the same way as NER.
The results show a much smaller improvement, less than 1%, over M-BERT
compared to NER tagging. This is probably due to F1 being very high to begin
with since a lot of words have a static mapping to a corresponding POS tag.

Tag | KB-BERT | BERT-M |  Tag | KB-BERT | BERT-M
---|---|---|---|---|---
AB | 0.9495 | 0.9421 |  PC | 0.9526 | 0.9327
DT | 0.9905 | 0.9876 |  PL | 0.9131 | 0.8917
HA | 0.9274 | 0.9249 |  PM | 0.9785 | 0.9689
HD | 1.0000 | 1.0000 |  PN | 0.9899 | 0.9874
HP | 0.9727 | 0.9700 |  PP | 0.9923 | 0.9889
HS | 1.0000 | 1.0000 |  PS | 0.9976 | 0.9961
IE | 0.9992 | 0.9984 |  RG | 0.9696 | 0.9664
IN | 0.9574 | 0.9355 |  RO | 0.9610 | 0.9554
JJ | 0.9687 | 0.9579 |  SN | 0.9844 | 0.9804
KN | 0.9876 | 0.9860 |  UO | 0.8094 | 0.7197
NN | 0.9941 | 0.9896 |  VB | 0.9987 | 0.9954
AVG | 0.9863 | 0.9819 |  |  |
Table 3: FB1 for POS tags

###  5.3. Downstream tasks evaluation during pre-training

To get a sense of when to stop pre-training, evaluation of the NER task was
continually made at certain checkpoints as a measure of model performance.
Results show a fairly slow rise after an initial jump somewhere around ten
thousand steps with only marginal improvements after a few hundred thousand
steps. This is in line with the results from Deepset’s German BERT(Chan et
al., 2019).

Task | 10k | 50k | 150k | 350k | 700k | 1M | 2M
---|---|---|---|---|---|---|---
NER (F1) | 0.8687 | 0.912 | 0.918 | 0.926 | 0.925 | 0.923 | 0.927
Table 4: FB1 for NER during pre-training

###  5.4. Results for harder NER tasks

While our KB-BERT scores objectively high and outperforms existing models,
some specific tests designed to stress language understanding in downstream
tasks work only intermittently. For example, the latter part of the sentence
“Pelle och Kalle startar företaget Pelle och Kalle” [“Pelle and Kalle start
the company Pelle and Kalle”] gets correctly tagged as ORG in some runs but
not in others. This does not seem to stabilize with additional pre-training
above 1-2 million steps. However, there does seem to be a difference for more
complex NER tasks in the range between 350k and 1M pre-training steps. Perhaps
we are pushing the boundaries of the model, or perhaps more fine-tuning data
is needed. Further research is required to clarify this.

##  6\. Conclusion

In this paper, we have described how we used KB’s considerable textual
resources to train a new Swedish BERT model. Looking at the performance of KB-
BERT, we have demonstrated that it outperforms both Google’s multilingual
model, M-BERT, and existing Swedish models, such as Arbetsförmedlingen’s,
across a range of NLP tasks.

What conclusions can be drawn from this? Firstly, and most obviously, it
highlights the importance of data volume: the more data, the better the
performance. That a BERT model trained specifically for Swedish outperforms a
multilingual model is hardly surprising; this has been shown with comparable
results in other languages such as German and Finnish. However, the clear
improvement in language understanding of KB-BERT in relation to other Swedish
models is interesting, since it demonstrates the significance of the large
amount of data required if a model is to perform well.

Secondly, this study suggests the clear value of including a broad range of
language types in creating a truly robust model. Perhaps most surprising here
was that the use of apparently messy data, such as the formally unorthodox and
colloquial language used on social media that incorporates emojis, seems to
result in a BERT-model with a better capacity for understanding many different
types of language (not just specifically that used on social media). In
general terms, we believe that this supports our contention that a
democratically-inclined model trained on data derived from many social domains
leads to improved performance on NLP tasks. We intend to develop our analysis
of this issue further in future research at KBLab.

Thirdly, our work in this project has highlighted the continued problem of a
lack of testbeds for the Swedish language. In contrast to well-resourced
languages such as English, this lack presents particular problems for the
future development and testing of new Swedish language models. To help address
this issue, we are going to be participating in a new project aiming to
develop new and improved testbeds in Swedish over the coming year, in
collaboration with researchers at the Swedish Research Institute (RISE), the
Swedish Language Bank at Gothenburg University, and AI Innovation of Sweden
(“SuperLim: en svensk testmängd för språkmodeller”, funded by Vinnova,
2020-2021). Our hope is that this will lead to significantly improved
possibilities for the future evaluation of Swedish models.

##  7\. Acknowledgements

  * •

National Library of Sweden / KBLab

  * •

Resources from Stockholm University, Umeå University and Swedish Language Bank
at Gothenburg University were used when fine-tuning KB-BERT for NER.

  * •

Research supported with Cloud TPUs from Google’s TensorFlow Research Cloud
(TFRC)

  * •

Models are hosted on S3 by Huggingface

##  8\. Bibliographical References

## References

  * Ahearn, 2017 Ahearn, L. M.  (2017).  Living Language: An Introduction to Linguistic Anthropology. John Wiley & Sons, Chichester, second edition edition.
  * Anderson, 2016 Anderson, B.  (2016).  Imagined Communities: Reflections on the Origin and Spread of Nationalism.  Verso, London, revised edition edition.
  * Borin et al., 2013 Borin, L., Forsberg, M., and Lonngren, L.  (2013).  SALDO: a touch of yin to WordNet’s yang.  Language Resources and Evaluation, 47(4):1191–1211, December.  Place: Dordrecht Publisher: Springer WOS:000328213000011.
  * Branden et al., 2019 Branden, C., Möller, T., Pietsch, M., Soni, T., and Yeung, C. M.  (2019).  Open Sourcing German BERT: Insights into pre-training BERT from scratch.
  * Devlin et al., 2019 Devlin, J., Chang, M.-W., Lee, K., and Toutanova, K.  (2019).  BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding.  arXiv:1810.04805 [cs], May.  arXiv: 1810.04805.
  * Honnibal and Montani, 2017 Honnibal, M. and Montani, I.  (2017).  spaCy 2: Natural language understanding with Bloom embeddings, convolutional neural networks and incremental parsing.  To appear.
  * Kokkinakis et al., 2014 Kokkinakis, D., Niemi, J., Hardwick, S., Linden, K., and Borin, L.  (2014).  HFST-SweNER - A New NER Resource for Swedish.  European Language Resources Assoc-Elra, Paris.  Pages: 2537-2543 Publication Title: Lrec 2014 - Ninth International Conference on Language Resources and Evaluation WOS:000355611004025.
  * Kudo and Richardson, 2018 Kudo, T. and Richardson, J.  (2018).  Sentencepiece: A simple and language independent subword tokenizer and detokenizer for neural text processing.  CoRR, abs/1808.06226.
  * Linder et al., 2020 Linder, L., Jungo, M., Hennebert, J., Musat, C., and Fischer, A.  (2020).  Automatic Creation of Text Corpora for Low-Resource Languages from the Internet: The Case of Swiss German.  arXiv:1912.00159 [cs], June.  arXiv: 1912.00159.
  * Martin et al., 2020 Martin, L., Muller, B., Suárez, P. J. O., Dupont, Y., Romary, L., de la Clergerie, E. V., Seddah, D., and Sagot, B.  (2020).  CamemBERT: a Tasty French Language Model.  arXiv:1911.03894 [cs], May.  arXiv: 1911.03894.
  * Pires et al., 2019 Pires, T., Schlinger, E., and Garrette, D.  (2019).  How multilingual is Multilingual BERT?  arXiv:1906.01502 [cs], June.  arXiv: 1906.01502.
  * Rönnqvist et al., 2019 Rönnqvist, S., Kanerva, J., Salakoski, T., and Ginter, F.  (2019).  Is Multilingual BERT Fluent in Language Generation?  arXiv:1910.03806 [cs], October.  arXiv: 1910.03806.
  * Virtanen et al., 2019 Virtanen, A., Kanerva, J., Ilo, R., Luoma, J., Luotolahti, J., Salakoski, T., Ginter, F., and Pyysalo, S.  (2019).  Multilingual is not enough: BERT for Finnish.  arXiv:1912.07076 [cs], December.  arXiv: 1912.07076.
  * Wolf et al., 2019 Wolf, T., Debut, L., Sanh, V., Chaumond, J., Delangue, C., Moi, A., Cistac, P., Rault, T., Louf, R., Funtowicz, M., and Brew, J.  (2019).  Huggingface’s transformers: State-of-the-art natural language processing.  ArXiv, abs/1910.03771.
  * Wu and Dredze, 2019 Wu, S. and Dredze, M.  (2019).  Beto, Bentz, Becas: The Surprising Cross-Lingual Effectiveness of BERT.  arXiv:1904.09077 [cs], October.  arXiv: 1904.09077.

##  9\. Language Resource References

swedish_bert languageresource


# Numerical and quantum simulation of a quantum disentangled liquid

E. Abbasgholinejad  Department of Physics, Sharif University of Technology,
P.O.Box 11155-9161, Tehran, Iran  Department of Electrical and Computer
Engineering, University of Washington, Seattle, WA 98195, USA  S. Raeisi
Department of Physics, Sharif University of Technology, P.O.Box 11155-9161,
Tehran, Iran  email A. Langari  Department of Physics, Sharif University of
Technology, P.O.Box 11155-9161, Tehran, Iran<EMAIL_ADDRESS>
###### Abstract

The illustrative wave function for a quantum disentangled liquid (QDL)
composed of light and heavy particles is examined within numerical
simulations. Initial measurement on light particles gives rise to the volume
law of the entanglement entropy of the heavy particles subsystem. The entropy
reaches its maximum value as the ratio of the system to subsystem sizes
increases. The standard deviation of entanglement entropy from its
thermodynamic limit due to the initial configuration of the light particle is
diminished within ensemble averaging. We have introduced a quantum circuit to
simulate the underlying QDL state. The results of the quantum simulation are
in agreement with the numerical simulations which confirms that the introduced
circuit realizes a QDL state.

###### pacs:

75.10.Jm

##  I Introduction

Thermalization or the approach to thermal equilibrium in isolated many-body
quantum systems is a fundamental problem in quantum statistical physics
Deutsch (1991); Srednicki (1994); Rigol  _et al._ (2008). There are various
quantum systems, whose properties can be described on the basis of equilibrium
statistical mechanics, while there are also systems in which the thermal
equilibrium is not reached. For instance, integrable models or systems which
exhibit many-body localization Anderson (1958); Kramer and MacKinnon (1993);
Gornyi  _et al._ (2005); Basko D.M.  _et al._ (2006); Oganesyan and Huse
(2007); Pal and Huse (2010); Bauer and Nayak (2013); Imbrie (2016) due to
strong disorder, represent non-thermalized phase D’Alessio  _et al._ (2016);
Abanin and Papić (2017); Alet and Laflorencie (2018). Recent theoretical
studies have revealed that integrability and the existence of a static
disorder are not necessary conditions Schiulaz  _et al._ (2015); Yao  _et al._
(2016a); Smith  _et al._ (2017); Yarloo  _et al._ (2018, 2019) for the
violation of thermalization, while there are translationally invariant and
non-integrable systems that do not thermalize in the conventional sense, among
them is the Quantum Disentangled Liquid (QDL) Grover and Fisher (2014).

QDL is a finite temperature phase of translationally invariant quantum liquid
consisting of two (or more) species of indistinguishable particles with a
large mass ratio, namely heavy and light ones. The basic characteristic of the
QDL phase is that, while the heavy degrees of freedom are fully thermalized,
the light ones satisfy an area law of entanglement entropy for a typical fixed
configuration of the heavy particles. This suggests that the light particles
are “localized” by their heavy counterparts. Thus, in a QDL phase, thermal
equilibration is incomplete and it is not a fully ergodic phase, where we are
faced with the partial breakdown of thermalization in a translationally
invariant system.

The playground is a closed quantum system, which evolves in time and its
dynamics is given by Hamiltonian via the Schrödinger equation. The time
evolution operator is unitary, which guarantees no loss of information for the
whole state of the system. However, the underlying information is distributed
to all parts of the system according to the time evolution process, which
changes the entanglement between different parts of the system. In this
respect, part of a closed quantum system is immersed in a heat bath of the
rest of the system, which opens the opportunity to localize a subsystem while
the rest is in a thermalized phase. A qualitative diagnostic to identify the
QDL phase is the bipartite entanglement entropy after a projective measurement
of the heavy/light degrees of freedom. The entanglement entropy scales with
the volume of the system (volume law) in a thermalized phase, while it grows
at most proportional to the area of the system (area law) for a non-
thermalized phase. Proposals for QDL state have been introduced in bosonic
model Schiulaz  _et al._ (2015), spin-charge degrees of freedom in a Hubbard-
like model Garrison  _et al._ (2017) and also anyons of the Kitaev ladder
Yarloo  _et al._ (2018).

In this work, we examine the QDL state introduced in Ref.Grover and Fisher
(2014). We present a quantum simulation scheme for QDL state Cirac and Zoller
(2012); Georgescu  _et al._ (2014). Specifically, we introduce a quantum
circuit for implementing the QDL state on a quantum simulator. First, we use
numerical simulation to justify the QDL properties of the introduced state and
investigate its characteristics on a finite-size system. We observe that the
measurement of light particles leaves the heavy particles in a thermalized
phase, where the entanglement entropy obeys the volume law. However, in a
finite-size system entanglement entropy is attained for a large ratio of the
size of the heat bath to its subsystem. An ensemble average justifies that our
results do not depend on the initial configuration of light particles. The
inverse process, measurement of heavy particles, leads to zero entanglement
entropy of light particles, which verifies the area law for the light species.
We compare the results of numerical simulations with their counterparts on a
quantum simulator. Although the number of qubits in the quantum simulator is
limited, our results justify that the introduced quantum circuit truly produce
a QDL state. A QDL state can be considered as a pragmatic system for quantum
bits, which should be interesting to realize a quantum memory Grover and
Fisher (2014); Yao  _et al._ (2016b).

In the next section, we review the notion of QDL and establish the state
introduced in Ref.Grover and Fisher (2014). Then, our proposed quantum circuit
is introduced. The numerical and quantum simulations are presented in Sec.III.
Finally, we summarize and conclude our results.

##  II Quantum Disentangled Liquid

###  II.1 Illustrative wavefunction

Grover and Fisher introduced an illustrative wave function for a lattice model
consisting of two types of heavy and light particles, which represents the
state of quantum disentangled liquid Grover and Fisher (2014). The proposed
state can be understood in terms of the notion of many-body localization. In
an extreme limit, where the position of heavy particles are fixed without any
dynamics, the pattern of heavy particle positions can be considered as a
random background for the light ones. This randomness is a disordered space
for the light particles and could lead to the localization of light particles.
On the other hand, if we only think about heavy particles a random
configuration of heavy particles would be a result of thermalized state. The
system is composed of ${\cal N}$ unit cells, each containing a pair of heavy
and light particles as shown in Fig.1-(a). In the framework of occupation
number representation, the wavefunction is written as:

$\ket{\psi}_{QDL}=\sum_{\\{N_{j}\\},\\{n_{j}\\}}\psi(\\{N_{j}\\},\\{n_{j}\\})\,\ket{\\{N_{j}\\}}\ket{\\{n_{j}\\}},$
(1)

where $\\{N_{j}=0,1;j=1,\dots,{\cal N}\\}$ , $\\{n_{j}=0,1;j=1,\dots,{\cal
N}\\}$ are the occupation number of heavy and light particles, respectively.
The coefficients are given by the following expression

$\psi(\\{N_{j}\\},\\{n_{j}\\})=\psi(\\{N_{j}\\})\prod_{j=1}^{{\cal
N}}\frac{1}{\sqrt{2}}[\delta_{n_{j},0}+e^{i\pi N_{j}}\delta_{n_{j},1}].$ (2)

The summation in Eq. (1) is over all possible $4^{\cal N}$ configurations. The
first term on the right hand side of Eq.(2), namely random phase, is a
function of heavy particles alone,

$\psi(\\{N_{j}\\})=\frac{1}{\sqrt{2^{{\cal N}}}}sgn(\\{N_{j}\\}),$ (3)

where $sgn(\\{N_{j}\\})=\pm 1$ is a random sign, chosen independently for all
of $2^{\cal N}$ heavy particles configurations. Moreover, the second part of
Eq. (2), namely occupation phase, could be categorized into four different
configurations as explained below,

$\begin{cases}N_{j}=\mathbb{0},n_{j}=0\,:\,\delta_{n_{j},0}+e^{i\pi
N_{j}}\delta_{n_{j},1}=1,\\\
N_{j}=\mathbb{0},n_{j}=1\,:\,\delta_{n_{j},0}+e^{i\pi
N_{j}}\delta_{n_{j},1}=1,\\\
N_{j}=\mathbb{1},n_{j}=0\,:\,\delta_{n_{j},0}+e^{i\pi
N_{j}}\delta_{n_{j},1}=1,\\\
N_{j}=\mathbb{1},n_{j}=1\,:\,\delta_{n_{j},0}+e^{i\pi
N_{j}}\delta_{n_{j},1}=-1.\\\ \end{cases}$ (4)

As a result there is only a minus sign due to occupation of both heavy and
light particles, simultaneously.

Consider a system consisting of $2{\cal N}$ qubits, the first ${\cal N}$
qubits represent the heavy particles and the rest represents light particles
at different sites. As an example, for ${\cal N}=2$ the Hilbert space bases
are given below

$\displaystyle\\{\ket{\\{N_{j}\\}}\\}$
$\displaystyle=\\{\ket{\mathbb{00}},\ket{\mathbb{01}},\ket{\mathbb{10}},\ket{\mathbb{11}}\\},$
(5) $\displaystyle\\{\ket{\\{n_{j}\\}}\\}$
$\displaystyle=\\{\ket{00},\ket{01},\ket{10},\ket{11}\\}.$

Hence, using the configurations given in Eq. (4), the QDL wavefunction is
written in the following form

$\displaystyle\ket{\psi}_{QDL}=\frac{1}{2^{2}}$ $\displaystyle\Big{[}$ (6)
$\displaystyle sgn(\\{\mathbb{00}\\})\big{(}$
$\displaystyle\ket{\mathbb{00},00}+\ket{\mathbb{00},01}+\ket{\mathbb{00},10}+\ket{\mathbb{00},11}\big{)}$
$\displaystyle+sgn(\\{\mathbb{01}\\})\big{(}$
$\displaystyle\ket{\mathbb{01},00}-\ket{\mathbb{01},01}+\ket{\mathbb{01},10}-\ket{\mathbb{01},11}\big{)}$
$\displaystyle+sgn(\\{\mathbb{10}\\})\big{(}$
$\displaystyle\ket{\mathbb{10},00}+\ket{\mathbb{10},01}-\ket{\mathbb{10},10}-\ket{\mathbb{10},11}\big{)}$
$\displaystyle+sgn(\\{\mathbb{11}\\})\big{(}$
$\displaystyle\ket{\mathbb{11},00}-\ket{\mathbb{11},01}-\ket{\mathbb{11},10}+\ket{\mathbb{11},11}\big{)}\Big{]}.$

We emphasize that the first two qubits are heavy particles at each site and
the rest shows light particles. For example $\ket{\mathbb{01},01}$ means that
the first site has no particles and the second one is occupied with both a
heavy and a light particle, which leads to a minus sign as a consequence of
the second site.

---
(a)

---
(b)

---
(c)

---
(d)

Figure 1: (a) Schematic representation of Light (small circles) and heavy
(large circles) particles in the system. (b) Measurement is done on light
particles. (c) Tracing out a subsystem of heavy particles after measurement on
light particles. (d) Tracing out light particles in the initial system.

###  II.2 Quantum circuit

Figure 2: The six-qubit quantum circuit for preparing the QDL state with
${\cal N}=3$. The first three qubits, named ($q_{0}$, $q_{1}$, $q_{2}$),
represent the occupation of heavy particles in each site and the qubits
($q_{3}$, $q_{4}$, $q_{5}$) are used for light particles. All the qubits are
initialized in $\ket{0}$ state. Then, the Hadamard gates make a uniform
distribution of all the possible basis states. The next step is to set the
phases in Eq. (2). The phase is -1 when both the heavy and light sites are
occupied simultaneously. The qubits that are connected by controlled-Z gates
share the same site on the lattice. So when both sites are occupied, i.e.
$\mid\mathbb{1}1\rangle$, the states pick up a negative sign, otherwise, it
stays unchanged. The last step is to apply the random signs. The sign is
chosen independently for each of the $2^{\cal N}$ heavy particles
configuration. This part of the circuit assigns a random sign to the states of
the heavy sites. For more details on this part of the circuit, see Appendix A.

In this section we introduce an algorithm for preparing the wavefunction
explained in the previous section (Eq. (1)). First, we prepare a uniform
superposition of all the elements of the computational basis. Each element of
the superposition is affected by two sets of phases. The first one is the
phase that is set by the occupation of light and heavy particles on each site.
The second one is the random phase in Eq. 3.

In our algorithm, we start with all qubits initialized in $\ket{0}$. Then, we
apply Hadamard gates to all of them to prepare the uniform superposition of
all elements of the computational basis. This gives

$H^{\otimes 2\cal N}\ket{0}^{\otimes 2\cal N}=\frac{1}{2^{\cal
N}}\prod_{1}^{2\cal N}\big{(}\ket{0}+\ket{1}\big{)}=\frac{1}{2^{\cal
N}}\sum_{i=0}^{2^{2\cal N}-1}\ket{i}.$ (7)

Next we need to add the occupation phases. Specifically, for a given site,
when there is both a heavy and a light particle, the state picks up a negative
sign. This can be implemented using Controlled-Z gates between qubits
corresponding to the light and heavy particle of each site. This indicates
that we need $\cal N$ Controlled-Z gates for this stage of the algorithm.

For the last part, we need to add the random phase to each element of the
superposition. Specifically, we select a random array of $\pm 1$ of length $N$
and apply these random signs to different configurations of heavy particles.
Note that the phases only depend on the heavy particles and as a result, this
part should only act on the first $\cal N$ qubits. This can be done by a
diagonal unitary gate (in the computational basis) acting only on the qubits
corresponding to the heavy particles. The diagonal of the unitary operator
would be the random array of $\pm 1$ signs.

$U_{\vec{\theta}}=diag\\{e^{-i\theta_{0}},\cdots,e^{-i\theta_{m-1}}\\}$ (8)

To break this operator into single and two-qubit gates, we use the method
introduced by Schuch and Siewert in Schuch and Siewert (2003), which is
described fully in Appendix A. It utilizes the one-qubit gate, phase shift
(z-rotation), and the CNOT gate to build any arbitrary $m=2^{\cal N}$
dimensional controlled phase shift. Note that this prepares any random phase
and is not limited to $\pm 1$. The optimal circuit, at worst, is made up of
$2^{\cal N}-1$ one-qubit z-rotation gates and $2^{\cal N}-2$ CNOT. Therefore,
the depth of the circuit for this part grows exponentially with $\cal N$.

The quantum circuit for the full algorithm, for three sites, ${\cal N}=3$, is
shown in Fig.2. The occupation phases are implemented with the Controlled-Z
gates in part I and the random phases are implemented in part II.

##  III Numerical simulation of QDL state

Here, we examine the properties of a QDL state using numerical simulation at
finite system size, ${\cal N}$. The QDL state given in Eq.(1) is realized
numerically, assuming a random distribution of signs defined in Eq.(3). Next,
we will average over an ensemble. The details of the averaging will be further
discussed.

Firstly, we would like to see the scaling of entanglement under the
partitioning between the two components of particles. In this respect, the
light particles’ ($m$) degrees of freedom are integrated out from the density
matrix ($\rho_{Q}$), which is constructed from the pure state of Eq.(1). The
process is shown schematically in Fig.1-(d). The resulting reduced density
matrix ($\rho_{M}$) is given by

$\displaystyle\rho_{Q}$ $\displaystyle=$
$\displaystyle\ket{\psi}_{QDL\;\;QDL}\bra{\psi},$ (9) $\displaystyle\rho_{M}$
$\displaystyle=$ $\displaystyle\tr_{m}(\rho_{Q}).$ (10)

The Rényi entanglement entropies can be used to quantify the amount of
entanglement

$S(\alpha)=\frac{1}{1-\alpha}\ln\tr_{M}[(\rho_{M})^{\alpha}],$ (11)

where it reduces to the von Nuemann entropy in the limit of $\alpha\rightarrow
1$. Here, we use the second Rényi entropy as a measure of entanglement

$S_{M}=-\ln[\tr_{M}(\rho_{M}^{2})],$ (12)

and call it Rényi-2 for simplicity. The von Neumann entropy is greater than or
equal to the Rényi-2 entropy, hence, the volume law scaling of Rényi-2 entropy
results in the volume law of the von Neumann entropy. The numerical simulation
of the Rényi-2 entropy of the reduced density matrix $S_{M}$, is plotted in
Fig.3, for several values of system sizes. It shows a linear dependence on the
number of sites for the entanglement (see the dotted line in Fig.3). The
Rényi-2 entropy of heavy particles scales as the volume law, which indicates
that the underlying state belongs to an ergodic phase (for heavy particles),

$S_{M}=-\ln[\tr_{M}(\rho_{M}^{2})]={\cal N}\ \ln 2\propto L^{d}.$ (13)

Figure 3: Rényi-2 entropy of all heavy particles is plotted versus the number
of sites. It shows that the inter-component entanglement scales with the
volume of the whole system.

Now, we intend to see the effect of measurement on the QDL state. We split the
procedure into two subsections, measuring light and heavy particles.

Figure 4: Rényi-2 entropy of heavy particles after measurement of light
particles, versus system size (${\cal N}$) for A-subsystem with a fixed number
of sites $N_{A}=2$. At high ${\cal N}$, the heavy particles’ entropy reaches
its maximum,$S^{M/m}_{A}=2\ln 2$, manifesting the volume law, which is a key
feature of QDL. Ensemble average is taken into account for 200 initial
configurations of light particles.

###  III.1 Measurement of light particles

We start with the state of system given by the density matrix, $\rho_{Q}$. By
measuring the position of light particles, which has been shown schematically
in Fig.1-(b), the reduced density matrix of heavy particles is obtained,

$\rho(\\{\tilde{n}_{j}\\})=\bra{\tilde{n}}\rho_{Q}\ket{\tilde{n}},$ (14)

which is a function of the measured values of $\\{\tilde{n}_{j}\\}$. Then, we
consider a bi-partitioning of system into two regions A(subsystem) and
B(environment) with number of sites $N_{A}$ and $N_{B}$, respectively, where
$N_{A}+N_{B}={\cal N}$. The reduced density matrix of region $A$ is obtained
by tracing out heavy particles of region $B$:

$\rho_{A}(\\{\tilde{n}_{j}\\})=\tr_{B}\rho(\\{\tilde{n}_{j}\\}),$ (15)

which has been depicted in Fig.1-(c). The Rényi-2 entropy of subsystem A is
given by

$S_{A}^{M/m}=-\ln[\tr_{A}\big{(}\rho_{A}(\\{\tilde{n}_{j}\\})^{2}\big{)}].$
(16)

It should be emphasized that the resulting entropy ($S_{A}^{M/m}$) is still a
function of the measured position of light particles, which has been randomly
chosen for the initial QDL state. To get an impression of the ensemble average
of the QDL state, the obtained entropy ($S_{A}^{M/m}$) is averaged over the
random locations of light particles, which are assumed to come from a
homogeneous distribution. Our result, which has been averaged over 200 samples
of random configurations of light particles, is plotted in Fig.4. In this
figure, the normalized entropy of heavy particles is plotted versus the system
size (${\cal N}$), where the size of the A-subsystem is fixed to $N_{A}=2$. As
Fig.4 shows, the entropy reaches its maximum for large system size
(thermodynamic limit), i.e.

$S_{A}^{M/m}=N_{A}\ \ln 2\sim L_{A}^{d}.$ (17)

This is consistent with the theoretical expectations in Ref.Grover and Fisher,
2014. The above equation (Eq.(17)) represents a volume law scaling, which
justifies that measuring light particles does not affect the thermal behavior
of heavy ones. This is a necessary condition of a QDL state.

We have also examined how the size of the A-subsystem affects the entropy of
heavy particles. Fig.5 shows the normalized entropy of heavy particles versus
the whole system size, ${\cal N}$, where different plots show fixed values of
A-subsystem size, $N_{A}=2,3,4,5$. Again, the ensemble average is taken over
200 configurations of light particles. All plots confirm that upon increasing
system size (${\cal N}$), the entropy reaches its maximum value. However, for
a larger A-subsystem, the maximum entropy is attained for larger size of the
whole system, while the error bar of ensemble average is smaller for larger
$N_{A}$. The qualitative behavior of all plots is the same and shows similar
trends.

The error bar that arises from the ensemble average is measured by the
standard deviation ($\sigma$) of our results. For a fixed size of A-subsystem,
$N_{A}=2$, we have plotted the standard deviation versus system size, ${\cal
N}$, in Fig.6-(a). The top panel also shows that the error bar decreases with
the system size. This might be due to the exponential growth of Hilbert space
on the system size, which results in less sensitivity to an initial
configuration of light particles. We have considered three sampling numbers,
$50,100,and200$, for comparison in the top panel, which justifies that
increasing the number of sampling does not change the final results anymore.
Hence, most of our results in this work have been plotted for $200$ number of
samples. It has been justified in Fig.6-(b), where the standard deviation is
plotted versus the number of sampling. The bottom panel shows almost a
constant value for the error bar versus the number of sampling as far as the
number of sampling is greater than $100$.

Figure 5: Normalized Rényi-2 entropy of heavy particles after measurement of
light particles versus system size for several A-subsystem, $N_{A}=2,3,4,5$.
The average on the random initial configurations of light particles is done on
200 samples.

###  III.2 Measurement of heavy particles

In the next numerical experiment, we measure on the locations of heavy
particles and obtain the reduced density matrix of light particles,

$\rho(\\{\tilde{N}_{j}\\})=\bra{\tilde{N}}\rho_{Q}\ket{\tilde{N}}.$ (18)

Similar to the previous case, we then trace over part of light particles,
namely tracing over B-subsystem, to get the reduced density matrix of light
particles.

$\rho_{A}(\\{\tilde{N}_{j}\\})=\tr_{B}\rho(\\{\tilde{N}_{j}\\}),$ (19)

The corresponding Rényi-2 entropy of subsystem A is given as the following

$S_{A}^{m/M}=-\ln[\tr_{A}\big{(}\rho_{A}(\\{\tilde{N}_{j}\\})^{2}\big{)}].$
(20)

For instance, consider an example of $N_{A}=2$ and ${\cal N}\geq 2$. After
measuring the heavy particles the system will result in a pure state and the
above entropy is zero,

$S_{A}^{m/M}=0\sim L_{A}^{d-1},$ (21)

which justifies the area law for A-subsystem. It means that light particles
are localized by the heavy ones although the heavy ones are in a thermal
state.

(a)
---

(b)
---

Figure 6: (a) The standard deviation ($\sigma$) of ensemble average versus
system size, for three different numbers of sampling, namely: 50, 100, and
200. (b) The standard deviation is plotted versus the number of sampling,
which justifies reaching the reliable results for sampling larger than
$\sim$100.

###  III.3 Quantum simulation on IBM systems

The quantum circuit of the QDL state, which has been introduced in Sec.II.2,
is implemented on the IBM quantum simulator IBM-Quantum (2022). However, there
is a limitation due to the number of available qubits and the quantum volume
of the system. As a result, we have done the simulation on the IBM system up
to ${\cal N}=4$. The number of sampling on the quantum simulator is set to
$20$ due to time constraints. Consequently, we have only considered $N_{A}=2$
as the size of the A-subsystem. For measuring the Rényi-2 entropy of the
circuit’s output, we estimate the final density matrix using the state
tomography Smolin  _et al._ (2012). As a consequence, there will be an error
with respect to the classical simulation even if the state preparation,
quantum gates, and measurements work ideally. In this paper, we used ibmq-
oslo, which is one of the IBM quantum resources. In this system, the CNOT gate
has some error that its median is on the order of $O(10^{-}3)$. Similarly the
readout error is on order of $O(10^{-2})$.

The results of the quantum simulations are plotted in Fig.7. In Fig.7-(a) the
entropy of heavy particles $S^{M/m}$ is shown versus ${\cal N}=2,3,4$, where
the location of light particles has been measured first, similar to the
procedure explained in Sec.III.1. We have also plotted the results of ideal
numerical simulations for comparison. Both plots indicate the trend of volume
law to reach its maximum as ${\cal N}$ increases. However, we observe the
difference between numerical and quantum simulations, which is beyond the
error bar of random configuration of signs (ensemble averaging).

In order to resolve the difference between numerical and quantum simulations,
we obtained the the entropy of light particles $S^{m/M}$ versus ${\cal
N}=2,3,4$, where the position of heavy particles has been measured first,
which is plotted in Fig.7-(b). As we know the quantum state for any subsystem
of light particles is a pure state after measurement on heavy particles, so
the entropy of light particles should be zero in an ideal simulation. The
results of numerical simulation, as expected, shows zero for the entropy of
the light particles. However, the results of quantum simulation show a non-
zero value that indicates the quantum states of the quantum simulator are not
pure. Note that for a mixed (not pure) quantum state we get positive entropy
in the measurement. This indicates that the results of the quantum simulation
are offset due to the lack of purity and presence of noise. The plot for
$S^{m/M}$ in Fig.7-(b) reflects the combined state tomography and quantum
hardware error. For an ideal quantum simulator, we would expect to get the
numerical simulation results.

However, the trend of numerical and quantum simulations for $S^{M/m}$ is the
same with a difference of less than 10 percent (within the error bar margin).
It justifies that the proposed quantum circuit (Fig.2) truly captures the
properties of a quantum disentangled liquid state.

(a)
---

(b)
---

Figure 7: The (a)((b)) panel shows the normalized Rényi-2 entropy of heavy
(light) particles subsystem after the measurement of light (heavy) particles
is done, versus system size (${\cal N}$) for both numerical and quantum
simulations. The size of the A-subsystem is fixed $N_{A}=2$. Quantum
simulations are obtained by making measurements on IBM quantum systems, which
show an offset bias due to the noise and lack of purity in the quantum
simulator. (a) Upon increasing ${\cal N}$, the entropy approaches its maximum,
$S^{M/m}_{A}=2\ln 2$, resembling the volume law for heavy particles. (b) This
figure shows the entropy of the light particles. For pure states, is seen in
the simulations, we expect this to be zero. However, from the quantum
simulations we get non-zero values which indicates that the states are not
pure.

##  IV Summary and Conclusions

We have simulated the quantum disentangled liquid state of Eq. (1) to examine
its entropy versus system/subsystem sizes, both numerically and on a quantum
simulator. The QDL state is composed of light and heavy particles, where heavy
ones are thermalized in a closed quantum system, while light particles are
localized that do not thermalize. Our numerical simulation shows that the
reduced density matrix of heavy particles obeys a volume law, which justifies
their ergodic behavior. If a measurement on the light particles is done and
the resulting entropy corresponding to a subsystem of heavy particles is
calculated, it still shows volume law behavior. However, a measurement on the
heavy particles leads to zero entropy (area law) of a subsystem of light
particles. It justifies that light particles are localized due to the
measurement of heavy ones, while heavy particles are always in a thermalized
phase. We have shown that reaching the maximum entropy of heavy particles on a
finite system depends on the ratio of the system to subsystem sizes as it
grows. Moreover, the error bar due to the initial configuration of light
particles diminishes as the size of the subsystem becomes large. The average
on the initial random configuration converges to its ensemble average after
$\sim 100$ sampling. We have also introduced a quantum circuit to realize the
QDL state on a quantum simulator. Our results on the quantum simulator are in
agreement with our numerical simulation, which justifies that our circuit
truly demonstrates a QDL state. However, there is an offset bias between
numerical and quantum simulations, which is due to the noise in the quantum
simulator.

It will be interesting to explore an adiabatic quantum simulation of these
systems Biamonte  _et al._ (2011). One potential candidate could be
optomechanical arraysLudwig and Marquardt (2013). It has been proposed to use
these systems for quantum simulation of topological insulators and quench
dynamics Raeisi and Marquardt (2020). The optical and mechanical modes in
these optomechanical arrays provide two fundamentally different degrees of
freedom that can be used for simulation of the light and heavy particles.

It will also be interesting to explore further optimization of the quantum
circuit. In particular, parallelization has been previously used for reducing
the cost of quantum simulations Raeisi  _et al._ (2012). Similar ideas might
find application for such simulations.

##  V Acknowledgements

We acknowledge the use of IBM Quantum services for this work. The views
expressed are those of the authors, and do not reflect the official policy or
position of IBM or the IBM Quantum team.

##  Appendix A Random Phase Circuit

In this appendix we review the algorithm to break an arbitrary diagonal phase
unitary gate, introduced by Schuch and Siewert Schuch and Siewert (2003)

$U_{\vec{\theta}}=diag\\{e^{-i\theta_{0}},\cdots,e^{-i\theta_{2^{\cal
N}-1}}\\}.$ (22)

The building block gates of this method are one qubit z-rotation and two qubit
CNOT.

$R_{Z}(\phi)=\begin{pmatrix}e^{-i\phi/2}&0\\\
0&e^{i\phi/2}\end{pmatrix},\text{CNOT}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\
0&0&0&1\\\ 0&0&1&0\\\ \end{pmatrix}.$ (23)

Suppose we prepare $\cal N$ qubits at initial computational state
$\ket{\mathbb{x}}=\ket{x_{1},\cdots,x_{\cal N}}$, $x_{i}\in\\{0,1\\}$.
Moreover, consider a binary string of length $\cal N$,
$\mathbb{y}=\\{y_{1},\cdots,y_{\cal N}\\}$. Let divide the indices that have
value $1$ in the $y$ string, $\\{m|m\in\mathbb{N},y_{m}=1\\}$, with set size
$M$. Then, use $M-1$ number of CNOTs with target qubit $y_{m_{M}}$ and
$\\{y_{1},\cdots,y_{m_{M-1}}\\}$ for control qubits. Finally, if we apply a
z-rotation on qubit $m_{M}$ with angle $\phi_{\mathbb{x}}$ and reuse the CNOT
gates exactly the same as before, the final wave function will be as the
following:

$\ket{\mathbb{x}}\to
e^{-i\theta_{\mathbb{x}}}\ket{\mathbb{x}},\quad\theta_{\mathbb{x}}=\frac{1}{2}(-1)^{\mathbb{x}.\mathbb{y}}\phi,$
(24)

where the definition of $\mathbb{x}$ and $\mathbb{y}$ product is
$\mathbb{x}.\mathbb{y}=x_{1}.y_{1}\oplus\cdots\oplus x_{\cal N}.y_{\cal N}$.
If we do the same algorithm for all $2^{\cal N}$ binary strings of $y$, the
rotation angle $\theta$ is

$\theta_{\mathbb{x}}=\frac{1}{2}\sum_{y=0}^{2^{{\cal
N}-1}}(-1)^{\mathbb{x}.\mathbb{y}}\phi_{\mathbb{x}}.$ (25)

More generally, since Eq. 25 is true for every computational state
$\mathbb{x}$, we can write

$\vec{\theta}=\frac{1}{2}{\cal H}\vec{\phi},$ (26)

where ${\cal H}$ is the $\cal N$-bit Hadamard transformation, and ${\cal
H}^{2}=2^{\cal N}\mathbb{I}$. The rotation angles
$\vec{\phi}=(\phi_{0},\cdots,\phi_{2^{\cal N}-1})$ for phase shift gates are
related to the phases $\vec{\theta}=(\theta_{0},\cdots,\theta_{2^{\cal N}-1})$
as given by

$\vec{\phi}=\frac{1}{2^{\cal N}-1}{\cal H}\vec{\theta}.$ (27)

Therefore, using CNOT and z-rotation gates, with rotation values of Eq. 27, we
are able to break a random unitary phase shift gate in Eq. 22.

## References

  * Deutsch (1991) J. M. Deutsch, Phys. Rev. A 43, 2046 (1991).
  * Srednicki (1994) M. Srednicki, Phys. Rev. E 50, 888 (1994).
  * Rigol  _et al._ (2008) M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854–858 (2008), 10.1038/nature06838.
  * Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Kramer and MacKinnon (1993) B. Kramer and A. MacKinnon, Reports on Progress in Physics 56, 1469 (1993).
  * Gornyi  _et al._ (2005) I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev. Lett. 95, 206603 (2005).
  * Basko D.M.  _et al._ (2006) Basko D.M., Aleiner I.L., and Altshuler B.L., Annals of Physics 321, 1126–1205 (2006).
  * Oganesyan and Huse (2007) V. Oganesyan and D. A. Huse, Phys. Rev. B 75, 155111 (2007).
  * Pal and Huse (2010) A. Pal and D. A. Huse, Phys. Rev. B 82, 174411 (2010).
  * Bauer and Nayak (2013) B. Bauer and C. Nayak, Journal of Statistical Mechanics: Theory and Experiment 2013, P09005 (2013).
  * Imbrie (2016) J. Z. Imbrie, Journal of Statistical Physics 163, 998–1048 (2016).
  * D’Alessio  _et al._ (2016) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Advances in Physics 65, 239 (2016).
  * Abanin and Papić (2017) D. A. Abanin and Z. Papić, Annalen der Physik 529, 1700169 (2017).
  * Alet and Laflorencie (2018) F. Alet and N. Laflorencie, Comptes Rendus Physique 19, 498 (2018), quantum simulation / Simulation quantique.
  * Schiulaz  _et al._ (2015) M. Schiulaz, A. Silva, and M. Müller, Phys. Rev. B 91, 184202 (2015).
  * Yao  _et al._ (2016a) N. Y. Yao, C. R. Laumann, J. I. Cirac, M. D. Lukin, and J. E. Moore, Phys. Rev. Lett. 117, 240601 (2016a).
  * Smith  _et al._ (2017) A. Smith, J. Knolle, D. L. Kovrizhin, and R. Moessner, Phys. Rev. Lett. 118, 266601 (2017).
  * Yarloo  _et al._ (2018) H. Yarloo, A. Langari, and A. Vaezi, Phys. Rev. B 97, 054304 (2018).
  * Yarloo  _et al._ (2019) H. Yarloo, M. Mohseni-Rajaee, and A. Langari, Phys. Rev. B 99, 054403 (2019).
  * Grover and Fisher (2014) T. Grover and M. P. A. Fisher, Journal of Statistical Mechanics: Theory and Experiment 2014, P10010 (2014).
  * Garrison  _et al._ (2017) J. R. Garrison, R. V. Mishmash, and M. P. A. Fisher, Phys. Rev. B 95, 054204 (2017).
  * Cirac and Zoller (2012) J. I. Cirac and P. Zoller, Nature physics 8, 264 (2012).
  * Georgescu  _et al._ (2014) I. M. Georgescu, S. Ashhab, and F. Nori, Rev. Mod. Phys. 86, 153 (2014).
  * Yao  _et al._ (2016b) N. Y. Yao, C. R. Laumann, J. I. Cirac, M. D. Lukin, and J. E. Moore, Phys. Rev. Lett. 117, 240601 (2016b).
  * Schuch and Siewert (2003) N. Schuch and J. Siewert, Phys. Rev. Lett. 91, 027902 (2003).
  * IBM-Quantum (2022) IBM-Quantum, “https://quantum-computing.ibm.com/,”  (2022).
  * Smolin  _et al._ (2012) J. A. Smolin, J. M. Gambetta, and G. Smith, Physical Review Letters 108 (2012), 10.1103/physrevlett.108.070502.
  * Biamonte  _et al._ (2011) J. D. Biamonte, V. Bergholm, J. D. Whitfield, J. Fitzsimons, and A. Aspuru-Guzik, AIP Advances 1, 022126 (2011).
  * Ludwig and Marquardt (2013) M. Ludwig and F. Marquardt, Phys. Rev. Lett. 111, 073603 (2013).
  * Raeisi and Marquardt (2020) S. Raeisi and F. Marquardt, Phys. Rev. A 101, 023814 (2020).
  * Raeisi  _et al._ (2012) S. Raeisi, N. Wiebe, and B. C. Sanders, New Journal of Physics 14, 103017 (2012).