4
Theory of Macro-Scale Molecular Communications

4.1 Introduction

T his chapter focuses on the mathematical modelling of the macro-scale molecular communication. As mentioned, macro-scale is seldom studied compared to micro-scale therefore mathematical models describing the macro-scale effects (i.e., advection) are rare. In this chapter, mathematical descriptions used throughout the study are given. Depending on the environment (i.e., boundary-less, open-air) models are revised. Each Cartesian dimension ( $x, y, z$ ) is considered in modelling the communication using the advective-diffusion equation (ADE).

4.2 Transmission of Molecules

A communication based on the transmission of particles (i.e., molecules, silt, pollutants etc.) can be explained using the mass transport phenomena. This can be briefly explained as the exchange of mass, energy, charge and momentum between observed systems. Since molecular communications utilises the mass of the particles in the exchange between systems, Fick’s $2^{nd}$ law can be used to describe the system when there is no additional flow present in the system [133]:

\frac{\partial c}{\partial t} = D \nabla^{2} c,

(4.1)

The derivation of the Fick’s $2^{n d}$ law can be seen in Appendix 8.2.5 . At micro scales (nm to $μ$ m) the above equation can be used to model and explain the whole system. However, at macro-scales (cm to m) the system needs external propulsion to transmit the chemicals. Therefore, an additional parameter, the advection flux, is needed to introduce the flow component into the system. To describe the system, a generalised expression of advection-diffusion (general scalar transport [343]) is used:

\frac{\partial c}{\partial t} = D \nabla^{2} c - \nabla \cdot (u c) + R .

(4.2)

The derivation of the ADE from the continuity equation is shown in Chapter 2. To simplify the equation, it is assumed that the diffusion coefficient is constant ( $\partial D ∕ \partial t = 0$ ), there are no sources or sinks ( $R$ = 0), and the velocity has zero divergence ( $u \equiv 0$ ). Based on these definitions, the equation can be rewritten as:

\frac{\partial c}{\partial t} = D \nabla^{2} c - \nabla \cdot (u c) .

(4.3)

The solution for this equation can be derived based on the following initial boundary conditions, known as “ thin-film solution” :

\begin{align} c (| x | > 0, | y | > 0, | z | > 0, t_{0}) & = 0, & (4.4a) \\ c (x = 0, y = 0, z = 0, t_{0}) & = M_{0} δ (x) δ (y) δ (z), & (4.4b) \\ c (| x | \to \infty, | y | \to \infty, | z | \to \infty, t) & = 0 . & (4.4c) \end{align}

where $M_{0}$ is the initial mass injected into the environment (kg), $x$ , $y$ , $z$ are the dimensions in which the propagation takes place, $t_{0}$ is the initial time (s) and $δ (x)$ , $δ (y)$ , $δ (z)$ are the Dirac delta functions of their respective dimensions. The mathematical definition of Dirac delta function is given below [344]:

δ (x) = {\begin{matrix} 1, & x = 0, \\ 0, & x \neq 0 . \end{matrix}

(4.5)

Based on these boundaries to the PDE , the solution with no physical boundaries applied to the propagation can be expressed as:

\begin{align} c (λ_{x_{0}}, t) & = \frac{M_{0}}{A_{y z} \sqrt{4 π D_{x} t}} \exp (- \frac{λ_{x_{0}}^{2}}{4 D_{x} t}), & (4.6a) \\ c (λ_{x_{0}}, λ_{y_{0}}, t) & = \frac{M_{0}}{L_{z} 4 π t \sqrt{D_{x} D_{y}}} \exp (- \frac{λ_{x_{0}}^{2}}{4 D_{x} t} - \frac{λ_{y_{0}}^{2}}{4 D_{y} t}), & (4.6b) \\ c (λ_{x_{0}}, λ_{y_{0}}, λ_{z_{0}}, t) & = \frac{M_{0}}{{(4 π t)}^{3 ∕ 2} \sqrt{D_{x} D_{y} D_{z}}} \exp (- \frac{λ_{x_{0}}^{2}}{4 D_{x} t} - \frac{λ_{y_{0}}^{2}}{4 D_{y} t} - \frac{λ_{z_{0}}^{2}}{4 D_{z} t}) . & (4.6c) \end{align}

where $D_{x}$ , $D_{y}$ and $D_{z}$ are the diffusive coefficients in $x$ , $y$ and $z$ coordinates respectively ( ${cm}^{2} ∕ s$ ), $A_{y z}$ and $L_{z}$ are the area-scale ( $m^{2}$ ) and length-scale (m) of the neglected dimensions. $t$ is the duration of the experiment (s), $c$ is the concentration of the chemical in 1D, 2D and 3D ( $kg ∕ m$ , $kg ∕ m^{2}$ and $kg ∕ m^{3}$ ) and $λ_{x_{0}}$ , $λ_{y_{0}}$ and $λ_{z_{0}}$ are the moving reference frames with the following descriptions:

\begin{align} λ_{x_{0}} & = x - (x_{0} + u_{x} t), & (4.7a) \\ λ_{y_{0}} & = y - (y_{0} + u_{y} t), & (4.7b) \\ λ_{z_{0}} & = z - (z_{0} + u_{z} t), & (4.7c) \end{align}

where $x_{0}$ , $y_{0}$ and $z_{0}$ are the injection points of the mass (m), $u_{x}$ , $u_{y}$ and $u_{z}$ are the mean flow velocity (m/s) and $u_{x} t$ , $u_{y} t$ and $u_{z} t$ are the distances the centre of mass of the cloud travelled (m) in a given time of $t$ (s).

In macro-scale molecular communications, the diffusion ( $D_{x}$ , $D_{y}$ and $D_{z}$ ) caused by the particles may be governed by a combination of momentum and turbulent diffusion.

The equations given in Eq. ( 4.6 ) quantify the concentration value of the sample in a given time ( $t$ ) and space $(x, y, z)$ . By integrating the concentration function with respect to distance ( $x, y, z$ ) the particles that are present in the environment ( $𝜃_{E}$ ) in a given time $t$ can be calculated:

\begin{align} 𝜃_{E} (x_{d}, t) & = \int_{- x_{𝜖}}^{+ x_{d}} c (x, t) d x, & (4.8a) \\ 𝜃_{E} (x_{d}, y_{d}, t) & = \int_{- x_{𝜖}}^{+ x_{d}} \int_{- y_{𝜖}}^{+ y_{d}} c (x, y, t) d x d y, & (4.8b) \\ 𝜃_{E} (x_{d}, y_{d}, z_{d}, t) & = \int_{- x_{𝜖}}^{+ x_{d}} \int_{- y_{𝜖}}^{+ y_{d}} \int_{- z_{𝜖}}^{+ z_{d}} c (x, y, z, t) d x d y d z, & (4.8c) \end{align}

where $x_{d}, y_{d}, z_{d}$ are the distances from the detector to the origin point (m) and $x_{𝜖}, y_{𝜖}, z_{𝜖}$ are the distance that particles travel against the flow (m) in Cartesian coordinates ( $x, y, z$ ). As the system has no sink nor source ( $R = 0$ ), the chemicals that are used in the transmission can either be in transmission ( $𝜃_{E}$ ) or have been absorbed by the detector ( $𝜃_{1}$ ). Therefore, both the aforementioned mass values must add up to the initial introduction of mass ( $M_{0}$ ) at the beginning of the transmission:

M_{0} = 𝜃_{E} + 𝜃_{1} .

(4.9)

The mass absorbed by the detector ( $𝜃_{1}$ ), therefore, can be calculated by simply subtracted from the initial mass ( $M_{0}$ ) [2], [3], [142]:

\begin{align} 𝜃_{1} (x_{d}, t) & = M_{0} - \int_{- x_{𝜖}}^{+ x_{d}} c (x, t) d x, & (4.10a) \\ 𝜃_{1} (x_{d}, y_{d}, t) & = M_{0} - \int_{- x_{𝜖}}^{+ x_{d}} \int_{- y_{𝜖}}^{+ y_{d}} c (x, y, t) d x d y, & (4.10b) \\ 𝜃_{1} (x_{d}, y_{d}, z_{d}, t) & = M_{0} - \int_{- x_{𝜖}}^{+ x_{d}} \int_{- y_{𝜖}}^{+ y_{d}} \int_{- z_{𝜖}}^{+ z_{d}} c (x, y, z, t) d x d y d z . & (4.10c) \end{align}

The closed-form solution for the mass absorption in 1D, 2D and 3D can be seen in Eq. ( 4.11 ) :

\begin{align} \begin{aligned} 𝜃_{1} (x_{d}, t) & = M_{0} - \frac{M_{0}}{2} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{x} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{x} t}})], \end{aligned} & (4.11a) \\ \begin{aligned} 𝜃_{1} (x_{d}, y_{d}, t) & = M_{0} - \frac{M_{0}}{4} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{x} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{x} t}})] \\ \times [e r f (\frac{y_{d} - u_{y} t}{\sqrt{4 D_{y} t}}) + e r f (\frac{y_{𝜖} + u_{y} t}{\sqrt{4 D_{y} t}})], \end{aligned} & (4.11b) \\ \begin{aligned} 𝜃_{1} (x_{d}, y_{d}, z_{d}, t) & = M_{0} - \frac{M_{0}}{8} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{x} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{x} t}})] \\ \times [e r f (\frac{y_{d} - u_{y} t}{\sqrt{4 D_{y} t}}) & + e r f (\frac{y_{𝜖} + u_{y} t}{\sqrt{4 D_{y} t}})] [e r f (\frac{z_{d} - u_{z} t}{\sqrt{4 D_{z} t}}) + e r f (\frac{z_{𝜖} + u_{z} t}{\sqrt{4 D_{z} t}})], \end{aligned} & (4.11c) \end{align}

where $e r f (\cdot)$ is the error function with the following definition [345]:

e r f (x) = \frac{1}{\sqrt{π}} \int_{- x}^{+ x} \exp (- t^{2}) d t .

(4.12)

The chemicals that are absorbed by the detector ( $𝜃_{1}$ ) in a given period of $T_{S}$ are given in Eq. ( 4.13 ):

\begin{align} M_{R} = 𝜃_{1} (x_{d}, t = T_{S}) & - 𝜃_{1} (x_{d}, t = 0), & (4.13a) \\ M_{R} = 𝜃_{1} (x_{d}, y_{d}, t = T_{S}) & - 𝜃_{1} (x_{d}, y_{d}, t = 0), & (4.13b) \\ M_{R} = 𝜃_{1} (x_{d}, y_{d}, z_{d}, t = T_{S}) & - 𝜃_{1} (x_{d}, y_{d}, z_{d}, t = 0) . & (4.13c) \end{align}

Therefore, the removal of chemicals from the detector ( $𝜃_{0}$ ) to the outside environment can be expressed by the following expression:

\begin{align} \begin{aligned} 𝜃_{0} (x_{d}, t) & = \frac{M_{R}}{2} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{x} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{x} t}})] \end{aligned}, & (4.14a) \\ \begin{aligned} 𝜃_{0} (x_{d}, y_{d}, t) & = \frac{M_{R}}{4} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{x} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{x} t}})] \\ \times [e r f (\frac{y_{d} - u_{y} t}{\sqrt{4 D_{y} t}}) + e r f (\frac{y_{𝜖} + u_{y} t}{\sqrt{4 D_{y} t}})], \end{aligned} & (4.14b) \\ \begin{aligned} 𝜃_{0} (x_{d}, y_{d}, z_{d}, t) & = \frac{M_{R}}{8} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{x} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{x} t}})] \\ \times [e r f (\frac{y_{d} - u_{y} t}{\sqrt{4 D_{y} t}}) & + e r f (\frac{y_{𝜖} + u_{y} t}{\sqrt{4 D_{y} t}})] [e r f (\frac{z_{d} - u_{z} t}{\sqrt{4 D_{z} t}}) + e r f (\frac{z_{𝜖} + u_{z} t}{\sqrt{4 D_{z} t}})] . \end{aligned} & (4.14c) \end{align}

As can be seen from Eq. ( 4.11 ) and Eq. ( 4.14 ) the mass parameter is different in each equation: former being the mass injected into the environment ( $M_{0}$ ) and the latter is the mass that is absorbed by the detector ( $M_{R}$ ). While the injected mass is a predefined value, $M_{R}$ dynamically changes as the transmission evolves [2], [3].

The model that represents the transmission can be seen in Figure 4.1 .

Figure 4.1: A diagram of the model used in the study. (A) At

t

= 0 s mass is injected into the environment. This is defined by the initial boundary condition

c (x_{0}, t_{0}) = M_{0} δ (x)

. (B) As the transmission evolves, the particles start propagating via advection and diffusion. (C) When the specified period is passed the transmitter stops releasing particles and transmits only advective flow. (D) The flow then forces the particles to be transferred from the detector to the outside environment. (E) After a finite amount has passed the particles are removed from the detector and transferred to the outside environment. However, some particles are left in the detector, which can cause inter-symbol interference for future transmissions (

𝜃_{I S I}

4.2.1 Proposed Transmission Model

Based on the interaction of chemicals with the environment and the detector mentioned in Section 4.2 , experimental work was carried out and was validated by the use of the aforementioned model in 3D. In the experimental setup, the flow occurs in the $x$ -direction only. Therefore ( $u_{y}$ ) and ( $u_{z}$ ) are assumed to be negligible :

u_{y} = 0, u_{z} = 0 .

(4.15)

The $y$ and $z$ coordinates of the environment are given the value of the radius of the detector ( $R_{D}$ ):

y_{d} = y_{𝜖} = R_{D}, z_{d} = z_{𝜖} = R_{D} .

(4.16)

In the propagation there are two types of diffusion, former being the longitudinal diffusion ( $D_{L}$ ) and the latter being the radial diffusion ( $D_{T}$ ). Because the advective flow ( $u_{x}$ ) occurs in the $x$ -direction, the particles propagate in the $x$ -direction (i.e., increases the diffusion in the $x$ -direction). A final note is that the diffusion occurring along the $x$ -axis is the longitudinal diffusion and the diffusion happening in the y and z axes are defined as radial diffusion ( $D_{T}$ ):

D_{L} = D_{x}, D_{T} = D_{y} = D_{z}, D_{L} ≫ D_{T} .

(4.17)

Based on these definitions the equation that gives the absorbed mass by the detector from the outside environment ( $𝜃_{1}$ ) and the mass removed from the detector to the outside environment ( $𝜃_{0}$ ) given can be shown as:

\begin{align} 𝜃_{1} (x_{d}, R_{D}, t) & = M_{0} - \frac{M_{0}}{2} e r f {(\frac{R_{D}}{\sqrt{4 D_{T} t}})}^{2} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{L} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{L} t}})], & (4.18a) \\ 𝜃_{0} (x_{d}, R_{D}, t) & = \frac{M_{R}}{2} e r f {(\frac{R_{D}}{\sqrt{4 D_{T} t}})}^{2} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{L} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{L} t}})] & (4.18b) \end{align}

4.2.2 Open-Air Transmission

To model the effect of open air transmission an additional parameter, decay ( $λ_{D}$ ), is required. Therefore, to include this parameter in the model, the ADE is revisited.

\frac{\partial c}{\partial t} = \nabla \cdot (D \nabla c) - \nabla \cdot (u c) + R .

(4.19)

The system has no defined boundary (i.e., pipe), therefore the amount introduced by the system will not be equal to the amount detected by the receiver since some amount of particles will stray away from the path of the detector. This property can be modelled by introducing a sink ( $R = - λ_{D} c$ ) to the equation. Inserting this parameter to Eq. ( 4.19 ) yields:

Figure 4.2: A Diagram showing the effect of open-air transmission.

\frac{\partial c}{\partial t} = \nabla \cdot (D \nabla c) - \nabla \cdot (u c) - λ_{D} c,

(4.20)

where $λ_{D}$ is the decay parameter of the function ( $s^{- 1})$ . The prototypical solution ( $c$ ( $x$ = 0, $t$ = 0) = $M_{0} δ$ ( $x$ )) to the given expression in Eq. ( 4.20 ) is given in Eq. ( 21 ):

\begin{matrix} c (x, y, z, t) = \frac{M_{0}}{{(4 π t)}^{3 ∕ 2} \sqrt{D_{x} D_{y} D_{z}}} \\ \times \exp (- \frac{{(x - u_{x} t)}^{2}}{4 D_{x} t} - \frac{{(y - u_{y} t)}^{2}}{4 D_{y} t} - \frac{{(z - u_{z} t)}^{2}}{4 D_{z} t} - λ_{D} t) . & (4.21) \end{matrix}

The equation above represents the concentration of the introduced sample in a given time and space. By integrating the concentration function with respect to distance the particles that are present in the environment can be calculated. To calculate the chemicals absorbed by the detector, the integration function is subtracted from the injected mass [142]:

𝜃_{1} (x_{d}, y_{d}, z_{d}, t) = M_{0} - \int_{- x_{𝜖}}^{+ x_{d}} \int_{- y_{𝜖}}^{+ y_{d}} \int_{- z_{𝜖}}^{+ z_{d}} c (x, y, z, t) d x d y d z .

(4.22)

The solution to the above Eq. ( 4.22 ) for open distance transmission with decay can be obtained by integrating the concentration function in 3D Cartesian dimensions and subtracting the original mass. This solution is given in Eq. ( 23 ):

\begin{matrix} 𝜃_{1} (x_{d}, R_{D}, t) = M_{0} \exp (- λ_{D} t) \\ \times {1 - \frac{1}{2} e r f {(\frac{R_{D}}{\sqrt{4 D_{T} t}})}^{2} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{L} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{L} t}})]} . & (4.23) \end{matrix}

The chemicals that are absorbed by the detector from the outside environment ( $𝜃_{1}$ ) in a given period of $T_{S}$ is given in Eq. ( 4.24 ):

M_{R} = 𝜃_{1} (x_{d}, R_{D}, t = T_{S}) - 𝜃_{1} (x_{d}, R_{D}, t = 0) .

(4.24)

Therefore, the removal of chemicals from the detector ( $𝜃_{0}$ ) to the outside environment can be expressed by the following Eq. ( 25 ):

\begin{matrix} 𝜃_{0} (x_{d}, r, t) = M_{R} \exp (- λ_{D} t) \\ \times {\frac{1}{2} e r f {(\frac{r}{\sqrt{4 D_{T} t}})}^{2} [e r f (\frac{x_{d} - u_{x} t}{\sqrt{4 D_{L} t}}) + e r f (\frac{x_{𝜖} + u_{x} t}{\sqrt{4 D_{L} t}})]} . & (4.25) \end{matrix}

As can be seen from Eq. ( 23 ) and Eq. ( 25 ) the mass parameter is different in each equation: former being the mass injected into the environment ( $M_{0}$ ) and the latter is the mass that is absorbed by the detector ( $M_{R}$ ). This process of introduction/removal of particles can be seen in detail in [2]. To model the detrimental effects of open air transmission, the decay parameter ( $λ_{D}$ ) with respect to transmission distance ( $x_{d}$ ) is approximated using a power equation with $a$ and $b$ being fitting parameters [346], [347]. The parameters of the decay equation ( $a$ , $b$ ) can be influenced by numerous parameters, such as the temperature and the pressure of the environment, particles and eddies present in the transmission medium.

λ_{D} (x_{d}) = a x_{d}^{b}, a = 6.743 \times 1 0^{- 5}, b = 2.616 .

(4.26)

While this decay function along with the equations developed are used in open-distance transmission, the problems faced with closed source are different a different approach needs to be implemented.

4.2.3 Closed-Boundary

To describe the closed boundary effect on the communication, the 3D solution to the PDE given in Eq. ( 4.6c ) is revisited:

\begin{matrix} c (x, y, z, t) = \frac{M_{0}}{\sqrt{{(4 π D_{x} D_{y} D_{z} t)}^{3}}} \\ \times \exp (- \frac{{(x - u_{x} t)}^{2}}{4 D_{x} t} - \frac{{(y - u_{y} t)}^{2}}{4 D_{y} t} - \frac{{(z - u_{z} t)}^{2}}{4 D_{z} t}), & (4.27) \end{matrix}

where $D_{x}$ , $D_{y}$ and $D_{z}$ are diffusion coefficients of their respective dimensions ( $m^{2} ∕ s$ ) and $u_{x}$ , $u_{y}$ and $u_{z}$ are the advective flow in $x$ , $y$ and $z$ dimensions respectively ( $m ∕ s$ ). The equation given in Eq. ( 27 ) shows the solution of system that has no boundaries, where the chemicals are free to disperse or diffuse in any direction possible. To create a boundary effect on the propagation, the Cartesian definition of the function needs to be converted to cylindrical coordinates which the following section focuses on.

Figure 4.3: A descriptive diagram of the model used in the study. At the initial stage of the experiment (t = 0

s

) a mass is injected into the environment. Once the mass is injected, to simulate the effect of physical boundary of the environment (

A^{(0)}

), additional gas pulses are generated in the

y

-coordinates with positive sides being

y = 2 L

y = 4 R

. . .

y = 2 n R

(

A^{(1, 2, . ., \infty)}

). This is also carried out in the

x

-coordinates. As transmission evolves, the gas pulses transcend the boundary of the propagation medium, at which point the mirror pulses are added to the actual transmission to create the effect of the boundary.

4.2.4 The Radial-Advective-Diffusion Equation

Eq. ( 27 ) represents the concentration function in 3D Cartesian space and to describe the cylindrical geometry of the transmission medium, the equation is converted to cylindrical coordinates with the following transformations:

\begin{align} x & = R_{M} \cos 𝜃, y = R_{M} \sin 𝜃 z = z, & (4.28a) \\ D_{x} & = D_{y} = D_{T}, D_{z} = D_{L}, & (4.28b) \\ u_{x} & = u_{y} = u_{r}, u_{z} = u_{z} . & (4.28c) \end{align}

where $R_{M}$ is the radius of the propagation medium (m). Following this conversion process, Eq. ( 27 ) can be written in its cylindrical form:

\begin{matrix} c (R_{M}, 𝜃, z, t) = \frac{M}{\sqrt{{(4 π D_{T}^{2} D_{L} t)}^{3}}} \\ \exp (- \frac{{(R_{M} \cos 𝜃 - u_{r} t)}^{2}}{4 D_{T} t} - \frac{{(R_{M} \sin 𝜃 - u_{r} t)}^{2}}{4 D_{T} t} - \frac{{(z - u_{z} t)}^{2}}{4 D_{L} t}) . & (4.29) \end{matrix}

This equation can be further simplified by using trigonometric identities (i.e., $R_{M}^{2} = {R_{M}}^{2} \cos^{2} 𝜃 + {R_{M}}^{2} \cos^{2} 𝜃$ ) and omitting the radial advective flow ( $u_{r} = 0$ ) to the following expression:

c (R_{M}, z, t) = \frac{M}{\sqrt{{(4 π D_{T}^{2} D_{L} t)}^{3}}} \exp (- \frac{{R_{M}}^{2}}{4 D_{T} t} - \frac{{(z - u_{z} t)}^{2}}{4 D_{z} t}) .

(4.30)

4.2.4.1 Boundary Conditions

To create a boundary condition for this function, the method of mirror images is used. This is a mathematical tool for solving PDEs by adding the mirror image of the function with respect to the symmetry hyperplane.

For example to have a boundary at $x = x_{0}$ in 1D, the same function is added at $x = 2 x_{0}$ . This ensures that the change of concentration at the defined boundary $x_{0}$ equal to zero (i.e., zero flux at the radial boundary of the pipe). However, as the transmission evolves, more images are needed to maintain the accuracy of the function. Therefore, continuing the example, mirror images are added at distances $x = 4 x_{0}$ , $x = 6 x_{0}$ , ... .

In the environment used in this study, there is only the radial no-flux boundary at $R_{M}$ , where $R_{M}$ is the radius of the boundary ( $m$ ). If the transmission of particles is assumed to be in $z$ -direction, the boundary condition can be expressed as:

{\frac{\partial c}{\partial r} |}_{r = R_{M}} = 0 .

(4.31)

4.2.4.2 Absorption/Desorption Process

Based on the boundary condition described in Eq. ( 4.31 ), the mirror functions can be implemented to the concentration function in 3D which is given as:

c (r, z, t) = \frac{M}{\sqrt{{(4 π t)}^{3} D_{T}^{2} D_{L}}} \sum_{n = 0}^{\infty} \exp (- \frac{{(r - 2 n R_{M})}^{2}}{4 D_{T} t} - \frac{{(z - u_{z} t)}^{2}}{4 D_{L} t}),

(4.32)

where $n$ is the summation index. As can be seen in the equation above, the function is independent from $𝜃$ as it possesses angular symmetry. By integrating the concentration function with respect to the cylindrical volume element the particles that are present in the environment, ( $𝜃_{E}$ ) can be calculated.

𝜃_{E} = \underset{V}{∭} c d z r d r d 𝜃 .

(4.33)

As the system has no sink/source ( $R = 0$ ), the chemicals that are used in the transmission can either be in the environment ( $𝜃_{E}$ ) or have been absorbed by the detector ( $𝜃_{1}$ ). Therefore, both the aforementioned mass values must add upto the initial introduction of mass in the beginning of the transmission,

M = 𝜃_{E} + 𝜃_{1} .

(4.34)

The mass absorbed by the detector ( $𝜃_{1}$ ), can be calculated by subtracting from the initial mass ( $M$ ) [2], [4],

\begin{align} 𝜃_{1} (r, 𝜃, z, t) & = M - 𝜃_{E} (r, 𝜃, z, t), & (4.35a) \\ 𝜃_{1} (r, 𝜃, z, t) & = M - \int_{0}^{2 π} \int_{0}^{R} \int_{0}^{L} c (r, 𝜃, z, t) d z r d r d 𝜃, & (4.35b) \end{align}

where $x_{d}$ is the distance between the transmitter and the detector ( $m$ ). The solution to this equation, which expresses the absorbed particles by the detector, can be expressed as:

\begin{matrix} 𝜃_{1} (R_{M}, x_{d}, t) = M - \frac{M}{i \sqrt{4 D_{T} t}} [erf (\frac{u_{z} t}{\sqrt{4 D_{L} t}}) + erf (\frac{x_{d} - u_{z} t}{\sqrt{4 D_{L} t}})] \\ \sum_{n = 0}^{\infty} \exp (- n^{2} \frac{{R_{M}}^{2}}{D_{T} t}) {i \sqrt{D_{T} t} [1 - \exp ((4 n - 1) \frac{{R_{M}}^{2}}{4 D_{T} t})] + n R \sqrt{π} \exp (n^{2} \frac{{R_{M}}^{2}}{D_{T} t}) \\ \times [e r fi (n \frac{i R_{M}}{\sqrt{D_{T} t}}) - e r fi (\frac{(2 n - 1)}{2} \frac{i R_{M}}{\sqrt{D_{T} t}})]}, & (4.36) \end{matrix}

where $i$ is the imaginary unit with the identity $i^{2} = - 1$ and $e r f i (\cdot)$ is the imaginary error function with the following identity [345]:

e r fi (x) = - i e r f (i x) = \frac{2}{\sqrt{π}} \int_{0}^{x} e^{t^{2}} d t .

(4.37)

Once the chemicals are absorbed by the detector, the removal process can be initiated. To begin the calculation of the desorption process the particles that have been absorbed need to be quantified. To achieve this, the travel time of the signal has to be taken into account. As the chemical travels long distances, the response time is also delayed considerably and therefore the removal of particles from the detector to the outside environment is also delayed by the same amount of time. Therefore, the particles that are absorbed by the detector ( $M_{R}$ ) can be calculated as:

M_{R} = 𝜃_{1} (R_{M}, L, T_{S} + t_{e m p}) - 𝜃_{1} (R_{M}, L, t_{e m p}),

(4.38)

where $T_{S}$ is the symbol period ( $s$ ) and $t_{e m p}$ is the empirically measured time for the detection of chemical with respect to distance. The emprical fitting of this equation is given below and the fitting process can be seen in Figure 4.4 .

Figure 4.4: Experimentally measured chemical detection with comparison to empirical fitting (linear

R^{2} = 0.9891

\begin{align} t_{e m p} (L) = p_{1} L + p_{2}, & (4.39a) \\ p_{1} = 27.6, p_{2} = 16.51, & (4.39b) \end{align}

where $p_{1}$ and $p_{2}$ are the fitting parameters to the empirical fitting function. This parameter would change depending on the chemical that is used for sending information and in this study the function is based on acetone being the signal chemical. Based on these preliminary definitions, the removal of particles from the detector to the outside environment can be defined by the following function:

\begin{matrix} 𝜃_{0} (R_{M}, L_{M}, t) = \frac{M_{R}}{i \sqrt{4 D_{T} t}} [erf (\frac{u_{z} t}{\sqrt{4 D_{L} t}}) + erf (\frac{L_{M} - u_{z} t}{\sqrt{4 D_{L} t}})] \\ \sum_{n = 0}^{\infty} \exp (- n^{2} \frac{{R_{M}}^{2}}{D_{T} t}) {i \sqrt{D_{T} t} [1 - \exp ((4 n - 1) \frac{{R_{M}}^{2}}{4 D_{T} t})] + n R_{M} \sqrt{π} \exp (n^{2} \frac{{R_{M}}^{2}}{D_{T} t}) \\ \times [e r fi (n \frac{i R_{M}}{\sqrt{D_{T} t}}) - e r fi (\frac{(2 n - 1)}{2} \frac{i R_{M}}{\sqrt{D_{T} t}})]}, & (4.40) \end{matrix}

where $L_{M}$ is the distance between the membrane and the detector ( $m$ ). It must be noted that unlike the absorption process, where chemicals travel long distances to reach the detector, in the desorption process the chemical propagation begins from the detector membrane and ends at the outside environment ( $L ≫ L_{M}$ ). A detailed description of introduction/removal of particles can be seen in [2], [4] and a diagram of the model used in the study is presented in Figure 4.3 .

4.2.5 Calculation of the coefficient of Diffusivity

To calculate the longitudinal diffusivity coefficient of the propagation ( $D_{L}$ ), which plays a pivotal role in this type of communication transmission, the characteristic properties of the fluid motion must be established.

In a communication where particles are propagated through a medium, the main propeller of these particles is the volumetric flow rate ( $Q$ ). This is defined as the amount of volume transported in a given amount of time ( $m^{3} ∕ s$ ) and the velocity parameter ( $u$ ) can be obtained by dividing the volumetric flow rate by the cross-sectional area of the tube ( $A$ ):

u_{m a x} = \frac{Q}{A} = \frac{Q}{4 π R^{2}}, u_{a v g} = \frac{1}{2} u_{m a x} .

(4.41)

After obtaining the velocity parameter, the next characteristic of a fluid motion to be established is whether the motion is laminar or turbulent. To calculate this value, the Reynolds number (Re) is used. The equation for Reynolds number can be seen below [348]:

Re = \frac{ρ u_{a v g} L_{D}}{η},

(4.42)

where $ρ$ is the density of the fluid ( $k g ∕ m^{3}$ ), $u_{a v g}$ is the mean velocity of the fluid (m/s), $L_{D}$ is the characteristic linear dimension of the environment ( $m$ ). For the system in question it is the diameter of the pipe, and $η$ is the dynamic viscosity of the fluid ( $N \cdot s$ / $m^{2}$ ).

4.2.5.1 Entrance Length

Entrance length is defined as the distance a flow travels after entering a pipe before the flow becomes fully developed. Since the Reynolds number is low (Re $<$ 2000) the flow can be considered laminar and the entrance length for the system is calculated by the the following equation [348]:

L_{E} = 0.1 R_{M} R e .

(4.43)

4.2.5.2 Longitudinal Diffusivity

Longitudinal diffusivity is defined as diffusion paralel to the advective vector. As the flow is laminar (Re $<$ 2000) the longitudinal diffusivity can be calculated as [349]:

D_{L} = D_{m} + (\frac{{u_{a v g}}^{2} R_{M}^{2}}{48 D_{m}}),

(4.44)

where $D_{m}$ is the molecular diffusion ( ${c m}^{2} ∕ s$ ). The derivation of this equation can be seen in Appendix B.

4.2.5.3 Transverse Diffusion

The presence of the membrane affects the transverse diffusion more profoundly than longitudinal since the main propagator of motion along the radial axis is diffusion rather that advection aided diffusion seen in longitudinal diffusion. To calculate the coefficient, mean squared displacement is used [350]:

D_{T} = \lim_{t \to \infty} \frac{d}{d t} \sum_{n = 1}^{N} \frac{1}{2 N} ((x_{n} (t) - x_{n} (0))),

(4.45)

where $x_{n} (0)$ is the initial position coordinate of the gas molecules and $x_{n} (t)$ is the position coordinate of the gas molecule after time, $t$ .

Figure 4.5: A diagram of the algorithm used to simulate on-off keying (OOK) for macro-scale molecular communications.

4.3 Simulation Framework

Based on the equation derived from Eq. ( 4.20 ), a simulation of the macro-scale molecular communications is developed. The flowchart for the simulation can be seen in Figure 4.5 and the description of the simulation is as follows;

1.

A state array of $Γ$ and a frequency array $k$ with a length of $j$ is generated. The state array describes the possible states (i.e., in OOK there are two states; 0 and 1) and the frequency array describes the duration of occurrence of each state before changing into a new state.

2.

An initial definition of the captured mass ( $𝜃_{A}$ ) is generated and given the value of $𝜃_{A} = 0$ before any transmission commences.

3.

The system checks the array to see if the first transmitted state in the array is 0. This is done since the definition of the removal of particles ( $𝜃_{0}$ ) relies on the mass absorbed by the detector ( $M_{R}$ ). Therefore, the removal of particles is not initiated and the captured mass value is stays $𝜃_{A} = 0$ and the simulation continues to the next state in the array. This process continues until the system detects a 1 in the state array $(Γ)$ and moves the transmission window by the period of $T_{S}$ times the frequency array with the same index. This check is carried out by the $m$ parameter and once the simulation switches from 0 to 1 for the first time the m parameter is assigned the value of $m = 1$ and stops this loop until the termination of the simulation.

𝜃_{A} = [0_{1, (x, k (1) \times T_{S})}]

(4.46)

4.

The system checks the states in the array $Γ (j)$ . If the bit value is 1, the system introduces the chemicals into the system. The duration of the bit-1 is dependent on the value of the frequency array $k (j)$ with the same index value. For example, for a symbol duration of $T_{S}$ with three bit-1s, the total mass accumulated by the system is $𝜃_{1} (x, 3 T_{S})$ with each bit-1 value having $𝜃_{1} (x, T_{S})$ , $𝜃_{1} (x, 2 T_{S})$ and $𝜃_{1} (x, 3 T_{S})$ respectively. Based on this behaviour the captured mass is updated accordingly.

𝜃_{A} = [0_{1, (x, k (1) \times T_{S})} 𝜃_{1} (x, k (2) \times T_{S})]

(4.47)

This process can be seen in detail in Figure 4.6 .

5.

After the system has accumulated particles from the bit-1 transmission ( $𝜃_{1}$ ), when the next state is 0 and based on the value the system flushes the chemicals relative to the particles already in the system ( $M_{R}$ ) and continues to remove the particles from the detector until state array moves to the next value.

𝜃_{A} = [0_{1, (x, k (1) \times T_{S})} 𝜃_{1} (x, k (2) \times T_{S}) 𝜃_{0} (x, k (3) \times T_{S})]

(4.48)

6.

When the next bit-1 is to be transmitted, the leftover chemicals are added to the particle introduction.

7.

The operation continues until the count operator ( $j$ ) reaches the length of the state and frequency arrays.

𝜃_{A} = [0_{1, (x, k (1) \times T_{S})} 𝜃_{1} (x, k (2) \times T_{S}) \dots 𝜃_{Γ (j)} (x, k (j) \times T_{S})]

(4.49)

8.

Before the simulation concludes, additive white Gaussian noise (AWGN) is added to the simulated signal [3].

The pseudo-code for an OOK transmission 3D environment is given overleaf. The first part of the algorithm converts a given bit array into two different arrays: (i) state array (ii) frequency array. The second part of the algorithm is to simulate the transmission based on the two generated arrays.

Figure 4.6: A representative diagram of how the transmission is simulated with an example transmission of a bit sequence of 011100111000 with states of the transmission

(Γ_{S})

shown above the transmission. The first part of the simulation is to analyse the sequence based on the states. In this context the states are defined as a bit value, which in this example are 0 and 1. In this example there are five states which are 0-1-0-1-0 with durations of

1 T_{S}

3 T_{S}

2 T_{S}

3 T_{S}

3 T_{S}

. After this assessment the system carries out the following procedures to initiate the simulation. In this example at time-point (1) , the detector starts absorbing particles with the absorbing function

𝜃_{1}

and this function continues until the time period of the state concludes at (4) in which the duration of the state is shown as the feedback loop to the state itself, with each bit-1 value having the absorbed mass value of

𝜃_{1} (x, T_{S})

𝜃_{1} (x, 2 T_{S})

and

𝜃_{1} (x, 3 T_{S})

, respectively. When the time duration passes the time mark (4) , the removal function (

𝜃_{0}

) initiates and starts removing the particles from the detector based on how many particles it has absorbed in the previous state.

4.4 Conclusions

In this chapter, the mathematical description that is used throughout the thesis is described. The process of particle transport across an environment can be explained by the use of the advective-diffusion equation (ADE). The model is based on the mass exchange between the detector and the environment. Depending on the environment (i.e., boundary, open space) the equation that describes the evolution can change, and are described for both open distance and closed distance. In addition, the diffusion parameter for closed-distance transmission is investigated in detail. The chapter concludes by describing the simulation framework used throughout the study.

Three chapters follow that use the experimental setup described in Chapter 3 and theoretical model developed in Chapter 4. The next chapter describes the fundamental parameters of macro-scale molecular communication. These parameters play a pivotal role in describing the limit of the communications. From these limitations, reliable communication methods can be developed and analysed which are studied in Chapter 6. Chapter 7 focuses on further improving the communication speed and the mutual information of the system by using multiple chemicals for simultaneous transmission.

4 Theory of Macro-Scale Molecular Communications