# Ponderomotive force driven mechanism for electrostatic wave excitation and
energy absorption of Electromagnetic waves in overdense magnetized plasma

Laxman Prasad Goswami∗1, Srimanta Maity∗2, Devshree Mandal†3, Ayushi
Vashistha†4 and Amita Das∗5 ∗Department of Physics, Indian Institute of
Technology Delhi, Hauz Khas, New Delhi-110016, India
†Institute for Plasma Research, HBNI, Bhat, Gandhinagar-382428, India
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The excitation of electrostatic waves in plasma by laser electromagnetic pulse
is important as it provides a scheme by which the power from the laser
electromagnetic (EM) field can be transferred into the plasma medium. The
paper presents a fundamentally new ponderomotive pressure-driven mechanism of
excitation of electrostatic waves in an overdense magnetized plasma by a
finite laser pulse. Particle-in-cell (PIC) simulations using the EPOCH-4.17.10
framework have been utilized for the study of a finite laser pulse interacting
with a magnetized overdense plasma medium. The external magnetic field is
chosen to be aligned parallel to the laser propagation direction. In this
geometry, the electromagnetic wave propagation inside the plasma is identified
as whistler or R and L waves. The group velocity of these waves being
different, a clear spatial separation of the R and L pulses are visible. In
addition, excitation of electrostatic perturbation associated with the EM
pulses propagating inside the plasma is also observed. These electrostatic
perturbations are important as they couple laser energy to the plasma medium.
The excitation of electrostatic oscillations are understood here by a
fundamentally new mechanism of charge separation created by the difference
between the ponderomotive force (of the electromagnetic pulse) felt by the two
plasma species, viz., the electrons and the ions in a magnetized plasma.

  * April 2021

##  1 Introduction

It is well known that the interaction of electromagnetic wave (EM) with
magnetized plasma is considerably rich compared to un-magnetized plasma [1, 2,
3, 4]. This regime, however, has not yet been explored in the laboratory
laser-plasma experiments due to the magnitude of the magnetic field required
to make even the lighter electron species of the plasma magnetized at the
laser frequency. With the recent technological advancements, magnetic fields
of the order of kilo Tesla have now been achieved [5]. This is typically a
factor of $10$ smaller than the required magnetic field strength for making
the lighter electron species magnetized at the $CO_{2}$ laser frequency.
Recently a theoretical mechanism has been proposed for the generation of much
higher magnetic fields [6]. It appears, therefore, that there is a good
likelihood for the technological advancements to catch up in near future to
have the regime of magnetized electron response accessible in laser-plasma
interaction experiments. The other regime where even the heavier ions become
magnetized appear considerably difficult to achieve in laboratory laser-plasma
experiments in near future. Both, these regimes, however, are of importance in
astrophysical scenarios where an EM wave can readily encounter a strongly
magnetized plasma [7, 8]. It is, therefore, important to understand the
underlying physics of the EM wave interaction with a strongly magnetized
plasma in detail. It is important to note that Particle - In - Cell (PIC)
simulations have in recent decades developed a sophistication to capture the
physics of laser (or EM wave ) plasma interaction phenomena in their pristine
glory [9, 10]. We have, therefore, chosen to carry out PIC studies to
investigate the interaction of laser with a magnetized plasma which cannot be
explored as yet in laboratory experiments. There exists many familiar schemes
that facilitate the laser energy absorption in the plasma [11, 12]. For
instance, $\vec{J}\times\vec{B}$, resonance absorption, Brunel heating schemes
etc. [13, 14, 15, 16, 17, 18, 19].

It has been shown in earlier studies [4, 20, 21] concerning the magnetized
laser plasma interaction that when the external magnetic field is directed
normal to both the propagation direction and the direction of laser electric
field (X mode geometry) excitation of electrostatic lower Hybrid wave are
observed as a result of the difference between the $\vec{E}\times\vec{B}$
drifts of ions and electrons. This electrostatic coupling was responsible for
a new mechanism of laser energy absorption. In the O mode geometry though the
EM wave is unable to penetrate the overdense plasma, it generates harmonics
that propagate inside the plasma if appropriate conditions are satisfied [22].
In this work, we restrict ourselves to a geometry for which the external
magnetic field is directed along the propagation direction of the laser pulse.

The paper has been organized as follows. Section 2 contains the simulation
details. In section 3 we present the observations wherein we observe an
electrostatic excitation in the plasma trailing the transmitted
electromagnetic signal inside the plasma. The results are analyzed in section
4. We identify the mechanism of the generation of electrostatic fluctuations
in the plasma to the difference in the ponderomotive force of the laser field
felt by the two species (electrons and ions) in the presence of an external
magnetic field [23, 24, 25, 26, 27]. Section 5 provides the summary and
conclusion.

##  2 Simulation details

We have employed EPOCH-4.17.10 for carrying out 1-D particle-in-cell (PIC)
simulations to study the interaction of a laser with plasma in the presence of
an external magnetic field [EPOCH]. The magnetic field is directed along the
laser propagation direction. The schematic of the simulation geometry (not to
scale) has been shown in Fig.1. The laser is propagating along $\hat{x}$
direction and the applied external magnetic field is also along $\hat{x}$. The
electric field of laser is directed along $\hat{y}$. A 1-D simulation box with
dimension $L_{x}=1000\mu m$ has been chosen. The number of particles per cell
are taken to be $200$. The plasma boundary starts from $x=0$. There is vacuum
between $x=-150$ to $0$ $\mu m$. At $x=0$, there is a sharp plasma vacuum
interface. The spatial resolution is taken as $30$ cells per electron skin
depth and it corresponds to the grid size $\Delta x=0.01\mu m$. The laser is
incident on the plasma target from left side. We consider a short-pulse laser
of wavelength $\lambda_{l}=9.42$ $\mu m$ (frequency $\omega_{l}=0.2\times
10^{15}rad/s$), which corresponds to $CO_{2}$ laser. The choice of long
wavelength $CO_{2}$ laser reduces the requirements of the externally applied
magnetic field for allowing magnetized response of plasma. The laser profile
is Gaussian with the peak intensity of $I=3.5\times 10^{19}Wm^{-2}$ . Boundary
conditions are taken as absorbing in the longitudinal direction. The number
density of the plasma is taken to be $n_{0}=3\times 10^{26}m^{-3}$ for which
the electron plasma frequency is equal to $10^{15}rad/s$. Thus, the plasma is
overdense for the incident laser pulse.

Figure 1: Schematic (not to scale) showing laser energy being coupled into
plasma via excitation of L and R waves in the system. The external magnetic
field $B_{0}$ is along laser propagation direction ($\hat{x}$). The laser
electric field is in along $\hat{y}$ while the laser magnetic field is in the
$\hat{z}$ direction. The laser frequency $\omega$ is taken less tha the ion
cyclotron frequency $\omega_{ci}$.

The dynamics of both electrons and ions has been followed. To reduce the
computational time, we carried out the simulations at a reduced mass of ions,
which is taken to be $25$ (or $100$ for a few selected runs) times heavier
than electrons (i.e. $m_{i}=25m_{e}$, where $m_{i}$ and $m_{e}$ represent the
rest mass of the ion and electron species respectively). The external magnetic
field is chosen in such a way that for simulations in case(A) the condition
$\omega_{ci}<\omega_{l}<\omega_{ce}$ is satisfied. To satisfy this condition
the value of magnetic field has been taken to be $2.5m_{e}c\omega_{pe}e^{-1}$.
For case (B) we have $\omega_{l}<\omega_{ci}<\omega_{ce}$ which corresponds to
magnetic field of $10m_{e}c\omega_{pe}e^{-1}$.

In this work, we identify a new mechanism of electrostatic wave excitation in
the plasma by laser which depends on the difference between the ponderomotive
force of the EM radiation felt by the electron and ions in the presence of an
external magnetic field. To distinguish it from all other mechanisms we choose
a specific simulation configuration. The laser is chosen to be incident on the
sharp plasma surface in the normal direction which ensures the absence of
resonance and vacuum heating processes. Also, the role of
$\vec{J}\times\vec{B}$ electron heating is made negligible by choosing the
laser intensity in the non-relativistic domain of
$a_{0}(=eE/m_{e}\omega_{l}c)<1$ ($a_{0}=0.6$ ). The difference in
$\vec{E}\times\vec{B}$ drift velocities felt by the electron and ions in the
presence of external magnetic field reported in [4, 20] will also not be
operative in the longitudinal direction as the applied magnetic field is also
directed along $\hat{x}$. We have provided the chosen values of the laser and
plasma simulation parameters in Table-1.

Table 1: Values of simulation parameters in normalized and corresponding standard units for $m_{i}=25$ Parameters | Values in SI unit | Values in SI unit
---|---|---
| Case A | Case B
Plasma Parameters
$n_{0}$ | $3\times 10^{20}cm^{-3}$ | $3\times 10^{20}cm^{-3}$
$\omega_{pe}$ | $10^{15}rad/s$ | $10^{15}rad/s$
$\omega_{pi}$ ($M/m=25$) | $0.2\times 10^{15}rad/s$ | $0.2\times 10^{15}rad/s$
Laser Parameters
$\omega_{l}$ | $0.2\times 10^{15}rad/s$ | $0.2\times 10^{15}rad/s$
$\lambda_{l}$ | $9.42\mu m$ | $9.42\mu m$
Intensity | $3.5\times 10^{19}W/m^{2}$ | $3.5\times 10^{19}W/m^{2}$
External Magnetic Field Parameters
$B_{0}$ | $2.5m_{e}c\omega_{pe}e^{-1}$ | $10m_{e}c\omega_{pe}e^{-1}$
$\omega_{ce}$ | $2.5\times 10^{15}rad/s$ | $10\times 10^{15}rad/s$
$\omega_{ci}$ | $0.1\times 10^{15}rad/s$ | $0.4\times 10^{15}rad/s$

##  3 Observations

We have considered two cases in our simulation. One for which the magnetic
field is such that only the lighter electron species is magnetized at the
laser frequency and another where both the electron and ion species are
magnetized. We refer to these two cases as case (A) and case (B) respectively.
Since the applied magnetic field is along the direction of laser propagation,
for case (A) one observes the EM wave to get converted to the whistler (R)
wave as it enters the plasma medium. In case(B) both L and R waves form and
are seen to get separated spatially as their group speeds differ in the plasma
medium. We also observe the formation of electrostatic disturbances in the
plasma medium.

###  3.1 Case (A)

In this case, the laser frequency satisfies the condition
($\omega_{ci}<\omega_{l}<\omega_{ce}$) such that only the lighter electron
species are magnetized. In Fig. 2 we show the evolution of various fields as a
function of $x$ at different times (viz., $t=0$, $t=1ps$, $t=2.0ps$ and
$t=4ps$). The laser magnetic field, shown in the Fig. 2(a), (b), (c) and (d)
shows the reflection and transmission of the laser field from the plasma
surface. The incident laser is plane-polarized with its electric field
directed along the $y$ axis. As the transmitted signal propagates inside the
plasma it leaves behind a trail of electrostatic disturbances as evident from
the subplots of $E_{x}$ and $J_{x}$ (the $x$ component of the electric field
and current density respectively). The unbalanced ion (red color) and electron
(blue) densities are also evident from the subplots of the Fig.7(f), (g) and
(h). We would now like to characterize these disturbances (electromagnetic as
well as electrostatic observed in the plasma).

Figure 2: Time evolution of spatial variation of $z$-component of magnetic
field $B_{z}$, perturbation of electron and ions density ($n$), $x$-component
of electric field $E_{x}$ and $x$-component of total current $J_{x}$ with
respect to $\hat{x}$ for case (A), showing that as the laser interacts with
magnetized plasma, there is finite magnitude of perturbed $B_{z},E_{x}$ and
$J_{x}$ inside the plasma. We can infer from these plots that a part of laser
energy has been coupled into plasma.

We observe that though the incident laser was plane-polarized with its
electric field component directed along $y$, the transmitted electromagnetic
signal has both $y$ and $z$ components (Fig.3). Furthermore, from the phase
lag between the $E_{y}$ and $E_{z}$ components (Fig.3 (a), (b)) we can
identify as the right circularly polarised EM radiation. The Fast Fourier
Transform (FFT) of the $B_{z}$ signal in time at $x=200\mu m$ and in space at
$t=2ps$ has been shown in Fig.4 (a), (b) respectively. The frequency peak
occurs at the laser frequency $0.2\times 10^{15}rad/sec=0.2\omega_{pe}$ and
the wavenumber from the spatial FFT is found to be $k=1.04\times 10^{6}m^{-1}$
(which is distinct from the laser wavenumber being $0.67\times 10^{6}m^{-1}$).

The R and L-wave dispersion relations [29] are as follows:

$n^{2}=\frac{c^{2}k^{2}}{\omega^{2}}=1-\frac{\omega_{pe}^{2}}{\left(\omega+\omega_{ci}\right)\left(\omega-\omega_{ce}\right)}$
(1)

$n^{2}=\frac{c^{2}k^{2}}{\omega^{2}}=1-\frac{\omega_{pe}^{2}}{\left(\omega-\omega_{ci}\right)\left(\omega+\omega_{ce}\right)}$
(2)

This frequency and wavenumber lies in the pass band for R-mode and stop band
for L-mode Fig.5 (a) and (b) respectively and satisfy the R wave dispersion
relation Eq.(1) of the EM wave inside the plasma. The EM disturbance in the
plasma thus clearly gets identified as an R wave.

Figure 3: Phase variation of electric field $E_{y}$ and $E_{z}$ at $x=200\mu
m$ (a) for case (A) and (c) for case (B). Phase variation of electric field
$E_{y}$ and $E_{z}$ at $t=2ps$ (b) for case (A) and (d) for case (B). It
evident that for case (A) the excited mode is Right Circularly Polarized
(R-mode) and for case (B) both Right and Left Circularly Polarized modes are
obtained. Also (c) and (d) indicate that the R and L modes are separated in
time and space. Figure 4: (a) and (c) indicate the Fast Fourier Transform
(FFT) for electromagnetic field $B_{z}$ in time at $x=200\mu m$. ((b) and (d))
indicate the FFT for $B_{z}$ in space at $t=2ps$. FFT in (a) and (b) are for
case (A). FFT in (c) and (d) are for case (B). The FFT is taken after the
laser has reflected from the vacuum plasma interface.

We now study the characteristic features of the electrostatic disturbance
which gets created in the plasma. The FFT in time at $x=200\mu m$ and in space
at $t=2ps$ have been shown in Fig.6 (a), (b) respectively. For the chosen
configuration ($k$ is parallel to the $B_{0}$), beside the R and L wave, there
is also an elctrostatic mode present with a frequency same as that of electron
plasma frequency. The wavenumber, however, shows two peaks at $3.04\times
10^{6}m^{-1}$ and $4.0\times 10^{6}m^{-1}$. We shall understand the physics
behind these observations in the next section.

Figure 5: Diagram of dispersion curves for case (A) ((a) R-mode (a) and (b)
L-mode) and for case (B) ((c) R-mode (c) and (d) L-mode). For R, L-mode
geometry, region-1 and region-3 are passbands while region-2 is stopband. For
case (A), electromagnetic wave frequency ($\omega=0.2\times 10^{15}rad/s$)
lies (a) in the passband of R-mode and (b) in the stop band of L-mode. For
case (B), electromagnetic wave frequency ($\omega=0.2\times 10^{15}rad/s$)
lies in the passband of both R-mode and L-mode. For case (A), only R-mode is
obtained, while for case (B), both R-mode and L-mode are obtained. Figure 6:
(a) and (c) indicate the Fast Fourier Transform (FFT) for electrostatic field
$E_{x}$ in time at $x=200\mu m$. ((b) and (d)) indicate the FFT for $E_{x}$ in
space at $t=2ps$. FFT in (a) and (b) are for case (A). FFT in (c) and (d) are
for case (B). The FFT is taken after the laser has reflected from the vacuum
plasma interface.

###  3.2 Case (B)

We now consider the laser frequency to satisfy the condition
($\omega_{l}<\omega_{ci}<\omega_{ce}$). For this case, both electrons and ions
would exhibit a magnetized response at the laser frequency. In Fig.7 we show
the evolution of various fields as a function of $x$ at different times
(viz.,$t=0$, $t=1.0ps$, $t=2.0ps$ and $t=3.0ps$) for this particular case. In
contrast to case(A), in this case, the transmitted magnetic field of the EM
wave shows two distinct pulses. The spatial separation between the two pulses
is observed to increase with time. In this case too as the transmitted EM
signal propagates inside the plasma it leaves behind a trail of electrostatic
disturbances as evident from the subplots of $E_{x}$ and $J_{x}$ (the $x$
component of the electric field and current density respectively).

To check the polarization of the transmitted EM signal we plot the $E_{y}$ and
$E_{z}$ component of the electric field at a particular $x$ location
($x=200\mu m$) as a function of time in Fig.3(c) and at particular time
($t=2ps$) as a function of space in Fig.3(d). It is clear from the plot that
the pulse which moves ahead has a right-hand circular polarization whereas the
trailing pulse has the Left-hand circular polarization. The FFT of the $B_{z}$
signal in space and time has been shown in Fig.4(c), (d) respectively. The
frequency peak occurs at the laser frequency $0.2\times
10^{15}rad/sec=0.2\omega_{pe}$ which lies in the passband of R-mode and L-mode
as evident from Fig.5 (c) and (d) respectively and the wavenumbers from the
FFT are found to be $k_{1}=0.73\times 10^{6}m^{-1}$ and $k_{2}=0.85\times
10^{6}m^{-1}$. This frequency and wavenumbers ($k_{1}$ and $k_{2}$) satisfy
the R wave and L wave dispersion relation Eq.(1) and Eq.(2) respectively of
the EM wave inside the magnetized plasma.

The FFT in space at $t=2ps$ and time at $x=200\mu m$ has been shown in Fig.6
(c), (d) for the electrostatic oscillations that get generated inside the
plasma. In this case too the frequency of the electrostatic disturbances is
observed to match with the electron plasma frequency. However, two wavenumber
peaks are obtained at $3.06\times 10^{6}m^{-1}$ and $3.43\times 10^{6}m^{-1}$.
We shall discuss the physics behind these observation in the next section.

Figure 7: Time evolution of spatial variation of electromagnetic field
$B_{z}$, density ($n$) perturbation of electron and ions, $x$-component of
electric field $E_{x}$ and $x$-component of total current $J_{x}$ with respect
to $\hat{x}$ for case (B), showing that as the laser interacts with magnetized
plasma, there is finite magnitude of perturbed $B_{z}$, $E_{x}$ and $J_{x}$.
We can infer from these plots that a part of laser energy has been coupled
into plasma.

##  4 The genesis of electrostatic excitation

It should be noted that for a pulsed EM field, a ponderomotive force,
proportional to the gradient of EM wave intensity is operative on the plasma
particles. In our 1-D problem with pulsed EM field, this force will be
operative along the direction of propagation of the EM pulse, and the
corresponding acceleration felt by the two species is given by the following
expression:

${{{\frac{\partial{u^{e,i}_{x}}}{\partial t}}=-\frac{e^{2}(1+\alpha^{2}\pm
2\alpha\omega_{c{e,i}}/\omega_{l})}{m_{e,i}^{2}\omega_{l}^{2}\left(1-\omega^{2}_{c{e,i}}/\omega_{l}^{2}\right)}\nabla
E^{2}(x)}}$ (3)

Here, the plus sign is for electrons, and the minus sign is for ions. Thus,
$\omega_{ce}=eB_{0}/m_{e}$ and $\omega_{ci}=eB_{0}/m_{i}$ are the magnitude of
the electron and ion cyclotron frequencies respectively. Where $m_{e}$ and
$m_{i}$ are respectively the electron and ion mass. Also, $\omega_{l}$ stands
for the laser frequency and $e$ represent the magnitude of electronic charge.
Here $\alpha$ takes the value of $0$ and $\pm 1$ for plane and right and left
circularly polarised EM (electromagnetic) waves respectively. The derivation
of Eq.(3) has been provided in A.2. In the absence of an applied magnetic
field, the above equation reduces to the standard form of (as given in A.1)

${{{\frac{\partial{u^{e,i}_{x}}}{\partial
t}}=-\frac{e^{2}(1+\alpha^{2})}{m_{e,i}^{2}\omega_{l}^{2}}\nabla E^{2}(x)}}$
(4)

which shows that the electrons experience a stronger ponderomotive
acceleration.

The physics of the electrostatic fluctuation driven by the difference in the
ponderomotive force experienced by electrons and ions can be understood as
follows. From Eq.(3) it is clear that the electrons and ions acquire an
additional velocity (other than the quiver velocity in the laser field) along
$x$ direction. The difference between the $x$ component of ion and electron
velocity generates $x$ component of current. The laser propagation direction
being along $x$, this current has a finite divergence. From the charge
continuity equation, viz.,

$\frac{\partial\rho}{\partial t}+\vec{\nabla}\cdot\vec{J}=0$ (5)

the current divergence leads to space charge fluctuations which have been
observed in our simulations (second row of Fig.2 and Fig.7). We have thus
identified a new ponderomotive force-driven mechanism of electrostatic wave
excitation in the plasma.

In the limit of strong magnetic field for which
$\omega_{l}<<\omega_{ci}<<\omega_{ce}$ (e.g. Case (B) of our simulations),
Eq.(3) reduces for the two species to

${{{\frac{\partial{u^{e,i}_{x}}}{\partial t}}=\frac{(1+\alpha^{2}\pm
2\alpha\omega_{c{e,i}}/\omega_{l})}{B_{0}^{2}}\nabla E^{2}(x)}}$ (6)

Clearly, in this limit of the strong magnetic field the acceleration of the
two species depends predominantly on the applied magnetic field and reduces as
$\sim 1/B_{0}^{2}$. The two species experience different acceleration only for
circularly polarised waves for which $\alpha=\pm 1$. For plane polarization
with $\alpha=0$ both the species experience identical ponderomotive
acceleration. In this case thus one would not generate any charge fluctuation.
When the laser frequency falls in the regime of
$\omega_{ci}<<\omega_{l}<<\omega_{ce}$ (case (A) of our simulations) the ion
would experience the unmagnetized ponderomotive expression provided by Eq.(4)
and electron would be governed by the Eq.(6).

Let us now try to analyze our observations from these considerations. For case
(A) since the ion and electron satisfy the unmagnetized and magnetized
expressions for the ponderomotive acceleration, for this case, there will
always be an unbalanced current (for the plane as well as circular
polarization). At the plasma vacuum surface the ponderomotive pressure of the
plane incident laser is operative and electrostatic fluctuations are evident
as clear from the zoomed view of the electrostatic field in Fig.2. This should
be contrasted with case (B) where there are no surface electrostatic
excitations. Both ion and electrons here satisfy Eq.(6) which is identical for
plane polarization of the incident laser pulse. The overall magnitude of the
electrostatic fluctuations in case (B) typically generated by the L and R mode
in the bulk of the plasma is also comparatively much low compared with case(A)
as evident from Fig.8.

Energy density plot in Fig.9(a) signifies the generation of electrostatic
field energy for case (A). It is also evident from Fig.9(b) that the
eelctrostatic field energy generation is more in case (A) than in case(B). We
have repeated our simulations with the increasing magnitude of the magnetic
field. In Fig.10 we provide a comparison of the evolution of the electrostatic
field energy in the plasma for different values of $B_{0}$. It can be seen
clearly, that the electrostatic energy in the plasma reduces with the
increasing magnetic field.

These studies thus show that the observations on electrostatic excitations and
absorption are consistent with the proposed mechanism driven by the difference
between the ponderomotive acceleration felt by the two species.

Figure 8: Time evolution of Electrostatic field variation for case (A) (red
desh line) and case (B) (blue line). It is evident from the plots that the
magnitude of electrostatic field generation due to ponderomotive force is more
in case (A) than in case (B)  Figure 9: The comparision for the time variation
of (a) for case (A) the mean of Electrostatic field energy density (red desh
line, left $y$ axis) and electromagnetic (blue line, right $y$ axis), and (b)
the mean of Electrostatic field energy density for case (A) (red desh line,
left $y$ axis) and case (B) (blue line, right $y$ axis). It is evident from
subplot (a) that electrostatic modes are generated along with the
electromagnetic field. Subplot (b) signifies that the electrostatic field
eenrgy density reduces with increase in external magnetic field and
electrostatic field energy is more in case(A) than that of in case (B). Figure
10: Figure shows the effect of varying external magnetic field on the
Electrostatic Energy. It can be observed from the graph that the energy
absorbed in the plasma medium reduces predominantly at higher magnetic fields.

##  5 Conclusion

We have carried out extensive PIC simulations with the help of EPOCH-4.17.10
for the propagation of laser pulse in a magnetized plasma. The direction of
the external magnetic field has been chosen along the laser propagation
direction. We demonstrate a novel mechanism of electrostatic wave excitation
leading to laser energy absorption in the plasma medium. The mechanism is
essentially driven by the difference between the ponderomotive acceleration
felt by electrons and ions.

##  6 Acknowledgment

The authors would like to acknowledge the EPOCH Consortium, for providing us
the access to the EPOCH-4.17.10 framework [EPOCH]. This research work has been
supported by the Core Research Grant No. CRG/2018/000624 of Department of
Scient and Technology (DST), Government of India. We also acknowledge support
from J. C. Bose Fellowship grant of AD (JCB-000055/2017) from the Science and
Engineering Research Board (SERB), Government of India. The authors thank IIT
Delhi HPC facility for computational resources. L.P.G also wishes to thank
Sathi Das for discussions.

##  7 References

## References

  * [1] VB Krasovitskii, VG Dorofeenko, VI Sotnikov, and BS Bauer.  Interaction of powerful laser pulse with magnetized plasma.  Physics of Plasmas, 11(2):724–742, 2004.
  * [2] Davoud Dorranian, Mikhail Starodubtsev, Hiromichi Kawakami, Hiroaki Ito, Noboru Yugami, and Yasushi Nishida.  Radiation from high-intensity ultrashort-laser-pulse and gas-jet magnetized plasma interaction.  Physical Review E, 68(2):026409, 2003.
  * [3] Pallavi Jha, Punit Kumar, Gaurav Raj, and Ajay K Upadhyaya.  Modulation instability of laser pulse in magnetized plasma.  Physics of plasmas, 12(12):123104, 2005.
  * [4] Ayushi Vashistha, Devshree Mandal, and Amita Das.  Excitation of lower hybrid and magneto-sonic perturbations in laser plasma interaction.  Nuclear Fusion, 61(2):026016, 2021.
  * [5] Daisuke Nakamura, A Ikeda, H Sawabe, YH Matsuda, and S Takeyama.  Record indoor magnetic field of 1200 t generated by electromagnetic flux-compression.  Review of Scientific Instruments, 89(9):095106, 2018.
  * [6] Ph Korneev, E d’Humières, and V Tikhonchuk.  Gigagauss-scale quasistatic magnetic field generation in a snail-shaped target.  Physical Review E, 91(4):043107, 2015.
  * [7] Christoph Federrath.  Magnetic field amplification in turbulent astrophysical plasmas.  Journal of Plasma Physics, 82(6), 2016.
  * [8] PK Shukla, G Feix, and L Stenflo.  Nonlinearly coupled electromagnetic ion-cyclotron and magnetosonic waves in astrophysical plasmas.  Astrophysics and space science, 147(2):383–388, 1988.
  * [9] John M. Dawson.  Particle simulation of plasmas.  Reviews of Modern Physics, 55(2):403–447, 1983.
  * [10] Charles K Birdsall.  Particle-in-cell charged-particle simulations, plus monte carlo collisions with neutral atoms, pic-mcc.  IEEE Transactions on plasma science, 19(2):65–85, 1991.
  * [11] PK Kaw.  Nonlinear laser–plasma interactions.  Reviews of Modern Plasma Physics, 1(1):1–42, 2017.
  * [12] Predhiman Krishan Kaw and JM Dawson.  Laser-induced anomalous heating of a plasma.  The Physics of Fluids, 12(12):2586–2591, 1969.
  * [13] JP Freidberg, RW Mitchell, Ro L Morse, and LI Rudsinski.  Resonant absorption of laser light by plasma targets.  Physical Review Letters, 28(13):795, 1972.
  * [14] Thomas H Stix.  Radiation and absorption via mode conversion in an inhomogeneous collision-free plasma.  Physical Review Letters, 15(23):878, 1965.
  * [15] SC Wilks, WL Kruer, M Tabak, and AB Langdon.  Absorption of ultra-intense laser pulses.  Physical review letters, 69(9):1383, 1992.
  * [16] L Chopineau, A Leblanc, G Blaclard, A Denoeud, M Thévenet, JL Vay, G Bonnaud, Ph Martin, H Vincenti, and F Quéré.  Identification of coupling mechanisms between ultraintense laser light and dense plasmas.  Physical Review X, 9(1):011050, 2019.
  * [17] F Brunel.  Not-so-resonant, resonant absorption.  Physical review letters, 59(1):52, 1987.
  * [18] WLl Kruer and Kent Estabrook.  J$\times$ b heating by very intense laser light.  The Physics of fluids, 28(1):430–432, 1985.
  * [19] Y Ping, R Shepherd, BF Lasinski, M Tabak, H Chen, HK Chung, KB Fournier, SB Hansen, A Kemp, DA Liedahl, et al.  Absorption of short laser pulses on solid targets in the ultrarelativistic regime.  Physical review letters, 100(8):085004, 2008.
  * [20] Ayushi Vashistha, Devshree Mandal, Atul Kumar, Chandrasekhar Shukla, and Amita Das.  A new mechanism of direct coupling of laser energy to ions.  New Journal of Physics, 22(6):063023, 2020.
  * [21] Atul Kumar, Chandrasekhar Shukla, Deepa Verma, Amita Das, and Predhiman Kaw.  Excitation of kdv magnetosonic solitons in plasma in the presence of an external magnetic field.  Plasma Physics and Controlled Fusion, 61(6):065009, 2019.
  * [22] Srimanta Maity, Devshree Mandal, Ayushi Vashistha, Laxman Prasad Goswami, and Amita Das.  Generation of higher harmonic by laser interacting with magnetized plasma.  arXiv preprint arXiv:2102.03104, 2021.
  * [23] R Singh and AK Sharma.  Ponderomotive acceleration of electrons by a whistler pulse.  Applied Physics B, 100(3):535–538, 2010.
  * [24] FV Hartemann, SN Fochs, GP Le Sage, NC Luhmann Jr, JG Woodworth, MD Perry, YJ Chen, and AK Kerman.  Nonlinear ponderomotive scattering of relativistic electrons by an intense laser field at focus.  Physical Review E, 51(5):4833, 1995.
  * [25] Anamika Sharma and VK Tripathi.  Ponderomotive acceleration of electrons by a laser pulse in magnetized plasma.  Physics of Plasmas, 16(4):043103, 2009.
  * [26] CS Liu and VK Tripathi.  Ponderomotive effect on electron acceleration by plasma wave and betatron resonance in short pulse laser.  Physics of plasmas, 12(4):043103, 2005.
  * [27] V Sazegari, M Mirzaie, and B Shokri.  Ponderomotive acceleration of electrons in the interaction of arbitrarily polarized laser pulse with a tenuous plasma.  Physics of plasmas, 13(3):033102, 2006.
  * [28] TD Arber, Keith Bennett, CS Brady, A Lawrence-Douglas, MG Ramsay, NJ Sircombe, P Gillies, RG Evans, Holger Schmitz, AR Bell, et al.  Contemporary particle-in-cell approach to laser-plasma modelling.  Plasma Physics and Controlled Fusion, 57(11):113001, 2015.
  * [29] TJ Boyd, TJM Boyd, and JJ Sanderson.  The physics of plasmas.  Cambridge University Press, 2003.
  * [30] K. Nishikawa and M. Wakatani.  Plasma physics.  Springer-Verlag,Berlin, 2000.

##  Appendix A General derivation for Ponderomotive Force in magnetized plasma
for elliptically polarized laser

To derive for the ponderomotive force, we consider the response of a
homogeneous magnetized plasma to a high frequency elliptically polarised field
whose amplitude is spatially varying along $\hat{x}$ given as:

$\vec{E}=E(x)\left(\hat{y}+i\alpha\hat{z}\right)e^{-i\omega t}$ (7)

Using the Maxwell equation

$\nabla\times\vec{E}=-\frac{\partial\vec{B}}{\partial t}$ (8)

We calculate the oscillating magnetic field as:

$\vec{B}_{1}=-\frac{i}{\omega}\frac{\partial E(x)}{\partial
x}\left(\hat{z}-i\alpha\hat{y}\right)e^{-i\omega t}$ (9)

We are considering a constant external magnetic field given by the following
expression:

$\vec{B}_{ext}=B_{x}\hat{x}+B_{y}\hat{y}+B_{z}\hat{z}=b_{x}B_{0}\hat{x}+b_{y}B_{0}\hat{y}+b_{z}B_{0}\hat{z}$
(10)

Where,

$b_{i}=\frac{B_{i}}{B_{0}}$

Thus the total magnetic field is summation of time varying oscillating
magnetic field and the constant external magnetic field given as:

$\vec{B}=\vec{B}_{ext}+\vec{B}_{1}$ (11)

Figure 11: Geometry for (a) RL mode (the external magnetic field $B_{0}$ is
applied along laser propagation direction $\hat{x}$ with oscillating laser
electric field along $\hat{y}$ direction), (b) X-mode (the external magnetic
field $B_{0}$ is applied along direction of laser magnetic field $\hat{z}$
with oscillating laser electric field along $\hat{y}$ direction), and (c)
O-mode configurations (the external magnetic field $B_{0}$ is applied along
direction of laser electric field $\hat{y}$).

If we neglect the electron - ion pressure terms, the force equation is given
as:

$\frac{\partial\vec{u}}{\partial
t}+\vec{u}\cdot\nabla\vec{u}=\frac{q_{p}}{m_{p}}\left(\vec{E}+\vec{u}\times\vec{B}\right)$
(12)

Where $p=e,i$ respectively indicate the electron and ion species. So $q_{e}$
and $q_{i}$ are the electron and ion charges. Similarly $m_{e}$ and $m_{i}$
are the electron and ion masses.

To lowest order in $|\vec{E}|$, $\vec{u}=\vec{u}^{h}$ we can ignore the non-
linear terms from the momentum equation

$\frac{\partial\vec{u}^{h}}{\partial
t}=\frac{q_{p}}{m_{p}}\left(\vec{E}+\vec{u}^{h}\times\vec{B}{{}_{2}}\right)$
(13)

After taking fourier transform of the above Eq.(13), and solving for
$u_{x}^{h}$, $u_{y}^{h}$ and $u_{z}^{h}$

$u_{x}^{h}=i\frac{\omega_{cp}v_{E}}{\beta\omega^{2}}\left[\left(b_{y}\alpha-
b_{x}b_{y}\frac{\omega_{cp}}{\omega}\right)+i\left(b_{z}-\alpha
b_{x}b_{z}\frac{\omega_{cp}}{\omega}\right)\right]$ (14)

$u_{y}^{h}=\frac{v_{E}}{\beta\omega}\left[\alpha
b_{y}b_{z}\frac{\omega^{2}_{cp}}{\omega^{2}}+i\left(1-b_{y}^{2}\frac{\omega^{2}_{cp}}{\omega^{2}}-b_{x}\alpha\frac{\omega_{cp}}{\omega}\right)\right]$
(15)

$u_{z}^{h}=-\frac{v_{E}}{\beta\omega}\left[\alpha\left(1-b_{z}^{2}\frac{\omega^{2}_{cp}}{\omega^{2}}\right)-b_{x}\frac{\omega_{cp}}{\omega}+ib_{y}c_{z}\frac{\omega^{2}_{cp}}{\omega^{2}}\right]$
(16)

Where,

$v_{E}=\frac{q_{p}E(x)}{m_{p}}e^{-i\omega t}$ (17)

$\omega_{cp}=\frac{q_{p}B_{0}}{m_{p}}$ (18)

$\beta=1-\frac{\omega_{cp}^{2}}{\omega^{2}}\left(b_{x}^{2}+b_{y}^{2}+b_{z}^{2}\right)=1-\frac{\omega_{cp}^{2}}{\omega^{2}}\frac{B_{ext}^{2}}{B_{0}^{2}}$
(19)

The electrons and ions are simply oscillating in the oscillating electric and
magnetic fields. By averaging the force equation over these rapid
oscillations, we obtain

$\frac{\partial\vec{u}^{s}}{\partial
t}=\frac{q_{p}}{m_{p}}\left<\vec{u}^{h}\times\vec{B}_{1}\right>_{t}-\left<\vec{u}\cdot\nabla\vec{u}\right>_{t}$
(20)

Where $\left<\right>_{t}$ denotes an average over high frequency oscillations.
The second term in Eq.(20) can be written as:

Substituting the known values in above equation, we have:

$\frac{q_{p}}{m_{p}}\left<\vec{u}^{h}\times\vec{B}_{1}\right>_{t}=-\frac{q_{p}^{2}}{2m_{p}^{2}\beta\omega^{2}}\left[\left(1-b_{y}^{2}\frac{\omega_{cp}^{2}}{\omega^{2}}\right)+\alpha^{2}\left(1-b_{z}^{2}\frac{\omega_{cp}^{2}}{\omega^{2}}\right)-2\alpha
b_{x}\frac{\omega_{cp}}{\omega}\right]\nabla|E|^{2}$ (21)

$\left<\vec{u}\cdot\nabla\vec{u}\right>_{t}=\frac{q^{2}_{p}\omega^{2}_{cp}}{2m^{2}_{p}\beta^{2}\omega^{4}}\left[\left(\alpha
b_{y}-b_{x}b_{y}\frac{\omega_{cp}}{\omega}\right)^{2}+\left(b_{z}-\alpha
b_{x}b_{z}\frac{\omega_{cp}}{\omega}\right)^{2}\right]\nabla|E|^{2}$ (22)

Using Eq.(21) and Eq.(22), we can write Eq.(20) as:

$\frac{\partial\vec{u}^{s}}{\partial
t}=-\frac{q_{p}^{2}\nabla|E|^{2}}{2m_{p}^{2}\beta\omega^{2}}\left[A_{1}+\alpha^{2}A_{2}-2\alpha
a\frac{\omega_{cp}}{\omega}+\frac{\omega_{cp}^{2}}{\beta\omega^{2}}A_{3}\right]$
(23)

Where,

$A_{1}=\left(1-b_{y}^{2}\frac{\omega_{cp}^{2}}{\omega^{2}}\right)$
$A_{2}=\left(1-b_{z}^{2}\frac{\omega_{cp}^{2}}{\omega^{2}}\right)$
$A_{3}=\left(\alpha
b_{y}-b_{x}b_{y}\frac{\omega_{cp}}{\omega}\right)^{2}+\left(b_{z}-\alpha
b_{x}b_{z}\frac{\omega_{cp}}{\omega}\right)^{2}$

###  A.1 Un-magnetized case

For unmagnetized case, we have $b_{x}=b_{y}=b_{z}=0$, thus $\beta=1$. So the
ponderomotive force in this case will be [30]:

$\frac{\partial\vec{u}^{s}}{\partial
t}=-\frac{q_{p}^{2}\nabla|E|^{2}}{2m_{p}^{2}\omega^{2}}\left(1+\alpha^{2}\right)$
(24)

###  A.2 RL-mode geometry

For RL-mode case, we have $b_{x}=1$ and $b_{y}=b_{z}=0$, thus
$\beta=1-\frac{\omega_{cp}^{2}}{\omega^{2}}$. So the ponderomotive force in
this geomrety is will be:

$\frac{\partial\vec{u}^{s}}{\partial
t}=-\frac{q_{p}^{2}\nabla|E|^{2}}{2m_{p}^{2}\omega^{2}}\frac{\left[1+\alpha^{2}-2\alpha\frac{\omega_{cp}}{\omega}\right]}{\left(1-\frac{\omega_{cp}^{2}}{\omega^{2}}\right)}$
(25)

###  A.3 X-mode geometry

For X-mode case, we have $b_{x}=b_{y}=0$ and $b_{z}=1$, thus
$\beta=1-\frac{\omega_{cp}^{2}}{\omega^{2}}$. So the ponderomotive force in
this geomrety is will be:

$\frac{\partial\vec{u}^{s}}{\partial
t}=-\frac{q_{p}^{2}\nabla|E|^{2}}{2m_{p}^{2}\omega^{2}}\frac{1}{\left(1-\frac{\omega_{cp}^{2}}{\omega^{2}}\right)}\left[1+\alpha^{2}\left(1-\frac{\omega_{cp}^{2}}{\omega^{2}}\right)+\frac{\omega_{cp}^{2}}{\omega^{2}}\frac{1}{\left(1-\frac{\omega_{cp}^{2}}{\omega^{2}}\right)}\right]$
(26)

###  A.4 O-mode geometry

For O-mode case, we have $b_{x}=b_{z}=0$ and $b_{y}=1$, thus
$\beta=1-\frac{\omega_{cp}^{2}}{\omega^{2}}$. So the ponderomotive force in
this geomrety is will be:

$\frac{\partial\vec{u}^{s}}{\partial
t}=-\frac{q_{p}^{2}\nabla|E|^{2}}{2m_{p}^{2}\omega^{2}}\frac{1}{\left(1-\frac{\omega_{cp}^{2}}{\omega^{2}}\right)}\left[\left(1-\frac{\omega_{cp}^{2}}{\omega^{2}}\right)+\alpha^{2}+\frac{\omega_{cp}^{2}}{\omega^{2}}\frac{\alpha^{2}}{\left(1-\frac{\omega_{cp}^{2}}{\omega^{2}}\right)}\right]$
(27)


# Properties of Hagen-Poiseuille flows in channel networks

Ana Filipa Valente  Rodrigo Almeida  Rui Dilão<EMAIL_ADDRESS>University of Lisbon, Instituto Superior Técnico, Nonlinear Dynamics Group,
Av. Rovisco Pais, 1049-001 Lisbon, Portugal.

###### Abstract

We derive the main properties of adaptive Hagen-Poiseuille flows in elastic
microchannel networks akin to biological veins in organisms. We show that
adaptive Hagen-Poiseuille flows successfully simulate key features of Physarum
polycephalum networks, replicating physiological out-of-equilibrium phenomena
like peristalsis and shuttle streaming, associated with the mechanism of
nutrient transport in Physarum. A new topological steady state has been
identified for asynchronous adaptation, supporting out-of-equilibrium laminar
fluxes. Adaptive Hagen-Poiseuille flows show saturation effects on the fluxes
in contractile veins, as observed in animal and artificial contractile veins.

Hagen-Poiseuille flows in networks, Biological networks, Shuttle streaming,
Peristalsis, Physarum.

###### pacs:

87.16.Uv, 87.17.Jj, 47.60.+i

##  I Introduction

Channel networks are ubiquitous transportation and communication structures in
natureTakahashi, (2014). This is the case of organ and tumour irrigation veins
Rubinow, (1973); Chang et al., (2019), nutrient transporting channels in plant
leaves Katifori et al., (2010), microorganisms with adaptive channel networks
like Physarum polycephalum Baumgarten et al., (2010); Oettmeier et al.,
(2017), transportation systems Tero et al., (2006); Hu and Cai, (2013) or
electrical networks Bohn and Magnasco, (2007). These systems can be described
using graphs and compared to formal systems like the travelling salesman
problem Zhu et al., (2018) or optimal Steiner tree geometries Tero et al.,
(2010). These are difficult network optimisation problems in graph theory and
have important applications in physics, biology and engineering.

Physarum polycephalum, also known as slime mould, is a protist which exhibits
complex behaviour. Its body consists of several nuclei (syncytium) and an
adaptive network of veins transporting endoplasmatic fluid from food sources
to different body regions. This network’s geometry and adaptability are
responsible for Physarum’s growth and motility Dove et al., (1985). The vein
walls of the network channels exhibit periodic transversal oscillations
regulated by food sources and other external stimuli Wohlfarth, (1979). Even
though slime mould lacks any neural system or centralised control, it shows
high-level behaviour, such as choosing the shortest path between food sources
in a maze Nakagaki et al., (2000).

To recreate and understand the dynamics of growth and movement of Physarum,
several different models have been developed to study the formation and
optimisation of network paths Tero et al., (2010); Almeida, (2021); Almeida
and Dilão, (2022); Valente, (2023). These models were based on the
phenomenology of Hagen-Poiseuille flows on graphs.

From all the features of the formation and adaptation of channel networks in
Physarum, some physical phenomena are poorly explained. This is the case of
the validation of Murray’s law at the branching points of channelsMurray,
(1926); Sherman, (1981), oscillations of channel radius and inversion of
stream velocities (peristalsis and shuttle streaming), conductivity saturation
in elastic vein networks, non-unicity of steady-state solutions of adaptive
Hagen-Poiseuille (H-P) flows, and suboptimal and optimal channel lengths in
steady-state solutions of H-P network flows. In these flows, adaptability is
associated with the transverse elasticity of channels.

This paper focuses on the properties of the adaptive H-P flows described
above, identifying the origin of these phenomena. We will follow the Hagen-
Poiseuille formalism, a Navier-Stokes type approach with a low Reynolds
number, imposing fluid mass conservation and energy minimisation as developed
by Almeida et alAlmeida, (2021); Almeida and Dilão, (2022). In the next
section, we describe the exact results of the temporal evolution of adaptive
H-P flows in channel networks embedded in $n$-dimensional Euclidian spaces. In
section III, we derive the properties of the steady states of the channel
flows and their spatial geometries, and we validate the Murray law for
adaptive H-P flows at steady state. In section IV, we compare the steady-state
H-P flow patterns with Physarum patterns. In section V, we analyse shuttle
streaming, and in section VI, we analyse peristalsis, both observed in out-of-
equilibrium flows. Finally, in section VII, we summarise the paper’s results.

##  II Adaptive Hagen-Poiseuille flows on graphs: exact results

The fluid flow in elastic straight channel networks of several systems, such
as in Physarum polycephalum or blood vessels in organisms, can be described as
the flow of viscous and incompressible fluids, with several sources and sinks
and adaptive conductivities.

A network of veins is described by an undirected graph
$\mathcal{G}=(\mathcal{V},E)$, embedded in a $n$ dimensional Euclidean space,
where $\mathcal{V}$ is the set of $N$ nodes or vertices with coordinates
$(x_{i}^{1},...,x_{i}^{n})$, with $i=1,...,N$, and $E$ is the set of $M$
straight edges $(i,j)$ connecting the nodes $i$ and $j$. We assume that the
graph $\mathcal{G}$ is connected, having only one component.

As we will describe out-of-equilibrium flows, the graph $\mathcal{G}$ is
generated through the Delaunay triangulation of a square lattice of size
$m\times m$, embedded in a unit square, where a white Gaussian stochastic
algorithm modified the coordinates of the interior points. The edges of
$\mathcal{G}$ represent elastic cylindrical channels that can expand or
contract transversely to the flow. This network has $K$ sources and $R$ sinks
at fixed nodes. The input or output fluid flux at node $j$ is represented by
$S_{j}$. If the network has a source at node $j$, then $S_{j}>0$. If the
network has a sink at node $j$, $S_{j}<0$. Otherwise, $S_{j}=0$. As the fluid
is incompressible,

$\sum_{j=1}^{N}S_{j}=\sum_{j:\text{sources}}S_{j}+\sum_{j:\text{sinks}}S_{j}=0,$
(1)

and the volume of fluid in the network is constant.

The edge connecting node $i$ to node $j$ has fixed length $L_{ij}$ and
(variable) radius $r_{ij}$. Edges have (variable) conductivities $D_{ij}=\pi
r_{ij}^{4}/8\eta$, where $\eta$ is the dynamic viscosity of the fluid, and
$\sqrt{D_{ij}}\simeq r_{ij}^{2}$ is proportional to the area of the channel.
The volume of the channel connecting node $i$ to node $j$ is $V_{ij}=\pi
r_{ij}^{2}L_{ij}=\sqrt{8\pi\eta}L_{ij}\sqrt{D_{ij}}$. As the fluid is
incompressible, the volume of fluid in the network is constant and is given by

$V=\sum_{(i,j)\in E}V_{ij}=\beta\sum_{(i,j)\in E}L_{ij}\sqrt{D_{ij}},$ (2)

where $\beta=\sqrt{8\pi\eta}$.

The Kirchhoff law determines the fluxes at every node $j$. As fluid flux is
conserved, we have

$\sum_{j:(i,j)\in E}Q_{ij}=S_{i}\quad,i=1,...,N.$ (3)

The H-P flux law states that the flux of fluid in a cylindric channel
connecting node $i$ to node $j$ depends on the pressures at the nodes and is
determined by the relation

$Q_{ij}=D_{ij}\frac{(p_{i}-p_{j})}{L_{ij}}.$ (4)

Biological elastic channels show transverse oscillations regulated by their
fluxRubinow and Keller, (1972); Wohlfarth, (1979). To describe these fluxes,
we consider the network of channels $\mathcal{G}$ filled with an
incompressible fluid, where the channel conductivities have a positive value.
Assuming that all the fluxes $S_{i}$ are known, by the H-P relationship (4),
the pressures $p_{i}$ at the vertices are obtained by solving the Kirchhoff
laws (3). Further assuming that the veins are elastic and the conductivity can
change in time, we have:

######  Theorem II.1 (Almeida and Dilão, (2022))

Consider a H-P flow in a channel network described by a connected graph
$\mathcal{G}$. Assuming that the volume $V$ of fluid in the network is
constant over time, then the conductivity of the elastic channels change in
time according to the adaptation law

$\frac{d}{dt}\sqrt{D_{ij}}=\alpha\frac{g(Q_{ij})}{\sum_{(k,m)\in
E}L_{km}g(Q_{km})}-\sqrt{D_{ij}}\ ,\quad(i,j)\in E,$ (5)

where $g(\cdot)$ is some arbitrary function describing the elasticity of
channels, $\alpha=V/\beta$ is a constant, and $t$ is the adimensional time.

Let $\mathcal{P}=\sum_{(i,j)\in E}Q_{ij}^{2}L_{ij}/D_{ij}$ be the dissipated
power of the H-P flow. To minimise the dissipated power while also keeping the
total volume of fluid constant, we consider the Lagrangian
$\mathcal{L}=\mathcal{P}-\lambda\left(V-\beta\sum_{(k,m)\in
E}L_{km}\sqrt{D_{km}}\right)$, where $\lambda$ is a Lagrange multiplier. Then
we have:

######  Theorem II.2 (Almeida and Dilão, (2022))

Assuming that the dissipated energy per unit time of the H-P flow is minimum,
then

$g(Q_{ij})=Q_{ij}^{2/3}.$ (6)

From these two theorems, the adaptive H-P flow has steady states obeying the
conditions

$D_{ij}^{*}=\alpha^{2}\frac{Q_{ij}^{4/3}}{\left(\sum_{(k,m)\in
E}L_{km}Q_{ij}^{2/3}\right)^{2}}\ \ \hbox{or}\ \ D_{ij}^{**}=0.$ (7)

The two theorem II.1 and II.2 lead to the generalised Murray’s law:

######  Corollary II.3 (Almeida and Dilão, (2022))

Assuming that the dissipated energy per unit time of the H-P flow is minimum,
then

$\sum_{j:(i,j)\in E,Q_{ij}>0}r_{ij}^{3}=\sum_{j:(i,j)\in
E,Q_{ij}<0}r_{ij}^{3},$ (8)

for every node $j$ of the graph $\mathcal{G}$.

Under the above construction, we can simulate out-of-equilibrium and steady
laminar fluxes on networks of channels with fixed sources and sinks. The
transverse adaptability of the network veins will determine a steady-state
graph that connects sources and sinks.

##  III Properties of the adaptive Hagen-Poiseuille flows

To study the properties of the adaptive H-P flow described by the above
theorems, we restrict our simulations to channel networks described by planar
graphs embedded in the two-dimensional Euclidean space. This approach easily
extended to $n$-dimensional Euclidian space through a graph embedding.

As biological networks evolve in time and adapt in response to the intensity
of flows, we start with graphs with small length edges resulting from a
Delaunay triangulation. The triangulation is generated from a distribution of
nodes in a square lattice with their positions randomly perturbed with
Gaussian noise, with a standard deviation $\sigma=0.5$. The lengths of each
edge $L_{ij}$ are calculated based on the coordinates of nodes $i$ and $j$.
Some of these nodes can be sources or sinks, and equation (1) is always
verified. The flow is initialised by randomly choosing positive arbitrary
conductivities $D_{ij}(t=0)>0$.

As equation (2) shows, the initial distribution of conductivities determines
the total fluid volume in the network of channels, thus affecting the value of
the parameter $\alpha=V/\beta$ in equation (5). To maintain the volume of
fluid constant for different initial conditions, we first choose the volume of
fluid $V$. As for any initial choice of conductivities,
$V_{0}=\beta\sum_{(i,j)\in E}L_{ij}\sqrt{D_{ij}(t_{0})}$, we renormalize the
initial conductivities to the new values

$D_{ij}^{\prime}(t=0)=\left(\frac{V}{V_{0}}\right)^{2}D_{ij}(t=0).$ (9)

Therefore, simulations for different initial conductivities and volumes can be
compared. This renormalisation condition is important as steady-state flows
depend on the distribution of initial conductivities.

To describe the temporal evolution of the adapted radii and fluid flows in all
the network channels, we start by calculating the initial pressures at all the
graph nodes. Then, channel conductivities are adapted according to equations
(5), in time increments according to the following rules: After each time step
$\Delta t$, new conductivities and new pressures and fluxes are calculated. To
find the pressure at each node $(i,j)$, we define $N$-dimensional vectors
$\mathbf{p}$ and $\mathbf{S}$, whose $i$th entries are $p_{i}$ and $S_{i}$,
respectively. The system of equations (3) is rewritten in the matrix form
$\mathbf{Wp}=\mathbf{S}$, where $\mathbf{W}$ is an $N\times N$ symmetric
matrix with entries $w_{ij}=\left(\sum_{(i,j\in
E)}C_{ik}\right)\delta_{ij}-C_{ij}$, and $C_{ij}=D_{ij}/L_{ij}$. We obtain the
pressures at each graph node after solving this linear system of equations. By
equation (4), the newly adapted fluxes $Q_{ij}$ are then obtained. To simplify
these calculations, without loss of generality, we assume that out-flow node
pressures are zero.

We used the Euler method to numerically integrate equations equations (5),
with $\Delta t=0.1$. Some parameters are kept constant for all the
simulations, namely $\beta=1$, $V=100$ and
$\sum_{j\in\\{\text{sources}\\}}S_{j}=1$.

These steps are repeated over time until a steady state of the channel
conductivities is reached. The steady state is reached when the condition

$\text{max}_{(i,j)\in E}\left|D_{ij}(n\Delta t)-D_{ij}\left((n-1)\Delta
t\right)\right|<D_{\text{thr}}$ (10)

is fulfilled. In all our simulations, we made the choice
$D_{\text{thr}}=5\times 10^{-6}$, where $n$ is the first integer for which the
inequality (10) is verified.

###  III.1 Length distributions of the steady states of the H-P flows

In all simulations, the initial distribution of nodes is inscribed in a square
of side $1$. The initial graph network $\mathcal{G}$ is shown in grey, and the
thickness of the black lines representing conducting channels is proportional
to the radius of the edges $\sim D_{ij}^{1/4}$. Source and sink nodes are
represented by circles and triangles, respectively. Different choices of
initial conductivities $\\{D_{ij}(t=0),(i,j)\in E\\}$ may lead to different
steady-state channel topologies connecting sources and sinks.

To quantify the different steady states of the adaptive H-P flows, the lengths
of the steady state graphs are $L=\sum_{(i,j)\in E^{\prime}}L_{ij}$, where
$E^{\prime}=\\{(i,j)\in E:D_{ij}>D_{\text{thr}}\\}$, where $E^{\prime}$ is the
set of conducting edges. For any $t\geq 0$ and finite $t$, the graph
$\mathcal{G}$ has only one connected component with conductivities partitioned
in two sets of low and high values.

Figures 1 and 2 show different steady-states obtained with one source and one
sink inscribed in a square of side $1$. In all simulations, we have obtained
different steady states generated by 300 different random initial conductivity
conditions. The statistics of the distribution of the lengths of the steady-
state graphs are also shown.

In figure 1, these simulations show that steady-state graphs are suboptimal
relative to the Euclidian distance between the source and the sink. Decreasing
the mean lengths of the sides of Delaunay triangles (figure 1), the
distribution of lengths approaches the minimum distance between the source and
the sink. However, in the case of figure 2, this tendency in not evident.

Figure 1: Steady states of minimum and maximal length of the adaptive H-P
flows for one source and one sink, in two graph lattices of sizes (a)
$25\times 25$, and (b) $40\times 40$. The steady-state graphs’ length
distribution is also shown in each case. In a) The mean of the length
distribution is $\overline{L}=1.312$, and the standard deviation is
$\sigma=0.032$. In b), $\overline{L}=1.245$, and $\sigma=0.017$. The Euclidian
distance between source and sink is $L_{\text{lin}}=4\sqrt{2}/5\approx 1.131$.
Figure 2: Steady states of minimum and maximal length of the adaptive H-P
flows for five sources and four sinks, in two graph lattices of sizes (a)
$25\times 25$, and (b) $40\times 40$. The steady-state graphs’ length
distribution is also shown in each case. In a) The mean of the length
distribution is $\overline{L}=3.39$, and the standard deviation is
$\sigma=0.28$. In b), $\overline{L}=3.49$, and $\sigma=0.33$. Figure 3: Two
possible steady-state solutions of the adaptive H-P flows for two sources and
two sinks in a $25\times 25$ graph lattice. The trees were generated using the
same initial conditions. In (a), $S_{\text{sources}}=-S_{\text{sinks}}=0.5$.
In (b), $S_{\text{left source}}=-S_{\text{right sink}}=0.3$, and
$S_{\text{right source}}=-S_{\text{left sink}}=0.7$.

In figure 3, we show two different steady-states of the adaptive H-P flows for
the same initial conductivities. The source and sink nodes were the same in
both cases, but the intensities of sources and sinks differed. In this case,
one of the steady states is a connected graph tree, and the other is a
disconnected graph tree.

One of the remarkable features of the H-P flow is that all the simulations
reached a steady state graph topology, even if we introduce stochasticity in
the intensity of sources and sinksAlmeida, (2021); Valente, (2023). In all the
cases, the steady-state flows obey the conditions (7).

###  III.2 Validation of Murray’s law

The generalised Murray’s law (8) is a direct consequence of the set of
equations used to determine the steady state of the adaptive H-P flow, but it
is unknown whether or not it is verified during out-of-equilibrium states.

To verify the transient behaviour of the H-P flow, the source and sink setup
of figure 2 were used, with initial conductivities $D_{ij}(t=0)=1$, for every
pair $(i,j)\in E$. We have chosen two branching points of the steady state,
and we have analysed the temporal evolution of the quantity

$\left|R_{\text{out}}(t)-R_{\text{in}}(t)\right|:=\left|\sum_{\begin{subarray}{c}j:(i,j)\in
E\\\ Q_{ij}<0\end{subarray}}r_{ij}^{3}(t)-\sum_{\begin{subarray}{c}j:(i,j)\in
E\\\ Q_{ij}>0\end{subarray}}r_{ij}^{3}(t)\right|$

The Murray’s law is verified if the quantity
$\left|R_{\text{out}}(t)-R_{\text{in}}(t)\right|$ is zero.

In figure 4, we show the localisation of two-channel branching nodes of a
steady state H-P flow and the time evolution of
$\left|R_{\text{out}}-R_{\text{in}}\right|$ at the two branching nodes. This
shows that the generalised Murray’s law (8) is only verified for steady-state
flows. Experiments with Physarum suggested the validity of Murray’s lawAkita
et al., (2017).

Figure 4: Converge to Murray’s law at two different branching nodes. a) Square
boxes show branching nodes at a steady state of the H-P flow for the source
and sink setup of figure 2. In b) and c), it is shown the variation of
$\left|R_{\text{out}}(t)-R_{\text{in}}(t)\right|/\left|R_{\text{out}}(0)-R_{\text{in}}(0)\right|$
over time for leftmost and rightmost branching nodes.

###  III.3 Properties of simple networks

The formalisms of adaptive H-P flows can be applied to simple networks with
fixed sources and sinks, and the fluid flows through all the channels.

####  III.3.1 One channel networks

We apply the conductivity adaptation algorithm to a single channel of length
$L_{12}$, with one source and sink. Assuming that the source is located at the
node number $1$ and the sink at node number $2$, the adaptation equation is

$\frac{d}{dt}\sqrt{D_{12}}=\alpha\frac{1}{L_{12}}-\sqrt{D_{12}}=\sqrt{D_{12}(0)}\frac{1}{L_{12}}-\sqrt{D_{12}},$
(11)

where we have used the relations $\alpha=V/\beta$, $\beta=\sqrt{8\pi\eta}$,
$V=\beta L_{12}\sqrt{D_{12}(0)}$ and $\eta$ is the dynamic viscosity of the
fluid. The fluxes at the nodes $1$ and $2$ are $Q_{12}=I_{1}>0$ and
$Q_{21}=-I_{1}$. Equation (11) is linear and its general solution is

$D_{12}(t)=\left(\frac{\alpha}{L_{12}}-\left(\frac{\alpha}{L_{12}}-\sqrt{D_{12}(0)}\right)e^{-t}\right)^{2}=D_{12}(0).$
(12)

The steady state is
$D_{12}^{*}=\lim_{t\to\infty}D_{12}(t)=\frac{\alpha^{2}}{L_{12}^{2}}=D_{12}(0)$.

As $D_{12}=\pi r_{12}^{4}/(8\eta)$, and the flux is
$Q_{12}=D_{12}(p_{1}-p_{2})/L_{12}$, at the steady state we have

$Q_{12}^{*}=D_{12}^{*}\frac{1}{L_{12}}(p_{1}(\infty)-p_{2})=\frac{\alpha^{2}}{L_{12}^{3}}p_{1}(\infty)=I_{1},$
(13)

where $p_{1}(\infty)$ is determined by the flux at the source $Q_{12}=I_{1}$,
and we made the choice $p_{2}=0$.

The steady-state solution $Q_{12}^{*}$ in (13) shows that, for fixed volume
$V$ ($\alpha$ constant), the flux at the source increases linearly with the
pressure. The radius of the channel is constant for every $t\geq 0$ and is
determined by $D_{12}(0)$.

In figure 5, we show the relationship between the pressure and channel length
and source fluxes.

Figure 5: Relationship between the pressure and channel length and source
fluxes in a one channel H-P flows, with constant volume. The parameters are
$\alpha=1$, and $I_{1}=1$ in a), and $L_{12}=1$ in b).

####  III.3.2 Two channels networks

In the case of two channels with two sinks connected to a single source with
an input flux $I_{1}$, the adaptation equations for the Hagen-Poiseuille flow
are

$\begin{array}[]{l}\displaystyle\frac{d}{dt}\sqrt{D_{12}}=\alpha\frac{Q_{12}^{2/3}}{L_{12}Q_{12}^{2/3}+L_{13}Q_{13}^{2/3}}-\sqrt{D_{12}}\\\
\displaystyle\frac{d}{dt}\sqrt{D_{13}}=\alpha\frac{Q_{13}^{2/3}}{L_{12}Q_{12}^{2/3}+L_{13}Q_{13}^{2/3}}-\sqrt{D_{13}}\end{array}$
(14)

with node fluxes

$Q_{12}+Q_{13}=I_{1}>0,\ Q_{21}=I_{2}<0,\ \hbox{and}\ Q_{31}=I_{3}<0$ (15)

and the conservation law $I_{1}+I_{2}+I_{3}=0$.

Introducing the relations (15), into (14), these equations simplify and become
linear. Therefore, the steady-state conductivities are

$\begin{array}[]{l}\displaystyle
D_{12}^{*}=\alpha^{2}\frac{|I_{2}|^{4/3}}{(L_{12}|I_{2}|^{2/3}+L_{13}|I_{3}|^{2/3})^{2}}\\\
\displaystyle
D_{13}^{*}=\alpha^{2}\frac{|I_{3}|^{4/3}}{(L_{12}|I_{2}|^{2/3}+L_{13}|I_{3}|^{2/3})^{2}},\end{array}$
(16)

where $\alpha=V/\beta=(L_{12}\sqrt{D_{12}(0)}+L_{13}\sqrt{D_{13}(0)})$.
Therefore, as $D_{12}^{*}\leq\alpha^{2}/L_{12}^{2}$ and
$D_{13}^{*}\leq\alpha^{2}/L_{13}^{2}$, we have

$\left(D_{12}^{*}\frac{1}{L_{12}}+D_{13}^{*}\frac{1}{L_{13}}\right)p_{1}(\infty)=I_{1}.$
(17)

Therefore, for fixed volume $V$, at the steady state, the pressure at the
source is linear as a function of the flux at the source. However, during time
evolution, the flow will adapt for the conductivity values $D_{12}^{*}$ and
$D_{13}^{*}$.

In figure 6, we show the time evolution of the adapted conductivities and the
pressure at the source for the two channels with fixed lengths and constant
volume. Equation (17) shows that the steady-state pressure at the source
depends on the lengths of the channels.

Figure 6: Channel conductivities and pressure at the source in two channels
H-P flows, with constant volume. The parameters of the simulation are
$\alpha=1$, $I_{1}=1$, $I_{2}=-0.3$, $I_{3}=-0.7$, $L_{12}=0.4$, $L_{13}=0.6$,
$D_{12}(0)=0.8$, $D_{13}(0)=1.14572$.

In figure 7, we have calculated the steady-state pressures for different
channel lengths and several fluxes at the source. Data shows the increase in
pressure at the source does not increase the flux at the source. This
saturation phenomenon has been observed experimentally for contractile
veinsRubinow and Keller, (1972).

Figure 7: Fluxes and pressures at the source in two channels H-P flows. The
parameters of the simulations are $\alpha=1$, $I_{1}\in[1,3]$,
$I_{2}=-I_{1}/2$, $I_{3}=-I_{1}/2$, $L_{12}\in[0.5,0.8]$, $L_{12}+L_{13}=1$,
$D_{12}(0)=0.8$, $D_{13}(0)=1.14572$.

##  IV Comparison between steady-state H-P flow patterns with Physarum
patterns

Physarum polycephalum’s growth patterns are one of the biological phenomena
used to develop the class of models described by adaptive H-P flows on graphs.

To test the development of channels in a Physarum-like configuration of
sources and sinks, we consider one source in the middle of a square grid
surrounded by $20$ sinks evenly placed along a circumference of radius $0.4$.
The fluxes are $S_{\text{source}}=1$ and $S_{j}=-1/20$, for every
$j\in\\{\text{sinks}\\}$. All sources and sinks are active at all times. This
configuration describes the body of a Physarum organism with a food source in
the middle. The sinks represent the reproducing Physarum peripheric nuclei
undergoing mitosis and therefore responsible for growth.

The shortest and longest channel graphs at steady state obtained after
adaptation are shown in figure 8a). In figure 8b), the steady-state graph
length distributions and the steady-state graph length $L$ as a function of
the pressure at the source node $p_{\text{source}}$ are shown. In these
simulations, $p_{j}=0$, for every $j\in\\{\text{sinks}\\}$. The patterns
obtained with the adaptive H-P flow are similar to the ones found in
PhysarumAlmeida and Dilão, (2022); Jones, (2013).

Figure 8: a) Steady-state configurations with minimum and maximum lengths for
a Physarum-like configuration of sources and sinks, for $V=100$,
$S_{\text{source}}=1$, $S_{j}=-1/20$, for every $j\in\\{\text{sinks}\\}$,
obtained from 300 runs of the simulation with random starting conductivities.
The simulations were done in a $25\times 25$ graph lattice. b) The steady-
state graph length distributions and the steady-state graph length $L$ as a
function of the pressure at the source node $p_{\text{source}}$. The mean of
the length distribution is $\overline{L}=6.50$ with the standard deviation
$\sigma=0.22$.

We now test Physarum-like configuration, but we consider that sinks are not
simultaneously active (asynchronous adaptation). In particular, we have
considered that at each iteration, a random number of sinks
$N_{\text{sinks}}$, between 1 and 20, is chosen to become active. The other
sinks are inactive. In this case, $S_{j}=-1/N_{\text{sinks}}$, for every
$j\in\\{\text{active sinks}\\}$, otherwise, $S_{j}=0$, for every
$j\in\\{\text{inactive sinks}\\}$. The source is always active with
$S_{\text{source}}=1$. This configuration also mimic the outward radial
distribution of channels in Physarum.

Figure 9 shows six final states after $3\times 10^{4}$ iterations for this
configuration and their respective lengths $L$, obtained using different
random initial conditions. The active sinks at the end of the simulation are
shown as grey triangles, while inactive sinks are shown as white triangles.

The states of figure 9 are not formally steady states as the active sinks
change every time step, the values of the conductivities $D_{ij}$ never reach
a fixed value, changing over time. However, the graph’s topology does reach a
constant form, and the length of the graph has a constant value, where $L$
decreases over time to a minimum stable length. For this reason, the final
channel configurations are called topological steady-states. However, during
the asynchronous adaptation, the fluxes vary over time, but the topology of
channels remains constant.

Figure 9: Six topological steady-state channels resulting from the
asynchronous adaptation of Physarum-like configuration of sources and sinks,
for $V=100$ and $S_{\text{source}}=1$, in a $25\times 25$ graph lattice. The
number of active sinks $N_{\text{sinks}}$ is chosen randomly in the range
$1-20$. So, $S_{j}=-1/N_{\text{sinks}}$, for every $j\in\\{\text{active
sinks}\\}$, otherwise, $S_{j}=0$, for every $j\in\\{\text{inactive sinks}\\}$.
Initial conductivities were chosen randomly. The topological steady-states
stabilised after $3\times 10^{4}$ iterations.

Comparing the simulations of figures 8 and 9, in the latter there are fewer
channels emanating from the source, explaining the growth pattern selection
during Physarum growthJones, (2013).

One of the features that are not explained in this approach is the development
of channel loops in the stabilised regions of the Physarum’s channel growth.

##  V Shuttle streaming

To analyse shuttle streaming in Physarum, we started with the asynchronous
algorithm of the previous section. The source in the middle is always active
with $S_{\text{source}}=1$ and, out of all the $20$ possible sinks circling
the source, at each iteration, $N$ sinks are randomly chosen to be active
sinks. We have run the simulations for $10^{4}$ iterations.

Figure 10: Topological steady-state configuration for the asynchronous model
of figure 9, obtained after $10^{4}$ iterations of the adaptation algorithm.
To analyse the shuttle stream we have chosen $7$ different edges and we have
indexed them with the capital letters A-G.

After reaching the topological steady-state, $7$ different edges were chosen.
In figure 10, we show the chosen edges surrounded by squared boxes. Then we
have calculated the fluxes $Q_{ij}$ as a function of time for these edges.

In figure 11, we show $Q_{ij}(t)$ for edges E and F. In edge E, we observe the
inversion of the direction of flux – shuttle streaming – during all the
iterations. In edge F (and B, D, G), we observe a transient shuttle streaming
effect, and after some iteration threshold, the sign of $Q_{ij}(t)$ becomes
constant. Edges A and C don’t show shuttle streaming for all the $10^{4}$
iterations.

Due to the nature of the asynchronous adaptation, shuttle streaming appears in
out-of-equilibrium regions of the channel network domain.

Figure 11: Time evolutions of the fluxes $Q_{ij}$ during the adaptation of H-P
flow, for the edges E and F shown in figure 10. These simulations prove the
presence of shuttle streaming during asynchronous adaptation.

##  VI Peristalsis

the adaptive H-P flow on elastic channel networks contracts and relaxes the
vein walls of the network over time.

In Physarum growth, some of the elastic channels contract and relax over time.
This phenomenon is referred to as peristalsis. To evaluate the occurrence of
peristalsis, we evaluate the channel radius $r_{ij}(t)=(8\eta
D_{ij}(t)/\pi)^{1/4}$ of the edges shown in figure 10 over time, during the
first $10^{3}$ iterations of the H-P adaptation algorithm. In figure 12, we
show the results for nodes A and E of figure 10. Therefore, the peristalsis
phenomenon is observed during the transient process of the asynchronous
adaptation in the H-P flow.

Figure 12: Time evolutions of the channel radius $r_{ij}$ during the
adaptation of H-P flow, at edges A and E of figure 10. These simulations prove
the presence of peristalsis during asynchronous adaptation.

Peristalsis and shuttle streaming are only shown in simulations with
asynchronous adaptation.

##  VII Conclusions

We have derived the main properties of adaptive Hagen-Poiseuille flows in
channel networks, simulating the transverse elasticity of vein networks in
organisms.

We have shown that the model tested produces results consistent with physical
observations, namely the tendency of the adaptation to select channels close
to the shortest path connecting sources and sinks. The lengths of the steady-
state adapted veins connecting sources and sinks are suboptimal, at least for
finite triangular geometries, refuting the hypothesis of the possibility of
solving travelling salesman or Steiner geometry problems in graphs.

Additionally, Murray’s law, experimentally observed in several biological
networks, is verified only for steady-state H-P flow regimes.

The adaptive H-P flow description applied to simple one-channel and two-
channel networks predicts the behaviour of the pressure in contractile veins
as observed in several organisms. This process is associated with conductivity
saturation, which results in a threshold in the fluid fluxes at the source.
Increasing the pressure at the source the flux remains constant.

We have compared the Physarum’s radial growth pattern of veins for fixed
sources and sinks, and sinks activated randomly in time – asynchronous
adaptation. In the latter, the pattern characteristics of Physarum networks
were evident, showing fewer thicker veins closer to food sources, branching
out into thinner veins as the periphery approached. A new topological steady
state has been identified, supporting out-of-equilibrium fluxes. In these
topological steady states, shuttle streaming and peristalsis are observed.
These out-of-equilibrium effects are also observed in Physarum. These
simulations show that these out-of-equilibrium effects observed in Physarum
may have a hydrodynamic origin.

## References

  * Takahashi, (2014) T. Takahashi, Microcirculation in Fractal Branching Networks, Springer, Tokyo, 2014.
  * Rubinow, (1973) S. I. Rubinow, Mathematical Problems in the Biological Sciences, SIAM, Philadelphia, 1973.
  * Chang et al., (2019) S. Chang and M. Roper, Microvascular networks with uniform flow, J. Theo. Bio. 462 (2019) 48-64.
  * Katifori et al., (2010) E. Katifori, G. J. Szollosi and M. O. Magnasco, Damage and fluctuations induce loops in optimal transport networks, Phys. Rev. Lett. 104 (2010) 048704.
  * Baumgarten et al., (2010) W. Baumgarten, T. Ueda and M. J. B. Hauser, Plasmodial vein networks of the slime mold Physarum polycephalum form regular graphs, Phys. Rev. E 82 (2010) 046113.
  * Oettmeier et al., (2017) C. Oettmeier, K. Brix, H.-G. Döbereiner, Physarum polycephalum — a new take on a classic model system, J. of Phys. D: Appl. Phys. 50 (2017) 413001.
  * Tero et al., (2006) A. Tero, R. Kobayashi and T. Nakagaki, Physarum solver: A biologically inspired method of road-network navigation. Physica A 363 (2006) 115-119.
  * Hu and Cai, (2013) D. Hu and D. Cai, Adaptation and optimization of biological transport networks, Phys. Rev. Lett. 111 (2013) 138701.
  * Bohn and Magnasco, (2007) S. Bohn and M. O. Magnasco, Structure, scaling, and phase transition in the optimal transport network, Phys. Rev. Lett. 98 (2007) 088702.
  * Zhu et al., (2018) L. Zhu, S.-J. Kim, M. Hara and M. Aono, Remarkable problem-solving ability of unicellular amoeboid organism and its mechanism, R. Soc. open sci. 5 (2018) 180396.
  * Tero et al., (2010) A. Tero, T. Nakagaki, K. Toyabe, K. Yumiki, and R. Kobayashi, A method inspired by Physarum for solving the Steiner problem, Int. J. of Unconventional Computing 6 109-123 (2010).
  * Dove et al., (1985) W. F. Dove, J. Dee, S. Hatano, F. B. Aaugli and K.-E. Wohlfarth-BottermannThe Molecular Biology of Physarum polycephalum, NATO ASI series. Series A, Life Sciences, v. 106, Plenum Press, 1985.
  * Wohlfarth, (1979) K. E. Wohlfarth-Bottermann, Oscillatory contraction activity in Physarum, J. Exp. Biol. 81(1) (1979) 15-32.
  * Nakagaki et al., (2000) T. Nakagaki, H. Yamada, and Á. Tóth, Maze-solving by an amoeboid organism, Nature, 407 (2000) 470.
  * Tero et al., (2010) A. Tero, S. Takagi, T. Saigusa, K. Ito, D. P. Bebber, M. D. Fricker, K. Yumiki, R. Kobayashi, and T. Nakagaki, Rules for biologically inspired adaptive network design, Science 327 (2010) 439-442.
  * Almeida, (2021) R. Almeida, Formation and Optimisation of Vein Networks in Physarum, Master thesis, University of Lisbon, IST, 2021, arXiv:2305.12244v1.
  * Almeida and Dilão, (2022) R. Almeida and R. Dilão, Adaptive Hagen-Poiseuille flows on graphs, Physica D 426 (2022) 133322.
  * Valente, (2023) A. F. Valente, Applications to Biological Networks of Adaptive Hagen-Poiseuille Flow on Graphs, Master thesis, University of Lisbon, IST, 2023, arXiv:2309.10577.
  * Murray, (1926) C. D. Murray, The Physiological Principle of Minimum Work, PNAS 12 (1926) 207-214.
  * Sherman, (1981) T. F. Sherman, On connecting large vessels to small. The meaning of Murray’s law. The Journal of General Physiology, 78(4) (1981) 431-453.
  * Rubinow and Keller, (1972) S. I. Rubinow and J. B. Keller, Flow of a viscous fluid through an elastic tube with applications to blood flow, J. Theor. Biol. 35 (1972) 299-313.
  * Jones, (2013) J. Jones and A. Adamatzky, Slime mould inspired generalised Voronoi diagrams with repulsive fields, Int. J. Bifur. Chaos, arXiv:1503.06973v1, 2015.
  * Akita et al., (2017) D. Akita, I. Kunita, M. D. Fricker, S. Kuroda, K. Sato and T. Nakagaki, Experimental models for Murray´s law, J. Phys. D: Appl. Phys. 50 (2017) 024001.


# The Mathematics of Chlorophyll Photoexcitation

Archit Chaturvedi
University of California  Berkeley
<EMAIL_ADDRESS>
(March 2024)

###### Abstract

Photosynthesis is a fundamental process for plants to produce energy and
survive. It is a well known fact that the light reactions in photosynthesis
are a significant part of the overall process, and are carried out by
chlorophyll molecules. This paper outlines the quantum chemical and
mathematical description concerning electronic state transitions that take
place in chlorophyll molecules during photosynthesis. In particular, the
electronic transition of chlorophylls is assumed as a diatomic rovibrational
$0\longrightarrow 0$ transition, and is mathematically described accordingly.
Then, the correlation of the Franck-Condon Principle in accordance to the
photoexcitation of chlorophylls is delineated. Finally, existing research
concerning the concepts laid out in this paper is discussed, and future
implications of quantum chemical theory concerning chlorophyll photoexcitation
is discussed.

##  1 Introduction

Photosynthesis is a fundamental process in plants. It is responsible for
producing chemical energy in plants to drive their survival. Photosynthesis is
divided into two parts: 1. Light Reactions and 2. Calvin Cycle ([6]). In the
context of this paper, light reactions are of focus. Chlorophylls are the most
significant molecule when discussing light reactions in photosynthesis.
Chlorophylls are located in the thylakoid membranes of a plant’s chloroplast.
These chlorophylls are photoexcited through absorbing photons, and it is this
photoexcitation of chlorophylls that ultimately drives photosynthesis. Of the
chlorophyll molecules, chlorophyll a is the most significant molecule, as it
is the key pigment in a plant’s light reactions. There are other types of
chlorophylls, such as chlorophyll b, which are discussed towards the end of
the paper. Ultimately, theoretical formulations in quantum chemistry can be
used to mathematically describe the photoexcitation of chlorophyll molecules.
This paper aims to elucidate the mathematics behind such electronics
transitions that take place in a chlorophyll molecule. One model is to assume
the electronic transitions of chlorophylls as diatomic rovibrational
transitions ([1, 2, 3]). This can give the total energy of a chloroplast’s
electronic state, and can also describe the transition energy required for the
photoexcitation of chloroplasts. Furthermore, the Franck-Condon principle is
significant when discussing electronic transitions ([1, 4, 5]). It describes
electronic transitions under a fixed nuclear position, and can give the
probability that an electronic transition occurs based on nuclear overlap.

##  2 Electronic Transitions in Chlorophylls

We can view the electronic transition of a chlorophyll molecule as a diatomic
rovibrational transition; that is, a transition with both rotations and
vibrations. Approximating based on an anharmonic oscillator model with a
nonrigid rotator, the energy of the chlorophyll’s ground state can be written
as:

$\displaystyle E_{ground}=G(g)+F(J)$ (1)
$\displaystyle=\nu_{g}\left(g+\frac{1}{2}\right)-\chi_{g}\nu_{g}\left(g+\frac{1}{2}\right)^{2}+B_{g}J_{g}(J_{g}+1)-D_{g}J_{g}^{2}(J_{g}+1)^{2}$
(2)

In this equation, $\nu_{g}$ denotes the vibration frequency at ground state,
$\chi_{g}$ denotes the anharmonic constant at ground state, $g$ denotes the
vibrational quantum number at ground state, $B_{g}$ denotes the rotational
constant at ground state, $J_{g}$ denotes the rotational quantum number of the
ground state, and $D_{g}$ denotes the centrifugal distortion constant at the
ground state. Ultimately, $G(g)$ marks the vibrational energy of the molecule,
and $F(J)$ marks the rotational energy of the molecule. Denoting the total
energy at the excited state requires simply adding the change in electronic
transition energy, denoted as $\nu_{el}$:

$E_{excited}=\nu_{el}+\nu_{g}\left(g+\frac{1}{2}\right)-\chi_{g}\nu_{g}\left(g+\frac{1}{2}\right)^{2}+B_{g}J_{g}(J_{g}+1)-D_{g}J_{g}^{2}(J_{g}+1)^{2}$
(3)

The eigenstate to eigenstate transition between the ground state and excited
state can be written as:

$E_{transition}=E_{excited}-E_{ground}$ (4)

Assuming a $0\longrightarrow 0$ transition between the ground state and
excited state would be the most reasonable in the case of chlorophylls, as
this is the most common and significant type of electronic transition.
Ultimately, in such a case, there would be no vibrational change, and so,
letting $\nu_{e}$ denote the vibrational frequency of the excited state and
letting $\chi_{e}$ denote the anharmonic constant of the excited state, the
total transition frequency can be written as:

$E_{transition}=T_{el}+\left(\frac{1}{2}\nu_{g}-\frac{1}{4}\chi_{g}\nu_{g}\right)-\left(\frac{1}{2}\nu_{e}-\frac{1}{4}\chi_{e}\nu_{e}\right)$
(5)

In the equation above, it should be noted that $E_{transition}$ represents the
transition frequency, or the overall difference in energy between the two
electronic and vibration states. And $T_{el}$ represents the electronic
transition energy. The rest of the equation marks the vibrational contribution
to the transition as a result of anharmonicity. Ultimately, $E_{transition}$
denotes the smallest possible energy to observe an electronic transition
between the ground state and excited state in a chlorophyll molecule. And the
electronic transition energy can be found by subtracting the vibrational
component of equation 5 from $E_{transition}$.

##  3 Franck-Condon Principle and Chlorophylls

It would be rather vacuous to discuss electronic transitions in chlorophylls
without including the Franck-Condon Principle. As described, light reactions
within a chlorophyll molecule involve the absorption of a photon, which is
ultimately what causes the electronic transition. Franck-Condon Principle
describes the intensities of vibronic transitions by assuming a fixed nuclear
position during such transitions. The probability of a transition occurring
from the ground state to the excited state can be given as:

$\left\langle\psi_{p,e}|\mu|\psi_{p,g}\right\rangle=\left\langle\psi_{nuc,e}^{*}|\psi_{nuc,g}\right\rangle\left\langle\psi_{el,e}^{*}|\mu|\psi_{el,g}\right\rangle$
(6)

The first operator represents the nuclear overlap of the transition, and the
second operator represents the electronic component of the transition.
Ultimately, if the nuclear overlap is 0, then the electronic transition will
not occur, despite the electronic component of the equation. A harmonic
oscillator model is perhaps the simplest of models used to denote nuclear
overlap in electronic transitions based on the Franck-Condon principle. In
such a model, the wavefunction of the ground electronic state can be written
as:

$\left|\psi_{g}(Q)\right\rangle=\left|\left(\frac{\alpha}{\pi}\right)^{\frac{1}{4}}\textbf{exp}\left(\frac{-\alpha(Q-Q_{g})^{2}}{2}\right)\right\rangle$
(7)

In the equation above, $\alpha$ is a parameter concerning the vibrational
motion of the chlorophyll molecule:

$\alpha=\frac{\sqrt{mk}}{\hbar},$ (8)

with $m$ being the reduced mass of the system and $k$ denoting the spring
constant. Furthermore, $Q$ denotes the position of the nuclei in the molecule,
and $Q_{g}$ denotes the equilibrium bond length of the ground electronic
state. Based on $\left|\psi_{g}(Q)\right\rangle$, the wavefunction of the
excited electronic state can also be written as:

$\left|\psi_{e}(Q)\right\rangle=\left|\left(\frac{\alpha}{\pi}\right)^{\frac{1}{4}}\textbf{exp}\left(\frac{-\alpha(Q-Q_{e})^{2}}{2}\right)\right\rangle$
(9)

From this, the nuclear overlap integral can be written as:

$S_{ge}=\sqrt{\frac{\alpha}{\pi}}\int_{-\infty}^{\infty}\textbf{exp}\left(\frac{-\alpha\left(2Q^{2}-QQ_{g}-2QQ_{e}+Q_{g}^{2}+Q_{e}^{2}\right)}{2}\right)dQ$
(10)

Ultimately, since:

$Q_{g}^{2}+Q_{e}^{2}=\frac{1}{2}\left[(Q_{g}+Q_{e})^{2}+(Q_{g}-Q_{e})^{2}\right],$
(11)

the integral turns out to be a Gaussian integral, and so therefore, $S_{ge}$
can finally be expressed as:

$S_{ge}=\textbf{exp}\left(\frac{-\alpha(Q_{g}-Q_{e})^{2}}{4}\right)$ (12)

This represents the value of the nuclear overlap of an electronic transition
in a chlorophyll molecule and can be used to calculate the probability of
chlorophyll photoexcitation given the electronic components of the chlorophyll
molecule.

##  4 Existing Research and Future Implications

There has been quite some research done on the quantum chemical properties of
chlorophylls. The vibrational properties of chlorophylls have been explored by
numerous studies. For example, source [7] calculated the vibrational
frequencies for chlorophyll a-4, chlorophyll a4+, chlorophyll a5, and
chlorophyll a5+. Furthermore, source [8] calculated the vibrational
frequencies associated witht he excited state of chlorophyll f and chlorophyll
a based on the Franck-Condon Principle. Sources [9] and [10] explore further
chemical properties of chlorophyll a and bacteriochlorophyll a in accordance
to the Franck-Condon principle. However, while much research has been done
concerning the electronic transitions and quantum chemical properties of
chlorophyll molecules, further research is needed to be done to fully
understand the overall picture of photoexcitation in chlorophylls during
photosynthesis. For example, using the models discussed in this paper, while
vibrational frequencies for chlorophylls have been calculated, the anharmonic
constant of chlorophyll molecules still needs to be researched upon further.
Also, the rotational component of electronic transitions in chlorophyll
molecules needs further experimental research. Research concerning the
equilibrium bond length of the ground/excited states of different types of
chlorophyll molecules, as well as the orientation of atomic nuclei in
different chlorophyll molecules may also be useful when discussing the
photoexcitation of chlorophyll molecules. Finally, since there are so many
different types of chlorophyll molecules out there, research concerning the
electronic transitions in such chlorophylls, as well as the Franck-Condon
components of these chlorophyll molecules, would open up exciting discoveries
and opportunities concerning the photoexcitation of chlorophylls during light
reactions.

## References

  * [1] Mark Tuckerman, Quantum Chemistry, 2024
  * [2] Tully, John C. ”Molecular dynamics with electronic transitions.” The Journal of Chemical Physics 93.2 (1990): 1061-1071.
  * [3] Kasha, Michael. ”Characterization of electronic transitions in complex molecules.” Discussions of the Faraday society 9 (1950): 14-19.
  * [4] Coolidge, Albert Sprague, Hubert M. James, and Richard D. Present. ”A study of the Franck‐Condon principle.” The Journal of Chemical Physics 4.3 (1936): 193-211.
  * [5] Condon, Edward U. ”The Franck-Condon principle and related topics.” American journal of physics 15.5 (1947): 365-374.
  * [6] Reece, Jane B., Lisa A. Urry, and Michael L. Cain. Campbell biology. Pearson, 2017.
  * [7] Parameswaran, Sreeja, Ruili Wang, and Gary Hastings. ”Calculation of the vibrational properties of chlorophyll a in solution.” The Journal of Physical Chemistry B 112.44 (2008): 14056-14062.
  * [8] Zamzam, Noura, and Jasper J. van Thor. ”Excited state frequencies of chlorophyll f and chlorophyll a and evaluation of displacement through franck-condon progression calculations.” Molecules 24.7 (2019): 1326.
  * [9] Zazubovich, V., I. Tibe, and G. J. Small. ”Bacteriochlorophyll a Franck- Condon Factors for the S0→ S1 (Q y) Transition.” The Journal of Physical Chemistry B 105.49 (2001): 12410-12417.
  * [10] Pieper, J., J. Voigt, and G. J. Small. ”Chlorophyll a Franck- Condon factors and excitation energy transfer.” The Journal of Physical Chemistry B 103.13 (1999): 2319-2322.