# Tellurium emission line in kilonova AT 2017gfo

Kenta Hotokezaka,1 Masaomi Tanaka,2 Daiji Kato,3,4 Gediminas Gaigalas5
1Research Center for the Early Universe, Graduate School of Science,
University of Tokyo, Bunkyo, Tokyo 113-0033, Japan
2Astronomical Institute, Tohoku University, Aoba, Sendai 980-8578, Japan
3National Institute for Fusion Science, 322-6 Oroshi-cho, Toki 509-5292, Japan
4Department of Advanced Energy Engineering Science, Kyushu University, Kasuga,
Fukuoka 816-8580, Japan
5Institute of Theoretical Physics and Astronomy, Vilnius University,
Saulėtekio Ave. 3, Vilnius, Lithuania E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The late-time spectra of the kilonova AT 2017gfo associated with GW170817
exhibit a strong emission line feature at $2.1\,{\rm\mu m}$. The line
structure develops with time and there is no apparent blue-shifted absorption
feature in the spectra, suggesting that this emission line feature is produced
by electron collision excitation. We attribute the emission line to a fine
structure line of Tellurium (Te) III, which is one of the most abundant
elements in the second r-process peak. By using a synthetic spectral modeling
including fine structure emission lines with the solar r-process abundance
pattern beyond the first r-process peak, i.e., atomic mass numbers $A\gtrsim
88$, we demonstrate that [Te III] $2.10\,\rm\mu m$ is indeed expected to be
the strongest emission line in the near infrared region. We estimate that the
required mass of Te III is $\sim 10^{-3}M_{\odot}$, corresponding to the
merger ejecta of $0.05M_{\odot}$, which is in agreement with the mass
estimated from the kilonova light curve.

###### keywords:

transients: neutron star mergers

††pubyear: 2022

##  1 Introduction

The origin of r-process elements is a long-standing problem in astrophysics
(Burbidge et al., 1957; Cameron, 1957). Neutron star mergers have been
considered as promising sites of r-process nucleosynthesis (Lattimer &
Schramm, 1974). A neutron star merger, GW170817, was accompanied by an uv-
optical-infrared counterpart, a kilonova (or macronova) AT 2017gfo, which
provides strong evidence that r-process nucleosynthesis occurs in neutron star
merger ejecta (see, e.g., Metzger, 2017; Nakar, 2020; Margutti & Chornock,
2021, for reviews).

A series of spectral data of the kilonova AT 2017gfo was obtained in the the
optical and near infrared bands from 0.5 to 10 days after the merger (Andreoni
et al., 2017; Chornock et al., 2017; Kasliwal et al., 2017; Pian et al., 2017;
Smartt et al., 2017; Tanvir et al., 2017; Troja et al., 2017). The kilonova AT
2017gfo is dominated by the photospheric emission at the early times. The
photospheric emission around a few days after the merger peaks in the near
infrared band, indicating the existence of lanthanides, which have strong
absorption at optical to near infrared wavelengths (Barnes & Kasen, 2013;
Kasen et al., 2013; Tanaka & Hotokezaka, 2013; Fontes et al., 2020; Tanaka et
al., 2020; Kawaguchi et al., 2018; Barnes et al., 2021). The early spectra
also exhibit several absorption structures including: (i) the $0.8\,{\rm\mu
m}$ feature attributed to Sr II or He I (Watson et al., 2019; Domoto et al.,
2021; Gillanders et al., 2022; Perego et al., 2022; Tarumi et al., 2023) and
(ii) the $1.3\,{\rm\mu m}$ and $1.5\,{\rm\mu m}$ features attributed to La III
and Ce III, respectively (Domoto et al., 2022). In addition to the elemental
identification, Sneppen et al. (2023) demonstrated that the spectra in the
photospheric phase are useful to study the geometry of the outer part of the
kilonova ejecta, $\gtrsim 0.2c$.

After the photospheric phase, kilonovae enter the nebula phase, where the
ejecta is heated by charged decay products of the radioactivity of r-process
nuclei and the heat is radiated through atomic emission lines. Examining
kilonova nebular spectra provides opportunities to identify atomic species
synthesized in the merger ejecta that may not appear as absorption lines
during the photospheric phase. For instance, Hotokezaka et al. (2022)
interpreted the detection of Spitzer (Villar et al., 2018; Kasliwal et al.,
2022) at $4.5\,{\rm\mu m}$ at $43$ and $74$ days after the merger as emission
lines of selenium (Se) or tungsten (W). In the early nebular phase, $\sim
10\,{\rm days}$, the infrared emission is of particular interest because the
absorption opacity due to atomic transitions is lower compared to the optical
region (e.g., Tanaka et al., 2020), and thus, the emission lines are expected
to appear as early as $\lesssim 10$ days. Most of infrared emission lines are
expected to arise from fine-structure transitions in the ground terms of heavy
elements, for which the line wavelengths and transition rates can be obtained
with reasonably high accuracy from the experimentally calibrated atomic energy
levels. Furthermore, such emission lines can be used to estimate the mass
distribution of the emitting ions from the emission line spectra. In fact, the
mass distributions of ions in type Ia supernova ejecta have been derived from
the infrared nebular spectra (e.g., Kwok et al., 2023; DerKacy et al., 2023)

In section 2, we study an emission line feature at $2.1\,{\rm\mu m}$ in the
kilonova AT 2017gfo spectra from 7.5 to 10.5 days. We attribute this line to a
fine-structure line of doubly ionized Tellurium (Te III, atomic number 52).
The Te III mass that is required to explain the observed data is estimated as
$\sim 10^{-3}M_{\odot}$. With a synthetic spectral modeling with the solar
r-process abundance pattern, we show that [Te III] $2.10\,{\rm\mu m}$ is the
strongest fine structure emission line in the near infrared region. In section
3, we conclude the results and discuss the uncertainties and implications.

##  2 Te III line in kilonova

The emission lines produced through radiative de-excitation of atoms emerge
from the optically thin region of the ejecta. The optical depth of the
kilonova ejecta with an expansion velocity of $v_{\rm ej}$ and a mass of
$M_{\rm ej}$ is

$\displaystyle\tau$ $\displaystyle\approx\frac{\kappa M_{\rm ej}}{4\pi(v_{\rm
ej}t)^{2}},$ (1) $\displaystyle\approx 1\left(\frac{\kappa}{1\,{\rm
cm^{2}g^{-1}}}\right)\left(\frac{M_{\rm
ej}}{0.05M_{\odot}}\right)\left(\frac{v_{\rm
ej}}{0.1c}\right)^{-2}\left(\frac{t}{10\,{\rm day}}\right)^{-2},$ (2)

where $\kappa$ is the opacity and $t$ is the time since merger. The opacity is
dominated by bound-bound transitions of heavy elements and depends on the
composition and wavelengths. Tanaka et al. (2020) show that the expansion
opacity decreases with wavelength, e.g., $\sim 10$ – $100\,{\rm cm^{2}g^{-1}}$
around $0.5\,{\rm\mu m}$ and $\lesssim 1\,{\rm cm^{2}g^{-1}}$ around
$2\,{\rm\mu m}$. Therefore, infrared emission lines are expected to emerge at
the earlier time than optical lines. With the ejecta parameters of AT 2017gfo,
we expect emission lines to dominate over the photospheric emission as early
as $\sim 10\,{\rm days}$ around $2\,{\rm\mu m}$.

Figure 1 shows the spectral series of the kilonova AT 2017gfo from $7.5$ to
$10.5$ days after the merger taken by X-shooter on the Very Large Telescope
(Pian et al., 2017). The observed spectra are composed of several line
features and a continuum component extending from the optical to near infrared
bands. We model the underlying continuum spectrum by blackbody radiation,
where the photospheric velocity and temperature for $7.5$–$10.5$ days are
$0.06$ – $0.08c$ and $1700$ – $2400\,{\rm K}$, respectively. The observed
spectra clearly show an emission line at $2.1\,{\rm\mu m}$ (see Gillanders et
al. 2023 for a detailed analysis). The expansion velocity of the line emitting
region is $\sim 0.07c$ derived from Doppler broadening of the line, which is
consistent with the picture where the emission line is produced outside the
photosphere. The line flux remains roughly constant with time while the
continuum flux declines, and thus, the line-to-continuum ratio increases from
$\sim 1$ at 7.5 days to $\sim 1.5$ at 10.5 days. This development of the
emission line without a blue-shifted absorption feature indicates that the
emission at $2.1\,{\rm\mu m}$ is a forbidden line driven by electron collision
rather than an emission line associated with an absorption line, e.g., a
P-Cygni line or a fluorescence line. The wavelength of the peak of the
emission line feature indeed coincides with a fine structure line, [Te III]
$2.10\,{\rm\mu m}$, arising from the transition between the ground level
${}^{3}{\rm P}_{0}$ and the first excited level ${}^{3}{\rm P}_{1}$. It is
worth noting that [Te III] $2.10\,{\rm\mu m}$ has been detected in planetary
nebulae (Madonna et al., 2018). Note that the transition between the ground
level ${}^{3}{\rm P}_{2}$ and the second excited level ${}^{3}{\rm P}_{1}$ of
Te I also produces an emission line at 2.1 ${\rm\mu m}$. As discussed later,
the contribution of Te I line is weaker than Te III line.

It may not be surprising that Te III produces the strongest emission lines
because Te is among the most abundant elements in the second r-process peak.
Figure 2 shows the mass fraction of each atom at 10 days after the merger.
Here we assume that the final abundance pattern matches the solar r-process
residual with atomic numbers $A\geq 88$ (Hotokezaka & Nakar, 2020), i.e., the
elements beyond the first r-process peak. With this assumption, the most
abundant element is Sr and the second most is Te at 10 days. Note also that
[Te III] $2.10\,{\rm\mu m}$ is particularly expected to be strong as long as
Te III is abundant outside the photosphere because this line is produced by
radiative decay of the first fine structure transition level, which is easily
excited by electron collision. For the iron peak elements, [Co III]
$11.89\,{\rm\mu m}$ and [Co II] $10.52\,{\rm\mu m}$ represent lines of the
same nature. Indeed, these are among the most prominent mid-IR lines observed
in SNe Ia and SN 1987A, respectively (Kwok et al., 2023; Wooden et al., 1993).

Figure 1: Spectral series of the kilonova AT 2017gfo 7.5–10.5 days after the
merger. The observed data were taken by X-shooter on VLT (Pian et al., 2017).
The synthetic spectra are composed of fine structure emission lines (dashed
curve) and a continuum (dotted curve), where the continuum emission is
approximate by a blackbody with temperatures $T_{\rm
BB}=2400,\,2100,\,1800,\,1700\,{\rm K}$ at 7.5, 8.5, 9.5 and 10.5 days,
respectively. The electron temperature is fixed to be $2000\,{\rm K}$. The
ejecta model assumes $M_{\rm ej}=0.05M_{\odot}$, $v_{\rm exp}=0.07c$, and
$n_{e}=10^{7}\,{\rm cm^{-3}}(t/9.5{\rm d})^{-3}$. We use ionization fractions
of $(Y^{+0},\,Y^{+1},\,Y^{+2},\,Y^{\geq+3})=(0.25\,\,0.4,\,0.25,\,0.1$) for
all the atomic species for simplicity. The composition is assumed to be the
solar r-process abundance pattern with $A\geq 88$ (figure 2). The shape of
each emission line is assumed to be a Gaussian profile with a broadening
parameter of $0.07c$. The distance to the source is set to $=40\,{\rm Mpc}$.
The wavelength of [Te I] $2.10\,{\rm\mu m}$ and [Te III] $2.10\,{\rm\mu m}$ is
shown as a vertical dashed line. Also shown as vertical lines are possibly
strong emission lines at 10.5 day. The gray shaded vertical regions depict the
wavelength ranges between the atmospheric windows. The wavelength of [Te I]
$2.10\,{\rm\mu m}$ is shown with an offset of $+0.04\,{\rm\mu m}$.

Let us first give an estimate of the amount of Te III from the observed line
flux assuming that the observed line flux is predominantly produced by Te III
and the ejecta is optically thin to the [Te III] $2.10\,{\rm\mu m}$ line. The
total line luminosity is given by

$\displaystyle L\sim h\nu_{10}A_{10}f_{1}N({\rm Te\,III}),$ (3)

where $h\nu_{10}$, $A_{10}\approx 2\,{\rm s^{-1}}$, and $f_{1}$ are the
excitation energy, the radiative decay rate, and the fraction of Te III ions
in the ${}^{3}{\rm P}_{1}$ level, respectively, and $N({\rm Te\,III})$ is the
total number of Te III ions in the ejecta (see equations (1) and (2) in
Hotokezaka et al., 2022, for the formula of M1 transition probabilities). The
observed flux at $2.1\,{\rm\mu m}$ after subtracting the underlying continuum
is $\sim 0.1\,{\rm mJy}$, corresponding to the observed line luminosity of
$L_{\rm obs,line}\sim 2\cdot 10^{39}\,{\rm erg\,s^{-1}}$ with $D=40\,{\rm
Mpc}$. Equation (3) leads to a total mass of Te III:

$\displaystyle M({\rm Te~{}III})\sim
10^{-3}M_{\odot}\left(\frac{f_{1}}{0.1}\right)^{-1}\left(\frac{L_{\rm
obs,line}}{2\cdot 10^{39}~{}{\rm erg~{}s^{-1}}}\right).$ (4)

Note that the electron density at $t\sim 10.5$ day may be comparable to the
critical density of Te III ${}^{3}{\rm P}_{1}$ (Madonna et al., 2018), and
therefore, the level fraction $f_{1}$ is comparable to or slightly less than
that expected from the thermal distribution, i.e., $f_{1}\approx 0.1$ in
thermal equilibrium at $T_{e}=2000\,{\rm K}$. Given the total ejecta mass of
$\sim 0.05M_{\odot}$, we estimate that the mass fraction of Te is greater than
a few per cent.

Figure 2: Elemental mass fraction used in the synthetic modeling. The
abundance pattern at 10 days is determined such that the final abundance
pattern at $\sim 5\,{\rm Gyr}$ matches the solar r-process residual for
$88\leq A\leq 209$. Tellurium is the second most abundant element in this
model.

We now turn to the comparison of the observed spectra with a synthetic
spectral model. The synthetic spectrum is composed of fine structure emission
lines and a continuum component, where the continuum emission is approximated
by blackbody radiation. The blackbody temperature and radius at a given epoch
are determined such that the synthetic spectrum roughly agrees with the
observed one at the near IR region $\lesssim 2\,{\rm\mu m}$. The emission line
spectrum is computed by the one-zone modeling presented in Hotokezaka et al.
(2022), where the energy level populations are solved by balance between
collision and radiative decay for a given electron density, temperature, and
ionization state. We use the collision strengths of the fine structure
transitions of the ground term of Te III derived by Madonna et al. (2018). The
collision strengths of other elements that are relevant for the nebular
modeling at $\lambda\lesssim 3.5\,{\rm\mu m}$ are obtained by using an atomic
structure code Hullac (Bar-Shalom et al. 2001, see also Hotokezaka et al.
2022) and the M1 line list is constructed by using the NIST database (Kramida
et al., 2021) and the LS selection rules with the single-configuration
approximation (Hotokezaka et al. in prep). Note that the wavelengths and
radiative transition rates of the M1 lines in the list are sufficiently
accurate for our purpose.

In the modeling, the ejecta composition is assumed to be the solar r-process
abundance pattern with $A\geq 88$ (figure 2), which is the same as the second
and third peak model used in Hotokezaka et al. (2022). The model also assumes
the electron temperature, $T_{e}=2000\,{\rm K}$ and the ionization fractions
$(Y^{+0},Y^{+1},Y^{+2},Y^{\geq+3})=(0.25,0.4,0.25,0.1)$111We neglect the
emission lines of ions in $Y^{\geq+3}$.. These quantities are assumed to be
constant with time for simplicity. This choice of the ionization fraction is
somewhat motivated by Hotokezaka et al. (2021), where the ionization fractions
of Nd atoms in the kilonova nebular phase are studied. With these parameters,
the ejecta mass of $0.05M_{\odot}$ and the expansion velocity of $0.07c$, [Te
III] $2.10\,{\rm\mu m}$ is the strongest emission line among M1 transitions of
all the heavy elements beyond the first r-process peak and the synthetic
spectra can roughly reproduce the emission line structure around $2.1\,{\rm\mu
m}$. However, one should keep in mind that the ionization fraction varies
among different atomic species. The ejecta mass of $0.05M_{\odot}$ agrees with
the ejecta mass estimated from the energy budget of the bolometric light curve
(Waxman et al., 2019; Kasen & Barnes, 2019; Hotokezaka & Nakar, 2020).

If this interpretation is correct, we expect that the $2.1\,{\rm\mu m}$ line
remains at the later times while the continuum flux keeps declining. It is
worth noting that Te III may produce another emission line at $2.93\,{\rm\mu
m}$ arising from the transition between the first and second excited levels
(${}^{3}{\rm P}_{1}$-${}^{3}{\rm P}_{2}$) at the later times because the
electron temperature is expected to gradually increase with time (Hotokezaka
et al., 2021; Pognan et al., 2022a). Although this line may be hidden by
several other lines of Os II, III, and Pd III, detecting the two lines of Te
III in future events can provide solid confirmation of the Te III production
in mergers. Furthermore, the ratio of these line fluxes can be used to
diagnose the electron temperature.

##  3 Conclusion and discussion

The observed spectra of the kilonova AT 2017gfo exhibit a strong emission line
at $2.1\,{\rm\mu m}$. The emission line with the lack of an apparent blue-
shifted absorption feature suggests that the emission feature is a forbidden
line excited through electron collision. We attribute this line to the fine
structure line, [Te III] $2.10\,{\rm\mu m}$, which has also been detected in
planetary nebulae (Madonna et al., 2018). Note that Te is one of the most
abundant elements in the second r-process peak. We estimate that the mass of
Te III is roughly $10^{-3}M_{\odot}$ to account for the observed line flux.

We compare the observed spectra with a synthetic model, where the spectrum is
composed of fine structure emission lines and a continuum component
approximated by blackbody radiation. The spectrum of fine structure emission
lines is computed with the one-zone model presented in Hotokezaka et al.
(2022). With the solar r-process abundance beyond the first r-process peak,
$T_{e}\sim 2000\,{\rm K}$, and $Y^{+2}\sim 0.3$, we show that [Te III]
$2.10\,{\rm\mu m}$ is the strongest emission line among M1 transitions of all
the heavy elements around 10 days after merger. Our model agrees with the
observed spectra for the ejecta of $0.05M_{\odot}$ with the solar r-process
abundance pattern with $A\geq 88$, an expansion velocity of $0.07c$, and
electron temperature of $\sim 2000\,{\rm K}$, and the ionization fraction of
$Y^{+2}\sim 0.3$. Because blackbody radiation may be a poor approximation to
the continuum flux around $2\,{\rm\mu m}$ we should keep in mind that the
amount of Te III in our analysis may be affected by the continuum flux model.
It is also interesting to note that the same abundance pattern can reproduce
the Spitzer $4.5\,{\rm\mu m}$ detection at 40 days (Villar et al., 2018;
Kasliwal et al., 2022), in which the emission is attributed to a fine
structure line of W III (Hotokezaka et al., 2022). Note that, if the lighter
r-process elements are abundant, they are expected to produce emission lines
around $2\,{\rm\mu m}$ such as [Kr III] $2.20\,{\rm\mu m}$ and [Se IV]
$2.29\,{\rm\mu m}$. However, the observed spectra peak at $2.1\,{\rm\mu m}$,
suggesting that these ions are less abundant compared to Te III in the line
emitting region of the ejecta.

Our model does not include electric dipole (E1) lines, which may produce
strong absorption and emission lines. Recently, Gillanders et al. (2023)
suggested that $2.1\,{\rm\mu m}$ may be composed of two lines and an E1 line
of Ce III is the best candidate producing this line feature. Although we
cannot quantify the contamination of E1 lines to the $2.1\,{\rm\mu m}$
feature, we emphasize that the M1 emission of Te can account the observed line
flux with reasonable parameters. To verify this hypothesis we need spectral
modelings with E1 lines. We also note that the M1 lines in our list cannot
account for the observed feature at $1.6\,{\rm\mu m}$. The flux in this line
declines with time as the continuum flux declines, indicating that this
emission feature may be produced by E1 lines. Interestingly, Domoto et al.
(2022) show that Ce III has several strong E1 lines around $1.6\,{\rm\mu m}$.

From the early blue emission in the photospheric phase, it is suggested that
the emission is dominated by the ejecta composed of light r-process elements
in order to avoid significant absorption in the optical band by lanthanides.
Furthermore, the analyses of the kilonova spectra in the photospheric phase
lead to the similar conclusion. The absorption feature around $0.8\,{\rm\mu
m}$ is likely caused by one of light r-process elements, Sr ($Z=38$), or even
He (Watson et al., 2019; Gillanders et al., 2022; Tarumi et al., 2023). Domoto
et al. (2022) propose that La ($Z=57$) and Ce ($Z=58$) produce the absorption
lines around $1.2$ and $1.5\,{\rm\mu m}$, respectively. But the abundances of
La and Ce inferred from the spectral analysis are lower than the solar
r-process residuals by factor of $\sim 10$. These indicate that the outer part
of the ejecta ($v\gtrsim 0.2c$) is predominantly composed of light r-process
elements. In contrast to the early emission, our analysis implies that heavier
elements, i.e., the second r-process peak, are likely more abundant in the
slower part of the ejecta.

In order to obtain better constraints on the elemental abundances and ejecta
parameters, the spectral modelings should be improved by developing non-LTE
radiation transfer modelings (e.g., Pognan et al., 2022b) and by improving
atomic data such as the radiative transition rates (e.g., Gaigalas et al.,
2019), collision strengths, and recombination rate coefficients. For future
kilonova events, the spectroscopic observations with the JWST as well as
ground-based telescopes will be useful to identify more elements in the
nebular spectra with a wider wavelength range.

## Acknowledgments

We thank Nanae Domoto and Yuta Tarumi for useful discussion. This research was
supported by JST FOREST Program (Grant Number JPMJFR212Y, JPMJFR2136), NIFS
Collaborative Research Program (NIFS22KIIF005), the JSPS Grant-in-Aid for
Scientific Research (19H00694, 20H00158, 21H04997, 20K14513, 20H05639,
22JJ22810).

## DATA AVAILABILITY

The data presented this article will be shared on request to the corresponding
author.

## References

  * Andreoni et al. (2017) Andreoni I., et al., 2017, Publ. Astron. Soc. Australia, 34, e069
  * Bar-Shalom et al. (2001) Bar-Shalom A., Klapisch M., Oreg J., 2001, J. Quant. Spectrosc. Radiative Transfer, 71, 169
  * Barnes & Kasen (2013) Barnes J., Kasen D., 2013, ApJ, 775, 18
  * Barnes et al. (2021) Barnes J., Zhu Y. L., Lund K. A., Sprouse T. M., Vassh N., McLaughlin G. C., Mumpower M. R., Surman R., 2021, ApJ, 918, 44
  * Burbidge et al. (1957) Burbidge E. M., Burbidge G. R., Fowler W. A., Hoyle F., 1957, Reviews of Modern Physics, 29, 547
  * Cameron (1957) Cameron A. G. W., 1957, PASP, 69, 201
  * Chornock et al. (2017) Chornock R., et al., 2017, ApJ, 848, L19
  * DerKacy et al. (2023) DerKacy J. M., et al., 2023, ApJ, 945, L2
  * Domoto et al. (2021) Domoto N., Tanaka M., Wanajo S., Kawaguchi K., 2021, ApJ, 913, 26
  * Domoto et al. (2022) Domoto N., Tanaka M., Kato D., Kawaguchi K., Hotokezaka K., Wanajo S., 2022, arXiv e-prints, p. arXiv:2206.04232
  * Fontes et al. (2020) Fontes C. J., Fryer C. L., Hungerford A. L., Wollaeger R. T., Korobkin O., 2020, MNRAS, 493, 4143
  * Gaigalas et al. (2019) Gaigalas G., Kato D., Rynkun P., Radžiūtė L., Tanaka M., 2019, ApJS, 240, 29
  * Gillanders et al. (2022) Gillanders J. H., Smartt S. J., Sim S. A., Bauswein A., Goriely S., 2022, MNRAS, 515, 631
  * Gillanders et al. (2023) Gillanders J. H., Sim S. A., Smartt S. J., Goriely S., Bauswein A., 2023, arXiv e-prints, p. arXiv:2306.15055
  * Hotokezaka & Nakar (2020) Hotokezaka K., Nakar E., 2020, ApJ, 891, 152
  * Hotokezaka et al. (2021) Hotokezaka K., Tanaka M., Kato D., Gaigalas G., 2021, MNRAS, 506, 5863
  * Hotokezaka et al. (2022) Hotokezaka K., Tanaka M., Kato D., Gaigalas G., 2022, MNRAS, 515, L89
  * Kasen & Barnes (2019) Kasen D., Barnes J., 2019, ApJ, 876, 128
  * Kasen et al. (2013) Kasen D., Badnell N. R., Barnes J., 2013, ApJ, 774, 25
  * Kasliwal et al. (2017) Kasliwal M. M., et al., 2017, Science, 358, 1559
  * Kasliwal et al. (2022) Kasliwal M. M., et al., 2022, MNRAS, 510, L7
  * Kawaguchi et al. (2018) Kawaguchi K., Shibata M., Tanaka M., 2018, ApJ, 865, L21
  * Kramida et al. (2021) Kramida A., Yu. Ralchenko Reader J., and NIST ASD Team 2021, NIST Atomic Spectra Database (ver. 5.9), [Online]. Available: https://physics.nist.gov/asd [2017, April 9]. National Institute of Standards and Technology, Gaithersburg, MD.
  * Kwok et al. (2023) Kwok L. A., et al., 2023, ApJ, 944, L3
  * Lattimer & Schramm (1974) Lattimer J. M., Schramm D. N., 1974, ApJ, 192, L145
  * Madonna et al. (2018) Madonna S., et al., 2018, ApJ, 861, L8
  * Margutti & Chornock (2021) Margutti R., Chornock R., 2021, ARA&A, 59
  * Metzger (2017) Metzger B. D., 2017, Living Reviews in Relativity, 20, 3
  * Nakar (2020) Nakar E., 2020, Phys. Rep., 886, 1
  * Perego et al. (2022) Perego A., et al., 2022, ApJ, 925, 22
  * Pian et al. (2017) Pian E., et al., 2017, Nature, 551, 67
  * Pognan et al. (2022a) Pognan Q., Jerkstrand A., Grumer J., 2022a, MNRAS, 510, 3806
  * Pognan et al. (2022b) Pognan Q., Jerkstrand A., Grumer J., 2022b, MNRAS, 513, 5174
  * Smartt et al. (2017) Smartt S. J., et al., 2017, Nature, 551, 75
  * Sneppen et al. (2023) Sneppen A., Watson D., Bauswein A., Just O., Kotak R., Nakar E., Poznanski D., Sim S., 2023, Nature, 614, 436
  * Tanaka & Hotokezaka (2013) Tanaka M., Hotokezaka K., 2013, ApJ, 775, 113
  * Tanaka et al. (2020) Tanaka M., Kato D., Gaigalas G., Kawaguchi K., 2020, MNRAS, 496, 1369
  * Tanvir et al. (2017) Tanvir N. R., et al., 2017, Astrophys. J., 848, L27
  * Tarumi et al. (2023) Tarumi Y., Hotokezaka K., Domoto N., Tanaka M., 2023, arXiv e-prints, p. arXiv:2302.13061
  * Troja et al. (2017) Troja E., et al., 2017, Nature, 551, 71
  * Villar et al. (2018) Villar V. A., et al., 2018, ApJ, 862, L11
  * Watson et al. (2019) Watson D., et al., 2019, Nature, 574, 497
  * Waxman et al. (2019) Waxman E., Ofek E. O., Kushnir D., 2019, ApJ, 878, 93
  * Wooden et al. (1993) Wooden D. H., Rank D. M., Bregman J. D., Witteborn F. C., Tielens A. G. G. M., Cohen M., Pinto P. A., Axelrod T. S., 1993, ApJS, 88, 477


#  Power law decay of entanglement quantifiers in a single agent to a many
body system coupling

Ohad Shpielberg<EMAIL_ADDRESS>Haifa Research Center for Theoretical
Physics and Astrophysics, University of Haifa, Mt. Carmel, Haifa 31905, Israel

###### Abstract

Manipulating many body quantum systems is a challenge. A useful way to achieve
it would be to entangle the system to a diluted system, with a small particle
number. Preparation of such entangled states can be facilitated as ground
state of a many body Hamiltonian or the steady state of a many body open
quantum system. Here we study two-site lattice models with a conserved boson
number, biased to display a large occupancy in one of the sites. The Von
Neumann entanglement entropy as well as the Logarithmic negativity show a
typical power law decay in $R$, the occupancy ratio between the two sites.
These results imply that it is feasible to entangle a large many body system
to a single atom, as recently reported experimentally.

##  I Introduction

Entanglement is a key resource in quantum information Horodecki  _et al._
(2009); Amico  _et al._ (2008); Vedral (2018), quantum computing Nielsen and
Chuang (2002); Preskill (2018) and quantum metrology Pezze  _et al._ (2018).
Recently, there has been significant advancement in generating, manipulating
and measuring entangled many body states; both experimentally and
theoretically Omran  _et al._ (2019); Islam  _et al._ (2015); Kaufman  _et
al._ (2016); Kallush  _et al._ (2021). Both preparation of the entangled state
and its validation using, e.g. entanglement witnesses Chruściński and Sarbicki
(2014), are challenging aspects in many body systems and are the focus of
ongoing research Lacroix (2020); Weimer  _et al._ (2010); Carr and Saffman
(2013); Verstraete  _et al._ (2009).

Preparation of a desired entangled state can be realized as the ground state
of a carefully designed Hamiltonian. Therefore, understanding the entanglement
properties of ground states is of practical importance. Significant effort has
been directed for extended systems, especially in $1D$. The ground state of a
Hamiltonian with local interactions typically exhibits an area law in the
bipartite Von Neumann entanglement entropy Eisert  _et al._ (2008), contrary
to the generic volume law of typical quantum states. However, the area law
does not characterise a system composed of a few sites, where each site can
occupy a large number of particles (see Fig. 1). For two sites, the average
Von Neumann entanglement entropy is known Page (1993); Sen (1996), but the
characteristic properties of the ground state entanglement remain largely
unexplored.

A particularly appealing case is for a single agent in system $B$ to be
entangled to a large number of particles in system $A$. Such entanglement
allows to manipulate the many body system via the single agent. This setup was
experimentally demonstrated in McConnell, R., Zhang, H., Hu, J. et al. (2015),
where a single photon was entangled with roughly $3,000$ atoms. At this point,
it is unknown whether there is a limit to the number of particles that could
realistically be entangled with a single agent. It is further unknown whether
our setup leads to typical ground state entanglement properties.

To answer the above questions, the entanglement needs to be quantified Friis
_et al._ (2019); Plenio and Virmani (2007). Choosing an appropriate
entanglement quantifier depends on the intended application, e.g. entanglement
distillation to produce Bell states. However, the entanglement quantifier can
alternatively be chosen to accommodate fast calculations of known density
matrices, e.g. the Logarithmic negativity $E_{ln}$. For pure states, the Von
Neumann entanglement entropy $E_{vn}$ serves both purposes.

Figure 1: The systems $A$ and $B$ are made to interact such that we find a
large occupancy of particles at system $A$ and a low occupancy at system $B$.

Keeping in mind the motivation of entangling a single agent to a many body
system, we turn to a simpler theoretical setup.  In this work we consider $N$
bosons, occupying a two-site system. We study the ground state entanglement
properties when the system is tuned to display an overwhelming majority of
bosons occupying site $A$ (see Fig. 1). To make this statement precise, let
$\hat{n}_{A,B}$ be the corresponding number operators. Then,
$R=\langle\hat{n}_{A}\rangle/\langle\hat{n}_{B}\rangle$ is the ratio of the
particle occupancies. We study the large $R$ behaviour of both the ground
state Von Neumann entanglement entropy and Logarithmic negativity of the Bose-
Hubbard Hamiltonian. For open quantum systems, the same setup can be
considered, where the steady state takes the role of the ground state. We
study the steady state Logarithmic negativity of a Lindblad super-operator
model – the quantum asymmetric inclusion process at large $R$ values Eisler
(2011); Bernard  _et al._ (2018).

Both models are studied at different scaling regimes. Nevertheless, they all
consistently lead to a power law decay in the entanglement quantifiers.
Quantitatively, the Von Neumann entanglement entropy $E_{vn}\sim\frac{\log
R}{R^{\alpha}}$ and the Logarithmic negativity
$E_{ln}\sim\frac{1}{R^{\alpha}}$ for $R\gg 1$. See Tables 1 and 2 for a
summary of the results.

We argue that the power law decay is typical, as we have considered two
disparate models and different scaling regimes for each model. The slow power
law decay, contrasting with an exponential decay, answers in a quantifiable
way how realistic it is to entangle a single, or a few atoms to a highly
occupied many body state. The exponent $\alpha$ is non-universal. Therefore,
interacting systems that result in small $\alpha$ values are favorable to
facilitate entanglement between the diluted system to the large occupancy
system.

The structure of this paper is as follows. In Sec. II we present the
Hamiltonian and Lindblad models and summarise the main results. Sec. III
presents in full the analytical and numerical treatment of the systems under
study. Finally, Sec. IV recaps the main findings and their physical relevance
and suggests future directions.

##  II Models and results

The aim of this work is to quantify the bipartite entanglement of a composite
$AB$ system at large $R$ values. Therefore, it is natural to study lattice
models, where the distinction between the two subsystems is clear cut. In
particular, we study lattice models with two sites, $A$ and $B$.

To demonstrate the power law behaviour is typical, two disparate lattice
models are considered. First, the ground state entanglement of the two-site
Bose-Hubbard model is extensively studied. Second, we consider a
generalization of the asymmetric inclusion process Grosskinsky  _et al._
(2011) to the quantum realm via a Lindblad equation, dubbed here the quantum
asymmetric inclusion process (QASIP). We then study the entanglement
properties of the steady state at large $R$.

###  II.1 The two-site Bose-Hubbard model

The Bose-Hubbard model is a simple yet rich many body lattice model of spin-
less bosons. It allows studying the superfluid-insulator transition Fisher
_et al._ (1989) and can be experimentally implemented using optical lattices
Jaksch  _et al._ (1998); Greiner  _et al._ (2002). The particular case of the
two-site Bose-Hubbard model was extensively used in the literature to study
tunneling effects between potential wells Links  _et al._ (2006) as well as
fragmentation Spekkens and Sipe (1999). Importantly for our purposes, the two
site Bose-Hubbard model is expected to be both analytically tractable and to
present typical physical behaviour in terms of the ground state entanglement.
Hence, it serves as the starting point of our analysis.

The two-site Bose-Hubbard Hamiltonian is given by

$H_{BH}=-J(\hat{b}^{\dagger}_{A}\hat{b}_{B}+\hat{b}^{\dagger}_{B}\hat{b}_{A})-\mu\hat{n}_{A}+\frac{U}{2}\sum_{i=A,B}\hat{n}_{i}-\hat{n}^{2}_{i},$
(1)

where $J$ is the hopping matrix element between neighboring sites and $U$
determines the strength of the on-site interaction. The operators
$\hat{b}_{i},\hat{b}^{\dagger}_{i}$ are the site-dependent bosonic creation
and annihilation operators and $\hat{n}_{i}=\hat{b}^{\dagger}_{i}\hat{b}_{i}$
is the number operator for $i=A,B$. In an optical lattice, the potential wells
are represented by the two sites Jaksch  _et al._ (1998). The potential offset
between the two asymmetric potential wells is given by $\mu$. It furthermore
allows to imbalance the system towards large $R$ values.

It is useful to note that the total particle number,
$\hat{N}=\hat{n}_{A}+\hat{n}_{B}$ is conserved. Therefore, we analyze the
ground state with $N$ bosons. In what follows, we consider two scaling schemes
leading to large $R$ values.

First, taking large $\mu$ values and keeping $N,J$ and $U$ fixed, a
perturbative treatment leads to $R=\mu^{2}/J^{2}+O(\mu)$. See Sec. III. In
this limit, and as long as $\frac{J\sqrt{N}}{\mu},\frac{UN^{2}}{\mu}\ll 1$, we
find analytically the power law behaviour described in Table 1. These results
are also numerically corroborated in Fig. 2. Note that the logarithmic
correction in the Von Neumann entanglement hardly changes the behaviour from a
clean power law.

Figure 2:  Corroborating the perturbation theory analytical predictions
presented in Table 1. The entanglement quantifiers are evaluated numerically
for the two-site Bose-Hubbard model and compared with perturbation theory.
Different $U,J,N$ values are considered (see legend), in the range
$\mu\in\left[10,10^{3}\right]$ to facilitate large $R$ values. No fitting
parameters are required to observe the collapse onto the expected power law
behaviour.

Second, we consider the large $N$ limit, with fixed $\mu,J$ and $U$. A
perturbative approach is harder in this case as the Hilbert space of the
effective Hamiltonian depends on the particle number $N$ (see Sec. III). See
Fig. 3, Table 1 and the appendix C for the numerical analysis of the large $N$
limit.

Figure 3:  Numerical fitting of the entanglement quantifiers at the large $N$ limit of the two-site Bose-Hubbard model. Both $\frac{E_{vn}}{\log R}$ and $E_{ln}$ exhibit power law behaviour at large $R$ values, albeit with different exponents. Numerical convergence of arbitrary parameters to the same exponent is obtained when sufficiently large $R$ values are reached. (a) The Von Neumann entanglement entropy is shown to scale like $E_{vn}\propto\frac{\ln R}{R^{\alpha}}$, where $\alpha\approx 0.518$. The range $N\in\left[10,400\right]$ was tested. (b) The logarithmic negativity is shown to scale like $E_{ln}\propto R^{-\alpha}$, where $\alpha\approx 0.22$ is numerically fitted in the range $N\in\left[10,40\right]$. See the Appendix C for a detailed discussion on the reduced range of $N$ for the logarithmic negativity.  | $R$ | $E_{vn}=\Gamma\frac{\ln R}{R^{\alpha}}$ | $E_{ln}=\frac{\Gamma}{R^{\alpha}}$
---|---|---|---
$\mu\gg 1$ | $\mu^{2}/J^{2}$ | $\Gamma=N,\alpha=1$ | $\Gamma=\frac{2\sqrt{N}}{\ln 2},\alpha=\frac{1}{2}$
$N\gg 1$ | $\propto N^{\beta},\,\beta\approx 2.15$ | $\alpha\approx 0.518$ | $\alpha\approx 0.220$
Table 1: The power law behaviour of the ground state entanglement of the Bose-
Hubbard model at large $R$ values. For the large $N$ limit, the exponents are
evaluated numerically.

In conclusion, the ground state of the two-site Bose-Hubbard model leads to a
power law behaviour of both the Von Neumann entanglement entropy and the
Logarithmic negativity at large $R$ values in two different scaling schemes.

Next, we explore the large $R$ steady state entanglement in an open quantum
system setup: the quantum asymmetric inclusion process.

###  II.2 Open quantum system

Realistically, physical systems are never truly isolated. Interestingly, the
coupling of a quantum system to an environment does not always lead to
complete loss of entanglement in the system. Indeed, there are examples where
the environment can be engineered to produce a desirable entangled state
Verstraete  _et al._ (2009); Shpielberg (2020).

Here, we aim to study whether the large $R$ power law behaviour persists for
steady states of open quantum systems as well. To that end, we focus on a
quantum analog of the asymmetric inclusion process given in terms of the
Lindblad equation.

In this setup, one usually assumes that a quantum system is coupled to an
environment with fast relaxation times. This in turn, allows to discard non-
Markovian contributions to the evolution of the density matrix and results in
the Lindblad equation Petruccione, F. and Breuer, H.-P. (2002); Gorini  _et
al._ (1976); Lindblad (1976)

$\displaystyle\partial_{t}\rho$ $\displaystyle=$
$\displaystyle\mathcal{H}(\rho)+\sum_{k}\mathcal{D}_{\hat{L}_{k}}(\rho)$ (2)
$\displaystyle\mathcal{H}(\rho)$ $\displaystyle=$
$\displaystyle-i\left[H,\rho\right]$
$\displaystyle\mathcal{D}_{\hat{L}}(\rho)$ $\displaystyle=$
$\displaystyle\hat{L}\rho\hat{L}^{\dagger}-\frac{1}{2}\\{\hat{L}^{\dagger}\hat{L},\rho\\}.$

In eq.(2), $H$ is a Hermitian operator and
$\left[\bullet,\bullet\right],\\{\bullet,\bullet\\}$ are the commutation and
anti-commutation relations correspondingly. Despite the restriction to
Markovian dynamics, enough quantumness remains in the Lindblad equation Spohn,
Herbert and Lebowitz, Joel L. (1978); Kosloff, Ronnie (2013); Dutta and Cooper
(2019). At this point, we turn to study a particular Lindblad model, allowing
to facilitate the large $R$ limit analytically.

The QASIP describes the dynamics of bosons on a two-site lattice, where the
boson interactions are environment assisted Eisler (2011); Bernard  _et al._
(2018). The evolution of the density matrix is given by

$\displaystyle\partial_{t}\hat{\rho}$
$\displaystyle=\mathcal{L}_{QASIP}(\hat{\rho})$ (3)
$\displaystyle\mathcal{L}_{QASIP}(\hat{\rho})$
$\displaystyle=\mathcal{H}_{\mathrm{tb}}(\hat{\rho})+\mathcal{L}_{\mathrm{D}}(\hat{\rho})+\mathcal{L}_{E}(\hat{\rho})$
$\displaystyle H_{\mathrm{tb}}(\hat{\rho})$
$\displaystyle=\varepsilon\left(\hat{b}^{\dagger}_{A}\hat{b}_{B}+\hat{b}^{\dagger}_{B}\hat{b}_{A}\right)$
$\displaystyle\mathcal{L}_{\mathrm{D}}(\hat{\rho})$
$\displaystyle=\eta\sum_{k=A,B}\mathcal{D}_{\hat{n}_{k}}(\hat{\rho})$
$\displaystyle\mathcal{L}_{\mathrm{E}}(\hat{\rho})$
$\displaystyle=\gamma\mathcal{D}_{\hat{b}^{\dagger}_{A}\hat{b}_{B}}(\hat{\rho}).$

See eq.(2) for the definition of $\mathcal{H},\mathcal{D}_{\hat{L}}$ for any
operator $\hat{L}$. Here, the Hermitian $H_{\rm{tb}}$ is a tight-binding
Hamiltonian inducing particle jumps between the two sites.
$\mathcal{L}_{\mathrm{D}}$ is responsible for dephasing at each site and
$\mathcal{L}_{\mathrm{E}}$ explicitly breaks the symmetry between the two
sites and induces the occupation imbalance. Later on, it will be shown that
controlling $\gamma$ allows to induce large $R$ values.

The relation of eq.(3) to the classical asymmetric inclusion process (ASIP)
Grosskinsky  _et al._ (2011) is as follows. The dephasing and bias terms alone
acting on the diagonal terms of the density matrix in the number basis lead to
the ASIP master equation. Namly, the quantum master equation is split to the
diagonal terms and the coherent terms, each having a closed set of equations
(in the number operator basis). The tight-binding Hamiltonian mixes between
the cohrent terms and the diagonal terms, hence a quantum ASIP.

Three comments are in order before we present the results. First, note that
here one cannot assume a priory that the quantum system is coupled to a series
of thermalized baths as we have not performed a microscopic derivation of the
Lindblad equation. Since we are not interested in studying thermalization, the
lack of a microscopic derivation is of no importance. Second, the steady state
density matrix is not pure. Hence, we will only use the Logarithmic negativity
as an entanglement quantifier. Third, in eq.(3) we have set $\hbar$ to unity.
Furthermore, we will assume $\varepsilon,\eta,\gamma$ and the time $t$ to be
dimensionless for convenience. When presenting the different scaling schemes,
the inverse time dimensions of $\varepsilon,\eta,\gamma$ could be restored.

The QASIP, like the Bose-Hubbard model can be shown to conserve the particle
number $\hat{N}=\hat{n}_{A}+\hat{n}_{B}$ (see the appendix B). However, a
related but more general property exists for the QASIP. In the number operator
basis, we can write the density matrix as

$\displaystyle\hat{\rho}$ $\displaystyle=$
$\displaystyle\sum_{S}a_{S}\hat{\rho}_{S},\quad\textrm{where}$ (4)
$\displaystyle\hat{\rho}_{S}$ $\displaystyle=$
$\displaystyle\sum^{S}_{x,y=0}\varrho_{S}(x,y)\ket{x,S-x}\bra{y,S-y}.$

Here, $S$ takes non-negative integer values and $a_{S}$ are non-negative
prefactors that sum to $1$. Note that the hermitianity of the density matrix
$\hat{\rho}_{S}$ implies $\varrho^{*}_{S}(x,y)=\varrho_{S}(y,x)$ and unity
trace implies $\sum^{S}_{x=0}\varrho_{S}(x,x)=1$.

The dynamics of eq.(3) is restricted to the subspace of $\hat{\rho}_{S}$:

$\displaystyle\partial_{t}\varrho_{S}(x,y)$ $\displaystyle=$
$\displaystyle-\eta(x-y)^{2}\varrho_{S}(x,y)$
$\displaystyle+\gamma\,x_{-}y_{-}\varrho_{S}(x-1,y-1)$
$\displaystyle-\frac{1}{2}\gamma\left(x^{2}_{+}+y^{2}_{+}\right)\varrho_{S}(x,y)$
$\displaystyle-i\varepsilon\sum_{z=\pm
1}x_{z}\varrho_{S}(x+z,y)-y_{z}\varrho_{S}(x,y+z)$ $\displaystyle X_{+}$
$\displaystyle=$ $\displaystyle\sqrt{(X+1)(S-X)}$ $\displaystyle X_{-}$
$\displaystyle=$ $\displaystyle\sqrt{X(S-X+1)},\quad\textrm{where }X=x,y.$

Namely, we have replaced the treatment of the infinite dimensional density
matrix $\hat{\rho}$, with a treatment of finite dimensional $(S+1)^{2}$
density matrices $\hat{\rho}_{S}$ at fixed $S$.

The conserved number $S$ of $\hat{\rho}_{S}$ equals the number of particles in
the system as

$\Tr\hat{\rho}_{S}\hat{N}=\sqrt{\Tr\hat{\rho}_{S}\hat{N}^{2}}=S.$ (6)

Eq.(II.2) implies therefore that the process conserves the particle number.
From hereon out, we may replace $S$ by $N$.

In eq.(3), there are different scaling schemes leading to the large $R$ limit
at the steady state. The Table 2 summarises the power law behaviour of the
Logarithmic negativity in three different scaling regimes. The results in
Table 2 were verified both analytically and numerically. See Fig. 4, for the
large $\gamma$ limit and Fig. 5 for the large $\eta$ limit. For the large $N$
only numerical evidence is currently present, see Fig. 6 and 8.

Figure 4:  Corroborating the perturbation theory analytical predictions presented in Table 2. The Logarithmic negativity is evaluated numerically for the QASIP model and compared with perturbation theory at the large $\gamma$ limit in the range $\gamma\in\left[1,200\right]$. Different $\varepsilon,\eta,N$ values are considered, all showing the expected collapse onto the $1/\sqrt{R}$ plot at large $R$ values, with no fitting parameters.  Figure 5:  Corroborating the perturbation theory analytical predictions presented in Table 2. The Logarithmic negativity is evaluated numerically for the QASIP model and compared with perturbation theory at the large $\eta$ limit in the range $\eta\in\left[1,500\right]$. Different $\varepsilon,\gamma,N$ values are considered, all showing the expected collapse onto the $1/R$ plot at large $R$ values, with no fitting parameters.  | $R$ | $E_{ln}(R)=\Gamma R^{-\alpha}$
---|---|---
$\gamma\gg 1$ | $\frac{\gamma^{2}N^{2}}{4\varepsilon^{2}}$ | $\Gamma=\frac{2\sqrt{N}}{\ln 2},\quad\alpha=1/2$
$\eta\gg 1$ | $\eta\frac{N\gamma}{2\varepsilon^{2}}$ | $\Gamma=\frac{2N^{3/2}\gamma}{\varepsilon\ln 2},\quad\alpha=1$
$N\gg 1$ | $\propto N^{\beta},\beta\approx 1.988$ | $\alpha\approx 0.236$
Table 2: The power law behaviour of the Logarithmic negativity in the steady
state QASIP. For the large $N$ limit, the exponents are evaluated numerically.

In conclusion, the steady state QASIP in eq.(3) exhibits a power law decay in
the Logarithmic negativity, similarly to the Bose-Hubbard dynamics.

In the next section, we provide a detailed derivation of the results for the
Bose-Hubbard model and for the QASIP.

Figure 6:  Numerical fitting of the Logarithmic negativity at the large $N$
limit of the QASIP model. Different $\eta,\gamma$ values were tested, leading
to the large $R$ limit at $\varepsilon=1,N\in[1,25]$. Numerically, it is shown
that $E_{ln}\propto R^{-\alpha}$, where $\alpha\approx 0.236$.

##  III Analytical and numerical analysis

In Sec. II, we have introduced two lattice models: the Bose-Hubbard model and
the quantum asymmetric inclusion process. The power law behaviour of the Von
Neumann entanglement and the Logarithmic negativity was summarized in Tables 1
and 2. In this section, we describe the analysis of these results in detail.

###  III.1 Bose-Hubbard model

To analyze the entanglement properties of the two-site Bose-Hubbard model in
eq.(1), we need to find the $N$-particle ground state of the Hamiltonian.
Given a description of the ground state, finding the ratio $R$ and the
entanglement quantifiers $E_{vn},E_{ln}$ becomes straight-forward, but
sometimes technically cumbersome. See the appendix A.

The Hilbert space of $N$ particles for the Bose-Hubbard Hamiltonian eq.(1) is
spanned by the $N+1$ Fock states
$\ket{n_{A},n_{B}}=\frac{1}{n_{A}!n_{B}!}(\hat{b}^{\dagger}_{A})^{n_{A}}(\hat{b}^{\dagger}_{B})^{n_{B}}\ket{0,0}$.
Namely,

$\ket{\psi_{N}}=\sum^{N}_{k=0}a_{k}\ket{k,N-k},\quad\textrm{where
}\sum^{N}_{k=0}|a_{k}|^{2}=1.$ (7)

Then, the $N$ particle Bose-Hubbard Hamiltonian can be written as a $N+1\times
N+1$ matrix. Finding the ground state can be done analytically and for all
$\mu,U,J$ values when we set $N=1$. This simple case will reveal the intuitive
scaling limits leading to the large $R$ behaviour, both for the Bose-Hubbard
model as well as for the quantum asymmetric inclusion process in the next
subsection.

Indeed, for $N=1$, the Bose-Hubbard Hamiltonian can be represented by the
$2\times 2$ matrix

$H^{(N=1)}_{BH}=\begin{pmatrix}-\mu&-J\\\ -J&0,\end{pmatrix}$ (8)

where the wavefunction is in its most general form

$\ket{\psi_{N=1}}=\cos\zeta\ket{10}+{\rm
e}^{i\phi_{1}}\sin\zeta\ket{01}=\begin{pmatrix}\cos\zeta\\\ {\rm
e}^{i\phi_{1}}\sin\zeta\end{pmatrix}$ (9)

for real $\zeta,\phi_{1}$ values and the right hand side of eq.(9) is in a
vector notation corresponding to the matrix in eq.(8). Clearly, the ratio is
$R=\cot^{2}\zeta$. The lowest eigenvalue of eq.(8) is
$\epsilon=-\frac{1}{2}(\mu+\sqrt{\mu^{2}+4J^{2}})$ with the ground state
$\ket{\psi_{N=1}}=\frac{1}{\mathcal{N}}(-\epsilon\ket{10}+J\ket{01})$ and
$\mathcal{N}^{2}=\epsilon^{2}+J^{2}$ is a normalization constant. This implies
that $R=\epsilon^{2}/J^{2}$. Therefore, in this particular case, $R\gg 1$ only
if $\mu/J\gg 1$ and leads to $R=(\mu/J)^{2}+O(\mu)$ asymptotically.

The Von Neumann entanglement for the ground state in eq.(9) is

$E_{vn}=-2\cos^{2}\zeta\ln(\cos\zeta)-2\sin^{2}\zeta\ln(\sin\zeta).$ (10)

Using trigonometric identities, we recover $\cos^{2}\zeta=\frac{R}{1+R}$ and
$\sin^{2}\zeta=\frac{1}{1+R}$. Asymptotically for large $R$, $E_{vn}=\frac{\ln
R}{R}$ to leading order as reported in Table 1.

To find the Logarithmic negativity, we need to write the partially transposed
density matrix

$\displaystyle\rho^{PT}$ $\displaystyle=$
$\displaystyle\cos^{2}\zeta\ket{10}\bra{10}+\sin^{2}\zeta\ket{01}\bra{01}$
$\displaystyle+\cos\zeta\sin\zeta\left({\rm
e}^{i\phi_{1}}\ket{00}\bra{11}+{\rm e}^{-i\phi_{1}}\ket{11}\bra{00}\right).$

The eigenvalues of the partially transposed density matrix are
$\cos^{2}\zeta,\sin^{2}\zeta$ and $\pm\frac{\sqrt{1-\cos 4\zeta}}{2\sqrt{2}}$.
Therefore, the Logarithmic negativity is $E_{ln}=\log_{2}(1+\frac{\sqrt{1-\cos
4\zeta}}{\sqrt{2}})$. For large $R$ and to leading order, we find
$E_{ln}=\frac{2}{\ln 2}R^{-1/2}+O(R^{-1})$ as reported in Table 1.

For $N>1$, the ground state solution becomes cumbersome, but still requires
dealing with a $N+1\times N+1$ Bose-Hubbard matrix. The numerical code that
produced Fig. 2,3 and 7 finds the ground state of the Bose-Hubbard matrix at
some finite $N$. Then, it calculates the Von Neumann entanglement and the
Logarithmic negativity.

From the $N=1$ example , we have seen that the large $\mu$ limit leads to a
large bias $R\gg 1$. This happens when the $\mu$ term dominates the energy of
the ground state, i.e. for $\mu\gg\sqrt{N}J,UN^{2}$. Another large $R$ limit
is recovered for finite $\mu,U,J>0$ and large $N$. In this limit, the
particles condense due to the strong attractive energy $\sim UN^{2}$. The
symmetry is broken by the potential offset $\mu$, leading to condensation of
the particles in site $A$ and to large $R$ values.

A perturbative approach for the large $N$ limit is non-trivial due to as the
dimension of the effective Hilbert space changes with $N$. However, a direct
numerical analysis clearly reveals the power law behaviour in this limit. See
Table 1, Sec.II and the appendix C for more details. In what follows, we
consider the large $\mu$ limit and evaluate the ground state and the
entanglement quantifiers using a perturbative approach.

Let us develop a standard perturbation theory for the Hamiltonian
$H=H_{0}+\frac{1}{\mu}H_{1}{\color[rgb]{0,0,0}=\frac{1}{\mu}H_{BH}}$

$H_{0}=-\hat{n}_{A}\quad H_{1}=H_{BH}-\mu H_{0},$ (12)

at large $\mu$.

The eigenstates of $H_{0}$ are $\ket{\phi^{(0)}_{n}}=\ket{n,N-n}$ with
energies $\epsilon^{(0)}_{n}=-n$. The first order correction to the ground
state is

$\ket{\phi_{N}}=\ket{N,0}+\lambda\ket{N-1,1}+O(\lambda^{2})$ (13)

where $\lambda=\frac{J\sqrt{N}}{\mu}$ assumed to be small as well as
$UN^{2}/\mu\ll 1$. At this limit we find
$R=\langle\hat{n}_{A}\rangle/\langle\hat{n}_{B}\rangle=N/\lambda^{2}.$ So, we
can approximate at small $\lambda$, $R=\mu^{2}/J^{2}+O(\mu)$. The Von Neumann
entanglement can be calculated for the ground state $\ket{\phi_{N}}$

$\displaystyle
E_{vn}=-\frac{1}{1+\lambda^{2}}\ln\frac{1}{1+\lambda^{2}}-\frac{\lambda^{2}}{1+\lambda^{2}}\ln\frac{\lambda^{2}}{1+\lambda^{2}}$
(14) $\displaystyle=-\lambda^{2}\ln\lambda^{2}+O(\lambda^{3}).$

At this limit, we find the reported scaling

$E_{vn}=\frac{N}{R}\ln\frac{N}{R}$ (15)

to leading order. This is corroborated numerically in Fig. 2 and summarised in
Table 1. Recall that in this scaling, $N$ may be large, but $R\gg N$.

To find the Logarithmic negativity, we need to calculate the eigenvalues of
the partially transposed density matrix of $\ket{\phi_{N}}$. For $N>1$, the
only non-zero eigenvalues are
$\frac{1}{1+\lambda^{2}},\frac{\lambda^{2}}{1+\lambda^{2}},\frac{\pm\lambda}{1+\lambda^{2}}$.
This leads to

$E_{ln}=\log_{2}\left(1+\frac{2\lambda}{1+\lambda^{2}}\right)=\frac{2}{\ln
2}\sqrt{\frac{N}{R}}+O(1/R).$ (16)

For large $R$ values, where the perturbation theory applies, the Logarithmic
negativity dominates the Von Neumann entanglement as it should Plenio and
Virmani (2007); Plenio  _et al._ (2005). Again, we refer to Fig. 2 to see the
excellent agreement with the numerical evaluation.

Other scaling schemes, leading to large $R$ values can exist. Nevertheless,
the power law behaviour of the entanglement quantifiers is believed to
persist, based on the $N=1$ exactly solvable cases.

We turn to study the large $R$ entanglement properties of a completely
different setup – the quantum asymmetric inclusion process.

###  III.2 The QASIP

To analyze the steady state entanglement properties of the QASIP at large $R$,
we need to find the steady state density matrix with a fixed $S$, i.e.
$\mathcal{L}_{QASIP}(\hat{\rho}_{S})=0$. Namely, we wish to find
$\varrho_{S}(x,y)$ such that the right hand side of eq.(II.2) vanishes.

As in the Bose-Hubbard model, it is useful to study first the simple case of
$S=1$. Here, $x,y=\\{0,1\\}$ and demanding a steady state in eq.(II.2) leads
to

$\begin{pmatrix}-\gamma&i\varepsilon&-i\varepsilon&0\\\
i\varepsilon&-\frac{1}{2}\gamma-\eta&0&-i\varepsilon\\\
-i\varepsilon&0&-\frac{1}{2}\gamma-\eta&i\varepsilon\\\
\gamma&-i\varepsilon&i\varepsilon&0\end{pmatrix}\begin{pmatrix}\varrho_{1}(0,0)\\\
\varrho_{1}(0,1)\\\ \varrho_{1}(1,0)\\\
\varrho_{1}(1,1)\end{pmatrix}=\begin{pmatrix}0\\\ 0\\\ 0\\\ 0\end{pmatrix}.$
(17)

Solving eq.(17), we find the steady state solution for $S=1$

$\displaystyle\mathcal{N}_{1}\hat{\rho}_{1}$ $\displaystyle=$ $\displaystyle
4\varepsilon^{2}\ket{0,1}\bra{0,1}$
$\displaystyle+2i\gamma\varepsilon\left(\ket{1,0}\bra{0,1}-\ket{0,1}\bra{1,0}\right)$
$\displaystyle+(4\varepsilon^{2}+\gamma^{2}+2\gamma\eta)\ket{1,0}\bra{1,0}.$

where $\mathcal{N}_{1}=\gamma^{2}+2\gamma\eta+8\varepsilon^{2}$ ensures unity
trace of $\hat{\rho}_{1}$. From the steady state solution we recover
$R=1+\frac{\gamma(\gamma+2\eta)}{4\varepsilon^{2}}$ and
$E_{ln}=\log_{2}\left(1+\frac{4\gamma\varepsilon}{\gamma^{2}+8\varepsilon^{2}+2\gamma\eta}\right)$.
A few observations can already be made. At $\eta\rightarrow\infty$ the
entanglement vanishes as can be expected in the large dephasing limit Bernard
_et al._ (2018); Fischer  _et al._ (2016). Moreover, the Logarithmic
negativity becomes positive due to a combination of biasing and coherent
hopping, i.e. $\gamma\varepsilon>0$. We identify two limits where $R$ becomes
large. For $\gamma\rightarrow\infty$ and finite $\eta,\varepsilon$ we recover
$R\propto\gamma^{2}$ and $E_{ln}\propto 1/\sqrt{R}$. Similarly, for finite
$\gamma,\varepsilon$ and large $\eta$ we recover $R\propto\eta$ and
$E_{ln}\propto 1/R$. Already for the $S=1$ case we find that the entanglement
power law behavior persists for large $R$ values. However, the exponent is
non-universal and depends on the scaling scheme.

Exact solution of the steady state for $S>1$ is at best tedious. Instead, we
will find the steady state solution at the two limits noted above, using a
perturbative approach; The limits at finite $S$ will be shown to agree with
the $S=1$ exact analysis carried out in the above.

####  III.2.1 Large asymmetry between the sites

For $\gamma\gg 1$, and at finite $S,\eta,\varepsilon$, we develop the steady
state density matrix as a perturbative sum

$\hat{\rho}_{S}=\frac{1}{\mathcal{N}_{\gamma}}(\hat{\rho}^{(0)}_{S}+\frac{1}{\gamma}\hat{\rho}^{(1)}_{S}+\frac{1}{\gamma^{2}}\hat{\rho}^{(2)}_{S})+O(1/\gamma^{3}),$
(19)

where $\mathcal{N}_{\gamma}$ is a normalization constant to assure trace one
of the truncated density matrix. This perturbative approach implies the order
by order steady state solutions

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\mathcal{L}_{\mathrm{E}}(\hat{\rho}^{(0)}_{S})$ (20)
$\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{1}{\gamma}\mathcal{L}_{\mathrm{E}}(\hat{\rho}^{(1)}_{S})+(\mathcal{L}_{\mathrm{D}}+\mathcal{H}_{tb})\hat{\rho}^{(0)}_{S}$
(21) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\frac{1}{\gamma}\mathcal{L}_{\mathrm{E}}(\hat{\rho}^{(2)}_{S})+(\mathcal{L}_{\mathrm{D}}+\mathcal{H}_{tb})\hat{\rho}^{(1)}_{S}$
(22)

Eq.(20) admits a uniqe solution $\hat{\rho}^{(0)}_{S}=\ket{S,0}\bra{S,0}$,
namely to leading order site $A$ is maximally occupied and site $B$ is
depleted. Using the leading order solution, we find
$\hat{\rho}^{(1)}_{S}=\frac{2i\varepsilon}{\sqrt{S}}(\ket{S,0}\bra{S-1,1}+\ket{S-1,1}\bra{S,0})$.
Note that to first order, there are yet no corrections to the occupancies.
Hence, we solve to second order in $1/\gamma$ obtaining

$\displaystyle\hat{\rho}^{(2)}_{S}$ $\displaystyle=$
$\displaystyle\frac{4\varepsilon^{2}}{\gamma^{2}S}\ket{S-1,1}\bra{S-1,1}$
$\displaystyle+$
$\displaystyle\frac{4i\varepsilon\eta}{S^{3/2}}\ket{S-1,1}\bra{S,0}$
$\displaystyle-$
$\displaystyle\frac{2\sqrt{2}\varepsilon^{2}}{\sqrt{S(S-1)}}\ket{S,0}\bra{S-2,2}+\textrm{h.c.}.$

To second order, we find
$\mathcal{N}_{\gamma}=\frac{4\varepsilon^{2}}{\gamma^{2}S}+1$. From the
perturbative solution of eq.(19), we find that

$R=S-1+\frac{S^{2}\gamma^{2}}{4\varepsilon^{2}}+O(\gamma)\approx\frac{S^{2}\gamma^{2}}{4\varepsilon^{2}}.$
(24)

This approximation also implies the assumption $\varepsilon^{2}\ll
S\gamma^{2}$. Also, the Logarithmic negativity can be calculated as there are
at most four non-zero eigenvalue for partially transposed density matrix for
any $S$ value. We find to leading order

$E_{ln}=\frac{4\varepsilon}{\gamma\sqrt{S}\ln
2}+O(\frac{1}{\gamma^{2}})\approx\frac{2\sqrt{S}}{\ln
2}\frac{1}{\sqrt{R}}+O(\frac{1}{R}).$ (25)

As noted in Sec. II, the $E_{ln}(R)$ power law behaviour was verified
numerically in Fig. 4.

####  III.2.2 Large dephasing limit

Here we consider the large $\eta$ limit with fixed $S,\gamma,\varepsilon$. We
write the density matrix as a perturbative series in $1/\eta$

$\hat{\rho}_{S}=\frac{1}{\mathcal{N}_{\eta}}(\hat{\rho}^{(0)}_{S}+\frac{1}{\eta}\hat{\rho}^{(1)}_{S}+\frac{1}{\eta^{2}}\hat{\rho}^{(2)}_{S})+O(1/\eta^{3}),$
(26)

where here $\mathcal{N}_{\eta}$ is a normalization constant ensuring the
truncated density matrix has trace $1$. Again, the perturbative series implies
the order by order steady state solutions

$\displaystyle 0=$
$\displaystyle\mathcal{L}_{\mathrm{D}}(\hat{\rho}^{(0)}_{S})$ (27)
$\displaystyle 0=$
$\displaystyle\frac{1}{\eta}\mathcal{L}_{\mathrm{D}}(\hat{\rho}^{(1)}_{S})+(\mathcal{L}_{\mathrm{E}}+\mathcal{H}_{tb})\hat{\rho}^{(0)}_{S}$
(28) $\displaystyle 0=$
$\displaystyle\frac{1}{\eta}\mathcal{L}_{\mathrm{D}}(\hat{\rho}^{(2)}_{S})+(\mathcal{L}_{\mathrm{E}}+\mathcal{H}_{tb})\hat{\rho}^{(1)}_{S}$
(29)

Eq.(27) admits a degenerate solution

$\hat{\rho}^{(0)}_{S}=\sum^{S}_{k=0}a_{k}\ket{k,S-k}\bra{k,S-k}$ (30)

with $a_{k}$ non-negative coefficients. This degeneracy is broken in the next
order, i.e. eq.(28). We find $\hat{\rho}^{(0)}_{S}=\ket{S,0}\bra{S,0}$,
however the degeneracy moves to the next order

$\displaystyle\hat{\rho}^{(1)}_{S}=$
$\displaystyle\sum^{S}_{k=0}b_{k}\ket{k,S-k}\bra{k,S-k}$ (31) $\displaystyle+$
$\displaystyle
i\varepsilon\sqrt{S}(\ket{S,0}\bra{S-1,1}-\ket{S-1,1}\bra{S,0}),$

where $b_{k}$ are again non-negative coefficients. To evaluate to leading
order $R$, we have to break the degeneracy in $b_{k}$. This breaking is
obtained at the next order, i.e. eq.(29), where we find
$b_{k}=\delta_{k,S-1}\frac{2\varepsilon^{2}}{\gamma}$ and

$\displaystyle\hat{\rho}^{(2)}_{S}=$
$\displaystyle-\frac{\varepsilon^{2}\sqrt{S-1}\sqrt{S}}{2\sqrt{2}}\ket{S-2,2}\bra{S,0}$
$\displaystyle+$ $\displaystyle i\frac{\gamma^{2}\varepsilon
S^{3/2}+4\varepsilon^{3}\sqrt{S}}{2\gamma}\ket{S-1,1}\bra{S,0}$
$\displaystyle-$ $\displaystyle
i\frac{2\sqrt{2}\varepsilon^{3}\sqrt{S-1}}{\gamma}\ket{S-2,2}\bra{S-1,1}$
$\displaystyle+$ $\displaystyle\sum^{S}_{k=0}c_{k}\ket{k,S-k}\bra{k,S-k}$
$\displaystyle+\textrm{h.c.}.$ (32)

The degeneracy in the non-negative terms $c_{k}$ is broken at the third order
of the expansion. To leading order in $\eta$, we find
$\mathcal{N}_{\eta}=1+\frac{2\varepsilon^{2}}{\gamma\eta}$. Therefore, to
leading order

$R=\frac{\gamma\eta S}{2\varepsilon^{2}}+S-1\approx\frac{\gamma\eta
S}{2\varepsilon^{2}}.$ (33)

Again, the spectrum of the partially transposed density matrix is composed of
only four non-zero eigenvalues for any $S>1$:
$(-\frac{\gamma\varepsilon\sqrt{S}}{\gamma\eta+2\varepsilon^{2}},\frac{\gamma\varepsilon\sqrt{S}}{\gamma\eta+2\varepsilon^{2}},\frac{2\varepsilon^{2}}{\gamma\eta+2\varepsilon^{2}},\frac{\gamma\eta}{\gamma\eta+2\varepsilon^{2}})$.
The Logarithmic negativity is thus given by

$E_{\ln}=\log_{2}\left(1+\frac{2\gamma\varepsilon\sqrt{S}}{\gamma\eta+2\varepsilon^{2}}\right)\approx\frac{2\varepsilon\sqrt{S}}{\eta\ln
2}=\frac{\gamma S^{3/2}}{\varepsilon\ln 2}\frac{1}{R}.$ (34)

As noted in Sec. II, the $E_{ln}(R)$ power law behaviour was verified
numerically in Fig. 5.

###  III.3 Large number of particles

We also studied the scaling limit $S\gg 1$ and finite
$\eta,\gamma,\varepsilon$. Analytically, a perturbative solution in this case
becomes hard due to the change in the state space, similarly to the Bose-
Hubbard case. Nevertheless, it is possible to numerically find the steady
state and calculate the Logarithmic negativity even for large $S$ values. This
was carried out numerically (Fig. 6) and reported on in Sec. II.

##  IV Discussion

State of the art experimental techniques allows to entangle a single agent to
thousands of atoms McConnell, R., Zhang, H., Hu, J. et al. (2015). However, it
was unclear whether one could push the experimental techniques to
significantly increase the number of atoms entangled to the agent.

Here, we have explored the theoretical bounds on entangling one or a few
agents to a many body system. The ground state of a two-site Bose-Hubbard
model, with an occupancy bias $R\gg 1$ leads to a power law decay in the
Logarithmic negativity and the Von Neumann entanglement entropy in different
scaling limits. Furthermore, the steady state of the QASIP biased to large $R$
values also exhibits a power law decay in the Logarithmic negativity. We
stress that while the power law behaviour is typical, the exponent depends on
the scaling limits, see Tables 1 and 2.

From the slow decay of the entanglement, it is now clear it is typically
possible to entangle thousands of atoms to a single agent. Furthermore,
designing systems with slow entanglement decay (small $\alpha$) allows to
entangle more particles in the many body system to the one agent (or a few).
It would be particularly appealing to develop a perturbative approach to the
large $N$ limit in both models. Such an approach would allow to extract the
exponents analytically, explore their range and dependence on the model
parameters. Furthermore, it would suggest how best to tune the parameters to
entangle the diluted system $B$ to the highly occupied system $A$.

The average Von Neumann entanglement entropy over the random pure state of
Hilbert space $N\times N$ is $E_{vn}\sim\log N$ Page (1993); Sen (1996).
Therefore, the Von Neumann entanglement entropy of the ground state is
fundamentally different than that of the average. This is not too surprising
when one relates to the area law of ground states in extended systems compared
to the typical volume law. In turn, the low Von Neumann entanglement entropy
of the ground state suggests that ground states in the large $R$ limit could
be susceptible to analytical and numerical techniques, even for large many
body systems.

From the analysis so far, it may seem that the two-site lattice model is
paramount to achieve the power law behaviour. We have carried out preliminary
tests in a three site Hubbard model. Taking sites $A,B$ to occupy most of the
particles in the system, namely
$R=\langle\hat{n}_{A}+\hat{n}_{B}\rangle/\langle\hat{n}_{C}\rangle$. The Von
Neumann entanglement entropy between site $C$ and the subsystem $AB$ still
exhibits a power law in large $R$. The analysis is beyond the scope of this
work and will be presented elsewhere.

Another question that comes to mind is whether the power law behaviour
persists also in continuum models, and not only in lattice models. We believe
this is not the case. After coarse graining a lattice model into a continuum
model, an increase is expected in the Von Neumann entanglement entropy due to
loss of information. This increase does not depend on the occupancies and
hence adds a constant to the Von Neumann entanglement entropy. Therefore, in
the large $R$ limit we expect to observe a saturation to a constant with a
power law correction. Naively, that should be the same power law of the
lattice model. It would be interesting to test this conjecture in future
works.

Acknowledgments: I would like to thank Guy Cohen, Shahaf Asban and Ofir E.
Alon for stimulating talks on the subject.

##  Appendix A Entanglement quantifiers

In this section, we provide a brief introduction to the entanglement
quantifiers used in this text: the Von Neumann entanglement entropy and the
Logarithmic negativity. The purpose of quantifiers is to distinguish between
entangled to non-entangled states (separable) and furthermore to suggest a
hierarchy of values for entangled states. Here we do not aim to give an
exhaustive account of quantum quantifiers, but to motivate the usage of the
Von Neumann entanglement entropy and the Logarithmic negativity in the case at
hand.

For pure states, all entanglement measures are defined to correspond to the
Von Neumann entanglement entropy Plenio and Virmani (2007). In bipartite
system $AB$,

$E_{vn}(\rho_{AB})=-\Tr\rho_{A}\ln\rho_{A}=-\Tr\rho_{B}\ln\rho_{B},$ (35)

where $\rho_{A}=\Tr_{B}\rho_{AB}$ is the reduced density matrix. $E_{vn}>0$
only for non-separable pure states.

In terms of wave functions (which are pure states), the Schmidt decomposition
using orthonormal states imply
$\ket{\psi}=\sum_{i}\alpha_{i}\ket{u_{i}}_{A}\otimes\ket{v_{i}}_{B}$. Then, we
find $E_{vn}(\ket{\psi})=-\sum_{i}|\alpha_{i}|^{2}\log|\alpha_{i}|^{2}$.

Entanglement is harder to quantify for mixed states. Many different measures
for the entanglement exists. Typically, entanglement measures are given in the
form of some minimization problem, making them hard to calculate. Instead, we
will use the Logarithmic negativity which is an entanglement monotone and not
a measure. Namely, for pure state the Logarithmic negativity does not
correspond to the Von Neumann entanglement entropy (except for specific
cases). However, it is straight-forward to calculate the Logarithmic
negativity, making it a favorable entanglement quantifier.

The Logarithmic Negativity is given by

$E_{ln}(\rho)=\log_{2}\lVert\rho^{PT}\rVert_{1},$ (36)

where $\rho^{PT}$ is the partially transposed density matrix, and $\lVert
A\rVert_{1}\equiv\Tr\sqrt{AA^{\dagger}}$. Intuitively speaking, the
Logarithmic negativity counts the amount of negative eigenvalues in the
partially transposed density matrix relating it to the Peres–Horodecki
criterion Peres (1996); Horodecki  _et al._ (1996). We note that a positive
Logarithmic negativity values insures non-separability, but a vanishing value
does not guarantees separability.

The Logarithmic Negativity is an entanglement monotone Plenio (2005); Plenio
and Virmani (2007); Friis  _et al._ (2019), which implies that on average,
under locally quantum operations and classical communication (LOCC), the
Logarithmic Negativity does not increase. Furthermore, the Logarithmic
negativity was shown to be an upper bound for the distillation entanglement,
connecting it to useful quantum operations using maximally entangled states
Audenaert  _et al._ (2003). Since the distillation entanglement is an
entanglement measure, it is evident that for pure states, the Von Neumann
entanglement entropy is bounded by the Logarithmic negativity. This fact
provides a consistency check in our numerical assessment.

##  Appendix B The Lindblad adjoint dynamics

The purpose of this section is to introduce the Heisenberg operator evolution
picture for the Lindblad dynamics.

For an observable $\hat{O}$ (explicitly time-independent), we have the
expectation value $\langle\hat{O}\rangle=\Tr\hat{O}\rho$. Therefore,

$\partial_{t}\langle\hat{O}\rangle=\Tr\hat{O}\partial_{t}\rho=\Tr\hat{O}\mathcal{L}(\rho),$
(37)

where $\mathcal{L}(\rho)$ is a Lindblad super-operator of eq.(2). Then, the
formal adjoint $\mathcal{L}^{\dagger}$ is defined such that

$\partial_{t}\langle\hat{O}\rangle=\Tr\mathcal{L}^{\dagger}(\hat{O})\rho.$
(38)

For the Lindblad super-operator in eq.(2), it implies the Heisenberg picture

$\partial_{t}\hat{O}=\mathcal{L}^{\dagger}(\hat{O})=-\mathcal{H}(\hat{O})+\sum_{k}\hat{L}^{\dagger}_{k}\hat{O}\hat{L}_{k}-\frac{1}{2}\\{\hat{L}^{\dagger}_{k}\hat{L}_{k},\hat{O}\\}.$
(39)

It is rather straight-forward to see that $\partial_{t}\hat{N}=0$ for the
QASIP.

##  Appendix C Additional numerical data for the Bose-Hubbard model

Here we present further technical details on the numerical analysis of the
Bose-Hubbard model.

For large $N$ values, the Von Neumann entanglement entropy is easier to obtain
than the logarithmic negativity. In the Appendix A, it was shown that the pure
state Von Neumann entanglement entropy can be obtained from the ket state.
However, the logarithmic negativity requires finding the spectrum of the
partially transposed density matrix. The complexity of handling density
matrices is certainly higher than that of handling ket states, hence the lower
values that were reached for the logarithmic negativity.

In Fig. 7 the scaling $R\propto N^{\beta}$ is presented in the large $N$ limit
for three values of the parameters $(J,U,\mu)$. The values are picked to
produce large $R$ values and to span over a few length scales in $R$,
providing a reliable prediction for the exponents.

Figure 7:  In the large $N$ limit of the two-site Bose-Hubbard model,
numerical evaluation show that $R\propto N^{\beta}$ with $\beta\approx 2.15$.
The three parameters in the range $N\in\left[1,400\right]$, are picked to
allow for a sufficiently large $R$ values.

##  Appendix D Additional numerical data for the QASIP model

Here we present further technical details on the numerical analysis of the
QASIP model.

In Fig. 8 the large $N$ scaling of $R\propto N^{\beta}$ is plotted for
different $\gamma$ values. Since testing the logarithmic negativity in large
$N$ values is numerically challenging, the parameters were chosen to
facilitate as large $R$ as possible which improves exponent fitting.

Figure 8:  In the large $N$ limit of the QASIP model, it is numerically
corroborated that $R\propto N^{\beta}$ where $\beta\approx 1.988$. Different
$\eta,\gamma$ values were tested, leading to the large $R$ limit at
$\varepsilon=1,N\in[1,25]$.

## References

  * Horodecki  _et al._ (2009) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Rev. Mod. Phys. 81, 865 (2009).
  * Amico  _et al._ (2008) L. Amico, R. Fazio, A. Osterloh, and V. Vedral, Reviews of modern physics 80, 517 (2008).
  * Vedral (2018) V. Vedral,  _Decoding reality: the universe as quantum information_ (Oxford University Press, 2018).
  * Nielsen and Chuang (2002) M. A. Nielsen and I. Chuang, “Quantum computation and quantum information,” (2002).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Pezze  _et al._ (2018) L. Pezze, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Reviews of Modern Physics 90, 035005 (2018).
  * Omran  _et al._ (2019) A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, _et al._ , Science 365, 570 (2019).
  * Islam  _et al._ (2015) R. Islam, R. Ma, P. M. Preiss, M. E. Tai, A. Lukin, M. Rispoli, and M. Greiner, Nature 528, 77 (2015).
  * Kaufman  _et al._ (2016) A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Science 353, 794 (2016).
  * Kallush  _et al._ (2021) S. Kallush, R. Dann, and R. Kosloff, arXiv preprint arXiv:2107.11767 (2021).
  * Chruściński and Sarbicki (2014) D. Chruściński and G. Sarbicki, Journal of Physics A: Mathematical and Theoretical 47, 483001 (2014).
  * Lacroix (2020) D. Lacroix, Physical Review Letters 125, 230502 (2020).
  * Weimer  _et al._ (2010) H. Weimer, M. Müller, I. Lesanovsky, P. Zoller, and H. P. Büchler, Nature Physics 6, 382 (2010).
  * Carr and Saffman (2013) A. Carr and M. Saffman, Physical Review Letters 111, 033607 (2013).
  * Verstraete  _et al._ (2009) F. Verstraete, M. M. Wolf, and J. I. Cirac, Nature physics 5, 633 (2009).
  * Eisert  _et al._ (2008) J. Eisert, M. Cramer, and M. B. Plenio, arXiv preprint arXiv:0808.3773 (2008).
  * Page (1993) D. N. Page, Physical Review Letters 71, 1291 (1993).
  * Sen (1996) S. Sen, Physical Review Letters 77, 1 (1996).
  * McConnell, R., Zhang, H., Hu, J. et al. (2015) McConnell, R., Zhang, H., Hu, J. et al. , Nature 519, 439–442 (2015).
  * Friis  _et al._ (2019) N. Friis, G. Vitagliano, M. Malik, and M. Huber, Nat. Rev. Phys. 1, 72 (2019).
  * Plenio and Virmani (2007) M. B. Plenio and S. Virmani, Quant.Inf.Comput 7 (2007).
  * Eisler (2011) V. Eisler, Journal of Statistical Mechanics: Theory and Experiment 2011, P06007 (2011).
  * Bernard  _et al._ (2018) D. Bernard, T. Jin, and O. Shpielberg, EPL (Europhysics Letters) 121 (6), 60006 (2018).
  * Grosskinsky  _et al._ (2011) S. Grosskinsky, F. Redig, and K. Vafayi, J. Stat. Phys. 142, 952 (2011).
  * Fisher  _et al._ (1989) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
  * Jaksch  _et al._ (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Physical Review Letters 81, 3108 (1998).
  * Greiner  _et al._ (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, nature 415, 39 (2002).
  * Links  _et al._ (2006) J. Links, A. Foerster, A. P. Tonel, and G. Santos, in  _Annales Henri Poincaré_ , Vol. 7 (Springer, 2006) pp. 1591–1600.
  * Spekkens and Sipe (1999) R. Spekkens and J. Sipe, Physical Review A 59, 3868 (1999).
  * Shpielberg (2020) O. Shpielberg, EPL (Europhysics Letters) 129, 60005 (2020).
  * Petruccione, F. and Breuer, H.-P. (2002) Petruccione, F. and Breuer, H.-P.,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Gorini  _et al._ (1976) V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Journal of Mathematical Physics 17, 821 (1976).
  * Lindblad (1976) G. Lindblad, Communications in Mathematical Physics 48, 119 (1976).
  * Spohn, Herbert and Lebowitz, Joel L. (1978) Spohn, Herbert and Lebowitz, Joel L., J. Adv. Chem. Phys 38, 109 (1978).
  * Kosloff, Ronnie (2013) Kosloff, Ronnie, Entropy 15 (6), 2100 (2013).
  * Dutta and Cooper (2019) S. Dutta and N. R. Cooper, Phys. Rev. Lett. 123, 250401 (2019).
  * Plenio  _et al._ (2005) M. B. Plenio, J. Eisert, J. Dreissig, and M. Cramer, Physical review letters 94, 060503 (2005).
  * Fischer  _et al._ (2016) M. H. Fischer, M. Maksymenko, and E. Altman, Phys. Rev. Lett. 116, 160401 (2016).
  * Peres (1996) A. Peres, Phys.Rev.Lett. 77, 1413 (1996).
  * Horodecki  _et al._ (1996) M. Horodecki, P. Horodecki, and R. Horodecki, Physics Letters A 223, 1 (1996).
  * Plenio (2005) M. B. Plenio, Phys. Rev. Lett. 95, 090503 (2005).
  * Audenaert  _et al._ (2003) K. Audenaert, M. B. Plenio, and J. Eisert, Phys. Rev. Lett. 90, 027901 (2003).