# Sharing quantum nonlocality and genuine nonlocality with independent
observables

Tinggui Zhang1,† and Shao-Ming Fei2,3,♯ ${1}$ School of Mathematics and
Statistics, Hainan Normal University, Haikou, 571158, China
$2$ School of Mathematical Sciences, Capital Normal University, Beijing
100048, China
$3$ Shenzhen Institute for Quantum Science and Engineering, Southern
University of Science and Technology, Shenzhen 518055, China
† Correspondence to<EMAIL_ADDRESS>
♯ Correspondence to<EMAIL_ADDRESS>
###### Abstract

Recently the authors in [Phys. Rev. Lett. 125, 090401 (2020)] considered the
following scenario: Alice and Bob each have half of a pair of entangled qubit
state. Bob measures his half and then passes his part to a second Bob who
measures again and so on. The goal is to maximize the number of Bobs that can
have an expected violation of the Clauser-Horne-Shimony-Holt (CHSH) inequality
with the single Alice. By taking the maximally entangled pure two-qubit state
$|\phi\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$ as an example, it has
been constructively proved that arbitrarily many independent Bobs can share
the nonlocality with the single Alice. Here we demonstrate that arbitrarily
many independent observers can share the nonlocality of a single arbitrary
dimensional bipartite entangled but not necessary two-qubit entangled state.
Further, taking the generalized GHZ states as an example, we show that at most
two Charlies can share the genuine nonlocality of a single generalized GHZ
state with an Alice and a Bob.

###### pacs:

03.67.-a, 02.20.Hj, 03.65.-w

##  I Introduction

Quantum nonlocality, such as that revealed by violating the Bell inequalities
of quantum entangled states jsbe , is one of the most startling predictions of
quantum mechanics. Recently, as confirmed in loophole-free experiments bmlk ,
nonlocality has been proven to be useful in many quantum tasks such as device-
independent cryptography aanb  and randomness certification spaa ; rcol ; lkyz
; wlml . A successful and secure quantum network relies on quantum
correlations distributed and shared among many sites hjki . Different kinds of
multipartite quantum correlations have been considered as valuable resources
for various applications in quantum communication tasks. A key property is
that such quantum correlations cannot be freely shared among the multipartite
systems, see e.g. monogamy  and references therein.

Recently, in Ref rsng ; smam ; assd ; dass ; pjbr , the authors studied the
fundamental limits on nonlocality, asking whether a single pair of entangled
qubits could generate a long sequence of nonlocal correlations. This
sequential scenario (see FIG. 1) was introduced first in rsng . With the same
sharpness of the two measurements applied by each Bob, in Ref. smam  the
authors shown that at most two Bobs can achieve an expected CHSH chsh
violation with a single Alice, in line with the numerical evidence from rsng .
Equal sharpness two-outcome measurements were also adopted in assd  to show
that at most two Bobs can share the Bell nonlocality of a maximally entangled
state with a single Alice in this scenario. It has been shown that at most two
Bobs can exhibit bipartite nonlocality with a single Alice by using local
realist inequalities with three and four dichotomic measurements pre observer
dass . More recently, in pjbr  the authors studied such scenario and shown
that if the Bobs’ apply different measurements, then arbitrarily many
independent Bobs can share the nonlocality of the maximally entangled pure
two-qubit state $|\phi\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$ with
the single Alice.

Figure 1: A quantum state $\rho_{AB}^{(1)}$ is initially shared by Alice and
Bob(1). Bob(1) performs first a measurement on his part and then passes it to
Bob(2). The post-measurement state is $\rho_{AB}^{(2)}$. Bob(2) measures
$\rho_{AB}^{(2)}$ on his part and passes it to Bob(3) and so on.

High-dimensional quantum entanglement and nonlocality provide a playground for
fundamental research and also lead to technological advances, with stronger
locality violations that can be exploited to tolerate larger amounts of noise
in quantum communication protocols. Various physical and technical approaches
on how to manipulate multilevel quantum states in different degrees of freedom
have been presented, inspiring new synergies that create new technologies such
as teleporting the complete quantum information stored in a single “photon”
zeilinger . Therefore, in this article we study such nonlocal correlation
sharing scenario for arbitrary high dimensional bipartite entangled pure
states. We show that arbitrarily many independent observers can share the
nonlocality of any single arbitrary dimensional bipartite entangled states.
Furthermore, we investigate the genuine nonlocality sharing among Alice, Bob
and Charlies, see FIG. 2. By using the Svetlichny inequality svet , we show
that at least two Charlies can share the genuine nonlocality of a single
generalized GHZ state with Alice and Bob.

Figure 2:  A quantum state $\rho_{ABC}^{(1)}$ is initially shared by Alice,
Bob and Charlie(1). Charlie(1) performs a measurement on her part and passes
it to Charlie(2). The post-measurement state is $\rho_{ABC}^{(2)}$. Charlie(2)
measures $\rho_{ABC}^{(2)}$ on her part and passes it to Charlie(3) and so on.

##  II Nonlocal sharing of bipartite high-dimensional pure states

We first consider the measurement scenario such that Alice attempts to share
the nonlocal correlations of an entangled pure state with $n$ independent Bobs
(FIG. 1). Any bipartite pure state $|\psi\rangle\in H_{A}\otimes H_{B}$ with
$dim(H_{A})=s$ and $dim(H_{B})=t$ ($s\leq t$) has Schmidt decomposition form,
$|\psi\rangle=\sum_{i=1}^{s}c_{i}|i_{A}\rangle|i_{B}\rangle$, where
$c_{i}\in[0,1]$, $\sum_{i}^{s}c_{i}^{2}=1$, $\\{i_{A}\\}_{1}^{s}$ and
$\\{i_{B}\\}_{1}^{t}$ are the orthonormal bases of $H_{A}$ and $H_{B}$,
respectively. $|\psi\rangle$ is entangled if and only if at least two $c_{i}$s
are nonzero. Without loss of generality, below we assume that $c_{i}$ are
arranged in descending order.

To begin with, Alice and Bob(1) share an arbitrary entangled bipartite pure
state
$\rho_{AB}^{(1)}=|\psi\rangle\langle\psi|=\Sigma_{i,j}c_{i}c_{j}|ii\rangle\langle
jj|$. Bob(1) proceeds by choosing a uniformly random input, performing the
corresponding measurement and recording the outcome. Denote the binary input
and output of Alice (Bob(k)) by $X$ ($Y^{(k)}$) and $A$ ($B^{(k)}$),
respectively. Suppose Bob(1) performs the measurement according to $Y^{(1)}=y$
with the outcome $B^{(1)}=b$. Averaged over the inputs and outputs of Bob(1),
the unnormalized state shared between Alice and Bob(2) is given by

$\rho_{AB}^{(2)}=\frac{1}{2}\Sigma_{b,y}(I_{s}\otimes\sqrt{B_{b|y}^{(1)}})\rho_{AB}^{(1)}(I_{s}\otimes\sqrt{B_{b|y}^{(1)}})$

where $B^{(1)}_{b|y}$ is the positive operator-valued measure (POVM) effect
corresponding to outcome $b$ of Bob(1)’s measurement for input $y$, $I_{s}$ is
the $s\times s$ identity matrix. Repeating this process, one gets the state
$\rho_{AB}^{(k)}$ shared between Alice and Bob(k).

To detect the nonlocality we employ the CHSH inequality chsh ,
$I_{CHSH}=\langle\mathbb{B}\rangle\leq 2$, where
$\langle\mathbb{B}\rangle=Tr(\mathbb{B}\rho)$, $\mathbb{B}=A_{0}\otimes
B_{0}+A_{0}\otimes B_{1}+A_{1}\otimes B_{0}-A_{1}\otimes B_{1}$, $A_{i}$ and
$B_{i}$, $i=0,1$, are Hermitian operators with eigenvalues $\in[-1,1]$. If for
some binary observables $A_{i}$ and $B_{i}^{(k)}$, $i=0,1$,
$I_{CHSH}^{(k)}\equiv Tr(\mathbb{B}\rho_{AB}^{(k)})>2$, then the state
$\rho_{AB}^{(k)}$ is nonlocally correlated.

For the case that $s$ and $t$ are even, we employ the POVMs with measurement
operators $\\{E,I-E\\}$, where $E$ has the form
$E=\frac{1}{2}(I_{m}+\gamma(I_{\frac{m}{2}}\otimes{\vec{r}}\cdot\vec{\sigma}))$,
$\vec{r}\in R^{3}$ with $\|\vec{r}\|=1$,
$\vec{r}\cdot\vec{\sigma}=r_{1}\sigma_{1}+r_{2}\sigma_{2}+r_{3}\sigma_{3}$,
$\sigma_{i}$, $i=1,2,3$, are the standard Pauli matrices, $\gamma\in[0,1]$ is
the sharpness of the measurement, $I_{m}$ stands for the $m\times m$ identity
matrix, $m=s,t$. We set the Alice’s POVMs to be

$A_{0|0}=\frac{1}{2}(I_{s}+(I_{\frac{s}{2}}\otimes(\cos\theta\sigma_{3}+\sin\theta\sigma_{1})),$
(1)
$A_{0|1}=\frac{1}{2}(I_{s}+(I_{\frac{s}{2}}\otimes(\cos\theta\sigma_{3}-\sin\theta\sigma_{1})),$
(2)

for some $\theta\in(0,\frac{\pi}{4}]$. For each $k=1,2,\cdots,n$, Bob(k)’s
POVMs are defined as

$B^{(k)}_{0|0}=\frac{1}{2}(I_{t}+(I_{\frac{t}{2}}\otimes(\cos\theta\sigma_{3})),$
(3)
$B^{(k)}_{0|1}=\frac{1}{2}(I_{t}+(I_{\frac{t}{2}}\otimes(\gamma_{k}\sin\theta\sigma_{1})),$
(4)

$k=1,2,\cdots,n$.

When $s$ and $t$ are odd, we employ the POVM measurement operators
$\\{E,I-E\\}$, with

$E=\frac{1}{2}\left[I_{m}+\gamma\left(\begin{array}[]{cc}I_{[\frac{m}{2}]}\otimes{\vec{r}}\cdot\vec{\sigma}&0\\\
0&1\end{array}\right)\right],$

where $m=s,t$, $[p]$ represents the integer less or equal to $p$. The Alice’s
POVMs are defined as

$A_{0|0}=\frac{1}{2}\left[I_{s}+\left(\begin{array}[]{cc}I_{[\frac{s}{2}]}\otimes(\cos\theta\sigma_{3}+\sin\theta\sigma_{1})&0\\\
0&1\\\ \end{array}\right)\right],$ (5)
$A_{0|1}=\frac{1}{2}\left[I_{s}+\left(\begin{array}[]{cc}I_{[\frac{s}{2}]}\otimes(\cos\theta\sigma_{3}-\sin\theta\sigma_{1})&0\\\
0&1\\\ \end{array}\right)\right]$ (6)

for some $\theta\in(0,\frac{\pi}{4}]$. The Bob(k)’s POVMs set to be

$B^{(k)}_{0|0}=\frac{1}{2}\left[I_{t}+\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma_{3}&0\\\
0&1\\\ \end{array}\right)\right],$ (7)
$B^{(k)}_{0|1}=\frac{1}{2}\left[I_{t}+\gamma_{k}\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma_{1}&0\\\
0&1\\\ \end{array}\right)\right],$ (8)

$k=1,2,\cdots,n$.

The observables are then given by $A_{x}=A_{0|x}-A_{1|x}$ and
$B_{y}^{(k)}=B_{0|y}-B_{1|y}$, $x,y=0,1$, for both even and odd $t,s$. We have
the following conclusion for the expected CHSH value for Alice and Bob(k), see
proof in Appendix.

######  Theorem 1.

For any initial entangled bipartite pure quantum state $|\psi\rangle\in
H_{A}\otimes H_{B}$ with Schmidt decomposition
$|\psi\rangle=\sum^{s}_{i=1}c_{i}|i_{A}\rangle|i_{B}\rangle$, the expected
CHSH value of $\rho_{AB}^{(k)}$ is given by

$I^{(k)}_{CHSH}\geq
2^{2-k}\left[\gamma_{k}L\sin\theta+\cos\theta\prod_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}^{2}})\right],$
(9)

where
$L=2(c_{1}c_{2}+c_{3}c_{4}+\cdots+c_{2[\frac{s}{2}]-1}c_{2[\frac{s}{2}]})$.

Next we show that there exist suitable parameters $\gamma_{k}$ and $\theta$
such that $I_{CHSH}^{(k)}>2$ for arbitrary $k$. From Theorem 1 we require that

$\displaystyle\gamma_{k}>\frac{2^{k-1}-\cos\theta\prod_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}^{2}})}{L\sin\theta}.$
(10)

Set $\gamma_{1}(\theta)=(1+\epsilon)\frac{1-\cos(\theta)}{L\sin(\theta)}$ for
$\varepsilon>0$. One has

$\displaystyle\gamma_{k}(\theta)=(1+\epsilon)\frac{2^{k-1}-\cos\theta\prod_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}^{2}})}{L\sin\theta}$
(11)

for $\gamma_{k-1}(\theta)\in(0,1)$, $k\in\\{2,\cdots,n\\}$.

######  Theorem 2.

For each $n\in\mathbb{N}$, there exists a sequence $\\{\gamma_{k}\\}_{1}^{n}$
and $\theta_{n}\in(0,\frac{\pi}{4}]$ such that $I_{CHSH}^{(k)}>2$ for all
$k=1,2,\cdots,n$, with $\theta\in(0,\theta_{n})$ and $\gamma_{k}(\theta)<1$
for all $k\leq n$.

The proof is given in Appendix. Theorem 2 shows that arbitrarily many
independent observers can share the nonlocality of a single arbitrary
dimensional bipartite entangled state, as long as at least two Schmidt
coefficients $c_{i}$ of the state are nonzero. Moreover, the state is not
necessary to be maximally entangled.

##  III Genuine nonlocal sharing of three-qubit states

Multipartite nonlocal correlations have not only foundational implications
ngis  but also novel applications in quantum communication and quantum
computation msju ; msee ; pzet ; cylu , as well as in phase transitions and
criticality in many-body systems pzet . Essentially different from the
bipartite case, one has so called genuine multipartite nonlocality for
multipartite systems. In the following, we consider the sharing ability of the
genuine nonlocality of the three-qubit generalized GHZ state,
$|\psi_{\alpha}\rangle=\cos\alpha|000\rangle+\sin\alpha|111\rangle.$

Let $A_{i}=\vec{a}_{i}\cdot\vec{\sigma}$, $B_{i}=\vec{b}_{i}\cdot\vec{\sigma}$
and $C_{i}=\vec{c}_{i}\cdot\vec{\sigma}$, $i=0,1$, be the measurement
observables on the first, second and third qubit, respectively, with
$\vec{a}_{i}$, $\vec{b}_{i}$ and $\vec{c}_{i}$ the real unit vectors. The
Svetlichny operator is defined by svet ,

$\displaystyle S=A_{0}(B_{0}+B_{1})C_{0}+A_{0}(B_{0}-B_{1})C_{1}$
$\displaystyle+A_{1}(B_{0}-B_{1})C_{0}-A_{1}(B_{0}+B_{1})C_{1}.$ (12)

If a state $|\psi\rangle$ admits bi-local hidden variable model, then the
expectation value of the Svetlichny operator satisfies the Svetlichny’s
inequality,

$\displaystyle\mathbb{S}(\psi)=\langle\psi|S|\psi\rangle\leq 4.$ (13)

If (13) is violated, $|\psi\rangle$ must be a genuine three-qubit nonlocally
correlated state. The maximal violation of the Svetlichny inequality (13) for
the state $|\psi_{g}\rangle$ has been studied in sgns . It has been shown that
when $\sin^{2}2\alpha>\frac{1}{2}$, the state $|\psi_{g}\rangle$ is genuine
three-qubit nonlocal.

We consider the following measurement scenario, see FIG. 2: Alice and Bob want
to share the genuine three-qubit nonlocality of single three-qubit state
$|\psi_{g}\rangle$ with possible $n$ independent Charlies. Denote the binary
input and output of Alice (Bob) by $X$ ($Y$) and $A$ ($B$), respectively. For
each $k\in{N}$ we denote the binary input and output of Charlie(k) by
$Z^{(k)}$ and $C^{(k)}$, respectively. At the beginning, the three-qubit state
$|\psi_{g}\rangle$ is shared among Alice, Bob and Charlie(1). Charlie(1)
proceeds by choosing a uniformly random input, performing the corresponding
measurement and recording the outcome. The postmeasurement qubit is then sent
to Charlie(2). Suppose Charlie(1) performs the measurement according to
$Z^{(1)}=z$ and receives the outcome $C^{(1)}=c$. The postmeasurement state
can be described by the Lüders rule. Averaged over the inputs and outputs of
Charlie(1)’s, the postmeasurement unnormalized state $\rho_{ABC}^{(2)}$ shared
among Alice, Bob and Charlie(2) is given by

$\displaystyle\rho_{ABC}^{(2)}=\frac{1}{2}\sum_{c,z}(I_{2}\otimes
I_{2}\otimes\sqrt{C_{c|z}^{(1)}})\rho_{ABC}^{(1)}(I_{2}\otimes
I_{2}\otimes\sqrt{C_{c|z}^{(1)}}),$

where $C^{(1)}_{c|z}$ is the POVM effect corresponding to the outcome $c$ of
Charlie(1)’s measurement for input $z$. Repeating this process, one can
compute the state $\rho_{ABC}^{(k)}$ shared among Alice, Bob and Charlie(k).
The expected value of Svetlichny operator associated with the state
$\rho_{ABC}^{(k)}$ and the binary random observables $A_{i}$, $B_{i}$ and
$C_{i}^{(k)}$, $i=0,1$, is given by $S^{(k)}\equiv Tr(S\rho_{ABC}^{(k)})$.

To see the maximal $k$ such that $S^{(k)}>4$, we consider again two-outcome
POVMs $\\{E,I-E\\}$. Let Alice’s POVMs be given by

$A_{0|0}=\frac{1}{2}(I+\sigma_{1}),~{}~{}~{}A_{0|1}=\frac{1}{2}(I+\sigma_{2})$

and Bob’s POVMs by

$\begin{array}[]{rcl}B_{0|0}&=&\displaystyle\frac{1}{2}(I+\cos\theta\sigma_{1}-\sin\theta\sigma_{2}),\\\\[5.69054pt]
B_{0|1}&=&\displaystyle\frac{1}{2}(I+\cos\theta\sigma_{1}+\sin\theta\sigma_{2})\end{array}$

for some $\theta\in(0,\frac{\pi}{4}]$. The Charlie(k)’s POVMs are defined by

$C^{(k)}_{0|0}=\frac{1}{2}(I+\sigma_{1}),~{}~{}~{}C^{(k)}_{0|1}=\frac{1}{2}(I+\gamma_{k}\sigma_{2})$

for $k=1,2,\cdots,n$.

Set $A_{x}=A_{0|x}-A_{1|x}$, $B_{y}=B_{0|y}-B_{1|y}$ and
$C^{(k)}_{z}=C^{(k)}_{0|z}-C^{(k)}_{1|z}$, $x,y,z=0,1$. We have the following
conclusion, see detailed proof in Appendix.

######  Theorem 3.

For the initially shared generalized GHZ state $|\psi_{\alpha}\rangle$, the
expected value of the Svetlichny operator with respect to the state
$\rho_{ABC}^{(k)}$ is given by

$S^{(k)}=2^{2-k}\sin
2\alpha(\cos\theta+\sin\theta)(\gamma_{k}+\prod_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}^{2}})).$
(14)

For $k=1$, we have $S^{(1)}=2\sin
2\alpha(\cos\theta+\sin\theta)(\gamma_{k}+1)$. If $\gamma_{k}=1$ and
$\theta=\frac{\pi}{4}$, then $S^{(1)}>4$ as long as
$\sin^{2}(2\alpha)>\frac{1}{2}$. By detailed analysis, we have, see proof in
Appendix,

######  Theorem 4.

There are at most two Charlies sharing the genuine nonlocality with Alice and
Bob, $S^{(k)}>4$, $k=1,2$, when $\sin^{2}2\alpha\in(\frac{8}{9},1]$ for the
initially shared generalized GHZ state $|\psi_{\alpha}\rangle$.

In sdsd  the authors considered this problem from the formalism of weak or
unsharp measurements for the GHZ stat,
$|GHZ\rangle=\frac{1}{\sqrt{2}}(|000\rangle+|111\rangle)$. Our conclusion is
for the generalized GHZ states: when $\sin^{2}2\alpha\in(\frac{8}{9},1]$, at
most two Charlies can share the three-qubit genuine nonlocality of a single
generalized GHZ state with an Alice and a Bob, which coincides with the result
derived in [24] as a special case ($\alpha=\frac{\pi}{4}$).

##  IV Conclusions and discussions

Quantum nonlocality is a fundamental feature in quantum mechanics. We have
demonstrated that it is possible for arbitrarily many independent Bobs to
violate the CHSH inequality with a single Alice by using any bipartite high-
dimensional pure states that are either maximally or non-maximally entangled.
As the approach used for qubit case can not be used for high dimensional case
since relations like
$Tr(\rho(\sigma_{\vec{a}}\otimes\sigma_{\vec{b}}))=(\vec{a},T(\rho)\vec{b})$
pjbr  is only correct for two-qubit states, where $T(\rho)$ is a matrix with
entries given by $T_{i,j}(\rho)=Tr[\rho(\sigma_{i}\otimes\sigma_{j})]$, we
have presented a new approach in derivations. Our innovation also lies in
choosing the POVM measurement operators and in calculating the expected CHSH
values. Our approach can be also extended to the case of mixed states.

We have also investigated the shareability of genuine tripartite nonlocality.
For the generalized GHZ state, it has been shown that from our measurement
schemes two Charlies can share the genuine nonlocality with Alice and Bob. Our
results may also highlight researches on sharing general multipartite quantum
nonlocalities and other quantum correlations such as quantum steerability sdsa
; sgag , entanglement asau ; adaa  and coherence sasm .

Finally, in this article we have constructed the higher-dimensional dichotomic
POVM measurement operators in terms of the Pauli operators. As the Pauli
operators are easily implemented in experiments, the POVM operators we
constructed may have potential advantages in some specific experimental
implementations mzxc ; mlmg . Our approach may also highlight the related
applications in randomness generation fmrm , quantum teleportation sasa  and
random access codes kanb . It would be also interesting to give some insight
into larger multipartite scenarios and explore the relationship between our
methodology and, for instance, nonlocality depth sharing fjan  and monogamy
relations rmmj .

Acknowledgments: This work is supported by the National Natural Science
Foundation of China under Grant Nos. 11861031, 11675113 and 12075159, Beijing
Municipal Commission of Education under grant No. KZ201810028042, Beijing
Natural Science Foundation (Z190005), the Education Department of Hainan
Province, project number Hnky2020ZD-10, Academy for Multidisciplinary Studies,
Capital Normal University and Shenzhen Institute for Quantum Science and
Engineering, Southern University of Science and Technology (Grant No.
SIQSE202005).

## References

  * (1) J. S. Bell, Physics, 1, 195 (1964).
  * (2) B. Hensen, et al, Nature, 526, 682 (2015);
M. Giustina, et al, Phys. Rev. Lett 115, 250401, (2015);
L. K. Shalm, et al, Phys. Rev. Lett. 115, 250402, (2015).

  * (3) A. Acín, N. Brunner, N. Gisin, S. Massar, S. Pironio, and V. Scarani, Phys. Rev. Lett. 98, 230501 (2007).
  * (4) S. Pironio, et al, Nature 464, 1021 (2010).
  * (5) R. Colbeck, arXiv:0911.3814 [quant-ph] (2019).
  * (6) L. K. Shalm, et al, arXiv:1912.11158 [quant-ph] (2019).
  * (7) W.Z. Liu, et al, arXiv:1912.11159 [quant-ph] (2019).
  * (8) H.J. Kimble, Nature (London) 453, 1023 (2008).
  * (9) Z.X. Jin, S.M. Fei, Phys. Rev. A 99, 032343 (2019).
  * (10) R. Silva, N. Gisin, Y. Guryanova, and S. Popescu, Phys. Rev. Lett. 114, 250401 (2015).
  * (11) S. Mal, A. Majumdar, and D. Home, Mathematics 4, 48 (2016).
  * (12) A. Shenoy H., S. Designolle, F. Hirsch, R. Silva, N. Gisin, and N. Brunner, Phys. Rev. A 99, 022317 (2019).
  * (13) D. Das, A. Ghosal, S. Sasmal, S. Mal, and A. S. Majumdar, Phys. Rev. A 99, 022305 (2019).
  * (14) P. J. Brown and R. Colbeck, Phys. Rev. Lett. 125, 090401 (2020).
  * (15) J. F. Clauser, M.A. Horne, A. Shimony and R.A. Holt, Phys. Rev. Lett. 23, 880 (1969).
  * (16) M. Erhard, M. Krenn, A. Zeilinger, Nature Rev. Phys. 2, 365-381 (2020).
  * (17) G. Svetlichny, Phys. Rev. D 35, 3066 (1987).
  * (18) N. Gisin, arXiv:quant-ph/0702021v2, (2007).
  * (19) M. Seevinck and J. Uffink, Phys. Rev. A 78, 032101 (2008); M. Hein, J. Eisert, and H. J. Briegel, Phys. Rev. A 69, 062311 (2004); K. Chen and H.-K. Lo, Quantum Inf. Comput. 7, 689 (2007); V. Scarani and N. Gisin, Phys. Rev. Lett. 87, 117901 (2001).
  * (20) M. Seevinck, Ph.D. thesis, Utrecht University, arXiv:0811.1027v2 [quant-ph].
  * (21) P. Zoller et al., Eur. Phys. J. D 36, 203 (2005).
  * (22) C. Y. Lu et al., Nature Phys. 3, 91 (2007).
  * (23) S. Ghose, N. Sinclair, S. Debnath, P. Rungta, and R. Stock Phys. Rev. Lett. 102, 250404 (2009);
V. Gupta, arxiv:1408.3073v1 [quant-ph].

  * (24) S. Saha, D. Das, S. Sasmal, D. Sarkar, K. Mukherjee, A. Roy, S. S. Bhattacharya, Quant. Inf. Process. 18, 42 (2019).
  * (25) S. Sasmal, D. Das, S. Mal, and A. S. Majumdar, Phys. Rev. A 98, 012305 (2018).
  * (26) S. Gupta, A. G. Maity, D. Das, A. Roy, A. S. Majumdar, arXiv:2007.03552 [quant-ph].
  * (27) A. Bera, S. Mal, A. Sen De, and U. Sen, Phys. Rev. A 98, 062304 (2018).
  * (28) A. G. Maity, D. Das, A. Ghosal, A. Roy, and A. S. Majumdar, Phys. Rev. A 101, 042340 (2020).
  * (29) S. Datta, and A. S. Majumdar, Phys. Rev. A 98, 042311 (2018).
  * (30) M. J. Hu, Z. Y. Zhou, X. M. Hu, C. F. Li, G. C. Guo, and Y. S. Zhang, npj Quantum Inf. 4, 63 (2018).
  * (31) M. Schiavon, L. Calderaro, M. Pittaluga, G. Vallone, and P. Villoresi, Quantum Sci. Technol. 2, 015010 (2017).
  * (32) F. J. Curchod, M. Johansson, R. Augusiak, M. J. Hoban, P. Wittek, and A. Acin, Phys. Rev. A 95, 020102 (2017).
  * (33) S. Roy, A. Bera, S. Mal, A. Sen De, and U. Sen, arXiv:1905.04164 [quant-ph].
  * (34) K. Mohan, A. Tavakoli, and N. Brunner, New J. Phys. 21, 083034 (2019).
  * (35) F. Baccari, J. Tura, M. Fadel, A. Aloy, J.-D. Bancal, N. Sangouard, M. Lewenstein, A. Acín, and R. Augusiak, Phys. Rev. A 100, 022121 (2019).
  * (36) R. Augusiak, M. Demianowicz, M. Pawlowski, J. Tura, and A. Acín, Phys. Rev. A 90, 052323 (2014).

## APPENDIX

###  IV.1 Proof of Theorem 1

We first consider the simplest case that both $t$ and $s$ are even. From (1),
(2), (3) and (4) we have

$A_{0}=A_{0|0}-A_{1|0}=I_{\frac{s}{2}}\otimes(\cos\theta\sigma_{3}+\sin\theta\sigma_{1}),$
$A_{1}=A_{0|1}-A_{1|1}=I_{\frac{s}{2}}\otimes(\cos\theta\sigma_{3}-\sin\theta\sigma_{1})$

and

$B_{0}^{(k)}=B_{0|0}-B_{1|0}=I_{\frac{t}{2}}\otimes\sigma_{3},$
$B_{1}^{(k)}=B_{0|1}-B_{1|1}=I_{\frac{t}{2}}\otimes\gamma_{k}\sigma_{1},$

we have

$\displaystyle I_{CHSH}^{(k)}$ (A1) $\displaystyle=$ $\displaystyle
Tr[\rho_{AB}^{(k)}((A_{0}+A_{1})\otimes B_{0}^{(k)})]$
$\displaystyle+Tr[\rho_{AB}^{(k)}((A_{0}-A_{1})\otimes B_{1}^{(k)})]$
$\displaystyle=$ $\displaystyle 2\cos\theta\
Tr[\rho_{AB}^{(k)}(I_{\frac{s}{2}}\otimes\sigma_{3})\otimes(I_{\frac{t}{2}}\otimes\sigma_{3})]$
$\displaystyle+2\gamma_{k}\sin\theta\
Tr[\rho_{AB}^{(k)}(I_{\frac{s}{2}}\otimes\sigma_{1})\otimes(I_{\frac{t}{2}}\otimes\sigma_{1})].$

For the case $k=1$, we get

$\displaystyle I_{CHSH}^{(1)}$ (A2) $\displaystyle=$ $\displaystyle
2\cos\theta\
Tr[\rho_{AB}^{(1)}(I_{\frac{s}{2}}\otimes\sigma_{3})\otimes(I_{\frac{t}{2}}\otimes\sigma_{3})]$
$\displaystyle+2\gamma_{1}\sin\theta\
Tr[\rho_{AB}^{(1)}(I_{\frac{s}{2}}\otimes\sigma_{1})\otimes(I_{\frac{t}{2}}\otimes\sigma_{1})]$
$\displaystyle=$ $\displaystyle
2\cos\theta(\Sigma_{i}c_{i}^{2})+2\gamma_{1}(2(c_{1}c_{2}+c_{3}c_{4}+\cdots))\sin\theta$
$\displaystyle=$ $\displaystyle 2\cos\theta+2\gamma_{1}L\sin\theta,$

where we have used $\rho_{AB}^{(1)}=\Sigma_{i,j}c_{i}c_{j}|ii\rangle\langle
jj|$ and $L=2(c_{1}c_{2}+c_{3}c_{4}+\cdots)$.

To obtain the relation between $\rho_{AB}^{(1)}$ and$\rho_{AB}^{(k)}$ we use
the following equation,

$\displaystyle\sqrt{\frac{1}{2}(I_{t}\pm\gamma_{k}I_{\frac{m}{2}}\otimes\sigma)}=$
$\displaystyle\frac{(\sqrt{1+\gamma_{k}}+\sqrt{1-\gamma_{k}})I_{t}\pm(\sqrt{1+\gamma_{k}}-\sqrt{1-\gamma_{k}})(I_{\frac{m}{2}}\otimes\sigma)}{2\sqrt{2}}$

for $\sigma=\sigma_{1},\sigma_{3}$. By using the Lüders update rule we have

$\displaystyle\rho_{AB}^{(k)}=\frac{1}{2}\Sigma_{b,y}(I_{s}\otimes\sqrt{B_{b|y}^{(k-1)}})\rho_{AB}^{(k-1)}(I_{s}\otimes\sqrt{B_{b|y}^{(k-1)}})$
(A3) $\displaystyle=$
$\displaystyle\frac{1}{4}(I_{s}\otimes(I_{\frac{t}{2}}\otimes\sigma_{3}))\rho_{AB}^{(k-1)}(I_{s}\otimes(I_{\frac{t}{2}}\otimes\sigma_{3}))$
$\displaystyle+$
$\displaystyle\frac{1-\sqrt{1-\gamma_{k-1}^{2}}}{4}(I_{s}\otimes(I_{\frac{t}{2}}\otimes\sigma_{1}))\rho_{AB}^{(k-1)}(I_{s}\otimes(I_{\frac{t}{2}}\otimes\sigma_{1}))$
$\displaystyle+$
$\displaystyle\frac{2+\sqrt{1-\gamma_{k-1}^{2}}}{4}\rho_{AB}^{(k-1)}.$

Substituting (A3) into (A1) and taking into account the relations
$\sigma_{3}\sigma_{3}\sigma_{3}=\sigma_{3}$ and
$\sigma_{1}\sigma_{3}\sigma_{1}=-\sigma_{3}$, we get

$\displaystyle
Tr[\rho_{AB}^{(k)}(I_{\frac{s}{2}}\otimes\sigma_{3})\otimes(I_{\frac{t}{2}}\otimes\sigma_{3})]$
$\displaystyle=$
$\displaystyle(\frac{1}{4}-\frac{1-\sqrt{1-\gamma_{k-1}^{2}}}{4}+\frac{2+\sqrt{1-\gamma_{k-1}^{2}}}{4})$
$\displaystyle\cdot$ $\displaystyle
Tr[\rho_{AB}^{(k-1)}(I_{\frac{s}{2}}\otimes\sigma_{3})\otimes(I_{\frac{t}{2}}\otimes\sigma_{3})]$
$\displaystyle=$ $\displaystyle\frac{1+\sqrt{1-\gamma_{k-1}^{2}}}{2}$
$\displaystyle\cdot$ $\displaystyle
Tr[\rho_{AB}^{(k-1)}(I_{\frac{s}{2}}\otimes\sigma_{3})\otimes(I_{\frac{t}{2}}\otimes\sigma_{3})].$

Similarly we can obtain

$\displaystyle
Tr[\rho_{AB}^{(k)}(I_{\frac{s}{2}}\otimes\sigma_{1})\otimes(I_{\frac{t}{2}}\otimes\sigma_{1})]$
$\displaystyle=$
$\displaystyle\frac{1}{2}Tr[\rho_{AB}^{(k-1)}(I_{\frac{s}{2}}\otimes\sigma_{1})\otimes(I_{\frac{t}{2}}\otimes\sigma_{1})].$

By recursion, we get

$\displaystyle
Tr[\rho_{AB}^{(k)}(I_{\frac{s}{2}}\otimes\sigma_{3})\otimes(I_{\frac{t}{2}}\otimes\sigma_{3})]=$
$\displaystyle
2^{1-k}\prod_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}^{2}})Tr[\rho_{AB}^{(1)}(I_{\frac{s}{2}}\otimes\sigma_{3})$
$\displaystyle\otimes(I_{\frac{t}{2}}\otimes\sigma_{3})]$ (A4)

and

$\displaystyle
Tr[\rho_{AB}^{(k)}(I_{\frac{s}{2}}\otimes\sigma_{1})\otimes(I_{\frac{t}{2}}\otimes\sigma_{1})]=$
$\displaystyle
2^{1-k}Tr[\rho_{AB}^{(1)}(I_{\frac{s}{2}}\otimes\sigma_{1})\otimes(I_{\frac{t}{2}}\otimes\sigma_{1})].$
(A5)

Substituting (IV.1) and (IV.1) into (A1) and using (A2), we obtain (9).

Now we prove the most complex case that both $s$ and $t$ are odd. From (5),
(6), (7) and (8) we have

$A_{0}=A_{0|0}-A_{1|0}=\left(\begin{array}[]{cc}I_{[\frac{s}{2}]}\otimes(\cos\theta\sigma_{3}+\sin\theta\sigma_{1})&0\\\
0&1\end{array}\right),$
$A_{1}=A_{0|1}-A_{1|1}=\left(\begin{array}[]{cc}I_{[\frac{s}{2}]}\otimes(\cos\theta\sigma_{3}-\sin\theta\sigma_{1})&0\\\
0&1\end{array}\right),$
$B_{0}^{(k)}=B_{0|0}-B_{1|0}=\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma_{3}&0\\\
0&1\end{array}\right)$

and

$B_{1}^{(k)}=B_{0|1}-B_{1|1}=\gamma_{k}\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma_{1}&0\\\
0&1\end{array}\right).$

Correspondingly, we have

$\displaystyle I_{CHSH}^{(k)}$ (A15) $\displaystyle=$ $\displaystyle
Tr[\rho_{AB}^{(k)}((A_{0}+A_{1})\otimes B_{0}^{(k)})]$
$\displaystyle+Tr[\rho^{(k)}((A_{0}-A_{1})\otimes B_{1}^{(k)})]$
$\displaystyle=$ $\displaystyle
2Tr\left[\rho_{AB}^{(k)}\left(\begin{array}[]{cc}I_{[\frac{s}{2}]}\otimes(\cos\theta\sigma_{3})&0\\\
0&1\\\
\end{array}\right)\otimes\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma_{3}&0\\\
0&1\\\ \end{array}\right)\right]$ $\displaystyle+2\gamma_{k}\sin\theta
Tr\left[\rho_{AB}^{(k)}\left(\begin{array}[]{cc}I_{[\frac{s}{2}]}\otimes\sigma_{1}&0\\\
0&0\\\
\end{array}\right)\otimes\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma_{1}&0\\\
0&1\\\ \end{array}\right)\right].$

For the case $k=1$, one has

$\displaystyle I_{CHSH}^{(1)}$ (A25) $\displaystyle=$ $\displaystyle
2Tr[\rho_{AB}^{(1)}\left(\begin{array}[]{cc}I_{[\frac{s}{2}]}\otimes(\cos\theta\sigma_{3})&0\\\
0&1\\\
\end{array}\right)\otimes\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma_{3}&0\\\
0&1\\\ \end{array}\right)]$ $\displaystyle+2\gamma_{k}\sin\theta
Tr[\rho_{AB}^{(1)}\left(\begin{array}[]{cc}I_{[\frac{s}{2}]}\otimes\sigma_{1}&0\\\
0&0\\\
\end{array}\right)\otimes\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma_{1}&0\\\
0&1\\\ \end{array}\right)]$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{lr}2\cos\theta(1-c_{s}^{2})+2c_{s}^{2}+2\gamma_{1}L\sin\theta,\
\ \mbox{if}\ s=t\ \mbox{and}\ c_{s}\neq 0,\\\\[2.84526pt]
2\cos\theta+2\gamma_{1}L\sin\theta,\ \ \mbox{otherwise}\end{array}\right.$
$\displaystyle\geq$ $\displaystyle 2\cos\theta+2\gamma_{1}L\sin\theta,$

where we have used
$\rho_{AB}^{(1)}=\Sigma_{i,j=1}^{s}c_{i}c_{j}|ii\rangle\langle jj|$ and
$L=2(c_{1}c_{2}+c_{3}c_{4}+\cdots+c_{2[\frac{s}{2}]-1}c_{2[\frac{s}{2}]})$.

Using the following identity

$\displaystyle\sqrt{\frac{1}{2}(I_{t}+\gamma_{k}\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma&0\\\
0&1\\\ \end{array}\right))}$ $\displaystyle=$
$\displaystyle\frac{(\sqrt{1+\gamma_{k}}+\sqrt{1-\gamma_{k}})I_{t}}{2\sqrt{2}}$
$\displaystyle\pm$
$\displaystyle\frac{(\sqrt{1+\gamma_{k}}-\sqrt{1-\gamma_{k}})\left(\begin{array}[]{cc}I_{[\frac{t}{2}]}\otimes\sigma&0\\\
0&1\\\ \end{array}\right)}{2\sqrt{2}}$ (A34)

and repeating the similar process for even $t$ and $s$, we get the relations
similar to Eq. (IV.1) and Eq.(IV.1). Moreover, we get

$I^{(k)}_{CHSH}\geq
2^{2-k}(\gamma_{k}L\sin\theta+\cos\theta\prod_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}^{2}})).$

It is straightforward to prove that the above inequality holds also for the
cases of even (odd) $t$ and odd (even) $s$.

###  IV.2 Proof of Theorem 2

The Theorem can be proved by using two lemmas in pjbr . Since the sequence
$\gamma_{k}(\theta)\\}_{k\in\mathbb{N}}$ defined by (11) reduces to the ones
given in pjbr  when $L=\lambda_{1}$, we have that
$\\{\gamma_{1}(\theta)=(1+\epsilon)\frac{1-\cos(\theta)}{L\sin(\theta)}$ is
positive and an increasing function of $\theta$ for
$\theta\in(0,\frac{\pi}{4}]$ and $\varepsilon>0$. Moreover, the subsequence
$\\{\gamma_{k}(\theta)\\}_{k\in\mathbb{N}}$ consisting of all finite terms is
a strictly increasing sequence. In addition to the sequence being
monotonically increasing, each term in the sequence also has a vanishing limit
as $\theta$ approaches $0$. For any $n\in\mathbb{N}$ there exists some
$\theta_{n}\in(0,\frac{\pi}{4}]$ such that for all $k\leq n$ and
$\theta\in(0,\theta_{n})$, $\gamma_{k}(\theta)<1$. Moreover, we have
$\lim_{\theta\rightarrow 0^{+}}\gamma_{n}(\theta)=0$ for all $n\in\mathbb{N}$.

Therefore, there exists some $\theta_{n}\in(0,\frac{\pi}{4}]$ such that
$\gamma_{k}(\theta)<1$ and
$0<\gamma_{1}(\theta_{n})<\gamma_{2}(\theta_{n})<\cdots<\gamma_{n}(\theta_{n})$,
with each $\gamma_{k}(\theta_{n})$ satisfying the condition (10) and thus
giving rise to Bell violations.

###  IV.3 Proof of Theorem 3

First of all, we have

$\begin{array}[]{l}A_{0}=A_{0|0}-A_{1|0}=\sigma_{1},\\\
A_{1}=A_{0|1}-A_{1|1}=\sigma_{2},\end{array}$
$\begin{array}[]{l}B_{0}=B_{0|0}-B_{1|0}=\cos\theta\sigma_{1}-\sin\theta\sigma_{2},\\\
B_{1}=B_{0|1}-B_{1|1}=\cos\theta\sigma_{1}+\sin\theta\sigma_{2}\end{array}$

and

$\begin{array}[]{l}C_{0}^{(k)}=C_{0|0}-C_{1|0}=\sigma_{1},\\\
C_{1}^{(k)}=C_{0|1}-C_{1|1}=\gamma_{k}\sigma_{2}.\end{array}$

Then

$\displaystyle S^{(k)}=Tr[\mathbb{S}\rho_{ABC}^{(k)}]$ (A35) $\displaystyle=$
$\displaystyle 2\cos\theta
Tr[\rho_{ABC}^{(k)}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}]$
$\displaystyle-$ $\displaystyle 2\sin\theta
Tr[\rho_{ABC}^{(k)}\sigma_{1}\otimes\sigma_{2}\otimes\sigma_{2}]$
$\displaystyle-$ $\displaystyle
2\sin\theta\gamma_{k}Tr[\rho_{ABC}^{(k)}\sigma_{2}\otimes\sigma_{2}\otimes\sigma_{1}]$
$\displaystyle-$ $\displaystyle
2\cos\theta\gamma_{k}Tr[\rho_{ABC}^{(k)}\sigma_{2}\otimes\sigma_{1}\otimes\sigma_{2}].$

For the case $k=1$, substituting
$\rho_{ABC}^{1}=|\psi_{g}\rangle\langle\psi_{g}|$ into (A35) we get

$\displaystyle S^{(1)}$ $\displaystyle=$ $\displaystyle 2\cos\theta\sin
2\alpha+2\sin\theta\sin 2\alpha$ $\displaystyle+2\sin\theta\gamma_{1}\sin
2\alpha+2\cos\theta\gamma_{1}\sin 2\alpha$ $\displaystyle=$ $\displaystyle
2\sin 2\alpha(\cos\theta+\sin\theta)(1+\gamma_{1}).$

By using the following identity

$\displaystyle\sqrt{\frac{1}{2}(I\otimes I+\gamma_{k}\sigma)}$
$\displaystyle=$
$\displaystyle\frac{(\sqrt{1+\gamma_{k}}+\sqrt{1-\gamma_{k}})I\otimes
I\pm(\sqrt{1+\gamma_{k}}-\sqrt{1-\gamma_{k}})\sigma}{2\sqrt{2}}$

for $\sigma=\sigma_{1},\sigma_{2}$ and repeating the same processes in the
proof of Theorem 1, we obtain

$\displaystyle
Tr[\rho_{ABC}^{(k)}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}]$
$\displaystyle=$
$\displaystyle(\frac{1}{4}-\frac{1-\sqrt{1-\gamma^{2}_{k-1}}}{4}+\frac{2+\sqrt{1-\gamma^{2}_{k-1}}}{4})$
$\displaystyle\cdot$ $\displaystyle
Tr[\rho_{ABC}^{(k-1)}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}]$
$\displaystyle=$
$\displaystyle(\frac{1+\sqrt{1+\gamma^{2}_{k-1}}}{2})Tr[\rho_{ABC}^{(k-1)}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}].$

Then

$\displaystyle
Tr[\rho_{ABC}^{(k)}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}]=$
$\displaystyle
2^{1-k}Tr[\rho_{ABC}^{(1)}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}]\Pi_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}^{2}}).$

At last we have

$S^{(k)}=2^{2-k}\sin
2\alpha(\cos\theta+\sin\theta)(\gamma_{k}+\prod_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}^{2}})).$

###  IV.4 Proof of Theorem 4

From (14) Alice, Bob and Charlie share a genuine nonlocally correlated state
if

$\displaystyle\gamma_{k}>\frac{2^{k}}{\sin
2\alpha(\cos\theta+\sin\theta)}-\prod_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}^{2}}).$

Let $\varepsilon>0$ and $\gamma_{1}(\theta)=(1+\epsilon)(\frac{2}{\sin
2\alpha(\cos\theta+\sin\theta)}-1)$. For $k\in\\{2,\cdots,n\\}$ recursively
set

$\gamma_{k}(\theta)=\left\\{\begin{array}[]{lr}(1+\epsilon)(\frac{2^{k}}{\sin
2\alpha(\cos\theta+\sin\theta)}-P_{k}),~{}\mbox{if}~{}\gamma_{k-1}(\theta)\in(0,1),\\\
\infty,\ \ \mbox{otherwise},\end{array}\right.$

where $P_{k}=\prod_{j=1}^{k-1}(1+\sqrt{1-\gamma_{j}(\theta)^{2}})$.

Suppose there is a finite integer number $m\geq 2$, such that
$0<\gamma_{j}(\theta)<1$ for all $j=1,2,\cdots,m$. Then
$1<1+\sqrt{1-\gamma_{j}(\theta)^{2}}<2$. The bound
$1+\sqrt{1-\gamma_{j}^{2}(\theta)}<2$ implies that
$\gamma_{j}(\theta)/\gamma_{j-1}(\theta)>2$.

Take $\sin^{2}2\alpha=\frac{8}{9}$. We have

$\gamma_{1}(\theta)>\gamma_{1}=\frac{3\sqrt{2}}{4}\frac{2}{\cos\theta+\sin\theta}-1.$

Because $1\leq\cos\theta+\sin\theta\leq\sqrt{2}$,
$\gamma_{1}(\theta)\in[\frac{1}{2},1]$ for some $\theta\in[0,\frac{\pi}{4}]$.
In this case, $\gamma_{2}(\theta)>2\gamma_{1}(\theta)>1$. Thus, there does not
exist $\theta$ such that $\gamma_{2}(\theta)<1$. As $\gamma_{1}$ is a
monotonically decreasing function of $\sin^{2}2\alpha$, when
$\sin^{2}2\alpha>\frac{8}{9}$, $\gamma_{1}(\theta)$ could be less than
$\frac{1}{2}$. The smallest value of $\gamma_{1}(\theta)$ is $\sqrt{2}-1$. In
this case there exits some $\theta$ such that
$1>\gamma_{2}(\theta)>2\gamma_{1}(\theta)$, which completes the proof of
Theorem 4.


Appendix

This appendix describes the model of robot power in a typical circulation
loop. This description consists of three parts: 1) a model of flow in a
circulation loop and variation in hematocrit through that loop that treats
cells and plasma as separate compartments, 2) how oxygen transfer from blood
cells to plasma affects the concentration in these vessel compartments, and 3)
the rates at which blood cells release oxygen in response to consumption by
robots and tissue.

##  Appendix A Vessel Circuit Model

The circuit we consider starts as blood leaves a lung capillary, where its
oxygen concentration is set by that from the air in the lung. After leaving
the lung, oxygen concentration decreases due to consumption by robots and
tissue. This circuit continues through the body and back to the lung. The
circuit ends as the blood is about to enter a lung capillary. The time to
complete a circuit depends on where the blood goes, e.g., shorter times for
transit through the head than through the feet. We consider an average
circulation time of about one minute, i.e., the time required for the resting
heart rate, $5\,\mbox{L}/\mbox{min}$, to pump the entire blood volume,
$V_{\textnormal{blood}}$.

###  A.1 Hematocrit and Vessel Size

The concentration of glucose in blood is much larger than that of oxygen [29],
so oxygen is the chemical that limits the power available to robots. Blood
carries most of its oxygen bound to hemoglobin in red blood cells. The number
of red cells in blood is commonly expressed by the fraction of the blood
volume occupied by cells, i.e., the hematocrit, which is typically around
$45\%$ [29]. This is the average value over the entire blood volume.

In small vessels the fluid flow tends to push cells toward the center of the
vessel where they move faster than the average speed of the flow. This
Fahraeus effect decreases the hematocrit in small vessels. While there is
considerable variation among such vessels, on average the hematocrit $h$ in a
vessel of diameter $d$ (measured in microns) is approximately [66]

$\frac{h}{h_{\textnormal{full}}}=h_{\textnormal{full}}+(1-h_{\textnormal{full}})\left(1+1.7e^{-0.35d}-0.6e^{-0.01d}\right)$
(A.1)

where $h_{\textnormal{full}}$ is the hematocrit of the entire blood volume.
For example, when overall hematocrit is $h_{\textnormal{full}}=0.45$, Eq. A.1
gives $h=0.34$ in a capillary with diameter $d=8\,\mbox{$\mu$m}$. The
circulation model uses Eq. A.1 to determine hematocrit, thereby treating
hematocrit as only depending on vessel diameter.

###  A.2 Parts of a Circulation Loop

Oxygen concentration in the plasma arises from the balance of consumption, by
robots and tissue, and replenishment, from red blood cells. Replenishment
depends on the number of cells in the plasma, i.e., the hematocrit.

From Eq. A.1, hematocrit only deviates significantly from
$h_{\textnormal{full}}$ in vessels whose diameters are less than about a
millimeter. Thus it is not necessary to distinguish diameters larger than this
to determine hematocrit. Instead, we only need to explicitly account for
vessel diameter during the portion of the circuit through small vessels. This
observation allows grouping the transport through large-diameter vessels and
the heart, which all have the same hematocrit, to produce a circuit with the
following parts:

  1. 1.

a sequence of small veins of increasing diameters starting from the end of a
lung capillary

  2. 2.

large pulmonary veins, the heart, large arteries

  3. 3.

a sequence of small arteries of decreasing diameters to the start of a body
capillary

  4. 4.

a body capillary

  5. 5.

a sequence of small veins of increasing diameters from the end of a body
capillary

  6. 6.

large veins, the heart, large pulmonary arteries

  7. 7.

a sequence of small arteries of decreasing diameters to the start of a lung
capillary

Passage through the lung capillary at the end of the circuit contributes to
the total circulation time but is not explicitly included in the model.
Instead, the concentration at the start of the circuit, i.e., just after
passing through a lung capillary, is a boundary condition because the transit
time through lung capillaries is more than sufficient to saturate red cells
with oxygen [25]. Moreover, the large spacing between robots in capillaries,
even at the largest number of robots considered here (Table 1), means that
robots do not significantly alter the oxygen available in lung capillaries.

The change in oxygen concentration in each part of the circuit depends on its
transit time and hematocrit. Eq. A.1 gives hematocrit equal to
$h_{\textnormal{full}}$ in the large-vessel parts of the circuit, i.e.,
circuit parts 2 and 6. We set the transit time for these parts so the total
circuit time equals one minute when combined with the estimates of transit
times through small vessels discussed below.

For the capillary part of the circuit, we use typical capillary diameter and
transit time. Eq. A.1 gives the hematocrit in the body capillary based on its
typical diameter.

For the remaining parts of the circuit, consisting of small branching vessels,
we estimate hematocrit from the vessel diameter, $d$, via Eq. A.1, and transit
time from vessel geometry: diameter, $d$, length, $\ell$, and number of such
vessel segments, $N_{\textnormal{vessel}}$. There is considerable variation in
these values. For the purpose of this model, we use average values, analogous
to using the average relation between vessel diameter and hematocrit in Eq.
A.1.

We relate these geometric parameters to transit time in the vessels of a given
type (i.e., artery, capillary or vein) and diameter. In aggregate, small
vessels have larger cross section than large vessels, which leads to slower
speeds in those vessels [48]. Nevertheless, the short length of the small
vessels more than offsets their slower flow speed, so blood spends most of the
circuit time in large vessels. The aggregate cross section of vessels with
diameter $d$ is $\pi(d/2)^{2}N_{\textnormal{vessel}}$. The entire blood volume
passes through this aggregate cross section, when neglecting the relatively
small portion of the flow through portal systems, so that

$N_{\textnormal{vessel}}\;\frac{\pi}{4}d^{2}\;v=\frac{1}{T}V_{\textnormal{blood}}$
(A.2)

where $v$ is the average flow speed in the vessels and $T$ is average transit
time for the blood volume $V_{\textnormal{blood}}$, i.e, about one minute.
This relation gives $v$ in terms of $d$ and $N_{\textnormal{vessel}}$. This
velocity determines the segment’s transit time as $t=\ell/v$.

###  A.3 Geometry of Branching Vessels

As described above, we use geometric parameters of small branching vessels to
quantify those parts of the circuit. For evaluating robot power in an average
circuit, it is sufficient to use representative average branch geometry. One
approach to quantifying vessel geometry uses scaling laws that relate flow,
transit time, morphology (e.g., vessel diameter and length) and topology (i.e,
connectivity) [68]. A second approach uses empirical measurements of vessel
geometry. We apply this approach using the branching of vessels in the lung
for which data on complete vascular trees is available [84, 74, 42, 80].
Although the structure of capillary networks in the lung differs to some
extent from that in other parts of the body [80], we take the branching of
small arteries and veins in the lung as representative of such vessels for a
typical circuit through the body.

Specifically, we use the measured geometry of arterial and venous trees in the
lungs [42], which indicates that most of the flow is through vessels of
successive branch orders. Flow that skips a few branch orders corresponds to
blood that reaches capillaries through fewer branchings than blood that goes
through all orders. The model of average flow neglects this variability and
instead focuses on the main flow through successive branching orders. This
gives a sequence of vessel lengths, diameters and number of branches at each
branching order for both arterial and venous trees [42, Tables 2 and 5].
Doubling the number of branches to account for flow through both lungs, Eq.
A.2 determines transit speed, and hence transit time, for each level of
branching.

###  A.4 Transit Through a Circulation Loop

Large arteries branch into successively smaller vessels until they reach
capillaries, and then merge into increasingly large veins. For circulation
through a sequence of vessels $i=1,2,\ldots$, the total length is
$\sum_{i}\ell_{i}$, and total passage time is $\sum_{i}\ell_{i}/v_{i}$, which
equals the typical total transit time $T$.

As described above, for small vessels we use geometric parameters of branching
in lungs. Subtracting the transit time through small vessels from the total
time $T$ gives the time in the circuit parts corresponding to large vessels,
i.e., circuit parts 2 and 6 in Appendix A.2. To partition this time between
these two large-vessel parts, i.e., arteries and veins to and from a
capillary, respectively, we take the time in veins to be $1.5$ times larger
that in arteries, corresponding to somewhat slower flow speed in the veins.

Figure A.1: A circulation loop starting and ending in lung capillaries. The
shaded vertical bars indicate passage through the heart, and the vertical
lines near the center indicate passage through a capillary. The dashed
sections are interpolations for large vessels to and from the small-vessel
branches into and out of a body capillary. Red and blue portions of the curve
correspond to arteries and veins, respectively. Figure A.2: Diameters of
vessels, on a log scale, as a function of time through the vessel circuit,
highlighting the diameters of small vessels, where hematocrit deviates from
its overall value.

Combining these vessel properties, Fig. A.1 illustrates the circuit model. For
the sake of illustration, the diagram shows the transit through each side of
the heart as taking $1\,\mbox{s}$ to transverse $50\,\mbox{mm}$. Since
hematocrit is the same in the heart and in large vessels, the precise time
spent in the heart has no effect on the model results. The diagram
schematically illustrates transit through large vessels as an interpolation
(dashed curves) between the heart and the branching through small vessels to
and from a capillary. This interpolation matches typical average flow speed of
blood leaving and entering the heart. Specifically, flow speed in the aorta
averaged over a heart beat is around $110\,\mbox{mm}/\mbox{s}$ and speed in
the vena cava is around $135\,\mbox{mm}/\mbox{s}$ [31]. This interpolation
illustrates the flow through large vessels. However, the model of oxygen
consumption only depends on the transit time through those vessels, not the
specific shape of the dashed curves in the figure.

Fig. A.2 shows the vessel diameters for the vessel circuit. For small vessels,
the diameters and transit times are the values described in Appendix A.3. The
diameters for large vessels are not shown in the figure, and are not relevant
for this model: such vessels are large enough that hematocrit equals the
overall value $h_{\textnormal{full}}$ (see Eq. A.1).

###  A.5 Flow of Plasma and Cells

|
---|---
(a) | (b)
Figure A.3: (a) Schematic vessel branching. (b) Aggregated vessel cross
section, with larger cross section corresponding to smaller vessels. The
fraction of volume occupied by cells, i.e., hematocrit, is smaller in the
smaller vessels.

Fig. A.3a is a schematic of vessel branching. Smaller vessels have larger
total cross section and lower hematocrit than larger vessels. The vessel
properties relevant to robot oxygen consumption are the flow speed and
hematocrit, not the detailed vessel branching. This allows using the
simplified aggregated model of the flow illustrated in Fig. A.3b. In this
aggregated model, the circulation consists of a single vessel with varying
cross section and hematocrit.

To evaluate robot oxygen consumption over a minute or so, we consider a one-
dimensional model of vessel flow and average over the variation in speed due
to heart contractions. This gives flow speed $v(x)$ depending on the location
$x$ along the aggregated vessel, but not on time or how close the fluid is to
the vessel wall. With this time averaging, the total cross section $A(x)$ is
independent of time.

Consider the flow in a vessel with cross section area $A(x)$ and flow speed
$v(x)$ at position $x$. Fluid flows through the cross section at $x$ at a rate
$\rho v(x)A(x)$ where $\rho$ is the fluid density. This rate is constant
throughout the vessel for incompressible fluid, as given in Eq. A.2. Thus
speed is inversely proportional to cross section area:

$v(x)=v_{0}A_{0}/A(x)$ (A.3)

where $v_{0}$ and $A_{0}$ are the speed and cross section at an arbitrarily
specified position along the circuit.

The total cross section $A(x)$ varies with position, as does hematocrit:
larger cross sections correspond to the aggregation of smaller vessels, which
have lower hematocrit as described in Appendix A.1. Thus the fraction of the
total cross section corresponding to plasma and cells varies in the aggregate
vessel, as illustrated in Fig. A.3b. This observation leads to a model of
aggregated vessels with two compartments, plasma and cells, whose relative
cross sections change with the changing hematocrit in the vessels.

Evaluating how much cells replenish oxygen in blood plasma requires specifying
the average speed of cells and plasma, $v_{\textnormal{cell}}$ and
$v_{\textnormal{plasma}}$, respectively. Due to changing hematocrit, these
speeds are not the same and vary with the changing cross section of the
vessel. One relation among these speeds is the average speed $v$ in terms of
hematocrit $h$:

$v=(1-h)v_{\textnormal{plasma}}+hv_{\textnormal{cell}}$ (A.4)

In this expression, $v$ is the average speed in the vessel, taken from the
parameters of the circuit in Fig. A.1. This relation assumes the addition of
robots to the blood does not noticeably alter the speed, which is reasonable
with the small fraction of the volume occupied by robots (Table 1) [29]. This
expression splits the blood volume between plasma and cells, ignoring the tiny
volume occupied by robots.

Another relation among these speeds arises from conservation of flow. To see
this, consider a small volume of blood with hematocrit $h$ containing plasma
and cells in a vessel with cross section area $A$. In a small time $\Delta t$,
the volume of plasma and cells moving across that area are
$(1-h)Av_{\textnormal{plasma}}$ and $hAv_{\textnormal{cell}}$, respectively.
Over a circulation time, e.g., about a minute, neither plasma nor cell volumes
change significantly. Thus the flow rates of plasma and cells must be the same
throughout the circuit, i.e., equal to some constants
$\alpha_{\textnormal{plasma}}$ and $\alpha_{\textnormal{cell}}$. That is,
$(1-h)Av_{\textnormal{plasma}}=\alpha_{\textnormal{plasma}}$ and
$hAv_{\textnormal{cell}}=\alpha_{\textnormal{cell}}$. These relations imply
that the ratio of cell to plasma speed is independent of the total cross
section:

$\frac{v_{\textnormal{cell}}}{v_{\textnormal{plasma}}}=\frac{1-h}{h}\frac{\alpha_{\textnormal{cell}}}{\alpha_{\textnormal{plasma}}}$
(A.5)

The Fahraeus effect does not apply to large vessels, where cells and plasma
move together with the average flow of the blood, i.e,
$v_{\textnormal{cell}}=v_{\textnormal{plasma}}$ and hematocrit is
$h_{\textnormal{full}}$. For this case, Eq. A.5 gives
$\alpha_{\textnormal{cell}}/\alpha_{\textnormal{plasma}}=h_{\textnormal{full}}/(1-h_{\textnormal{full}})$.
Thus Eq. A.5 becomes

$\frac{v_{\textnormal{cell}}}{v_{\textnormal{plasma}}}=\frac{1-h}{h}{\Bigg{/}}\frac{1-h_{\textnormal{full}}}{h_{\textnormal{full}}}$
(A.6)

Combined with Eq. A.1, this gives the velocity ratio as a function of vessel
diameter. For example, with the range of hematocrits used here for large and
small vessels, Eq. A.6 gives a ratio around $1.7$ in small vessels, which
matches reported values [5].

The vessel diameters illustrated in Fig. A.2 and Eq. A.1 give the variation in
$h(x)$ in the aggregated vessels. This value, combined with Eq. A.4 and A.6
gives the speeds $v_{\textnormal{plasma}}(x)$ and $v_{\textnormal{cell}}(x)$
as a function of position in the circuit.

##  Appendix B Oxygen Concentration in Vessels

The aggregate vessel model shown in Fig. A.3b consists of two compartments,
plasma and cells. The cross sections of these compartments vary along the
length of the aggregate vessel. Robots are a small portion of the blood and
not treated as a separate compartment for the discussion of this section.
These compartments exchange oxygen with each other, and with robots and
tissue. This section describes how oxygen concentration changes in vessels
with variable cross section: first for a single vessel, then for two such
vessels treated as two compartments exchanging oxygen. These behaviors
determine how concentration changes in a volume of fluid moving with the flow
in one of these vessels. This discussion differs from behavior in vessels of a
fixed cross section due to the changing cross section and fraction of the
total cross section occupied by the two compartments.

###  B.1 Concentration in a Vessel with Changing Cross Section

A chemical in the fluid moves by convection with the fluid’s motion as well as
diffusion. For the scales and flow speeds considered here, diffusion is a
minor contribution to changing concentration. In particular, the Peclet number
characterizes the relative importance of convection and diffusion [76]:

$\mbox{Pe}=\frac{vd}{D_{\textnormal{$\mbox{O}_{2}$}}}$ (B.1)

where $v$ is the flow speed, $d$ a characteristic distance and
$D_{\textnormal{$\mbox{O}_{2}$}}$ the diffusion coefficient. For flow through
a vessel of diameter $d$, Pe roughly corresponds to the number of vessel
diameters required for diffusion to spread the chemical across the vessel. For
motion along the vessel, the distance at which $\mbox{Pe}\approx 1$, i.e.,
$d=D_{\textnormal{$\mbox{O}_{2}$}}/v$, is the distance at which diffusion and
convection have about the same effect on mass transport in a moving fluid. At
significantly longer distances, convection is the dominant effect.

The distance over which diffusion is important is largest in the vessels with
slowest flow, i.e., the capillaries. Capillary flow speeds are around
$1\,\mbox{mm}/\mbox{s}$ for which
$D_{\textnormal{$\mbox{O}_{2}$}}/v=2\,\mbox{$\mu$m}$ (see Table D.1). This
distance, comparable to the size of the robots, is considerably smaller than a
typical capillary length, i.e., a millimeter, and the length of the full
circulation loop (e.g., as indicated in Fig. A.1) relevant for evaluating
systemic effects of robot oxygen consumption. Moreover, the typical distance
between neighboring robots is larger than $D_{\textnormal{$\mbox{O}_{2}$}}/v$
in the scenarios considered here (see Table 1). Thus neighboring robots do not
directly compete with each other for oxygen. This is unlike the case of
closely spaced robots, such as aggregates on vessel walls, where robots
significantly reduce oxygen available to their neighbors [40]. In light of
these observations, we neglect diffusive transport in modeling the change in
oxygen on the scale of the circulation loop.

For convective flow, the chemical flux at position $x$ along the vessel is
$J(x)=v(x)c(x)$ where $c(x)$ is the chemical’s concentration. In a small
section of vessel between $x$ and $x+\Delta x$, in time $\Delta t$,
$J(x)A(x)\Delta t$ molecules enter that section of vessel, and $J(x+\Delta
x)A(x+\Delta x)\Delta t$ leave it. In addition, reactions, such as release of
oxygen by cells or consumption by tissue, change the concentration at a rate
$R$. Combining these contributions to concentration change gives

$\frac{\partial c}{\partial t}=-\frac{1}{A}\frac{\partial}{\partial x}(JA)+R$

From Eq. A.3, at position $x$, $JA=v_{0}A_{0}c$ so the rate of change of
concentration is

$\frac{\partial c}{\partial t}=-\frac{A_{0}}{A}v_{0}\frac{\partial c}{\partial
x}+R=-v\frac{\partial c}{\partial x}+R$ (B.2)

with $v$ given by Eq. A.3.

As described in Section C.2, the time constant for robot oxygen consumption is
less than a minute. Proposed applications for the robots involve operation
over at least tens of minutes, corresponding to many circulations. Thus a
reasonable simplification is to focus on the steady-state concentration in the
vessels, so Eq. B.2 becomes

$v\frac{\partial c}{\partial x}=R$ (B.3)

For transient operations, e.g., for a few seconds after robots start consuming
oxygen, robots will have more oxygen than indicated by the steady-state
analysis considered here.

###  B.2 Concentration in Flowing Compartments that Exchange Oxygen

As illustrated in Fig. A.3, the aggregated vessel consists of two main
compartments: plasma and cells. These compartments have separate flow speeds,
indicated by Eq. A.6. Each compartment acts as a separate vessel in terms of
flow, but they can exchange oxygen. Generalizing Eq. B.3 to account for this
exchange gives the behavior of the concentrations, $c_{1}$ and $c_{2}$, in the
two compartments:

$\begin{split}v_{1}\frac{\partial c_{1}}{\partial
x}&=R_{1}+R_{\textnormal{from 2}}\\\ v_{2}\frac{\partial c_{2}}{\partial
x}&=R_{2}-R_{\textnormal{to 1}}\end{split}$ (B.4)

where $R_{\textnormal{from 2}}$ is the change of concentration in compartment
1 due to chemicals from compartment 2, $R_{\textnormal{to 1}}$ is the decrease
of concentration in compartment 2 from chemicals that move to compartment 1,
and $R_{i}$ is the rate concentration changes in compartment $i$ due to
production of chemical in that compartment (with a negative value for chemical
consumed).

The two compartments share the total cross section of the aggregated vessel
model, $A(x)$. With hematocrit $h(x)$, the cross sections of the two
compartments are $A_{1}=(1-h(x))A(x)$ and $A_{2}=h(x)A(x)$. These expressions
assume robots occupy a negligible fraction of the vessel volume so there is no
need to subtract the fraction of volume occupied by robots. Conservation of
flow means $v_{1}(x)A_{1}(x)$ and $v_{2}(x)A_{2}(x)$ are independent of $x$.

$R_{\textnormal{from 2}}$ and $R_{\textnormal{to 1}}$ are rates of
concentration change in the two compartments from oxygen moving from
compartment 2 to compartment 1. These rates must account for different
volumes, i.e., a given amount of oxygen makes a larger change to concentration
in a smaller volume. The number of molecules in volume element $i$, extending
from $x$ to $x+\Delta x$, is $A_{i}(x)\Delta xc_{i}$. Thus the rate molecules
move from compartment 2 to compartment 1 is both $R_{\textnormal{to
1}}A_{2}\Delta x$ and $R_{\textnormal{from 2}}A_{1}\Delta x$. These must be
the same, so

$R_{\textnormal{to 1}}=\frac{1-h}{h}R_{\textnormal{from 2}}$ (B.5)

For example, when $h$ is small, the transfer of a given amount of oxygen has a
much larger effect on the concentration in compartment 2 than it does on that
of compartment 1.

###  B.3 Concentration in a Volume Moving with the Fluid

Consider a robot moving in a small volume of fluid in compartment 2 of the
two-compartment vessel model. Eq. B.4 describes the steady-state concentration
in the two compartments. These relations determine how the concentration in
that fluid volume changes as it moves through a circuit, such as illustrated
in Fig. A.1. In time $\Delta t$, the fluid volume moves from position $x$ to
$x+v_{2}\Delta t$. Thus the time rate of change of concentration in the volume
element of compartment $i$, due to motion with the speed $v_{2}$ is
$dc_{i}/dt=v_{2}\partial c_{i}/\partial x$. Thus Eq. B.4 gives

$\begin{split}\frac{dc_{1}}{dt}&=\left(R_{1}+R_{\textnormal{from
2}}\right)\frac{v_{2}}{v_{1}}\\\ \frac{dc_{2}}{dt}&=R_{2}-R_{\textnormal{to
1}}\end{split}$ (B.6)

for the rate that concentration changes in the two compartments from the
viewpoint of a volume moving with the flow in compartment 2.

The circulation shown in Fig. A.1 starts with blood leaving lung capillaries.
Thus the initial condition for Eq. B.6 is the concentration in the lung.

##  Appendix C Changing Oxygen Concentration

For this study, we assume robots move with blood cells rather than plasma.
Specifically, this assumes the flow pushes robots away from vessel walls in a
manner similar to the behavior of blood cells [85]. This is reasonable in
capillaries where cells move through single-file: robots are too large to fit
in the gap between cells and the vessel wall, so robots move between
successive blood cells, and hence at similar speed. In somewhat larger
vessels, if robots are pushed closer to the vessel wall than the cells, robots
would move somewhat more slowly. The difference in speeds between cells and
the average flow rate is largest when hematocrit differs most from its overall
value, i.e., in small vessels. Even then, the difference is relatively minor
due to the short time (a few seconds out of the one minute circulation) robots
spend in those vessels. Thus the model results are not very sensitive to the
accuracy of this assumption. With this assumption, robots travel with the
speed of compartment 2 discussed in Section B.3.

Completing the model of concentration change requires specifying the reaction
rates appearing in Eq. B.6. During the circulation, robots and tissue consume
oxygen from the plasma, and red cells replenish it. Fig. 3 illustrates the
processes that change oxygen concentration. In terms of Eq. B.6,
$R_{1}=R_{\textnormal{robot}}+R_{\textnormal{tissue}}$ is the rate, per unit
volume, that robots and tissue remove oxygen from the plasma. There is no
consumption within cells so $R_{2}=0$. The rates $R_{\textnormal{from 2}}$ and
$R_{\textnormal{to 1}}$ are the transfer rates from cells to plasma that
maintains the equilibrium between red cells and plasma given by Eq. C.1. The
remainder of this section quantifies these processes.

###  C.1 Blood Cells

Oxygen removed from blood plasma is partially replaced by oxygen released from
nearby red blood cells. The time scale for this process is less than
$100\,\mbox{ms}$ [16], which is much shorter than the one minute circulation
time considered here, and even the one second capillary transit time during
which tissue extracts oxygen from the blood. Thus a reasonable approximation
is that oxygen bound inside red cells is in equilibrium with the concentration
in the surrounding plasma.

Oxygen bound in red cells is characterized by the hemoglobin saturation $S$:
the fraction of hemoglobin capacity in a cell which has bound oxygen [55]. The
oxygen concentration in the cell is $C_{\textnormal{$\mbox{O}_{2}$}}^{\rm
max}S$, where $C_{\textnormal{$\mbox{O}_{2}$}}^{\rm max}$ is the concentration
in the cell when all the hemoglobin has bound oxygen.

Quantitatively, the equilibrium saturation, conventionally expressed in terms
of the equivalent partial pressure $p$ of $\mbox{O}_{2}$ in the fluid around
the cell, is described by the Hill equation [63, 34]:

$S_{\textnormal{equib}}(a)=\frac{a^{n}}{1+a^{n}}$ (C.1)

where $a=p/p_{50}$ is the partial pressure ratio, $p_{50}$ is the partial
pressure at which half the hemoglobin is bound to oxygen and $n$ characterizes
the steepness of the change from low to high saturation. The saturation ranges
from near $1$ in the lungs to around $1/3$ in tissues.

Henry’s Law relates the partial pressure to the oxygen concentration in the
plasma around the cell:
$p=H_{\textnormal{$\mbox{O}_{2}$}}C_{\textnormal{$\mbox{O}_{2}$}}$ with the
proportionality constant $H_{\textnormal{$\mbox{O}_{2}$}}$ depending on the
temperature. Thus,
$a=(H_{\textnormal{$\mbox{O}_{2}$}}/p_{50})C_{\textnormal{$\mbox{O}_{2}$}}=C_{\textnormal{$\mbox{O}_{2}$}}/C_{\textnormal{half}}$
where $C_{\textnormal{half}}=p_{50}/H_{\textnormal{$\mbox{O}_{2}$}}$, which is
about $2.2\times 10^{22}\,\mbox{molecule}/\mbox{m}^{3}$, using the values from
Table D.1. This is comparable to lower range of oxygen concentration in
tissue.

This model does not consider deviations from Eq. C.1, which mainly occur at
low saturations [34, 63], and variations in its parameters with changing blood
chemistry, such as pH and carbon dioxide concentration.

###  C.2 Robots

Consider a small volume $\Delta V$ of blood with oxygen concentration $c$ in
the plasma. This volume contains $\nu_{\textnormal{robot}}\Delta V$ robots,
where $\nu_{\textnormal{robot}}$ is the robot number density in the blood.
With each robot absorbing at the rate given by Eq. 1, the total rate robots
remove oxygen per unit volume of plasma is

$R_{\textnormal{robot}}=\gamma_{\textnormal{robot}}c$ (C.2)

with rate constant $\gamma_{\textnormal{robot}}=4\pi
D_{\textnormal{$\mbox{O}_{2}$}}r_{\textnormal{robot}}\nu_{\textnormal{robot}}/(1-h)$.

This absorption rate assumes each robot draws oxygen independently from a
fluid with concentration $c$, that is robots are sufficiently far apart that
they do not compete with their neighbors. Such competition would reduce the
local oxygen concentration and hence robot absorption rate [40]. This
competition is insignificant for robots separated by at least about ten times
their size [9], which is the case for the numbers of robots considered here
(see Table 1).

The time constant for robots removing oxygen from plasma is
$\tau_{\textnormal{robot}}=1/\gamma_{\textnormal{robot}}$. From Table D.1,
$\tau_{\textnormal{robot}}\approx 10\,\mbox{s}$ for $10^{10}$ robots.
$\tau_{\textnormal{robot}}$ is correspondingly smaller for the scenarios with
larger numbers of robots described in Table 1. Thus concentration changes due
to robots will reach steady-state within a single circulation time.

###  C.3 Tissue

Tissue extracts oxygen from blood passing through capillaries. In terms of the
aggregated vessel model, consumption by tissue occurs over a short distance
around the peak in total cross section shown in Fig. A.3b, which corresponds
to capillaries.

A simple model of tissue oxygen consumption is using a cylindrical region of
tissue around each capillary [46, 63]. The radius of this tissue cylinder,
$r_{\textnormal{tissue}}$, corresponds to a typical maximum distance of tissue
supplied from a capillary, which is a few cell diameters.

In this model, the volume of tissue receiving oxygen from a length $\Delta x$
of capillary, with radius $r_{\textnormal{cap}}$ is
$\pi(r_{\textnormal{tissue}}^{2}-r_{\textnormal{cap}}^{2})\Delta x$. So the
ratio of tissue to vessel volume is
$(r_{\textnormal{tissue}}/r_{\textnormal{cap}})^{2}-1$.

A consistency check on this model is that multiplying the total capillary
volume by the tissue to capillary volume ratio should equal the total body
volume. The total capillary volume is about $4\times 10^{5}\,\mbox{mm}^{3}$
[29]. From Table D.1, the ratio of tissue to capillary volume is about 100,
giving corresponding tissue volume of about $0.04\,\mbox{m}^{3}$, which is
comparable to body volume.

Models of oxygen use and power generation in tissues can account for tissue
structure [63]. A simpler approach [56], adopted in this paper, treats the
tissue surrounding the vessel as homogeneous and metabolizing oxygen at the
rate that produces power

$P_{\textnormal{tissue}}=P^{\textnormal{max}}_{\textnormal{tissue}}\frac{C_{\textnormal{$\mbox{O}_{2}$}}}{K_{\textnormal{tissue}}+C_{\textnormal{$\mbox{O}_{2}$}}}$
(C.3)

where $P^{\textnormal{max}}_{\textnormal{tissue}}$ is the nominal demanded
power density (i.e, power per unit volume) of the tissue and
$K_{\textnormal{tissue}}$ is the concentration of $\mbox{O}_{2}$ giving half
the maximum reaction rate. When oxygen concentration is substantially larger
than $K_{\textnormal{tissue}}$, tissue power is nearly independent of oxygen
concentration. In that case, tissue metabolic demand rather than available
oxygen limit tissue power.

Dividing $P_{\textnormal{tissue}}$ by the reaction energy per oxygen molecule
consumed gives the rate of oxygen consumption per unit volume of tissue.
Oxidizing a single glucose molecule consumes six oxygen molecules. Thus the
energy per oxygen molecule is $e/6$, with $e$ the energy from oxidizing one
glucose molecule, given in Table D.1.

For the resting tissue power demand considered here, oxygen concentration in
tissue is nearly constant as a function of distance from the capillary [46],
even with additional consumption from robots [40]. Thus, for evaluating tissue
oxygen consumption from capillaries, we consider the oxygen concentration in
the tissue to be the same as that of the plasma in the capillary. That is, in
Eq. C.3 we set $C_{\textnormal{$\mbox{O}_{2}$}}$ equal to the oxygen
concentration in the capillary surrounded by the tissue. In this case,
$P_{\textnormal{tissue}}$ is constant within the tissue cylinder and total
tissue consumption is $P_{\textnormal{tissue}}$ multiplied by the tissue
volume around the capillary. Tissue takes oxygen from plasma in the capillary,
so the rate it reduces concentration in the plasma is enhanced by the ratio of
tissue to capillary volume given above, and the fraction of the vessel volume
that is plasma, i.e., $1-h$. Combining these factors, the rate tissue reduces
oxygen concentration in the plasma is

$R_{\textnormal{tissue}}=\frac{P_{\textnormal{tissue}}}{e/6}\left(\left(\frac{r_{\textnormal{tissue}}}{r_{\textnormal{cap}}}\right)^{2}-1\right)\frac{1}{1-h}$
(C.4)

in the capillaries, and zero elsewhere.

##  Appendix D Model Parameters

geometry
---
capillary radius | $r_{\textnormal{cap}}$ | $4\,\mbox{$\mu$m}$
tissue cylinder radius | $r_{\textnormal{tissue}}$ | $40\,\mbox{$\mu$m}$
robot radius | $r_{\textnormal{robot}}$ | $1\,\mbox{$\mu$m}$
fluid
density | $\rho$ | $10^{3}\,\mbox{kg}/\mbox{m}^{3}$
hematocrit | $h_{\textnormal{full}}$ | $45\%$
blood volume | $V_{\textnormal{blood}}$ | $5.4\,\mbox{L}$
tissue
power demand | $P^{\textnormal{max}}_{\textnormal{tissue}}$ | $4\,\mbox{kW}/\mbox{m}^{3}$
$\mbox{O}_{2}$ concentration for half power | $K_{\textnormal{tissue}}$ | $10^{21}\,\mbox{molecule}/\mbox{m}^{3}$
power
reaction energy from one glucose molecule | $e$ | $4\times 10^{-18}\,\mbox{J}$
robot fuel cell efficiency | $f_{\textnormal{robot}}$ | $50\%$
red blood cells
partial pressure for 50% $\mbox{O}_{2}$ saturation | $p_{50}$ | $3500\,\mbox{Pa}$
$\mbox{O}_{2}$ saturation exponent | $n$ | $2.7$
maximum $\mbox{O}_{2}$ concentration in cell | $C_{\textnormal{$\mbox{O}_{2}$}}^{\rm max}$ | $10^{25}\,\mbox{molecule}/\mbox{m}^{3}$
oxygen
$\mbox{O}_{2}$ diffusion coefficient | $D_{\textnormal{$\mbox{O}_{2}$}}$ | $2\times 10^{-9}\,\mbox{m}^{2}/\mbox{s}$
$\mbox{O}_{2}$ concentration in lung | $C_{\textnormal{$\mbox{O}_{2}$}}^{\rm lung}$ | $7\times 10^{22}\,\mbox{molecule}/\mbox{m}^{3}$
$\mbox{O}_{2}$ partial pressure to concentration ratio | $H_{\textnormal{$\mbox{O}_{2}$}}$ | $1.6\times 10^{-19}\,\mbox{Pa}/(\mbox{molecule}/\mbox{m}^{3})$
Table D.1: Parameters for fluids, oxygen and microscopic robots considered
here. Sources for these values are described in the text.

Table D.1 gives the parameter values used to evaluate robot power in the
circulation.

For tissue, $K_{\textnormal{tissue}}$ is from Ref. [56], and the blood cell
parameters are from Refs. [16] and [63]. The oxygen concentration in the lung
corresponds to arterial concentration [29]. Concentrations of glucose in blood
plasma are in the millimolar range (about
$10^{24}\,\mbox{molecule}/\mbox{m}^{3}$), far larger than the oxygen
concentrations [29].

For evaluating the rate oxygen diffuses to the robot surface (e.g., in Eq. 1),
it is convenient to express oxygen concentration in terms of molecules per
unit volume, as shown in Fig. 6. By contrast, macroscopic studies usually
express concentrations in more readily measurable quantities. These include
moles of chemical per liter of fluid (i.e., molar, M) and grams of chemical
per cubic centimeter. Discussions of gases dissolved in blood often specify
concentration indirectly via the corresponding partial pressure of the gas
under standard conditions. As an example of these units, oxygen concentration
$C_{\textnormal{$\mbox{O}_{2}$}}=10^{22}\,\mbox{molecule}/\mbox{m}^{3}$
corresponds to a $17\,\mbox{$\mu$M}$ solution,
$0.53\,\mbox{$\mu$g}/\mbox{cm}^{3}$ and to a partial pressure of
$1600\,\mbox{Pa}$ or $12\,\mbox{mmHg}$. This concentration corresponds to
$0.037\mbox{cm}^{3}\mbox{O}_{2}/100\mbox{cm}^{3}\,\mbox{tissue}$ with oxygen
volume measured at standard temperature and pressure.

Tissue power demands vary considerably, depending on the tissue type and
overall activity level. We set $P^{\textnormal{max}}_{\textnormal{tissue}}$ to
a typical resting power demand [29]. For comparison, peak metabolic rate in
human tissue can be as high as $200\,\mbox{kW}/\mbox{m}^{3}$ [56].

Conventional fuel cell efficiencies are around $50\%$ [29]. While the
efficiency of the fuel cells required for micron-size robots remains to be
determined, for definiteness we use fuel cell efficiency
$f_{\textnormal{robot}}=50\%$.

## References

  * [1] Achraf Ben Amar, Ammar B. Kouki, and Hung Cao.  Power approaches for implantable medical devices.  Sensors, 15:28889–28914, 2015.
  * [2] Liang An et al.  Alkaline direct oxidation fuel cell with non-platinum catalysts capable of converting glucose to electricity at high power output.  J. of Power Sources, 196:186–190, 2011.
  * [3] Hellmut G. Augustin and Gou Young Koh.  Organotypic vasculature: From descriptive heterogeneity to functional pathophysiology.  Science, 357:eaal2379, 2017.
  * [4] Prosenjit Bagchi.  Mesoscale simulation of blood flow in small vessels.  Biophysical Journal, 92:1858–1877, 2007.
  * [5] A. C. L. Barnard, L. Lopez, and J. D. Hellums.  Basic theory of blood flow in capillaries.  Microvascular Research, 1:23–34, 1968.
  * [6] Robert H. Bartlett.  Alice in intensiveland: Being an essay on nonsense and common sense in the ICU, after the manner of Lewis Carroll.  Chest, 108:1129–1139, 1995.
  * [7] S. C. Barton, J. Gallaway, and P. Atanassov.  Enzymatic biofuel cells for implantable and microscale devices.  Chemical Reviews, 104:4867–4886, 2004.
  * [8] Kateryna Bazaka and Mohan V. Jacob.  Implantable devices: Issues and challenges.  Electronics, 2:1–34, 2013.
  * [9] Howard C. Berg.  Random Walks in Biology.  Princeton Univ. Press, 2nd edition, 1993.
  * [10] Eric Bonabeau, Marco Dorigo, and Guy Theraulaz.  Swarm Intelligence: From Natural to Artificial Systems.  Oxford University Press, Oxford, 1999.
  * [11] Allan M. Brooks and Michael S. Strano.  A conceptual advance that gives microrobots legs.  Nature, 584:530–531, 2020.
  * [12] Fei Cao et al.  Photoacoustic imaging in oxygen detection.  Applied Sciences, 7:1262, 2017.
  * [13] Hakan Ceylan et al.  3D printed personalized magnetic micromachines from patient blood–derived biomaterials.  Science Advances, 7:eabh0273, 2021.
  * [14] C. K. Chapler and S. M. Cain.  The physiologic reserve in oxygen carrying capacity: Studies in experimental hemodilution.  Canadian J. of Physiology and Pharmacology, 64:7–12, 1986.
  * [15] Swades K. Chaudhuri and Derek R. Lovley.  Electricity generation by direct oxidation of glucose in mediatorless microbial fuel cells.  Nature Biotechnology, 21:1229–1232, 2003.
  * [16] Alfred Clark, Jr. et al.  Oxygen delivery from red cells.  Biophysical Journal, 47:171–181, 1985.
  * [17] Alexis R. Cole et al.  Safety of prolonged inhalation of hydrogen gas in air in healthy adults.  Critical Care Explorations, 8:e543, 2021.
  * [18] Colin A. Cook et al.  Oxygen delivery from hyperbarically loaded microtanks extends cell viability in anoxic environments.  Biomaterials, 52:376–384, 2015.
  * [19] K. A. Cook-Chennault, N. Thambi, and A. M. Sastry.  Powering MEMS portable devices – a review of non-regenerative and regenerative power supply systems with special emphasis on piezoelectric energy harvesting systems.  Smart Materials and Structures, 17:043001, 2008.
  * [20] Frank Davis and Seamus P. J. Higson.  Biofuel cells: Recent advances and applications.  Biosensors and Bioelectronics, 22:1224–1235, 2007.
  * [21] Lixin Dong and Bradley J. Nelson.  Robotics in the small. part II: Nanorobotics.  IEEE Robotics & Automation Magazine, 14:111–121, 2007.
  * [22] Shawn M. Douglas, Ido Bachelet, and George M. Church.  A logic-gated nanorobot for targeted transport of molecular payloads.  Science, 335:831–834, 2012.
  * [23] J. Ernsting.  The effect of brief profound hypoxia upon the arterial and venous oxygen tensions in man.  J. of Physiology, 169:292–311, 1963.
  * [24] Xinjian Fan et al.  Scale-reconfigurable miniature ferrofluidic robots for negotiating sharply variable spaces.  Science Advances, 8:eabq1677, 2022.
  * [25] Joseph Feher.  Quantitative Human Physiology.  Academic Press, 2nd edition, 2017.
  * [26] Dario Floreano and Hod Lipson.  From individual robots to robot societies.  Science Robotics, 6:eabk2787, 2021.
  * [27] Darrel P. Francis et al.  Quantitative general theory for periodic breathing in chronic heart failure and its clinical implications.  Circulation, 102:2214–2221, 2000.
  * [28] Robert A. Freitas Jr.  Exploratory design in medical nanotechnology: A mechanical artificial red cell.  Artificial Cells, Blood Substitutes and Immobilization Biotechnology, 26:411–430, 1998.
  * [29] Robert A. Freitas Jr.  Nanomedicine, volume I: Basic Capabilities.  Landes Bioscience, Georgetown, TX, 1999.  Available at www.nanomedicine.com/NMI.htm.
  * [30] Robert A. Freitas Jr.  Nanomedicine, volume IIA: Biocompatibility.  Landes Bioscience, Georgetown, TX, 2003.  Available at www.nanomedicine.com/NMIIA.htm.
  * [31] Ivor T. Gabe et al.  Measurement of instantaneous blood flow velocity and pressure in conscious man with a catheter-tip velocity probe.  Circulation, 40:603–614, 1969.
  * [32] Mohamed A. Ghanem et al.  Noninvasive acoustic manipulation of objects in a living body.  Proc. of the National Academy of Sciences, 117:16848–16855, 2020\.
  * [33] Zuzand Gogova, Jiri Hanika, and Jozef Markos.  Optimal design of a multifunctional reactor for catalytic oxidation of glucose with fast catalyst deactivation.  In Alisson V. Brito, editor, Dynamic Modelling, chapter 12, pages 209–232. InTech, Rijeka, Croatia, 2010.
  * [34] Daniel Goldman.  Theoretical models of microvascular oxygen transport to tissue.  Microcirculation, 15:795–811, 2008.
  * [35] John Happel and Howard Brenner.  Low Reynolds Number Hydrodynamics.  Kluwer, The Hague, 2nd edition, 1983.
  * [36] Juan P. Hernandez-Ortiz, Christopher G. Stoltz, and Michael D. Graham.  Transport and collective dynamics in suspensions of confined swimming particles.  Physical Review Letters, 95:204501, 2005.
  * [37] Tad Hogg.  Coordinating microscopic robots in viscous fluids.  Autonomous Agents and Multi-Agent Systems, 14(3):271–305, 2007\.
  * [38] Tad Hogg.  Stress-based navigation for microscopic robots in viscous fluids.  J. of Micro-Bio Robotics, 15:59–67, 2018.
  * [39] Tad Hogg.  Acoustic power management by swarms of microscopic robots.  J. of Micro-Bio Robotics, 17:93–102, 2021.
  * [40] Tad Hogg and Robert A. Freitas Jr.  Chemical power for microscopic robots in capillaries.  Nanomedicine: Nanotechnology, Biology, and Medicine, 6:298–317, 2010.
  * [41] Tad Hogg and Robert A. Freitas Jr.  Acoustic communication for medical nanorobots.  Nano Communication Networks, 3:83–102, 2012.
  * [42] W. Huang et al.  Morphometry of the human pulmonary vasculature.  J. of Applied Physiology, 81:2123–2133, 1996.
  * [43] Edwin W. H. Jager, Olle Inganas, and Ingemar Lundstrom.  Microrobots for micrometer-size objects in aqueous media: Potential tools for single-cell manipulation.  Science, 288:2335–2338, 2000.
  * [44] James Keener and James Sneyd.  Mathematical Physiology, volume II: Systems Physiology.  Springer, 2nd edition, 2009.
  * [45] Muhammad Saad Khan et al.  Oxygen-carrying micro/nanobubbles: Composition, synthesis techniques and potential prospects in photo-triggered theranostics.  Molecules, 23:2210, 2018.
  * [46] August Krogh.  The number and distribution of capillaries in muscles with calculations of the oxygen pressure head necessary for supplying the tissue.  J. of Physiology, 52:409–415, 1919.
  * [47] Philip W. Kuchel et al.  Surface model of the human red blood cell simulating changes in membrane curvature under strain.  Scientific Reports, 11:13712, 2021.
  * [48] Michael LaBarbera.  Principles of design of fluid transport systems in zoology.  Science, 249:992–1000, 1990.
  * [49] Marc J. Lajoie et al.  Designed protein logic to target cells with precise combinations of surface antigens.  Science, 369:1637–1643, 2020.
  * [50] Rienhard Lanbenbacher, James P. Sluka, and James A. Glazier.  Using digital twins in viral infection.  Science, 371:1105–1106, 2021.
  * [51] Jinxing Li, Berta Esteban-Fernandez de Avila, Wei Gao, Liangfang Zhang, and Joseph Wang.  Micro/nanorobots for biomedicine: Delivery, surgery, sensing, and detoxification.  Science Robotics, 2:eaam6431, 2017.
  * [52] C. S. Lim et al.  Venous hypoxia: A poorly studied etiological factor of varicose veins.  J. of Vascular Research, 48:185–194, 2011.
  * [53] Wendell A. Lim.  The emerging era of cell engineering: Harnessing the modularity of cells to program complex biological function.  Science, 378:848–852, 2022.
  * [54] Sylvain Martel.  The coming invasion of the medical nanorobots.  Nanotechnology Perceptions, 3:165–173, 2007.
  * [55] Benjamin Mauroy.  Following red blood cells in a pulmonary capillary.  ESAIM Proceedings, 23:48–65, 2008.
  * [56] B. J. McGuire and T. W. Secomb.  A theoretical model for oxygen transport in skeletal muscle under conditions of high oxygen demand.  Journal of Applied Physiology, 91:2255–2265, 2001.
  * [57] George McHedlishvili and Manana Varazashvili.  Hematocrit in cerebral capillaries and veins under control and ischemic conditions.  J. of Cerebral Blood Flow & Metabolism, 7:739–744, 1987.
  * [58] Sumit Mohanty, Islam S. M. Khalil, and Sarthak Misra.  Contactless acoustic micro/nano manipulation: a paradigm for next generation applications in life sciences.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 476:20200621, 2020.
  * [59] Kelly Morris.  Macrodoctor, come meet the nanodoctors.  The Lancet, 357:778, March 10 2001.
  * [60] Leslie E. Morris and Herrman L. Blumgart.  Velocity of blood flow in health and disease.  Circulation, XV:448–460, 1957.
  * [61] Bradley J. Nelson, Ioannis K. Kaliakatsos, and Jake J. Abbott.  Microrobots for minimally invasive medicine.  Annual Review of Biomedical Engineering, 12:55–85, 2010.
  * [62] Meir Nitzan et al.  Pulse oximetry: fundamentals and technology update.  Medical Devices: Evidence and Research, 7:231–239, 2014.
  * [63] Aleksander S. Popel.  Theory of oxygen transport to tissue.  Critical Reviews in Biomedical Engineering, 17:257–321, 1989.
  * [64] Aleksander S. Popel and Paul C. Johnson.  Microcirculation and hemorheology.  Annual Review of Fluid Mechanics, 37:43–69, 2005.
  * [65] Manu Prakash and Neil Gershenfeld.  Microfluidic bubble logic.  Science, 315:832–835, 2007.
  * [66] A. R. Pries, D. Neuhaus, and P. Gaehtgens.  Blood viscosity in tube flow: dependence on diameter and hematocrit.  American J of Physiology, 263:H1770–H1778, 1992.
  * [67] Benjamin I. Rapoport, Jakub T. Kedzierski, and Rahul Sarpeshkar.  A glucose fuel cell for implantable brain-machine interfaces.  PLoS ONE, 7(6):e38436, 2012.
  * [68] Mohammad S. Razavi, Ebrahim Shirani, and Ghassan S. Kassab.  Scaling laws of flow rate, vessel blood volume, lengths, and transit times with number of capillaries.  Frontiers in Physiology, 9:581, 2018.
  * [69] Erik L. Ritman and Amir Lerman.  The dynamic vasa vasorum.  Cadiovascular Research, 75:649–658, 2007.
  * [70] Michael Rubenstein, Alejandro Cornejo, and Radhika Nagpal.  Programmable self-assembly in a thousand-robot swarm.  Science, 345:795–799, 2014.
  * [71] S. K. Sahoo, S. Parveen, and J. J. Panda.  The present and future of nanotechnology in human health care.  Nanomedicine: Nanotechnology, Biology and Medicine, 3:20–31, 2007\.
  * [72] Simone Schuerle et al.  Synthetic and living micropropellers for convection-enhanced nanoparticle transport.  Science Advances, 5:eaav4803, 2019.
  * [73] Mark J. Schulz, Vesselin N. Shanov, and Yeoheung Yun, editors.  Nanomedicine Design of Particles, Sensors, Motors, Implants, Robots, and Devices.  Engineering in Medicine and Biology. Artech House, Boston, 2009.
  * [74] Siam Singhal et al.  Morphometry of the human pulmonary arterial tree.  Circulation Research, 33:190–197, 1973.
  * [75] Ambika Somasundar and Ayusman Sen.  Chemically propelled nano and micromotors in the body: Quo vadis?  Small, 17:2007102, 2021.
  * [76] Todd M. Squires and Stephen R. Quake.  Microfluidics: Fluid physics at the nanoliter scale.  Reviews of Modern Physics, 77:977–1026, 2005.
  * [77] Hagit Stauber et al.  Red blood cell dynamics in biomimetic microfluidic networks of pulmonary alveolar capillaries.  Biomicrofluidics, 11:014103, 2017.
  * [78] Soren Strauss et al.  Phyllotaxis: is the golden angle optimal for light capture?  New Phytologist, 225:499–510, 2019.
  * [79] C. A. Taylor and C. A. Figueroa.  Patient-specific modeling of cardiovascular mechanics.  Annual Review of Biomedical Engineering, 11:109–134, 2009.
  * [80] Mary I. Townsley.  Structure and composition of pulmonary arteries, capillaries and veins.  Comprehensive Physiology, 2:675–709, 2012.
  * [81] Tamas Vicsek et al.  Novel type of phase transition in a system of self-driven particles.  Physical Review Letters, 75:1226–1229, 1995.
  * [82] Xudong Wang et al.  Direct-current nanogenerator driven by ultrasonic waves.  Science, 316:102–105, 2007.
  * [83] Barbara Webb.  Robots with insect brains.  Science, 368:244–245, 2020.
  * [84] Ewald R. Weibel and Domingo M. Gomez.  Architecture of the human lung.  Science, 137:577–585, 1962.
  * [85] R. L. Whitmore.  A theory of blood flow in small vessels.  J. of Applied Physiology, 22:767–771, 1967.
  * [86] Inbal Wiesel-Kapah, Gal A. Kaminka, Guy Hachmon, Noa Agmon, and Ido Bachelet.  Rule-based programming of molecular robot swarms for biomedical applications.  In G. Brewka, editor, Proc. of the 25th Intl. Joint Conf. on Artificial Intelligence (IJCAI16), pages 3505–3512. AAAI Press, 2016.
  * [87] David F. Wilson and Franz M. Matschinsky.  Hyperbaric oxygen toxicity in brain: A case of hyperoxia induced hypoglycemic brain syndrome.  Medical Hypotheses, 132:109375, 2019.
  * [88] Hui Xie et al.  Reconfigurable magnetic microrobot swarm: Multimode transformation, locomotion, and manipulation.  Science Robotics, 4:eaav8006, 2019.
  * [89] A. Zebda et al.  Single glucose biofuel cells implanted in rats power electronic devices.  Scientific Reports, 3:1516, 2013.