##  Appendix A Numerical method

For the numerical simulations of section 2, we use the formulation introduced
in [26], which readily adapts to the viscoelastic setting. The formulation is
based on a combination of works by Moreau, et al. [23] and Maxian, et al.
[22]. For convenience, we recall the original formulation (1)-(4) of the
viscoelastic resistive force theory equations:

$\displaystyle\frac{\partial\bm{X}}{\partial t}(s,t)$
$\displaystyle=-(1+\mu)\big{(}{\bf I}+\gamma\bm{X}_{s}\bm{X}_{s}^{\rm
T}\big{)}\big{(}\bm{X}_{sss}-\tau\bm{X}_{s}-(\kappa_{0})_{s}\bm{e}_{\rm
n}-\frac{\mu}{1+\mu}\xi_{s}\bm{e}_{\rm n}\big{)}_{s}$ (95)
$\displaystyle\delta\dot{\xi}$ $\displaystyle=-\xi+\kappa-\kappa_{0}$
$\displaystyle\left\lvert\bm{X}_{s}\right\rvert^{2}$ $\displaystyle=1$
$\displaystyle(\bm{X}_{ss}-\kappa_{0}\bm{e}_{\rm n})\big{|}_{s=0,1}$
$\displaystyle=0\,,\quad(\bm{X}_{sss}-\tau\bm{X}_{s}-(\kappa_{0})_{s}\bm{e}_{\rm
n})\big{|}_{s=0,1}=0\,,\quad\xi\big{|}_{s=0,1}=\xi_{s}\big{|}_{s=0,1}=0\,.$

The formulation used in numerical simulations will be derived from (95). We
begin by parameterizing the filament using the tangent angle description (11).
In particular, we write

$\bm{X}(s,t)=\bm{X}_{0}(t)+\int_{0}^{s}\bm{e}_{\rm
t}(s^{\prime},t)ds^{\prime}\,,\quad\bm{e}_{\rm t}=\begin{pmatrix}\cos\theta\\\
\sin\theta\end{pmatrix}\,,\quad\bm{e}_{\rm n}=\begin{pmatrix}-\sin\theta\\\
\cos\theta\end{pmatrix}$ (96)

As in [26, 23, 22], we exploit that, due to the inextensibility constraint,
only the normal components of the hydrodynamic force along the filament
actually contribute to the fiber motion. In particular, using the
parameterization (96), we may rewrite the first three equations of (95) as a
closed system:

$\displaystyle\dot{\bm{X}}_{0}+\int_{0}^{s}\dot{\bm{e}}_{\rm
t}(s^{\prime})\,ds^{\prime}$ $\displaystyle=-(1+\mu)({\bf I}+\gamma\bm{e}_{\rm
t}\bm{e}_{\rm t}^{\rm T})\bm{h}(s)$ (97) $\displaystyle\big{(}{\bf
I}-\bm{e}_{\rm t}(s)\bm{e}_{\rm t}(s)^{\rm
T}\big{)}\int_{0}^{s}\bm{h}(s^{\prime})ds^{\prime}$ $\displaystyle=\big{(}{\bf
I}-\bm{e}_{\rm t}(s)\bm{e}_{\rm t}(s)^{\rm
T}\big{)}\bigg{(}\bm{X}_{sss}-(\kappa_{0})_{s}\bm{e}_{\rm
n}-\frac{\mu}{1+\mu}\xi_{s}\bm{e}_{\rm n}\bigg{)}$ (98)
$\displaystyle\delta\dot{\xi}$ $\displaystyle=-\xi+\kappa-\kappa_{0}\,.$ (99)

Note that in (98), by projecting away from the tangential direction along the
filament, we have eliminated the need to solve for the unknown fiber tension.
Instead, inextensibility is enforced directly via the parameterization (96).

Solving (97) directly for $\bm{h}$ and inserting this expression in (98), we
obtain the system

$\displaystyle\bm{e}_{\rm n}(s,t)\cdot\int_{0}^{s}({\bf
I}-\frac{\gamma}{1+\gamma}\bm{e}_{\rm t}\bm{e}_{\rm t}^{\rm
T})\bigg{(}\dot{\bm{X}}_{0}+\int_{0}^{s^{\prime}}\dot{\bm{e}}_{\rm
t}(\overline{s})\,d\overline{s}\bigg{)}ds^{\prime}$
$\displaystyle=-(1+\mu)\theta_{ss}+(1+\mu)(\kappa_{0})_{s}+\mu\xi_{s}$ (100)
$\displaystyle\delta\dot{\xi}$ $\displaystyle=-\xi+\theta_{s}-\kappa_{0}$
(101)

for unknowns $\bm{X}_{0}(t)$, $\theta(s,t)$, and $\xi(s,t)$. Equations (100)
and (101) serve as the basis for our numerical method. The boundary conditions
$(\theta_{s}-\kappa_{0})\big{|}_{s=0,1}=0$ are enforced directly in the
discretization of $\theta_{ss}$ on the right hand side of (100), while
$\xi\big{|}_{s=0,1}=\xi_{s}\big{|}_{s=0,1}=0$ is enforced in the
discretization of $\xi$ in (101). To enforce the boundary condition
$(-\theta_{ss}+(\kappa_{0})_{s})\big{|}_{s=1}=0$, we will also need to require

$\int_{0}^{1}({\bf I}-\frac{\gamma}{1+\gamma}\bm{e}_{\rm t}\bm{e}_{\rm t}^{\rm
T})\bigg{(}\dot{\bm{X}}_{0}+\int_{0}^{s}\dot{\bm{e}}_{\rm
t}(s^{\prime})\,ds^{\prime}\bigg{)}\,ds=0.$ (102)

The analogous condition at $s=0$ is then satisfied automatically via the
formulation (100).

We discretize the arclength coordinate $s\in[0,1]$ into $N+1$ equally spaced
points $s_{i}$, $i=0,\dots,N$ and define $\bm{X}_{i}=\bm{X}(s_{i})$. We
consider the fiber as $N$ straight segments between each $\bm{X}_{i}$, and
define $\theta_{i}$, $i=1,\dots,N$, to be the angle between segment $i$ and
the $x$-axis.

The evolution equation (100) is enforced at the midpoint of each segment
$\bm{X}_{i-\frac{1}{2}}:=\frac{\bm{X}_{i-1}+\bm{X}_{i}}{2}$, $i=1,\dots,N$. In
particular, we parameterize the evolution $\dot{\bm{X}}_{i-\frac{1}{2}}$ of
each fiber segment as

$\displaystyle\dot{\bm{X}}_{i-\frac{1}{2}}=\begin{pmatrix}\dot{x}_{0}\\\
\dot{y}_{0}\end{pmatrix}+\frac{1}{2N}\begin{pmatrix}-\sin\theta_{i}\\\
\cos\theta_{i}\end{pmatrix}\dot{\theta}_{i}+\frac{1}{N}\sum_{k=1}^{i}\begin{pmatrix}-\sin\theta_{k}\\\
\cos\theta_{k}\end{pmatrix}\dot{\theta}_{k}\,,\quad i=1,\dots,N\,.$

We also define $\kappa_{0,i}=\kappa_{0}(s_{i-\frac{1}{2}})$ where
$s_{i-\frac{1}{2}}=\frac{s_{i-1}+s_{i}}{2}$, $i=1,\dots,N$, and we define
$\xi_{i}$ similarly.

For the middle segments $j=2,\dots,N-1$, we obtain $2(N-2)$ equations from the
discretization of (100) and (101):

$\displaystyle\frac{1}{N}\begin{pmatrix}-\sin\theta_{j}\\\
\cos\theta_{j}\end{pmatrix}\cdot\sum_{i=1}^{j}\bm{M}_{\rm
RFT}(\theta_{i})\dot{\bm{X}}_{i-\frac{1}{2}}$
$\displaystyle=-N^{2}(1+\mu)(\theta_{j-1}-2\theta_{j}+\theta_{j+1})+(1+\mu)(\kappa_{0})_{s,j}$
(103) $\displaystyle\qquad+\mu\frac{N}{2}(\xi_{j+1}-\xi_{j-1})\,,\quad
j=2,\dots,N-1$ $\displaystyle\delta\dot{\xi}_{j}$
$\displaystyle=-\xi_{j}+2N(\theta_{j+1}-\theta_{j-1})-\kappa_{0,j}\,,\quad
j=2,\dots,N-1\,.$ (104)

Here the $2N\times 2N$ matrix $\bm{M}_{\rm RFT}(\theta_{i})$ is given by

$\displaystyle\bm{M}_{\rm
RFT}(\theta_{i})=\begin{pmatrix}1-\frac{\gamma}{1+\gamma}\cos^{2}\theta_{i}&-\frac{\gamma}{1+\gamma}\cos\theta_{i}\sin\theta_{i}\\\
-\frac{\gamma}{1+\gamma}\cos\theta_{i}\sin\theta_{i}&1-\frac{\gamma}{1+\gamma}\sin^{2}\theta_{i}\end{pmatrix}\,.$

At the fiber endpoints, we set $\xi_{1}=\xi_{N}=0$ and enforce the boundary
conditions $(-\theta_{s}+\kappa_{0})_{s}\big{|}_{s=0,1}=0$ via the following
two equations:

$\displaystyle\frac{1}{N}\begin{pmatrix}-\sin\theta_{1}\\\
\cos\theta_{1}\end{pmatrix}\cdot\bm{M}_{\rm
RFT}(\theta_{1})\dot{\bm{X}}_{\frac{1}{2}}$
$\displaystyle=-N^{2}(2\theta_{2}-2\theta_{1})+2N\kappa_{0,1}$ (105)
$\displaystyle\frac{1}{N}\sum_{i=1}^{N-1}\bm{M}_{\rm
RFT}(\theta_{i})\dot{\bm{X}}_{i-\frac{1}{2}}$ $\displaystyle=0\,,$ (106)

To enforce $(\theta_{s}-\kappa_{0})\big{|}_{s=0,1}=0$, we discretize
$\theta_{ss}$ near the fiber endpoints as

$\theta_{ss}\big{|}_{s=0}\approx
N^{2}(2\theta_{2}-2\theta_{1})-2N\kappa_{0,1}\,,\qquad\theta_{ss}\big{|}_{s=1}\approx
N^{2}(2\theta_{N-1}-2\theta_{N})+2N\kappa_{0,N}\,.$ (107)

At the $s=1$ endpoint, equations (100) and (102) coincide to give the boundary
condition $\theta_{ss}\big{|}_{s=1}=(\kappa_{0})_{s}\big{|}_{s=1}$, which,
using (107), becomes an equation for $\theta_{N}$:

$\theta_{N}=\theta_{N-1}+\frac{1}{N}\kappa_{0,N}-\frac{1}{2N^{2}}(\kappa_{0})_{s,N}\,.$
(108)

Counting equations, we have $2(N-2)$ equations from (103) and (104), 1
equation from (105), 2 equations from (106), and 1 equation from (108) for a
total of $2N$ equations. These uniquely determine the $2N$ unknowns $x_{0}$,
$y_{0}$, $\theta_{1}$, …, $\theta_{N}$, $\xi_{2}$, …, $\xi_{N-1}$.

The equations (103)-(108) are evolved in time using a built-in ODE solver in
Matlab.

Acknowledgments. L.O. acknowledges support from NSF Postdoctoral Fellowship
DMS-2001959 and thanks Becca Thomases, Yoichiro Mori, and Dallas Albritton for
helpful discussion.

## References

  * [1] V. Angeles, F. A. Godínez, J. A. Puente-Velazquez, R. Mendez-Rojano, E. Lauga, and R. Zenit.  Front-back asymmetry controls the impact of viscoelasticity on helical swimming.  Physical Review Fluids, 6(4):043102, 2021.
  * [2] S. Camalet and F. Jülicher.  Generic aspects of axonemal beating.  New Journal of Physics, 2(1):24, 2000.
  * [3] G. J. Elfring and G. Goyal.  The effect of gait on swimming in viscoelastic fluids.  Journal of Non-Newtonian Fluid Mechanics, 234:8–14, 2016.
  * [4] J. Espinosa-Garcia, E. Lauga, and R. Zenit.  Fluid elasticity increases the locomotion of flexible swimmers.  Physics of Fluids, 25(3):031701, 2013.
  * [5] H. C. Fu, T. R. Powers, and C. W. Wolgemuth.  Theory of swimming filaments in viscoelastic media.  Physical review letters, 99(25):258101, 2007.
  * [6] H. C. Fu, C. W. Wolgemuth, and T. R. Powers.  Beating patterns of filaments in viscoelastic fluids.  Physical Review E, 78(4):041913, 2008.
  * [7] H. C. Fu, C. W. Wolgemuth, and T. R. Powers.  Swimming speeds of filaments in nonlinearly viscoelastic fluids.  Physics of Fluids, 21(3):033102, 2009.
  * [8] G. R. Fulford, D. F. Katz, and R. L. Powell.  Swimming of spermatozoa in a linear viscoelastic fluid.  Biorheology, 35(4-5):295–309, 1998.
  * [9] F. A. Godínez, L. Koens, T. D. Montenegro-Johnson, R. Zenit, and E. Lauga.  Complex fluids affect low-reynolds number locomotion in a kinematic-dependent manner.  Experiments in Fluids, 56(5):1–10, 2015.
  * [10] R. E. Goldstein and S. A. Langer.  Nonlinear dynamics of stiff polymers.  Physical review letters, 75(6):1094, 1995.
  * [11] J. Gray and G. Hancock.  The propulsion of sea-urchin spermatozoa.  Journal of Experimental Biology, 32(4):802–814, 1955.
  * [12] M. Hines and J. Blum.  Bend propagation in flagella. i. derivation of equations of motion and their simulation.  Biophysical Journal, 23(1):41–57, 1978.
  * [13] N. C. Keim, M. Garcia, and P. E. Arratia.  Fluid elasticity can enable propulsion at low reynolds number.  Physics of Fluids, 24(8):081703, 2012.
  * [14] L. Landau and E. Lifschitz.  Theory of elasticity—course of theoretical physics, vol. 7 oxford.  UK: Butterwoth-Heinemann, 1986.
  * [15] E. Lauga.  Propulsion in a viscoelastic fluid.  Physics of Fluids, 19(8):083104, 2007.
  * [16] E. Lauga.  The fluid dynamics of cell motility, volume 62.  Cambridge University Press, 2020.
  * [17] C. Li, B. Qin, A. Gopinath, P. E. Arratia, B. Thomases, and R. D. Guy.  Flagellar swimming in viscoelastic fluids: role of fluid elastic stress revealed by simulations based on experimental data.  Journal of The Royal Society Interface, 14(135):20170289, 2017.
  * [18] C. Li, B. Thomases, and R. D. Guy.  Orientation dependent elastic stress concentration at tips of slender objects translating in viscoelastic fluids.  Physical Review Fluids, 4(3):031301, 2019.
  * [19] G. Li, E. Lauga, and A. M. Ardekani.  Microswimming in viscoelastic fluids.  Journal of Non-Newtonian Fluid Mechanics, 297:104655, 2021.
  * [20] B. Liu, T. R. Powers, and K. S. Breuer.  Force-free swimming of a model helical flagellum in viscoelastic fluids.  Proceedings of the National Academy of Sciences, 108(49):19516–19520, 2011.
  * [21] V. A. Martinez, J. Schwarz-Linek, M. Reufer, L. G. Wilson, A. N. Morozov, and W. C. Poon.  Flagellated bacterial motility in polymer solutions.  Proceedings of the National Academy of Sciences, 111(50):17771–17776, 2014.
  * [22] O. Maxian, A. Mogilner, and A. Donev.  Integral-based spectral method for inextensible slender fibers in stokes flow.  Physical Review Fluids, 6(1):014102, 2021.
  * [23] C. Moreau, L. Giraldi, and H. Gadêlha.  The asymptotic coarse-graining formulation of slender-rods, bio-filaments and flagella.  Journal of the Royal Society Interface, 15(144):20180235, 2018.
  * [24] Y. Mori and L. Ohm.  An error bound for the slender body approximation of a thin, rigid fiber sedimenting in stokes flow.  Research in the Mathematical Sciences, 7(2):1–27, 2020.
  * [25] Y. Mori and L. Ohm.  Accuracy of slender body theory in approximating force exerted by thin fiber on viscous fluid.  Studies in Applied Mathematics, 2021.
  * [26] Y. Mori and L. Ohm.  Well-posedness and applications of classical elastohydrodynamics for a swimming filament.  arXiv preprint arXiv:2208.03350, 2022.
  * [27] Y. Mori, L. Ohm, and D. Spirn.  Theoretical justification and error analysis for slender body theory.  Communications on Pure and Applied Mathematics, 73(6):1245–1314, 2020.
  * [28] Y. Mori, L. Ohm, and D. Spirn.  Theoretical justification and error analysis for slender body theory with free ends.  Archive for Rational Mechanics and Analysis, 235(3):1905–1978, 2020\.
  * [29] L. Ohm.  Remarks on regularized stokeslets in slender body theory.  Fluids, 6(8):283, 2021.
  * [30] O. Pironneau and D. Katz.  Optimal swimming of flagellated micro-organisms.  Journal of Fluid Mechanics, 66(2):391–415, 1974.
  * [31] E. E. Riley and E. Lauga.  Enhanced active swimming in viscoelastic fluids.  EPL (Europhysics Letters), 108(3):34003, 2014.
  * [32] E. E. Riley and E. Lauga.  Small-amplitude swimmers can self-propel faster in viscoelastic fluids.  Journal of theoretical biology, 382:345–355, 2015.
  * [33] D. Salazar, A. M. Roma, and H. D. Ceniceros.  Numerical study of an inextensible, finite swimmer in stokesian viscoelastic flow.  Physics of Fluids, 28(6):063101, 2016.
  * [34] X. Shen and P. E. Arratia.  Undulatory swimming in viscoelastic fluids.  Physical review letters, 106(20):208101, 2011.
  * [35] S. E. Spagnolie, B. Liu, and T. R. Powers.  Locomotion of helical bodies in viscoelastic fluids: enhanced swimming at large helical amplitudes.  Physical review letters, 111(6):068101, 2013.
  * [36] J. Teran, L. Fauci, and M. Shelley.  Viscoelastic fluid response can increase the speed and efficiency of a free swimmer.  Physical review letters, 104(3):038101, 2010.
  * [37] B. Thomases and R. D. Guy.  Mechanisms of elastic enhancement and hindrance for finite-length undulatory swimmers in viscoelastic fluids.  Physical review letters, 113(9):098102, 2014.
  * [38] B. Thomases and R. D. Guy.  The role of body flexibility in stroke enhancements for finite-length undulatory swimmers in viscoelastic fluids.  Journal of Fluid Mechanics, 825:109–132, 2017.
  * [39] A.-K. Tornberg and M. J. Shelley.  Simulating the dynamics and interactions of flexible fibers in stokes flows.  Journal of Computational Physics, 196(1):8–40, 2004.
  * [40] C. H. Wiggins and R. E. Goldstein.  Flexive and propulsive dynamics of elastica at low reynolds number.  Physical Review Letters, 80(17):3879, 1998.
  * [41] C. H. Wiggins, D. Riveline, A. Ott, and R. E. Goldstein.  Trapping and wiggling: elastohydrodynamics of driven microfilaments.  Biophysical Journal, 74(2):1043–1060, 1998.


# Entanglement entropy in holographic moving mirror and Page curve

Ibrahim Akala, Yuya Kusukia, Noburo Shibaa, Tadashi Takayanagia,b,c and Zixia
Weia aYukawa Institute for Theoretical Physics, Kyoto University,
Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan  bInamori Research
Institute for Science, 620 Suiginya-cho, Shimogyo-ku, Kyoto 600-8411 Japan
cKavli Institute for the Physics and Mathematics of the Universe,
University of Tokyo, Kashiwa, Chiba 277-8582, Japan

###### Abstract

We calculate the time evolution of entanglement entropy in two dimensional
conformal field theory with a moving mirror. For a setup modeling Hawking
radiation, we obtain a linear growth of entanglement entropy and show that
this can be interpreted as the production of entangled pairs. For a setup,
which mimics black hole formation and evaporation, we find that the evolution
follows the ideal Page curve. We perform these computations by constructing
the gravity dual of the moving mirror model via holography. We also argue that
our holographic setup provides a concrete model to derive the Page curve for
black hole radiation in the strong coupling regime of gravity.

††preprint: YITP-20-149; IPMU20-0121

##  I Introduction

Moving mirrors have been known for a while as a class of instructive models
that mimic Hawking radiation Hawking:1974sw  based on quantum field theory
Birrell:1982ix ; Davies:1976hi  where unitarity is manifest. On the other
hand, in the case of black hole evaporation, it has been a significant problem
to understand whether unitarity is maintained in the gravitational theory. One
manifestation of unitary black hole evaporation is the Page curve for the
entropy of Hawking radiation Page:1993df . Based on the fine grained entropy
formula RT ; HRT ; Faulkner:2013ana ; Engelhardt:2014gca , this has been
derived semiclassically for field theories coupled to gravity
Penington:2019npb ; Almheiri:2019psf ; Almheiri:2019hni  and confirmed by
direct gravity replica computations Penington:2019kki ; Almheiri:2019qdq . See
e.g. Akers:2019nfi ; Almheiri:2019yqk ; Rozali:2019day ; Chen:2019uhq ;
Bousso:2019ykv ; Kusuki:2019hcg ; Pollack:2020gfa ; Liu:2020gnp ;
Marolf:2020xie ; Piroli:2020dlx ; Balasubramanian:2020hfs ; Verlinde:2020upt ;
Chen:2020wiq ; Gautason:2020tmk ; Anegawa:2020ezn ; Giddings:2020yes ;
Hashimoto:2020cas ; Hartman:2020swn ; Agon:2020fqs ; Hollowood:2020cou ;
Krishnan:2020oun ; Alishahiha:2020qza ; Geng:2020qvw ; Chen:2020uac ;
Almheiri:2020cfm ; Li:2020ceg ; Chandrasekaran:2020qtn ; Bak:2020enw ;
Bousso:2020kmy ; Dong:2020uxp ; Akal:2020wfl ; Engelhardt:2020qpv ;
Karlsson:2020uga ; Chen:2020jvn ; Chen:2020tes ; Hartman:2020khs ;
Murdia:2020iac ; Altland:2020ccq ; Balasubramanian:2020xqf ; Sybesma:2020fxg ;
Stanford:2020wkf ; Ling:2020laa ; Chakravarty:2020wdm ; Bhattacharya:2020uun ;
Harlow:2020bee ; Chen:2020ojn ; Kirklin:2020zic ; Goto:2020wnk ; Hsin:2020mfa
for further progress along this direction. For recent related works refer to
Nomura:2019qps ; Laddha:2020kvp ; Akal:2020ujg .

In this article, we first present concrete calculations of entanglement
entropy in moving mirror setups and show that this leads to an ideal Page
curve. This itself provides a novel nonequilibrium setup, where quantum
entanglement evolves rapidly. Moreover, we present a close connection between
moving mirror models and black hole radiation via a particular version of the
anti de-Sitter (AdS)/conformal field theory (CFT) correspondence Ma , namely,
in the case when the CFT is defined on a manifold with a boundary Karch:2000gx
; Ta . For earlier studies of entanglement entropy in moving mirror models
refer to Bianchi:2014qua ; Hotta:2015huj ; Good:2016atu ; Chen:2017lum ;
Good:2019tnf .

Figure 1: The moving mirror setup (left) and its conformal transformation into
a static mirror (right). The Mirror trajectory is depicted by the thick curve.
The shaded region shown in the right panel corresponds to an inside horizon
region, which is missing in the left picture.

##  II Moving mirror from conformal maps

A moving mirror setup in two dimensions is specified by the trajectory of a
mirror profile $x=Z(t)$. We consider a CFT which lives on the right region,
i.e. $x\geq Z(t)$. A conformal transformation (here, we set $u=t-x$ and
$v=t+x$) Birrell:1982ix ; Davies:1976hi ,

$\displaystyle\tilde{u}=p(u),\ \ \ \tilde{v}=v,$ (1)

maps this into a simple setup with a static mirror $\tilde{u}-\tilde{v}=0$, as
depicted in Fig. 1. Here, we choose the function $p(u)$ such that the mirror
trajectory is given by $v=p(u)$, i.e.

$\displaystyle t+Z(t)=p\left(t-Z(t)\right).$ (2)

For example, we can calculate the energy stress tensor from the conformal
anomaly via the map (1), such that

$\displaystyle T_{uu}={c\over 24\pi}\left({3\over
2}\left({p^{\prime\prime}(u)\over
p^{\prime}(u)}\right)^{2}-{p^{\prime\prime\prime}(u)\over
p^{\prime}(u)}\right),$ (3)

where the components $T_{uv}$ and $T_{vv}$ are vanishing.

As an example of a CFT, consider a massless free scalar $\phi$. We impose the
Dirichlet boundary condition $\phi(t,Z(t))=0$ along the mirror trajectory. A
complete set of positive frequency solutions to the equations of motion
$\partial_{u}\partial_{v}\phi=0$, which satisfy the latter boundary condition,
reads

$\begin{split}\phi_{\omega}(t,x)=i(4\pi\omega)^{-1/2}\left(e^{-i\omega
v}-e^{-i\omega p(u)}\right).\end{split}$ (4)

Then, $\phi$ can be expanded in terms of these modes Birrell:1982ix  as

$\begin{split}\phi(t,x)=\int_{0}^{\infty}d\omega\left[a_{\omega}^{\text{in}}\phi_{\omega}+a_{\omega}^{\text{in}\dagger}\phi_{\omega}^{*}\right],\end{split}$
(5)

where $a_{\omega}^{\text{in}}$ and $a_{\omega}^{\text{in}\dagger}$ are the
annihilation and creation operators, respectively. The in-coming vacuum
$|0_{\text{in}}\rangle$ is defined by the state annihilated by
$a_{\omega}^{\text{in}}$ and the out-going vacuum $|0_{\text{out}}\rangle$ is
given by a Bogoliubov transformation of $|0_{\text{in}}\rangle$. The
expectation value of the energy stress tensor, i.e. $\langle
0_{\text{in}}|T_{uu}|0_{\text{in}}\rangle$, reproduces (3) for $c=1$.

Figure 2: The graphs in the left figure show the time evolution of
entanglement entropy $S_{A}$. We choose the end point of $A$ to be
$x_{0}=-t+\xi_{0}$, with $\xi_{0}=1$ (thick), $\xi_{0}=0.1$ (dashed) and
$\xi_{0}=0.01$ (dotted). We set $\beta=1$, $\epsilon=0.1$ and
$S_{\text{bdy}}=0$. The right figure shows a quasi-particle picture of
entanglement growth for the moving mirror. The black thick curve represents
the mirror trajectory $x=Z(t)$. The purple dotted curve describes a spacelike
curve defined by $v+p(u)=0$. The red line corresponds to the null line. The
entangled pair production occurs on the purple dashed curve.

##  III Moving mirror with horizon

For a typical example which models Hawking radiation from a black hole, we
choose

$\displaystyle p(u)=-\beta\log(1+e^{-{u\over\beta}}),$ (6)

where the parameter $\beta$ plays the role of an inverse temperature. Its
profile is depicted in the left picture of Fig. 1. In the early time limit
$t\to-\infty$, we have $Z(t)\simeq 0$, while in the late time limit, the
mirror trajectory gets almost lightlike, $Z(t)\simeq-t-\beta e^{-2t/\beta}$.
As depicted in Fig. 1, the region $\tilde{u}\geq 0$ in the extended
coordinates is missing for the original coordinates. This region is analogous
to the inside horizon region in the black hole formation process.

The energy stress tensor (3) reads

$\displaystyle T_{uu}={c\over
48\pi\beta^{2}}\left(1-{1\over(1+e^{u/\beta})^{2}}\right),$ (7)

which vanishes at early time $u\to-\infty$, and becomes a constant thermal
flux, $T_{uu}\simeq{c\over 48\pi\beta^{2}}$, at late time $u\to\infty$.

Let us calculate the entanglement entropy $S_{A}$ for a semi infinite
subsystem $A$ given by $[x_{0},\infty]$ at time $t$. We can calculate $S_{A}$
from the one point function of the twist operator CC04 ; HLW ; Casini:2009sr
on the upper half plane $\tilde{x}>0$ by using the conformal map (1) via the
replica method. We find

$\displaystyle S_{A}={c\over
6}\log{t+x_{0}-p(t-x_{0})\over\epsilon\sqrt{p^{\prime}(t-x_{0})}}+S_{\text{bdy}}.$
(8)

Here, $\epsilon$ is the UV cutoff (lattice spacing) of the CFT and
$S_{\text{bdy}}$ is the boundary entropy Affleck:1991tk ; CC04 ; Ta . Note
that this formula holds for any CFT.

If we fix the end point of the subsystem $A$, i.e. $x_{0}$, we can approximate
(8) at late time $t\to\infty$, and find

$\displaystyle S_{A}\simeq{c\over 12\beta}(t-x_{0})+{c\over
6}\log{t\over\epsilon}+S_{\text{bdy}}.$ (9)

The first term linear in $t$ arises from entangled pair production due to the
moving mirror, while the second $\log t$ term comes from the standard vacuum
entanglement as the length of the complement of $A$ grows linearly.

To study the first contribution in more detail, we allow changing the value of
$x_{0}$ time dependently as

$\displaystyle x_{0}(t)=-t+\xi_{0}.$ (10)

In the late time limit $u\to\infty$, we obtain

$\displaystyle S_{A}=\\!{c\over 6}\\!\log\\!{\left(\xi_{0}\\!+\\!\beta
e^{-(t\\!-\\!x_{0}(t))/\beta}\\!\right)\sqrt{1+e^{(t\\!-\\!x_{0}(t))/\beta}}\over\epsilon}+S_{\text{bdy}}.$

We choose $\xi_{0}$ to be positive, but sufficiently small. If the left end
point of $A$, given by $(u,v)=(2t-\xi_{0},\xi_{0})$, satisfies $v+p(u)>0$, we
get the linear growth (see the left panel in Fig. 2),

$\displaystyle S_{A}\simeq{c\over 6\beta}t+{c\over
6}\log{\xi_{0}\over\epsilon}+S_{\text{bdy}}.$ (12)

In this way, we may conclude that the entangled pair production occurs along
the spacelike curve $v+p(u)=0$, and the propagation of the entangled pairs
gives the linear growth of the entanglement entropy (12). This is sketched in
the right panel of Fig. 2. We can also confirm this from the free scalar
example (5), where the spacial distribution of the pair production looks like

$\displaystyle\langle 0_{\text{in}}|\phi(u_{1},v_{1})\phi(u_{2},v_{2})\int
d\omega
a_{\omega}^{\text{in}\dagger}a_{\omega}^{\text{in}\dagger}|0_{\text{in}}\rangle\\!\propto$
(13)
$\displaystyle\int\\!{d\omega\over\omega}\\!\left[\\!e^{\\!-i\omega(\\!v_{1}\\!+\\!p_{2}\\!)}\\!+\\!e^{\\!-i\omega(\\!v_{2}\\!+\\!p_{1}\\!)}\\!-\\!e^{\\!-i\omega(\\!v_{1}\\!+\\!v_{2}\\!)}\\!-\\!e^{\\!-i\omega(\\!p_{1}\\!+\\!p_{2}\\!)}\\!\right],$

where we have defined $p_{i}:=p(u_{i})$ for brevity. The first two terms are
divergent at $v_{1}+p_{2}=0$ and $v_{2}+p_{1}=0$. This shows that the
entangled pairs are produced along the curve $v+p(u)=0$, and they propagate in
opposite directions at the speed of light.

##  IV AdS/BCFT and entanglement entropy

To compute $S_{A}$ for generic subsystems, we need to specify the target CFT.
For holographic CFTs, we can calculate $S_{A}$ via the gravity dual of a CFT
defined on a manifold $M$ with a boundary $\partial M$, i.e. boundary CFT
(BCFT) Cardy2004 , known as AdS/BCFT Ta . In this description, the dual
geometry is given by extending the boundary $\partial M$ into the bulk AdS,
which leads to a codimension one surface $Q$, called the end of the world
brane. This surface $Q$ obeys the Neumann boundary condition

$\displaystyle K_{ab}-h_{ab}K+{\cal T}h_{ab}=0,$ (14)

where $h_{ab}$ is the induced metric and $K_{ab}$ is the extrinsic curvature.
The parameter ${\cal T}$ is the tension of the brane $Q$ and depends on the
boundary condition of the CFT at $\partial M$. The condition (14) implies the
presence of boundary conformal invariance. Refer to Takayanagi:2020njm  for an
equivalent formulation using Chern-Simons gravity, and to Sully:2020pza  for
comparisons with CFT calculations.

We can find a gravity dual by applying the following coordinate
transformation, which is a special case of Banados:1998gg ; Roberts:2012aq
and is a bulk extension of the map (1),

$\displaystyle U=p(u),\quad V=v+{p^{\prime\prime}(u)\over
2p^{\prime}(u)}z^{2},\quad\eta=z\sqrt{p^{\prime}(u)}$ (15)

on Poincare AdS3

$\displaystyle ds^{2}={d\eta^{2}-dUdV\over\eta^{2}}.$ (16)

Using (3), this leads to the metric

$\displaystyle ds^{2}={dz^{2}\over z^{2}}+{12\pi\over
c}T_{uu}(u)(du)^{2}-{1\over z^{2}}dudv.$ (17)

In Poincare AdS3 (16), by solving the boundary condition (14), the profile of
$Q$ is given by $X=-\lambda\eta$, where we have defined $\lambda={{\cal
T}\over\sqrt{1-{\cal{T}}^{2}}}$, and introduced new coordinates $U=T-X$ and
$V=T+X$. The metric on $Q$ is given by that of Poincare AdS2 (see the left
panel in Fig. 3),

$\displaystyle ds^{2}={(1+\lambda^{2})d\eta^{2}-dT^{2}\over\eta^{2}}.$ (18)

Thus, the gravity dual in terms of the $(U,V,\eta)$ coordinates is given by a
part of Poincare AdS3 defined by $X+\lambda\eta>0$. Note that the surface $Q$
at the boundary $z=0$ coincides with the mirror trajectory $v=p(u)$ via the
map (15).

The gravity dual in terms of the coordinates $(u,v,z)$, which is given by the
metric (17) and is sketched in the right panel of Fig. 3, only covers the
region $U<0$ as $U=p(u)$ is always negative for any $u$. This is the bulk
extension of the mentioned inner horizon region shown in the right panel in
Fig. 1. In the coordinates $(u,v,z)$, the metric of the brane $Q$ reads

$\displaystyle ds^{2}\\!\\!$ $\displaystyle=\\!\\!{dz^{2}\over
z^{2}}\\!\\!+\\!\\!\left(\\!{p^{\prime\prime}\over
zp^{\prime}}\\!\\!+\\!\\!{2\lambda\sqrt{p^{\prime}}\over
z^{2}}\\!\right)\\!dudz+\\!\\!\left(\\!{p^{\prime\prime 2}\over 4p^{\prime
2}}\\!\\!-\\!\\!{p^{\prime}\over z^{2}}\\!\\!+\\!\\!{\lambda
p^{\prime\prime}\over z\sqrt{p^{\prime}}}\\!\right)\\!du^{2}.$ (19)

This covers only the part $T<-\lambda\eta$ of (18). In this way, the gravity
dual of the moving mirror has a horizon, analogous to a single sided AdS black
hole.

Figure 3: Gravity dual of a moving mirror in the coordinates $(U,V,\eta)$
(left) and $(u,v,z)$ (right). We set ${\cal T}=0$. We also show the
computation of holographic entanglement entropy.

In general, the holographic entanglement entropy RT ; HRT  in AdS/BCFT can be
computed Ta  as

$\displaystyle S_{A}={1\over
4G_{N}}\mbox{Min}_{\Gamma_{A}}\left[A(\Gamma_{A})\right],$ (20)

where $A(\Gamma_{A})$ is the length of $\Gamma_{A}$ which satisfies
$\partial\Gamma_{A}=\partial A\cup\partial I_{s}$, where $I_{s}$ (i.e. island)
is a region on the surface $Q$. The three dimensional Newton constant is
denoted by $G_{N}$. The minimum in (20) is taken over all possible choices of
$I_{s}$ and $\Gamma_{A}$. When $A$ is an interval $[x_{0},x_{1}]$ at time $t$,
there appear to be two candidates for $\Gamma_{A}$. One is a connected
geodesic between $x=x_{0}$ and $x=x_{1}$. The other one is a union of two
disconnected ones, each of which departs from $x=x_{0}$ (or $x=x_{1}$) and
ends on $Q$, respectively.

Figure 4: The time evolution of entanglement entropy for moving mirror (6) in
free Dirac fermion CFT (left) and in holographic CFT (right). Here, we set the
subsystem to be $A=[Z(t)+0.1,Z(t)+10]$ with $\beta=0.1$, $\epsilon=0.1$ and
$S_{\text{bdy}}=0$. In the right, the thick line and the dashed line show the
disconnected and connected entanglement entropy, respectively.

In our setup, they are explicitly given by
$S_{A}=\mbox{Min}[S^{\text{con}}_{A},S^{\text{dis}}_{A}]$, where the
disconnected and connected geodesic contributions $S^{\text{dis}}_{A}$ and
$S^{\text{con}}_{A}$ read

$\displaystyle S^{\text{dis}}_{A}$ $\displaystyle={c\over
6}\log{t\\!+\\!x_{0}\\!-\\!p(t\\!-\\!x_{0})\over\epsilon\sqrt{p^{\prime}(t\\!-\\!x_{0})}}\\!+{c\over
6}\log{t\\!+\\!x_{1}\\!-\\!p(t\\!-\\!x_{1})\over\epsilon\sqrt{p^{\prime}(t\\!-\\!x_{1})}}$
$\displaystyle+2S_{\text{bdy}},$ $\displaystyle S^{\text{con}}_{A}$
$\displaystyle={c\over
6}\log{(x_{1}-x_{0})\left(p(t-x_{0})-p(t-x_{1})\right)\over\epsilon^{2}\sqrt{p^{\prime}(t-x_{0})p^{\prime}(t-x_{1})}}.$
(21)

The boundary entropy is a function of the tension and is given by
$S_{\text{bdy}}={c\over 6}\log\sqrt{(1+{\cal T})/(1-{\cal T})}$. When $A$ is
semi infinite, i.e. $x_{1}\to\infty$, we always have
$S_{A}=S^{\text{dis}}_{A}$, and this reproduces (8). When $A$ is a finite
interval, $S^{\text{dis}}_{A}$ is initially favored and this gives the linear
growth as in (12). At later time, $S^{\text{con}}_{A}$ is favored and this
leads to a saturation as depicted in the right panel of Fig. 4. We have also
plotted $S_{A}$ for the massless free Dirac fermion case, which is shown in
the left panel of Fig. 4. Refer to Appendix A for detailed computations of
$S_{A}$ in the free Dirac fermion and holographic CFT case.

##  V Page curve from moving mirror

A typical moving mirror model which mimics an evaporating black hole is found
by setting

$\displaystyle
p(u)=-\beta\log(1+e^{-{u\over\beta}})+\beta\log(1+e^{{u-u_{0}\over\beta}}),$
(22)

whose mirror trajectory $x=Z(t)$ and energy flux $T_{uu}$ are depicted in Fig.
5. When $\beta$ is small, we can approximate the trajectory as $Z(t)\simeq 0$
for $t<0$, $Z(t)\simeq-t$ for $0<t<u_{0}/2$, and $Z(t)=-u_{0}/2$ for
$t>u_{0}/2$. The energy flux is nonvanishing, $T_{uu}\simeq{c\over
48\pi\beta^{2}}$, namely, only for the period $0<u<u_{0}$.

Figure 5: The profile of moving mirror, with the creation of entangled pairs
and their reflection at the mirror, is depicted in the left picture. The right
graph is the energy density $T_{uu}$ plotted as a function of $u$. We set
$\beta=0.1$ and $u_{0}=5$.

We can again calculate the holographic entanglement entropy
$S_{A}=\mbox{Min}[S^{\text{con}}_{A},S^{\text{dis}}_{A}]$, using (21) as
plotted in Fig. 6. In particular, when $A$ is a semi infinite line, $S_{A}$
takes the form of the Page curve. For $0<t<u_{0}/4$, $S_{A}$ grows linearly
${dS_{A}\over dt}\simeq{c\over 6\beta}$, and for $u_{0}/4<t<u_{0}/2$, it
decreases linearly, i.e. ${dS_{A}\over dt}\simeq-{c\over 6\beta}$. Note that
as is clear from Fig. 6, the disconnected result $S^{\text{dis}}$ gives the
dominant contribution (refer to Hotta:2015huj  for an earlier calculation of
the connected result $S^{\text{con}}$).

The initial linear growth of $S_{A}$ can be understood as in the previous
example (12) by considering entangled pair production along the curve
$v+p(u)=0$. Moreover, the linear decay of $S_{A}$ is explained by reflections
of the left moving partner, as shown in the left picture of Fig. 5. When $A$
is a finite interval, we have two Page peaks. The first peak occurs when only
the originally right moving particles are crossing $A$. The second peak
appears when only the reflected particles are crossing $A$.

##  VI Brane world gravity and island

The gravity dual of our moving mirror setup can be interpreted in an
alternative way by regarding the surface $Q$ as an end of the world brane in
the brane world setup Randall:1999ee ; Randall:1999vf ; Karch:2000ct . This
situation is depicted in the left panel of Fig. 7. According to this
interpretation, the CFT defined in the region $x\geq Z(t)$ will be coupled to
a two dimensional gravity theory on $Q$. By estimating the effective Newton
constant on $Q$ via Kaluza-Klein reduction Randall:1999ee ; Randall:1999vf ;
Karch:2000ct ; Akal:2020wfl , which we denote by $G^{(Q)}_{N}$, we can show

$\displaystyle{1\over 4G^{(Q)}_{N}}=S_{\text{bdy}}.$ (23)

The gravitational entropy of AdS2, i.e. the brane $Q$, will thus be equal to
the boundary entropy $S_{\text{bdy}}$.

We can regard $S_{A}$ in (8) as the entanglement entropy of the subregion $A$
in a system consisting of a CFT on $x\geq Z(t)$ and a gravitational theory on
$Q$, glued along the moving mirror. Then, we can interpret the first and
second term in (8) as the bulk entropy contribution $S^{\text{bulk}}_{A\cup
I_{s}}$ and the area term ${\text{Area}(\partial I_{s})\over 4G_{N}}$,
respectively, in the island formula Penington:2019npb ; Almheiri:2019psf ;
Almheiri:2019hni ; Penington:2019kki ; Almheiri:2019qdq . Note that here we
use the standard formula for computing holographic entanglement entropy
without invoking the quantum extremal surface prescription. The density matrix
under consideration is pure and radiation is manifestly unitary.

In our moving mirror model, the entropy (8) shows linear growth in the region
defined by the equation $v+p(u)>0$. This part in the entanglement entropy
would arise from the island. This region is not covered by the coordinate
patch (19), as sketched in the left panel of Fig. 7.

Figure 6: The time evolution of holographic entanglement entropy follows the
Page curve. We choose $A$ to be a semi infinite line $A=[Z(t)+0.1,\infty]$ in
the left and a finite interval $A=[Z(t)+0.1,Z(t)+10]$ in the right. The thick
and dashed curves describe $S^{\text{dis}}_{A}$ and $S^{\text{con}}_{A}$,
respectively. We set $\beta=0.1,u_{0}=5$, $\epsilon=0.1$ and
$S_{\text{bdy}}=0$.

An interesting feature in our setup is the entanglement between the
gravitational theory on $Q$ and the CFT by the amount (23). Also it is
important to note that the presence of energy flux from the boundary is
different from standard BCFTs (Cardy states CardyS ) which have no energy flux
condition $T(w)-\bar{T}(\bar{w})=0$ at the boundary.

The radiation in the CFT looks similar to the setups Penington:2019npb ;
Almheiri:2019psf ; Almheiri:2019hni ; Penington:2019kki ; Almheiri:2019qdq ,
where the Page curve was derived. However, unlike these, in our BCFT model, we
find that there is no radiation present in the gravitational system on $Q$.
That is, the flux does not come from the gravitational system, but is created
on the boundary. Indeed, the holographic energy stress tensor
Balasubramanian:1999re  on $Q$ is proportional to $h_{ab}$ as it follows from
(14). Hence, it can just be regarded as a negative cosmological constant. The
absence of radiation from $Q$ is obviously consistent with the fact that the
mirror is completely reflective. We can regard our analysis as a derivation of
the Page curve in the strong coupling regime of gravity, while
Penington:2019npb ; Almheiri:2019psf ; Almheiri:2019hni ; Penington:2019kki ;
Almheiri:2019qdq  focus on the weak coupling regime.

Indeed, by changing the profile of the brane $Q$, we can deform our setup of
the moving mirror such that it incorporates radiation resulting from the
gravitational sector on $Q$, see right panel of Fig. 7. If we modify the
surface $Q$ to make it close to the standard AdS boundary located at
$z=\epsilon$, the matter energy stress tensor on $Q$ will be approximated by
(7). This provides a special and concrete example of the setup considered in
Almheiri:2019hni , see also Almheiri:2019yqk ; Rozali:2019day ; Chen:2019uhq ;
Chen:2020uac ; Bousso:2020kmy ; Akal:2020wfl ; Chen:2020tes  for related
works. The holographic dual of this modified setup is a two dimensional CFT
coupled to two dimensional gravity. One advantage of our procedure is that our
calculation based on the BCFT analysis is much easier than the one based on
the conformal welding problem Almheiri:2019psf ; Almheiri:2019qdq . An
interesting future direction will be relating the two realizations in an
explicit way Future .

Figure 7: The left picture shows the global spacetime of a CFT (left triangle)
and the brane world gravity on $Q$ (right triangle), which are attached along
the mirror trajectory. The island (doubled green line) is shown as well. The
right picture sketches the deformation of the gravity dual of the moving
mirror (top) into that with black hole radiation (bottom).

##  VII Conclusion

In this article, we have presented a gravity dual of two dimensional CFT with
a moving mirror, which mimics black hole formation and evaporation. We have
explicitly calculated the time evolution of entanglement entropy in the
presence of the mirror. We have found that it follows the ideal Page curve.
This can be explained by the creation of entangled particles, their
propagation, and reflection from the mirror. We have also discussed that
modifying the profile of the end of the world brane in the gravity dual
results in a model for two dimensional black hole radiation. In order to
understand unitary evolution for realistic black hole evaporation, we will
have to incorporate the singularity. We expect that the presence of spacelike
boundaries in the CFT and its gravity dual Akal:2020wfl ; Chen:2020tes  will
be relevant Future .

Acknowledgements We are grateful to Tatsuma Nishioka, Kotaro Tamaoka and
Tomonori Ugajin for useful comments on a draft of this paper. IA, YK, and ZW
are supported by the Japan Society for the Promotion of Science (JSPS). IA is
supportted by the Alexander von Humboldt (AvH) foundation. IA and TT are
supported by Grant-in-Aid for JSPS Fellows No. 19F19813. YK is supported by
Grant-in-Aid for JSPS Fellows No. 18J22495. NS is also supported by JSPS
KAKENHI Grant No. JP19K14721. TT is supported by the Simons Foundation through
the “It from Qubit” collaboration. TT is supported by Inamori Research
Institute for Science and World Premier International Research Center
Initiative (WPI Initiative) from the Japan Ministry of Education, Culture,
Sports, Science and Technology (MEXT). SN and TT are supported by JSPS Grant-
in-Aid for Scientific Research (A) No. 16H02182. TT is also supported by JSPS
Grant-in-Aid for Challenging Research (Exploratory) 18K18766. ZW is supported
by the ANRI Fellowship and Grant-in-Aid for JSPS Fellows No. 20J23116.

## References

  * (1) S. W. Hawking, “Particle Creation by Black Holes,” Commun. Math. Phys. 43 (1975), 199-220 [erratum: Commun. Math. Phys. 46 (1976), 206]
  * (2) N. Birrell and P. Davies, _Quantum Fields in Curved Space_ , Cambridge Monographs on Mathematical Physics. Cambridge Univ. Press, Cambridge, UK, 2, 1984\.
  * (3) P. Davies and S. Fulling, _Radiation from a moving mirror in two-dimensional space-time conformal anomaly_ , Proc. Roy. Soc. Lond. A 348 (1976) 393.
  * (4) D. N. Page, “Average entropy of a subsystem,” Phys. Rev. Lett. 71 (1993), 1291-1294 [arXiv:gr-qc/9305007 [gr-qc]]; “Information in black hole radiation,” Phys. Rev. Lett. 71 (1993), 3743-3746 [arXiv:hep-th/9306083 [hep-th]].
  * (5) S. Ryu and T. Takayanagi, “Holographic derivation of entanglement entropy from AdS/CFT,” Phys. Rev. Lett. 96 (2006) 181602 [hep-th/0603001].
  * (6) V. E. Hubeny, M. Rangamani and T. Takayanagi, “A Covariant holographic entanglement entropy proposal,” JHEP 0707 (2007) 062 [arXiv:0705.0016 [hep-th]].
  * (7) T. Faulkner, A. Lewkowycz and J. Maldacena, “Quantum corrections to holographic entanglement entropy,” JHEP 11, 074 (2013) [arXiv:1307.2892 [hep-th]].
  * (8) N. Engelhardt and A. C. Wall, “Quantum Extremal Surfaces: Holographic Entanglement Entropy beyond the Classical Regime,” JHEP 01, 073 (2015) [arXiv:1408.3203 [hep-th]].
  * (9) G. Penington, “Entanglement Wedge Reconstruction and the Information Paradox,” JHEP 09 (2020), 002 [arXiv:1905.08255 [hep-th]].
  * (10) A. Almheiri, N. Engelhardt, D. Marolf and H. Maxfield, “The entropy of bulk quantum fields and the entanglement wedge of an evaporating black hole,” JHEP 12 (2019), 063 [arXiv:1905.08762 [hep-th]].
  * (11) A. Almheiri, R. Mahajan, J. Maldacena and Y. Zhao, “The Page curve of Hawking radiation from semiclassical geometry,” JHEP 03 (2020), 149 [arXiv:1908.10996 [hep-th]].
  * (12) G. Penington, S. H. Shenker, D. Stanford and Z. Yang, “Replica wormholes and the black hole interior,” [arXiv:1911.11977 [hep-th]].
  * (13) A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian and A. Tajdini, “Replica Wormholes and the Entropy of Hawking Radiation,” JHEP 05 (2020), 013 [arXiv:1911.12333 [hep-th]].
  * (14) C. Akers, N. Engelhardt and D. Harlow, “Simple holographic models of black hole evaporation,” JHEP 08 (2020), 032 [arXiv:1910.00972 [hep-th]].
  * (15) A. Almheiri, R. Mahajan and J. Maldacena, “Islands outside the horizon,” [arXiv:1910.11077 [hep-th]]; A. Almheiri, R. Mahajan and J. E. Santos, “Entanglement islands in higher dimensions,” SciPost Phys. 9 (2020) no.1, 001 [arXiv:1911.09666 [hep-th]].
  * (16) M. Rozali, J. Sully, M. Van Raamsdonk, C. Waddell and D. Wakeham, “Information radiation in BCFT models of black holes,” JHEP 05 (2020), 004 [arXiv:1910.12836 [hep-th]].
  * (17) H. Z. Chen, Z. Fisher, J. Hernandez, R. C. Myers and S. M. Ruan, “Information Flow in Black Hole Evaporation,” JHEP 03 (2020), 152 [arXiv:1911.03402 [hep-th]].
  * (18) R. Bousso and M. Tomašević, “Unitarity From a Smooth Horizon?,” Phys. Rev. D 102, no.10, 106019 (2020) [arXiv:1911.06305 [hep-th]].
  * (19) Y. Kusuki, Y. Suzuki, T. Takayanagi and K. Umemoto, “Looking at Shadows of Entanglement Wedges,” [arXiv:1912.08423 [hep-th]].
  * (20) J. Pollack, M. Rozali, J. Sully and D. Wakeham, ‘Eigenstate Thermalization and Disorder Averaging in Gravity,” Phys. Rev. Lett. 125 (2020) no.2, 021601 [arXiv:2002.02971 [hep-th]].
  * (21) H. Liu and S. Vardhan, “A dynamical mechanism for the Page curve from quantum chaos,” [arXiv:2002.05734 [hep-th]].
  * (22) D. Marolf and H. Maxfield, “Transcending the ensemble: baby universes, spacetime wormholes, and the order and disorder of black hole information,” JHEP 08 (2020), 044 [arXiv:2002.08950 [hep-th]]; “Observations of Hawking radiation: the Page curve and baby universes,” [arXiv:2010.06602 [hep-th]].
  * (23) L. Piroli, C. Sünderhauf and X. L. Qi, “A Random Unitary Circuit Model for Black Hole Evaporation,” JHEP 04 (2020), 063 [arXiv:2002.09236 [hep-th]].
  * (24) V. Balasubramanian, A. Kar, O. Parrikar, G. Sárosi and T. Ugajin, “Geometric secret sharing in a model of Hawking radiation,” [arXiv:2003.05448 [hep-th]].
  * (25) H. Verlinde, “ER = EPR revisited: On the Entropy of an Einstein-Rosen Bridge,” [arXiv:2003.13117 [hep-th]].
  * (26) Y. Chen, X. L. Qi and P. Zhang, “Replica wormhole and information retrieval in the SYK model coupled to Majorana chains,” JHEP 06 (2020), 121 [arXiv:2003.13147 [hep-th]].
  * (27) F. F. Gautason, L. Schneiderbauer, W. Sybesma and L. Thorlacius, “Page Curve for an Evaporating Black Hole,” JHEP 05 (2020), 091 [arXiv:2004.00598 [hep-th]].
  * (28) T. Anegawa and N. Iizuka, “Notes on islands in asymptotically flat 2d dilaton black holes,” JHEP 07 (2020), 036 [arXiv:2004.01601 [hep-th]].
  * (29) S. B. Giddings and G. J. Turiaci, “Wormhole calculus, replicas, and entropies,” JHEP 09 (2020), 194 [arXiv:2004.02900 [hep-th]].
  * (30) K. Hashimoto, N. Iizuka and Y. Matsuo, “Islands in Schwarzschild black holes,” JHEP 06 (2020), 085 [arXiv:2004.05863 [hep-th]].
  * (31) T. Hartman, E. Shaghoulian and A. Strominger, “Islands in Asymptotically Flat 2D Gravity,” JHEP 07 (2020), 022 [arXiv:2004.13857 [hep-th]].
  * (32) C. A. Agón, S. F. Lokhande and J. F. Pedraza, “Local quenches, bulk entanglement entropy and a unitary Page curve,” JHEP 08 (2020), 152 [arXiv:2004.15010 [hep-th]].
  * (33) T. J. Hollowood and S. P. Kumar, “Islands and Page Curves for Evaporating Black Holes in JT Gravity,” JHEP 08 (2020), 094 [arXiv:2004.14944 [hep-th]].
  * (34) C. Krishnan, V. Patil and J. Pereira, “Page Curve and the Information Paradox in Flat Space,” [arXiv:2005.02993 [hep-th]].
  * (35) M. Alishahiha, A. Faraji Astaneh and A. Naseh, “Island in the Presence of Higher Derivative Terms,” [arXiv:2005.08715 [hep-th]].
  * (36) H. Geng and A. Karch, “Massive islands,” JHEP 09 (2020), 121 [arXiv:2006.02438 [hep-th]].
  * (37) H. Z. Chen, R. C. Myers, D. Neuenfeld, I. A. Reyes and J. Sandor, “Quantum Extremal Islands Made Easy, Part I: Entanglement on the Brane,” JHEP 10 (2020), 166 [arXiv:2006.04851 [hep-th]]; H. Z. Chen, R. C. Myers, D. Neuenfeld, I. A. Reyes and J. Sandor, “Quantum Extremal Islands Made Easy, Part II: Black Holes on the Brane,” [arXiv:2010.00018 [hep-th]]; J. Hernandez, R. C. Myers and S. M. Ruan, “Quantum Extremal Islands Made Easy, PartIII: Complexity on the Brane,” [arXiv:2010.16398 [hep-th].
  * (38) A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian and A. Tajdini, “The entropy of Hawking radiation,” [arXiv:2006.06872 [hep-th]].
  * (39) T. Li, J. Chu and Y. Zhou, “Reflected Entropy for an Evaporating Black Hole,” [arXiv:2006.10846 [hep-th]].
  * (40) V. Chandrasekaran, M. Miyaji and P. Rath, “Including contributions from entanglement islands to the reflected entropy,” Phys. Rev. D 102 (2020) no.8, 086009 [arXiv:2006.10754 [hep-th]].
  * (41) D. Bak, C. Kim, S. H. Yi and J. Yoon, “Unitarity of Entanglement and Islands in Two-Sided Janus Black Holes,” [arXiv:2006.11717 [hep-th]].
  * (42) R. Bousso and E. Wildenhain, “Gravity/ensemble duality,” Phys. Rev. D 102 (2020) no.6, 066005 [arXiv:2006.16289 [hep-th]].
  * (43) X. Dong, X. L. Qi, Z. Shangnan and Z. Yang, “Effective entropy of quantum fields coupled with gravity,” JHEP 10 (2020), 052 [arXiv:2007.02987 [hep-th]].
  * (44) I. Akal, Y. Kusuki, T. Takayanagi and Z. Wei, “Codimension two holography for wedges,” [arXiv:2007.06800 [hep-th]].
  * (45) N. Engelhardt, S. Fischetti and A. Maloney, “Free Energy from Replica Wormholes,” [arXiv:2007.07444 [hep-th]].
  * (46) A. Karlsson, “Replica wormhole and island incompatibility with monogamy of entanglement,” [arXiv:2007.10523 [hep-th]].
  * (47) H. Z. Chen, Z. Fisher, J. Hernandez, R. C. Myers and S. M. Ruan, “Evaporating Black Holes Coupled to a Thermal Bath,” [arXiv:2007.11658 [hep-th]].
  * (48) Y. Chen, V. Gorbenko and J. Maldacena, “Bra-ket wormholes in gravitationally prepared states,” [arXiv:2007.16091 [hep-th]].
  * (49) T. Hartman, Y. Jiang and E. Shaghoulian, “Islands in cosmology,” [arXiv:2008.01022 [hep-th]].
  * (50) C. Murdia, Y. Nomura and P. Rath, “Coarse-Graining Holographic States: A Semiclassical Flow in General Spacetimes,” Phys. Rev. D 102 (2020) no.8, 086001 [arXiv:2008.01755 [hep-th]].
  * (51) A. Altland and J. Sonner, “Late time physics of holographic quantum chaos,” [arXiv:2008.02271 [hep-th]].
  * (52) V. Balasubramanian, A. Kar and T. Ugajin, “Islands in de Sitter space,” [arXiv:2008.05275 [hep-th]]; “Entanglement between two disjoint universes,” [arXiv:2008.05274 [hep-th]].
  * (53) W. Sybesma, “Pure de Sitter space and the island moving back in time,” [arXiv:2008.07994 [hep-th]].
  * (54) D. Stanford, “More quantum noise from wormholes,” [arXiv:2008.08570 [hep-th]].
  * (55) Y. Ling, Y. Liu and Z. Y. Xian, “Island in Charged Black Holes,” [arXiv:2010.00037 [hep-th]].
  * (56) J. Chakravarty, “Overcounting of interior excitations: A resolution to the bags of gold paradox in AdS,” [arXiv:2010.03575 [hep-th]].
  * (57) A. Bhattacharya, A. Chanda, S. Maulik, C. Northe and S. Roy, “Topological shadows and complexity of islands in multiboundary wormholes,” [arXiv:2010.04134 [hep-th]].
  * (58) D. Harlow and E. Shaghoulian, “Global symmetry, Euclidean gravity, and the black hole information problem,” [arXiv:2010.10539 [hep-th]].
  * (59) Y. Chen and H. W. Lin, “Signatures of global symmetry violation in relative entropies and replica wormholes,” [arXiv:2011.06005 [hep-th]].
  * (60) J. Kirklin, “Islands and Uhlmann phase: Explicit recovery of classical information from evaporating black holes,” [arXiv:2011.07086 [hep-th]].
  * (61) K. Goto, T. Hartman and A. Tajdini, “Replica wormholes for an evaporating 2D black hole,” [arXiv:2011.09043 [hep-th]].
  * (62) P. S. Hsin, L. V. Iliesiu and Z. Yang, “A violation of global symmetries from replica wormholes and the fate of black hole remnants,” [arXiv:2011.09444 [hep-th]].
  * (63) Y. Nomura, “Spacetime and Universal Soft Modes — Black Holes and Beyond,” Phys. Rev. D 101 (2020) no.6, 066024 [arXiv:1908.05728 [hep-th]]; “Black Hole Interior in Unitary Gauge Construction,” [arXiv:2010.15827 [hep-th]].
  * (64) A. Laddha, S. G. Prabhu, S. Raju and P. Shrivastava, “The Holographic Nature of Null Infinity,” [arXiv:2002.02448 [hep-th]].
  * (65) I. Akal, “Universality, intertwiners and black hole information,” [arXiv:2010.12565 [hep-th]].
  * (66) J. M. Maldacena, Adv. Theor. Math. Phys. 2 (1998) 231 [Int. J. Theor. Phys. 38 (1999) 1113] [arXiv:hep-th/9711200];
  * (67) A. Karch and L. Randall, “Open and closed string interpretation of SUSY CFT’s on branes with boundaries,” JHEP 06 (2001), 063 [arXiv:hep-th/0105132 [hep-th]].
  * (68) T. Takayanagi, “Holographic Dual of BCFT,” Phys. Rev. Lett. 107 (2011) 101602 [arXiv:1105.5165 [hep-th]]; M. Fujita, T. Takayanagi and E. Tonni, “Aspects of AdS/BCFT,” JHEP 11 (2011), 043 [arXiv:1108.5152 [hep-th]].
  * (69) E. Bianchi and M. Smerlak, “Entanglement entropy and negative energy in two dimensions,” Phys. Rev. D 90, no.4, 041904 (2014) [arXiv:1404.0602 [gr-qc]].
  * (70) M. Hotta and A. Sugita, “The Fall of Black Hole Firewall: Natural Nonmaximal Entanglement for Page Curve,” PTEP 2015 (2015) no.12, 123B04 [arXiv:1505.05870 [gr-qc]].
  * (71) M. R. R. Good, K. Yelshibekov and Y. C. Ong, “On Horizonless Temperature with an Accelerating Mirror,” JHEP 03 (2017), 013 [arXiv:1611.00809 [gr-qc]].
  * (72) P. Chen and D. h. Yeom, “Entropy evolution of moving mirrors and the information loss problem,” Phys. Rev. D 96 (2017) no.2, 025016 [arXiv:1704.08613 [hep-th]].
  * (73) M. R. R. Good, E. V. Linder and F. Wilczek, “Moving mirror model for quasithermal radiation fields,” Phys. Rev. D 101, no.2, 025012 (2020) doi:10.1103/PhysRevD.101.025012 [arXiv:1909.01129 [gr-qc]].
  * (74) P. Calabrese and J. L. Cardy, “Entanglement entropy and quantum field theory,” J. Stat. Mech. 0406 (2004) P06002 [arXiv:hep-th/0405152 [hep-th]]; “Entanglement entropy and conformal field theory,” J. Phys. A 42 (2009), 504005 [arXiv:0905.4013 [cond-mat.stat-mech]].
  * (75) C. Holzhey, F. Larsen and F. Wilczek, “Geometric and renormalized entropy in conformal field theory,” Nucl. Phys. B 424 (1994) 443 [arXiv:hep-th/9403108 [hep-th]].
  * (76) H. Casini and M. Huerta, “Entanglement entropy in free quantum field theory,” J. Phys. A 42 (2009), 504007 [arXiv:0905.2562 [hep-th]].
  * (77) I. Affleck and A. W. W. Ludwig, “Universal noninteger ’ground state degeneracy’ in critical quantum systems,” Phys. Rev. Lett. 67 (1991), 161-164
  * (78) J. L. Cardy, _Boundary conformal field theory_ , arXiv: hep-th/0411189.
  * (79) T. Takayanagi and T. Uetoko, “Chern-Simons Gravity Dual of BCFT,” [arXiv:2011.02513 [hep-th]].
  * (80) J. Sully, M. Van Raamsdonk and D. Wakeham, “BCFT entanglement entropy at large central charge and the black hole interior,” [arXiv:2004.13088 [hep-th]].
  * (81) M. Banados, “Three-dimensional quantum geometry and black holes,” AIP Conf. Proc. 484 (1999) no.1, 147-169 [arXiv:hep-th/9901148 [hep-th]].
  * (82) M. M. Roberts, _Time evolution of entanglement entropy from a pulse_ , JHEP 12 (2012) 027 [arXiv:1204.1982].
  * (83) L. Randall and R. Sundrum, “A Large mass hierarchy from a small extra dimension,” Phys. Rev. Lett. 83 (1999), 3370-3373 [arXiv:hep-ph/9905221 [hep-ph]].
  * (84) L. Randall and R. Sundrum, “An Alternative to compactification,” Phys. Rev. Lett. 83 (1999), 4690-4693 [arXiv:hep-th/9906064 [hep-th]].
  * (85) A. Karch and L. Randall, “Locally localized gravity,” JHEP 05 (2001), 008 [arXiv:hep-th/0011156 [hep-th]].
  * (86) J. L. Cardy, “Boundary Conditions, Fusion Rules and the Verlinde Formula,” Nucl. Phys. B 324 (1989), 581-596.
  * (87) V. Balasubramanian and P. Kraus, “A Stress tensor for Anti-de Sitter gravity,” Commun. Math. Phys. 208 (1999), 413-428 [arXiv:hep-th/9902121 [hep-th]].
  * (88) I. Akal, Y. Kusuki, N. Shiba, T. Takayanagi and Z. Wei, work in progress.

##  Appendix A Appendix A: Calculating entanglement entropy in two dimensional
CFTs with a moving mirror

In the following, we present details of our computations of entanglement
entropy in the presence of a moving mirror for CFTs in two dimensions. The
replica trick is often used to compute the entanglement entropy in quantum
field theories CC04 . For two dimensional CFTs, in particular, the $n$-th
Rényi entropy of a single interval, $A=[x_{0},x_{1}]$, can be evaluated as

$\displaystyle S^{(n)}_{A}={1\over
1-n}\log\langle\sigma_{n}(t,x_{0})\bar{\sigma}_{n}(t,x_{1})\rangle.$ (24)

Here, $\sigma_{n}(t,x_{0})$ and $\bar{\sigma}_{n}(t,x_{1})$ are twist
operators which act as primaries with conformal weights

$\displaystyle h_{n}=\bar{h}_{n}={c\over 24}\left(n-{1\over n}\right).$ (25)

The expectation value $\langle\cdots\rangle$ is evaluated by the path integral
on the manifold with a moving mirror. Once an analytical form of the Rényi
entropy is found, we can perform the analytic continuation to obtain the
entanglement entropy. To evaluate the 2-point function of twist operators, we
use the following conformal map

$\displaystyle\tilde{u}=p(u),\quad\tilde{v}=v$ (26)

to map the moving mirror setup into the right half plane (RHP) of
$\mathbb{R}^{1,1}$, which is parameterized by

$\displaystyle(\tilde{t},\tilde{x})=\left({\tilde{u}+\tilde{v}\over
2},{-\tilde{u}+\tilde{v}\over 2}\right).$ (27)

Under this conformal transformation, we have

$\displaystyle\langle\sigma_{n}(t,x_{0})\bar{\sigma}_{n}(t,x_{1})\rangle=$
(28)
$\displaystyle\left(p^{\prime}(u_{0})p^{\prime}(u_{1})\right)^{h_{n}}\cdot\langle\tilde{\sigma}_{n}\left(\tilde{t}_{0},\tilde{x}_{0}\right)\tilde{\bar{\sigma}}_{n}\left(\tilde{t}_{1},\tilde{x}_{1}\right)\rangle_{\rm
RHP}.$

Thus, our task is reduced to evaluating a 2-point function of twist operators
on the RHP. As concrete examples, in what follows, we will do so for free
massless Dirac fermion CFT and holographic CFT.

###  A.1 Dirac fermion CFT

For the massless Dirac fermion CFT, after performing doubling trick and
bosonization (see e.g. Casini:2009sr ), we have

$\displaystyle\langle\tilde{\sigma}_{n}\left(\tilde{t}_{0},\tilde{x}_{0}\right)\tilde{\bar{\sigma}}_{n}\left(\tilde{t}_{1},\tilde{x}_{1}\right)\rangle_{\rm
RHP}$ $\displaystyle=$
$\displaystyle\sqrt{\langle\tilde{{\sigma}}_{n}\left(\tilde{t}_{1},-\tilde{x}_{1}\right)\tilde{\bar{\sigma}}_{n}\left(\tilde{t}_{0},-\tilde{x}_{0}\right)\tilde{\sigma}_{n}\left(\tilde{t}_{0},\tilde{x}_{0}\right)\tilde{\bar{\sigma}}_{n}\left(\tilde{t}_{1},\tilde{x}_{1}\right)\rangle_{\mathbb{R}^{1,1}}}$
$\displaystyle\propto$
$\displaystyle\left({(\tilde{u}_{0}-\tilde{v}_{1})(\tilde{v}_{0}-\tilde{u}_{1})\over(\tilde{u}_{0}-\tilde{u}_{1})(\tilde{v}_{0}-\tilde{v}_{1})(\tilde{u}_{0}-\tilde{v}_{0})(\tilde{u}_{1}-\tilde{v}_{1})}\right)^{2h_{n}}.$
(29)

By plugging this into (A) and (24) with proper $p(u)$, and taking the limit
$n\rightarrow 1$, we will get the entanglement entropy for the moving mirror
setup in Dirac fermion CFT.

###  A.2 Holographic CFT

In holographic CFT, the 2-point function can be evaluated as (see e.g.
Sully:2020pza )

$\displaystyle\langle\tilde{\sigma}_{n}\left(\tilde{t}_{0},\tilde{x}_{0}\right)\tilde{\bar{\sigma}}_{n}\left(\tilde{t}_{1},\tilde{x}_{1}\right)\rangle_{\rm
RHP}$ (30) $\displaystyle=$
$\displaystyle\max\begin{cases}\langle\tilde{\sigma}_{n}\left(\tilde{t}_{0},\tilde{x}_{0}\right)\tilde{\bar{\sigma}}_{n}\left(\tilde{t}_{1},\tilde{x}_{1}\right)\rangle_{\mathbb{R}^{1,1}}\\\
g^{2(1-n)}\prod_{i\in\\{0,1\\}}\langle\tilde{\sigma}_{n}\left(\tilde{t}_{i},\tilde{x}_{i}\right)\tilde{\bar{\sigma}}_{n}\left(\tilde{t}_{i},-\tilde{x}_{i}\right)\rangle_{\mathbb{R}^{1,1}}^{1/2}\end{cases}.$

Here, $g$ is a constant which depends on the details of the boundary condition
on the moving mirror and results in the boundary entropy $S_{\text{bdy}}=\log
g$. These two cases correspond to the connected channel and the disconnected
channel, respectively. We can obtain (A.2) by assuming the simple Wick
contraction rule (so called generalized free field prescription) by factoring
the correlation functions into two point functions of the twist operators.
This is justified due to the large $c$ factorization property of holographic
CFTs. The disconnected channel arises by considering the mirror operator
across the boundary. Therefore, we call the von Neumann entropy obtained from
both channels the connected entanglement entropy, $S_{A}^{\text{con}}$, and
the disconnected entanglement entropy, $S_{A}^{\text{dis}}$, respectively. The
physical entanglement entropy is

$\displaystyle
S_{A}=\mathrm{Min}\left[S_{A}^{\text{con}},S_{A}^{\text{dis}}\right].$ (31)

Remembering that the 2-point function on $\mathbb{R}^{1,1}$ is given by

$\displaystyle\langle\tilde{\sigma}_{n}\left(\tilde{t},\tilde{x}\right)\tilde{\bar{\sigma}}_{n}\left(\tilde{t}^{\prime},\tilde{x}^{\prime}\right)\rangle_{\mathbb{R}^{1,1}}={1\over|(\tilde{t}-\tilde{t}^{\prime})^{2}-(\tilde{x}-\tilde{x}^{\prime})^{2}|^{2h_{n}}},$
(32)

plugging these into (A) and (24) with proper $p(u)$ and taking the limit
$n\rightarrow 1$, we get the connected and the disconnected entanglement
entropy for the moving mirror setup in holographic CFT. This leads to the
results in (21). As discussed in the main text, for deriving (21) we have
performed a holographic calculation by evaluating the length of geodesics in
the bulk. The computation here is performed on the CFT side. In fact, we can
easily find out that evaluating the 2-point function in the (dis)connected
channel on the CFT side is exactly equivalent to evaluating the length of the
(dis)connected geodesic on the AdS side.


11institutetext: BTU Cottbus-Senftenberg, Lehrstuhl Theoretische Informatik,
Platz der Deutschen Einheit 1, 03046 Cottbus, Germany, https://www.b-tu.de/
11email<EMAIL_ADDRESS>22institutetext: Nantes Université, École Centrale
Nantes, CNRS, LS2N, UMR 6004, F-44000 Nantes, France

# On Hybrid Gene Regulatory Networks

Adrian Wurm ID 1 Honglu Sun2

###### Abstract

In this work, we study a class of hybrid dynamical systems called hybrid gene
regulatory networks (HGRNs) which was proposed to model gene regulatory
networks. In HGRNs, there exist well-behaved trajectories that reach a fixed
point or converge to a limit cycle, as well as chaotic trajectories that
behave non-periodic or indeterministic. In our work, we investigate these
irregular behaviors of HGRNs and present theoretical results about the
decidability of the reachability problem, the probability of indeterministic
behavior of HGRNs, and chaos especially in 2-dimensional HGRNs.

###### Keywords:

Hybrid system Decidability Chaos Bisimulation.

##  1 Introduction

Hybrid systems are dynamical systems that admit both continuous and discrete
behavior. They were initially proposed to model the interactions between the
discrete digital world and the continuous physical world [2]. In recent years,
hybrid systems were also used in system biology to study biological
networks[10, 8, 11].

In this work, we study a class of hybrid systems called hybrid gene regulatory
networks (HGRN) [9, 7], which is an extension of Thomas’ discrete modeling
framework [19, 20]. These hybrid systems were proposed to model gene
regulatory networks, a framework to describe the regulation influences between
genes.

Two widely used formalisms to study gene regulatory networks are discrete
models (like Boolean networks [14]) and continuous models (differential
equations [6, 1] and stochastic models [13]). The major advantage of HGRNs
compared to discrete models is that they can provide and process temporal
information of the systems. Compared to purely continuous models, the
dynamical properties of HGRNs are easier to analyze, which makes them a more
reasonable choice to apply on large and complex instances.

HGRNs are similar to piecewise-constant derivative systems (PCD systems) [4]
which are a special case of hybrid automata [3]. The major difference between
HGRNs and PCD systems, as for example in [4, 5, 17], is the existence of
sliding modes, which means that when a trajectory reaches a black wall (a
boundary of the discrete state which can be reached but cannot be crossed by
trajectories), the trajectory does not end there but is instead forced to move
along this black wall. There exist other methods to define the behavior of
trajectories on a black wall [12, 16] which are different from the sliding
modes in HGRNs.

One important research topic of dynamical models of gene regulatory networks
is how to analyze the reachability of states, for instance, whether a
trajectory starting in a certain state can reach another given state. These
questions arise in cell reprogramming [15]. In [18], the reachability problem
of HGRNs is investigated and methods are proposed to analyze the reachability
of regular trajectories, meaning those that reach a fixed point or converge to
a limit cycle.

In HGRNs however, there also exists so-called chaotic trajectories that cross
an infinite sequence of discrete regions so that this sequence does not end up
in a cycle. For now, our understanding of this irregular behaviour is quite
limited. Further understanding of these trajectories can potentially lead to
the development of new methods to analyze the reachability problem for them.
In this work, we want to fill the gaps in this research field by investigating
irregular behavior occuring in HGRNs. Our major contributions are:

  * •

We prove that the reachability problem for HGRNs is undecidable and that this
enforces the existence of chaotic trajectories.

  * •

We prove that under reasonable conditions, the probability for a random
trajectory to admit indeterminism tends to 1 as the dimension of the system
tends to infinity.

  * •

We argue that in two-dimensional HGRNs, chaos is not possible without
indeterminism, and conclude that the probability for a random two-dimensional
HGRN to admit chaos with indeterminism is 0.

The paper is organized as follows: In Section 2 we collect basic notions and
recall the definition of HGRNs as well as their properties. In Section 3 we
prove that HGRNs are Turing-powerful and conclude that chaotic trajectories
must exist. In Section 4 we give some structural results on indeterminism and
provide conditions as well as probability estimations for the bifurcation of
trajectories.

The paper ends with some open questions.

##  2 Preliminaries and Network Reachability Problems

We start by defining the problems we are interested in; here, we follow the
definitions and notions of [18] for everything related to hybrid gene
regulatory networks.

Intuitively, a HGRN is a an $N$-dimensional rectangular cuboid/hyperrectangle
with integer edge lengths decomposed into hypercubes with unit edge length.
Each such hypercube $d$ is called a discrete state and is associated with a
celerity $c_{d}$ enforcing the direction in which a trajectory is allowed to
move inside $d$, so a trajectory may only change its direction in the moment
of transition between two discrete states.

######  Definition 1 (Hybrid gene regulatory network (HGRN))

Let $N\in\mathbb{N}$ and let $n_{i}\in\mathbb{N}$ for each
$i\in\\{1,...,N\\}$. A _Hybrid gene regulatory network_ (HGRN) on $N$ genes
with $n_{i}$ levels is a pair $\mathcal{H}=(E,c)$, where

$E:=\left\\{d\in\mathbb{N}^{N}\mid\forall i\in\left\\{1,2,...,N\right\\}:0\leq
d_{i}\leq n_{i}\right\\}$

is called the set of discrete states, and $c:E\rightarrow\mathbb{R}^{N}$ is
called the celerity function.

######  Definition 2 (Hybrid state of HGRN)

A _hybrid state_ of a HGRN is a pair $h=(\pi,d)$ containing a _fractional
part_ $\pi\in[0,1]^{N}$, and a discrete state $d\in E$. $E_{h}$ is the set of
all hybrid states and $E_{sh}$ is the set of all finite or infinite sequences
of hybrid states:

$E_{sh}=\left\\{(h_{0},h_{1},...,h_{m})\in(E_{h})^{m+1}\mid
m\in\mathbb{N}\cup\\{\infty\\}\right\\}.$

Note that $E_{h}=:E_{h}^{1}\subseteq E_{sh}$.

######  Definition 3 (Boundary)

A boundary of a discrete state $d$ is a set of hybrid states defined by
$e(i,\pi_{0},d):=\left\\{(\pi,d)\in E_{h}\mid\pi^{i}=\pi_{0},\right\\}$, where
$i\in\left\\{1,2,...,N\right\\}$ and $\pi_{0}\in\left\\{0,1\right\\}$.

######  Definition 4 (Hybrid trajectory)

Let $t_{0}\in\mathbb{R}^{+}\cup\\{\infty\\}$ and $[0,t_{0}]$ a time interval.
For a function $\tau:[0,t_{0}]\rightarrow E_{sh}$, we denote by abuse of
notation by $\pi:[0,t_{0}]\rightarrow E_{sh}$ the projection of $\tau$ on the
second component.

For $\tau(t)=h$, we denote by
$\frac{d\tau(t^{+})}{dt}:=\frac{d\pi(t^{+})}{dt}$ the temporal derivative of
the fractional part of $\tau$, $\frac{d\tau(t^{-})}{dt}$ and
$\frac{d\tau(t)}{dt}$ are defined analogous.

We call $\tau$ a hybrid trajectory, if for every time $t\in[0,t_{0}]$ with
$\tau(t)=h$ and $h=(h_{0},...,h_{m})\in E_{h}^{m+1}$ a sequence of $m+1$
hybrid states $h_{j}=(\pi_{j},d_{j})$, the following conditions are satisfied:

  * •

For $t=0$ we have $m=0$, and $\tau(0)\in E_{h}$ is not on any boundary.

  * •

If $h_{0}$ does not belong to any boundary, then $m=0$ and
$\frac{d\tau(t)}{dt}=c(d_{0})$ for $t>0$ and
$\frac{d\tau(t^{+})}{dt}=c(d_{0})$ for $t=0$

  * •

If $h_{0}$ only belongs to one boundary $e$, let $e$ be the upper boundary in
the $i^{th}$ dimension (the following is analogous when $e$ is the lower
boundary). In case $d_{0}^{i}$ is not the maximal discrete level of gene $i$,
let the discrete state on the other side of $e$ be noted as $d_{r}$, where
$d_{0}^{k}=d_{r}^{k}$ for all $k\neq i$, and $d_{0}^{i}+1=d_{r}^{i}$. There
are four possible cases:

    * –

If $c(d_{0})^{i}=0$, then $m=0$ and the trajectory from the current state will
slide along the boundary $e$ with speed
$\frac{d\tau(t^{+})}{dt}=\frac{d\tau(t^{-})}{dt}=c_{0}$, which is then called
a _neutral boundary_ of $d_{0}$.

    * –

If $c(d_{0})^{i}>0$, and either $d_{0}^{i}$ is the maximal discrete level of
the $i^{th}$ gene, or $d_{0}^{i}$ is not the maximal discrete level of the
$i^{th}$ gene but the $i^{th}$ component of $c(d_{r})$ is negative, then $m=0$
and the trajectory from the current state will slide along the boundary $e$,
which is called an _attractive boundary_ of $d_{0}$. We then have that
$\frac{d\tau(t)}{dt}^{k}=c(d_{0})^{k}$ for all $k\neq i$,
$\frac{d\tau(t^{+})}{dt}^{i}=0$ and
$\frac{d\tau(t^{-})}{dt}^{i}\in\\{0,c(d_{0})^{i}\\}$.

    * –

If $c(d_{0})^{i}>0$, $d_{0}^{i}$ is not the maximal discrete level of the
$i^{th}$ gene, and the $i^{th}$ component of $c_{r}$ is positive, then the
trajectory from the current state instantly crosses the boundary $e$ and enter
the discrete state $d_{r}$. We then have $m=1$ and
$\tau(t)=((\pi_{0},d_{0}),(\pi_{1},d_{1}))$, where $d_{1}=d_{r}$,
$\pi_{0}^{j}=\pi_{1}^{j}$ for all if $j\neq i$, and $\pi_{1}^{i}=0$. Such a
boundary is called an _output boundary_ and the counterpart of $d_{0}$ is
called an _input boundary_ of $d_{r}$. We then demand that
$\frac{d\tau(t^{+})}{dt}=c(d_{r})$ and $\frac{d\tau(t^{-})}{dt}=c(d_{0})$.

  * •

If $h$ belongs to several boundaries, then the previous cases can be mixed:

    * –

If in these boundaries there is no output boundary, then the trajectory from
the current state will exit all input boundaries and slide along all
attractive or neutral boundaries.

    * –

If in these boundaries there is only one output boundary, then the trajectory
from the current state will cross this output boundary.

    * –

If in these boundaries there are several output boundaries, then the
trajectory from the current state can cross any of them, but can only cross
one boundary at one time, which causes indeterministic behavior.

A hybrid trajectory $\tau$ is called

  * •

_finite_ , if $t_{0}\neq\infty$ and the sequence of states that are reached is
finite,

  * •

_periodic_ , if it is defined on $t_{0}=\infty$ and $\exists T,t_{1}>0,\forall
t\in[t_{1},\infty[,\tau(t)=\tau(t+T)$,

  * •

_chaotic_ , if the sequence of discrete states that it passes is infinite and
non-periodic, meaning it does not end up in a cycle and

  * •

_regular_ , if it is not chaotic.

######  Definition 5 (Reachability Problem)

The reachability problem for HGRNs is the following: Given a HGRN and two sets
$S,T$ of hybrid states in form of linear inequalities, the question is whether
there exists a trajectory $\tau$ starting in a hybrid state $s\in S$ and
reaching at some time $t$ a state $f\in F$, meaning

$\exists\tau\exists s\in S,f\in F,t>0:\tau(0)=s\land\tau(t)=f.$

If such a trajectory exists, we say that $F$ is reachable from $T$.

Note the special cases that $S,T$ might contain only one hybrid state or that
they may be a whole discrete state. It was shown in [18] that the reachability
problem is decidable for regular trajectories.

######  Definition 6 (Area of attraction)

Consider a subset $A$ of hybrid states of a HGRN. The area of attraction
$\varphi(A)$ is the set of all hybrid states with the property that every
trajectory starting in them will eventually intersect $A$. The extended area
of attraction $\Phi(A)$ is the set of all hybrid states with the property that
there exists a trajectory starting in them that will at some point intersect
$A$.

###### Example 1

Fig. 1 is a sketch of a HGRN on three genes with levels
$n_{1}=2,n_{2}=4,n_{3}=2$, the arrows indicate the celerities where an
encircled circle indicates an upward movement and an encircled X a downward
movement, missing arrows indicate that the celerity is the zero-vector. The
left side of the figure illustrates level 1 of gene 3, the right side level 2.

Figure 1: Binary expansion

The network computes the binary expansion of a number in $[0,1]$ in the
following way: If a trajectory passes $A=e(2,1,(1,2,1))$ (blue area),
interpret its distance to the rightmost part $A\cap B,B=e(1,1,(1,2,1))$ as its
current value $x$. If $x<\frac{1}{2}$, the trajectory will enter the small
circle of discrete states $(1,2,1),(2,2,1),(2,1,1)$ and $(1,1,1)$ and finally
end up in $A$ again with value $2x$. If $x>\frac{1}{2}$, it will pass the big
circle of discrete states $(1,2,1)$ (where $\frac{1}{2}$ is subtracted from
the value)$,(1,3,1),(1,3,2),(1,4,2)$ (where its value is multiplied with two)
$,(2,4,2),(2,3,2),(2,2,2),(2,1,1)$ and $(1,1,1)$ and again reach $A$, this
time with value $2(x-\frac{1}{2})$. If one follows the trajectory and
interprets the first case as 0 and the second case as 1, one obtains the
binary expansion of the initial $x$.

If for example we start in $a=((\frac{3}{16},\frac{3}{16}),(1,1,1))$, we
obtain the trajectory drawn in red. It passes $A$ at value $\frac{3}{8}$, so
it runs through the small circle and enters $A$ again at value $\frac{3}{4}$,
therefore entering the big circle ending up in $A$ again at value
$\frac{1}{2}$. So far, we obtained the binary expansion 0.01, but now the
trajectory becomes indeterministic because the celerities allow a bifurcation
at $((1,1),(1,2,1))$. This corresponds to the fact that $\frac{3}{8}$ has the
two binary expansions $0.011$ and $0.010\bar{1}$.

Observe that we obtain a periodic trajectory for every rational value $x$ and
a chaotic trajectory for $x$ irrational.

##  3 Turing completeness and chaos

Turing-completeness was already shown for PCD systems in [4], a model related
to HGRNs. We would now like to prove it for HGRNs, for they are more
restrictive in the sense that they only allow a change of derivatives at
integer coordinate changes and admit indeterminism with positive probability.

###### Theorem 3.1

HGRNs are Turing powerful.

###### Proof

Let $M=(Q,\Sigma,\Gamma,\delta,q_{1},\square,F)$ be a deterministic Turing-
Machine starting on an empty tape, $Q=\\{q_{1},...,q_{n}\\},F=\\{q_{n}\\}$.
Assume that all symbols are encoded in binary, so $\Sigma=\Gamma=\\{0,1\\}$
and that the encoding of $\square$ consists only of zeroes. Also assume that
for each two states $q_{i}$ and $q_{j}$ there are transitions in only one
direction, either $q_{i}\rightarrow q_{j}$ or $q_{j}\rightarrow q_{i}$ but not
both, a property that can be guaranteed by introducing auxiliary states. We
construct an equivalent HGRN $H$ with the following genes:

$\bullet$

One with two discrete levels storing the current value called $G_{c}$

$\bullet$

$G_{\ell,1}$ with three discrete levels and $G_{\ell,2}$ with two discrete
levels for the tape values left of the head, where the current value is stored
as the binary expansion with the most significant bits corresponding to the
values closest to the head in the fractional part of the upmost discrete state
of $G_{\ell,1}$.

$\bullet$

$G_{r,1}$ and $G_{r,2}$ analogous for the tape values right of the head.

$\bullet$

One called $G_{i}$ for every $q_{i}\in Q$ with two discrete levels each.

$\bullet$

Two called $G_{\alpha}$ and $G_{\beta}$ with two discrete levels each, keeping
track of the current calculation step as in Fig. 2.

The whole network is therefore of shape
$[0,2]\times[[0,3]\times[0,2]]^{2}\times[0,2]^{n}\times[0,2]^{2}$ and our
trajectory starts in
$((0,0),[(2,0),(0,0)]^{2},(1,1),[(0,0)]^{n-1},(0,1),(0,0))$.

Let $(q_{i},a)\rightarrow(q_{j},b,d)$ be a transition in $\delta$, without
loss of generality $d=r$ stating that the head moves to the right ($d=\ell$ is
analogous and $d=n$ is easily derived). We choose the celerities as follows:

Figure 2: Calculation step accounting

In the first calculation step $\alpha=\beta=0$ we put $q_{j}$ into the upmost
position. Note that we do not loose the information on which state we started
in because we assumed that there is no transition from $q_{j}$ to $q_{i}$.

The celerities describing how to update the left side tape start with making
the bits less significant by one position because the head recedes by one
position, which corresponds to dividing the numerical value of $G_{\ell,1}$ by
two, see Fig 3, left. This is done in the first calculation step
$\alpha=\beta=0$ along with adding half of $b$ to it, see Fig. 3, middle,
because $b$ is on the position left of the head after that transition.

In the second calculation step $\alpha=0,\beta=1$, this value is then copied
back from $G_{\ell,2}$ to the original position, see Fig. 3, right,
$G_{\ell,1}$, both $G_{\ell,1}$ and $G_{\ell,2}$ remain inactive for the last
two steps.

Figure 3: Left side tape update

The celerities describing how to update the right side tape work similar. We
start by multiplying the value by two, see Fig. 4 left, because each bit now
becomes more significant by one position. This is done in the first
calculation step $\alpha=\beta=0$, where we also catch the trajectory on the
border between $G_{r,1}=1$ and $G_{r,1}=2$ if it is in $G_{r,2}=1$ and let it
pass until the lower boundary of $G_{r,1}=1$ if $G_{r,1}=2$, in order to keep
track of what our new current value $G_{c}$ will be in the next Turing step.

Figure 4: Right side tape update

Note that if the initial value was $\frac{1}{2}$, we have an indeterminism
when hitting the lower boundary of $G_{r,1}=3$ and the upper boundary of
$G_{r,2}=1$. In the second calculation step $\alpha=0,\beta=1$, we therefore
adjust the celerities to make sure that we end in $G_{r,2}=2$ whenever the
initial value was greater than or equal $\frac{1}{2}$, see Fig. 4 middle (the
case that all positions on the right side tape are 1’s is not possible).

In the third and fourth calculation step $\alpha=1,\beta=1$ and
$\alpha=1,\beta=0$, we pop the most significant bit of the right side tape to
be the new $G_{c}$, which means that if the value of $G_{r,1}$ is greater than
1 we subtract 1 from $G_{r,1}$ and add move $G_{c}$ to the upmost position,
otherwise to the lowest, see Fig. 4 right.

Note that we loose the information on whether we moved left or right, so we
actually have the celerities of the fourth step on the right side of the tape
for the left side of the tape as well. However, the left side is not effected
by that because it stays in $G_{\ell,2}=3$ from the end of step two on and
therefore never enters a hypercube with non-zero celerities induced be the
update of $G_{c}$.

In the fourth step we also bring $G_{i}$ to the lowest possible point. By
assumption of anti-symmetry, $i$ can still be distinguished from $j$ at that
point without the knowledge of $d$.$\hfill\blacksquare$

###### Corollary 1

The reachability problem for HGRNs is undecidable.

###### Proof

The halting problem can be reformulated as the reachability problem of state
$q_{n}$.$\hfill\blacksquare$

Note that the trajectories that are obtained by the Turing machines that
behave chaotic automatically lead to chaotic trajectories, for the sequence of
reached states cannot have a repeating pattern.

##  4 Indeterminism

Most attempts to solve verification tasks such as the reachability problem
become a lot more difficult as soon as trajectories are not regular any more.
It is therefore desirable to estimate how likely a trajectory is to be chaotic
or indeterministic.

###### Example 2

In opposition to what is claimed in [18], indeterminism will happen with non-
zero probability in networks with at least 3 genes having at least two
discrete levels each when all celerities and the starting point are chosen
uniformly random in $(0,1)$. Fig. 5 is a sketch of a HGRN with 3 Genes
$G_{1},G_{2}$ and $G_{3}$ with two levels 0 and 1 each that allows the
indeterministic red trajectory starting in point $a$ in $(0,0,0)$.

Figure 5: Indeterminism with positive probability

The trajectory reaches the attractive boundary to discrete state 1 of $G_{1}$,
slides along until it reaches the attractive boundary to discrete state 1 of
$G_{2}$ and slides "upwards" along this intersection of discrete domains until
it reaches the hybrid state $((1,1,1),(0,0,0))$. The border to the discrete
state 1 of $G_{3}$ is an output boundary, so it is crossed and the hybrid
state $((1,1,0),(0,0,1))$ is immediately reached. Now the boundaries to the
discrete levels 1 of $G_{1}$ and $G_{2}$ both are output boundaries, so the
trajectory may pass any of them, leading to indeterminism. Note that we did
not rely on probability-zero-events; The probability of the celerity in the
first hypercube to point into the middle in all three dimensions is
$(\frac{1}{2})^{3}$, the first two reached boundaries are then attractive with
probability $\frac{1}{2}$ each as well as the third one being and output
boundary and the two new boundaries after the crossing are output boundaries
with probability $\frac{1}{4}$ each. The accumulated probability is then
$\frac{1}{8}\cdot(\frac{1}{2})^{2}\cdot\frac{1}{2}\cdot\frac{1}{4}^{2}=(\frac{1}{2})^{10}$.
This probability significantly rises once we increase the size of the network,
meaning the number of genes as well as their dimensions.

We can be more general:

###### Theorem 4.1

Let all celerities of a HGRN with $n$ genes and at least two discrete levels
per gene as well as the starting point of a trajectory be randomly chosen by
the uniform distribution in $[0,1]^{n}$. Then the probability for this
trajectory to lead to indeterminism tends to 1 as $n$ tends to $\infty$.

###### Proof

First, note that for each gene we have that an inner boundary is reached at
some point with probability at least $\frac{1}{2}$. These events are
independent, so it follows by the central limit theorem that the probability
of at least one fourth of all genes reaching an inner boundary tends to 1.
Each of these $\frac{n}{4}$ boundaries is attractive with probability
$\frac{1}{2}$, if we group them into $\frac{n}{12}$ triples of boundaries,
each triple starts with two attractive boundaries followed by an output
boundary with probability $\frac{1}{8}$ as in Example 2. The probability that
after the crossing of such a third boundary, the two other boundaries become
output boundaries is $(\frac{1}{2})^{2}$ each, so the accumulated probability
for each triple to lead to indeterminism is at least
$\frac{1}{8}\cdot(\frac{1}{2})^{2}\cdot(\frac{1}{2})^{2}=(\frac{1}{2})^{7}$.
The probability that none of the $\frac{n}{12}$ triples leads to indeterminism
is therefore at most
$(\frac{2^{7}-1}{2^{7}})^{\frac{n}{12}}\xrightarrow{n\rightarrow\infty}0$.$\hfill\blacksquare$

Note that Example 2 relied on the third dimension, for hitting two edges at
the same time without sliding before is highly unlikely:

###### Theorem 4.2

Random HGRNs in two dimensions admit chaos only with probability 0.

###### Proof

Since a trajectory cannot cross itself, it is obvious that chaos can only be
obtained by bifurcation.

A trajectory starts in the inner of a hypercube with probability 1. In order
to bifurcate, it needs to reach two output boundaries at the same time,
meaning it has to hit one of finitely many points. Reaching such an
intersection point $q$ by sliding towards it along an attractive boundary will
not lead to indeterminism at $q$ because the boundary which one is sliding
along is not an output boundary. A point $p$ in a two-dimensional HGRN that
has an area of attraction with non-zero measure must either have integer
coordinates connected to an attractive boundary or be reachable from such a
point, which is a zero-set.

The area of attraction of a point that admits two output boundaries is
therefore one-dimensional with probability 1, meaning it is a zero-set.
$\hfill\blacksquare$

Note that an indeterministic network such as in Example 2 inherently enables
chaotic trajectories if both choices lead back into the area of attraction of
the point $a$ at which the path was bifurcating, therefore forming two cycles
$T_{1},T_{2}$ with the same starting point $a$.

Figure 6: Guaranteed indeterminism with positive probability

The idea is sketched in Fig. 6 for two dimensions.

Take an irrational number, for example $\pi$, and when reaching $a$ for the
$n$-th time, follow the cycle $T_{i}$ iff the $n$-th position in the binary
expansion of $\pi$ is $i$. The picture shows indeterministic behavior in two
dimensions, which only happens with probability zero, it is however obvious
that this construction can be lifted into three dimensions using Example 2 as
a gadget to obtain a similar chaotic path with positive probability.

Fig. 6 is also an example for a two-dimensional HGRN that leads to
indeterminism independent by the starting point of the trajectory, and
therefore to chaos, because every possible trajectory in $\mathcal{N}$ will at
some point pass $a$, meaning $\varphi(a)=\Phi(a)=\mathcal{N}$.

This kind of chaos differs from the chaos that makes HGRNs Turing-complete. It
depends on the choices made on the way rather than the restrictions of the
network, and another trajectory starting in the same point might be periodic.
To examine more precisely how the reachability problem and the identification
of chaos relate to those, from now on we call a chaotic path that reaches a
certain point an infinite number of times recurrent, and otherwise say that
its chaos is proper.

The proof of Theorem 3.1 relies on proper chaos, so reachability could still
be decidable in recurrent chaos.

The following is a direct consequence of [4]:

###### Lemma 1

A two-dimensional HGRN never admits proper chaos.

We now want to define a property that most HGRNs of natural origin share:

###### Definition 7

For a HGRN $\mathcal{N}$ with two discrete levels for every gene and
$H=(V,E,\rho)$ a complete weighted directed graph with $V=\\{1,...,N\\}$ the
genes of $\mathcal{N}$ and $\rho:E\rightarrow\mathbb{Q}$, we say that
$\mathcal{N}$ is induced by $G$, iff for every $\alpha\in\\{0,1\\}^{N}$ we
have that celerities are the weighted sums over the incoming edges

$c(\alpha)_{j}=\sum\limits_{\begin{subarray}{c}i=1\\\
\alpha_{i}=1\end{subarray}}^{N}\rho(i,j)-\sum\limits_{\begin{subarray}{c}i=1\\\
\alpha_{i}=0\end{subarray}}^{N}\rho(i,j)$

We call $H$ the influence graph of $\mathcal{N}$.

We want to pose the question whether our problems become easier when
inducedness is given.

Note that induced HGRNs form a zero-set among all possible HGRNs.

###### Theorem 4.3

It can be decided whether a HGRN is induced or not in P. Moreover, the
influence graph is unique and can be computed in polynomial time.

###### Proof

When given all celerity vectors, the decisional as well as the computational
problem are linear equation systems on rational domains and therefore solved
efficiently. It remains to show that there can only be one solution.

Suppose there were two graphs $H_{1}$ and $H_{2}$ inducing the same HGRN, then
for all dimensions $i,j\in\\{1,...,n\\}$ we have that

$\rho_{1}(i,j)=\frac{c([0]^{i-1},1,[0]^{n-i})_{j}-c(0,...,0)_{j}}{2}=\rho_{2}(i,j)$

which gives us uniqueness.$\hfill\blacksquare$

The constellation in Example 2 can be obtained by an induced HGRN with
positive probability, for example by

##  5 Conclusion and Further Questions

We examined the reachability problem for HGRNs and showed in what way the size
and dimension of the net influence chaos and indeterminism.

Further open questions are:

1.)

Are induced HGRNs still Turing-powerful?

2.)

Does Theorem 4.1 still hold for induced HGRNs?

3.)

For $H_{\mathcal{N}}$ the set of all possible trajectories in $\mathcal{N}$,
how do the sets $R_{t}:=\\{h(t)\mid h\in H_{\mathcal{N}}\\}$ behave for large
$t$? With what probability and how fast do they approach a zero-set?

4.)

Can these sets $R_{t}$ be computed efficiently?

5.)

Can HGRNs be expressed by induced HGRNs with the help of auxiliary dimensions?

Acknowledgment: We want to thank Klaus Meer for helpful discussion.

## References

  * [1] Sofia Almeida, Madalena Chaves, and Franck Delaunay.  Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids.  Royal Society Open Science, 7(2):192054, 2020.
  * [2] Rajeev Alur, Costas Courcoubetis, Thomas A Henzinger, and Pei-Hsin Ho.  Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems.  In International Hybrid Systems Workshop, pages 209–229. Springer, 1991.
  * [3] Rajeev Alur, Costas Courcoubetis, Thomas A Henzinger, and Pei-Hsin Ho.  Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems.  Technical report, Cornell University, 1993.
  * [4] Eugene Asarin, Oded Maler, and Amir Pnueli.  Reachability analysis of dynamical systems having piecewise-constant derivatives.  Theoretical Computer Science, 138(1):35–65, 1995.
  * [5] Eugene Asarin, Venkatesh P Mysore, Amir Pnueli, and Gerardo Schneider.  Low dimensional hybrid systems–decidable, undecidable, don’t know.  Information and Computation, 211:138–159, 2012.
  * [6] Debashis Barik, William T Baumann, Mark R Paul, Bela Novak, and John J Tyson.  A model of yeast cell-cycle regulation based on multisite phosphorylation.  Molecular systems biology, 6(1):405, 2010.
  * [7] Jonathan Behaegel, Jean-Paul Comet, Gilles Bernot, Emilien Cornillon, and Franck Delaunay.  A hybrid model of cell cycle in mammals.  Journal of bioinformatics and computational biology, 14(01):1640001, 2016.
  * [8] Madalena Chaves and Miguel Preto.  Hierarchy of models: From qualitative to quantitative analysis of circadian rhythms in cyanobacteria.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 23(2):025113, 2013.
  * [9] Emilien Cornillon, Jean-Paul Comet, Gilles Bernot, and Gilles Enée.  Hybrid gene networks: a new framework and a software environment.  advances in Systems and Synthetic Biology, 2016.
  * [10] Etienne Farcot and Jean-Luc Gouzé.  Periodic solutions of piecewise affine gene network models with non uniform decay rates: the case of a negative feedback loop.  Acta Biotheoretica, 57(4):429–455, 2009.
  * [11] Elena Firippi and Madalena Chaves.  Topology-induced dynamics in a network of synthetic oscillators with piecewise affine approximation.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(11):113128, 2020.
  * [12] Jean-Luc Gouzé and Tewfik Sari.  A class of piecewise linear differential equations arising in biological models.  Dynamical systems, 17(4):299–316, 2002.
  * [13] Guy Karlebach and Ron Shamir.  Modelling and analysis of gene regulatory networks.  Nature reviews Molecular cell biology, 9(10):770–780, 2008.
  * [14] Stuart A Kauffman.  Metabolic stability and epigenesis in randomly constructed genetic nets.  Journal of theoretical biology, 22(3):437–467, 1969.
  * [15] Hugues Mandon, Cui Su, Stefan Haar, Jun Pang, and Loïc Paulevé.  Sequential reprogramming of boolean networks made practical.  In Computational Methods in Systems Biology: 17th International Conference, CMSB 2019, Trieste, Italy, September 18–20, 2019, Proceedings 17, pages 3–19. Springer, 2019.
  * [16] Erik Plahte and Sissel Kjøglum.  Analysis and generic properties of gene regulatory networks with graded response functions.  Physica D: Nonlinear Phenomena, 201(1-2):150–176, 2005.
  * [17] Andrei Sandler and Olga Tveretina.  Deciding reachability for piecewise constant derivative systems on orientable manifolds.  In Reachability Problems: 13th International Conference, RP 2019, Brussels, Belgium, September 11–13, 2019, Proceedings 13, pages 178–192. Springer, 2019.
  * [18] Honglu Sun, Maxime Folschette, and Morgan Magnin.  Reachability analysis of a class of hybrid gene regulatory networks.  In International Conference on Reachability Problems, pages 56–69. Springer, 2023.
  * [19] René Thomas.  Boolean formalization of genetic control circuits.  Journal of theoretical biology, 42(3):563–585, 1973.
  * [20] René Thomas.  Regulatory networks seen as asynchronous automata: a logical description.  Journal of theoretical biology, 153(1):1–23, 1991.