# Instabilities, Anomalous Heating and Stochastic Acceleration of Electrons in
Colliding Plasmas

M.A. Malkov1 and V.I. Sotnikov2 1CASS and Department of Physics, University of
California, San Diego, La Jolla, CA 92093
2Air Force Research Laboratory, WPAFB, OH

###### Abstract

The collision of two expanding plasma clouds is investigated, emphasizing
instabilities and electron energization in the plasma mixing layer. This work
is directly relevant to laboratory experiments with explosively-created laser
or z-pinch plasmas but may also elucidate naturally occurring plasma
collisions in astrophysical or space physics contexts. In the previous
publications [MalkovSotnikov2018PhPl, Sotnikov2020], we have studied,
analytically and numerically, the flow emerging from interpenetrating coronas
launched by two parallel wires vaporized in a vacuum chamber. The main foci of
the studies have been on the general flow pattern and lower-hybrid and thin-
shell instabilities that under certain conditions develop in the collision
layer. The present paper centers around the initial phase of the
interpenetration of the two plasmas. A two-stream ion-ion instability,
efficient electron heating, and stochastic acceleration dominate plasma mixing
at this phase. Both the adiabatic (reversible) electron heating and stochastic
(irreversible) acceleration and heating mechanisms, powered by unstably driven
electric fields, are considered. The irreversibility results from a
combination of electron runaway acceleration in the wave electric fields and
pitch-angle scattering on ions and neutrals.

##  I Introduction

Macroscopic plasma motions often harbor microscopic phenomena that are
challenging to understand or even observe. The hi-speed collision of two
plasma clouds is an excellent example of how microinstabilities control the
macroscopic flow that results from the collision. The collision of ordinary
gases results in a contact surface and a pair of shocks propagating into each
gas. No mixing other than a relatively slow diffusion across the contact
surface is expected, at least in a perfectly symmetric stable flow. Plasma
collisions are different in several ways. First, colliding plasma clouds may
go through each other without much effect. This regime is expected when the
relative plasma motion is collisionless, owing to its high velocity, whereas
the individual plasmas may or may not be. At the same time, instabilities
associated with the relative motion (two-stream instability) may have no time
to develop during the plasma interaction or be suppressed by the internal
pressure of the streams [VVS_1961SvPhUsp]. The difference in collisionality
regimes may occur when the relative velocity between the clouds is much larger
than the thermal velocity in each of them, owing to the steep growth of the
collision mean free path with the relative particle velocity, $\lambda\propto
v^{4}$. Sustainable counterstreaming ion flows have also been observed in
numerical simulations with a contact discontinuity-type of initial state
[Moreno_2019]. The interpenetration of adjacent plasmas of different densities
results in a phase space structure where mixed populations intersperse
residual areas of counterstreaming protons. So the phase-space develops
density holes.

This paper will consider microscopic instabilities and electron heating in
interpenetrating plasmas using unmagnetized plasma approximation near a
symmetry plane of collision. The collision regime of magnetized electrons and
unmagnetized ions has been addressed in [MalkovSotnikov2018PhPl] along with
the macroscopic two-dimensional flow emerging from the collision. Macroscopic
instabilities of this flow have been studied numerically in [Sotnikov2020] in
hydrodynamic approximation. When focusing on the microscopic instabilities and
associated heating processes, further progress can be made by considering
plasma interpenetration along the centerline in one dimension, which we pursue
in this paper.

Plasma collisions are best studied in experiments with laser and z-pinch
plasmas. The present work is motivated by laboratory experiments with two
ohmically-exploding parallel wires. Being vaporized, they launch plasma
coronas toward each other. The interpenetration and mixing of the hot coronas
are followed by a collision and mixing of denser phases of the melted wire
material. They are also much colder than the coronas, weakly ionized, and move
at an order of magnitude lower speed (typical at 3-5 km/s vs. 50-100 km/s,
[Sarkisov2005]).

Nevertheless, this phase is characterized by protracted light emission from
the collision layer, which is considerably less pronounced and shorter in time
in the case of a single wire explosion. The difference speaks to the
importance of the corona interpenetration stage of the collision process.
Otherwise, two closely spaced parallel wires would exhibit the same flow
morphology as a single wire with a similar total energy deposition. On the
contrary, the single- and double-wire explosions show different emission
patterns in intensity and morphology. Namely, the double-wire explosion
develops a thin glowing layer between the two wires. It extends symmetrically
to a chamber scale, remains narrower than the gap between the wires, and
persists for times until their core material engulfs the entire chamber.

Turbulent electron heating and acceleration during the early stage of corona
interpenetration may enhance the line emission from the colliding neutral core
materials that trail coronas. This premise and the above considerations have
motivated the present study. Besides, we are interested in characterizing the
entire collision in terms of its scales, growth rates, and thresholds of
underlying instabilities. Spatio-temporal characteristics of the instabilities
and associated turbulent motion is another objective of this study.

The analytic and numerical studies in [MalkovSotnikov2018PhPl, Sotnikov2020]
have indicated that plasmas that stream against each other from two ohmically
vaporized wires may indeed mix in a relatively thin layer, in accord with the
preceding experiments [Sarkisov2005]. The shocked corona layer is separated
from the inflowing material by a pair of termination shocks (for supersonic
flows), and the plasma pressure inside the layer reaches its maximum at the
flow stagnation point. The emerging pressure gradient drives two symmetric
outflows perpendicular to the plane of the wires. These outflows are squeezed
between the inflowing coronas of the exploded wires (see, e.g., Figs.4-9 in
[Sotnikov2020]). A similar flow pattern forms a dense material of the wires
that follows the corona flow (Fig.10 in the above reference). The flow pattern
is qualitatively similar for the collisionless and hydrodynamic regimes. The
role of pressure in the former case plays electrostatic and pondermotive
potentials. The electrostatically mediated flow in the case of Boltzmannian
electrons is formally equivalent to the hydrodynamic flow with an adiabatic
index $\gamma=1$. This analogy justifies hydrodynamic simulations of colliding
plasmas [Sotnikov2020] as a guide for studying the microscopic processes to
which we turn in the present paper. Besides, the $\gamma=1$ value indicates
that the shocks are radiative, which is consistent with the enhanced emission
from a surprisingly thin shocked plasma layer observed in the experiments. The
thin layer points to high shock compression, characteristic of radiative
shocks.

The remainder of the paper is organized as follows. In the next section, we
consider instabilities that we expect to develop when two coronas from
exploded wires begin to penetrate each other. In Sec.II, we first evaluate the
electron response to the unstable waves and obtain their velocity distribution
established in the turbulent electric field of these waves in an adiabatic
(non-resonant) wave-particle interaction regime. Next, in Sec.III.2, we
include the electron-ion and electron-neutral collisions, thus making electron
dynamics irreversible. In particular, we will obtain the steady-state analytic
solutions to the problem and time-dependent, numerical ones. We conclude the
paper in Sec.IV with a summary of the results and a brief discussion.

##  II Two-Stream Instability in Counterstreaming Plasmas

When two coronas expanding from electrically exploded wires collide
supersonically, the collision region is bound by a pair of termination shocks
that propagate into the inflowing plasmas on each side of the shocked layer.
In both high- and low-Mach-number cases, an electrostatic potential builds up
between interpenetrating coronas and slows their ions down before the mixed
flow deflects off the plane of the two wires. To produce such deflection, the
electrostatic potential, $\Phi$, must reach the level of

$e\Phi\sim Mu^{2}$

Figure 1: Growth rate of the two-stream instability as described by eq.(4).

where $u$ is the speed of plasma upstream of the collision layer, and $M$ is
the ion mass. In more detail, a two-dimensional flow that emerges from the
collision of two concentric outflows has been studied using the hodograph
transformation in [MalkovSotnikov2018PhPl]. The magnetic field carried to the
collision layer with the flows is assumed to have no substantial effect on the
ion motion since the flow ram pressure is typically about an order of
magnitude higher than the magnetic pressure [Sarkisov2005]. This condition is
equivalent to the ion Larmor radius $\rho_{i}\sim u/\omega_{ci}$ being much
larger than the ion skin depth, which is the basic scale of a magnetized
collisionless shock. However, electrons are likely to be magnetized, at least
before being heated, so that the electric field produces an
$\boldsymbol{E}\times\boldsymbol{B}$ drift of electrons with the drift
velocity:

$V_{E}\sim c\Phi/LB$

where $L$ is the characteristic scale of the potential $\Phi$ that slows the
ions down. Therefore, the electron drift velocity $V_{E}$ can be estimated as
follows:

$V_{E}\sim\frac{\rho_{i}}{L}u.$

The collision region, where the counterstreaming flows slow down as they
approach the mid-plane, is expected to undergo the following. Initially, one
of the rapidly growing electron-scale instabilities sets in, depending on the
velocity of the relative motion and thermal velocities of the interpenetrating
plasmas. The typical speed of the relative motion is of the order of a several
$10$km/s [Sarkisov2005], whereas the electron thermal velocity is considerably
higher, up to $10^{3}$km/s. The Buneman-type two-stream electron instability
requires $u>V_{Te},$and probably does not start. Even if it does, it will
quickly thermalize and mix up the electron populations before the much heavier
ions start responding to this instability. The modified two-stream
instability, in which electrons ExB drift in the wire direction, has been
considered [MalkovSotnikov2018PhPl]. It should produce a similar effect of
electron preheating. After the saturation of these fast instabilities, a
slower instability of the counterstreaming ions will develop at their (slower)
time scale. A simple dispersion relation for this instability, written down in
the limit of monoenergetic counterstreaming ion beams, has the following form

$1+\chi_{e}-\frac{\omega_{pi}^{2}}{\left(\omega-
ku\right)^{2}}-\frac{\omega_{pi}^{2}}{\left(\omega+ku\right)^{2}}=0$ (1)

The last two terms come from the counterstreaming ions, moving at the bulk
speeds, $\pm u$. The electron susceptibility, $\chi_{e},$ can be taken in an
adiabatic approximation form

$\chi_{e}=\frac{2}{k^{2}\lambda_{De}^{2}},$ (2)

which requires $\omega\ll k_{z}V_{Te}$. This requirement is consistent with
the electron preheating during the early stage of collision. The electron
contribution in eq.(2) to the dielectric function corresponds to their density
perturbations related to that of the electrostatic potential as
$n_{\text{e}}^{\prime}/n_{0}\approx e$$\phi^{\prime}/T_{\text{e}}$. The
coefficient ’2’ in the above expression stands for the two electron
populations that neutralize the respective ion beams of a number density
$n_{0}$.

The role of the magnetic field is difficult to foresee, as it is of opposite
polarity in the inflowing coronas and may undergo a significant turbulent
reconnection/annihilation as the plasmas mix together. It is possible that as
the two plasmas expand adiabatically before they collide, the regime of hot
electrons, implied in eqs.(1) and (2), is established only after significant
interpenetration and heating of the coronas. During the initial collision
phase, an opposite regime of cold electrons, $k_{z}V_{Te}\ll\omega$, is more
likely to occur. Assuming that the unstable waves have a small wave vector
component along the magnetic field at this stage, $\chi_{e}$ takes a different
form. Indeed, if the electrons are magnetized in their motion across the local
magnetic field, $\omega_{ci}\ll\omega\ll\omega_{ce}$, instead of eq.(2) we can
write

$\chi_{e}=2\frac{\omega_{pe}^{2}}{\omega_{ce}^{2}}-2\frac{k_{z}^{2}}{k^{2}}\frac{\omega_{pe}^{2}}{\omega^{2}}$
(3)

The significance of this form of $\chi_{e}$ is a lower instability threshold
for certain wavenumbers. In this case, however, the instability saturates more
easily by electron heating. We will therefore override the phase of initial
electron heating and return to the adiabatic regime, which we justify below.
But, first, let us introduce the following dimensionless parameters:

$\xi=\frac{k^{2}u^{2}}{\omega_{pi}^{2}},\,\,\,\,\,\,\,\,\,\,\mathcal{M}^{2}=\frac{u^{2}}{C_{s}^{2}}\equiv
M\frac{u^{2}}{T_{e}}.$

Using this notation, eq.(1) rewrites

$1+\frac{2\mathcal{M}^{2}}{\xi}-\frac{1}{\left(\Omega-\sqrt{\xi}\right)^{2}}-\frac{1}{\left(\Omega+\sqrt{\xi}\right)^{2}}=0.$

One of its solutions describes an aperiodically growing mode, Fig.1:

$\frac{\omega^{2}}{\omega_{pi}^{2}}\equiv\Omega^{2}=-\frac{\xi}{\xi+2\mathcal{M}^{2}}\left(\sqrt{8\mathcal{M}^{2}+4\xi+1}-1-2\mathcal{M}^{2}-\xi\right)$
(4)

The aperiodic instability ($\omega^{2}<0$) develops when the term in
parentheses in eq.(4) is positive. The condition for that can be written in
the following short form:

$\mathcal{M}^{2}+\frac{\xi}{2}<1$ (5)

We observe that this condition can be met only in subsonic flows. Such flows
are expected in the plasma collision region, even if the flows are initially
highly supersonic, such as $u\gtrsim V_{Te}$. In this case, as we discussed,
the electrons should be heated via Buneman instability. Likewise, the
interpenetrating flows will also be subsonic when standing termination shocks
are formed on the two sides of the collision zone. We should note here,
though, that a situation in which the flows are supersonic but not to the
extent sufficient for triggering the Buneman instability is also possible. In
this case, electrons can be preheated by other instabilities, such as the
modified two-stream instability, leading to the generation of lower-hybrid
waves. This case has already been considered in Ref.MalkovSotnikov2018PhPl. We
may conclude from the above considerations that at least one of the discussed
instabilities will result in electron preheating. We will, therefore, continue
our analysis of a slower, two-stream ion instability under the assumption
$\omega\ll k_{z}V_{Te}$.

The above ion-driven instability can be suppressed if ions in the beam are
also heated by generated waves or arrive in the collision zone with a
significant velocity dispersion, carried over from the wire evaporation. We
assume the finite ion thermal speed in each ion beam, $V_{Ti}$, to examine
this possibility. Eq.(1) generalizes to the following [VVS_1961Usp] (see
Appendix)

$1+\frac{2}{k^{2}\lambda_{De}^{2}}-\frac{\omega_{pi}^{2}}{\left(\omega-
ku\right)^{2}-k^{2}V_{Ti}^{2}}-\frac{\omega_{pi}^{2}}{\left(\omega+ku\right)^{2}-k^{2}V_{Ti}^{2}}=0$
(6)

Again, coefficient 2 in the electron contribution accounts for the two
electron components that neutralize the two counterstreaming ion beams. These
are represented by the last two terms on the l.h.s. of the above equation.
Using the dimensionless variables introduced earlier, we rewrite this equation
as follows

$1+\frac{2}{k^{2}\lambda_{De}^{2}}-\frac{1}{\left(\Omega-\sqrt{\xi}\right)^{2}-k^{2}\lambda_{Di}^{2}}-\frac{1}{\left(\Omega+\sqrt{\xi}\right)^{2}-k^{2}\lambda_{Di}^{2}}=0$
(7)

A compact expression for the solution can be obtained using the parameters
$\mathcal{M}$ and $\zeta\equiv 2+\xi/\mathcal{M}^{2}$:

$\Omega^{2}\equiv\frac{\omega^{2}}{\omega_{pi}^{2}}=\left(\zeta-2\right)\left[\mathcal{M}^{2}+\frac{1}{\zeta}+\frac{T_{i}}{T_{e}}-\frac{1}{\zeta}\sqrt{4\zeta\mathcal{M}^{2}\left(1+\frac{T_{i}}{T_{e}}\zeta\right)+1}\right]$
(8)

For the onset of instability, the r.h.s of this expression should be negative.
This condition can be manipulated into the following inequality

$0<\zeta\left(\mathcal{M}^{2}-\frac{T_{i}}{T_{e}}\right)<2,$

which, in turn, can be rewritten as the following constraint on the Mach
number of the counterstreaming flows:

$\frac{T_{i}}{T_{e}}<\mathcal{M}^{2}<\frac{T_{i}}{T_{e}}+\frac{1}{1+k^{2}C_{s}^{2}/2\omega_{pi}^{2}},$

or on their velocity:

$\frac{T_{i}}{M}<u^{2}<\frac{T_{i}}{M}+\frac{C_{s}^{2}}{1+k^{2}C_{s}^{2}/2\omega_{pi}^{2}}.$
(9)

The growth rate given in eq.(2) is shown in Fig.2 as a function of $\xi$ and
$T_{i}/T_{e}$ for $\mathcal{M}=0.4$. The instability obviously, disfavors
short waves, for which even insignificant change of ion temperature or the
flow velocity $u$, or both, may produce a stabilizing effect. The electron
heating, on the contrary, widens the range of the flow velocities where the
instability may develop. Under these circumstances, it seems plausible to
investigate electron heating in more detail. Due to the low mass ratio,
electrons will be heated more efficiently than ions as the instability
develops. Electron heating broadens the velocity range of unstable modes and
impedes their stabilization, based on the above criterion.

##  III Electron Heating and Acceleration

At least in a linear approximation, the instability considered above excites
electrostatic perturbations that are non-propagating in the laboratory frame.
Physically, the instability may be thought of as intersecting backward-
propagating ion-acoustic modes in each of the two plasmas. They have phase
velocities $\left|\omega_{\pm}/k\right|=u$, where $\omega_{\pm}$ denote the
wave frequencies in the frames of the two counterstreaming plasmas. They
resonate with each other in the laboratory frame, in which they do not
propagate. So, the mode grows aperiodically.

###  III.1 Rapid Electron Preheating

The counterstreaming ion instability discussed in the previous section has a
growth rate $\gamma_{k}<\omega_{\text{pi}}$, so the electrons respond
adiabatically, as we assumed in deriving the dispersion relation, that is
$kV_{\text{Te}}\sim\omega_{\text{pi}}V_{\text{Te}}/u>\gamma_{k}$. Representing
then the fluctuating electric field as a superposition of growing modes,
$E=\sum_{k}E_{k}\left(t\right)\exp\left(i\boldsymbol{kr}\right)$, we can write
the following quasilinear equation for the heating process, e.g., [SagdGal69]

$\frac{\partial f_{e}}{\partial
t}=\frac{e^{2}}{m^{2}}\frac{\partial}{\partial\boldsymbol{v}}\sum_{k}\frac{\gamma_{k}\left|E_{k}\right|^{2}}{k^{2}v^{2}+\gamma_{k}^{2}}\frac{\partial
f_{e}}{\partial\boldsymbol{v}}$ (10)

This equation describes nonresonant interaction of electrons with the unstable
waves. Based on our assumption of electron adiabaticity, $kV_{Te}>\gamma_{k}$,
we can neglect the $\gamma_{k}^{2}$ in the denominator of the above
expression. Next, we introduce the electron velocity diffusion coefficient and
represent it in the form

$D_{e}=\frac{e^{2}}{2m^{2}}\frac{\partial}{\partial
t}\sum_{k}\left|\phi_{k}\right|^{2}=\frac{e^{2}}{m^{2}}\sum_{k}\gamma_{k}\left|\phi_{k}\right|^{2},$

Figure 2: Instability growth rate given by eq.(8) for _$\mathcal{M}=0.4$_ as a
function of $\xi$ and $g\equiv T_{i}/T_{e}$.

where we have denoted the spectral components of the wave electric potential
by $\left|\phi_{k}\right|=\left|E_{k}/k\right|$. Introducing then a new
variable

$\tau=\frac{e^{2}}{2m^{2}}\sum_{k}\left|\phi_{k}\left(t\right)\right|^{2},$

instead of $t$ and assuming first a one-dimensional spectrum of electrostatic
perturbations oriented along the flow, eq.(10) can be written in the following
short form

$\frac{\partial f_{e}}{\partial\tau}=\frac{\partial}{\partial
v}\frac{1}{v^{2}}\frac{\partial f_{e}}{\partial v}$

The last equation has a simple self-similar solution:

$f_{e}=\frac{C}{\tau^{1/4}}\exp\left(-\frac{v^{4}}{16\tau}\right)$

where C is a normalization constant. This solution describes reversible
electron heating that disappears together with the field oscillations. The
perceived electron “temperature”, $T_{e}^{\text{*}}=2\sqrt{\tau}$, in the
above formula merely reflects the amplitude of coherent electron oscillations
in an ensemble of waves with random phases.

In the case of three-dimensional isotropic fluctuations, the above equation
takes a slightly different form

$\frac{\partial
f_{e}}{\partial\tau}=\frac{1}{v^{2}}\frac{\partial^{2}f_{e}}{\partial v^{2}}$

but its solution has a similar dependence on the particle velocity:

$f_{e}=\frac{C}{\tau^{3/4}}\exp\left(-\frac{v^{4}}{4\tau}\right),$

as well as the same physical meaning as in the one-dimensional case.

While the adiabatic heating of electrons may impart significant amount of
kinetic energy, the energized electrons will cool down as soon as the
oscillations are damped or disappear otherwise, e.g., propagate away from the
excitation region. On the other hand, it is not difficult to imagine a
situation in which the entropy of electron gas increases, thus rendering the
heating process irreversible with no substantial energy input. It might come
about in two different ways or as a combination thereof. In one way, the
random but fixed phases of the waves begin to evolve dynamically with a
certain degree of stochasticity, which should randomize the electron
trajectories. In a second way, the electron orbits are randomized
independently of the oscillating field that sets them in motion. Such a
randomizing factor can be the scattering of electrons on ions and neutrals in
the shocked plasma layer. Since we expect such scattering in the mixing layer,
we include it in the electron kinetics below, which will make the adiabatic
heating irreversible.

###  III.2 Stochastic Heating and Runaway Acceleration of Electrons in The
Shocked Plasma Layer

Observations of enhanced emission near the flow stagnation region in
experiments, e.g., [Sarkisov2005], suggest that the shocked plasma thermalizes
most rapidly where the counterstreaming flows collide head-on. Due to their
rapid response to the instability, considered earlier, the two
counterstreaming electron distributions will couple and form essentially one
distribution with a zero average drift speed across the layer. At the same
time, ions will continue counterstreaming, thus providing free energy for the
instability considered in the preceding section. The instability should,
however, eventually saturate. The most straightforward saturation mechanism is
the convection of electrostatic perturbations with plasma outflow from the
instability region in the layer direction.

Meanwhile, the merger of electron distributions cannot be completed everywhere
in the mixing layer, especially at its edges. s. Therefore, we first describe
electrons that inflow into the mixing region from its opposite sides as
separate components and denote their distributions as $f_{\pm}$, so that the
full electron distribution is simply $f=f_{\text{+}}+f_{-}$. The distributions
$f_{\pm}$ enter the collision layer with drift velocities $\pm u$,
respectively. The equations for $f_{\pm}$ can be written as follows

$\frac{\partial f_{\pm}}{\partial t}+\left(v\mu\pm u\right)\frac{\partial
f_{\pm}}{\partial z}-\frac{e}{m}E\left(\mu\frac{\partial f_{\pm}}{\partial
v}+\frac{1-\mu^{2}}{v}\frac{\partial
f_{\pm}}{\partial\mu}\right)=\frac{1}{v^{2}}\frac{\partial}{\partial
v}\left[v^{2}\nu_{\text{e}}\left(V_{\text{Te}}^{2}\frac{\partial
f_{\pm}}{\partial
v}+vf_{\pm}\right)\right]+\frac{\nu}{2}\frac{\partial}{\partial\mu}\left(1-\mu^{2}\right)\frac{\partial
f_{\pm}}{\partial\mu}$ (11)

where we have denoted

$\nu_{\text{e,i}}=\frac{4\pi
e^{4}N_{\text{e,i}}}{m^{2}v^{3}}\ln\left(\frac{mv^{2}\lambda_{\text{De}}}{e^{2}}\right)$
(12)

and

$\nu=\nu_{\text{i}}+\nu_{\text{e}}\left(1-\frac{V_{\text{Te}}^{2}}{v^{2}}\right)$
(13)

Below, explain the two equations for $f_{\pm}$ combined in eq.(11) for
simplicity. Formally, they are written in separate reference frames moving
with velocities $\pm u$ in $z-$direction, so that the velocity projection on
$z$ direction is $v\mu\pm u$ in each of these equations. The cosine of the
particle pitch-angle to the shock normal ($z$\- direction), $\mu$, is
introduced in a standard form $v_{z}=v\mu$, where $v_{z}$ is the $z-$
component of electron velocity $\boldsymbol{v}$. Apart from the terms with
$\pm u$ on the lhs of eqs.(11), the equations for $f_{\pm}$ are identical.
Moreover, since $v$ is at least $v\sim V_{\text{Te}}\gg V_{\text{Ti}}\sim u$,
or even $v\gg V_{\text{Te}}$, since we expect significant electron heating and
acceleration by turbulent electric field, $E$, we assume that
$f_{\text{+}}\approx f_{-}$ inside the shocked layer. The role of the $\pm u$
term is then limited to the edges of the layer where initially cold electrons
enter the layer from opposite sides. Likewise, the pitch-angle scattering
operator (the last term on the rhs) can be regarded as acting in the
laboratory rather than separately in each of the comoving reference frames, as
originally required in eq.(11). Except for $z-$dependence (introduced
straightforwardly), separation of the two groups of electrons $f_{\pm}$, and
an oscillatory rather than the constant electric field, $E$, eq.(11) is
similar to that used in studies of electron runaway, e.g., [Gurevich61RunAway,
KruskalBernstein64].

####  III.2.1 Equation for runaway electrons

Taking the above discussion into account, introducing pitch-angle averaged
quantities

$\left\langle\cdot\right\rangle=\frac{1}{2}\int_{-1}^{1}\left(\cdot\right)d\mu,$

and averaging eq.(11) over $\mu$, we arrive at the following equation for
$f=f_{\text{+}}+f_{-}$:

$\frac{\partial\left\langle f\right\rangle}{\partial
t}+\frac{1}{2}\left(\frac{eE}{mv^{2}}\frac{\partial}{\partial
v}v^{2}-v\frac{\partial}{\partial
z}\right)\left\langle\left(1-\mu^{2}\right)\frac{\partial
f}{\partial\mu}\right\rangle=\frac{1}{v^{2}}\frac{\partial}{\partial
v}\left[v^{2}\nu_{\text{e}}\left(V_{\text{Te}}^{2}\frac{\partial\left\langle
f\right\rangle}{\partial v}+v\left\langle
f\right\rangle\right)\right]+\left\langle S\right\rangle$ (14)

Here the term $\left\langle S\right\rangle$ arose from summing up the fluxes
$\pm uf_{\pm}$ through the boundaries of the layer. They are not accounted for
automatically when we add the equations for $f_{\pm}$ together in eq.(11) and
assume $f_{\text{+}}=f_{-}$. The reason for this extra term is that
$f_{\text{+}}$$\neq f_{-}$ near the boundaries, which we mentioned earlier. We
may include in the term $\left\langle S\right\rangle$ other transport
channels, such as particle convection along $x$ and y axes. To obtain an
evolution equation for $\left\langle f\right\rangle$, we need to express the
$\left\langle\left(1-\mu^{2}\right)\partial f/\partial\mu\right\rangle$ in the
last equation through $\left\langle f\right\rangle$. This can be done by
applying Chapman-Enskog method to eq.(11). This method will provide an
expansion of equation for $\left\langle f\right\rangle$ in powers of
$\overline{E^{2}}$\- an average value of $E^{2}$ over the distance and time
that exceeds the mfp $\lambda=v/\nu$ and collision time $1/\nu$, respectively.
In particular, if the electric field comes from the instability of
counterstreaming ion beams considered earlier,
$E=\left(1/\sqrt{2}\right)\sum_{k}E_{k}\left(t\right)\exp\left(ikz\right)+c.c.$,
then $\overline{E^{2}}=\sum_{k}\left|E_{k}\right|^{2}$.

The Chapman-Enskog method selects the zeroth order approximation that makes
the rhs (particle collision part) vanish. This is the assumption for a system
being close to thermal equilibrium. However, the collision operator on the rhs
of eq.(11) acts on the angular part and the absolute value of particle
velocity separately. Now we consider the second term (angular part) on the rhs
of eq.(11) to be dominant. Although it is formally of the same order as the
first term, it must vanish in the zeroth order by virtue of eq.(14), as we can
set $\left\langle S\right\rangle\approx 0$ inside the shocked plasma, as
discussed earlier. Also note that for a Maxwellian distribution,
$f\propto\exp\left(-v^{2}/V_{\text{Te}}^{2}\right)$, the collision term
vanishes. As we shall see, for higher velocities the runaway due the secular
action of electric field dominates the thermal equilibtration represented by
the rhs in eq.(14). This dominance justifies neglecting the first term on the
rhs of eq.(11) for the purpose of calculation of anisotropic correction to the
equilibrium state.

Turning to the lhs of eq.(11) we note that the first term ($\partial
f/\partial t$) needs to be dropped to this order of approximation as it is
responsible for a short-time relaxation of the anisotropic component, e.g.
[MS2015]. From eq.(11), to the first order, we thus write

$\frac{\partial
f}{\partial\mu}\approx\frac{1}{\nu}\left(\frac{eE}{m}\frac{\partial\left\langle
f\right\rangle}{\partial v}-v\frac{\partial\left\langle
f\right\rangle}{\partial z}\right)$

On substituting this expression into eq.(14), the latter takes the following
form

$\frac{\partial\left\langle f\right\rangle}{\partial
t}=\frac{1}{v^{2}}\frac{\partial}{\partial
v}v^{2}\left[\left(\nu_{\text{e}}V_{\text{Te}}^{2}+\frac{e^{2}\overline{E^{2}}}{3m^{2}\nu}\right)\frac{\partial\left\langle
f\right\rangle}{\partial v}+\nu_{\text{e}}v\left\langle
f\right\rangle\right]+\left\langle
S\right\rangle+\frac{v^{2}}{3\nu}\frac{\partial^{2}\left\langle
f\right\rangle}{\partial z^{2}}$ (15)

The last term on the rhs governs the diffusive losses from the shocked layer.
They might extend to the other coordinates, which can symbolically be included
into the $\left\langle S\right\rangle$-term, as we mentioned earlier. The
losses become important for large $v$, where they equate the acceleration
$\propto\overline{E^{2}}$. The acceleration occurs in the plasma mixing layer
that is bound by a pair of termination shocks at $z=\pm l/2$. A plausible
assumption about the distribution of accelerated electrons in $z$\- direction
is that their density vanishes at $z=\pm l/2$. The reason for this assumption
is twofold. First, the shocked layer has an elevated electrostatic potential
that slows down the inflowing ions in $z$ direction. Hence, this region is a
potential well for electrons, so they tend to have a higher concentration
toward the middle of this region. Second, the high-energy electrons, for which
the losses become significant, most likely never return to the shocked layer
if they cross one of its edges. Indeed, the flow upstream of the termination
shocks is relatively laminar. This means that hot electrons accelerated in the
mixing layer by the turbulent electric field are not deflected strongly from
their paths. Their transport regime changes at the edge from diffusive to
ballistic. Depending on the electron energy and magnetic field strength,
electrons leave the layer by either following rectilinear trajectories or by
ExB drifting away from the flow stagnation region across the magnetic field
convected with the upstream flow. Therefore, we may use the following
expansion for $\left\langle f\right\rangle$ in $z$\- direction

$\left\langle
f\right\rangle=\sum_{n=1}^{\infty}f_{n}\left(t,v\right)\cos\left(\frac{\pi
nz}{l}\right)$ (16)

While this representation implies instantaneous losses of hot electrons from
the shocked layer as soon as they reach its boundaries at $z=\pm l/2$, there
is also an influx of cold electrons from the upstream. By keeping only the
leading term in the last Fourier expansion, eq.(15) can be rewritten as
follows

$\frac{\partial f_{1}}{\partial t}=\frac{1}{v^{2}}\frac{\partial}{\partial
v}v^{2}\left[\left(\nu_{\text{e}}V_{\text{Te}}^{2}+\frac{e^{2}\overline{E^{2}}}{3m^{2}\nu}\right)\frac{\partial
f_{1}}{\partial
v}+\nu_{\text{e}}vf_{1}\right]+S_{1}-\frac{\pi^{2}v^{2}}{3l^{2}\nu}f_{1}$ (17)

Here, we have included the input of the $\left\langle
S\left(z\right)\right\rangle$-term in eq.(15) and denoted it by $S_{1}$. Note
that both $V_{\text{Te}}$ and $\overline{E^{2}}$ in the last equation,
generally speaking, depend on time. We start, however, from considering a
steady state solution of the above equation, assuming that these quantities
have already attained their time-asymptotic magnitudes.

####  III.2.2 Steady State Fermi acceleration of Runaway Electrons

Since the source term $S_{1}$ in eq.(17) is associated with the influx of
relatively cold electrons into the shocked plasma layer, it is plausible to
place it at the origin in the velocity space. Thus, we set
$S_{1}=S\delta\left(v\right)/v^{2}$. Then, a steady state solution for $v>0$
corresponds to a constant electron flux $S$ between the source and the sink
(the last term in eq.[17] controls the particle losses) in velocity space:

$S=-v^{2}\left[\left(\nu_{\text{e}}V_{\text{Te}}^{2}+\frac{e^{2}\overline{E^{2}}}{3m^{2}\nu}\right)\frac{\partial\left\langle
f\right\rangle}{\partial v}+\nu_{\text{e}}v\left\langle
f\right\rangle\right]=cosnt$ (18)

Figure 3: Left panel: steady state solutions of eq.(21) plotted for different
$\vartheta$. To emphasize the development of the runaway tail, the
normalization constant $\sigma$ in each curve is chosen to satisfy the
condition $\left\langle f\left(0\right)\right\rangle=1$, for all $\vartheta$
rather than fixing the total number of particles. Right panel: the same
solution shown as a function of $\vartheta$ for different particle energies,
$w$.

This spectrum will cut off at sufficiently large $v$ due to losses, which we
have neglected for now because they grow very sharply at high velocities (see
eq.(13)) and are thus not important in the velocity range considered. For
relatively low values of $\overline{E^{2}}$ and not very large $v$ (for which
$\nu$ and $\nu_{\text{e}}$ would become small), the term with $E^{2}$ can be
neglected and the solution for $\left\langle f\right\rangle$ is close to a
Maxwellian (because $S=0$ as well for $v>0$). The electron runaway effect of
the electric field becomes important for $v\gg V_{\text{Te}}$, and we can
substitute $\nu\approx 2\nu_{e}$ in eq.(13). Furthermore, the expression under
the logarithm in eq.(12) is typically so large that the velocity variation
under the logarithm can be neglected, and $v$ can be set to a value
$V_{\text{*}}$ at which the acceleration and dynamical friction terms in
eq.(17) equate, that is

$\frac{4\pi
e^{4}N_{\text{e,i}}}{m^{2}v^{2}}\ln\left(\frac{mV_{\text{*}}^{2}\lambda_{\text{De}}}{e^{2}}\right)\sim\frac{e\sqrt{\overline{E^{2}}}}{m}$

By defining then the Dreicer’s field, familiar from the theory of runaway
electrons [Gurevich61RunAway, KruskalBernstein64] as

$E_{\text{D}}=\frac{4\pi
e^{3}N_{\text{e,i}}}{T_{\text{e}}}\ln\left(\frac{mV_{\text{*}}^{2}\lambda_{\text{De}}}{e^{2}}\right),$

the critical velocity $V_{\text{*}}$ above which the runaway process dominates
dynamical friction can be expressed through $E_{\text{D}}$ as follows:

$V_{\text{*}}=V_{\text{Te}}\frac{E_{\text{D}}^{1/2}}{\left(\overline{E^{2}}\right)^{1/4}}$
(19)

It is now convenient to use dimensionless electron energy
$w=v^{2}/6^{1/3}V_{\text{Te}}^{2/3}V_{\text{*}}^{4/3}$ and field strength
parameter

$\vartheta=6^{1/3}V_{\text{*}}^{4/3}/2V_{\text{Te}}^{4/3}.$ (20)

Using these quantities, eq.(18) can be rewritten as follows

$\left(1+w^{3}\right)\frac{\partial\left\langle f\right\rangle}{\partial
w}+\vartheta\left\langle
f\right\rangle=-S\vartheta\frac{m^{2}}{eT_{\text{e}}E_{\text{D}}}\equiv-\sigma$
(21)

This equation can be solved in quadratures:

$\left\langle
f\right\rangle=\sigma\exp\left[-\frac{\vartheta}{\sqrt{3}}\tan^{-1}\left(\frac{2w-1}{\sqrt{3}}\right)\right]\frac{\left(1-w+w^{2}\right)^{\vartheta/6}}{\left(1+w\right)^{\vartheta/3}}\int_{w}^{\infty}\frac{\left(1+w\right)^{\vartheta/3-1}dw}{\left(1-w+w^{2}\right)^{\vartheta/6+1}}\exp\left[\frac{\vartheta}{\sqrt{3}}\tan^{-1}\left(\frac{2w-1}{\sqrt{3}}\right)\right]$
(22)

Fig.3 (left panel) shows the above solution for different values of field-
strength parameter $\vartheta$. It is seen that for strong electric fields,
$\vartheta\sim 1$, the power-law tail, $\left\langle f\right\rangle\propto
w^{-2}\propto v^{-4}$, grows immediately from the thermal velocity. For weak
fields ($\vartheta\gg 1)$ the Maxwellian core is distinct from the power-law
tail of runaway electrons. On the right panel of Fig.3, the electron
distribution function $\left\langle f\right\rangle$ is shown as a function of
$\vartheta$ for several values of the fixed electron energies.

The asymptotic power-law energy spectrum $\left\langle f\right\rangle\propto
w^{-2}$ may be interpreted as a classical Fermi acceleration at work. Fermi
[Fermi49] showed that under competing particle acceleration-loss events, a
power-law energy spectrum, $w^{-q}$, is naturally established. Its index $q$,
which is in our case, $q=2$, is determined by the ratio of particle
confinement time, $\tau_{\text{loss}}$, and its acceleration time,
$\tau_{\text{acc}}$. Fermi originally devised this mechanism for galactic
cosmic rays, in which case $\tau_{\text{loss}}$ is an average time spent by a
CR particle in the galaxy before it escapes into the intergalactic space. The
acceleration time, in turn, is simply the energy e-folding time,
$dE/dt=E/\tau_{\text{acc}}$, where $E$ the particle energy. According to
Fermi, the power-law index, $q$, in the particle distribution, $f\propto
E^{-q}$, can be calculated as follows:

$q=1+\frac{\tau_{\text{acc}}}{\tau_{\text{loss}}}$

The Fermi acceleration has been recognized later as a general mechanism
pertinent to many other systems where a nonthermal spectrum emerges from the
balance between the particle energy gain and its loss or particle escape. In
our case, accelerated particles loose energy as they scatter on slowly moving
particles (thermal electrons, ions, and neutrals). Between each scattering
they gain energy from the electric field. So, the energy gain and loss occur
during the same time interval, $\tau_{\text{loss}}=\tau_{\text{acc}}$. The
result $q=2$ immediately follows from the Fermi’s celebrated formula.

The asymptotic $q=2$ acceleration regime is pertinent to high-energy
particles, for which the contribution of thermal core into the flux $S=const$
in eq.(18) can be neglected (all terms other than the acceleration term,
$\propto\overline{E^{2}}$). The acceleration term, however, includes both the
actual energy gain between two consecutive collisions and the particle
scattering that results from these collisions. This conclusion can be inferred
from the derivation of eq.(17) in Sec.III.2.1. The full steady state solution
in eq.(22) does include the thermal core part of the spectrum and it has a
somewhat complex form. For the high-energy particles, the asymptotic spectrum
can be most easily obtained from the relation $w^{3}\partial\left\langle
f\right\rangle/\partial w=-\sigma=const$, according to the insignificant role
of the thermal core for particles with $w\gg 1$.

####  III.2.3 Time-dependent runaway

It is convenient to study time-dependent solutions of eq.(17) using a new
variable $U$ instead of the particle energy $w$ or velocity $v$. These
variables are related in the following way:

$U=w^{3/2}=\frac{v^{3}}{\sqrt{6}V_{\text{Te}}V_{\text{*}}^{2}}.$

By measuring time $t$ in the units of electron acceleration time,
$\tau_{\text{a}}$ ($t/\tau_{\text{a}}\to t$), and introducing the particle
loss parameter, $\lambda$, owing to their diffusion in $z-$ direction,

$\tau_{\text{a}}^{-1}=6^{1/6}\frac{e}{m}E_{D}\frac{V_{\text{Te}}^{7/3}}{V_{\text{*}}^{10/3}},\;\;\;\;\;\lambda=\frac{3^{4/3}m^{2}}{4E_{D}^{2}e^{2}}\frac{V_{\text{*}}^{20/3}}{l^{2}V_{Te}^{8/3}}$

eq.(17) rewrites as follows

$\frac{\partial f_{1}}{\partial t}=\frac{\partial}{\partial
U}\left[\frac{3}{2}U^{1/3}\left(1+U^{2}\right)\frac{\partial f_{1}}{\partial
U}+\vartheta f_{1}\right]+\tau_{a}S_{1}-\lambda U^{5/3}f_{1}$ (23)

This equation contains two parameters, $\vartheta$ and $\lambda$; these are
the rates of collisional energy losses of runaway electrons and their escape
through the boundaries of the mixing layer, respectively. Besides, the
parameter $\vartheta$ determines the thermal core temperature that establishes
through a balance between dynamical friction (term proportional to
$\vartheta$) and velocity diffusion (term proportional to $1+U^{2}$, where the
first term in the sum corresponds to the thermal velocity diffusion due to
Column collisions, while $U^{2}$ stands for the anomalous velocity diffusion).
The part of the velocity diffusion $\propto U^{2}$ stems from the electric
field generated by the ion beam instability discussed in Sec.II.

According to our discussion in the preceding section, the source supplies
particles at lower energies, that is at $U\ll 1$. To keep the treatment
simple, we therefore draw equivalence between the source term and the boundary
condition at the lower end of the integration domain. In other words, instead
of specifying the source profile $S_{1}\left(U\right),$ we keep the value
$f\left(U_{\text{min}}\right)$ fixed as a boundary condition. It is clear that
acceleration is negligible near this end as long as $U_{\text{min}}\ll 1$.
Also, to circumvent an insignificant problem with the algebraic singularity,
$U^{1/3}$, at $U=0$ in the above equation, we set $U_{\text{min }}=0.1$
instead of $U_{\text{min }}=0$. This simplification is justified by the upper
end of the integration interval being at $U_{\text{max}}\gg 1$. We then impose
$f_{1}\left(U_{\text{min}}\right)=1$ and $f_{1}\left(U_{\text{max}}\right)=0$
as a complete set of boundary conditions for eq.(23). The latter condition
implies that particles accelerated to $U_{\text{max}}$ escape the mixing layer
regardless of the value of $\lambda$, but the effect of finite $\lambda$ is
still present in the solution.

Figure 4: Time evolution of the spectrum of runaway electrons, as described by
eq.(23), starting from a relatively cold initial distribution with
$U_{\text{th}}=0.01$. Four combinations of equation parameters $\vartheta$ and
$\lambda$ are shown in each plot. The left two panels show the loss-free
solutions, $\lambda=0$, except for a particle instantaneous escape upon
reaching the maximum energy at $U_{\text{max}}.$

Fig.4 illustrates the electron runaway acceleration, as it evolves in time,
shown for four combinations of parameters $\vartheta$ and $\lambda$. All four
runs start from a relatively cold Maxwellian (if expressed in terms of
particle velocity rather than $U$):

$f_{1}(0,U)=\exp\left[-\frac{U^{2/3}-U_{\text{min}}^{2/3}}{U_{\text{th}}^{2/3}}\right]$
(24)

where $U_{\text{th }}$ characterizes the initial electron temperature.
Although the electron adiabatic preheating, as described in Sec.III.1, should
result in a steeper decay of electron distribution, we assume the above
Maxwellian form. The reason for that choice is twofold. First, the initial
distribution is too narrow to significantly affect its evolution beyond a
short initial phase. Second, the collisional relaxation of any initial
distribution localized at $U\ll 1$ dominates the evolution of
$f_{1}\left(t,U\right)$ at early times. It quickly relaxes to
$f_{1}\propto\exp\left(-\vartheta U^{2/3}\right)$, as can be seen from
eq.(23). This initial dynamics is seen in the two left panels of Fig.4.
Indeed, for the chosen two values of dynamical friction parameter,
$\vartheta=12$ and $\vartheta=18$, the electron distribution broadens to the
respective Maxwellian, $\exp\left(-\vartheta U^{2/3}\right)$ at $U\lesssim 1$
very quickly and proceeds then to higher energies. At $U\gtrsim 1$, the
particle distribution develops a suprathermal tail that also very quickly
“forgets” the initial velocity distribution.

While the initial distribution of electrons is not significant for the
development of the suprathermal tail, the tail amplitude is very sensitive to
the dynamical friction parameter, $\vartheta$ (see also Fig.3). As seen from
Fig.4, for a 50% increase in $\vartheta$ from the top panel to the bottom
ones, the tail height drops by about two orders of magnitude. This is similar
to the conventional runaway process (in a constant electric field) as the
suprathermal tail is detached from the thermal Maxwellian at a critical
velocity, $V_{\text{*}}$, which in our case is defined in eq.(19) and is
related to the dynamical friction, $\vartheta$. Because of the exponential
decay of Maxwellian distribution at $v\sim V_{\text{*}}$, small variations of
$V_{\text{*}}$ result in large changes in the tail height, connected to the
Maxwellian at $v\approx V_{\text{*}}$.

##  IV Conclusions and Discussion

In this paper, we have studied microscopic processes occurring in a mixing
layer between two colliding plasma clouds. Our major findings are as follows:

  1. 1.

Initial relative velocity of interpenetrating plasmas needs to be in the range
between the ion thermal velocity, $V_{\text{Ti}}=\sqrt{T_{\text{i}}/M}$, and
ion-sound velocity, $C_{\text{s}}=\sqrt{T_{\text{e}}/M}$, for triggering a
strong anomalous coupling of the plasma streams that is supported by the two-
stream instability.

  2. 2.

Subsequent nonresonant (reversible) heating of electrons then broadens the
range of unstable wavenumbers.

  3. 3.

Inclusion of electron-ion and electron-neutral collisions makes the electron
energization irreversible.

  4. 4.

In addition to the heating of an electron core distribution, a suprathermal
power-law tail is produced in wave-particle interactions. It decays with
energy as $w^{-2}$, thus showing a notable analogy with the stochastic Fermi
acceleration.

  5. 5.

The amplitude of the power-law electron tail strongly depends on the ratio of
the turbulent electric field to the electron collision frequency.

While the instability condition in (1) may appear to be restrictive, the
following arguments can be advanced in its favor. Since each of the colliding
plasmas undergoes a transition from supersonic to subsonic speed in the
collision region, the instability condition is likely to be met at least in a
limited part of the collision region. The subsequent electron heating will
then broaden the unstable range of relative velocity.

Under favorable instability conditions, specified above and in more detail in
Sec.II, a significant fraction of kinetic energy of counterstreaming corona
ions is likely to be converted in ion-acoustic waves. When also the electron-
ion and electron-neutral collision frequencies are limited by the condition
$\vartheta=6^{1/3}V_{\text{*}}^{4/3}/2V_{\text{Te}}^{4/3}\lesssim 10$ (here
$V_{\text{*}}$, eq.[19], is a critical electron velocity beyond which they
enter a runaway regime), the wave energy is efficiently transferred to thermal
and nonthermal electrons. The anomalous (wave-driven) electron heating and
acceleration will likely to have a defining impact on the next phase of plasma
mixing.

The next phase starts when heavier wire-core materials, which follow the light
corona plasmas, begin to interpenetrate. As the materials are largely in a
neutral or weakly ionized state at the moment of touching each other, the
presence of hot and suprathermal coronal electrons in the mixing layer will
affect the interpenetration process by primarily ionizing the neutrals and
producing line emission. Although this phase of the plasma collision was out
of the scope of this paper, focused on the coronal phase, it is worthwhile to
briefly discuss possible scenarios that may be derived from the presence of
hot and suprathermal electrons.

As we mentioned, energetic electrons will ionize counterstreaming wire-core
atoms, thus partially converting them into counterstreaming ion beams. The
plasma stability analysis then becomes similar to that of the coronal
collision studied earlier. However, several differences are evident. First,
the beams will likely become multi-charged as the ionization potential of wire
material (e.g., aluminum, tungsten) grows significantly with the charge
number. Second, as the energetic electron population left behind after the
coronal phase of plasma collision has a much lower concentration than the
wire-core atoms, the ionization is likely to occur in an avalanche regime
[GurevichUFN2001]. The neutral gas electric breakdown supported by the
oscillatory electric field generated by the two-stream instability of freshly
ionized beams may then naturally explain the protracted line emission observed
in experiments with double-wire electric explosion [Sarkisov2005]. Note, that
this phenomenon is not observed in a single-wire experimental setup. Another
important aspect of these experiments is that the voltage collapse in the
vacuum chamber occurs when the light emission starts. Therefore, the electric
discharge is not driven by the external electric field, which also speaks to
an internal mechanism of the observed glow in the mixing layer.

Strong radiation from the mixing layer suggests that its bounding shocks may
be radiative. The removal of the shocked plasma internal energy by line
emission softens the equation of state thus driving the adiabatic index to
$\gamma\approx 1$, which makes the shocks more compressive. The plasma layer
between them becomes thinner. Our hydrodynamic simulation [Sotnikov2020] of
colliding concentric wire outflow have shown that the low value of
$\gamma$$\approx 1$ is required to keep the mixing layer thin enough and not
to expand under the internal pressure of shocked plasma. The layer remains way
inside the gap between the original wire positions, as observed in experiments
[Sarkisov2005].

###### Acknowledgements.

This work was supported by the Air Force Office of Scientific Research under
the grant LRIR No. 19RYCOR062. Public release approval record: AFRL-2021-4552

##  Appendix A Dispersion Equation for Interpenetrating Plasma Flows

Thermal corrections to the plasma dispersion function are often given in an
additive form. It is more restrictive to the magnitude of these corrections
than those we used in eq.(6). As a reference, we provide below its simple
derivation, following the work [VVS_1961Usp]. The equations of motion for each
of the ion beams moving with the bulk velocities $\pm u,$ with the respective
perturbation, $u_{\pm}^{\prime},$ is

$\frac{\partial u^{\prime}}{\partial t}\pm u\frac{\partial
u^{\prime}}{\partial x}=-\frac{1}{Mn_{i}}\frac{\partial p^{\prime}}{\partial
x}-\frac{e}{M}\frac{\partial\phi}{\partial x}$ (25)

To lighten the notation, we have omitted the subscripts $\pm$ in $u^{\prime}$,
$p^{\prime}$, $n^{\prime}$ for either beam, as it should not cause any
confusion. The respective equations for $u_{\pm}^{\prime}$, etc., are marked
by the sign in front of the beam velocity, $\pm u$. In the above equations, we
have denoted the pressure perturbation as $p^{\prime}$, while $n_{i}$ below
denotes the unperturbed ion density of each beam and $\phi$ is the
perturbation of electrostatic potential. The continuity equation for the
density perturbation, $n_{i}^{\prime}$

$\frac{\partial n_{i}^{\prime}}{\partial t}\pm u\frac{\partial
n_{i}^{\prime}}{\partial x}+n_{i}\frac{\partial u^{\prime}}{\partial x}=0$
(26)

while the Poisson equation

$\frac{\partial^{2}\phi}{\partial x^{2}}=4\pi
e\sum_{\pm}\left(n_{e}^{\prime}-n_{i}^{\prime}\right)$ (27)

where $n_{e}^{\prime}$ denotes the total electron density perturbation, while
the sum is taken over the two ion beams and the respective electron
contributions that are initially moving at the same speed as the ions. We
relate the ion pressure perturbation with that of their density, assuming a
constant temperature, $T_{i},$ $p^{\prime}=T_{i}n_{i}^{\prime},$while a
Boltzmann distribution for each electron component simply yields
$n_{e}^{\prime}/n_{i}=e\phi/T_{e}.$ Assuming that the perturbed quantities
depend on $x$ and $t$ as $\exp\left(ikx-i\omega t\right)$, expressing
$n_{i}^{\prime}$ through $\phi$ using eqs.(25-26), and substituting in
eq.(27), we arrive at eq.(6)


# Branching Exponents of Synthetic Vascular Trees under Different Optimality
Principles

Etienne Jessen
Institute for Mechanics, Computational Mechanics Group
Technical University of Darmstadt
64287 Darmstadt, Germany
<EMAIL_ADDRESS>
& Marc C. Steinbach
Institute of Applied Mathematics
Leibniz Universität Hannover
30167 Hannover, Germany
<EMAIL_ADDRESS>
& Charlotte Debbaut
IBiTech – Biommeda
Ghent University
Ghent, Belgium
<EMAIL_ADDRESS>
& Dominik Schillinger
Institute for Mechanics, Computational Mechanics Group
Technical University of Darmstadt
64287 Darmstadt, Germany
<EMAIL_ADDRESS>

###### Abstract

Objective: The branching behavior of vascular trees is often characterized
using Murray’s law. We investigate its validity using synthetic vascular trees
generated under global optimization criteria. Methods: Our synthetic tree
model does not incorporate Murray’s law explicitly. Instead, we assume it
holds implicitly and investigate the effects of different physical constraints
and optimization goals on the branching exponent that is now allowed to vary
locally. In particular, we include variable blood viscosity due to the
Fåhræus–Lindqvist effect and enforce an equal pressure drop between inflow and
the micro-circulation. Using our global optimization framework, we generate
vascular trees with over one million terminal vessels and compare them against
a detailed corrosion cast of the portal venous tree of a human liver. Results:
Murray’s law is implicitly fulfilled when no additional constraints are
enforced, indicating its validity in this setting. Variable blood viscosity or
equal pressure drop leads to deviations from this optimum, but with the
branching exponent inside the experimentally predicted range between 2.0 and
3.0. The validation against the corrosion cast shows good agreement from the
portal vein down to the venules. Conclusion: Not enforcing Murray’s law
explicitly reduces the computational cost and increases the predictive
capabilities of synthetic vascular trees. Significance: The ability to study
optimal branching exponents across different scales can improve the functional
assessment of organs

_Keywords_ branching exponents $\cdot$ Fåhræus–Lindqvist effect $\cdot$ human
liver $\cdot$ Murray’s law $\cdot$ synthetic vascular trees $\cdot$ vascular
corrosion cast

##  1 Introduction

The cardiovascular system is responsible for transporting blood to and from
all cells in the human body, leading to hierarchical networks of vessels,
called vascular trees, inside each organ. According to Murray [1], this
hierarchy obeys scaling relations based on the minimization of the total
energy expenditure of the system. Many factors influence and constrain this
minimization process, such as the type and shape of the organ supplied, the
demand for the organ’s cells, and the existence of vascular diseases. The
goals and constraints guiding their structural development and influence on
the vascular system have yet to be entirely understood, even though extensive
work has been carried out for over a century. Thus the analysis of vascular
diseases based on the anatomy and physiology of the vascular structure remains
a challenge.

Murray first described a minimization problem for vascular segments in 1926
[1, 2]. Here, a tree is approximated as a bifurcating network consisting of
rigid tubes, and the physical principles for fluid flow follow Poiseuille’s
law. The goal is to minimize the total power of the tree network with its
minimum characterized by _Murray’s law_. It describes the relationship of the
radius of a parent vessel $r_{0}$ against the radii of its children’s vessels
$\left(r_{1},r_{2}\right)$ as a power law with

$\displaystyle r_{0}^{\gamma}=r_{1}^{\gamma}+r_{2}^{\gamma}.$ (1)

The branching exponent $\gamma$ became an essential parameter for
characterizing the branching behavior of vascular trees. In Murray’s original
formulation, $\gamma=3.0$ is constant across the entire network.

An extensive number of studies have been conducted to investigate Murray’s law
experimentally [3, 4, 5]. In general, exponents between $2.0$ and $3.0$ were
measured. In [6], exponents were observed even going over the theoretical
limit of 3.0 with $\gamma=3.2$. Multiple theoretical studies have analyzed the
possible factors contributing to these branching behaviors. An extension to
Murray’s law was proposed by Uylings [7], which incorporated the effects of
turbulent flow into the minimization problem. Results show branching exponents
as low as $2.33$ for turbulent flow. In [8], a vascular model was
investigated, which considered the role of elastic tubes. Compared to rigid
tubes, the effect of pulsatile flow lowered the optimal value to $2.3$. Zhou,
Kassab, and Molloi [9, 10] generalized Murray’s law hypothesis to an entire
coronary arterial tree by defining a vessel segment as a stem and the tree
distal to the stem as a crown. They showed that $\gamma$ deviates from $3.0$
even for steady-state flow and depends on the ratio between metabolic demand
and viscous power dissipation.

An alternative approach to investigate branching exponents is to construct
vascular trees synthetically. The most well-known generation method here is
constrained constructive optimization (CCO) [11]. The local optimization
approach is directly based on Murray’s minimization principles and allows to
investigate different, albeit constant, values for $\gamma$, like $2.55$ [12]
or $3.0$ [13]. Another approach, based on Simulated Annealing (SA), included
the branching exponent as an optimization parameter [14]. Results show that
the vascular topology and the metabolic demand significantly influence the
value of the branching exponent. Recently, the authors extended the CCO
approach to finding a synthetic tree optimal both in (global) geometry and
topology [15]. Finding the optimal geometry is cast into a nonlinear
optimization problem (NLP), which allows the investigation of various possible
goal functions and constraints.

In this paper, we utilize this flexibility and go beyond previous studies by
allowing the branching exponent $\gamma$ to vary locally. Furthermore, we
include a blood viscosity law based on the Fåhræus–Lindqvist effect [16] and
enforce equal pressure drop to terminal vessels. The goal is to investigate
the change in branching exponents under these influences. We start by
introducing the relevant definitions and assumptions to generate synthetic
trees. We then cast our goals and constraints into NLPs and introduce our
optimization framework in more detail. Finally, we generate full portal venous
trees of the human liver with up to one million terminal vessels and compare
them against a vascular corrosion cast of a human liver [17, 18].

##  2 Methods

###  2.1 Definitions and assumptions

We represent a vascular tree as a directed branching network
$\mathbb{T}=(\mathbb{V},\mathbb{A})$ with nodes $u\in\mathbb{V}$ and segments
$a\in\mathbb{A}$. Each segment $a=uv$ connects a proximal node $x_{u}$ with a
distal node $x_{v}$. It approximates a vessel as a rigid and straight
cylindrical tube and is defined by its radius $r_{a}$, length $\ell_{a}=\lVert
x_{u}-x_{v}\rVert$, volumetric flow $Q_{a}$ and apparent blood viscosity
$\eta_{a}$. The distal nodes of _terminal segments_ are terminal nodes
(_leaves_) $v\in\mathbb{L}$, and the proximal node of the (single) _root
segment_ is the root node $x_{0}$. A synthetic vascular tree perfuses blood at
a steady state from the root segment down to the terminal segments inside a
given (non-convex) perfusion volume $\Omega\subset\mathbb{R}^{3}$,
schematically shown in Fig. 1

Figure 1: Schematic of a vascular tree and its relation to nodes and segments.
Red circles denote a node, and white rectangles a segment. This tree has a
given inflow $Q_{1}=Q_{\textup{perf}}$ and equal terminal outflow
$Q_{2}=Q_{3}=Q_{\textup{term}}$ through each of the outlets (leaves)

As in Murray’s original paper [1], we assume laminar flow and approximate
blood as an incompressible homogeneous Newtonian fluid. We express the
hydrodynamic resistance $R_{a}$ of segment $a$ by Poiseuille’s law with

$R_{a}=\frac{8\eta_{a}}{\pi}\frac{\ell_{a}}{r_{a}^{4}}\quad\forall
a\in\mathbb{A}.$ (2)

The pressure drop over a segment can now be computed as

$\Delta p_{a}=R_{a}Q_{a}\quad\forall a\in\mathbb{A},$ (3)

and the pressure at a node $v$ follows with

$p_{v}=p_{u}-\Delta p_{a}\quad\forall uv\in\mathbb{A}.$ (4)

We further assume that the (known) perfusion flow $Q_{\textup{perf}}$ is
homogeneously distributed among all $N$ terminal segments, leading to a
terminal flow value $Q_{\textup{term}}=Q_{\textup{perf}}/N$. All remaining
flow values can then be computed using Kirchhoff’s law with
$Q_{uv}=\sum_{vw\in\mathbb{A}}Q_{vw}\forall
v\in\mathbb{V}\setminus\left({0}\cup\mathbb{L}\right)$.

We aim at generating vessels down to the smallest arterioles/venules with
typical radii in the range of $0.015\text{\,}\mathrm{m}\mathrm{m}$ to
$0.1\text{\,}\mathrm{m}\mathrm{m}$. The Fåhræus–Lindqvist effect [19] should
be accounted for at this scale. It describes how the blood viscosity decreases
as the vessel diameter decreases. The tendency of red blood cells to migrate
toward the vessel center is largely responsible for this effect. In turn, this
forces plasma toward the walls and decreases peripheral friction. At the
smallest vessels with radii approaching the radii of red blood cells, the
viscosity sharply rises again. Pries et al. [16] derived an empirical
relationship for this behavior with

$\displaystyle\eta(r_{a})$
$\displaystyle=\eta_{p}\bigl{(}\kappa+\kappa^{2}\bigl{(}\eta_{45}-1\bigr{)}\bigr{)},$
(5) $\displaystyle\eta_{45}$ $\displaystyle=6\exp(-170r_{a}/\text{mm})-$ (6)
$\displaystyle 2.44\exp(-8.08(r_{a}/\text{mm})^{0.645})+3.2,$ (7)
$\displaystyle\kappa$
$\displaystyle=\frac{r_{a}^{2}}{(r_{a}-0.00055\text{mm})^{2}},$ (8)

where $\eta_{p}$ is the viscosity of the plasma, which we set to
$\eta_{p}=$1.125\text{\,}\mathrm{c}\mathrm{P}$$. $\eta_{45}$ is the relative
apparent blood viscosity for a discharge hematocrit of $0.45$. This
relationship is depicted in Fig. 2 for the relevant radii between
$0.015\text{\,}\mathrm{m}\mathrm{m}$ and $10\text{\,}\mathrm{m}\mathrm{m}$.

Figure 2: Change in apparent blood viscosity due to the Fåhræus–Lindqvist
effect as approximated by Pries et al. [16]

###  2.2 Design goals and constraints

####  2.2.1 Murray’s minimization problem

The original minimization formulation by Murray states that the total power of
a vascular tree consists of the metabolic power required to sustain blood
$P_{\textup{vol}}$ and the viscous power $P_{\textup{vis}}$ required to pump
blood from the root down to the micro-circulation. The cost function of a tree
is then defined as

$\displaystyle
f_{\mathbb{T}}=P_{\textup{vol}}+P_{\textup{vis}}=\sum_{a\in\mathbb{A}}m_{b}\pi\ell_{a}r_{a}^{2}+\frac{8\eta_{a}}{\pi}\frac{\ell_{a}}{r_{a}^{4}}Q_{a}^{2},$
(9)

where $m_{b}$ is the metabolic demand of blood in $\mathrm{\SIUnitSymbolMicro
W}\text{\,}{\mathrm{mm}}^{-3}$. As described in [15], we can now include the
nodal positions $x$, the length $\ell$ and radii $r$ as well as the blood
viscosity $\eta$ in the vector of optimization variables $y$, leading to
$y_{1}=(x,\ell,r,\eta)$. We add physical lower bounds $\ell^{-},\ r^{-}$ and
$\eta^{-}$ and, for numerical efficiency, upper bounds $\ell^{+},\ r^{+}$ and
$\eta^{+}$. The best geometry is then found in the rectangle defined as

$\displaystyle
Y_{1}=\mathbb{R}^{3|\mathbb{V}|}\times[\ell^{-},\ell^{+}]^{\mathbb{A}}\times[r^{-},r^{+}]^{\mathbb{A}}\times[\eta^{-},\eta^{+}]^{\mathbb{A}}.$
(10)

Our NLP ” _Power minimization_ ” finally reads:

$\displaystyle\min_{y_{1}\in Y_{1}}\quad$
$\displaystyle\sum_{a\in\mathbb{A}}m_{b}\pi\ell_{a}r_{a}^{2}+8\eta/\pi
Q_{a}^{2}\ell_{a}/r_{a}^{4}$ (11) s.t. $\displaystyle
0=x_{u}-\leavevmode\bar{x}_{u},$ $\displaystyle u$
$\displaystyle\in\mathbb{V}_{0}\cup\mathbb{L}$ (12) $\displaystyle
0=\ell_{uv}^{2}-\lVert x_{u}-x_{v}\rVert^{2},$ $\displaystyle uv$
$\displaystyle\in\mathbb{A}$ (13) $\displaystyle 0=\eta_{a}-\eta(r_{a})$
$\displaystyle a$ $\displaystyle\in\mathbb{A}$ (14)

Eq. 12 fixes the position of terminal nodes and Eq. 13 ensures consistency
between nodal positions and segment length. The third constraint in Eq. 14
enforces the Fåhræus–Lindqvist effect as defined in Eq. 5.

####  2.2.2 Enforcing equal pressure drop

In Murray’s original formulation, no consideration of the resulting pressure
values at terminal segments was given. This terminal pressure is a crucial
parameter for the regulation of blood flow and blood velocity at the
microcirculatory domains. Since we assume these domains are roughly
homogeneous across the organ and have equal demand, the pressure should not
differ significantly. Therefore we enforce equal pressure at each terminal
segment by adding the pressure $p_{v}$ at each node as a new unknown in our
(second) NLP with variables $y_{2}=(x,\ell,r,\eta,p)$, leading to

$\displaystyle
Y_{2}=\mathbb{R}^{4|\mathbb{V}|}\times[\ell^{-},\ell^{+}]^{\mathbb{A}}\times[r^{-},r^{+}]^{\mathbb{A}}\times[\eta^{-},\eta^{+}]^{\mathbb{A}}.$
(15)

Secondly, we constrain the pressure drop between the root and the terminal
nodes as a prescribed constant value $\Delta_{p}$. Since the viscous power at
each segment $a$ is directly proportional to the pressure drop by a factor of
$Q_{a}$, the total viscous power $P_{\textup{vis}}$ becomes a constant. Thus,
we can remove it from the cost function. Finally, we can drop the constant
factor $m_{b}$, leading to a minimization goal proportional to the tree volume
$V_{\mathbb{T}}$. This formulation is used in most synthetic tree studies,
e.g., [11, 12, 20, 14, 15]. Our NLP ” _Volume minimization_ ” then reads:

$\displaystyle\min_{y_{2}\in Y_{2}}\quad$
$\displaystyle\sum_{a\in\mathbb{A}}\pi\ell_{a}r_{a}^{2}$ (16) s.t.
$\displaystyle 0=x_{u}-\leavevmode\bar{x}_{u},$ $\displaystyle u$
$\displaystyle\in\mathbb{V}_{0}\cup\mathbb{L}$ (17) $\displaystyle
0=\ell_{uv}^{2}-\lVert x_{u}-x_{v}\rVert^{2},$ $\displaystyle uv$
$\displaystyle\in\mathbb{A}$ (18) $\displaystyle
0=p_{u}-p_{v}-(8\eta/\pi)Q_{uv}\ell_{uv}/r_{uv}^{4},$ $\displaystyle uv$
$\displaystyle\in\mathbb{A}$ (19) $\displaystyle 0=p_{u},$ $\displaystyle u$
$\displaystyle\in\mathbb{L}$ (20) $\displaystyle 0=p_{0}-\Delta_{p},$ (21)
$\displaystyle 0=\eta_{a}-\eta\left(r_{a}\right)$ $\displaystyle a$
$\displaystyle\in\mathbb{A}$ (22)

Remark 1: Murray’s law, as stated in Eq. 1, is not incorporated explicitly.
Instead, we assume it holds implicitly and compute the associated branching
exponent $\gamma$ using the Newton-Raphson method after the optimization is
finished.

####  2.2.3 Additional optimization variants

To better isolate the individual influence of different factors, we define
additional variants of our two minimization problems. Firstly, we simplify
both problems to a constant apparent viscosity
$\eta_{\textup{const}}=$3.6\text{\,}\mathrm{c}\mathrm{P}$$, removing $\eta$
from the vector of optimization variables and dropping the corresponding
constraints Eq. 14 and Eq. 22. Secondly, we investigate the influence of the
metabolic demand $m_{b}$ and the total pressure drop $\Delta_{p}$. We consider
values between $0.1\text{\,}\mathrm{\SIUnitSymbolMicro
W}\text{\,}{\mathrm{mm}}^{-3}$ and $1.0\text{\,}\mathrm{\SIUnitSymbolMicro
W}\text{\,}{\mathrm{mm}}^{-3}$ for the metabolic demand, but note that
estimates of this parameter vary significantly [21]. For the total pressure
drop, we set the terminal pressure to
$p_{\textup{term}}=$6\text{\,}\mathrm{mmHg}$$ and vary the root pressure
between $10\text{\,}\mathrm{mmHg}$ and $14\text{\,}\mathrm{mmHg}$. Finally,
for each variant, a separate tree is generated, where Murray’s law (see Eq. 1)
is enforced directly with a single exponent $\gamma_{\textup{opt}}$, included
in the optimization variables. All variants include the same root flow
$Q_{\textup{perf}}=$1.1\text{\,}\mathrm{l}\mathrm{/}\mathrm{m}\mathrm{i}\mathrm{n}$$.

###  2.3 Generation framework

We generate our synthetic trees using the framework introduced in [15], which
we summarize in the following. First, we generate $N_{\textup{topo}}$ terminal
nodes on a regular cubic grid inside our organ’s volume. The root position is
manually set and connected to the geometric center of the volume, which in
turn is connected to all terminal nodes. We swap segments to explore new
topologies from this initial (fan-shaped) tree. A _swap_ detaches a segment
from its parent and connects it with another segment. After each swap, the
geometry is optimized by solving the corresponding NLP. The newly created
topology is accepted on the basis of an SA approach. After topology
optimization, we grow the tree using a modified CCO approach. Here, we
optimize the global geometry each time after adding $N_{\textup{geo}}$ new
terminal nodes and then increase $N_{\textup{geo}}$ heuristically based on the
current density of the tree. Notably, we drop the local optimization of
branching positions and set them to their flow-weighted mean, similar to [22].
Due to our repeated global geometry optimization, this simplification had no
significant impact on the final tree structure. In the last step of the
optimization, we delete all segments that reached the lower bound $\ell^{-}$
(_degenerate segments_), possibly creating $n$-furcations $\left(n\geq
3\right)$.

We then classify the hierarchy throughout the finished tree by assigning each
segment an order number corresponding to the Strahler ordering method [23].
Continuous segments of the same Strahler order correspond to one _vessel_.
Additionally, we employ an ordering scheme based on [4] to allow direct
comparison (in reverse order) to the _generation_ notation used for the
vascular corrosion cast in [18].

Remark 2: The complete optimization framework is only applied once to obtain a
common topology. We solve each NLP variant with this topology to get the
corresponding global geometry. We choose this method to focus on the geometry
changes and to allow a direct comparison of branching exponents and radii at
the same branch types.

##  3 Results

###  3.1 Overall structure of the synthetic portal vein

The full synthetic portal vein tree for the case of volume minimization with
variable viscosity is depicted in Fig. 3. The left side shows the complete
tree inside the non-convex liver domain with two zoom levels. The tree splits
into four major branches, which further split into $8$ main branches. These
results align with previous results of a sparser tree in [15].

For a detailed comparison between the different optimization variants
discussed in Section 2.2.3, we summarized the results in Table 1. Here, a
column represents the results of a single variant, with the first four columns
corresponding to the NLP ”Power minimization” and the last four columns
corresponding to the NLP ”Volume minimization”. For ”Power minimization”, we
included the results for metabolic demands $m_{b}$ of
$0.1\text{\,}\mathrm{\SIUnitSymbolMicro W}\text{\,}{\mathrm{mm}}^{-3}$ and
$1.0\text{\,}\mathrm{\SIUnitSymbolMicro W}\text{\,}{\mathrm{mm}}^{-3}$.
Similarly, for ”Volume minimization”, we included the results for root
pressures $p_{\textup{root}}$ of $10\text{\,}\mathrm{mmHg}$ and
$14\text{\,}\mathrm{mmHg}$. The last row indicates the results of each variant
after enforcing a single (optimal) branching exponent $\gamma_{\textup{opt}}$.

For ”Power minimization”, an increase in metabolic demand $m_{b}$ from
$0.1\text{\,}\mathrm{\SIUnitSymbolMicro W}\text{\,}{\mathrm{mm}}^{-3}$ to
$1.0\text{\,}\mathrm{\SIUnitSymbolMicro W}\text{\,}{\mathrm{mm}}^{-3}$ leads
to an increase in viscous power $P_{\textup{vis}}$, shown in the first row of
Table 1, by around $364\text{\,}\mathrm{\char 37\relax}$ and a reduction in
volume $V_{\mathbb{T}}$, shown in the second row Table 1, by around
$53\text{\,}\mathrm{\char 37\relax}$. Similarly, for the ”Volume
minimization”, an increase in the root pressure from
$10\text{\,}\mathrm{mmHg}$ to $14\text{\,}\mathrm{mmHg}$ leads to a
$320\text{\,}\mathrm{\char 37\relax}$ increase in viscous power and a
$51\text{\,}\mathrm{\char 37\relax}$ reduction in volume. The
Fåhræus–Lindqvist effect had a minor influence. It decreased the viscous power
by around $1\text{\,}\mathrm{\char 37\relax}$ and the total volume by around
$4\text{\,}\mathrm{\char 37\relax}$ in all cases.

Figure 3: Complete synthetic vascular tree of the portal vein of a human liver
with $1,000,000$ terminal vessels (volume minimization with variable viscosity
and $p_{\textup{root}}=$10\text{\,}\mathrm{mmHg}$$). Two zoom levels highlight
the hierarchical structure at different scales. The radii are between
$5.1\text{\,}\mathrm{m}\mathrm{m}$ (root vessel) and
$0.017\text{\,}\mathrm{m}\mathrm{m}$ (smallest terminal vessel)

###  3.2 Vessel radii

The root radius $r_{\textup{root}}$ is shown in the third row of Table 1. It
decreased by $31\text{\,}\mathrm{\char 37\relax}$ after the metabolic demand
$m_{b}$ was increased for ”Power minimization”. Similarly, in the case of
”Volume minimization”, it decreased by $30\text{\,}\mathrm{\char 37\relax}$
after the root pressure $p_{\textup{root}}$ was increased to
$14\text{\,}\mathrm{mmHg}$. In all cases, the Fåhræus–Lindqvist effect had a
limited influence on the root radius with changes less than
$0.1\text{\,}\mathrm{\char 37\relax}$. However, a significant decrease in
radius can be observed for vessels between Strahler order $1$ and $6$, shown
in Fig. 4. In both NLP cases, this decrease was highest at the terminal
vessels with $2\text{\,}\mathrm{\char 37\relax}$ for ”Power minimization”
(Fig. 4a) and $1.5\text{\,}\mathrm{\char 37\relax}$ for ”Volume minimization”
(Fig. 4b). A notable difference between both cases is the variance of radii at
each Strahler order. In the case of ”Power minimization”, the highest variance
is observed at Strahler order $6$, whereas the terminal radii are constant. In
the case of ”Volume minimization”, the highest variance is at the terminal
segments and decreases as the Strahler order increases.

Table 1: Comparison of the results between the different variants of our two minimization problems as introduced in Section 2.2.3. The branching exponents $\gamma_{\textup{opt}}$ in the last row (marked with *) are the results of separate runs for each variant, where a single constant branching exponent was enforced. | Power minimization | Volume minimization
---|---|---
| $m_{b}=$0.1\text{\,}\mathrm{\SIUnitSymbolMicro W}\text{\,}{\mathrm{mm}}^{-3}$$ | $m_{b}=$1.0\text{\,}\mathrm{\SIUnitSymbolMicro W}\text{\,}{\mathrm{mm}}^{-3}$$ | $p_{\textup{root}}=$10.0\text{\,}\mathrm{mmHg}$$ |  $p_{\textup{root}}=$14.0\text{\,}\mathrm{mmHg}$$
Parameter | $\eta_{\textup{const}}$ | $\eta(r)$ | $\eta_{\textup{const}}$ | $\eta(r)$ | $\eta_{\textup{const}}$ | $\eta(r)$ | $\eta_{\textup{const}}$ |  $\eta(r)$
$P_{\textup{vis}}$ in $\mathrm{m}\mathrm{W}$ | 2.67 | 2.65 | 12.34 | 12.23 | 3.33 | 3.33 | 13.99 | 13.99
$V_{\mathbb{T}}$ in $\mathrm{m}\mathrm{m}^{3}$ | 53,400.60 | 52,283.94 | 24,786.36 | 24,042.08 | 50,413.84 | 49,158.58 | 24,587.02 | 23,154.13
$r_{\textup{root}}$ in $\mathrm{m}\mathrm{m}$ | 5.35 | 5.35 | 3.64 | 3.63 | 5.10 | 5.09 | 3.56 | 3.54
$p_{\textup{term}}$ in $\mathrm{m}\mathrm{m}\mathrm{H}\mathrm{g}$ | $[10.34,11.68]$ | $[10.41,11.73]$ | $[1.15,10.67]$ | $[1.74,10.98]$ | 6.00 | 6.00 | 6.00 | 6.00
$\gamma$ | 3.00 | $[2.90,3.00]$ | 3.00 | $[2.90,3.00]$ | $[1.75,3.01]$ | $[1.76,3.17]$ | $[1.43,3.00]$ | $[1.46,3.08]$
$\gamma_{\textup{opt}}$ (constant) | 3.00* | 2.91* | 3.00* | 2.92* | 2.84* | 2.76* | 2.82* | 2.74*

(a) Power minimization ($m_{b}=$0.1\text{\,}\mathrm{\SIUnitSymbolMicro
W}\text{\,}{\mathrm{mm}}^{-3}$$)

(b) Volume minimization ($p_{\textup{root}}=$10\text{\,}\mathrm{mmHg}$$)

Figure 4: Influence of Fåhræus–Lindqvist effect on vessel radii for different
Strahler orders

###  3.3 Pressure drop

After ”Power minimization”, the terminal pressures are not constant across the
tree but exhibit a wide range of values, see Table 1 row 4 column 1–4. This
range widens further for higher metabolic demands and also increases the mean
total pressure drop from the root to the terminal segments, shown in Fig. 5.
In Fig. 6a, the pressure values at different Strahler orders are shown for
$m_{b}=$0.1\text{\,}\mathrm{\SIUnitSymbolMicro
W}\text{\,}{\mathrm{mm}}^{-3}$$. Pressure values drop and variances increase
with decreasing Strahler order. Including the Fåhræus–Lindqvist effect leads
to slightly higher pressure values for the Strahler orders $1$ to $6$. The
terminal pressures after ”Volume minimization” are fixed to
$p_{\textup{term}}=$6\text{\,}\mathrm{mmHg}$$ as enforced by Eq. 19 – Eq. 21.
The effect of these constraints is highlighted in Fig. 6b for root pressure
$p_{\textup{root}}=$10\text{\,}\mathrm{mmHg}$$. In contrast to ”Power
minimization”, variances are significantly higher at the intermediate Strahler
orders $3$ to $11$. Furthermore, the influence of the Fåhræus–Lindqvist effect
is more pronounced, decreasing the pressure values between Strahler order $2$
to $10$.

Figure 5: Density plot of the total pressure drop from root vessel to terminal
vessels for the power minimization under different metabolic demand

(a) Power minimization ($m_{b}=$0.1\text{\,}\mathrm{\SIUnitSymbolMicro
W}\text{\,}{\mathrm{mm}}^{-3}$$)

(b) Volume minimization ($p_{\textup{root}}=$10\text{\,}\mathrm{mmHg}$$)

Figure 6: Influence of Fåhræus–Lindqvist effect on pressure values of the
distal nodes of vessels for different Strahler orders

###  3.4 Branching behavior

The resulting branching exponents of all variants are summarized in row 5 of
Table 1. For ”Power minimization” with constant viscosity
$\eta_{\textup{const}}$, the exponents are constant with $\gamma=3.0$ across
all branches regardless of metabolic demand $m_{b}$. In contrast, the
inclusion of the Fåhræus–Lindqvist effect leads to deviations from $3.0$, with
branching exponents reaching a minimum of $2.9$. For ”Volume minimization,”
exponents are not constant even for constant blood viscosity. Instead, most
values fall between $2.0$ and $3.0$, with the smallest outliers having values
of $1.43$.

For a more detailed comparison, the probability density function of the
branching exponents for both optimization cases is compared in Fig. 7. The
influence of the Fåhræus–Lindqvist effect shifts most branching exponents from
a constant $3.0$ to $2.9$ during ”Power minimization” (Fig. 7a). During
”Volume minimization” with constant blood viscosity, most exponents are at
$3.0$ (Fig. 7b in red) and are shifted to $2.9$ when including the
Fåhræus–Lindqvist effect (Fig. 7b in green).

Fig. 8 highlights the distribution of mean branching exponents across
different branch types. Each cell $(i,j)$ corresponds to a branch with child
segments of Strahler order $i$ and $j$. In Fig. 8(a) the effect of variable
blood viscosity on ”Power minimization” is depicted. The branching exponents
decrease if the Strahler order of either child decreases, leading to the
smallest branching exponent of $2.9$ at branches with two terminal segments.
If both child segments have a Strahler order over $8$, the mean branching
exponent is at its maximum of $3.0$. The effect of enforcing equal terminal
pressure is shown in Fig. 8(b). Here, the higher the difference between the
Strahler orders of both children is, the smaller the branching exponent is.
Again, the smallest mean branching exponents are observed at branches
connecting two terminal segments with a value of $2.76$. Fig. 8(c) shows the
accumulated effect of both constraints with a minimum mean exponent of $2.7$,
again at terminal branches.

(a) Power minimization ($m_{b}=$0.1\text{\,}\mathrm{\SIUnitSymbolMicro
W}\text{\,}{\mathrm{mm}}^{-3}$$)

(b) Volume minimization ($p_{\textup{root}}=$10\text{\,}\mathrm{mmHg}$$)

Figure 7: Effect of enforcing variable viscosity and equal pressure onto the
branching exponents Figure 8: Mean values of branching exponent $\gamma$ for
different types of branches, e.g., the cell $(1,2)$ corresponds to branches,
where child vessels of Strahler order $1$ and $2$ meet. The results are
symmetric. White cells correspond to branch types, which do not occur in the
tree topology. (a) Power minimization with variable viscosity, (b) Volume
minimization with constant viscosity, (c) Volume minimization with variable
viscosity

###  3.5 Comparison to vascular corrosion cast

We now directly compare our synthetic trees against a corrosion cast of the
portal vein of the human liver [18]. In Fig. 9a, the radii per generation
(radius-adjusted Strahler order in reverse) are compared between the ”Power
minimization” with $m_{b}=$0.1\text{\,}\mathrm{\SIUnitSymbolMicro
W}\text{\,}{\mathrm{mm}}^{-3}$$, the volume minimization with
$p_{\textup{root}}=$10\text{\,}\mathrm{mmHg}$$ and the corrosion cast data,
including measurements and a best-fit trend line, based on the least sum of
square errors. The synthetic trees’ radii fit the data and trend line well for
the generations $5$ to $15$. Notably, they significantly deviate for the first
four generations, especially against the measurements with errors of around
25%. In contrast, the number of vessels in Fig. 9b of both synthetic trees fit
the data of the corrosion cast well for all generations.

(a) Mean radii

(b) Number of vessels

Figure 9: Comparison of synthetic trees (with variable viscosity) against
corrosion cast data [18]

##  4 Discussion

The exclusion of Murray’s law simplifies our two optimization problems and
automatically allows the branching exponents to vary locally. Both problems
can generate synthetic trees with similar geometry and topology using
appropriate parameters. Power minimization with a metabolic factor of
$m_{b}=$0.1\text{\,}\mathrm{\SIUnitSymbolMicro W}\text{\,}{\mathrm{mm}}^{-3}$$
leads to more realistic radii than a value of
$1.0\text{\,}\mathrm{\SIUnitSymbolMicro W}\text{\,}{\mathrm{mm}}^{-3}$, which
is in line with estimated values around $0.3-0.4$ $\mathrm{\SIUnitSymbolMicro
W}\text{\,}{\mathrm{mm}}^{-3}$ for venous trees [24, 21]. Similarly, a root
pressure of $p_{\textup{root}}=$14\text{\,}\mathrm{mmHg}$$ results in a
pressure drop of $\Delta_{p}=$8\text{\,}\mathrm{mmHg}$$, which is related to
portal hypertension. Such a high pressure drop leads to unrealistic small
radii in comparison to a more realistic root pressure of
$p_{\textup{root}}=$10\text{\,}\mathrm{mmHg}$$.

Under power minimization, the radii are at their individual local optima,
i.e., the radius of each segment can be solved independently as the minimum of
the metabolic demand and the viscous power dissipation. This observation is in
line with the findings of Murray and explains the constant branching exponent
of $3$ in Table 1. The variations in radii for Strahler orders $2$ to $12$ in
Fig. 4a are based entirely on the branching type, which is completely defined
by the flow values of the child segments.

For volume minimization, no such simplification can be made, as the constraint
of equal pressure creates dependencies between segments on the same path to
the root. This constraint also forces radii to deviate from their local
optima, which means a deviation from the branching exponent $\gamma=3.0$. The
highest deviations are found at branches between two terminal segments,
because they can be adjusted to a given pressure drop without significantly
increasing the tree’s volume. Given the same length, they also constitute a
higher pressure drop than segments with bigger radius.

The inclusion of the Fåhræus–Lindqvist effect reduces the viscosity of the
smaller vessels and, in turn, increases their pressure drop. The corresponding
branching exponents also deviate from $\gamma=3.0$, as can be observed in Fig.
8. In contrast, the effect on bigger vessels is negligible and results in
constant exponents $\gamma=3.0$, as seen in Fig. 8(a), for branches where both
children have Strahler orders over $8$.

While the generated trees generally fit the corrosion cast data well, the
underestimation of the largest radii (generation $1$ to $4$) is significant.
Since the vessels of the portal venous tree are more elliptical than circular,
the radii of the vascular corrosion cast were estimated based on their cross-
sectional area. Perhaps the synthetic trees’ radii may correspond better with
estimations based on other criteria such as a maximal inscribed circle.
Another explanation is that the total energy dissipation is likely higher than
our models predict. This underestimation could be due to ignoring the effect
of turbulent flow [8] and simplified geometric modelling of the branching of
vessels. The effect of pulsatile flow, however, can be neglected because the
blood comes not directly from the heart but first flows through the digestive
tract, leaving only a limited amount of pulsatility inside the portal venous
tree.

##  5 Conclusion

Our optimization framework can handle complex constraints and goal functions
while generating synthetic trees up to but not including the capillary level
of the microcirculation. We used our framework to investigate the local
branching behavior for different constraints and goal functions. Branching
exponents automatically lie in the experimentally predicted range between
$2.0$ and $3.0$. Even small changes to Murray’s original optimization problem,
like the inclusion of variable blood viscosity, significantly affect the
optimal branching exponents of vessels. The topology and geometry of our
synthetic trees closely follow the vascular corrosion cast of a human portal
venous tree, with significant deviations only in the largest vessels.

Without enforcing any pressure constraint, terminal pressures vary by up to
$1.0\text{\,}\mathrm{mmHg}$, leading to highly heterogeneous boundary
conditions for the microcirculation. In contrast, enforcing equal terminal
pressure in its current form leads to pressure variations up to
$2.0\text{\,}\mathrm{mmHg}$ in the intermediate vessels of the
mesocirculation. Both these results need to be critically evaluated and
compared against measurements of real vascular trees. One approach to improve
the pressure constraint could be to prescribe a certain physical range using
inequality constraints rather than a fixed value. In the future, we plan to
include pulsatile and turbulent flow effects in a closed form, improving our
framework’s computational ability. Furthermore, shear stress, a critical
parameter for vascular growth, must also be incorporated into the model.

A more mature version of this model, specifically the ability to predict
optimal branching exponents under different constraints, could have many
potential applications in the medical field. An example would be to predict
and relate the branching behavior across the scales to vascular diseases.
These predictions could improve the interpretation of medical images by giving
valuable input to the functional assessment of organs.

## Acknowledgment

The results presented in this paper were obtained as part of the ERC Starting
Grant project “ImageToSim” that has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (Grant agreement No. 759001). The authors gratefully
acknowledge this support.

## References

  * [1] C. D. Murray, “The physiological principle of minimum work: I. The vascular system and the cost of blood volume,” Proceedings of the National Academy of Sciences of the United States of America, vol. 12, no. 3, p. 207, 1926.
  * [2] C. D. Murray, “The physiological principle of minimum work applied to the angle of branching of arteries,” The Journal of General Physiology, vol. 9, no. 6, pp. 835–841, 1926.
  * [3] K. Horsfield and M. J. Woldenberg, “Diameters and cross-sectional areas of branches in the human pulmonary arterial tree,” The Anatomical Record, vol. 223, no. 3, pp. 245–251, 1989.
  * [4] Z. Jiang, G. Kassab, and Y. Fung, “Diameter-defined Strahler system and connectivity matrix of the pulmonary arterial tree,” Journal of Applied Physiology, vol. 76, no. 2, pp. 882–892, 1994.
  * [5] E. VanBavel and J. Spaan, “Branching patterns in the porcine coronary arterial tree. Estimation of flow heterogeneity.,” Circulation Research, vol. 71, no. 5, pp. 1200–1212, 1992.
  * [6] Y. Nakamura and S. Awa, “Radius exponent in elastic and rigid arterial models optimized by the least energy principle,” Physiological reports, vol. 2, no. 2, p. e00236, 2014.
  * [7] H. Uylings, “Optimization of diameters and bifurcation angles in lung and vascular tree structures,” Bulletin of mathematical biology, vol. 39, no. 5, pp. 509–520, 1977.
  * [8] P. R. Painter, P. Edén, and H.-U. Bengtsson, “Pulsatile blood flow, shear force, energy dissipation and murray’s law,” Theoretical Biology and Medical Modelling, vol. 3, no. 1, pp. 1–10, 2006.
  * [9] Y. Zhou, G. S. Kassab, and S. Molloi, “On the design of the coronary arterial tree: A generalization of Murray’s law,” Physics in Medicine & Biology, vol. 44, no. 12, p. 2929, 1999.
  * [10] G. S. Kassab, “Scaling laws of vascular trees: of form and function,” American Journal of Physiology-Heart and Circulatory Physiology, vol. 290, no. 2, pp. H894–H903, 2006.
  * [11] W. Schreiner and P. F. Buxbaum, “Computer-optimization of vascular trees,” IEEE Transactions on Biomedical Engineering, vol. 40, no. 5, pp. 482–491, 1993.
  * [12] R. Karch, F. Neumann, M. Neumann, and W. Schreiner, “A three-dimensional model for arterial tree representation, generated by constrained constructive optimization,” Computers in Biology and Medicine, vol. 29, no. 1, pp. 19–38, 1999.
  * [13] W. Schreiner, R. Karch, M. Neumann, F. Neumann, S. M. Roedler, and G. Heinze, “Heterogeneous perfusion is a consequence of uniform shear stress in optimized arterial tree models,” Journal of theoretical biology, vol. 220, no. 3, pp. 285–301, 2003.
  * [14] J. Keelan and J. P. Hague, “The role of vascular complexity on optimal junction exponents,” Scientific Reports, vol. 11, no. 1, pp. 1–14, 2021\.
  * [15] E. Jessen, M. C. Steinbach, C. Debbaut, and D. Schillinger, “Rigorous mathematical optimization of synthetic hepatic vascular trees,” Journal of the Royal Society Interface, vol. 19, no. 191, p. 20220087, 2022.
  * [16] A. Pries, T. Secomb, T. Gessner, M. Sperandio, J. Gross, and P. Gaehtgens, “Resistance to blood flow in microvessels in vivo.,” Circulation Research, vol. 75, no. 5, pp. 904–915, 1994.
  * [17] C. Debbaut, D. Monbaliu, C. Casteleyn, P. Cornillie, D. Van Loo, B. Masschaele, J. Pirenne, P. Simoens, L. Van Hoorebeke, and P. Segers, “From vascular corrosion cast to electrical analog model for the study of human liver hemodynamics and perfusion,” IEEE Transactions on Biomedical Engineering, vol. 58, no. 1, pp. 25–35, 2010.
  * [18] C. Debbaut, P. Segers, P. Cornillie, C. Casteleyn, M. Dierick, W. Laleman, and D. Monbaliu, “Analyzing the human liver vascular architecture by combining vascular corrosion casting and micro-CT scanning: A feasibility study,” Journal of Anatomy, vol. 224, no. 4, pp. 509–517, 2014.
  * [19] R. Fahraeus and T. Lindqvist, “The viscosity of the blood in narrow capillary tubes,” American Journal of Physiology-Legacy Content, vol. 96, no. 3, pp. 562–568, 1931.
  * [20] M. Kretowski, Y. Rolland, J. Bézy-Wendling, and J.-L. Coatrieux, “Physiologically based modeling of 3-D vascular networks and CT scan angiography,” IEEE Transactions on Medical Imaging, vol. 22, no. 2, pp. 248–257, 2003.
  * [21] Y. Liu and G. S. Kassab, “Vascular metabolic dissipation in murray’s law,” American Journal of Physiology-Heart and Circulatory Physiology, vol. 292, no. 3, pp. H1336–H1339, 2007.
  * [22] A. A. Guy, A. W. Justin, D. M. Aguilar-Garza, and A. E. Markaki, “3d printable vascular networks generated by accelerated constrained constructive optimization for tissue engineering,” IEEE Transactions on Biomedical Engineering, vol. 67, no. 6, pp. 1650–1663, 2019.
  * [23] A. N. Strahler, “Quantitative analysis of watershed geomorphology,” Eos, Transactions American Geophysical Union, vol. 38, no. 6, pp. 913–920, 1957.
  * [24] L. A. Taber, “An optimization principle for vascular radius including the effects of smooth muscle tone,” Biophysical journal, vol. 74, no. 1, pp. 109–114, 1998.