# A Quantum Optimal Control Problem with State Constrained Preserving
Coherence

Nahid Binandeh Dehaghani and Fernando Lobo Pereira SYSTEC - Research Center
for Systems and Technologies, FEUP - Faculty of Engineering, Porto University,
Rua Dr. Roberto Frias sn, i219, 4200-465 Porto, Portugal
nahid<EMAIL_ADDRESS>
###### Abstract

In this work, we address the problem of maximizing fidelity in a quantum state
transformation process controlled in such a way as to keep decoherence within
given bounds. We consider a three-level $\Lambda$-type atom subjected to
Markovian decoherence characterized by non-unital decoherence channels. We
introduce fidelity as the performance index for the quantum state
transformation process since the goal is to maximize the similarity of the
final state density operator with the one of the desired target state. We
formulate the quantum optimal control problem with state constraints where the
later reflect the fact that the decoherence level remains within a pre-defined
bound. Optimality conditions in the form of a Maximum Principle of Pontryagin
in Gamkrelidze’s form are given. These provide a complete set of relations
enabling the computation of the optimal control strategy. This is a novel
approach in the context of quantum systems.

##  I INTRODUCTION

The article is focused on the application of Maximum Principle of Pontryagin
in the Gamkrelidze’s form, [1] for state constrained optimal control problems
arising in the control of quantum systems. To the best of our knowledge, this
is the first instance of a formulation of a optimal quantum control problem
with state constraints, and, also, for which, this type of necessary
conditions of optimality have been applied. In this work, a three-level
$\Lambda$-type atom subjected to Markovian decoherence characterized by non-
unital decoherence channels is considered for which we preserve the quantum
coherence within some bound specified by the state constraints. Then, by
considering in the quantum state transformation process, as performance index,
the quantum fidelity of density operators, we aim at maximizing the similarity
of the desired target density operator with the final one.

With the increasing progress of quantum technology, the control of quantum
phenomena has become a significant issue in a variety of areas such as quantum
computing, quantum information processing, quantum optics, quantum chemistry,
atomic and molecular physics, and nuclear magnetic resonance (NMR). Quantum
control theory investigates how to drives quantum system dynamics by means of
external controls so that tasks are performed, [2, 3, 4, 5],e.g., the
application of quantum optimal control for quantum sensing, [6], and quantum
simulation, [7].

Due to the extra noise and disturbances of closed-loop control methods in
fragile quantum systems, the configuration of control law is considered
without experimental feedback or measurement, as open-loop, in most quantum
control protocols. Hence, optimal control theory (OCT) can be used as a strong
tool to control quantum processes while minimizing a cost function according
to the properties of quantum system, e.g., time optimal, [8], energy optimal,
[9], etc. Moreover, in comparison with other control methods, OCT [2] has the
key superiority of the very high versatility in the control problem
formulation, which may include a wide variety of constraints, such as state,
control, and mixed state-and-control and constraints. Necessary conditions of
optimality for optimal control problem in the form of the Pontryagin’s Maximum
Principle (PMP) has been already considered for a long time for quantum
systems, [10, 2]. However, to the best of our knowledge, state constraints
have not been considered. The reason for this relies on the fact that, in the
presence of state constraints, the multipliers, in general, include Borel
measures which make virtually impossible to use them in computationally
efficient indirect methods based these optimality conditions. However, in [1,
11], a Maximum Principle of Pontryagin in the Gamkrelidze form was derived
which under some regularity conditions, involve nontrivial multipliers which
exhibit a degree of regularity that make indirect methods computationally
effective.

The control of quantum decoherence, [12], has attracted much attention as an
important issue which has been remained unsolved yet. Loss of coherence means
non-unitary evolutions, [13, 14], leading to the loss of information and
probability leakage towards the environment, which causes the lack of
preservation of quantum superposition and entanglement, [15]. Under perfect
isolation, a quantum system would indefinitely keep its coherence, however,
the manipulation and investigation of an isolated quantum system is not
possible. For instance, during the process of testing or manipulating a
physical system in order to obtain numerical results, coherence is shared with
the environment, and will be lost over periods of time. Until now, quantum
presservation has been studied in several schemes such as quantum error-
correction codes, [16], error-avoiding codes, [17], decoherence-free subspace,
[18], dynamical decoupling, [19], tracking control, [20]. In addition, a
feedback control strategy based on quantum weak measurements, [21], has been
investigated, however, the strategy was restricted to single-qubit systems and
unital decoherence channels, while high-dimensional quantum systems have
attracted much more attention in the applications of quantum control.
Moreover, spontaneous emission, which can be modelled through non-unital
decoherence channels, demonstrate properties different from unital decoherence
channels and is typically used for decoherence modelling. Here, we take three-
level quantum systems and spontaneous emission into account.

Amongst the main three quantum control problems, including quantum state
preparation, transformation, and gate synthesis, the evolution of a quantum
state from some given initial state into a desired target state is the basis
of quantum computing, [22], coherent control, [23], etc. The performance of
such applications usually require a criterion, named fidelity, in order to be
quantified. By means of comparison, fidelity provides a primary measure of
similarity of two quantum states, for instance, in order to understand how
accurately some given approximate calculation replicates some desired density
operator, [24]. However, fidelity can also be used for analysing quantum
channels instead of pure states or density operators, [25]. The concept of
fidelity has been extensively used in quantum memories, [26], quantum phase-
space simulations, [27], quantum teleportation, [28], etc. In quantum
teleportation, fidelity between the initial state and the final state is the
standard measure used to quantify the effectiveness of the process, [28].

This work has been organized as follows: first we express the master equation
indicating the mathematical model of a Markovian open quantum system and the
system Hamiltonian. Then, we formulate the quantum coherence for non-unital
decoherence channel and the ways to quantifying fidelity. We then present a
novel formulation of the optimal control problem in the context of quantum
physics with state constraints. In the next section, we present the solution
to the proposed quantum control problem by means of Maximum Principle of
Pontryagin in Gamkrelidze form for the first time in quantum context.

##  II System Dynamics

According to Liouville equation we have

$\dot{\rho}(t)=-\frac{i}{\hbar}\left[H\left(t\right),\rho(t)\right]$ (1)

where the quantum state is represented by the density operator $\rho$, and
$\hbar$ is the reduced Planck’s constant. The system Hamiltonian
$\displaystyle H(t)=H_{0}+H_{C}(u(t))=H_{0}+\sum_{k}u_{k}(t)H_{k}$ determines
the controlled evolution, where $H_{0}$ and $H_{k}$ indicate the free and
interaction Hamiltonians, respectively. The system control can be realized by
a set of control functions $u_{k}(t)$ in $\mathbb{R}$, which is coupled to the
quantum system via time independent interaction Hamiltonians $H_{k}$,
$k=0,1,2,\dots$. The evolution of the Hamiltonian system is unitary for closed
quantum systems, so the spectrum of the quantum state is preserved. However,
in many practical applications, the quantum system is open, and its evolution
cannot be described by unitary transformations. The mathematical model of a
Markovian open quantum system can be expressed by the master equation, [29,
30, 31, 32],

$\dot{\rho}(t)=-\frac{i}{\hbar}\left[H(t),\rho(t)\right]+{\mathcal{D}}(\rho(t))$
(2)

where the dissipator $\mathcal{D}{{\rho}\left(t\right)}$ reflects the
interactions between the system and the environment through different
decoherence channels and is represented by

$\\!\\!\\!\\!\\!{\mathcal{D}}(\rho(t))\\!=\\!\sum_{k}\\!\gamma_{k}\\!\left(\\!L_{k}\rho(t)L_{k}^{\dagger}\\!-\\!\frac{1}{2}L_{k}^{\dagger}L_{k}\rho(t)\\!-\\!\frac{1}{2}\rho(t)L_{k}^{\dagger}L_{k}\\!\right)$
(3)

in which $L_{k}$ indicates the set of Lindblad operators, [29], and
$\gamma_{k}$ is the damping rate. By using the notation
$vec\left(\left|\psi\right\rangle\left\langle\xi\right|\right)=\left|\xi\right\rangle\otimes\left|\psi\right\rangle$,
we can write (2) in vector space as

$\frac{d}{dt}vec\left(\rho\right)=\hat{\mathcal{L}}vec\left(\rho\right)$ (4)

where

$\displaystyle\\!\\!\\!\hat{\mathcal{L}}=-i\left(I\otimes H-{{H}^{T}}\otimes
I\right)$ (5)
$\displaystyle\\!\\!\\!\\!+\sum\limits_{k}{{{\gamma}_{k}}\left(L_{k}^{*}\otimes{{L}_{k}}-\frac{1}{2}I\otimes
L_{k}^{\dagger}{{L}_{k}}-\frac{1}{2}{{\left(L_{k}^{\dagger}{{L}_{k}}\right)}^{T}}\otimes
I\right)}$

in which $I$ indicates the identity matrix.

##  III Specification of the quantum control problem

Consider a $\Lambda$-type atom subjected to Markovian decoherence. In
${\displaystyle\Lambda}$-type atoms, there are 3 states:
${\displaystyle|a\rangle}$, ${\displaystyle|b\rangle}$, and
${\displaystyle|e\rangle}$. In this type, ${\displaystyle|e\rangle}$ is at the
highest energy level, while ${\displaystyle|a\rangle}$ and
${\displaystyle|b\rangle}$ are at lower levels. The ground states
${\displaystyle|a\rangle}$ and ${\displaystyle|b\rangle}$ can couple to the
the excited state ${\displaystyle|e\rangle}$ with resonance frequencies
(coupling constants)
${{\omega}_{1}}={\left({{E}_{e}}-{{E}_{a}}\right)}/{\hbar}\;$ and
${{\omega}_{2}}={\left({{E}_{e}}-{{E}_{b}}\right)}/{\hbar}$, respectively, in
which ${{E}_{a}}$, ${{E}_{b}}$, and ${{E}_{e}}$ indicate the corresponding
eigenvalues. ${{\omega}_{3}}$ can also be defined in a similar way to
demonstrate the resonant frequency between ${\displaystyle|a\rangle}$ and
${\displaystyle|b\rangle}$ for the atom and laser field. The free Hamiltonian
of such a system can be expressed as

${{H}_{0}}=\frac{{{\omega}_{1}}}{3}\left({{\rho}_{e}}-{{\rho}_{a}}\right)+\frac{{{\omega}_{2}}}{3}\left({{\rho}_{e}}-{{\rho}_{b}}\right)+\frac{{{\omega}_{3}}}{3}\left({{\rho}_{a}}-{{\rho}_{b}}\right)$
(6)

where ${{\rho}_{j}}=\left|j\right\rangle\left\langle j\right|$, and the
constant energy term ${\left({{E}_{a}}+{{E}_{b}}+{{E}_{e}}\right)}/{3}$ is
ignored. Since The excited state is not stable, two electrons in state
${\displaystyle|e\rangle}$ decay to the two ground states
${\displaystyle|a\rangle}$ and ${\displaystyle|b\rangle}$ at rates
$\gamma_{a}$ and $\gamma_{b}$, and two kinds of emission are radiated, [33].
By assuming that the decay process is Markovian, the decoherence channel can
be characterized by the Lindblad operator
${{L}_{k}}=\left|k\right\rangle\left\langle e\right|,$ $k=a,b$. in the
dissipator equation (3). In order to inhibit the decay process, we can apply a
classical field such that drive the states between $|a\rangle$ and
$|e\rangle$. By assuming that the transition dipole moments for the linearly
polarized field are real, one can express the control field in the dipole
approximation as
$E\left(t\right)=u\left(t\right)\cos\left({{\omega}_{d}}t+{{\phi}_{d}}\right)$,
where $\omega_{d}$ and $\phi_{d}$ indicate the frequency and initial phase of
the driving field. Hence, the control Hamiltonian is expressed as

${{H}_{C}}\left(t\right)=u\left(t\right)\left({{e}^{i{\phi}}}\left|a\right\rangle\left\langle
e\right|+{{e}^{-i{{\phi}}}}\left|e\right\rangle\left\langle
a\right|\right)\cos\left({{\omega}}t\right)$ (7)

Here, decoherence channels are considered as non-unital, meaning that

$\mathcal{C}\left(\rho\right)=\sqrt{\sum\limits_{i=1}^{2}{{{\left(tr\left({{M}_{i}}\rho\right)\right)}^{2}}}}$
(8)

where ${M_{1}}=\left|e\right\rangle\left\langle
a\right|+\left|a\right\rangle\left\langle e\right|$ and
${M_{2}}=-i\left(\left|e\right\rangle\left\langle
a\right|-\left|a\right\rangle\left\langle e\right|\right)$. In order to
preserve the coherence during the whole evolution, we impose the following
constraint:

$\alpha{{c}_{0}}\leq\mathcal{C}\left(\rho\right)\leq{{c}_{0}},\quad
0<\alpha<1$ (9)

In practice, any experimental quantum state preparation or transformation is
exposed to environmental noise and imperfections. Hence, one needs to quantify
the degree of closeness for a pair of quantum states, i.e. the degree of
accuracy of an approximate calculation for replication of some intended target
density operator. Given two density operators $\rho$ and $\sigma$, the
fidelity is expressed as, [24],

${\mathcal{F}(\rho,\sigma)=\left(\operatorname{tr}{\sqrt{{\sqrt{\rho}}\sigma{\sqrt{\rho}}}}\right)^{2}}$
(10)

An equivalent expression for the fidelity may be written using the trace norm
$\mathcal{F}(\rho,\sigma)={\Big{(}}\operatorname{tr}|{\sqrt{\rho}}{\sqrt{\sigma}}|{\Big{)}}^{2}$
where the absolute value of an operator is here defined as
${\displaystyle|A|\equiv{\sqrt{A^{\dagger}A}}}$.

##  IV Optimal Control Formulation

Consider the following optimal control problem $(P_{1})$

$\begin{split}\text{Minimize}\quad&-\mathcal{F}(\rho,\sigma)\\\ \text{subject
to}\quad&\dot{\rho}(t)=F\left(\rho(t),u(t)\right)\\\
&\rho\left(0\right)={{\rho}_{0}}\in{{\mathbb{R}}^{n}}\times{{\mathbb{R}}^{n}}\\\
&u\left(t\right)\in\mathcal{U}:=\left\\{u\in{{L}_{\infty}}:u\left(t\right)\in\Omega\subset{{\mathbb{R}}^{m}}\right\\}\\\
&h\left(\rho\left(t\right)\right)\leq 0\quad\forall t\in\left[0,1\right]\\\
\end{split}$

where $\dot{\rho}:=\displaystyle\frac{\partial\rho}{dt}$,
$t\in\left[0,1\right]$ determines the time variable, $\rho$ is the state
variable in ${{\mathbb{R}}^{n}}\times{{\mathbb{R}}^{n}}$ supposed to satisfy
the differential constraints $\dot{\rho}=F\left(\rho,u\right)$, and $\rho_{0}$
is the so-called initial quantum state. $u(\cdot)$ is the measurable bounded
function termed as control. The inequality
$h\left(\rho\left(t\right)\right)\leq 0$ $\forall t\in\left[0,1\right]$
defines the state constrains for the density operator. Here, we aim to
maximize fidelity so that the density operator $\rho$ has the maximum overlap
with the target $\sigma$. The fidelity criteria signifies a security level for
quantum state transformation or the effectiveness of quantum gate synthesis.

Let $\left({{\rho}^{*}},{{u}^{*}}\right)$ be the solution to the problem
$(P_{1})$. The pair $\left({{\rho}^{*}},{{u}^{*}}\right)$ is termed as optimal
if the value of fidelity is the maximum possible over the set of all feasible
solutions. Then, there exist Lagrangian multipliers: $p$, $\mu$, and
$\lambda$, where $p$ is an absolutely continuous function
$p:\left[0,1\right]\to{{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}}$, $\mu$ is
bounded variation, non-increasing $\mu:\left[0,1\right]\to{{\mathbb{R}}^{k}}$,
such that $\mu(t)=cte$ on the time interval for which the state constraint is
non-active, i.e.,
$\left\\{t\in\left[0,1\right],h\left(\rho\left(t\right)\right)<0\right\\}$,
and $\lambda\geq 0$ is non-negative. We also consider the following
assumptions:

  * •

$h\left(\rho\left(t\right)\right)$ is $C^{2}$.

  * •

$\left({{\rho}^{*}},{{u}^{*}}\right)$ is a regular control process.

  * •

The multiplier $\mu$ is continuous.

Since Gamkrelidze set of necessary conditions is used here, [1, 34], the
second-order derivative for $h\left(\rho\left(t\right)\right)$, which
indicates state constraints, is required, [11].

##  V Maximum Principle

For the indicated set of conditions, the extended Hamilton-Pontryagin function
can be written as

$\mathcal{H}\left(\rho,u,\pi,\mu\right)=tr\left({{\left(\pi-\mu{{\nabla}_{\rho}}h\right)}^{\dagger}}F\left(\rho,u\right)\right)$
(11)

in which $tr$ indicates the trace of matrix. According to Pontryagin’s maximum
principle, the optimal state trajectory $\rho^{*}$, optimal control $u^{*}$,
and corresponding Lagrange multiplier matrix $\pi^{*}$ must maximize the
Hamiltonian $\mathcal{H}$ so that

$\mathcal{H}\left({\rho}^{*}(t),u,{\pi}(t),\mu(t)\right)\leq\mathcal{H}\left({\rho}^{*}(t),{u}^{*}(t),{\pi}(t),\mu(t)\right)$
(12)

for all time $t\in[0,T]$ and for all permissible control inputs
$u\in\mathcal{U}$.

Additionally, the adjoint equation, and its terminal conditions imply that,
Lebesgue a.e.,

$(-\dot{\pi}^{\dagger}(t),\dot{\rho}^{*}(t))=\nabla_{(\rho,\pi)}{\cal
H}(\rho^{*}(t),u^{*}(t),\pi(t),\mu(t))$ (13)
$\pi^{\dagger}(1)=\nabla_{\rho}\mathcal{F}(\rho^{*}(1),\sigma(1))$ (14)

To satisfy the state constraint indicated in (9), we defined
$h\left(\rho\left(t\right)\right)$ as the following

$h\left(\rho\right)=\left(\begin{matrix}{{h}_{1}}\left(\rho\right)\\\
{{h}_{{}_{2}}}\left(\rho\right)\\\
\end{matrix}\right)=\left(\begin{matrix}\mathcal{C}\left(\rho\right)-{{c}_{0}}\\\
\alpha{{c}_{0}}-\mathcal{C}\left(\rho\right)\\\ \end{matrix}\right)$ (15)

For the sake of simplicity, we define
$\bar{\mathcal{C}}\left(\rho\right)={{\mathcal{C}}^{2}}\left(\rho\right)$, and
rewrite the constraints as

$\bar{h}\left(\rho\right)=\left(\begin{matrix}{\bar{h}_{1}}\left(\rho\right)\\\
{\bar{h}_{{2}}}\left(\rho\right)\\\
\end{matrix}\right)=\left(\begin{matrix}\bar{\mathcal{C}}\left(\rho\right)-c_{0}^{2}\\\
{{\alpha}^{2}}c_{0}^{2}-\bar{\mathcal{C}}\left(\rho\right)\\\
\end{matrix}\right)\leq 0$ (16)

Hence,

${{\nabla}_{\rho}}\bar{h}\left(\rho\right)=\left(\begin{matrix}{{\nabla}_{\rho}}\bar{\mathcal{C}}\left(\rho\right)\\\
-{{\nabla}_{\rho}}\bar{\mathcal{C}}\left(\rho\right)\\\ \end{matrix}\right)$
(17)

in which
${{\nabla}_{\rho}}\bar{C}\left(\rho\right)=2\sum\limits_{i=1}^{2}{\left(tr\left({{M}_{i}}\rho\right)\right){{M}_{i}}}$.

Noting that $\mu=\left(\begin{matrix}{{\mu}_{1}}&{{\mu}_{2}}\\\
\end{matrix}\right)$, we can rewrite (11) as

$\\!\\!\\!\\!\mathcal{H}\\!=\\!tr\left(\\!{{\left(\\!\pi\\!-\\!\bar{\mu}\sum\limits_{i=1}^{2}{tr\left({{M}_{i}}\rho\\!\right)\\!{{M}_{i}}}\\!\right)}^{\dagger}}\\!\\!\left(\\!-i\left[H,\rho\right]\\!+\\!\mathcal{D}\left(\rho\right)\\!\right)\\!\right)$
(18)

where $\bar{\mu}=2(\mu_{1}-\mu_{2})$.

From (13) we have

$-{{\dot{\pi}}^{\dagger}}\\!=\\!\mathcal{G}\left({{\pi}^{\dagger}}\right)\\!-\\!\bar{\mu}\sum\limits_{i=1}^{2}{\\!\\!{{M}_{i}}tr\\!\left(\mathcal{G}\left(M_{i}^{\dagger}\right)\rho\right)\\!+\\!tr\left({{M}_{i}}\rho\right)\\!\mathcal{G}\left(M_{i}^{\dagger}\\!\right)\\!}$
(19)

where

$\mathcal{G}\left(\mathcal{X}\right)=-i\left[\mathcal{X},H\left(t\right)\right]+\mathcal{D}\left(\mathcal{X}\right)+{{\mathcal{D}}_{0}}\left(\mathcal{X}\right)$
(20)

in which

${{\mathcal{D}}_{0}}\left(\mathcal{X}\right)=\sum\limits_{k}{{{\gamma}_{k}}\left(L_{k}^{\dagger}\mathcal{X}{{L}_{k}}-{{L}_{k}}\mathcal{X}L_{k}^{\dagger}\right)}$
(21)

Hence, by vectorization we have

$-\frac{d}{dt}vec\left({{\pi}^{\dagger}}\right)={\hat{\mathcal{L}}}^{\prime}vec\left({{\pi}^{\dagger}}\right)+vec\left(cte\left(\rho(t),\bar{\mu}(t)\right)\right)$
(22)

where

${\hat{\mathcal{L}}}^{\prime}=\hat{\mathcal{L}}+\sum\limits_{k}{{\gamma}_{k}}\left(\left(L_{k}^{T}\otimes
L_{k}^{\dagger}\right)+\left(L_{k}^{*}\otimes{{L}_{k}}\right)\right)$ (23)

The boundary condition (14) at the final time is obtained by computing, [35],

$\pi^{\dagger}(1)=\nabla_{\rho}\left(tr\sqrt{\sqrt{\rho(1)}\sigma(1)\sqrt{\rho(1)}}\right)^{2}$
(24)

for which by means of Taylor expansion at the point $\rho=I$, we have

$\sqrt{\rho}=\sum_{k=0}^{\infty}\frac{1}{k!}\frac{d^{k}}{d\rho^{k}}\sqrt{\rho}_{|\rho=I}(\rho-I)^{k}$
(25)

which, in turn, under some conditions, the application of the Cayley–Hamilton
theorem yields, for a certain coefficients $\alpha_{k}$, $k=0,\ldots,n-1$,
being $n$ the dimension of the square matrix $\rho$,

$\sqrt{\rho}=\sum_{k=0}^{n-1}\alpha_{k}(\rho-I)^{k}.$

Thus, by using matricial calculus, we obtain

$\displaystyle\pi^{\dagger}(1)=2tr\sqrt{\rho(1)\sigma(1)}\nabla_{\rho}(tr\sqrt{\rho(1)\sigma(1)})$
(26)
$\displaystyle=2tr\sqrt{\rho(1)\sigma(1)}\sum\limits_{k=0}^{n-1}{\alpha}_{k}\sum_{i=0}^{k-1}\bar{\rho}(1)^{i}\sqrt{\sigma(1)}\bar{\rho}(1)^{k-i-1}$

where $\bar{\rho}(1)=\rho(1)-I$.

##  VI Application of the Pontryagin Maximum Principle

Let denote the excited state as

$\rho_{e}=\left(\begin{matrix}1&0&0\\\ 0&0&0\\\ 0&0&0\\\ \end{matrix}\right)$
(27)

and the ground states as

$\rho_{a}=\left(\begin{matrix}0&0&0\\\ 0&1&0\\\ 0&0&0\\\
\end{matrix}\right),\quad\rho_{b}=\left(\begin{matrix}0&0&0\\\ 0&0&0\\\
0&0&1\\\ \end{matrix}\right)$ (28)

The drift Hamiltonian for such system is expressed as

$\displaystyle{{H}_{0}}$
$\displaystyle=\sum\limits_{j=1}^{3}{{{E}_{j}}{{\rho}_{j}}}=0.8{{\rho}_{e}}+0.5{{\rho}_{a}}+0.4{{\rho}_{b}}$
(29) $\displaystyle=\left(\begin{matrix}0.8&0&0\\\ 0&0.5&0\\\ 0&0&0.4\\\
\end{matrix}\right)$

and the control Hamiltonian is

${{H}_{C}}\left(t\right)=u\left(t\right)\left(\begin{matrix}0&{{e}^{-i{\phi}}}&0\\\
{{e}^{i{{\phi}}}}&0&0\\\ 0&0&0\\\
\end{matrix}\right)\cos\left({{\omega}}t\right)$ (30)

From (3), we have

$\displaystyle\mathcal{D}\left(\rho\left(t\right)\right)\\!=\\!10^{-1}\\!\left(\\!{{L}_{1}}\rho\left(t\right)L_{1}^{\dagger}-\frac{1}{2}L_{1}^{\dagger}{{L}_{1}}\rho\left(t\right)-\frac{1}{2}\rho\left(t\right)L_{1}^{\dagger}{{L}_{1}}\\!\right)$
(31)
$\displaystyle+10^{-3}\left({{L}_{2}}\rho\left(t\right)L_{2}^{\dagger}-\frac{1}{2}L_{2}^{\dagger}{{L}_{2}}\rho\left(t\right)-\frac{1}{2}\rho\left(t\right)L_{2}^{\dagger}{{L}_{2}}\right)$

in which

${{L}_{1}}=\left(\begin{matrix}0&0&0\\\ 1&0&0\\\ 0&0&0\\\
\end{matrix}\right),\quad{{L}_{2}}=\left(\begin{matrix}0&0&0\\\ 0&0&0\\\
1&0&0\\\ \end{matrix}\right)\cdot$

In this example, the initial state reflects the fact that the initial
population distribution is mostly on the levels $\rho_{a}$ and $\rho_{e}$.
Hence, as it can be seen from (31), the relaxation rate $\gamma_{2}$ is
considerably smaller in comparison with $\gamma_{1}$, so the three-level
system can be considered as a two-level one. The coherence of the system is
expressed by

$\bar{\mathcal{C}}\left(\rho\right)={{\left(tr\left({{M}_{1}}\rho\right)\right)}^{2}}+{{\left(tr\left({{M}_{2}}\rho\right)\right)}^{2}}$
(32)

where ${{M}_{1}}=\left(\begin{matrix}0&1&0\\\ 1&0&0\\\ 0&0&0\\\
\end{matrix}\right)$ and ${{M}_{2}}=\left(\begin{matrix}0&-i&0\\\ i&0&0\\\
0&0&0\\\ \end{matrix}\right)$. Suppose that $\rho_{0}=\rho_{e}$ is the initial
state, meaning that the system remains on the excited state. Since we consider
the effectiveness of our method for any pure state, we present the desired
target state by

$\sigma\left(1\right)=\left(\begin{matrix}0&0&0\\\
0&{{\sin}^{2}}\beta&\sin\beta\cos\beta\\\
0&\sin\beta\cos\beta&{{\cos}^{2}}\beta\\\ \end{matrix}\right)\ $ (33)

for $\beta\in[0,2\pi)$.

Now, we present a time discretized computational scheme to solve the optimal
control problem associated with $(P_{1})$ by using an indirect method based on
the Maximum Principle. Let $N$ be the number of discrete time subintervals.
Given the smooth properties of the problems data, let us consider a uniform
discretization. Thus, we consider the $N$ points $\displaystyle
t_{m}=\frac{m}{N}$ for $m=0,\ldots,N-1$ and denote the value of any function
$f(t_{m})$ by $f_{m}$. We consider both the system dynamics, and the adjoint
differential equation approximated by a first order Euler approximation.
Higher order methods yield better approximations but our option is made to
keep the presentation simple. Let $j=0,\ldots$ be the iterations counter, and
we denote the $j^{th}$ iteration of the function $f$ at time $t_{m}$ by
$f_{m}^{j}$.
From now on the control is considered as
$V_{m}^{j}=\left(u_{m}^{j},{{\phi}^{j}},{{\omega}^{j}}\right)$, in which $u$
depends on time, however, $\phi$ and $\omega$ modify in each iteration.

The proposed algorithm is as follows:

  * Step 1 - Initialization.

Let $j=0$, for $m=1,\ldots,N-1$ in (30),

    * Initialize $V_{m}^{j}=\left(u_{m}^{j},{{\phi}^{j}},{{\omega}^{j}}\right)$

  * Step 2 - Computation of the state trajectory.
For $m=0,\dots,N$

    * Vectorize $\rho_{0}^{j}$.

    * Compute

$vec\left(\rho_{m}^{j}\right)={{e}^{\sum\limits_{k=0}^{m-1}{\frac{1}{N}\hat{\mathcal{L}}\left(V_{k}^{j}\right)}}}vec\left(\rho_{0}^{j}\right),$

    * and obtain

$\rho_{m}^{j}=ve{{c}^{-1}}\left(\rho_{m}^{j}\right)\cdot$

  * Step 3 - Computation of the adjoint trajectory.
Compute $\pi_{N}^{j}$ by using (26) with $\rho_{N}^{j}$ computed in Step 2.
Initialize $\mu^{j}_{1,N}=0$ and $\mu^{j}_{2,N}=0$.
Obtain $vec(\pi_{N}^{j})$ and $vec(cte(\rho_{m}^{j},\bar{\mu}_{m}^{j}))$.

For $m=N-1$ to $0$,

    * Compute

$\displaystyle vec\left(\pi_{m}^{j}\right)$
$\displaystyle={{e}^{\sum\limits_{k=m}^{N-1}{\frac{1}{N}{\hat{\mathcal{L}}}^{\prime}\left(V_{k}^{j}\right)}}}vec\left(\pi_{N}^{j}\right)$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!+\sum\limits_{k=m}^{N-1}{\frac{1}{N}{{e}^{\sum\limits_{l=k}^{N-1}{\frac{1}{N}{\hat{\mathcal{L}}}^{\prime}\left(V_{l}^{j}\right)}}}vec\left(cte\left(\rho_{m}^{j},\bar{\mu}_{m}^{j}\right)\right)}$

    * Obtain

$\pi_{m}^{j}=ve{{c}^{-1}}\left(vec\left(\pi_{m}^{j}\right)\right)$

    * Obtain

$\bar{\mu}_{m-1}^{j}=2\left(\mu_{1,m-1}^{j}-\mu_{2,m-1}^{j}\right)$

by checking

      * *

If $\mathcal{\bar{C}}(\rho_{m}^{j})<c_{0}^{2}$, then

$\mu^{j}_{1,m-1}=\mu^{j}_{1,m}$

and if $\mathcal{\bar{C}}(\rho_{m}^{j})>\alpha^{2}c_{0}^{2}$, then

$\mu^{j}_{2,m-1}=\mu^{j}_{2,m}$

or

      * *

If $\mathcal{\bar{C}}(\rho_{m}^{j})=c_{0}^{2}$, then

$\mu^{j}_{1,m-1}=\mu^{j}_{1,m}+\frac{1}{N}{{\nabla}_{\rho}}\bar{\mathcal{C}}\left(\rho_{m}^{j}\right)$

and if $\mathcal{\bar{C}}(\rho_{m}^{j})=\alpha^{2}c_{0}^{2}$, then

$\mu^{j}_{2,m-1}=\mu^{j}_{2,m}-\frac{1}{N}{{\nabla}_{\rho}}\bar{C}\left(\rho_{m}^{j}\right)$

  * Step 4 - Computation of the Pontryagin Hamilton function
For $m=0,\dots,N-1$, let $\mathcal{H}_{m}^{j}(u,\omega,\phi)$ as

$\\!\\!\\!\\!\\!tr\left(\\!{{\left(\\!\pi^{j}_{m}\\!-\\!{\bar{\mu}^{j}_{m}}\sum\limits_{i=1}^{2}{tr\left({{M}_{i}}\rho\right){{M}_{i}}}\\!\right)}^{\dagger}}\\!\\!\left(\\!-i\left[H^{j}_{m},\rho^{j}_{m}\right]\\!+\\!\mathcal{D}\left(\rho\right)\\!\right)\\!\right)$

where $\displaystyle H(u,\omega,\phi,t_{m})$, from (29) and (30), is given by

$\left(\begin{matrix}0.8&u{{e}^{-i{\phi}}}cos(\omega t_{m})&0\\\
u{{e}^{i{{\phi}}}}cos(\omega t_{m})&0.5&0\\\ 0&0&0.4\\\ \end{matrix}\right)$

  * Step 5: Update the control function.
For $m=0,\ldots,N-1$, compute the values $u_{m}^{j+1}$, $\omega^{j+1}$, and
$\phi^{j+1}$ that maximize the map

$(u,\omega,\phi)\to\mathcal{H}_{m}^{j}(u,\omega,\phi).$

  * Step 6: Stopping test.

For given small positive numbers $\varepsilon_{\omega}$, $\varepsilon_{\phi}$,
and $\varepsilon_{u}$ (tolerance errors), check whether all the following
inequalities hold:

$\displaystyle\hskip 42.67912pt|\omega^{j}-\omega^{j+1}|<\varepsilon_{\omega}$
$\displaystyle\hskip 42.67912pt|\phi^{j}-\phi^{j+1}|<\varepsilon_{\phi}$
$\displaystyle\hskip
14.22636pt\max_{m=0,\ldots,N-1}\\{|u_{m}^{j}-u_{m}^{j+1}|\\}<\varepsilon_{u}$

If yes, let

$\omega^{*}=\omega^{j+1}$ $\phi^{*}=\phi^{j+1}$

and for $m=0,\ldots,N-1$,

$u^{*}(t_{m})=u_{m}^{j+1}$

and exit the algorithm.

Otherwise, let $j=j+1$, go to Step 2.

##  VII Conclusion

In this article, we have shown the application of Maximum Principle of
Pontryagin in Gamkrelidze’s form in maximizing fidelity in a quantum state
transformation process controlled in such a way as to preserve coherence
within given bounds while, simultaneously, maximizing the fidelity in the
transfer of the state variable. The context of the matrix-valued probability
density function, between two given state values with the dynamics of the
system expressed by the Lindblad master equation was adapted. We introduced a
complete set of equations enabling the computation of optimal control strategy
in the presence of state constraints, which is a novel approach in the context
of quantum systems. As an example, we took into account a three-level
$\Lambda$-type atom subjected to Markovian decoherence characterized by non-
unital decoherence channels. We finally introduced an algorithm to solve the
proposed problem.

##  VIII ACKNOWLEDGMENTS

The authors acknowledge the support of FCT for the grant 2021.07608.BD, the
ARISE AL, Ref. LA/P/0112/2020, SYSTEC-Base, Ref. UIDB/00147/2020, and
Programmatic, Ref. UIDP/00147/2020, and also the support of projects SNAP,
Ref. NORTE-01-0145-FEDER-000085, and MAGIC, Ref. PTDC/EEI-AUT/32485/2017,
COMPETE2020-POCI and FCT/MCTES (PIDDAC).

## References

  * [1] A.V. Arutyunov, D.Y. Karamzin, and F. Lobo Pereira ”The maximum principle for optimal control problems with state constraints by RV Gamkrelidze: revisited”, Journal of Optimization Theory and Applications, vol. 149, 2011, pp 474-493.
  * [2] S.J. Glaser, U. Boscain, T. Calarco, C.P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, and T. Schulte-Herbrüggen, and others, ”Training Schrödinger’s cat: quantum optimal control”, The European Physical Journal D, vol. 69, 2015, pp 1-24.
  * [3] D. d’Alessandro, ”Introduction to quantum control and dynamics”, Chapman and hall/CRC, 2021.
  * [4] C. Brif, R. Chakrabarti, and H. Rabtiz, ”Control of quantum phenomena: past, present and future”, New Journal of Physics, vol. 12, 2010, pp 075008.
  * [5] D. Dong and I.R. Petersen, ”Quantum control theory and applications: a survey”, IET Control Theory & Applications, vol. 4, 2010, pp 2651-2671.
  * [6] P. Rembold, N. Oshnik, M.M Müller, S. Montangero, T. Calarco, and E. Neu, ”Introduction to quantum optimal control for quantum sensing with nitrogen-vacancy centers in diamond”, AVS Quantum Science, vol. 2, 2020, pp 024701.
  * [7] T. Schaetz, C.R. Monroe, and T. Esslinger, ”Focus on quantum simulation”, New Journal of Physics, vol. 15, 2013, pp 085009.
  * [8] N. Khaneja, R. Brockett, and S.J. Glaser, ”Time optimal control in spin systems”, Physical Review A, vol. 63, 2001, pp 032308.
  * [9] S. Grivopoulos and B. Bamieh, ”Optimal population transfers for a quantum system in the limit of large transfer time”, Proceedings of the 2004 American Control Conference, vol. 3, 2004, pp 2481-2486.
  * [10] A.P. Peirce, M.A. Dahleh, and H. Rabitz, ”Optimal control of quantum-mechanical systems: Existence, numerical approximation, and applications”, Physical Review A, vol. 37, 1988, pp 4950.
  * [11] A. Arutyunov, D. Karamzin, and F. Lobo Pereira, ”A remark on the continuity of the measure Lagrange multiplier in state constrained optimal control problems”, 2018 IEEE Conference on Decision and Control (CDC), 2018, pp 49-54.
  * [12] P. Facchi, S. Tasaki, S. Pascazio, and H. Nakazato, A. Tokuse, and DA. Lidar, ”Control of decoherence: Analysis and comparison of three different strategies”, Physical Review A, vol. 71, 2005, pp 022302.
  * [13] Q. Ai, Y. Li, and G. Long, ”Influences of gate operation errors in the quantum counting algorithm”, Journal of Computer Science and Technology, vol. 21, 2006, pp 927-932.
  * [14] J. Niwa, K. Matsumoto, and H. Imai, ”General-purpose parallel simulator for quantum computing”, International Conference on Unconventional Methods of Computation, 2002, pp 230-251.
  * [15] E. Joos, and H.D. Zeh, C. Kiefer, D.JW. Giulini, J. Kupsch, and I. Stamatescu, ”Decoherence and the appearance of a classical world in quantum theory”, Springer Science & Business Media, 2013.
  * [16] P.W. Shor, ”Scheme for reducing decoherence in quantum computer memory”, Physical review A, vol. 52, 1995, pp R2493.
  * [17] L.M. Duan and G.C Guo, ”Prevention of dissipation with two particles”, Physical review A, vol. 57, 1998, pp 2399.
  * [18] P. Cappellaro, J. Hodges, T. Havel, and D. Cory, ”Principles of control for decoherence-free subsystems”, The Journal of chemical physics, vol. 125, 2006, pp 044514.
  * [19] S. Tasaki, A. Tokuse, P. Facchi, and S. Pascazio, ”Control of decoherence: Dynamical decoupling versus quantum Zeno effect: A case study for trapped ions”, International journal of quantum chemistry, vol. 98, 2004, pp 160-172.
  * [20] D.A. Lidar, S. Schneider, ”Stabilizing qubit coherence via tracking-control”, arXiv preprint quant-ph/0410048, 2004.
  * [21] J. Zhang, R. Wu, C. Li, and T. Tarn, ”Protecting coherence and entanglement by quantum feedback controls”, IEEE Transactions on Automatic Control, vol. 55, 2010, pp 619-633.
  * [22] D.G. Cory, A.F. Fahmy, and T.F. Havel, ”Ensemble quantum computing by NMR spectroscopy”, Proceedings of the National Academy of Sciences, vol. 94, 1997, pp 1634-1639.
  * [23] W.S. Warren, H. Rabitz, and M. Dahleh, ”Coherent control of quantum dynamics: the dream is alive”, Science, vol. 259, 1993, pp 1581-1589.
  * [24] Y. Liang, Y. Yeh, P.E. Mendonça, R.Y. Teh, M.D. Reid, and P.D. Drummond, ”Quantum fidelity measures for mixed states”, Reports on Progress in Physics, vol. 82, 2019, pp 076001.
  * [25] A. Gilchrist, N.K. Langford, and M.A. Nielsen, ”Distance measures to compare real and ideal quantum processes”, Physical Review A, vol. 71, 2005, pp 062310.
  * [26] T. Chanelière, G. Hétet, and N. Sangouard, ”Quantum optical memory protocols in atomic ensembles”, Advances In Atomic, Molecular, and Optical Physics, vol. 67, 2018, pp 77-150.
  * [27] M. Hillery, R.F. O’Connell, M.O. Scully, and E.P. Wigner, ”Distribution functions in physics: Fundamentals”, Physics reports, vol. 106, 1984, pp 121-167.
  * [28] S. Popescu, ”Bell’s inequalities versus teleportation: What is nonlocality?”, Physical review letters, vol. 72, 1994, pp 797.
  * [29] H.P. Breuer, F. Petruccione, and others, ”The theory of open quantum systems”, Oxford University Press on Demand, 2002.
  * [30] E.B. Davies, ”Markovian master equations. II”, Mathematische Annalen, vol. 219, 1976, pp 147-158.
  * [31] G. Lindblad, ”On the generators of quantum dynamical semigroups”, Communications in Mathematical Physics, vol. 48, 1976, pp 119-130.
  * [32] R. Alicki and K. Lendi, ”Quantum dynamical semigroups and applications”, Springer, vol. 717, 2007.
  * [33] M.O. Scully and M.S. Zubairy, ”Quantum optics”, American Association of Physics Teachers, 1999.
  * [34] A.V. Arutyunov, D.Y. Karamzin, and F. Lobo Pereira, ”Investigation of controllability and regularity conditions for state constrained problems”, IFAC-PapersOnLine, vol. 50, 2017, pp 6295-6302.
  * [35] N.B. Dehaghani and F. Lobo Pereira, ”High Fidelity Quantum State Transfer by Pontryagin Maximum Principle”, arXiv preprint arXiv:2203.04361, 2022.


# Finite Model Properties for Residuated Semigroups

Szabolcs Mikulás<EMAIL_ADDRESS>
(Date: 21/11/2022)

###### Abstract.

We have a quick look at various finite model properties for residuated
semigroups. In particular, we solve Problem 19.17 from the monograph by Hirsch
and Hodkinson [HH02].

##  1\. Introduction

Residuated semigroups have been investigated on their own right, but also in
connection with susbstructural logics, like the Lambek Calculus (LC) Lambek
[La58], and relevance logics Anderson et al. [ABD92]. Indeed, relational
semantics, i.e. models based on binary relations, have been proposed for these
substructural logics, see van Benthem [vB91] and Routley and Meyer [RM73], and
the algebraic settings of these relational semantics give rise to residuated
semigroups (and their expansions to larger similarity types). For instance,
representable residuated semigroups (see below) provide sound and complete
semantics for LC, see Andréka and Mikulás [AM94]. The connection between
(fragments of) relevance logics and (variants of) residuated semigroups have
been investigated in e.g. Bimbó et al. [BDM09], Maddux [Ma10], Hirsch and
Mikulás [HM11] and Mikulás [Mi09],[Mi15].

In this note, we will concentrate on representable residuated semigroups with
finite bases (see below) and show that they have smaller expressive power than
representable semigroups in general. In particular, we show that there is a
finite representable residuated semigroup that is not isomorphic to any
representable residuated semigroup with a finite base, thus providing a
solution to Problem 19.17 of Hirsch and Hodkinson [HH02].

Next we define (representable) residuated semigroups. Strictly speaking, these
are not algebras but algebraic structures (because of the ordering $\leq$) but
the reader should not have any problem with applying the usual algebraic
techniques (like taking subalgebras) to algebraic structures.

######  Definition 1.1.

A _residuated semigroup_
${\mathcal{A}}=(A,\leq,\circ,\mathbin{\backslash},\mathbin{/})$ is an algebra
such that the following are satisfied:

  * •

$\leq$ is a reflexive ordering on $A$,

  * •

$(A,\circ)$ is a semigroup,

  * •

$\circ$ is monotonic w.r.t. $\leq$,

  * •

$\mathbin{\backslash},\mathbin{/}$ are the right and left residuals of
$\circ$:

(1) $y\leq x\mathbin{\backslash}z\text{ iff }x\circ y\leq z\text{ iff }x\leq
z\mathbin{/}y$

for every $x,y,z$ in $A$.

We denote the class of residuated semigroups by ${\sf RS}$.

Next we define subclasses of ${\sf RS}$ that are isomorphic to algebras of
binary relation and the relation $\leq$ and operations
$\circ,\mathbin{\backslash},\mathbin{/}$ have “relational” interpretations.

######  Definition 1.2.

A _representable residuated semigroup_
${\mathcal{A}}=(A,\leq,\circ,\mathbin{\backslash},\mathbin{/})$ is a
subalgebra of some
$(\wp(W),\subseteq,\circ_{W},\mathbin{\backslash}_{W},\mathbin{/}_{W})$ where
$W$ is a transitive relation, $\subseteq$ is the subset relation and the
operations are interpreted as follows:

(2) $\displaystyle x\circ_{W}y$ $\displaystyle=\\{(u,v)\in W:(u,w)\in x\mbox{
and }(w,v)\in y\mbox{ for some }w\\},$ (3) $\displaystyle
x\mathbin{\backslash}_{W}y$ $\displaystyle=\\{(u,v)\in W:\text{for every
}w,(w,u)\in x\text{ implies }(w,v)\in y\\},$ (4) $\displaystyle
x\mathbin{/}_{W}y$ $\displaystyle=\\{(u,v)\in W:\text{for every }w,(v,w)\in
y\text{ implies }(u,w)\in x\\},$

see Figure 1.

$\textstyle{{\exists
w}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{y}$$\textstyle{{\forall
w}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{x}$$\scriptstyle{y}$$\textstyle{\forall
w}$$\textstyle{\land}$$\textstyle{\rightarrow}$$\textstyle{\leftarrow}$$\textstyle{{u}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{x}$$\scriptstyle{x\circ_{W}y}$$\textstyle{v}$$\textstyle{{u}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{x\mathbin{\backslash}_{W}y}$$\textstyle{v}$$\textstyle{{u}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{x\mathbin{/}_{W}y}$$\scriptstyle{x}$$\textstyle{{v}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$$\scriptstyle{y}$
Figure 1. Interpretations of composition and its residuals

Let $U$ be the smallest set such that $W\subseteq U\times U$. We call $U$ the
_base_ and $W$ the _unit_ of $\mathcal{A}$. We denote the class (of isomorphs)
of representable residuated semigroups by ${\sf RRS}$.

Let ${\sf RRS}^{\Box}$ be the “square” subclass of ${\sf RRS}$, where we
additionally require that $W=U\times U$.

The subscript $W$ to the operations of representable algebras is to remind the
reader that their meanings depend on the choice of the unit $W$ of the
algebra. In particular, if the unit $W$ is an irreflexive relation, then
$(u,u)\notin x\backslash_{W}x$ for any element $u$ in the base and algebra
element $x$, in stark contrast to $(v,v)\in y\backslash_{V}y$ for any base
element $v$ and algebra element $y$ for a reflexive unit $V$.

######  Definition 1.3.

let $\sf K$ be a class of representable algebras (like ${\sf RRS}$ or ${\sf
RRS}^{\Box}$).

We say that ${\mathcal{A}}\in{\sf K}$ is _representable on a finite base_ , or
_finitely representable_ , if the base $U$ can be chosen to be finite in
Definition 1.2.

The class $\sf K$ has the _finite representation property_ (frp) if every
finite ${\mathcal{A}}\in{\sf K}$ is finitely representable.

##  2\. The main result

Next we show that ${\sf RRS}$ lacks the frp.

######  Theorem 2.1.

There is a finite ${\mathcal{A}}\in{\sf RRS}$ that cannot be represented on a
finite base.

###### Proof.

Let $U=({\mathbb{Q}}\times\\{0\\})\cup({\mathbb{Q}}\times\\{1\\})$ and

$W=\\{((q,0),(r,0)):q,r\in{\mathbb{Q}},q<r\\}\cup\\{((q,1),(r,0)):q,r\in{\mathbb{Q}},q\leq
r\\}$

where $\mathbb{Q}$ denotes the set of rational numbers. Define

(5) $\displaystyle a$
$\displaystyle=\\{((q,0),(r,0)):q,r\in{\mathbb{Q}},q<r\\},$ (6) $\displaystyle
b$ $\displaystyle=\\{((q,1),(q,0)):q\in{\mathbb{Q}}\\}.$

Let ${\mathcal{A}}=(A,\leq,\circ_{W},\backslash_{W},/_{W})$ be the subalgebra
of $(\wp(W),\subseteq,\circ_{W},\backslash_{W},/_{W})$ generated by $a,b$.
$\mathcal{A}$ has the following seven elements: bottom element $\emptyset$
(that we will denote by $\bot$), top element $W$ (that we will denote by
$\top$), $a$, $b$,

$b\circ_{W}a=\\{((q,1),(r,0)):q,r\in{\mathbb{Q}},q<r\\},$

$a\cup b\circ_{W}a$ (that we will denote by $a^{\prime}$) and $b\cup
b\circ_{W}a$ (that we will denote by $b^{\prime}$).111Note that $\cup$ is not
an operation of $\mathcal{A}$, but it makes sense to take the union of two
elements of an algebra of binary relations. We include the tables for the
operations (where we abbreviate $b\circ_{W}a$ as $ba$).

$\circ_{W}$ | $a$ | $b$ | $\bot$ | $ba$ | $\top$ | $b^{\prime}$ | $a^{\prime}$
---|---|---|---|---|---|---|---
$a$ | $a$ | $ba$ | $\bot$ | $ba$ | $a^{\prime}$ | $ba$ | $a^{\prime}$
$b$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$
$\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$
$ba$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$
$\top$ | $a$ | $ba$ | $\bot$ | $ba$ | $a^{\prime}$ | $ba$ | $a^{\prime}$
$b^{\prime}$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$ | $\bot$
$a^{\prime}$ | $a$ | $ba$ | $\bot$ | $ba$ | $a^{\prime}$ | $ba$ | $a^{\prime}$
Table 1. Table for composition

$\mathbin{\backslash}_{W}$ | $a$ | $b$ | $\bot$ | $ba$ | $\top$ | $b^{\prime}$ | $a^{\prime}$
---|---|---|---|---|---|---|---
$a$ | $\top$ | $b^{\prime}$ | $\top$ | $b^{\prime}$ | $b^{\prime}$ | $b^{\prime}$ | $b^{\prime}$
$b$ | $b^{\prime}$ | $b^{\prime}$ | $\top$ | $b^{\prime}$ | $b^{\prime}$ | $b^{\prime}$ | $b^{\prime}$
$\bot$ | $b^{\prime}$ | $b^{\prime}$ | $\top$ | $b^{\prime}$ | $b^{\prime}$ | $b^{\prime}$ | $b^{\prime}$
$ba$ | $b^{\prime}$ | $\top$ | $\top$ | $\top$ | $b^{\prime}$ | $\top$ | $b^{\prime}$
$\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$
$b^{\prime}$ | $b^{\prime}$ | $\top$ | $\top$ | $\top$ | $b^{\prime}$ | $\top$ | $b^{\prime}$
$a^{\prime}$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$

$\mathbin{/}_{W}$ | $a$ | $b$ | $\bot$ | $ba$ | $\top$ | $b^{\prime}$ | $a^{\prime}$
---|---|---|---|---|---|---|---
$a$ | $a$ | $\bot$ | $\bot$ | $b^{\prime}$ | $\top$ | $b^{\prime}$ | $\top$
$b$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$
$\bot$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$
$ba$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$
$\top$ | $a$ | $\bot$ | $\bot$ | $b^{\prime}$ | $\top$ | $b^{\prime}$ | $\top$
$b^{\prime}$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$ | $\top$
$a^{\prime}$ | $a$ | $\bot$ | $\bot$ | $b^{\prime}$ | $\top$ | $b^{\prime}$ | $\top$

Table 2. Tables for the residuals

Observe that

(7) $\displaystyle a$ $\displaystyle=a\circ_{W}a$ (8) $\displaystyle a$
$\displaystyle\not\leq a\circ_{W}(b\backslash_{W}b)\circ_{W}a$

since $a$ is a transitve and dense relation, and $b\backslash_{W}b=b^{\prime}$
and $a\circ_{W}b^{\prime}=\bot$ whence

$a\circ_{W}(b\backslash_{W}b)\circ_{W}a=a\circ_{W}b^{\prime}\circ_{W}a=\bot\circ_{W}a=\bot\not\geq
a$

in ${\mathcal{A}}$.

Next assume that ${\mathcal{A}}$ can be represented on a finite base. That is,
we have a finite set $U^{\prime}$ and a transitive relation $V\subseteq
U^{\prime}\times U^{\prime}$ such that $\mathcal{A}$ is isomorphic to
${\mathcal{B}}=(B,\leq,\circ_{V},\backslash_{V},/_{V})\subseteq(\wp(V),\subseteq,\circ_{V},\backslash_{V},/_{V})$.

Let $(x,y)\in a$ in $\mathcal{B}$. We show that $(x,y)\in
a\circ_{V}(b\backslash_{V}b)\circ_{V}a$. Since $a$ is dense ($a\leq
a\circ_{V}a$), we have $z\in U^{\prime}$ (not necessarily distinct from $x,y$)
such that $(x,z),(z,y)\in a$. Continuing in the same vein (using
$(x,z),(z,y)\in a\leq a\circ_{V}a$, etc.) we can create an $a$-path of
elements $(x=z_{0},z_{1},\dots,z_{i},\dots,z_{n-1},z_{n}=y)$ of unlimited
length (i.e. each $(z_{i},z_{i+1})\in a$). Since $U^{\prime}$ is finite, after
finitely many iterations of the process of expanding the path we get that
$z_{i}=z_{j}$ for two elements of the path. Since $a$ is transitive
($a\circ_{V}a\leq a$), we have $(x,z_{i}),(z_{i},z_{j}),(z_{j},y)\in a$. Thus
$(z_{i},z_{j})=(z_{i},z_{i})\in V$, whence $(z_{i},z_{i})\in
b\backslash_{V}b$. We get that $(x,y)\in
a\circ_{V}(b\backslash_{V}b)\circ_{V}a$, contradicting that $\mathcal{B}$
should satisfy (8). ∎

##  3\. Conclusions

The above proof shows that ${\sf RRS}$ does not have the _finite base
property_ (fbp) w.r.t. implications: the implication

(9) $(a\leq a\circ a\land a\circ a\leq a)\to a\leq a\circ b\backslash b\circ
a$

is not ${\sf RRS}$-valid but cannot be refuted in any ${\sf RRS}$ with finite
base.

The solution to the following problem would answer the question of the _finite
relational model property_ (frmp) for the LC: whether non-derivable sequents
of LC are refutable in relational models with finite bases.

######  Problem 3.1.

Does ${\sf RRS}$ have the finite base property w.r.t. atomic formulas $t\leq
s$ (where $t,s$ are $\\{\circ,\backslash,/\\}$-terms of variables).

We can ask the same problems for the square version ${\sf RRS}^{\Box}$ of
${\sf RRS}$. These problems are related to that version of the LC that allows
the empty word and sequents with empty antecedents [La61]. The corresponding
${\sf RRS}$ is a class with reflexive units $W$, thus every
$x\mathbin{\backslash}_{W}x$ and $x\mathbin{/}_{W}x$ subsume the identity
relation, see Andréka and Mikulás [AM94]. Thus the proof of Theorem 2.1 does
not work in this case.

######  Problem 3.2.

Does ${\sf RRS}^{\Box}$ have the finite representation and finite base
properties?

###  Acknowledgements: Thanks are due to Robin Hirsch and Ian Hodkinson for
useful discussions.

## References

  * [ABD92] A.R. Anderson, N.D. Belnap and J.M. Dunn, Entailment. The Logic of Relevance and Necessity. Vol. II., Princeton University Press, 1992.
  * [AM94] H. Andréka and Sz. Mikulás, “Lambek calculus and its relational semantics: completeness and incompleteness”, Journal of Logic, Language and Information, 3:1–37, 1994.
  * [vB91] J. van Benthem, Language in Action, North-Holland, 1991.
  * [BDM09] K. Bimbó, J.M. Dunn and R.D. Maddux, “Relevance logic and relation algebras”, Review of Symbolic Logic, 2(1):102–131, 2009.
  * [HH02] R. Hirsch and I. Hodkinson, Relation Algebras by Games, North-Holland, 2002.
  * [HM11] R. Hirsch and Sz. Mikulás, “Positive fragments of relevance logic and algebras of binary relations”, _Review of Symbolic Logic_ , 4(1):81–105, 2011.
  * [La58] J. Lambek, “The mathematics of sentence structure”, American Mathematical Monthly, 65:154–170, 1958.
  * [La61] J. Lambek, “On the calculus of syntactic types”, in R. Jakobson (editor), Structure of Language and Its Mathematical Aspects, pages 166–178, AMS, 1961.
  * [Ma10] R.D. Maddux, “Relevance logic and the calculus of relations”, Review of Symbolic Logic, 3(01):41–70, 2010.
  * [Mi09] Sz. Mikulás, “Algebras of relations and relevance logic”, Journal of Logic and Computation, 19:305–321, 2009.
  * [Mi15] Sz. Mikulás, “Lower Semilattice-Ordered Residuated Semigroups and Substructural Logics”, Studia Logica, 103:453-478, 2015.
  * [RM73] R. Routley and R.K. Meyer, “The semantics of entailment (I)”, in H. Leblanc (editor), Truth, Syntax and Modality, pages 199–243, North-Holland, 1973.


# Cosmological advection flows in the presence of primordial black holes as
dark matter and formation of first sources

A. Kashlinsky  Code 665, Observational Cosmology Lab, NASA Goddard Space
Flight Center, Greenbelt, MD 20771 and SSAI, Lanham, MD 20770; email:
<EMAIL_ADDRESS>
###### Abstract

In the inflation-based cosmology the dark matter (DM) density component starts
moving with respect to the universal expansion at $z_{\rm eq}\sim 3,200$ while
baryons remain frozen until $z_{\rm rec}\sim 1,100$. It has been suggested
that in this case post-linear corrections to the evolution of small
fluctuations would result, for the standard $\Lambda$-dominated cold DM (CDM)
model, in delayed formation of early objects as supersonic advection flows
develop after recombination, so baryons are not immediately captured by the DM
gravity on small scales. We develop the hydrodynamical description of such
two-component advection and show that, in the supersonic regime, the advection
within irrotational fluids is governed by the gradient of the difference of
the kinetic energies of the two (DM and baryonic here) components. We then
apply this formalism to the case where DM is made up of LIGO-type black holes
(BHs) and show that there the advection process on scales relevant for early
structure collapse will differ significantly from the earlier discussed (CDM)
case because of the additional granulation component to the density field
produced during inflation. The advection here will lead efficiently to the
common motion of the DM and baryon components on scales relevant for collapse
and formation of first luminous sources. This leads to early collapse making
easier to explain the existence of supermassive BHs observed in quasars at
high $z>7$. The resultant net advection rate reaches minimum around
${{}_{<}\atop{}^{\sim}}10^{9}M_{\odot}$ and subsequently rises to a secondary
maximum near the typical mass of $\sim 10^{12}M_{\odot}$, which may be an
important consideration for formation of galaxies at
$z{{}_{<}\atop{}^{\sim}}$(a few).

##  I Introduction

The linear Newtonian growth theory of cosmic structures is now established
(e.g. Harrison, 1967). However, Tseliakhovich and Hirata (2010) noted the
importance of advection flows in post-linear approximation for standard
cosmological model, where the dark matter (DM) density field approaches the
power $P\propto k^{-3}$ at large $k$, while preserving the initial Harrison-
Zeldovich shape $P\propto k$ on scales, $\sim k^{-1}$, exceeding the horizon
at matter-radiation equality, redshift $z_{\rm eq}\simeq 3200$, when the DM
component starts growing. This results in a coherent velocity field on scales
${{}_{<}\atop{}^{\sim}}$ a few Mpc corresponding to early collapsing
structures. At the same time baryons remain frozen into the comoving
(radiation) frame until recombination, $z_{\rm rec}\simeq 1100$, when DM is
already moving supersonically relative to baryon’s sound speed (Ricotti  _et
al._ , 2008). Baryons then find themselves moving supersonically in highly
coherent DM flows likely delaying their collapse to form luminous objects
Tseliakhovich and Hirata (2010). Much work followed on this potentially
important effect (e.g. Maio  _et al._ , 2011; Stacy  _et al._ , 2011; Greif
_et al._ , 2011; Ahn, 2016; Ahn and Smith, 2018; Blazek  _et al._ , 2016;
Schmidt, 2016; Hirata, 2018).

The advection hydrodynamics would be modified if LIGO-mass primordial black
holes (PBHs) make up DM. The possibility was proposed to explain their
apparent merger rate (Bird  _et al._ , 2016; Clesse and García-Bellido, 2017),
or the source-subtracted cosmic infrared background (CIB) (Kashlinsky, 2016),
where this conjecture naturally reproduces the amplitude and shape of the
earlier uncovered near-IR source-subtracted CIB fluctuations (Kashlinsky  _et
al._ , 2005) and their strong spatial coherence with cosmic X-ray background
implying populations containing substantial BH proportions (Cappelluti  _et
al._ , 2013). See review Kashlinsky  _et al._ (2018). Following the first LIGO
detection during its short engineering test run (Abbott  _et al._ , 2016) the
O1+O2 run uncovered 10 significant BH mergers of $\sim 10-50M_{\odot}$ masses
with low-to-zero spins (The LIGO Scientific Collaboration  _et al._ , 2019).
The O3/O4 LIGO runs at increasing sensitivity should provide critical insights
into the possible PBH-DM collusion. Theoretical mechanisms for such PBHs are
discussed in García-Bellido  _et al._ (1996); Jedamzik (1997); Musco  _et al._
(2005); Carr  _et al._ (2016); Sasaki  _et al._ (2018). This decade’s new EM-
based efforts should shed critical light on the PBH-DM linkage (Kashlinsky
_et al._ , 2019a), particularly in the CIB realm (Kashlinsky  _et al._ ,
2019b, c). Recent discussions of cosmogonical implications of this conjecture
include (Carr and Silk, 2018; Inman and Ali-Haïmoud, 2019; Hasinger, 2020).

If PBHs constitute DM their granulation produces an additional power component
(Meszaros, 1974, 1975) of shot-noise type on scales beyond the horizon scale
at the time of their formation (Kashlinsky, 2016), which modifies the
advection mechanism compared to Tseliakhovich and Hirata (2010). Here we show
that the advection in the presence of PBHs as DM would be efficient in early
equalizing the velocity fields of DM and baryons. This enables a sufficiently
early formation of the BH seeds to explain the existence of supermassive BHs
(SMBHs) observed in quasars at $z{{}_{>}\atop{}^{\sim}}7$. The advection rate
resulting from the $\Lambda$DM and PBH components is such that the equalizing
the two velocity fields becomes less efficient for total halo masses in the
$10^{9}-10^{11}M_{\odot}$ range and peaks up again for haloes around a few
times $10^{12}M_{\odot}$, the typical masses of modern galaxies potentially
affecting/determining formation of structures on galactic mass-scales.

We assume two primordial density field components: 1) from the standard
inflationary era ($\Lambda$DM), and 2) the PBH granulation fluctuation
component generated after inflation. We adopt DM fraction, $f_{\rm PBH}\leq
1$, for PBHs; the results are scaleable with $f_{\rm PBH}$. (Lacki and Beacom,
2010) argued that PBHs must either make all the DM or contribute (almost)
nothing since otherwise they accrete particle DM producing highly luminous
annihilation sources contradicting $\gamma$-ray observations.

##  II Cosmological advection flows

We use hydrodynamical descriptions for the common evolution of the PBH-DM and
baryon components. While this is obviously applicable to baryons, some caveats
are in order for the PBHs. If PBHs of mass $M_{\rm PBH}$ contribute the
fraction $f_{\rm PBH}$ to the the average density, their mean comoving
separation is $\bar{r}_{\rm PBH}\simeq 0.44f_{\rm PBH}^{1/3}(M_{\rm
PBH}/30M_{\odot})^{-1/3}$Kpc and the number of PBHs contained in comoving
scale $r$ is $N_{\rm PBH}\simeq(r/\bar{r}_{\rm PBH})^{3}\gg 1$ on scales of
relevance here. Provided we consider scales encompassing $N_{\rm PBH}\gg 1$
their evolution is described by the stellar-dynamical Jeans equations, which,
absent stellar-dynamical pressure (in virialized stellar systems), are
equivalent to Euler’s fluid equations; we call both the “Euler-Jeans”
equations as they are derived from moments of Boltzman’s equation. We consider
comoving scales $r{{}_{>}\atop{}^{\sim}}0.03$Mpc, so the Euler-Jeans equations
are valid provided $M_{\rm PBH}\ll 10^{7}f_{\rm PBH}(r/0.03{\rm
Mpc})^{3}M_{\odot}$; this range covers the LIGO-type PBH masses examined
throughout. The PBH mass-range considered in the analysis here is comfortably
within the $1\\!-\\!\sigma$ upper limits on the granulation power, discussed
below, from the Lyman-forest observations and simulations Afshordi  _et al._
(2003); Murgia  _et al._ (2019).

We follow standard notations for the Newtonian evolution of density
fluctuations and their flows for the two components, “d” (DM) and “b”
(baryons). After recombination, the Euler-Jeans equations for their evolution
are given by eqs 6 of Tseliakhovich and Hirata (2010):

$\dot{\boldsymbol{v}}_{\rm d,b}+H\boldsymbol{v}_{\rm
d,b}+a^{-1}(\boldsymbol{v}_{\rm
d,b}\cdot\boldsymbol{\nabla})\boldsymbol{v}_{\rm
d,b}=-a^{-1}\boldsymbol{\nabla}\phi-\frac{a^{-1}}{{\rho}_{\rm
d,b}}\boldsymbol{\nabla}p_{\rm d,b}$ (1)

where $a=(1+z)^{-1},H(z)=\dot{a}/a$, $\phi$ is the gravitational potential and
$p$ is the presure in each component. We use the Lamb transformation
$(\boldsymbol{v}\cdot\boldsymbol{\nabla})\boldsymbol{v}=\frac{1}{2}\boldsymbol{\nabla}(v^{2})-\boldsymbol{v}\times\boldsymbol{\omega}$
where vorticity $\boldsymbol{\omega}=\boldsymbol{\nabla}\times\boldsymbol{v}$
and $v=|\boldsymbol{v}|$ (Lamb, 1975). Gravity is a potential force inducing
irrotational flows so vorticity is small in linear regime (Peebles, 1969;
Efstathiou and Jones, 1979); thus
$(\boldsymbol{v}\cdot\boldsymbol{\nabla})\boldsymbol{v}=\frac{1}{2}\boldsymbol{\nabla}(v^{2})$.

The relative baryons–DM velocity is important for collapsed structure
formation. Once the two components move together, the gas can collapse in the
formed DM halos and fragment subject to its cooling and fragmentation
efficiency (Hoyle, 1953; Rees and Ostriker, 1977; Kashlinsky, 1982). Since
both flows are driven by the same gravitational potential $\phi$, subtracting
“d” from “b” components in (1) gives, for irrotational flows:

$\frac{\partial}{\partial t}(a{\boldsymbol{V}}_{\rm
bd})=-\boldsymbol{\nabla}({\cal K}_{\rm b}-{\cal K}_{\rm
d}+c_{s}^{2}\delta_{b})$ (2)

where $\boldsymbol{V}_{\rm bd}\equiv\boldsymbol{v}_{\rm b}-\boldsymbol{v}_{\rm
d}$, each component’s kinetic energy per unit mass ${\cal
K}=\frac{1}{2}v^{2}$, $c_{s}^{2}\equiv\frac{\partial p_{\rm
b}}{\partial\rho_{\rm b}}$ is the adiabatic sound speed squared, assumed
uniform, and DM is taken to be pressurerless. Eq. 2 with the continuity
equation for DM later are equivalent to equations (6) of Tseliakhovich and
Hirata (2010) when baryons comove with DM. DM dominates the peculiar gravity
and starts moving at matter-radiation equality, $z_{\rm eq}\simeq 3,200$,
whereas baryons start growing with $\boldsymbol{v}_{b}\simeq 0$ at $z_{\rm
rec}\simeq 1,090$. Both components move supersonically after recombination
Ricotti  _et al._ (2008); Tseliakhovich and Hirata (2010). Relevant
solutions/consequences of (2) are:

  1. 1.

Eq. 2 provides an exact description of the relative motions of the
irrotational flow components until the DM shell-crossing. For supersonic flows
under the same potential force (gravity) the advection is influenced by the
gradient of the difference in kinetic energies, ${\cal K}$, of the flow
components, driving their (kinematic) mixing.

  2. 2.

In steady-state $\boldsymbol{v}_{\rm b}=\boldsymbol{v}_{\rm d}$, i.e. baryons
and DM move coherently on all scales (to within $c_{s}$). This solution always
exists even at time/space-varying $\boldsymbol{v}$ in the linear
approximation, when the RHS of (2) vanishes. This regime is reached as the
result of advection.

  3. 3.

The general solution to eq.2 is the sum of two parts, the first of which is
$\boldsymbol{V}_{\rm bd}(z)=-\boldsymbol{v}_{\rm
d}(z\\!\\!=\\!\\!1000)(1+z)/1000$ if $v_{\rm b}=0$ initially. This solution
fully describes the component of $\boldsymbol{V}_{\rm bd}$ perpendicular to
$\boldsymbol{\nabla}({\cal K}_{\rm d})$. When it dominates advection flows may
suppress/delay the onset of baryonic collapse Tseliakhovich and Hirata (2010).

  4. 4.

Because $v_{\rm d}$ grows with time, the advection term on the RHS of eq. 2
may become important and also more efficient in decreasing
$\boldsymbol{V}_{\rm bd}$. The remaining relative baryon-DM velocity,
$\boldsymbol{V}_{\rm bd}$, is then dominated by the solution along the
direction defined by $\boldsymbol{\nabla}({\cal K}_{\rm d})$. If $r$ denotes
the coordinate along this direction, the baryons eventually catch up with the
motion of DM, i.e. $\dot{V}_{\rm bd}>0$, if $\partial(v_{\rm d}^{2})/\partial
r<0$ ($v_{\rm d}^{2}$ decreases with increasing distance) and vice versa. The
rate at which baryons catch up with the moving DM is given by the gradient of
the kinetic energy of the DM bulk flow, $\partial{\cal K}_{\rm d}/\partial r$.

  5. 5.

Evolution of $\boldsymbol{V}_{\rm bd}$, (2), along the line defined by
$\boldsymbol{\nabla}({\cal K}_{\rm d})$ starts at $\boldsymbol{V}_{\rm
bd}=-\boldsymbol{v}_{\rm d}$ proceeding with baryons catching up with DM
($\boldsymbol{v}_{\rm b}=\boldsymbol{v}_{\rm d}$) on advection timescale
$t_{\cal A}\sim L_{\cal A}/v_{\rm d}$ where $L_{\cal
A}=a{\color[rgb]{0,0,0}v_{\rm d}^{2}}/|\boldsymbol{\nabla}({\cal K}_{\rm
b}-{\cal K}_{\rm d})|$; the rate at which velocities are equalized, $v_{\rm
d}/t_{\cal A}$, being closely related to the advection rate defined below.
Consequently, for the $\Lambda$DM power spectrum the advection terms can
become comparable to or larger than the expansion time (Sec. II.C in
Tseliakhovich and Hirata, 2010). If the power spectrum (e.g. PBH-DM) is such
that $t_{\cal A}$ wins over the expansion time the balancing of the b,d flows
proceeds more efficiently, similar to analogous problems in atmospheric
advection Vasil’kov (1975).

  6. 6.

When $v_{\rm b}=0$ to within $c_{s}$ to begin with, but $v_{\rm d}\gg c_{s}$
so that $v_{\rm b}^{2}\ll v_{\rm d}^{2}$, we can Fourier transform velocities,
rewriting (2)

$\frac{\partial[a\boldsymbol{U}_{\rm bd}(\boldsymbol{k})]}{\partial
t}\propto\boldsymbol{k}\int U_{\rm d}(\boldsymbol{k}-\boldsymbol{k}_{1})U_{\rm
d}(\boldsymbol{k}_{1})d^{3}\boldsymbol{k}_{1}$ (3)

where $\boldsymbol{U}_{\rm bd}$ is the Fourier mode of $\boldsymbol{V}_{\rm
bd}$, and the RHS integral represents the convolution of
$U_{d}\equiv|\boldsymbol{U}_{{\rm d},\boldsymbol{k}}|$ with itself, $U_{{\rm
d},\boldsymbol{k}}\star U_{{\rm d},\boldsymbol{k}}$. Consequently, Fourier
harmonics no longer evolve independently (e.g. Efstathiou, 1990).

  7. 7.

In the presence of rotation/vorticity the RHS of eq. 2 contains the additional
term $(\boldsymbol{v}_{\rm b}\times\boldsymbol{\omega}_{\rm
b}-\boldsymbol{v}_{\rm d}\times\boldsymbol{\omega}_{\rm d})$.

  8. 8.

When $c_{s}\ll v_{d}$ and $v_{\rm b}=0$ at recombination the relative baryon-
DM velocity evolves as

$\boldsymbol{V}_{\rm bd}=-\frac{1+z}{1+z_{\rm rec}}\boldsymbol{v}_{\rm
d,rec}-a^{-1}\int_{t(z_{\rm rec})}^{t(z)}\boldsymbol{{\cal A}}_{\cal K}\;adt$
(4)

where $\boldsymbol{{\cal A}}_{\cal K}\equiv-a^{-1}\boldsymbol{\nabla}{\cal
K}_{\rm d}$ is the local advection rate. When the latter dominates the first
term on the RHS of eq. 4 the advection speeds up the equalizing of the baryon-
DM velocity. Eq. 4 is correct to $O(v_{\rm b}^{2}/v_{\rm d}^{2})$ leading to
steady-state at $V_{\rm bd}=0$.

Below we adopt: $h\equiv H_{0}/$(100 km/sec/Mpc)$=0.7,\Omega_{\rm
d}h^{2}=0.11,\Omega_{\rm b}h^{2}=0.023,\Omega_{{\rm m},0}=\Omega_{\rm
d}+\Omega_{\rm b}=0.3,\Omega_{\Lambda}+\Omega_{\rm d}+\Omega_{\rm
b}=1,\sigma_{8}=0.9$. At $z_{\rm rec}\gg z\gg 1$, $\Omega_{\rm
m}(z)=\Omega_{{\rm m},0}(1+z)^{3}/[\Omega_{{\rm
m},0}(1+z)^{3}+\Omega_{\Lambda}]\simeq 1$, the Hubble constant $H(z)\simeq
H_{0}\sqrt{\Omega_{\rm m,0}}(1+z)^{3/2}$, $t_{\rm
cosm}(z)\simeq\frac{2}{3}H^{-1}(z)\simeq 0.2(\frac{1+z}{20})^{-3/2}$ Gyr.

##  III Advection flows for PBH dark matter

Following Kashlinsky (2016) we assume that the power spectrum, $P_{\rm m}(k)$,
of matter fluctuations responsible for structure formation at recombination,
$z\simeq 10^{3}$, is made up of 1) the $\Lambda$DM component from the
inflationary period and 2) the component from LIGO-type PBHs contributing a
fraction $f_{\rm PBH}$ of the DM:

$P_{\rm m}(k)=P_{\Lambda{\rm DM}}(k)+P_{\rm PBH}=P_{\Lambda{\rm
DM}}(k)+1.2\times 10^{-8}f_{\rm PBH}\left(\frac{M_{\rm
PBH}}{30M_{\odot}}\right)\left(\frac{1+z}{1000}\right)^{-2}{\rm Mpc}^{3}$ (5)

The PBH component $P_{\rm PBH}\propto f_{\rm PBH}M_{\rm PBH}$. For an extended
PBH mass function (e.g. Clesse and García-Bellido, 2018; Carr and Silk, 2018;
Inman and Ali-Haïmoud, 2019) $M_{\rm PBH}$ represents the effective PBH mass,
it being the eigen-value after suitably averaging over the PBH mass-function.

The advection eq. 2, when starting at $z_{\rm rec}$ at $\boldsymbol{v}_{\rm
b}=0$, must be complemented with the continuity equations for each component,
$\left[\dot{\delta}_{\rm
d,b}+a^{-1}\boldsymbol{\nabla}\cdot\boldsymbol{v}_{\rm
d,b}\right]=-a^{-1}\boldsymbol{\nabla}\cdot(\delta_{\rm
d,b}\boldsymbol{v}_{\rm d,b})$. Its general solution is the sum of two terms:
$\boldsymbol{v}_{d}\equiv\boldsymbol{v}_{{\rm d,1}}+\Delta\boldsymbol{v}_{\rm
d}$. The first term satisfies, in this gauge, the sum in square brackets being
$\left[\ldots\right]=0$. Hence $\Delta\boldsymbol{v}_{\rm
d}=-\frac{\delta_{\rm d}}{1+\delta_{\rm d}}\boldsymbol{v}_{{\rm
d,1}}\simeq-\delta_{\rm d}\boldsymbol{v}_{{\rm d,1}}+O(\delta_{d}^{2}v_{{\rm
d,1}})$. Fourier-transforming $\delta_{\rm
d}(\boldsymbol{r}),\boldsymbol{v}_{\rm d}(\boldsymbol{r})$ into $\Delta_{\rm
d}(\boldsymbol{k}),\boldsymbol{U}_{\rm d}(\boldsymbol{k})$ for the
irrotational flow, gives $\boldsymbol{U}_{\rm d}=\boldsymbol{U}_{\rm
d,1}+(\boldsymbol{U}_{\rm d,1}\star\Delta_{\rm d})_{\boldsymbol{k}}$, where
${\boldsymbol{U}}_{{\rm d,1}}=-i\boldsymbol{k}a\dot{\Delta}_{\rm d}/k^{2}$ and
$\boldsymbol{U}_{{\rm d,1}}\star\Delta_{d}\equiv\int\boldsymbol{U}_{{\rm
d,1}}(\boldsymbol{k}_{1})\Delta_{\rm
d}(\boldsymbol{k}-\boldsymbol{k}_{1})d^{3}\boldsymbol{k}_{1}$. The small-to-
large scale mode-coupling is less pronounced here because of the absence of
the $k^{-1}$ weight for $\Delta_{\rm d}$ compared to (3). For the white-noise
PBH component dominating small-scale power, $\Delta(k)$=const, the convolution
integral $\propto\int\boldsymbol{U}_{{\rm
d,1}}(\boldsymbol{k}_{1})d^{3}\boldsymbol{k}_{1}=\boldsymbol{v}_{\rm d,1}(0)$
is constant with $k$, while $U_{\rm d,1}$ increases toward small scales.

Hence, like Tseliakhovich and Hirata (2010), we adopt the linear growth
velocity evolution, i.e.
$k^{2}P_{v}(k,z)=\left[\frac{a\dot{\Delta}(k,z)}{\Delta(k,z)}\right]^{2}P_{m}(k,z)$.
In linear approximation $\Delta(k)\propto a$, but to describe the full non-
linear evolution of density fields requires solving (1) Tseliakhovich and
Hirata (2010) or using the ”stable-clustering” approach (Peacock and Dodds,
1996). Such corrections are scale-dependent but small (Tseliakhovich and
Hirata (2010) show in Fig. 2 the ${{}_{<}\atop{}^{\sim}}$O(10%) non-linear
corrections to $P_{\Lambda{\rm DM}}$ at $z{{}_{>}\atop{}^{\sim}}40$ peaking
around $k\simeq 100-300$Mpc-1) so we take $P_{v}=H^{2}(z)[\Omega_{\rm
m}(z)]^{1.2}P_{\rm m}/k^{2}$ with the caveats discussed later. Fig. 1,left
plots the $\Lambda$DM and PBH parts of $kP_{\rm m}(k)$ at
$z=1000,900,700,100,20$ for standard cosmological parameters using
CAMB111https://lambda.gsfc.nasa.gov/toolbox/tb_camb_form.cfm. Fig. 1,right
shows the rms density fluctuation, $\delta_{\rm
m,rms}^{2}=\frac{1}{2\pi^{2}}\int P_{\rm m}W^{2}(kr)k^{2}dk$ from these
components, where $W(y)=3j_{1}(y)/y$.

Figure 1: Left: Logarithmic contribution to the DM velocity variance, $kP_{\rm
m}(k)$, times the power growth factor for Einstein-deSitter regime. Lines, for
the $\Lambda$DM component, are for $z=1,000$ (red), 900 (green), 700 (indigo),
100 (blue) and 20 (black). Thick dashes show DM-PBH contribution for $f_{\rm
PBH}M_{\rm PBH}=30M_{\odot}$ Right: RMS density fluctuation dispersion times
the Einstein-deSitter growth-factor vs the comoving radius subtending the
matter mass in the upper horizontal axis. Solid, dashed, dotted lines show the
$\Lambda$DM, DM-PBH term, and their sum.

The mean kinetic energy of the DM component is:

$\bar{\cal K}_{\rm
d}=\frac{1}{2}\sigma_{v}^{2}=\frac{1}{4\pi^{2}}\Omega_{m}^{1.2}(z)a^{2}H^{2}\int_{0}^{\infty}P_{\rm
m}(k)W^{2}(kr)dk\equiv\frac{1}{2}(\sigma_{v,\Lambda{\rm
DM}}^{2}+\sigma_{v,{\rm PBH}}^{2})$ (6)

where the two RHS terms arise from the two power components, (5). The 1-D
velocity variance, $\sigma_{v}^{2}$, is related to the ”dot” velocity
correlation function,
$\xi_{v}(\boldsymbol{r})=\langle\boldsymbol{v}(\boldsymbol{r}^{\prime})\cdot\boldsymbol{v}(\boldsymbol{r}^{\prime}+\boldsymbol{r})\rangle$(Peebles,
1980; Vittorio  _et al._ , 1986; Kashlinsky, 1992); the relative 1-D velocity
relevant for turn-around/collapse is $v_{\rm
rel}(r)=\sqrt{\sigma_{v}^{2}(0)-\sigma_{v}^{2}(r)}$. Fig. 2,left shows the
$\Lambda$DM component of the velocity variance, $\sigma_{v,\Lambda{\rm DM}}$.
The resultant flow of the $\Lambda$DM component is highly coherent out to
comoving $r\sim$ (a few)Mpc. We define the velocity dispersion slope,
$n_{v}\equiv\frac{\partial\ln\sigma_{v}(r)}{\partial\ln
r}=\frac{1}{2}\frac{\partial\ln\bar{\cal K}_{\rm d}(r)}{\partial\ln r}$. On
scales where $\sigma_{v}$ is highly coherent the relative 1-D velocity is
$v_{\rm rel}\simeq\sqrt{-2n_{v}}\sigma_{v}(0)$. Fig. 2,right plots $n_{v}$ vs
$r$ for the $\Lambda$DM component. For the PBH component $n_{v}=-0.5$. Eq. 2
of Tseliakhovich and Hirata (2010) gives $c_{s}\simeq 5.5-1.4$ km sec-1 at
$z=1000-100$ with $c_{s}<0.5$ km/sec at $z<30$ (until reionization/reheating).
The acoustic pressure term is subdominant compared to the DM kinetic energy in
(2) until first sources form and reheat the baryonic gas.

Figure 2: Velocity parameters for the $\Lambda$DM component. Left: Solid lines
show the velocity dispersion vs comoving scale at different $z$ in the color
notation of Fig. 1,left. Thick dashed line shows the 1-D relative velocity
$v_{\rm rel}(r)\equiv\sqrt{\sigma_{v}^{2}(0)-\sigma_{v}^{2}(r)}$ at $z=20$.
Right: Logarithmic slope of the $\Lambda$DM component’s velocity dispersion,
$n_{v}$ vs $r$, for the color lines on the left. For PBH component
$n_{v}=-0.5$.

The mean of $\boldsymbol{V}_{\rm bd}\cdot$(eq. 2) is zero over a finite volume
because of the conservation of relative energy. At each point the two
components get mixed by the instantaneous gradient of the (difference in)
kinetic energy with non-zero rms rate. The rms measure of the advection rate
due to the PBH component with the mean bulk kinetic energy $\bar{\cal K}$ is
(with $\int_{0}^{\infty}W^{2}(y)dy\simeq 2$):

$\bar{{\cal A}}_{\cal K}\equiv-a^{-1}\frac{\partial\bar{\cal K}_{\rm
d}}{\partial r}=-a^{-1}n_{v,\Lambda{\rm DM}}\frac{\sigma_{v,\Lambda{\rm
DM}}^{2}}{r}+\frac{\Omega_{m}^{1.2}(z)}{\pi^{2}}aH^{2}P_{\rm
PBH}(z)\frac{1}{r^{2}}$ (7)

We express the advection rate in km/sec/Gyr, so ${\cal A_{K}}=1$ equalizes
relative motions of 1 km/sec in 1 Gyr.

The advection rate $\bar{{\cal A}}_{\cal K}=-n_{v}\sigma_{v}^{2}/r$ is
independent of $z$ in the Einstein-deSitter regime and is added to the
reduction due to expansion given by the first term on the RHS of eq. 4. The
expansion term equalizes the baryon-DM velocities at the (reduction) rate of
${\cal R}_{V}\equiv\frac{(1+z)}{(1+z_{\rm rec})}v_{\rm
d,rec}/t_{{\color[rgb]{0,0,0}\rm cosm}}(z)$. The comparison between the two is
shown in Fig. 3 assuming, as an example, the initial $v_{\rm d}$ at
recombination to be given by $\sigma_{v}\simeq 18$km/sec at $z=1000$. One can
see that the advection rate from PBH-DM dominates for scales corresponding to
those relevant for the formation of first collapsed structures already at
$z{{}_{<}\atop{}^{\sim}}50$.  The precise numerology may be affected by the
following: 1) $n_{v}$ changes in the non-linear regime at scales below
$r(\delta_{\rm m,rms}\\!\\!=\\!\\!1)\simeq 10^{-3}(f_{\rm PBH}M_{\rm
PBH}/30M_{\odot})^{1/3}[(1+z)/1000]^{-2/3}h^{-1}$Mpc, and 2) scale-dependent
growth due to post-linear corrections; Figs.2,3 of Peacock and Dodds (1996)
show that these corrections are small for white-noise power at $\delta_{\rm
m,rms}<1$. The first of these would decrease the effective advection somewhat
by reducing the effective $\sigma_{v}$, while the second would increase it by
lowering $n_{v}$.

Figure 3: The advection rate (solid) for the $\Lambda$DM (red) and PBH (black)
components with $f_{\rm PBH}M_{\rm PBH}$ marked. Dotted lines show total
${\cal A_{K}}$. The advection from the PBH-DM component is more efficient in
equalizng DM and baryonic velocity components than the $\Lambda$DM alone at
scales where $\Lambda$DM component starts dominating the density fluctuation
per Fig. 1,right. The PBH advection rate continues to be higher than of the
$\Lambda$DM even when the PBHs contribute $f_{\rm PBH}M_{\rm PBH}\ll
30M_{\odot}$. Blue dash-triple-dotted lines show ${\cal R}_{V}$ for the marked
values of $v_{\rm d,rec},z$.

Fig. 3 shows ${\cal A_{K}}$ from the $\Lambda$DM, PBH DM power components
assuming $f_{\rm PBH}M_{\rm PBH}=3,30M_{\odot}$. As a consequence of the
larger coherence of $\Lambda$DM velocity field compared to its density field,
the advection rate, with ${\cal A_{K}}_{\rm,PBH}\propto f_{\rm PBH}M_{\rm
PBH}$, is controlled predominantly by the PBH power component even at scales
where the overall density field is already dominated by the $\Lambda$DM power.
The red line appears in agreement with the formalism in Tseliakhovich and
Hirata (2010), which, however, is not applicable in the presence of the PBH-DM
power component. The advection on small scales relevant to first source
formation is driven by the PBH component even when $f_{\rm PBH}M_{\rm PBH}\ll
30M_{\odot}$.

##  IV Formation of SMBHs and galaxies in the presence of advection

Fig. 3 shows that if PBHs make up the DM, the advection makes baryons comoving
with DM quickly after recombination with ${\cal
A_{K}}{{}_{>}\atop{}^{\sim}}10-100$km/sec/Gyr. Baryons and DM participate in
formation/evolution of the same DM halos at the epochs they separate from the
comoving frame and collapse without delay, i.e. the normal evolution for
growth/collapse of density fluctuations applies Press and Schechter (1974). At
the same time the same granulation component ensures an early collapse of
first haloes and potentially early formation of compact objects.

Of specific relevance is the existence of SMBHs implied by QSO observations
deep inside the reionization epoch. The following are noteworthy here in
increasing $z$: 1) an ultraluminous quasar with $M_{\rm SMBH}\simeq 1.3\times
10^{10}M_{\odot}$ at $z=6.3$ (Wu  _et al._ , 2015), 2) a QSO at $z=7.1$
implying $M_{\rm SMBH}\simeq 2\times 10^{9}M_{\odot}$ (Mortlock  _et al._ ,
2011), and 3) a QSO at $z=7.5$ with $M_{\rm SMBH}\simeq 8\times
10^{8}M_{\odot}$ (Bañados  _et al._ , 2018). Their implications are
significant for standard CDM model, although models have been proposed how to
form them at $z\sim 12-15$ while reproducing CIB constraints (Yue  _et al._
(2013), cf. Ricarte  _et al._ (2019)). The difficulty stems from the limited
power in the $\Lambda$DM component (Fig. 1,right) Efstathiou and Rees (1988);
Kashlinsky (1993) and the advection efficiency from this component (e.g.
Tanaka  _et al._ , 2013). Fig. 2 of Bañados  _et al._ (2018) shows that, for
the Eddington accretion rates, the required SMBH seed masses must be already
of order a few $10^{3}M_{\odot}$ at $z{{}_{>}\atop{}^{\sim}}50$.

Formation of compact objects in primordial composition haloes happens if
baryons there can cool and maintain certain temperatures (e.g. Bromm, 2013,
and refs therein). This marks the critical halo masses for gas to collapse at
$z$: in haloes where no H2 formed the temperature would be $T\sim T_{4}\equiv
10^{4}$K, if H2 formed it can reach $T\sim T_{3}\equiv 10^{3}$K. For the gas
to collapse, pressure gradients must be less than gravity, defining haloes
with masses $M(z){{}_{>}\atop{}^{\sim}}\left[\frac{4\pi(1+\delta_{\rm
col})}{3}\right]^{-1/2}\left(k_{\rm B}T/m_{\rm
p}G\right)^{3/2}[\bar{\rho}_{\rm m}(z)]^{-1/2}$, where $\delta_{\rm
col}=1.68$. (This assumed $T$=const and, given its approximate values in the
presence of the coolants, omitted factors of molecular weight and the slope of
the pressure gradient). This delineates two critical total halo masses,
$M_{4,3}$ corresponding to $T_{4,3}$ at $z$; $M_{3}(z)\simeq
10^{6}\left(\frac{1+z}{20}\right)^{-3/2}M_{\odot}$,
$M_{4}(z)=(T_{4}/T_{3})^{3/2}M_{3}(z)$. The halo gas mass will be a factor of
$\Omega_{\rm b}/\Omega_{{\rm m},0}\sim 0.15$ smaller. H2 can be destroyed by
the Lyman-Werner (LW) radiation from the LIGO-type PBHs as the gas accretion
onto them would lead to multi-temperature accretion disks emitting at $T_{\rm
acc.disk}\sim(M_{\rm PBH}/M_{\odot})^{-1/4}$keV, which may be important in any
modeling involving detailed structure and compact objects formation in the
case of the PBHs, with $M_{4}$ providing a more reliable estimate.

Figure 4: Number of standard deviations, $Q_{\rm col}$, for halos to collapse
at $z$ and cool their baryons to the given temperature. Black(green) lines are
for $f_{\rm PBH}M_{\rm PBH}=30(3)M_{\odot}$ and $Q_{\rm col}\propto[f_{\rm
PBH}M_{\rm PBH}]^{-1/2}$ over most of the range. Granularity would also
introduce non-Gaussianity Inman and Ali-Haïmoud (2019) lowering $Q_{\rm col}$.

We evaluate the halo collapse likelihoods, assuming it leads to compact
object’s formation from the baryonic gas provided it can cool to the required
$T$ and that the advection due to PBH is efficient at equalizing the DM and
baryon velocities. For a Gaussian density field the probability for halo of
total mass $M$ to collapse at $z$ is ${\cal
P}_{M}(z)=\frac{1}{2}$erfc$(\frac{Q_{\rm col}}{\sqrt{2}})$ with $Q_{\rm
col}\equiv\frac{\delta_{\rm col}}{\delta_{\rm m,rms}(M,z)}$ (Press and
Schechter, 1974). At $Q_{\rm col}^{2}\gg 1$ this reduces to ${\cal
P}_{M}(z)\simeq\frac{1}{\sqrt{2\pi}}Q_{\rm col}^{-1}\exp(-\frac{1}{2}Q_{\rm
col}^{2})$ (Kashlinsky, 1993); this last expression is already within $<15\%$
of the true ${\cal P}_{M}$ for $Q_{\rm col}\geq 2.25$. Fig. 4 shows $Q_{\rm
col}(z)$ for the two cooling regimes. In the PBH-DM paradigm one can have a
reasonable abundance of haloes, say of $Q_{\rm
col}(z){{}_{<}\atop{}^{\sim}}3-6$, with the gas collapsing on mass-scales up
to ${{}_{<}\atop{}^{\sim}}10^{5}M_{\odot}$ at $z{{}_{<}\atop{}^{\sim}}30-50$,
if only a subset of these systems subsequently forms seed BHs of $M_{\rm
BH,seed}\sim 10^{3-4}M_{\odot}$ via e.g. any of the numerous mechanisms
suggested (e.g. Latif and Ferrara, 2016, and refs therein). A number of
models, applicable here, have been developed for formation of SMBHs inside
collapsed halos at high $z$. The proposed mechanisms are a result of stellar
dynamical relaxation processes, typically involving evolution of dense stellar
systems, coupled with gas collapse and dynamic instabilities discussed by e.g.
Begelman and Rees (1978); Kashlinsky and Rees (1983); Begelman  _et al._
(2006). Since $Q_{\rm col}\propto[f_{\rm PBH}M_{\rm PBH}]^{-1/2}$, to form
collapsed halos with $T_{3}$ by $z\sim 40$ would require $f_{\rm PBH}M_{\rm
PBH}{{}_{>}\atop{}^{\sim}}3M_{\odot}$ setting $Q_{\rm
col}{{}_{<}\atop{}^{\sim}}6$. Because of the advection efficiency and the
extra power of the PBH-DM component the SMBHs appear to support the PBH-DM
connection conjecture.

In the PBH-DM paradigm the advection rate as shown in Fig. 3 for $f_{\rm
PBH}M_{\rm PBH}=30{\color[rgb]{0,0,0},3}M_{\odot}$ reaches minimum around dark
matter mass scales $M_{\rm d}\sim 10^{9}{\color[rgb]{0,0,0},5\times
10^{7}}M_{\odot}$ and then rises again. This may require modifications and
additional considerations in various discussions involving subsequent
structure formation (e.g. Carr and Silk, 2018). Naively this implies that
there may be a pause after the first collapse era and resurgence of collapse
and luminous source formation around masses of order modern galaxies. After
the pause the gas would collapse to form the presently observed galaxy
morphology depending on the halo spin (Kashlinsky, 1982). This can explain the
existence of early-type galaxies with already established morphology and
${{}_{>}\atop{}^{\sim}}3.5$Gyr-old stellar populations at $z\sim 1.5-2$ Dunlop
_et al._ (1996); Saracco  _et al._ (2009); Renzini (2006).

###### Acknowledgements.

Fernando Atrio-Barandela and Alexandre Vassilkov are thanked for discussions
on cosmological aspects and advection hydrodynamics respectively. Support from
NASA/12-EUCLID11-0003 “LIBRAE: Looking at Infrared Background Radiation
Anisotropies with Euclid” project is acknowledged.

## References

  * Harrison (1967) E. R. Harrison, Reviews of Modern Physics 39, 862 (1967).
  * Tseliakhovich and Hirata (2010) D. Tseliakhovich and C. Hirata, Phys. Rev. D 82, 083520 (2010), arXiv:1005.2416 [astro-ph.CO] .
  * Ricotti  _et al._ (2008) M. Ricotti, J. P. Ostriker, and K. J. Mack, Astrophys. J. 680, 829 (2008), arXiv:0709.0524 [astro-ph] .
  * Maio  _et al._ (2011) U. Maio, L. V. E. Koopmans, and B. Ciardi, MNRAS 412, L40 (2011), arXiv:1011.4006 [astro-ph.CO] .
  * Stacy  _et al._ (2011) A. Stacy, V. Bromm, and A. Loeb, ApJL 730, L1 (2011), arXiv:1011.4512 .
  * Greif  _et al._ (2011) T. H. Greif, S. D. M. White, R. S. Klessen, and V. Springel, Astrophys. J. 736, 147 (2011), arXiv:1101.5493 .
  * Ahn (2016) K. Ahn, Astrophys. J. 830, 68 (2016), arXiv:1603.09356 [astro-ph.CO] .
  * Ahn and Smith (2018) K. Ahn and B. D. Smith, Astrophys. J. 869, 76 (2018), arXiv:1807.04063 [astro-ph.CO] .
  * Blazek  _et al._ (2016) J. A. Blazek, J. E. McEwen, and C. M. Hirata, Phys. Rev. Lett. 116, 121303 (2016), arXiv:1510.03554 [astro-ph.CO] .
  * Schmidt (2016) F. Schmidt, Phys. Rev. D 94, 063508 (2016), arXiv:1602.09059 [astro-ph.CO] .
  * Hirata (2018) C. M. Hirata, MNRAS 474, 2173 (2018), arXiv:1707.03358 [astro-ph.CO] .
  * Bird  _et al._ (2016) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Clesse and García-Bellido (2017) S. Clesse and J. García-Bellido, Physics of the Dark Universe 15, 142 (2017), arXiv:1603.05234 [astro-ph.CO] .
  * Kashlinsky (2016) A. Kashlinsky, ApJL 823, L25 (2016), arXiv:1605.04023 .
  * Kashlinsky  _et al._ (2005) A. Kashlinsky, R. G. Arendt, J. Mather, and S. H. Moseley, Nature (London) 438, 45 (2005), astro-ph/0511105 .
  * Cappelluti  _et al._ (2013) N. Cappelluti, A. Kashlinsky, R. G. Arendt, A. Comastri, G. G. Fazio, A. Finoguenov, G. Hasinger, J. C. Mather, T. Miyaji, and S. H. Moseley, Astrophys. J. 769, 68 (2013), arXiv:1210.5302 [astro-ph.CO] .
  * Kashlinsky  _et al._ (2018) A. Kashlinsky, R. G. Arendt, F. Atrio-Barandela, N. Cappelluti, A. Ferrara, and G. Hasinger, Reviews of Modern Physics 90, 025006 (2018), arXiv:1802.07774 .
  * Abbott  _et al._ (2016) B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, and et al., Phys. Rev. D 93, 122003 (2016), arXiv:1602.03839 [gr-qc] .
  * The LIGO Scientific Collaboration  _et al._ (2019) The LIGO Scientific Collaboration, the Virgo Collaboration, B. P. Abbott, _et al._ , arXiv e-prints , arXiv:1901.03310 (2019), arXiv:1901.03310 [astro-ph.HE] .
  * García-Bellido  _et al._ (1996) J. García-Bellido, A. Linde, and D. Wands, Phys. Rev. D 54, 6040 (1996), arXiv:astro-ph/9605094 [astro-ph] .
  * Jedamzik (1997) K. Jedamzik, Phys. Rev. D 55, R5871 (1997), astro-ph/9605152 .
  * Musco  _et al._ (2005) I. Musco, J. C. Miller, and L. Rezzolla, Classical and Quantum Gravity 22, 1405 (2005), arXiv:gr-qc/0412063 [gr-qc] .
  * Carr  _et al._ (2016) B. Carr, F. Kühnel, and M. Sandstad, Phys. Rev. D 94, 083504 (2016), arXiv:1607.06077 [astro-ph.CO] .
  * Sasaki  _et al._ (2018) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Classical and Quantum Gravity 35, 063001 (2018), arXiv:1801.05235 [astro-ph.CO] .
  * Kashlinsky  _et al._ (2019a) A. Kashlinsky, Y. Ali-Haïmoud, S. Clesse, J. Garcia-Bellido, L. Amendola, L. Wyrzykowski, J. Annis, A. Arbey, R. G. Arendt, F. Atrio-Barandela, and et al., in  _Bull. Am. Astron. Soc_ , Vol. 51 (2019) p. 51.
  * Kashlinsky  _et al._ (2019b) A. Kashlinsky, R. G. Arendt, N. Cappelluti, A. Finoguenov, G. Hasinger, K. Helgason, and A. Merloni, Astrophys. J. 871, L6 (2019b), arXiv:1812.01535 [astro-ph.CO] .
  * Kashlinsky  _et al._ (2019c) A. Kashlinsky, R. G. Arendt, M. Ashby, F. Atrio-Barandela, V. Bromm, N. Cappelluti, S. Clesse, A. Comastri, G. Fazio, and S. Driver, in  _Bull. Am. Astron. Soc_ , Vol. 51 (2019) p. 37, arXiv:1903.04324 [astro-ph.CO] .
  * Carr and Silk (2018) B. Carr and J. Silk, MNRAS 478, 3756 (2018), arXiv:1801.00672 [astro-ph.CO] .
  * Inman and Ali-Haïmoud (2019) D. Inman and Y. Ali-Haïmoud, Phys. Rev. D 100, 083528 (2019).
  * Hasinger (2020) G. Hasinger, arXiv e-prints , arXiv:2003.05150 (2020), arXiv:2003.05150 [astro-ph.CO] .
  * Meszaros (1974) P. Meszaros, Astron. Astrophys. 37, 225 (1974).
  * Meszaros (1975) P. Meszaros, Astron. Astrophys. 38, 5 (1975).
  * Lacki and Beacom (2010) B. C. Lacki and J. F. Beacom, Astrophys. J. 720, L67 (2010), arXiv:1003.3466 [astro-ph.CO] .
  * Afshordi  _et al._ (2003) N. Afshordi, P. McDonald, and D. N. Spergel, Astrophys. J. 594, L71 (2003), arXiv:astro-ph/0302035 [astro-ph] .
  * Murgia  _et al._ (2019) R. Murgia, G. Scelfo, M. Viel, and A. Raccanelli, Phys. Rev. Lett. 123, 071102 (2019), arXiv:1903.10509 [astro-ph.CO] .
  * Lamb (1975) H. Lamb,  _Hydrodynamics_ (Cambridge University Press, 1975).
  * Peebles (1969) P. J. E. Peebles, Astrophys. J. 155, 393 (1969).
  * Efstathiou and Jones (1979) G. Efstathiou and B. J. T. Jones, MNRAS 186, 133 (1979).
  * Hoyle (1953) F. Hoyle, Astrophys. J. 118, 513 (1953).
  * Rees and Ostriker (1977) M. J. Rees and J. P. Ostriker, MNRAS 179, 541 (1977).
  * Kashlinsky (1982) A. Kashlinsky, MNRAS 200, 585 (1982).
  * Vasil’kov (1975) A. P. Vasil’kov, Fluid Dynamics 10, 806 (1975).
  * Efstathiou (1990) G. Efstathiou, in  _Physics of the Early Universe_ (1990) pp. 361–463.
  * Clesse and García-Bellido (2018) S. Clesse and J. García-Bellido, Physics of the Dark Universe 22, 137 (2018), arXiv:1711.10458 [astro-ph.CO] .
  * Peacock and Dodds (1996) J. A. Peacock and S. J. Dodds, MNRAS 280, L19 (1996), arXiv:astro-ph/9603031 [astro-ph] .
  * Peebles (1980) P. J. E. Peebles,  _The large-scale structure of the universe_ (Princeton University Press, 1980).
  * Vittorio  _et al._ (1986) N. Vittorio, R. Juszkiewicz, and M. Davis, Nature (London) 323, 132 (1986).
  * Kashlinsky (1992) A. Kashlinsky, ApJL 386, L37 (1992).
  * Press and Schechter (1974) W. H. Press and P. Schechter, Astrophys. J. 187, 425 (1974).
  * Wu  _et al._ (2015) X.-B. Wu, F. Wang, X. Fan, W. Yi, W. Zuo, F. Bian, L. Jiang, I. D. McGreer, R. Wang, J. Yang, Q. Yang, D. Thompson, and Y. Beletsky, Nature (London) 518, 512 (2015), arXiv:1502.07418 [astro-ph.GA] .
  * Mortlock  _et al._ (2011) D. J. Mortlock, S. J. Warren, B. P. Venemans, M. Patel, P. C. Hewett, R. G. McMahon, C. Simpson, T. Theuns, E. A. Gonzáles-Solares, A. Adamson, S. Dye, N. C. Hambly, P. Hirst, M. J. Irwin, E. Kuiper, A. Lawrence, and H. J. A. Röttgering, Nature (London) 474, 616 (2011), arXiv:1106.6088 [astro-ph.CO] .
  * Bañados  _et al._ (2018) E. Bañados, B. P. Venemans, C. Mazzucchelli, E. P. Farina, F. Walter, F. Wang, R. Decarli, D. Stern, X. Fan, and F. B. Davies, Nature (London) 553, 473 (2018), arXiv:1712.01860 [astro-ph.GA] .
  * Yue  _et al._ (2013) B. Yue, A. Ferrara, R. Salvaterra, Y. Xu, and X. Chen, MNRAS 433, 1556 (2013), arXiv:1305.5177 .
  * Ricarte  _et al._ (2019) A. Ricarte, F. Pacucci, N. Cappelluti, and P. Natarajan, arXiv e-prints , arXiv:1907.03675 (2019), arXiv:1907.03675 [astro-ph.GA] .
  * Efstathiou and Rees (1988) G. Efstathiou and M. J. Rees, MNRAS 230, 5p (1988).
  * Kashlinsky (1993) A. Kashlinsky, ApJL 406, L1 (1993).
  * Tanaka  _et al._ (2013) T. L. Tanaka, M. Li, and Z. Haiman, MNRAS 435, 3559 (2013), arXiv:1309.2301 [astro-ph.CO] .
  * Bromm (2013) V. Bromm, Reports on Progress in Physics 76, 112901 (2013), arXiv:1305.5178 [astro-ph.CO] .
  * Latif and Ferrara (2016) M. A. Latif and A. Ferrara, Pub. Astron. Soc. Aust. 33, e051 (2016), arXiv:1605.07391 [astro-ph.GA] .
  * Begelman and Rees (1978) M. C. Begelman and M. J. Rees, MNRAS 185, 847 (1978).
  * Kashlinsky and Rees (1983) A. Kashlinsky and M. J. Rees, MNRAS 205, 955 (1983).
  * Begelman  _et al._ (2006) M. C. Begelman, M. Volonteri, and M. J. Rees, MNRAS 370, 289 (2006), arXiv:astro-ph/0602363 [astro-ph] .
  * Dunlop  _et al._ (1996) J. Dunlop, J. Peacock, H. Spinrad, A. Dey, R. Jimenez, D. Stern, and R. Windhorst, Nature (London) 381, 581 (1996).
  * Saracco  _et al._ (2009) P. Saracco, M. Longhetti, and S. Andreon, MNRAS 392, 718 (2009), arXiv:0810.2795 [astro-ph] .
  * Renzini (2006) A. Renzini, Ann. Rev. Astron. Astropys. 44, 141 (2006), arXiv:astro-ph/0603479 [astro-ph] .