# Extreme sensitivity of the magnetic ground-state to halide composition in
FeCl3-xBrx

Andrew Cole  Department of Physics, Boston College, Chestnut Hill, MA 02467,
USA  Alenna Streeter  Department of Physics, Boston College, Chestnut Hill, MA
02467, USA  Adolfo O. Fumega  Department of Applied Physics, Aalto University,
Espoo, Finland  Xiaohan Yao  Department of Physics, Boston College, Chestnut
Hill, MA 02467, USA  Zhi-Cheng Wang  Department of Physics, Boston College,
Chestnut Hill, MA 02467, USA  Erxi Feng  Neutron Scattering Division, Oak
Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA  Huibo Cao  Neutron
Scattering Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee
37831, USA  Jose L. Lado  Department of Applied Physics, Aalto University,
Espoo, Finland  Stephen E. Nagler  Neutron Scattering Division, Oak Ridge
National Laboratory, Oak Ridge, Tennessee 37831, USA  Fazel Tafti  Department
of Physics, Boston College, Chestnut Hill, MA 02467, USA

###### Abstract

Mixed halide chemistry has recently been utilized to tune the intrinsic
magnetic properties of transition-metal halides – one of the largest families
of magnetic van der Waals materials. Prior studies have shown that the
strength of exchange interactions, hence the critical temperature, can be
tuned smoothly with halide composition for a given ground-state. Here we show
that the ground-state itself can be altered by a small change of halide
composition leading to a quantum phase transition in FeCl3-xBrx. Specifically,
we find a three-fold jump in the Néel temperature and a sign change in the
Weiss temperature at $x=0.08$ corresponding to only $3\%$ bromine doping.
Using neutron scattering, we reveal a change of the ground-state from spiral
order in FeCl3 to A-type antiferromagnetic order in FeBr3. Using first-
principles calculations, we show that a delicate balance between nearest and
next-nearest neighbor interactions is responsible for such a transition. These
results support the proximity of FeCl3 to a spiral spin liquid state, in which
competing interactions and nearly degenerate magnetic $k$-vectors may cause
large changes in response to small perturbations.

##  I Introduction

Magnetic frustration provides a fascinating playground for the realization of
exotic quantum statesSavary and Balents (2016). A curious example of
frustrated magnet is the spiral spin-liquid (SSL) phase produced by competing
interactions in a bipartite lattice such as diamond (3D) and honeycomb (2D)
structures. The SSL is characterized by degenerate spin spirals with
$\mathbf{k}$-vectors lying on a surface in momentum space Yao  _et al._
(2021). Weak thermal fluctuations can entropically lift this degeneracy and
establish order by disorder Villain  _et al._ (1980). Spinel materials such as
MnSc2S4 and CoAl2O4 are candidates of SSL in the 3D diamond lattice described
by a frustrated $J_{1}$-$J_{2}$ Heisenberg model Fritsch  _et al._ (2004);
MacDougall  _et al._ (2011); Gao  _et al._ (2017); Bergman  _et al._ (2007);
Niggemann  _et al._ (2019). In 2D, however, the experimental realization of
such effects in the honeycomb lattice Mulder  _et al._ (2010); Zhang and Lamas
(2013), prevalent in van der Waals (VdW) materials, has remained elusive.

Recent neutron scattering experiments on the VdW magnet FeCl3 with a honeycomb
lattice have shown a ring of degenerate $\mathbf{k}$-vectors just above
$T_{\textrm{N}}$$=8.5$ K, indicating a 2D SSL phase Gao  _et al._ (2022).
Below $T_{\textrm{N}}$, a spiral order with
$\mathbf{k}=(\frac{4}{15},\frac{1}{15},\frac{3}{2})$ is established Cable  _et
al._ (1962); Jones  _et al._ (1969) indicating the entropic selection of this
propagation vector by spin fluctuations, i.e. order by disorder.

In this letter, we reveal the extreme sensitivity of the ground-state of FeCl3
to tiny amounts of disorder by characterizing a series of FeCl3-xBrx crystals
(Fig. 1a–c). We find a three-fold jump in the Néel temperature
($T_{\textrm{N}}$) and a sign change in the Weiss temperature
($\Theta_{\textrm{W}}$) between $x=0$ and $0.08$ corresponding to only $3\%$
bromine doping. Such dramatic enhancement of the order due to a tiny amount of
disorder suggests that FeCl3 is a frustrated magnet on the verge of a quantum
phase transition (QPT) between the SSL phase and a different ordered state
(Fig. 1d). For $x\geq 0.08$, we identify the ordered state to be A-type AFM
which is ferromagnetic (FM) within the layers and AFM between them. This is
similar to the ground-state of FeBr3 and different from the spiral order in
FeCl3.

Our experimental results are corroborated by first-principles calculations on
a 2D $J_{1}$-$J_{2}$ Heisenberg model that predict a QPT between the SSL and
FM states in FeCl3-xBrx as observed experimentally. We discuss the importance
of including $p$-orbital correlations, in addition to $d$-orbital
correlations, for obtaining the correct magnetic ground-state in density
functional theory (DFT). To our knowledge, the impact of $p$-orbital
correlations in DFT calculations for VdW materials has not be discussed in the
literature before. Details of crystal growth, neutron scattering, scanning
electron microscopy (SEM), magnetization measurements, and DFT calculations
are explained in the Supplementary Information (SI).

##  II Results and Discussion

Figure 1:  (a) Schematic illustration of the honeycomb $ab$-planes in the
mixed halide system FeCl1.5Br1.5. The $J_{1}$ and $J_{2}$ exchange paths are
highlighted in orange and blue colors, respectively. (b) Layered (VdW)
structure of FeCl1.5Br1.5 viewed from the [210] direction. (c) SEM image and
EDX color maps confirming the homogeneous distribution of halides in an
FeCl1.33Br1.67 crystal. (d) Phase diagram of the spiral and FM states with
FeCl3 and FeBr3 across the QPT. Here, the $(q,q,0)$ $y$-label is the
degenerate magnetic wave vector for the SSL state, not the final spiral
ground-state Gao  _et al._ (2022). FeBr3 is not in the SSL phase, and its
order is FM in 2D.

Among different classes of VdW magnets, transition-metal halides offer a
special opportunity in that their intrinsic properties can be tuned by mixing
the halide species (Cl, Br, and I) Abramchuk  _et al._ (2018). With increasing
halide size, the orbital overlaps and ligand spin-orbit coupling (SOC) are
enhanced, which in turn tune all magneto-optical properties as reported
previously in CrCl3-xBrx, CrBr3-yIy, and CrCl3-x-yBrxIy alloys Abramchuk  _et
al._ (2018); Tartaglia  _et al._ (2020).

Following the recent report of a SSL ground-state and order by disorder in
FeCl3, Gao  _et al._ (2022) we grew the heteroanionic crystals of FeCl3-xBrx
with the goal of tuning the SSL ground-state. Since both FeCl3 and FeBr3
crystallize in the rhombohedral space group $R\bar{3}m$ with a layered
honeycomb structure, a solid solution of FeCl3-xBrx crystals can also be grown
with the same structure (Figs. 1a,b). The composition and uniform distribution
of elements were confirmed using energy dispersive x-ray spectroscopy (EDX).
For example, Fig. 1c shows the distribution of Fe $K$-edge (yellow), Cl
$K$-edge (green), and Br $L$-edge (red) absorption intensities in
FeCl1.33Br1.67.

The super-exchange interactions between Fe3+ ions are mediated by Cl and Br
ligands within the honeycomb layers of each FeCl3-xBrx sample. The nearest
neighbor (NN) and next-nearest neighbor (NNN) exchange paths are denoted by
$J_{1}$ and $J_{2}$ in Fig. 1a, respectively. Although $J_{2}$ is weaker than
$J_{1}$, there are six NNNs and three NNs, leading to a close competition
between the $J_{1}$ and $J_{2}$ terms in the model Hamiltonian

$H=J_{1}\sum_{\langle
ij\rangle}\vec{S}_{i}\cdot\vec{S}_{j}+J_{2}\sum_{\langle\langle
ij\rangle\rangle}\vec{S}_{i}\cdot\vec{S}_{j}$ (1)

where $S=5/2$ in the high-spin state of Fe3+ (${}^{6}S_{5/2}$). The
competition between different magnetic exchange paths (direct and super-
exchange) leads to effective $J_{1}<0$ (FM) and $J_{2}>0$ (AFM) for the 2D
spin model of Eq. 1. The ratio between $J_{1}$ and $J_{2}$,
$\left|J_{1}/J_{2}\right|$, determines the magnetic ground-state of the
honeycomb lattice as well as the degenerate $(q,q,0)$ wave vector of the SSL
phase above $T_{\textrm{N}}$ which is shown schematically in Fig. 1d. It is
known from prior neutron diffraction experiments Cable  _et al._ (1962) that
spin spiral is the ground-state of FeCl3. We will show here that the Br
substitution modifies $|J_{1}/J_{2}|$ ratio and drives the system to an A-type
AFM state (which is FM in 2D) where FeBr3 is located in Fig. 1d.

###  II.1 Magnetization Measurements

Figure 2:  (a) Magnetic susceptibility as a function of temperature measured
under ZFC (open circles) and FC (full circles) conditions with $H\|c$. (b)
Same as in (a) but with $H\perp c$. (c) Magnetization curves with $H\|c$
showing metamagnetic (MM) transitions in FeCl3-xBrx samples with $x\geq 0.08$.
(d) The MM transitions are absent when $H\perp c$. (e) $T_{\textrm{N}}$ as a
function of Br content ($x$) showing an initial jump followed by a linear
decrease. (f) $\Theta_{\textrm{W}}$ as a function of $x$ showing an initial
sign change followed by a linear decrease. (g) $\mu_{\text{eff}}$ estimated
from the Curie-Weiss analysis (Fig. S1). Error bars in panels (e,f,g) are
mainly due to the uncertainty in evaluating the mass of thin VdW crystals. (h)
The critical field of MM transitions as a function of $x$ showing an initial
jump followed by a smooth increase. Error bars reflect the width of the
transition (inset).

The experimental evidence of competing FM and AFM interactions appear in the
magnetic susceptibility ($\chi$) and magnetization ($M$) data. The raw data
for all samples with $H\|c$ and $H\perp c$ are shown in Figs. 2a-d, and the
analyzed results are presented in Figs. 2e-h. The empty and full circles in
Figs. 2a,b correspond to zero-field-cooled (ZFC) and field-cooled (FC)
measurements, respectively. In Fig. 2a, the peaks in $\chi(T)$ curves with
$H\|c$ and the minimal difference between the ZFC and FC curves are
characteristics of AFM transitions. Unlike $H\|c$, the $H\perp c$ curves in
Fig. 2b do not go to zero when $T\to 0$, indicating a finite FM component. A
combination of FM and AFM correlations also exists in CrCl3, which undergoes
an A-type AFM order (FM within the layers and AFM between them) Cable  _et
al._ (1961). In CrCl3, only the $\mathrm{t_{2g}}$ manifold of Cr3+ is at half-
filling, whereas both $\mathrm{t_{2g}}$ and $\mathrm{e_{g}}$ levels are at
half-filling in FeCl3, maximizing the competition between FM and AFM
correlations according to Goodenough-Kanamori rules Goodenough (1958);
Kanamori (1959).

We determined $T_{\textrm{N}}$ from $d\chi/dT$ curves (Supplementary Fig. S1),
and plotted it as a function of bromine content ($x$ in FeCl3-xBrx) for both
field directions in Fig. 2e. Using a Curie-Weiss (CW) analysis (Supplementary
Fig. S1), we extracted $\Theta_{\textrm{W}}$, which is a rough measure of the
magnetic correlations, and plotted it as a function of $x$ in Fig. 2f. The
central observation in Figs. 2e,f is a jump in both $T_{\textrm{N}}$ and
$\Theta_{\textrm{W}}$ when a tiny amount of Br is added to FeCl3, i.e. at
$x=0.08$ in FeCl3-xBrx corresponding to only $3\%$ Br doping. The effect is
dramatic with $T_{\textrm{N}}$ showing a three-fold jump and
$\Theta_{\textrm{W}}$ changing sign, indicating a change of magnetic ground-
state at $x=0.08$. From the CW analysis, we found the effective magnetic
moment of all FeCl3-xBrx samples to be close to 5.9 $\mu_{B}$ within
experimental errors (Fig. 2g) as expected for Fe3+ in the high-spin state.

After the initial three-fold jump of $T_{\textrm{N}}$ from $8.5(2)$ K at $x=0$
to $23.1(1)$ K at $x=0.08$, it is suppressed linearly to $14.5(5)$ K at $x=3$
(Fig. 2e). This behavior is the same for both field directions. Similarly,
after the initial jump of $\Theta_{\textrm{W}}$ from $-13$ K at $x=0$ to $+21$
K at $x=0.08$, it is suppressed linearly to $8$ K at $x=3$ for both field
directions (Fig. 2f). Two conclusions can be drawn from these observations.
(i) The initial jump in $T_{\textrm{N}}$ and sign change in
$\Theta_{\textrm{W}}$ suggest an abrupt change of the magnetic ground-state of
FeCl3 by a tiny amount in Br doping, indicating its proximity to a QPT. (ii)
The linear decrease of both $T_{\textrm{N}}$ and $\Theta_{\textrm{W}}$ from
$x=0.08$ to $3$ suggests that the ground-states of all FeCl3-xBrx samples with
$x\geq 0.08$ are similar to the ground-state of FeBr3 and different from that
of FeCl3.

Another evidence of the change of magnetic ground-state at $x=0.08$ comes from
the field dependence of magnetization. The $M(H)$ curves of FeCl3 in Figs.
2c,d are nearly linear in both field directions, consistent with the spiral
AFM order reported in prior neutron diffraction studies Cable  _et al._
(1962). The subtle kinks in the $M(H\perp c)$ curves at $H=0$ and $4$ T are
due to the alignment of spiral domains with the field. Unlike FeCl3, the
FeCl3-xBrx samples with $x\geq 0.08$ exhibit field-induced metamagnetic (MM)
transitions when $H\perp c$ and linear $M(H)$ when $H\|c$ (Figs. 2c,d).
Tracing the critical field ($H_{c}$) of the MM transition as a function of $x$
(Fig. 2h) reveals an initial jump at $x=0.08$ followed by a smooth increase of
$H_{c}$ for $x\geq 0.08$. The $H_{c}$ values in Fig. 2h correspond to the peak
fields in the $dM/dH$ curves shown in the inset. Similar to the behavior of
$T_{\textrm{N}}$ and $\Theta_{\textrm{W}}$ (Figs. 2e,f), the initial jump of
$H_{c}$ at $x=0.08$ in Fig. 2h indicates a change of the magnetic ground-
state, and its subsequent smooth increase indicates that the ground-states of
FeCl3-xBrx samples with $x\geq 0.08$ are similar to that of FeBr3 and
different from FeCl3.

###  II.2 Neutron Diffraction.

Figure 3:  (a) The 2D neutron diffraction scan at 4.8 K and zero-field in the
$H$-$L$ plane at $K=\bar{1}$, showing strong nuclear reflections at integer
$L$ and weak magnetic reflections at half-integer $L$. (b) A 1D cut through
the data in panel (a) showing nuclear (strong) and magnetic (weak) Bragg peaks
that identify the magnetic propagation vector $\mathbf{k}=(0,0,1.5)$. (c) 2D
scan in the $(H\bar{1}L)$ plane showing another set of nuclear and magnetic
reflections. (d) 1D cut through the data in panel (c) showing Bragg peaks
along $(1\bar{1}L)$ direction. (e) 2D scan in the $(H0L)$ plane showing the
absence of magnetic reflections along $(00L)$ direction. (f) Temperature
dependence of the intensity of the magnetic reflection
$\mathbf{Q}=(2,\bar{1},4.5)$ with a power-law fit to extract $T_{\textrm{N}}$
and the critical exponent $\beta$. (g,h,i) Candidate magnetic ground-states of
FeBr3 with $R_{l}$ symmetry. (j,k,l) Candidate ground-states with $P_{s}$
symmetry. The ground-state is determined as A-type AFM shown in panel (g).

Single crystal neutron diffraction has been performed on FeCl3 previously, and
the ground-state was determined as AFM with a spiral modulation parallel to
the $[14\bar{5}]$ direction with a wavelength of 15 $(14\bar{5})$ $d$ spacing
Cable  _et al._ (1962). We performed single-crystal neutron diffraction on
FeBr3 at 4.8 K and zero field to probe its magnetic ground-state. Since all
FeCl3-xBrx samples with $x\neq 0$ have similar $\chi(T)$ and $M(H)$ behaviors
as FeBr3 (Fig. 2), we assume that their magnetic ground-states are similar to
that of FeBr3.

The observed Bragg peaks in Figs. 3a,b at $\mathbf{Q}=(2,\bar{1},L)$ consist
of structural (nuclear) peaks at $L=3$ and $6$ as well as magnetic peaks at
$L=\frac{3}{2}$, $3+\frac{3}{2}$, and $6+\frac{3}{2}$, hence the magnetic
propagation vector $\mathbf{k}=\left(0,0,1.5\right)$. The systematic absence
of the nuclear peaks other than $L=3n$ ($n\in$ integer) in Figs. 3a,b is
expected in the space group $R\bar{3}$ (#148) of FeBr3. A second scan along
$\mathbf{Q}=(1,\bar{1},L)$ in Figs. 3c,d confirms the propagation vector
$\mathbf{k}=\left(0,0,1.5\right)$. The absence of magnetic reflections at
$\mathbf{Q}=(0,0,4.5)$ and $(0,0,7.5)$ in Fig. 3e suggests that the ordered
magnetic moments lie along the $c$-axis because neutrons only probe the moment
perpendicular to scattering vector ($\mathbf{M}\perp\mathbf{Q}$). Thus, we
identify FeBr3 as an out-of-plane Ising system with $M_{z}\neq 0$ and
$M_{x,y}=0$.

We construct an order parameter plot in Fig. 3f by tracing the intensity of
the $\mathbf{Q}=(2,\bar{1},4.5)$ peak as a function of temperature. From a
power-law fit, we extract $T_{\text{N}}=14.0(5)$ K consistent with $14.5(5)$ K
from the magnetization measurements, and $\beta=0.324(15)$ consistent with a
3D Ising system.

A magnetic symmetry analysis based on the ordering wave-vector
$\mathbf{k}=\left(0,0,1.5\right)$ for $S=5/2$ Fe3+ ions in the structural
space group $R\bar{3}$ of FeBr3 identifies six possible ground-states. These
states are labeled in Figs. 3g-l by their magnetic subgroups and irreducible
representations (irreps) as $R_{l}\bar{3}$ (irrep: $mT_{1-}$), $R_{l}\bar{3}$
($mT_{1+}$), $R_{l}3$ ($mT_{1-}$), $P_{s}\bar{1}$ ($mT_{1-}$), $P_{s}\bar{1}$
($mT_{1+}$), and $P_{s}1$ ($mT_{1-}$). The first two candidates in Figs. 3g
and 3h represent maximal symmetry subgroups corresponding to out-of-plane
Ising spins ($M_{z}\neq 0$ and $M_{x,y}=0$) with A-type and C-type AFM
ordering, respectively. The lower symmetry $R_{l}3$ similarly allows only
$M_{z}$ components but with two different moment sizes (Fig. 3i). The later
three subgroups (Figs. 3j,k,l) are primitive ($P$) instead of rhombohedral
($R$). Having a lower symmetry, they allow for each of the previous cases in
Figs. 3g,h,i to have moments canted toward the $ab$-plane, i.e. $M_{x,y}\neq
0$.

We identify the ground-state of FeBr3 as A-type AFM order (Fig. 3g) for the
following reasons. (i) The lack of magnetic reflections along the $00L$
direction in Fig. 3e indicates that the allowed magnetic moments lie along the
$c$-axis, which disqualifies the configurations in Figs. 3j,k,l with allowed
moments in the $ab$-plane. (ii) Our neutron scattering refinement of the
magnetic moment size (SI) gives 4.90(99) $\mu_{B}$ on each Fe3+ site,
consistent with the results of magnetization measurements and slightly larger
than 4.3 $\mu_{B}$ reported for FeCl3 in prior studies Cable  _et al._ (1962);
Jones  _et al._ (1969). Thus, the configuration in Fig. 3i with different
moment sizes on different Fe-sites is also disqualified. (iii) The remaining
candidates in Figs. 3g and 3h belong to the maximal subgroup $R_{l}\bar{3}$
corresponding to the A-type and C-type AFM ordering, respectively. These two
structures can be distinguished according to the intensity of magnetic
diffraction peaks. We simulated the structure factor ($F^{2}_{sim}$) for
several nuclear and magnetic Bragg peaks in both A-type and C-type structures
assuming $\mu=4.9~{}\mu_{B}$, and summarized them in the Supplementary Table I
(SI). According to Table I, $F^{2}_{sim}$ values for two magnetic Bragg peaks
are vanishingly small in the C-type AFM structure and considerably large in
the A-type AFM structure. Thus, the diffraction pattern is best described by
A-type AFM ordering (Fig. 3g).

An A-type AFM ground-state is also consistent with the field-scale of the MM
transition in Fig. 2c which is about 3 T in FeBr3 corresponding to 2 K,
considerably smaller than $T_{\text{N}}=14$ K. Thus, the MM transitions are
likely due to a spin flip between the layers instead of within the layers, so
the ground-state must be A-type instead of C-type AFM.

###  II.3 First-Principles Calculations.

Figure 4:  (a) The NN ($J_{1}$) and NNN ($J_{2}$) coupling constants computed
from first-principles and plotted as a function of in-plane lattice parameter
$a$ in stoichiometric FeCl3. (b) The ratio $|J_{1}/J_{2}|$ plotted as a
function of $a$. (c,d) Same as in (a,b) but for the stoichimetric FeBr3. (e)
$J_{1}$ and $J_{2}$ traced as a function of bromine content $x$ in FeCl3-xBrx
using virtual lattice approximation. (f) $|J_{1}/J_{2}|$ as a function of $x$.

Our experimental results suggest that FeCl3 must be very close to a QPT since
the magnetic ground-state changes from the spiral order in FeCl3 to an A-type
AFM order upon $3\%$ bromine doping ($x=0.08$ in FeCl3-xBrx). To understand
the nature of the transition, we performed all-electron first-principles
calculations based on DFT in FeCl3 and FeBr3 to extract $J_{1}$ and $J_{2}$
exchange couplings in the model Hamiltonian Eq. 1 for a 2D honeycomb lattice
of $S=5/2$ spins (Fig. 1a). We extracted $J_{1}$ and $J_{2}$ as a function of
the in-plane lattice parameters interpolating between the experimental values
of FeCl3 and FeBr3, $a=6.05$ and $6.40$ Å, respectively. As a benchmark to our
methodology, we have obtained $J_{1}$ and $J_{2}$ values in good agreement
with those reported for FeCl3 via inelastic neutron scattering Gao  _et al._
(2022) giving rise to a spiral spin liquid state with $|J_{1}/J_{2}|\approx
4$.

To reproduce the high-spin configuration and $J_{1}/J_{2}$ found
experimentally in FeCl3 and FeBr3, onsite Coulomb interactions must be
included in DFT. We found that $U_{dd}=7.0$, $U_{pp}=0.5$, and $U_{pp}=3.5$ eV
were required in the Fe, Cl, and Br atoms, respectively, to obtain the SSL
state in FeCl3 and in-plane FM state in FeBr3 (Supplementary Figs. S2,3). We
verified that this trend was robust against changes in the local interactions,
including spin-orbit coupling effects (SI). The larger value of $p$-orbital
potential in FeBr3 was necessary for establishing the experimentally observed
in-plane FM state. Smaller values of $U_{pp}$ would lead to a much larger
$J_{2}$ and establish an AFM (Néel) state within the honeycomb planes of FeBr3
due to the more covalent nature of Fe-Br bonds compared to Fe-Cl bonds. Note
that the in-plane correlations are FM within an A-type AFM order (Fig. 3g)
consistent with the positive $\Theta_{\textrm{W}}$ observed in FeBr3 (Fig.
2f). Our calculations are performed on a 2D lattice without considering the
inter-layer coupling $J_{c}$ that eventually establishes the 3D A-type AFM
order in FeBr3 (SI).

Figures 4a-d show the evolution of $J_{1}$, $J_{2}$ and $|J_{1}/J_{2}|$ as a
function of the lattice parameter for FeCl3 and FeBr3. We observe that for
stoichiometric FeCl3, an increase in the lattice parameter dramatically
impacts $|J_{1}|$ unlike $|J_{2}|$ (Fig. 4a), while the system remains in the
spiral state (Fig. 4b). In contrast, the stoichiometric FeBr3 remains in the
FM state for all lattice parameter values due to a much larger $|J_{1}|$
(Figs. 4c,d). The dashed lines in Figs. 4b,d mark the critical value
$|J_{1}/J_{2}|=2Z=6$ for the theoretical transition from the SSL to FM/Néel
state Bergman  _et al._ (2007); Gao  _et al._ (2017); Niggemann  _et al._
(2019). To account for the alloys FeCl3-xBrx, we use a virtual crystal
approximation (SI) to estimate $J_{1}$ and $J_{2}$ as a function of $x$ based
on the alloy-dependent lattice constant (Figs. 4e,f). Within this
approximation, doping FeCl3 with Br drives the 2D system from the spiral to FM
phase, as observed experimentally.

The theoretical phase diagram in Fig. 4f shows a QPT between the spiral and FM
states at $x=0.8$ in 2D layers of FeCl3-xBrx. This is consistent with the
experimental data in Figs. 2e,f,h that show a jump in $T_{\textrm{N}}$, sign
change in $\Theta_{\textrm{W}}$, and MM transition at $x>0.08$. However, the
theoretically predicted critical doping $x_{c}=0.8$ differs from the
experimentally observed $x_{c}=0.08$. Such a difference likely stems from
limitations of the DFT methodology that neglects the impact of disorder on
exchange interactions. Spin liquid phases (e.g. SSL) are particularly
sensitive to disorder Bergman  _et al._ (2007); Niggemann  _et al._ (2019);
Kao  _et al._ (2021); Dantas and Andrade (2022) and thus, it is conceivable
that the alloying procedure would change the critical doping via disorder
effects that are not captured in the pristine-limit DFT calculations.

##  III Conclusion

To summarize, we have demonstrated a QPT by tuning the halide composition in
the frustrated VdW system FeCl3-xBrx. Our results demonstrate the application
of halide engineering in tuning the $J_{1}/J_{2}$ ratio of the underlying
frustrated honeycomb spin model, in particular crossing the critical point of
the spin spiral liquid at $x_{c}=0.08$. The transition driven by the halide
composition from a spiral to FM state was demonstrated with magnetometry
measurements, and the ground-states were determined by neutron scattering
results presented here for FeBr3 and elsewhere Cable  _et al._ (1962) for
FeCl3. Our theoretical calculations further show that the transition from the
spiral to FM state is driven by competing magnetic exchanges with a sizable
contribution from the $p$-orbital correlations of halides. Our results
establish a new strategy for engineering frustrated VdW magnetic materials by
exploiting a continuous parameter realized by mixed halide chemistry.

## ACKNOWLEDGMENTS

The work at Boston College was supported by the National Science Foundation
under the grant number DMR-2203512. J.L.L. and A.F. acknowledge the
computational resources provided by the Aalto Science-IT project, and the
financial support from the Academy of Finland Projects No. 331342, No. 336243
and No 349696, and the Jane and Aatos Erkko Foundation. A portion of this
research used resources at the High Flux Isotope Reactor, a DOE Office of
Science User Facility operated by the Oak Ridge National Laboratory.

## References

  * Savary and Balents (2016) L. Savary and L. Balents, Reports on Progress in Physics 80, 016502 (2016).
  * Yao  _et al._ (2021) X.-P. Yao, J. Q. Liu, C.-J. Huang, X. Wang, and G. Chen, Frontiers of Physics 16, 53303 (2021).
  * Villain  _et al._ (1980) J. Villain, R. Bidaux, J.-P. Carton, and R. Conte, Journal de Physique 41, 1263 (1980).
  * Fritsch  _et al._ (2004) V. Fritsch, J. Hemberger, N. Büttgen, E.-W. Scheidt, H.-A. Krug von Nidda, A. Loidl, and V. Tsurkan, Physical Review Letters 92, 116401 (2004).
  * MacDougall  _et al._ (2011) G. J. MacDougall, D. Gout, J. L. Zarestky, G. Ehlers, A. Podlesnyak, M. A. McGuire, D. Mandrus, and S. E. Nagler, Proceedings of the National Academy of Sciences 108, 15693 (2011).
  * Gao  _et al._ (2017) S. Gao, O. Zaharko, V. Tsurkan, Y. Su, J. S. White, G. S. Tucker, B. Roessli, F. Bourdarot, R. Sibille, D. Chernyshov, T. Fennell, A. Loidl, and C. Rüegg, Nature Physics 13, 157 (2017).
  * Bergman  _et al._ (2007) D. Bergman, J. Alicea, E. Gull, S. Trebst, and L. Balents, Nature Physics 3, 487 (2007).
  * Niggemann  _et al._ (2019) N. Niggemann, M. Hering, and J. Reuther, Journal of Physics: Condensed Matter 32, 024001 (2019).
  * Mulder  _et al._ (2010) A. Mulder, R. Ganesh, L. Capriotti, and A. Paramekanti, Physical Review B 81, 214419 (2010).
  * Zhang and Lamas (2013) H. Zhang and C. A. Lamas, Physical Review B 87, 024415 (2013).
  * Gao  _et al._ (2022) S. Gao, M. A. McGuire, Y. Liu, D. L. Abernathy, C. d. Cruz, M. Frontzek, M. B. Stone, and A. D. Christianson, Physical Review Letters 128, 227201 (2022).
  * Cable  _et al._ (1962) J. W. Cable, M. K. Wilkinson, E. O. Wollan, and W. C. Koehler, Physical Review 127, 714 (1962).
  * Jones  _et al._ (1969) E. R. Jones, O. B. Morton, L. Cathey, T. Auel, and E. L. Amma, The Journal of Chemical Physics 50, 4755 (1969).
  * Abramchuk  _et al._ (2018) M. Abramchuk, S. Jaszewski, K. R. Metz, G. B. Osterhoudt, Y. Wang, K. S. Burch, and F. Tafti, Advanced Materials 30, 1801325 (2018).
  * Tartaglia  _et al._ (2020) T. A. Tartaglia, J. N. Tang, J. L. Lado, F. Bahrami, M. Abramchuk, G. T. McCandless, M. C. Doyle, K. S. Burch, Y. Ran, J. Y. Chan, and F. Tafti, Science Advances 6, eabb9379 (2020).
  * Cable  _et al._ (1961) J. W. Cable, M. K. Wilkinson, and E. O. Wollan, Journal of Physics and Chemistry of Solids 19, 29 (1961).
  * Goodenough (1958) J. B. Goodenough, Journal of Physics and Chemistry of Solids 6, 287 (1958).
  * Kanamori (1959) J. Kanamori, Journal of Physics and Chemistry of Solids 10, 87 (1959).
  * Kao  _et al._ (2021) W.-H. Kao, J. Knolle, G. B. Halász, R. Moessner, and N. B. Perkins, Physical Review X 11, 011034 (2021).
  * Dantas and Andrade (2022) V. Dantas and E. C. Andrade, Physical Review Letters 129, 037204 (2022).


Unitary (semi)causal quantum-circuit representation
of black hole evaporation
Bogusław Broda
Department of Theoretical Physics, Faculty of Physics and Applied
Informatics, University of Łódź, 90-236 Łódź,
Pomorska 149/153, Poland

ORCID: 0000-0001-9753-7207

A general structure of unitary evolution (evaporation) of the black
hole, respecting causality imposed by the event horizon (semicausality),
has been derived and presented in the language of quantum circuits.
The resulting consequences for the evolution of the corresponding
entanglement entropy and the entropy curve have been determined. As
an illustration of the general scheme, two families of qubit toy models
have been discussed: tensor product models and controlled non-product

§ INTRODUCTION

The black hole (BH) information paradox formulated in [Hawking, 1976]
is still an interesting, unsolved, and inspiring problem. Recent progress,
starting with the two famous papers made in 2019, [Almheiri et al., 2019]
and [Penington, 2020], consists in the computation of the von
Neumann entropy of the radiation emitted by BHs in the framework of
some specific models of quantum gravity. It appears that in these
models, the shape of the entropy curve is essentially given by the
Page curve [Page, 1993]. Instead, in the present work, we concentrate
on the description of the unitary evolution and evaporation process
of BHs directly in terms of qubits, according to the ideas reviewed
by e.g. [Avery, 2013]. Extensive literature with some 50 references
to various qubit models is given in [Osuga and Page, 2018] (see also proposals
in [Broda, 2020, Broda, 2021, Broda, 2021]). Despite these large
studies, it seems that the important issue of properly understood
causality in the process of transport of quantum information in the
presence of the BH event horizon has been largely ignored. The aim
of our work is to at least partially fill the gap.

The paper consists of two parts, and the whole analysis is performed
in the language of quantum circuits. In the first part (Section <ref>),
we implement the concept of causality that is, in our opinion, most
natural for the description of the process of evolution (and evaporation)
of a BH. Namely, taking into account the role played by the BH event
horizon, the kind of causality appropriate in this case is given by
the (already existent in literature) notion of semicausality.
The semicausality in turn imposes some restrictions on the structure
of unitary evolution and, in consequence, on the evolution of entropy
and (in a graphic presentation) on the entropy curve. In the second
part of the paper (Section <ref>), we present
and analyze two large families of unitary qubit toy models that explicitly
implement the idea of (semi)causality. The first family comprises
some tensor product models, whereas the second one is a non-product
generalization of the first family with a small number of control

§ THE GENERAL STRUCTURE OF UNITARY SEMICAUSAL CIRCUITS

In this section, we would like to establish a general structure of
unitary transformations governing the evolution of evaporating BHs,
which properly respect causality implied by the presence of the BH
event horizon. Because of our later needs (see Section <ref>),
the answer will be formulated in the language of quantum circuits.
Since unitary transformations (and their corresponding quantum circuits)
are supposed to act in regions both outside and inside the BH, restrictions
following from causality should be properly encoded in the structure
of the circuit. It appears that the version of causality appropriate
in the context of a BH is semicausality. In quantum mechanics, the
notion of semicausality has a short history. The concept of semicausality
first appeared in [Beckman et al., 2001], was analyzed in [Eggeling et al., 2002]
(see also [Piani et al., 2006]), and was later discussed by [Schumacher and Westmoreland, 2012]
(see also [Schumacher and Westmoreland, 2005]) in the circuit (diagrammatic) representation.
One should also note that semicausality has already explicitly appeared
in the context of BHs in [Braunstein and Pirandola, 2018].

Occasionally, we loosely use the term (semi)causality instead
of semicausality (e.g., in the title) or even causality, but referring
to causality in the present paper, we exclusively mean the notion
of semicausality (precisely, causality is something different—it
is semicausality “valid in both directions”).

§.§ Semicausality in quantum circuits

An inherent characteristic property of the BH is its one-way transmission
of information (and particles) through its event horizon. This property
perfectly harmonizes with the notion of semicausality introduced in
the context of quantum mechanics and is phrased in the following way
[Beckman et al., 2001]: If a bipartite operation $\mathcal{E}$ does
not enable superluminal signalling from Bob to Alice (denoted in literature
as $\textrm{B}\nrightarrow\textrm{A}$, or similarly), then we say
that $\mathcal{E}$ is semicausal. In the interesting context of the
event horizon, Bob should be placed inside the BH, whereas Alice should
reside outside the BH. The main result of [Eggeling et al., 2002] is
a proof of a conjecture by DiVincenco (mentioned by [Beckman et al., 2001]),
which (in physics terms) says that if a device allows no signalling
from Bob to Alice (semicausality), we can represent it explicitly
as a device involving possibly a particle sent from Alice to Bob but
none in the other direction. More precisely (in mathematical terms),
let the unitary evolution of the tripartite system $\textrm{BXA}$
be described by the operator $U^{\left(\textrm{BXA}\right)}$. Then
$U^{\left(\textrm{BXA}\right)}$ can be represented by the following
product containing unitary operators $U^{\left(\textrm{BX}\right)}$
and $U^{\left(\textrm{XA}\right)}$ [Schumacher and Westmoreland, 2005]:

\begin{equation}
U^{\left(\textrm{BXA}\right)}=\left(U^{\left(\textrm{BX}\right)}\varotimes1^{\left(\textrm{A}\right)}\right)\left(1^{\left(\textrm{B}\right)}\otimes U^{\left(\textrm{XA}\right)}\right),\label{eq:U(BXA)_tensored}
\end{equation}

where $\textrm{B}$ denotes Bob's subsystem, $\textrm{X}$ is a “particles'
sent” subsystem, and $\textrm{A}$ denotes Alice's subsystem. A
diagrammatic (circuit in block form) presentation of (<ref>)
is depicted in Figure <ref> (see [Schumacher and Westmoreland, 2012]).
We can also present another, essentially equivalent version of the
theorem, which is better suited for our further considerations:

\begin{equation}
U^{\left(\textrm{BXA}\right)}=\left(U^{\left(\textrm{BX}\right)}\varotimes U^{\left(\textrm{A}\right)}\right)\left(U^{\left(\textrm{B}\right)}\otimes U^{\left(\textrm{XA}\right)}\right),\label{eq:U(BXA)_tensored_an}
\end{equation}

where the additional unitary operators $U^{\left(\textrm{A}\right)}$
and $U^{\left(\textrm{B}\right)}$ only redefine the subsystems $\textrm{A}$
and $\textrm{B}$ after and before the evolution (<ref>),
respectively. A diagrammatic version of (<ref>)
is presented in Figure <ref>.

B         2U^(BX)

X     2U^(XA)

B     U^(B)              2U^(BX)

X     2U^(XA)

A               					    U^(A)

Diagrammatic representation of the operator $U^{\left(\textrm{BXA}\right)}$
defining the unitary evolution of the tripartite system $\textrm{BXA}$,
where no signalling from Bob (B) to Alice (A) is allowed (semicausality).
Transmission of information can only be related to particle(s) (X)
sent from Alice (A) to Bob (B), but not in the other direction. Here
(a) presents the original version of the operator, given in (<ref>),
whereas (b) presents its slightly generalized version (<ref>),
which is better suited for our further considerations.

In the process of BH evaporation (and BH evolution in general), we
deal with a continuous flux of particles (e.g., external real “matter”
particles and/or Hawking's virtual ones) constantly crossing the event
horizon and entering the BH. Therefore, the formula (<ref>)
or (<ref>) should be appropriately generalized,
i.e., iterated. Namely, the unitary evolution of the $N$-partite
system $X_{1}X_{2}\ldots X_{N}$ described by the operator $U^{\left(X_{1}X_{2}\ldots X_{N}\right)}$
should be a (huge, in the BH case) iteration of (<ref>),
assuming the form

\begin{equation}
\begin{aligned}U^{\left(X_{1}X_{2}\ldots X_{N}\right)}= & \left(U^{\left(X_{1}X_{2}\ldots X_{N-1}\right)}\otimes U^{\left(X_{N}\right)}\right)\\
& \ldots\left(U^{\left(X_{1}X_{2}\right)}\otimes U^{\left(X_{3}\ldots X_{N}\right)}\right)\left(U^{\left(X_{1}\right)}\otimes U^{\left(X_{2}\ldots X_{N}\right)}\right)
\end{aligned}
\label{eq:U(many_X)_tensored}
\end{equation}

(see Figure <ref>), where all the operators $U^{(\cdots)}$
are unitary.

$X_1$	   U^(X_1)	     2U^(X_1X_2)     …    4U^(X_1X_2…X_N-1)

$X_2$	   4U^(X_2…X_N)      		     …

$X_3$       		      3U^(X_3…X_N)     …

1.5ex$\vdots$.5em      	   	    …

$X_N$	   		      		     …    U^(X_N)

Diagrammatic representation of the operator $U^{\left(X_{1}X_{2}\ldots X_{N}\right)}$,
which is an arbitrarily large iteration of the operator $U^{\left(\textrm{BXA}\right)}$
defined in (<ref>) and depicted in Figure <ref>.

Going from $N$-partite (hugely multi-particle, in general) quantum
subsystems $X_{1},X_{2},\:$$\ldots,X_{N}$ of
the system $X$ (here the symbol $X$ should not be confused with
a “particles' sent” subsystem X from the beginning of this subsection)
to $n$ qubits $x_{1},x_{2},\ldots,x_{n}$ (usually $n\gg N$), we
can replace (<ref>) with its (most) fine-grained

\begin{equation}
\begin{aligned}U^{\left(x_{1}x_{2}\ldots x_{n}\right)}= & \left(U^{\left(x_{1}x_{2}\ldots x_{n-1}\right)}\otimes U^{\left(x_{n}\right)}\right)\\
& \ldots\left(U^{\left(x_{1}x_{2}\right)}\otimes U^{\left(x_{3}\ldots x_{n}\right)}\right)\left(U^{\left(x_{1}\right)}\otimes U^{\left(x_{2}\ldots x_{n}\right)}\right),
\end{aligned}
\label{eq:U(many_x)_tensored}
\end{equation}

or equivalently (in a shorter notation) by

\begin{equation}
U=\left(U_{n-1}^{\left(1\right)}\otimes U_{n-1}^{\left(2\right)}\right)\ldots\left(U_{2}^{\left(1\right)}\otimes U_{2}^{\left(2\right)}\right)\left(U_{1}^{\left(1\right)}\otimes U_{1}^{\left(2\right)}\right)\equiv U_{n-1}\ldots U_{2}U_{1},\label{eq:U(many_x)_tensored_shorter}
\end{equation}

where $U_{k}\equiv U_{k}^{\left(1\right)}\otimes U_{k}^{\left(2\right)}$
($k=1,\ldots n-1$). The expansion (<ref>)
or equivalently (<ref>) can be a bit
further generalized to

\begin{equation}
\begin{aligned}U= & U_{n}\left(U_{n-1}^{\left(1\right)}\otimes U_{n-1}^{\left(2\right)}\right)\ldots\left(U_{2}^{\left(1\right)}\otimes U_{2}^{\left(2\right)}\right)\left(U_{1}^{\left(1\right)}\otimes U_{1}^{\left(2\right)}\right)U_{0}\\
\equiv & U_{n}U_{n-1}\ldots U_{2}U_{1}U_{0},
\end{aligned}
\label{eq:U(many_x)_tensored_generalized}
\end{equation}

where the “external” unitary operators $U_{0}$ and $U_{n}$ have
been added in analogy to the generalization (<ref>),
and they only redefine the in (initial) state $\left|\psi_{\textrm{in}}\right\rangle $
and the out (final) state $\left|\psi_{\textrm{out}}\right\rangle $,
respectively (see Figure <ref>).

style=very thin,fill=blue!5,inner sep=0em,background]

[wires=3]$\ket{\psi_{\textrm{in}}}$   3U_0≫31$U_0$

…    2U_n-1^(1)≫31$U_{n-1}$

[]3U_n≫31$U_n$    [wires=3]$\ket{ \psi_{\textrm{out}}}$

2U_1^(2)    …

…    U_n-1^(2)

(Color online) Diagrammatic representation of the unitary evolution
operator $U$ defined in (<ref>),
which is a slight generalization of the operator $U^{\left(x_{1}x_{2}\ldots x_{n}\right)}$
in (<ref>), interpreted as the most fine-grained
particle version of the operator $U^{\left(X_{1}X_{2}\ldots X_{N}\right)}$
given in (<ref>) and depicted in Figure <ref>.
The vertical dashed left line and the right line denote the theoretically
earliest moment and the latest one, when the BH is allowed to exist,

In the context of BH evolution, the operators $U_{0}$ and $U_{n}$
are supposed to act before the BH forms and after its complete evaporation,

One should stress that according to (<ref>),
the decomposition $\mathcal{H}=\mathcal{H}_{k}^{\left(1\right)}\otimes\mathcal{H}_{k}^{\left(2\right)}$
($k=1,\ldots,n-1$) of the total Hilbert space $\mathcal{H}$ onto
a product of the two Hilbert spaces $\mathcal{H}_{k}^{\left(a\right)}$,
$a=1,2$, where the superscripts 1 and 2 correspond to the interior
and the exterior of the BH, respectively, is “dynamical”. In other
words, the decomposition of $\mathcal{H}$ gradually changes in the
course of the (indexed by $k$) evolution of the BH, as depicted in
Figures <ref> and <ref>. This
observation undermines the validity of the fixed global decomposition
of the total Hilbert space $\mathcal{H}$ in the context of BH evaporation.

§.§ Bounds on entanglement entropy

The general situation and the area of our interest between two arbitrary
neighbouring unitary blocks $U_{k-1}$ and $U_{k}$ have been presented
in Figure <ref>.

U_k-1                 U_k

$\ket{x_1}$    2U_k-1^(1)∈ℋ_k-1^(1)
X$\ket{\psi_k}$       X$\mathcal{S}_{k}^{+}\dbinom{\rho_{k}^{(1)+}}{\rho_{k}^{(2)+}}$    3    U_k^(1)∈ℋ_k^(1) [wires=2]B

$\ket{x_k}$    3U_k-1^(2)∈ℋ_k-1^(2)    [wires=1]X                 [wires=1]X

1.5ex$\vdots$.5em          [wires=2]A	        	               2    U_k^(2)∈ℋ_k^(2)[wires=2]A

(Color online) Description of the configuration between two arbitrary
britishneighbouring unitary blocks of Figure <ref>.
The (pure) state between the two blocks $U_{k-1}$ and $U_{k}$ is
denoted by $\ket{\psi_{k}}$. In our notation, the correspondence
between the reduced density matrices and the subsystems is as follows:
and $\rho_{k}^{(2)+}$—$\,\left(\textrm{A}\right)\equiv\left(x_{k+1},\ldots,x_{n}\right)$
(see (<ref>) and (<ref>)).
The entanglement entropy $\mathcal{S}_{k}^{\mp}$ is by definition
the von Neumann entropy for the reduced density matrix $\rho_{k}^{\left(a\right)\mp}$
($a=1,2$), where, as it turns out, $a$ can be arbitrary and thus
finally skipped.

The total (pure) state exiting the unitary block $U_{k-1}$ and next
entering $U_{k}$ has been denoted by $\left|\psi_{k}\right\rangle $.
As follows from Figure <ref>, the configuration in
the area between the two pairs of unitary operators, $U_{k-1}\equiv U_{k-1}^{\left(1\right)}\otimes U_{k-1}^{\left(2\right)}$
and $U_{k}\equiv U_{k}^{\left(1\right)}\otimes U_{k}^{\left(2\right)}$,
gives rise to the four density operators $\rho_{k}^{(a)\mp}$, $a=1,2$.
Here the signs “$-$” and “$+$” denote the LHS and the RHS
of the area in Figure <ref>, respectively, and the
values of the superscript $a$ correspond to the values $1$ and $2$
in the two decompositions of the Hilbert space $\mathcal{H}=\mathcal{H}_{k-1}^{\left(1\right)}\otimes\mathcal{H}_{k-1}^{\left(2\right)}$
and $\mathcal{H}=\mathcal{H}_{k}^{\left(1\right)}\otimes\mathcal{H}_{k}^{\left(2\right)}$
for the LHS and the RHS, respectively. The entanglement entropy $\mathcal{S}$
is the von Neumann entropy of the reduced density matrix calculated
for one of the four subsystems, which in the block notation of Figure
<ref> could be denoted by the superscripts $\textrm{B}$,
$\textrm{XA}$, $\textrm{BX}$, and $\textrm{A}$, and which in the
current notation of Figure <ref> are denoted by the
superscripts $\left(1\right)-$, $\left(2\right)-$, $\left(1\right)+$,
and $\left(2\right)+$, respectively. Utilizing a well-known equality
between the von Neumann entropies in a two-partite division (here
denoted by the superscript $a=1,2$) of a system in a pure state $\left|\psi_{k}\right\rangle $
we can skip (on both sides, i.e., for both superscripts “$-$”
and “$+$”) the superscript $a=1,2$, i.e.,

\begin{equation}
\begin{aligned}\mathcal{S}\left(\rho_{k}^{\left(1\right)\mp}\right)= & -\textrm{Tr \ensuremath{\left[\rho_{k}^{\left(1\right)\mp}\ln\rho_{k}^{\left(1\right)\mp}\right]}}\\
= & -\textrm{Tr \ensuremath{\left[\rho_{k}^{\left(2\right)\mp}\ln\rho_{k}^{\left(2\right)\mp}\right]}=\ensuremath{\mathcal{S}\left(\rho_{k}^{\left(2\right)\mp}\right)}\ensuremath{\mathcal{\equiv S}\left(\rho_{k}^{\mp}\right)}}\equiv\mathcal{S}_{k}^{\mp}.
\end{aligned}
\label{eq:entropy}
\end{equation}

By assumption, the unitary evolution (<ref>)
is entirely general and only restricted by a single kinematic condition:
semicausality. It appears that even though we should not expect too
much without any additional dynamic input, some very general results
can be derived. Namely, the global (semi)causal structure of the total
unitary evolution $U$ defined in (<ref>)
(see Figure <ref>) implies bounds on the evolution
of the entanglement entropy $\mathcal{S}$. The evolution of $\mathcal{S}$
itself can be characterized by the steps (here denoted by $\mathcal{S}_{k}^{\mp}$)

\begin{equation}
\mathcal{S}_{1}^{-}\left(=0\right)\rightarrow\mathcal{S}_{1}^{+}\rightarrow\mathcal{S}_{2}^{-}\rightarrow\mathcal{S}_{2}^{+}\rightarrow\cdots\rightarrow\mathcal{S}_{k}^{-}\rightarrow\mathcal{S}_{k}^{+}\rightarrow\cdots\rightarrow\mathcal{S}_{n-1}^{+}\rightarrow\mathcal{S}_{n}^{-}\rightarrow\mathcal{S}_{n}^{+}\left(=0\right),\label{eq:EntropySteps}
\end{equation}

where the presence of zeros ($\mathcal{S}_{1}^{-}=\mathcal{S}_{n}^{+}=0$)
should be obvious by virtue of Figure <ref>.
In turn, the bounds on entropy can be derived from the Araki–Lieb
version of the triangle inequality [Araki and Lieb, 1970]

\begin{equation}
\left|\mathcal{S}^{\left(\textrm{B}\right)}-\mathcal{S}^{\left(\textrm{X}\right)}\right|\leq\mathcal{S}^{\left(\textrm{BX}\right)}\leq\mathcal{S}^{\left(\textrm{B}\right)}+\mathcal{S}^{\left(\textrm{X}\right)},\label{eq:Trianle_ineq}
\end{equation}

here written in the notation of Figures <ref> or <ref>.
In particular, for a one-qubit $\textrm{X}$, the evolution of the
von Neumann entropy is constrained by the inequality (essentially,
a version of the triangle inequality)

\begin{equation}
\left|\Delta\mathcal{S}\right|\equiv\left|\mathcal{S}^{\left(\textrm{BX}\right)}-\mathcal{S}^{\left(\textrm{B}\right)}\right|\leq\ln2.\label{eq:DeltaSIneq}
\end{equation}

Explicitly, (<ref>) is a consequence of the following
three inequalities:

\begin{equation}
\begin{cases}
1^{\textrm{st}} & \mathcal{S}^{\left(\textrm{BX}\right)}-\mathcal{S}^{\left(\textrm{B}\right)}\leq\mathcal{S}^{\left(\textrm{X}\right)};\\
2^{\textrm{nd}} & -\mathcal{S}^{\left(\textrm{X}\right)}\leq\mathcal{S}^{\left(\textrm{BX}\right)}-\mathcal{S}^{\left(\textrm{B}\right)}\\
3^{\textrm{rd}} & \mathcal{S}^{\left(\textrm{X}\right)}\leq\ln2,
\end{cases};\label{eq:eqDeltaSIneq}
\end{equation}

where the $1^{\textrm{st}}$ inequality follows from the RHS of the
triangle inequality (<ref>), the $2^{\textrm{nd}}$
one is a consequence of the LHS of (<ref>), whereas
the $3^{\textrm{rd}}$ inequality is a bound for the entropy for two-dimensional
Hilbert space of the one-qubit subsystem $\textrm{X}$. Alternatively,
instead of the above reasoning in the part referring to the $1^{\textrm{st}}$
and $2^{\textrm{nd}}$ inequality in (<ref>) one
can ignore the RHS of (<ref>) and using arguments
of elementary geometry permute the superscripts in (<ref>)
according to $\left(\begin{array}{ccc}
{\scriptstyle \left(\textrm{B}\right)} & {\scriptstyle \left(\textrm{X}\right)} & {\scriptstyle \left(\textrm{BX}\right)}\\
{\scriptstyle \left(\textrm{BX}\right)} & {\scriptstyle \left(\textrm{B}\right)} & {\scriptstyle \left(\textrm{X}\right)}
\end{array}\right)$.

Translating (<ref>) onto “the indexed notation”
of (<ref>) (see Figure <ref>), for
each “time” step $k$ ($k=1,\ldots,n$), we get the inequality

\begin{equation}
\left|\Delta\mathcal{S}_{k}\right|\equiv\left|\mathcal{S}_{k}^{+}-\mathcal{S}_{k}^{-}\right|\leq\ln2.\label{eq:IndeqDeltaSIneq}
\end{equation}

The description of the evolution of entropy can be simplified a bit
by the fact that entropy is invariant (isometric invariance) with
respect to arbitrary unitary transformations (similarity transformations,
in general), denoted below by the symbol $\mathsf{U}^{(a)}$ and defined
by the operator $U^{\left(a\right)}$ ($a=1,2$), of the corresponding
density matrix, i.e.,

\begin{equation}
\mathsf{U}^{(a)}\left[\mathcal{S}\left(\rho_{k}^{(a)\mp}\right)\right]\equiv\mathcal{S}\left[U^{\left(a\right)}\rho_{k}^{(a)\mp}\left(U^{\left(a\right)}\right)^{-1}\right]=\mathcal{S}\left(\rho_{k}^{(a)\mp}\right).\label{eq:UnitaryInv}
\end{equation}

Therefore, in the series (<ref>), we can make the
following identifications

\begin{equation}
\mathcal{S}_{k}^{+}\equiv\mathcal{S}\left(\rho_{k}^{\left(a\right)+}\right)=\mathsf{U}_{k}^{\left(a\right)}\mathcal{S}\left(\rho_{k}^{(a)+}\right)=\mathcal{S}\left(\rho_{k+1}^{(a)-}\right)\equiv\mathcal{S}_{k+1}^{-}\equiv\mathcal{S}_{k},\qquad a=1,2,\label{eq:UniInvInStep}
\end{equation}

where the first equality follows from <ref> and the
second one follows from the fact that the unitary transformation $\mathsf{U}_{k}^{(a)}$
corresponds to the time evolution (shift in indices) operator $U_{k}^{(a)}$
(entropy is a constant of motion), i.e., in terms of indices alone,
$\mathsf{U}_{k}^{(a)}:\,_{k}^{(a)+}\longmapsto\,_{k+1}^{(a)-}$. In
consequence, (<ref>) allows us to a bit simplify
our notation by skipping the superscripts “$+$” and “$-$”
as indicated in the last identity in (<ref>). By
identifying neighbouring terms according to (<ref>),
the series (<ref>) can be substantially simplified
and rearranged as follows

\begin{equation}
\mathcal{S}_{0}\left(=0\right)\rightarrow\mathcal{S}_{1}\rightarrow\cdots\rightarrow\mathcal{S}_{k}\rightarrow\cdots\rightarrow\mathcal{S}_{n-1}\rightarrow\mathcal{S}_{n}\left(=0\right).\label{eq:EntStepsRe}
\end{equation}

Now also (<ref>) assumes a simpler form

\begin{equation}
\left|\Delta\mathcal{S}_{k}\right|=\left|\mathcal{S}_{k}-\mathcal{S}_{k-1}\right|\leq\ln2,\label{eq:ReIndeqDeltaSIneq}
\end{equation}

where $k=1,\ldots,n$. Iteration of (<ref>)
yields the solution

\begin{equation}
\mathcal{S}_{k}\leq k\ln2+\mathcal{S}_{0},\label{eq:IterEnIneq}
\end{equation}

or (taking into account the initial condition $\mathcal{S}_{0}=0$)

\begin{equation}
\mathcal{S}_{k}\leq k\ln2,\label{eq:SingSIneq}
\end{equation}

or (by virtue of the non-negativity of entropy)

\begin{equation}
0\leq\mathcal{S}_{k}\leq k\ln2\label{eq:DoubleSIneq}
\end{equation}

(see Figure <ref>).

(Color online) Bounds on entropy curves for unitarily evolving systems,
which follow from semicausality. In particular: (a) presents the (beginning
part of the) curve of the quickest possible growth of entanglement
entropy according to the limitation imposed by (<ref>);
(b) takes into account “symmetric” limitations (<ref>)
(here $n=11$); in the limit of a very large $n$, the (shaded) triangle
in (c) presents an allowed area for any entropy curve for semicausal
unitary evolution; here the two upper (dotted) edges of the triangle
form a Page-like curve.

Taking into account that according to (<ref>) $\mathcal{S}_{n}=0$,
utilizing (<ref>) and going backwards in “time”
($k$), we obtain another (twin) inequality

\begin{equation}
\end{equation}

Alternatively, making use of the fact that the time-inversed evolution
$U^{\intercal}$ of the given time evolution $U$ is also time evolution,
we could repeat our analysis for the new (time-inversed) evolution
$U^{\intercal}=U^{-1}$ obtaining (<ref>). The
logical conjunction of (<ref>) and (<ref>)
yields our final (double) inequality

\begin{equation}
0\leq\mathcal{S}_{k}\leq\min\left(k,n-k\right)\ln2,\quad k=0,1,\ldots,n,\label{eq:FinalSIneq}
\end{equation}

which is illustrated (for $n=11$) by Figure <ref>.

§.§.§ The entropy curve for (semi)causal

For huge $n$ according to (<ref>) the lower and the
upper bounds for the entropy $\mathcal{S}_{k}$ can be graphically
represented as the lower (dashed) edge and the two upper (dotted)
edges of the (shaded) triangle in Figure <ref>,
respectively. In other words, for any unitary semicausal evolution,
the entropy curve $\mathcal{S}_{k}$ is supposed to satisfy the following
semicausal entropy bound(s):

* The curve $\mathscr{C}$, a graphic representation of the function
$\mathcal{S}_{k}$, should lie inside the shaded triangle ($\vartriangle$)
depicted in Figure <ref> ($\mathscr{C}\subseteq\vartriangle$).

* Its slope (derivative, in an appropriate continuous limit) should
be bounded by the inequality

\begin{equation}
\left|\dfrac{d\mathcal{S}}{dk}\right|\leq\ln2.\label{eq:EnBound}
\end{equation}

(<ref>) is a differential version of (<ref>).

§ EXAMPLES OF SEMICAUSAL MODELS

To illustrate ideas concerning the (semi)causal structure of unitary
models of BH evolution (evaporation) introduced in the previous section,
we present two large multiparametric families of qubit toy models:
(1) a family of tensor product models defined in Subsection <ref>;
(2) a family of controlled non-product models introduced in Subsection
<ref>. In Subsubsection <ref>
for pedagogical reasons, we start from a very simple three-qubit (3Q)
model, which can serve as a building block for subsequent constructions.
In Subsubsection <ref>, we “tensor exponentiate”
the 3Q model, first yielding its tensor square ($3\textrm{Q}^{2}$),
and next, in Subsubsection <ref>, its
$m$-th tensor power ($3\textrm{Q}^{m}$). In Subsection <ref>
as a less straightforward generalization of the $3\textrm{Q}$ model,
we propose a family of non-product models. In particular, Subsubsection
<ref> introduces a four-qubit
controlled version of the 3Q model, which we can call the controlled
three-qubit (C3Q) model. A (non-product) generalization of the C3Q
model to an arbitrarily large number $n$ of qubits is proposed in
Subsubsection <ref>. Obviously,
by virtue of the construction, all the introduced models are unitary
and semicausal, and in consequence, they satisfy the entropy bounds
of Subsection <ref>.

Since all the proposed models of BH evaporation are formulated in
a rather formal language of qubits and quantum circuits, some explanations
in terms of physics should be given. To begin with, let us observe
that there is a direct correspondence between qubit models and two-state
spin models. In turn, at least in principle, two-state spin models
can be translated (e.g., via the Jordan–Wigner transformation)
onto fermion models. Keeping that observation in mind, for our discussion,
we can loosely, i.e., ignoring the actual Jordan–Wigner
transformation, assume the following convention: the qubit $\left|0\right\rangle $
can be associated to a vacuum state, whereas the qubit $\left|1\right\rangle $
can be associated to a one-particle state.

In physics terms, the picture of the BH evolution depicted in Figure
<ref> looks as follows. The vertical dashed
left line and the right line denote the theoretically earliest moment
and the latest one, respectively, when the BH is allowed to exist—no
causal restrictions (single vertical boxes $U_{0}$ and $U_{n}$)
and then no event horizons. Obviously, in more realistic many-particle
situations, the BH forms much later (because usually, a huge number
of particles is necessary for the actual BH creation) and decays earlier
(it seems reasonable to assume that a too-low number of particles
is unable to support the existence of the BH). Anyway, when a state
(for some $k$ close to $n$, i.e., for the final stages of BH evolution)
quitting the upper box $U_{k-1}^{(1)}$ (potentially belonging to
a BH) is the vacuum state $\left(\ket{\psi_{k}}^{(1)}=\ket{0}^{(1)}\right)$,
i.e., the total state $\ket{\psi_{k}}$ assumes the particular (product)
form $\ket{\psi_{k}}=\ket{0}^{(1)}\otimes\ket{\psi_{k}}^{(2)}$, certainly,
the BH already (at this stage $k$) does not exist (has evaporated)—the
vacuum state $\ket{0}^{(1)}$ cannot support the existence of the
BH. Consequently, the rest of the unitary operators (corresponding
to $i=k,k+1,\ldots,n-1$) ceases to be essential from the point of
view of the BH evolution (evaporation) and thus can be put in the
trivial form $U_{i}=\mathbb{I}$ or (equivalently) can be simply skipped.

§.§ Product models

The product models discussed in this subsection are tensor products
($m$-th tensor powers) of the 3Q model defined in Subsubsection <ref>.
Obviously, $m$-th tensor powers are algebraically unique for a given
$m$, but as will be observed in Subsubsection <ref>,
as far as the (semi)causal structure is concerned, there is a multitude
of different product models for each $m$. From combinatorics, it
follows that there are precisely $2^{-m}\left(2m\right)!$ models
for a given $m$, which in general can yield different entropy curves.

§.§.§ Three-qubit model

The three-qubit (3Q) model first appeared (rather sketchy and in “a
four-qubit incarnation”) in the context of firewalls in [Almheiri et al., 2013]
and then directly and explicitly in [Osuga and Page, 2018]. Another but
akin version of the 3Q model was proposed and discussed in [Broda, 2020]
(see also [Broda, 2021]), whereas its some four-qubit extension
in [Broda, 2021]. In this Subsubsection, the 3Q model appears
for pedagogical reasons as well as an elementary building block for
our further constructions. Our present approach closely follows [Broda, 2020].

style=very thin,fill=blue!5,inner sep=0em,background]
style=thin,rounded corners,color=blue,fill=yellow!5,inner sep=-.1em, background]

$1^{\textrm{st}}\textrm{\ qubit:\ }\ket{x}$    ≫33$U_0$         2$\overset{\textrm{\normalsize$t_1$}}{\md}$    ≫12$U_1$   $\overset{\textrm{\normalsize$t_2$}}{\md}$    1≫22$U_2$    𝚄^†    ≫32$U_3$      $\ket{0}$

$2^{\textrm{nd}}\textrm{\ qubit:\ }\ket{0}$    22               ≫22                    $\ket{0}$

$3^{\textrm{rd}}\textrm{\ qubit:\ }\ket{0}$    𝚄    -1              ≫12                 $\ket{x}$

(Color online) Quantum circuit defining the three-qubit (3Q) model.
The first qubit, in the state $\ket{x}$, is interpreted as a “matter”
qubit entering the BH, whereas the other two (vacuum) qubits are utilized
for Hawking particle production. The Hawking effect parametrized by
$\mathtt{U}$ takes place in the rounded box (see (<ref>))
in the left lower corner of the circuit.

The 3Q model can be concisely described in terms of a quantum circuit,
as presented in Figure <ref>. In the notation of (<ref>)
(see Figure <ref>), the unitary evolution $U$
of the 3Q model assumes the following block structure

\begin{equation}
U=U_{3}U_{2}U_{1}U_{0}\equiv U_{3}\left(U_{2}^{\left(1\right)}\otimes U_{2}^{\left(2\right)}\right)\left(U_{1}^{\left(1\right)}\otimes U_{1}^{\left(2\right)}\right)U_{0}.\label{eq:UforThreeQM}
\end{equation}

Here the first unitary block $U_{0}$ implements the Hawking process
on the second and third qubit, followed by a swap operation on the
first and third one, $U_{1}=U_{3}=\mathbb{I}$ (identity), and $U_{2}=U_{2}^{(1)}\otimes U_{2}^{(2)}$,
where $U_{2}^{(1)}$ can be interpreted as an inverse of the Hawking
process and $U_{2}^{(2)}=\mathbb{I}$.

(Color online) Example entropy curves for the simplest few-qubit product
models. (a) Possible evolution of entanglement entropy for the three-qubit
(3Q) model. In this particular example, in the time interval between
$t_{1}$ and $t_{2}$, the entropy attains its (maximum) value $S=\ln2$
($t_{1}$ and $t_{2}$ are defined in Figure <ref>).
In general, i.e., for arbitrary $\mathtt{U}$, $0\protect\leq S\protect\leq\ln2$.
(b) Possible evolution of entanglement entropy for the tensor square
of the three-qubit ($3\textrm{Q}^{2}$) model. The entropy curve belongs
to “the Page type” (the cases (1)–(4) on the list (<ref>)).
In particular, for $S_{1}=S_{2}=\ln2$ (the extremal case), we would
exactly obtain “the two-qubit Page curve”. $t_{1}$, $t_{2}$,
$t_{3}$, and $t_{4}$ are defined in Figure <ref>.
(c) Another type of possible evolution of entanglement entropy for
the $3\textrm{Q}^{2}$ model. This time the curve is “minimal”
rather than of “the Page type”—the cases (5)–(6)
on the list (<ref>).

In principle, the one-qubit operator $\mathtt{U}$ in Figure <ref>
could be an arbitrary unitary operator, but we can as well assume
that it is orthogonal, given in the form

\begin{equation}
\mathtt{U}=\mathtt{U}_{r_{p}}=\left(\begin{array}{cc}
\cos r_{p} & -\sin r_{p}\\
\sin r_{p} & \cos r_{p}
\end{array}\right),\label{eq:defU}
\end{equation}

because of the correspondence to the Hawking process for fermions,
described by orthogonal (Bogolyubov) transformations. More precisely,
the gate

\begin{equation}
\label{eq:HGate}
\def\gg#1#2{\gategroup[{#1},steps={#2},
style={thin,rounded corners,color=blue,fill=yellow!5,inner sep=-.1em}, background]} %sep=.25em
%### style={dotted,rounded corners,color=green,fill=magenta!5,inner sep=0ex}
\begin{quantikz}
\ls{$\ket{0}$}	&\qw\gg{2}{2}{}		 &\targ{}	&\qw\rstick[wires=2]%
{\ $\cos r_{p}\ket{00}+\sin r_{p}\ket{11}$}	\\
\ls{$\ket{0}$}	&\gt{\mathtt{U}_{r_p}}  &\ctrl{-1}  &\qw
\end{quantikz}
\end{equation}

implements the (one-mode $p$) Hawking process for fermions with mass
$m$ and momentum $p$, where

\begin{equation}
\cos r_{p}=\left[2\cosh\left(\pi\omega_{p}/\kappa\right)\right]^{-1/2}\exp\left(\pi\omega_{p}/2\kappa\right),\quad\omega_{p}=\sqrt{p^{2}+m^{2}},\label{eq:cosrp}
\end{equation}

with $\kappa$ denoting the surface gravity [Mann, 2015].

In the 3Q model, the entanglement entropy $\mathcal{S}$ only changes
at the instants $t_{1}$ (there is a jump of the numeric value of
$\mathcal{S}$ from 0 to some non-negative value $S\leq\ln2$) and
$t_{2}$ (a jump from the value $S$ to $0$ back), as depicted in
Figure <ref> ($t_{1}$ and $t_{2}$ are defined in Figure
<ref>). In terms of the four blocks $U_{0},\ldots,U_{3}$
of Figure <ref>, the evolution of entanglement entropy
proceeds according to the pattern

\begin{equation}
\underset{0}{0}\rightarrow\underset{1}{S}\rightarrow\underset{2}{0}\rightarrow\underset{3}{0},\label{eq:EEntropy43Q}
\end{equation}

where the underscripts number the unitary blocks. The evolution of
entropy can be called extremal, if the entropy attains its maximal
allowed value, i.e., for this model $S=\ln2$. According to (<ref>)
this happens in the limit of $\kappa\rightarrow\infty$ or $\omega_{p}\left(=\sqrt{p^{2}+m^{2}}\right)\rightarrow0$
(maximal entanglement, see (<ref>)).

§.§.§ Tensor square of the three-qubit

As an intermediate stage between the 3Q model and its arbitrary $m$-th
tensor power $\left(3\textrm{Q}^{m}\right)$, we will consider the
$3\textrm{Q}^{2}$ model, which is a tensor square of the 3Q model.
The $3\textrm{Q}^{2}$ model operates on six qubits, and its unitary
evolution $U$ can be semicausally decomposed into $1+6\left(\textrm{number of qubits}\right)=7$

\begin{equation}
\begin{aligned}U= & U_{6}U_{5}U_{4}U_{3}U_{2}U_{1}U_{0}\\
\equiv & U_{6}\left(U_{5}^{\left(1\right)}\otimes U_{n-1}^{\left(2\right)}\right)\ldots\left(U_{2}^{\left(1\right)}\otimes U_{2}^{\left(2\right)}\right)\left(U_{1}^{\left(1\right)}\otimes U_{1}^{\left(2\right)}\right)U_{0}.
\end{aligned}
\label{eq:UforTSThreeQM}
\end{equation}

One of the versions of the model (a preferred one from our perspective)
is explicitly presented in Figure <ref>.

style=very thin,fill=blue!5,inner sep=.15em, outer sep=0em,background]

x_1   ≫63$U_0$      $\overset{\textrm{\normalsize$t_1$}}{\md}$    ≫12$U_1$
≫31$U_3$$\overset{\textrm{\normalsize$t_4$}}{\md}$    1≫41$U_4$    ≫51$U_5$    ≫61$U_6$   0

x_2            ≫52               2       0

0                  ≫41                0

0                     ≫31             0

0   1                     ≫21          x_1

0         2                   ≫11       x_2

(Color online) Quantum circuit defining one of the 6 variants of the
$3\textrm{Q}^{2}$ model (the first case on the list of possible semicausal
embeddings (<ref>)). The first two qubits, in the states
$\ket{x_{1}}$ and $\ket{x_{2}}$, are interpreted as “matter”
qubits entering the BH, whereas the other four (vacuum) qubits are
utilized for Hawking particle production.

As already mentioned at the beginning of Subsection <ref>,
even though from a purely algebraic point of view the tensor square
(and more generally, the tensor power) is unique, the corresponding
models are not, because of various possibilities of embeddings in
the surrounding semicausal frame (Figure <ref>).
That follows from the fact that particular qubits can enter successive
upper subboxes $\left(U_{1}^{(1)},U_{2}^{(1)},\ldots,U_{n-1}^{(1)},U_{n}\right)$
in various orders, giving rise to (in principle) different entropy
curves. For the $3\textrm{Q}^{2}$ model, we have six possibilities,
which can be identified by the following entropy jumps:

\begin{equation}
\begin{alignedat}{1}(1)\qquad & \underset{0}{0}\rightarrow\underset{1}{S_{1}}\rightarrow\underset{2}{S_{1}+S_{2}}\rightarrow\underset{3}{S_{2}}\rightarrow\underset{4}{0}\rightarrow\underset{5}{0}\rightarrow\underset{6}{0};\\
(2)\qquad & \underset{}{0}\rightarrow\underset{}{S_{1}}\rightarrow\underset{}{S_{1}+S_{2}}\rightarrow\underset{}{S_{1}}\rightarrow\underset{}{0};\\
(3)\qquad & \underset{}{0}\rightarrow\underset{}{S_{2}}\rightarrow\underset{}{S_{1}+S_{2}}\rightarrow\underset{}{S_{1}}\rightarrow\underset{}{0};\\
(4)\qquad & \underset{}{0}\rightarrow\underset{}{S_{2}}\rightarrow\underset{}{S_{1}+S_{2}}\rightarrow\underset{}{S_{2}}\rightarrow\underset{}{0};\\
(5)\qquad & \underset{}{0}\rightarrow\underset{}{S_{1}}\rightarrow\underset{}{0}\rightarrow\underset{}{S_{2}}\rightarrow\underset{}{0};\\
(6)\qquad & \underset{}{0}\rightarrow\underset{}{S_{2}}\rightarrow\underset{}{0}\rightarrow\underset{}{S_{1}}\rightarrow\underset{}{0.}
\end{alignedat}
\label{eq:SixPoss}
\end{equation}

The first case (1) on the list (<ref>) (with explicitly
numbered blocks $U_{0},\ldots,U_{6}$ by the underscripts $0,1.\ldots,6$)
corresponds to the model presented in Figure <ref>,
and its entropy jumps are qualitatively illustrated in Figure <ref>;
the second case (2) is some variation of (1); the cases (3) and (4)
are “qubit transpositions” of (1) and (2), respectively; whereas
the cases (5) and (6) are “causally decoupled” models. The cases
(1)–(4) belong to “the Page type”, i.e., the shape
of their entropy curves resembles the shape of the curve depicted
in Figure <ref>; in particular, for $S_{1,2}=\ln2$ (the
extremal case), we exactly obtain “the two-qubit Page curve”.
The cases (5) and (6) schematically presented in Figure <ref>
and Figure <ref>, respectively, belong to the “minimal”
type, i.e., the shape of their entropy curves resembles the shape
of the curve depicted in Figure <ref>.

One should also note that $S_{1}$ and $S_{2}$ $\left(0\leq S_{1,2}\leq\ln2\right)$
are contributions to the entanglement entropy $\mathcal{S}$ coming
from the Hawking effect rather than from the “matter” qubits $\ket{x_{1}}$
and $\ket{x_{2}}$.

≫[2,steps=1,style=dotted,rounded corners,color=green,fill=magenta!5,inner sep=0ex,background]

1^st triple:

|x_1⟩,|0⟩,|0⟩        Figure <ref>

2^nd triple:

|x_2⟩,|0⟩,|0⟩       [alternate]
Figure <ref>

[]≫[2,steps=1,style=dotted,rounded corners,color=green,fill=magenta!5,inner sep=0ex,background]

Figure <ref>

𝚄=𝚄_1 [alternate]≫

Figure <ref>

𝚄=𝚄_2 [alternate]

(Color online) Schematic presentation of quantum circuits defining
the last two cases on the list (<ref>) of all possible
semicausal embeddings of the $3\textrm{Q}^{2}$ model. The cases (5)
and (6), depicted in (a) and (b),  respectively, belong to the “minimal”
(or “decoupled”) type, i.e., the shape of their entropy curves
is similar to the shape of the curve depicted in Figure <ref>.

§.§.§ Arbitrary tensor power of the
three-qubit model

From the point of view of entropy curves (or entropy evolution), various
$3\textrm{Q}^{m}$ models with the fixed tensor power $m$ can be
characterized by various series of entropy contributions $\pm S_{1},\pm S_{2},\ldots,\pm S_{m}$,
where the numeric values $S_{i}$ ($0\leq S_{i}\leq\ln2$, $i=1,2,\ldots,m$)
of the entanglement entropy $\mathcal{S}$ depend on the explicit
form of the corresponding $U(2)$ matrices $\mathtt{U}_{i}$. For
example, the subset of all $3\textrm{Q}^{m}$ models corresponding
to entropy jumps (analogous to the cases (1)–(4) on the
list (<ref>))

\begin{equation}
\begin{alignedat}{1}0\rightarrow & S_{\sigma_{+}\left(1\right)}\rightarrow S_{\sigma_{+}\left(1\right)}+S_{\sigma_{+}\left(2\right)}\rightarrow\cdots\rightarrow S_{\sigma_{+}\left(1\right)}+S_{\sigma_{+}\left(2\right)}+\cdots+S_{\sigma_{+}\left(m\right)}\\
\rightarrow & S_{\sigma_{+}\left(1\right)}+S_{\sigma_{+}\left(2\right)}+\cdots+S_{\sigma_{+}\left(m\right)}-S_{\sigma_{-}\left(1\right)}\\
\rightarrow & S_{\sigma_{+}\left(1\right)}+S_{\sigma_{+}\left(2\right)}+\cdots+S_{\sigma_{+}\left(m\right)}-S_{\sigma_{-}\left(1\right)}-S_{\sigma_{-}\left(2\right)}\rightarrow\cdots\rightarrow0,
\end{alignedat}
\label{eq:PagePoss}
\end{equation}

where $\sigma_{+}$ and $\sigma_{-}$ are two independent permutations
of $m$ elements, is of “the Page type”, i.e., for $S_{i}\approx\ln2$
the shape of the entropy curves for these models is similar to the
shape of the curve depicted in Figure <ref> or
Figure <ref> for $m\gg1$ (two upper edges of
the triangle). Instead, the subset of $3\textrm{Q}^{m}$ models corresponding
to entropy jumps (analogous to the cases (5)–(6) in (<ref>))

\begin{equation}
0\rightarrow S_{\sigma\left(1\right)}\rightarrow0\rightarrow S_{\sigma\left(2\right)}\rightarrow0\rightarrow\cdots\rightarrow0\rightarrow S_{\sigma\left(m\right)}\rightarrow0,\label{eq:ZeroPoss}
\end{equation}

where $\sigma$ is a (single) permutation of $m$ elements, belongs
to the “minimal” type, i.e., for $m\gg1$ their entropy curves
are very close to the lower (dashed) edge of the (shaded) triangle
in Figure <ref>. Let us also note that the number
of models belonging to the Page type (<ref>) is equal
to $\left(m!\right)^{2}$, and for $m\gg1$ this number is much larger
than the number of “minimal” models (<ref>), which
is equal to $m!$.

§.§ Non-product models

A lot of (product) models introduced in the previous Subsection <ref>
can yield entropy curves arbitrarily close to the Page curve, but
from the point of view of physics, the models can seem to be almost
trivial—interactions are bounded to mutually isolated and
almost identical (i.e., possibly differing only by one-qubit operators
$\mathtt{U}_{i}$ or equivalently by the parameters $r_{p_{i}}$ in
(<ref>)) three-qubit sectors. The Hawking process in these
models is determined by the fixed family of matrices $\left\{ \mathtt{U}_{i}\right\} _{i=1}^{m}$.
Therefore, the whole system looks like an open one defined
by externally prescribed $r_{p_{i}}$, $i=1,2,\ldots,m$. To remove
these limitations, in this subsection we introduce a family of controlled
non-product models, which are controlled with a relatively small set
of $l$ control qubits $\left\{ \left|y_{i}\right\rangle \right\} _{i=1}^{l}$.
The control qubits should be interpreted as the part of all “matter”
qubits that gradually enter the BH, but their modest number is unable
to support the existence of the BH in the final stages of its evolution.
Thus, the total number of qubits in these models amounts to $n=l+3m$
($l$ control qubits and $3m$ “controlled” qubits, including
$m$ “matter” qubits as well as $2m$ Hawking ones), where for
large $m$ (the case of BHs), we should assume that $l\ll m$. The
simplest case, $l=1$ and $m=1$, i.e., the controlled three-qubit
(C3Q) model, is discussed in the following Subsubsection <ref>,
whereas the case of arbitrary $l$ and $m$, i.e., the $l$-controlled
$3m$-qubit ($l\textrm{C}3m\textrm{Q}$) model, is introduced in Subsubsection

§.§.§ The controlled three-qubit

The controlled three-qubit (C3Q) model operates on four qubits (see
Figure <ref>):

style=very thin,fill=blue!2,inner sep=.25em,background]

$\ket{y}$   3≫41$U_0$    ≫12$U_1$       ≫22$U_2$       ≫32$U_3$    1    ≫42$U_4$       $\ket{y}$

$\ket{x}$      ≫32    2          1
𝚄^†          $\ket{0}$

$\ket{0}$             ≫22                    $\ket{0}$

$\ket{0}$   𝚄    -1              ≫12             $\ket{x}$

Quantum circuit defining the controlled three-qubit (C3Q) model. The
first two qubits, in the states $\ket{y}$ (a control qubit) and $\ket{x}$,
are “matter” qubits, whereas the other two vacuum states are utilized
for Hawking particle production. In the first unitary block ($U_{0}$)
the qubit $\ket{y}$ controls the Hawking effect by the CNOT-$\mathtt{U}$
gate determined by the $U(2)$ matrix $\mathtt{U}$.

two “matter” qubits in arbitrary states $\ket{y}$ (a control
qubit) and $\ket{x}$ as well as two “Hawking qubits”, initially
both in the vacuum state $\ket{0}$. The model is parameterized by
a single one-qubit $U(2)$ matrix $\mathtt{U}$, possibly limited
to the orthogonal (Bogolyubov–Hawking) form (<ref>).
The qubit $\ket{y}$ controls (according to the NOT-rule) the Hawking
effect by the CNOT-$\mathtt{U}$ gate in the block $U_{0}$ of Figure

The evolution of the entanglement entropy $\mathcal{S}$ proceeds
in consecutive blocks $U_{0},\ldots,U_{4}$, according to the pattern
(compare to (<ref>))

\begin{equation}
\underset{0}{0}\rightarrow\underset{1}{S'}\rightarrow\underset{2}{S''}\rightarrow\underset{3}{0}\rightarrow\underset{4}{0},\label{eq:EEntropy4C3Q}
\end{equation}

where the numeric values of $\mathcal{S}$, i.e., $S'$ and $S''$,
with $0\leq S',S''\leq\ln2$, are determined by the qubit $\ket{y}$
and the explicit form of the matrix $\mathtt{U}$. One should remark
that the contributions $S'$ and $S''$ are coming from the Hawking
particles as well as from the control “matter” qubit $\ket{y}$.

Let us note that some immediate and natural generalization of the
C3Q model is possible by the replacement of the CNOT-$\mathtt{U}$
gate (defined with a single $U(2)$ matrix $\mathtt{U}$) in the block
$U_{0}$ of Figure <ref> by a gate determined with a
pair of $U(2)$ matrices $\mathtt{U}_{1}$ and $\mathtt{U}_{2}$,
i.e., for the whole two-qubit unitary operator, we have the generalization

\begin{equation}
\begin{pmatrix}\mathbb{I} & 0\\
0 & \mathtt{U}
\end{pmatrix}\longrightarrow\mathcal{U}\equiv\begin{pmatrix}\mathtt{U}_{1} & 0\\
0 & \mathtt{U}_{2}
\end{pmatrix},\label{eq:gen.contr.u-gate}
\end{equation}

where the matrices $\mathtt{U}_{1}$ and $\mathtt{U}_{2}$ correspond
to the control qubits $\ket{0}$ and $\ket{1}$, respectively.

§.§.§ Controlled multiqubit

The $l$-controlled $3m$-qubit ($l\textrm{C}3m\textrm{Q}$) model
is a direct generalization of the C3Q model of the previous subsubsection,
where instead of a single control qubit $\ket{y}$ we have a (small)
set of $l$ control qubits $\ket{y_{1}},\ket{y_{2}},\ldots,\ket{y_{l}},$
and instead of a single triple of the controlled qubits (initially
in the states $\ket{x}$, $\ket{0}$, $\ket{0}$) we have a set of
$m$ triples (initially in the states $\ket{x_{1}},\ket{0},\ket{0},\ldots,\ket{x_{m}},\ket{0},\ket{0}$).
Then, the total number of qubits amounts to $n=l+3m$. Moreover, for
each triple $i$ ($i=1,2,\ldots,m$) there is defined a collection
of the $2^{l}$ $U(2)$ matrices $\left\{ \mathtt{U}_{j}\left(i\right)\right\} _{j=1}^{2^{l}}$
that determine an $l+1$-qubit unitary matrix $\mathcal{U}_{i}$ of
dimension $2^{l+1}\times2^{l+1}$—a generalization of (<ref>)
to $l$ control qubits. Explicitly,

\begin{equation}
\mathcal{U}_{i}=\left(\begin{array}{cc}
\begin{array}{cc}
\mathtt{U}_{1}\left(i\right)\\
& \mathtt{U}_{2}\left(i\right)
\end{array} & \raisebox{-1ex}{\text{\Huge0}}\\
\kern-1em \raisebox{-1ex}{\text{\Huge0}} & \begin{array}{cc}
\ddots\\
& \mathtt{U}_{2^{l}}\left(i\right)
\end{array}
\end{array}\right),\label{eq:gen.gen.contr.u-gate}
\end{equation}

where the $U(2)$ matrices $\mathtt{U}_{j}\left(i\right)$ operate
in the $i$-th qubit triple and binary presentation of the number
$j$ corresponds to the series of “control” digits given by $y_{1}\cdots y_{l}$.

The general structure of the $l\textrm{C}3m\textrm{Q}$ model is schematically
depicted in block form in Figure <ref>.

$l$ control qubits:

[]1-.5em.5ex$l$     …    3

1^st triple:

{ 𝚄_j(1)} _j=1^2^l
(see Figure <ref>)

[alternate]      … [alternate]    [alternate]         [alternate]

m^th triple:

… [alternate]     { 𝚄_j(m)} _j=1^2^l
(see Figure <ref>)

[alternate]     [alternate]

Schematic presentation (in block form) of the quantum circuit for
the $l$-controlled $3m$-qubit ($l\textrm{C}3m\textrm{Q}$) model
without explicitly imposed semicausal structure. The model operates
on $n\left(=l+3m\right)$ qubits, where for each qubit triple $\left\{ \ket{x_{i}},\ket{0},\ket{0}\right\} $
($i=1,2,\ldots,m$), there is an associated collection of the $2^{l}$
$U(2)$ matrices $\left\{ \mathtt{U}_{j}\left(i\right)\right\} _{j=1}^{2^{l}}$
that define an $l+1$-qubit unitary operator $\mathcal{U}_{i}$ of
the form (<ref>). The matrices $\mathtt{U}_{1}\left(i\right)$,
$\mathtt{U}_{2}\left(i\right)$, $\ldots$ , $\mathtt{U}_{2^{l}}\left(i\right)$
act on a single qubit in the box corresponding to the given $i$-th

The whole family of the models is parametrized by $m$ $2^{l}$-element
collections of $U(2)$ matrices $\left\{ \mathtt{U}_{j}\left(i\right)\right\} _{j=1}^{2^{l}}$,
$i=1,\ldots,m$, and by a multitude of various semicausal structures.
By virtue of the construction, entropy curves are still within the
semicausal entropy bound given in (<ref>) and described
in Subsubsection (<ref>), but they in a
more complicated way depend on the matrices $\mathtt{U}_{j}\left(i\right)$
as well as on the initial values of the control “matter” qubits

§ SUMMARY AND FINAL REMARKS

Utilizing the universal form of semicausal evolution (<ref>)
(see Figure <ref>), we have established a general
structure of the quantum circuit presenting unitary evolution and
respecting causality imposed by the BH event horizon (see Figures
<ref> and <ref>). Moreover,
making use of several standard properties of the von Neumann entropy
(the triangle inequality, isometric invariance, and others), we have
determined principles governing the evolution of the corresponding
BH entanglement entropy (see (<ref>) and (<ref>)),
and next visualized them in Figure <ref> as well
as additionally summarized in Subsubsection <ref>.

In the second part of the work (Section <ref>),
for illustration purposes, two multiparametric families of unitary
semicausal qubit toy models have been considered: tensor product models
($3\textrm{Q}^{m}$ models) and controlled non-product models ($l\textrm{C}3m\textrm{Q}$
models). By virtue of the construction of the $3\textrm{Q}$ model
(Figure <ref>), the role played by individual qubits
and consequently the direction of the flow of information are rather
clear. In particular, we can state that “external information”
carried by the $1^{\textrm{st}}$ qubit (a “matter” qubit in the
state $\ket{x}$) is exchanged by the swap gate with that in the outgoing
Hawking qubit (the $3^{\textrm{rd}}$ qubit in the box $U_{0}$).
Therefore, evidently, in the $3\textrm{Q}^{m}$ models, Hawking radiation
produces BH entropy. In turn, at the latter stages of BH evolution,
the accumulated BH entropy is gradually being cancelled by the partner
qubits of the Hawking pairs (see the $2^{\textrm{nd}}$ qubit in box
$U_{2}$). In other words, “truly external” (carried by the $1^{\textrm{st}}$
qubit) information does not cross the BH horizon and is (therefore)
not lost. In Subsection <ref>, to “dynamically
close” the system, we have introduced the controlled non-product
models with a limited number of $l$ (external) control qubits, which
serve to parameterize the Hawking effect. Since the control qubits
are supposed to enter the BH as well as to be later finally released,
in these models, we are forced to (additionally, but hopefully reasonably)
assume that a very small number of particles is unable to support
the existence of the BH, and therefore the BH has to decay and free
the rest of particles (exactly, the control matter particles). Presumably,
such a physics scenario could be conceivable for the total mass of
all $l$ control particles less than the Planck mass, i.e., $M_{l\textrm{\thinspace control}}<M_{\textrm{Planck}}$.
Consequently, in the controlled non-product models, a small contribution
to the BH entropy also comes from the control matter particles.

[Hawking, 1976]
S. W. Hawking.
Breakdown of predictability in gravitational collapse.
Physical Review D, 14:0 2460–2473, November 1976.
URL <https://link.aps.org/doi/10.1103/PhysRevD.14.2460>.

[Almheiri et al., 2019]
Ahmed Almheiri, Netta Engelhardt, Donald Marolf, and Henry Maxfield.
The entropy of bulk quantum fields and the entanglement wedge of an
evaporating black hole.
Journal of High Energy Physics, 20190 (12):0
63, December 2019.

[Penington, 2020]
Geoffrey Penington.
Entanglement wedge reconstruction and the information paradox.
Journal of High Energy Physics, 20200 (9):0 2,
September 2020.

[Page, 1993]
Don N. Page.
Information in black hole radiation.
Physical Review Letters, 71:0 3743–3746, December
URL <https://link.aps.org/doi/10.1103/PhysRevLett.71.3743>.

[Avery, 2013]
Steven G. Avery.
Qubit Models of Black Hole Evaporation.
Journal of High Energy Physics, 01:0 176, 2013.

[Osuga and Page, 2018]
Kento Osuga and Don N. Page.
Qubit transport model for unitary black hole evaporation without
Physical Review D, 97:0 066023, March 2018.
URL <https://link.aps.org/doi/10.1103/PhysRevD.97.066023>.

[Broda, 2020]
Bogusław Broda.
Unitary toy qubit transport model for black hole evaporation.
The European Physical Journal C, 800 (5):0
418, May 2020.

[Broda, 2021]
Bogusław Broda.
Fermionic model of unitary transport of qubits from a black hole.
Physical Review D, 103:0 025022, January 2021.
URL <https://link.aps.org/doi/10.1103/PhysRevD.103.025022>.

[Broda, 2021]
Bogusław Broda.
Causal unitary qubit model of black hole evaporation.
Physics Letters B, 820:0 136564, September 2021.

[Beckman et al., 2001]
David Beckman, Daniel Gottesman, M. A. Nielsen, and John Preskill.
Causal and localizable quantum operations.
Physical Review A - Atomic, Molecular, and Optical Physics,
64, 2001.
ISSN 10941622.

[Eggeling et al., 2002]
T. Eggeling, D. Schlingemann, and R. F. Werner.
Semicausal operations are semilocalizable.
Europhysics Letters, 57, 2002.
ISSN 02955075.

[Piani et al., 2006]
M. Piani, M. Horodecki, P. Horodecki, and R. Horodecki.
Properties of quantum nonsignaling boxes.
Physical Review A - Atomic, Molecular, and Optical Physics,
74, 2006.
ISSN 10502947.

[Schumacher and Westmoreland, 2012]
Benjamin Schumacher and Michael D. Westmoreland.
Isolation and Information Flow in Quantum Dynamics.
Foundations of Physics, 42, 2012.
ISSN 00159018.

[Schumacher and Westmoreland, 2005]
Benjamin Schumacher and Michael D. Westmoreland.
Locality and information transfer in quantum operations.
Quantum Information Processing, 4, 2005.
ISSN 15700755.

[Braunstein and Pirandola, 2018]
Samuel L. Braunstein and Stefano Pirandola.
Quantum information versus black hole physics: Deep firewalls from
narrow assumptions.
Philosophical Transactions of the Royal Society A:
Mathematical, Physical and Engineering Sciences, 376, 2018.
ISSN 1364503X.

[Araki and Lieb, 1970]
Huzihiro Araki and Elliott H Lieb.
Entropy inequalities.
Communications in Mathematical Physics, 180
(2):0 160–170, 1970.

[Almheiri et al., 2013]
Ahmed Almheiri, Donald Marolf, Joseph Polchinski, Douglas Stanford, and James
An Apologia for Firewalls.
Journal of High Energy Physics, 09:0 018, 2013.

[Mann, 2015]
Robert B. Mann.
Black Holes: Thermodynamics, Information, and Firewalls.
SpringerBriefs in Physics. Springer, 2015.
ISBN 978-3-319-14495-5, 978-3-319-14496-2.

This work was supported by the University of Łódź.