# Bath Engineering Enhanced Quantum Critical Engines

Revathy B. S  Department of Physics, Indian Institute of Technology Palakkad,
Palakkad, 678557, India  Victor Mukherjee  Department of Physical Sciences,
IISER Berhampur, Berhampur 760010, India  Uma Divakaran  Department of
Physics, Indian Institute of Technology Palakkad, Palakkad, 678557, India

###### Abstract

Driving a quantum system across quantum critical points leads to non-adiabatic
excitations in the system. This in turn may adversely affect the functioning
of a quantum machine which uses a quantum critical substance as its working
medium. Here we propose a bath-engineered quantum engine (BEQE), in which we
use the Kibble–Zurek mechanism and critical scaling laws to formulate a
protocol for enhancing the performance of finite-time quantum engines
operating close to quantum phase transitions. In the case of free fermionic
systems, BEQE enables finite-time engines to outperform engines operating in
the presence of shortcuts to adiabaticity, and even infinite-time engines
under suitable conditions, thus showing the remarkable advantages offered by
this technique. Open questions remain regarding the use of BEQE based on non-
integrable models.

##  I Introduction

The field of quantum thermodynamics aims to form a coherent understanding of
the thermodynamics of quantum systems Alicki and Kosloff (2018); Vinjanampathy
and Anders (2016); Gemmer  _et al._ (2009); Kosloff (2013); Binder  _et al._
(2018); Bhattacharjee and Dutta (2021); Myers  _et al._ (2022). In classical
thermodynamics, one can then use this knowledge to understand the limitations
on the performance of quantum machines. In this respect, quantum control can
play a significant role in enabling us to go beyond these limitations and
develop high-performing quantum machines Erdman and Noé (2022); Khait  _et
al._ (2022). This can be especially significant in the case of finite-time
quantum machines Quan  _et al._ (2007); Mukherjee and Divakaran (2021), as
non-adiabatic excitations can be detrimental to the performances of such
machines, thus necessitating the application of controls in order to boost
their outputs Gelbwaser-Klimovsky  _et al._ (2015); Mukherjee  _et al._
(2019).

Control techniques such as shortcuts to adiabaticity (STA) have been shown to
be highly successful in enhancing the output of finite-time quantum engines
del Campo  _et al._ (2014); Hartmann  _et al._ (2020); Deng  _et al._ (2013);
Beau  _et al._ (2016); del Campo  _et al._ (2018); Sels and Polkovnikov
(2017). However, the application of STA can be highly non-trivial in many-body
quantum engines, owing to the diverging dimensions of the associated Hilbert
spaces. This can be especially challenging in quantum engines operating close
to quantum critical points, where the diverging length and time scales can
demand STA protocols involving long-range interactions del Campo  _et al._
(2012); Kolodrubetz  _et al._ (2017). The above challenges motivated us to
search for control protocols beyond STA for application in quantum engines
operating close to quantum phase transitions.

In this work we propose a control protocol aimed at enhancing the efficiency
as well as the output work of quantum engines based on free fermionic working
mediums (WMs) operating close to quantum critical points Sachdev (2011).
Quantum phase transitions have proven to be beneficial for quantum heat
engines Campisi and Fazio (2016); Chen  _et al._ (2019); Ma  _et al._ (2017);
Piccitto  _et al._ (2022); Fogarty and Busch (2020). The universality in
quantum critical machines arising from the Kibble–Zurek mechanism (KZM) has
already been studied in Revathy  _et al._ (2020). Here, we construct a quantum
heat engine using a working medium that undergoes quantum phase transition.
The formation of excitations close to the critical point due to the divergence
of relaxation time results in the loss of adiabaticity, thus reducing the
performance of the quantum machine Dutta  _et al._ (2015); Polkovnikov  _et
al._ (2011); Dziarmaga (2010). Although conventional control techniques such
as STA involve complex calculations and non-trivial many-body interactions, we
propose the implementation of the bath-engineered quantum engine (BEQE), in
which the working of the engine can be improved significantly through the
simple control of bath spectral functions.

The present work is organized as follows. We describe a many-body quantum Otto
cycle in Section II. In Section III, the operation of the BEQE using a generic
free fermionic WM is explained in detail. We study a specific example of the
BEQE using the transverse Ising model in Section IV. Finally, we summarize our
results in Section V. Details of the calculations presented in this work are
included in the Appendix.

##  II Many-Body Quantum Otto Cycle

We consider an Otto cycle with the working medium (WM) described by the
Hamiltonian $H(\lambda(t))$, where $\lambda$ is a time-dependent parameter.
The four-stroke quantum Otto cycle consists of two non-unitary strokes and two
unitary strokes, as described below (Figure 1).

Figure 1: Schematic diagram of a quantum Otto cycle with a many-body spin
system as the working medium.

  * (i)

Non-unitary stroke A $\rightarrow$ B: The WM with parameter
$\lambda=\lambda_{1}$ is connected to an energizing bath $\mathcal{E}$ until
it reaches the corresponding steady state at B by receiving energy
$\mathcal{Q}_{in}$ from the bath.

  * (ii)

Unitary stroke B $\rightarrow$ C: The WM is decoupled from the energizing bath
and $\lambda$ is changed from $\lambda_{1}$ to $\lambda_{2}$ at a speed of
$1/\tau_{1}$.

This unitary evolution is described by the Von Neumann equation of motion:

$\frac{d\rho}{dt}=-i[H,\rho].$ (1)

  * (iii)

Non-unitary stroke C $\rightarrow$ D: The WM with $\lambda=\lambda_{2}$ is now
connected to a decaying bath $\mathcal{D}$ until it reaches the corresponding
steady state at D; energy $\mathcal{Q}_{out}$ flows from the WM to the bath
during this stroke.

  * (iv)

Unitary stroke D $\rightarrow$ A: After decoupling from the decaying bath, the
parameter $\lambda$ is changed back to $\lambda_{1}$ from $\lambda_{2}$ with a
speed of $1/\tau_{2}$.

The WM crosses a quantum critical point at $\lambda=\lambda_{c}$ during the
unitary strokes, such that $\lambda_{2}\leq\lambda_{c}<\lambda_{1}$. The
energy at the end of stroke $i$ is calculated using the equation

${E}_{i}=\text{Tr}(H^{i}\rho^{i})$ (2)

where $H^{i}$ and $\rho^{i}$ are the Hamiltonian and the density matrix at
$i={{A}},{{B}},{{C}},{{D}}$. The heat input ($\mathcal{Q}_{in}$) and heat
output ($\mathcal{Q}_{out}$) can be calculated using

$\displaystyle\mathcal{Q}_{in}$ $\displaystyle=$
$\displaystyle{E}_{B}-{E}_{A}$ (3) $\displaystyle\mathcal{Q}_{out}$
$\displaystyle=$ $\displaystyle{E}_{D}-{E}_{C}.$ (4)

The output work is given by
$\mathcal{W}=-(\mathcal{Q}_{in}+\mathcal{Q}_{out})$. The sign convention used
here is as follows: energy is taken to be positive (negative) if it enters
(leaves) the WM. The Otto cycle works as an engine when
$\mathcal{Q}_{in}>0,\mathcal{Q}_{out}<0$ and $\mathcal{|W|}<0$. The
performance of the engine is then characterized using the quantity of
efficiency ($\eta$), which is defined as

$\displaystyle\eta$ $\displaystyle=$
$\displaystyle-\frac{\mathcal{W}}{\mathcal{Q}_{in}}.$ (5)

We note that other regimes of operation may arise for different signs of
$\mathcal{Q}_{in}$ and $\mathcal{Q}_{out}$, as discussed in Mukherjee  _et
al._ (2016).

##  III Bath Spectral Form Engineering

We consider a free fermionic WM, described by a Hamiltonian of the form

$\displaystyle H=\sum_{k}\psi_{k}^{\dagger}H_{k}\psi_{k}$ $\displaystyle
H_{k}=\vec{f}(k).\vec{\sigma}_{k}$ (6)

where
$\vec{\sigma}_{k}=\left(\sigma_{k}^{x},\sigma_{k}^{y},\sigma_{k}^{z}\right)$
denotes the Pauli matrices corresponding to the $k$-th mode; $\vec{f}(k)$ is a
model-dependent function for the $k$-th mode; and
$\psi_{k}^{\dagger}=\left(c_{1k}^{\dagger}~{}c_{2k}^{\dagger}\right)$, where
$c_{jk}$ and $c_{jk}^{\dagger}$ (with $j=1,2$) denote the fermionic operators
corresponding to the $k$-th mode.

For non-interacting $k$ modes, the density matrix $\rho$ of the system can be
written as $\rho=\otimes_{k}\rho_{k}$. The WM undergoes unitary dynamics
during the strokes D $\to$ A and B $\to$ C, described by the Von Neumann
equation:

$\displaystyle\dot{\rho}_{k}=-i\left[H_{k},\rho_{k}\right]$ (7)

for each $k$ mode. Furthermore, we assume fermionic baths such that each $k$
mode evolves independently during the non-unitary strokes, described by the
master equation Keck  _et al._ (2017)

$\displaystyle\frac{d\rho_{k}}{dt}=\mathcal{G}_{\alpha}(\Delta_{k})\mathcal{L}_{k}\left[\rho_{k}(t)\right]+\mathcal{G}_{\alpha}(-\Delta_{k})\mathcal{L}_{k}^{\dagger}\left[\rho_{k}(t)\right]$
(8)

where, following the Kubo–Martin–Schwinger condition, we have

$\displaystyle\mathcal{G}_{\alpha}(-\Delta_{k})=\exp(-\Delta_{k}/T_{\alpha})\mathcal{G}_{\alpha}(\Delta_{k}).$
(9)

Here $\mathcal{G}_{\alpha}(\nu)$ denotes the spectral function of the
$\alpha=\mathcal{E},\mathcal{D}$ bath at frequency $\nu$, whereas $T_{\alpha}$
is the effective temperature of the $\alpha$ bath Breuer and Petruccione
(2002); Revathy  _et al._ (2020). The superoperator $\mathcal{L}_{k}$ and
$\mathcal{L}_{k}^{\dagger}$ are of the form

$\displaystyle\mathcal{L}_{k}$ $\displaystyle=$
$\displaystyle\left(L_{k}\rho_{k}L_{k}^{\dagger}-\frac{1}{2}\\{L_{k}^{\dagger}L_{k},\rho_{k}\\}\right)$
(10) $\displaystyle\mathcal{L}_{k}^{\dagger}$ $\displaystyle=$
$\displaystyle\left(L_{k}^{\dagger}\rho_{k}L_{k}-\frac{1}{2}\\{L_{k}L_{k}^{\dagger},\rho_{k}\\}\right)$

with $L_{k}$ being the Lindblad operators denoting jumps between the different
eigenenergy levels. The above dynamics given in Equations (8)–(10) ensures
that each $k$ mode thermalizes independently with the bath, such that the
steady state of the WM at the end of an isochoric stroke is given by
$\rho=\otimes_{k}\rho_{k}^{\rm th}$, where $\rho_{k}^{\rm th}$ is the Gibbs
state corresponding to the $k$-th mode Deng  _et al._ (2011). However, in
general, $\Delta_{k}$ can be expected to vary with $k$, which may result in
the global state $\rho=\otimes_{k}\rho_{k}^{\rm th}$ being non-thermal. We
emphasize that even though globally the WM remains in a non-thermal state, in
contrast to quantum engines powered by squeezed thermal baths, none of the $k$
modes receive any ergotropy from the $\mathcal{E}$ and $\mathcal{D}$ baths in
this setup Roßnagel  _et al._ (2014); Niedenzu  _et al._ (2018). Furthermore,
this global steady state becomes thermal for $\Delta_{k}$ becoming a positive
$k$-independent constant, which may happen far away from a quantum critical
point, or in the limits $T_{\alpha}\to\infty$ and $T_{\alpha}\to 0$. In
addition, we note that the bath considered here is local in $k$ space, and
therefore can be expected to be non-local in real space. However, this bath
can become local far away from the critical point, in which regime the system
may be composed of effectively non-interacting particles (see Section IV). A
detailed discussion regarding the bath considered here is given in Keck  _et
al._ (2017).

Non-adiabatic excitations are inevitable when a quantum system is driven
across quantum critical points Dutta  _et al._ (2015). This results in a
reduction in output work as well as efficiency when a quantum critical
substance is made the working medium of a quantum Otto engine Revathy  _et
al._ (2020). Here, we propose bath spectral form engineering to prevent these
excitations from reducing the performance of the engine, henceforth called the
bath-engineered quantum engine (BEQE).

In the bath engineering technique, we choose bath spectral functions
$\mathcal{G}_{\alpha}$ with appropriate cut-offs, such that the modes which
have higher probabilities of getting excited, and which are therefore
detrimental to the performance of the finite-time quantum engine, are not
allowed to participate in the dynamics. Although techniques such as shortcuts
to adiabaticity are applied in the unitary strokes, bath engineering is
performed during the non-unitary strokes (Figure 2a).

Figure 2: (a) Schematic diagram showing bath engineering being applied during
the non-unitary strokes, whereas shortcuts to adiabaticity are applied during
the unitary strokes. (b) A $k$ mode is coupled to the bath if
$\Delta_{k}>\Delta^{*}$ ($k>k^{*}$) and is not coupled to the bath if
$\Delta_{k}<\Delta^{*}$ ($k<k^{*}$), where $\Delta_{k^{*}}=\Delta^{*}$.

The physics of excitations generated in a system which is driven at a finite
rate across a quantum critical point (QCP) are well established and are
described by the Kibble–Zurek mechanism (KZM) Zurek  _et al._ (2005);
Polkovnikov (2005); Dziarmaga (2010), which was also experimentally
demonstrated in Cui  _et al._ (2020); Bando  _et al._ (2020). According to the
adiabatic-impulse approximation Damski and Zurek (2006), these excitations
occur due to vanishing energy gaps which are below a threshold value (say,
$\Delta^{*}$), the expression of which can be obtained using KZ arguments as
described below. Bath engineering is carried out such that the energy levels
having gaps less than $\Delta^{*}$ are not allowed to interact with the bath,
i.e, $\mathcal{G}(\Delta_{k}<\Delta^{*})\approx 0$, thus preventing them from
participating in the operation of the cycle (Figure 2b).

### Kibble–Zurek-Mechanism-Assisted BEQE

According to KZM, the response of a system driven across a quantum critical
point is determined by the inherent time scale (relaxation time $\xi_{\tau}$)
of the system, and the rate of change of the system Hamiltonian Polkovnikov
(2005); Zurek  _et al._ (2005); Damski and Zurek (2006). When the relaxation
time $\xi_{\tau}$ of the system is greater than the rate at which the
Hamiltonian parameter $\lambda$ is changed, the system stops evolving
adiabatically, thus resulting in non-adiabatic excitations. In order to arrive
at a more quantitative analysis, let us assume that $t^{*}$ is the time at
which the system looses adiabaticity and excitations begin to occur. The
energy gap $\Delta_{k_{c}}$ at the critical mode $k_{c}$ scales with the
distance from the critical point $\lambda_{c}$ as Sachdev (2011)

$\Delta_{k_{c}}\sim|\lambda-\lambda_{c}|^{\nu z},$ (11)

where $\nu$ and $z$ are the correlation length and dynamical critical
exponents, respectively. When the parameter $\lambda$ is varied using the
quench protocol $\lambda=\lambda_{2}+(\lambda_{1}-\lambda_{2})\frac{t}{\tau}$,
one can write

$\lambda-\lambda_{c}=\lambda_{2}-\lambda_{c}+(\lambda_{1}-\lambda_{2})\frac{t}{\tau}.$
(12)

According to the adiabatic-impulse approximation Damski and Zurek (2006), the
time $t^{*}$ is determined by the condition that the relaxation time
$\xi_{\tau}$ is of the order of the time scale with which $\lambda$ is
changed, i.e.,

$\displaystyle\frac{\Delta_{k_{c}}}{\dot{\Delta}_{k_{c}}}|_{t=t^{*}}\sim\xi_{\tau}.$
(13)

Furthermore, the relaxation time diverges according to the scaling

$\xi_{\tau}\sim\frac{1}{\Delta_{k_{c}}}\sim|\lambda-\lambda_{c}|^{-\nu z}.$
(14)

Using the expressions (11)–(14) one obtains

$\displaystyle\frac{\Delta_{k_{c}}}{\dot{\Delta}_{k_{c}}}|_{t=t^{*}}$
$\displaystyle\sim$
$\displaystyle\frac{(\lambda_{2}-\lambda_{c})+(\lambda_{1}-\lambda_{2})\frac{t^{*}}{\tau}}{\nu
z(\frac{\lambda_{1}-\lambda_{2}}{\tau})}$ (15) $\displaystyle\sim$
$\displaystyle[(\lambda_{2}-\lambda_{c})+(\lambda_{1}+\lambda_{2})\frac{t^{*}}{\tau}]^{-\nu
z}$ $\displaystyle\Rightarrow t^{*}$ $\displaystyle\sim$ $\displaystyle
t_{c}+\frac{\tau}{\lambda_{1}-\lambda_{2}}\left(\nu
z(\frac{\lambda_{1}-\lambda_{2}}{\tau})\right)^{\frac{1}{1+\nu z}}$ (16)

where
$t_{c}=\tau(\lambda_{c}-\lambda_{2})/\left(\lambda_{1}-\lambda_{2}\right)$ is
the time such that $\lambda(t_{c})=\lambda_{c}$, and we have assumed that
$\left(\lambda-\lambda_{c}\right)>0$ for simplicity.

Thus, the energy gap at which the excitations begin to happen for the critical
mode is given by

$\tilde{\Delta}^{*}=\Delta_{k_{c}}|_{t^{*}}\sim\left(\lambda(t^{*})-\lambda_{c}\right)^{\nu
z}\sim\left(\frac{\nu z(\lambda_{1}-\lambda_{2})}{\tau}\right)^{\frac{\nu
z}{1+\nu z}}.$ (17)

In the quantum Otto cycle, bath engineering is implemented during the non-
unitary strokes ${C}\to{D}$ and ${A}\to{B}$ by choosing

$\displaystyle\mathcal{G}_{\mathcal{D}}(\Delta_{k})$ $\displaystyle\approx 0$
$\displaystyle~{}~{}\text{for}~{}~{}\Delta_{k}<\Delta^{*}$
$\displaystyle\mathcal{G}_{\mathcal{E}}(\Delta_{k})$ $\displaystyle\approx 0$
$\displaystyle~{}~{}\text{for}~{}~{}\Delta_{k}<\gamma\Delta^{*}$
$\displaystyle\Delta^{*}$ $\displaystyle=$ $\displaystyle
C_{1}\tilde{\Delta}^{*}$ (18)

respectively, such that small energy gaps which have a higher probability of
getting excited do not participate in the dynamics. The scaling parameter
$C_{1}$ ($C_{1}>0$), along with Equations (17) and (III), determine the lower
cut-offs for the bath spectral functions; one can choose an appropriate
$C_{1}$ depending on the details of the setup and the constraints involved in
order to improve the performance of an engine. In the numerical results given
below, we have chosen $C_{1}=1$ for simplicity; $\gamma$ is the scaling factor
by which a typical energy gap changes in the ${D}\to{A}$ stroke, and we have
assumed that $\lambda_{2}$ is close to the quantum critical point, i.e.,

$\displaystyle|\lambda_{2}-\lambda_{c}|^{\nu z}\ll\tilde{\Delta}^{*}.$ (19)

For non-critical $\lambda_{2}$ (i.e., $|\lambda_{2}-\lambda_{c}|^{\nu
z}\gg\tilde{\Delta}^{*}$), the energy gaps of the system at $\lambda_{2}$ are
of the order of

$\displaystyle\Delta_{k}\approx C_{2}\left|\lambda_{2}-\lambda_{c}\right|^{\nu
z}+f\left(k,h_{2}\right),$ (20)

where $C_{2}$ is a model-dependent constant related to the minimum energy gap
of the system, whereas $f(k,h_{2})$ is a model-dependent function for the mode
$k$. For low-energy modes, one can expect $|f(k)|\ll
C_{2}\left|\lambda_{2}-\lambda_{c}\right|^{\nu z}$ Sachdev (2011).
Consequently, in this case we take

$\displaystyle\mathcal{G}_{\mathcal{D}}(\Delta_{k})$ $\displaystyle\approx
0~{}~{}$ $\displaystyle\text{for}~{}~{}\Delta_{k}<\Delta^{*}$
$\displaystyle\mathcal{G}_{\mathcal{E}}(\Delta_{k})$ $\displaystyle\approx
0~{}~{}$ $\displaystyle\text{for}~{}~{}\Delta_{k}<\gamma\Delta^{*}$
$\displaystyle\Delta^{*}$ $\displaystyle=$ $\displaystyle
C_{2}\left|\lambda_{2}-\lambda_{c}\right|^{\nu z}+C_{3}.$ (21)

As before, $C_{3}$ ($|C_{3}|\ll C_{2}\left|\lambda_{2}-\lambda_{c}\right|^{\nu
z}$) is a constant which we choose depending on the details of the WM and
constraints on the bath spectral functions. We note that ideally one should
consider $C_{3}$ to be a function of $\tau$; however, in contrast to Equations
(17) and (III), here we consider a $\tau$-independent $\Delta^{*}$ since
$C_{3}$ can be considered to be a small correction over the first term
$C_{2}\left|\lambda_{2}-\lambda_{c}\right|^{\nu z}$ (see Equation (20) and the
text below).

In this control protocol, the bath spectral functions of the modes with large
$\Delta_{k}$ (see Equations (III) and (21)), and therefore the thermalization
times for these modes, remain unchanged and finite. On the other hand, the
modes with small $\Delta_{k}$ do not evolve during the non-unitary strokes.
Consequently, the durations of the non-unitary strokes of a BEQE, and in turn
the total cycle period, remain the same as that of a finite-time engine
without controls. Furthermore, only the modes with large $\Delta_{k}$ values
go to their respective steady states at the end of a non-unitary stroke in a
BEQE, thereby in general giving rise to non-thermal global steady states at
$B$ and $D$.

Next we demonstrate the bath engineering technique using a free fermionic
model, which is described in the following section.

##  IV BEQE with Transverse Ising Model as a WM

A prototypical example of a free fermionic system undergoing quantum phase
transition is the one-dimensional transverse Ising model (TIM). It is an
exactly solvable model and is thus widely studied. The Hamiltonian of the
transverse Ising model is

$H(t)=-J\sum_{n}\sigma_{n}^{x}\sigma_{n+1}^{x}-h(t)\sum_{n}\sigma_{n}^{z}$
(22)

where $J$ is the nearest neighbor interaction strength; $h(t)$ is the
transverse field which is time-dependent, playing the role of $\lambda$ in the
previous section; and $n$ is the lattice site index. Here, $\sigma_{n}^{i}$
with $i=x,y,z$ are the Pauli matrices at each site $n$. This system shows a
zero temperature quantum phase transition from a paramagnetic to ferromagnetic
state at the quantum critical point $h=\pm J$ Lieb  _et al._ (1961); Pfeuty
(1970); Bunder and McKenzie (1999). We set $J=1$ throughout the paper so that
$h=\pm 1$ are the critical points.

After performing Jordan–Wigner fermionization and taking the Fourier
transform, the Hamiltonian $H_{k}$ takes the form Pfeuty (1970):

$\displaystyle H_{k}=-2\left(h-\cos k\right)\sigma^{z}+2\sin k\sigma^{x}.$
(23)

Even though unitary dynamics allows transitions only between $|0\rangle$ and
$|k,-k\rangle=c_{k}^{\dagger}c_{-k}^{\dagger}|0\rangle$, the system bath
interactions lead to transitions to the $|\pm k\rangle=c_{\pm
k}^{\dagger}|0\rangle$ states as well, resulting in the mixing of states Keck
_et al._ (2017); Bandyopadhyay  _et al._ (2018). Therefore, the Hamiltonian is
rewritten in the basis $|0\rangle,|k\rangle,|-k\rangle,|k,-k\rangle$ as

$H_{k}=\begin{bmatrix}-2(h-\cos k)&0&0&2\sin k\\\ 0&0&0&0\\\ 0&0&0&0\\\ 2\sin
k&0&0&2(h-\cos k)\end{bmatrix}$ (24)

with eigenenergies $-\epsilon_{k},0,0,\epsilon_{k}$ where
$\epsilon_{k}=2\sqrt{(h-\cos k)^{2}+{\sin k}^{2}}$.

We now focus on the strokes of the Otto cycle with the TIM as the WM. The
density matrix at B is given by $\rho^{B}$ ($=\otimes_{k}\rho_{k}^{B}$), where

$\rho_{k}^{B}=\begin{bmatrix}\frac{e^{\beta\epsilon_{k}}}{Z_{k}}&0&0&0\\\
0&\frac{1}{Z_{k}}&0&0\\\ 0&0&\frac{1}{Z_{k}}&0\\\
0&0&0&\frac{e^{-\beta\epsilon_{k}}}{Z_{k}}\end{bmatrix}$ (25)

is the thermal state for the mode $k$ corresponding to $T=T_{H}$ and
$h=h_{1}$.

Here $\beta=\frac{1}{k_{B}T_{H}}$ ($k_{B}$ is set to unity for the rest of the
paper) and $Z_{k}=2+e^{\beta\epsilon_{k}}+e^{-\beta\epsilon_{k}}$ are the
partition functions for each $k$ mode. In the unitary stroke (B $\rightarrow$
C), the transverse field is changed from $h_{1}$ to $h_{2}$ according to the
protocol,

$h(t)=h_{1}+(h_{2}-h_{1})(\frac{t}{\tau_{1}}),\hskip
28.45274ptt\in[0,\tau_{1}]$ (26)

in a time $\tau_{1}$ with $h_{1}\gg h_{2}$. During the non-unitary stroke (C
$\rightarrow$ D) the system again reaches a state
$\rho^{D}=\otimes_{k}\rho_{k}^{D}$, where $\rho_{k}^{D}$ is the thermal state
for the mode $k$ corresponding to $T_{C}$ and $h=h_{2}$ at D. The transverse
field $h_{2}$ is then changed back to $h_{1}$ using the same quench protocol
but in time $\tau_{2}$ in the unitary stroke D $\rightarrow$ A.

Now let us examine how bath engineering is implemented in TIM. We first focus
on the case of $h_{2}\to 1$. As discussed before, we make use of selective
coupling between the bath and the working medium so that some $k$ modes close
to the critical mode $k_{c}$, having an energy gap $\Delta_{k}$ lower than the
threshold value $\Delta{{}^{*}}$ ($=\left(\frac{\nu
z(h_{1}-h_{2})}{\tau}\right)^{\frac{\nu z}{1+\nu z}}$), are prohibited from
interacting with the bath, thereby preventing these modes from thermalizing.

The energy gaps between the adjacent non-degenerate eigenstates of the
Hamiltonian (24) are given by

$\displaystyle\Delta_{k}=\epsilon_{k}=2\sqrt{(h-\cos k)^{2}+\sin^{2}k}.$ (27)

For TIM, the critical exponents are $\nu=1$ and $z=1$ so that $\Delta^{*}$, as
obtained in Equations (17) and (III), is given by

$\Delta^{*}=\sqrt{\frac{h_{1}-h_{2}}{\tau_{2}}}$ (28)

Below we present the steps required to incorporate bath engineering into the
quantum Otto cycle.

  * (a)

As discussed above, we consider a lower cut-off for the decaying bath spectral
function, given by $\mathcal{G}_{\mathcal{D}}\left(\Delta\right)\approx 0$ for
$\Delta_{k}<\Delta^{*}$ (see Equations (III) and (28)). This choice of bath
spectral function ensures that modes with $\Delta_{k}<\Delta^{*}$ are not
allowed to interact with the decaying bath, so that
$\rho_{k}^{D}=\rho_{k}^{C}$ for these modes.

On the other hand, modes with $\Delta_{k}>\Delta^{*}$ thermalize with the
decaying bath and reach the state
$\rho_{k}^{D}=\frac{e^{-\beta_{C}H_{k}(h_{2})}}{Z_{k}}$ at ${D}$.

  * (b)

In the D $\rightarrow$ A stroke, the Hamiltonian is changed from $h_{2}$ to
$h_{1}$, starting from the state $\rho_{k}^{D}$ to reach $\rho_{k}^{A}$.

  * (c)

At A, the lower cutoff for the energizing bath is chosen to be
$\mathcal{G}_{\mathcal{E}}\approx 0$ for $\Delta_{k}<\gamma\Delta^{*}$ where
$\gamma$ is chosen in such a way that $\gamma\Delta^{*}$ is of the order of
the lower-energy gaps for $h=h_{1}$ which allows for some modes to be bath-
engineered in the energizing bath stroke. This results in
$\rho_{k}^{B}=\rho_{k}^{A}$ for such modes.

The modes with $\Delta_{k}>\gamma\Delta^{*}$ are allowed to interact with the
energizing bath, leading the system to the steady state given by Equaion (25).

  * (d)

From B to C, the system is quenched and all modes evolve to reach $\rho^{C}$.

The total heat input and output of the system are calculated using

$\displaystyle\mathcal{Q}_{in}$ $\displaystyle=$
$\displaystyle\sum_{k}\mathcal{Q}_{in}^{k}$ (29)
$\displaystyle\mathcal{Q}_{out}$ $\displaystyle=$
$\displaystyle\sum_{k}\mathcal{Q}_{out}^{k}$ (30)

and the work output and efficiency of the engine are obtained using

$\displaystyle\mathcal{W}$ $\displaystyle=$
$\displaystyle\sum_{k}\mathcal{W}_{k}$ $\displaystyle\mathcal{W}_{k}$
$\displaystyle=$
$\displaystyle-\left(\mathcal{Q}_{in}^{k}+\mathcal{Q}_{out}^{k}\right)$
$\displaystyle\eta$ $\displaystyle=$
$\displaystyle-\frac{\mathcal{W}}{\sum_{k}\mathcal{Q}_{in}^{k}}.$ (31)

As discussed in Section II, a mode $k$ acts as an engine with non-zero work
output for $\mathcal{Q}_{in}^{k}>0$, $\mathcal{Q}_{out}^{k}<0$ and
$\mathcal{W}_{k}<0$.

We depict the variation of the output work and the efficiency of the engine
after implementing bath engineering in Figure 3. To obtain a complete picture,
we compare BEQE with finite-time engines without any control, finite-time
engines with the presence of shortcuts to adiabaticity in the unitary strokes,
and engines operating in the adiabatic limit, i.e.,
$\tau_{1}=\tau_{2}=\tau\rightarrow\infty$ (or infinite-time engines). As shown
in del Campo  _et al._ (2012), the STA Hamiltonian involves long-range
interactions. However, one can truncate the control Hamiltonian to $M$-body
terms to obtain a physically realizable approximate STA protocol. In Figure 3,
we present a comparison of the output work as a function of $\tau$
($=\tau_{1}=\tau_{2})$. As expected, engines using STA always perform better
than the finite-time engines without controls. However, interestingly, the
BEQE outperforms the engines using STA, as well as the perfectly adiabatic
engine, for a wide range of $\tau$ values, thus exhibiting the remarkable
benefits offered by the bath engineering technique.

Similarly, we also plot the efficiency $\eta$ as a function of $\tau$ (see
inset of Figure 3) and compare BEQE with engines operating with different
techniques. Here also we find that BEQE outperforms all other engines for the
same range of $\tau$ values as in the work output analysis. The expressions
for $\mathcal{|W|}_{adia}$ and $\eta_{adia}$ are given in Appendix A.

The fact that BEQE outperforms other engines can be explained using Figure 4,
where the $\mathcal{Q}_{in}$ and $\mathcal{Q}_{out}$ values for a perfect
adiabatic engine are plotted as a function of individual $k$ modes. In Figure
4, it can be seen that even when the engine works in the adiabatic limit,
there are some $k$ modes close to the critical mode which do not function as
an ’engine’ ($\mathcal{Q}_{in}<0$). To understand this better, let us consider
the adiabatic limit where $\rho_{k}^{C}=\rho_{k}^{B}$ in the eigenbasis. In
this limit, $\mathcal{Q}_{in}$ for each $k$ mode is given by (see Appendix A
for details)

$\displaystyle\mathcal{Q}_{in}^{k}$ $\displaystyle=$
$\displaystyle\frac{\Delta_{k}(h_{1})}{2}\Big{[}\frac{\left(e^{-\frac{\beta_{H}\Delta_{k}(h_{1})}{2}}-e^{\frac{\beta_{H}\Delta_{k}(h_{1})}{2}}\right)}{Z(h_{1})}$
(32) $\displaystyle-$
$\displaystyle\frac{\left(e^{-\frac{\beta_{C}\Delta_{k}(h_{2})}{2}}-e^{\frac{\beta_{C}\Delta_{k}(h_{2})}{2}}\right)}{Z(h_{2})}\Big{]}$
(33)

For $\mathcal{Q}_{in}^{k}$ to be positive,

$\displaystyle\frac{\sinh(\frac{\beta_{H}\Delta_{k}(h_{1})}{2})}{2+\cosh(\frac{\beta_{H}\Delta_{k}(h_{1})}{2})}<\frac{\sinh(\frac{\beta_{C}\Delta_{k}(h_{2})}{2})}{2+\cosh(\frac{\beta_{C}\Delta_{k}(h_{2})}{2})}.$
(34)

There can be modes for which this condition is not satisfied, resulting in
’non-engine’ modes in the adiabatic limit. BEQE helps to remove these non-
engine modes from participating in the non-unitary strokes, thereby boosting
the performance of the engine compared to the perfectly adiabatic engine. It
can be noted that the presence of non-engine modes is essential for the BEQE
to outperform the adiabatic engine. For instance, when
$T_{H}\rightarrow\infty$ ($\beta_{H}\rightarrow 0$), Equation 34 may be
satisfied for all modes so that the technique of BEQE will not provide better
results compared to the adiabatic engine. However, we emphasize that although
Equation (34) and the discussion above are specific to TIM WM, BEQE can be
expected to perform better than generic finite-time free fermionic quantum
critical engines, following the arguments presented in Section III.

Figure 3: $|\mathcal{W}|$ is plotted as a function of $\tau$ for the critical
engine using different techniques. Inset: $\eta$ is plotted as a function of
$\tau$. The parameters used are
$L=1000,h_{1}=10,h_{2}=1,T_{H}=20,T_{C}=1,\gamma=6.5,\tau_{1}=\tau_{2}=\tau$.
Figure 4: $\mathcal{Q}_{in}^{k}$ and $\mathcal{Q}_{out}^{k}$ are plotted as
functions of $k$ modes for the critical engine in the adiabatic limit. The
parameters are $L=1000,h_{1}=10,h_{2}=1,T_{H}=20,T_{C}=1$. The green dashed
line represents the zero of heat.

We note that BEQE depends on the appropriate choice of bath spectral function
(III), which again depends on $\tau$ through Equation (17). However, in
experimental setups, it might be difficult to change the bath-spectral
function for every change of $\tau$. Consequently, we examine the robustness
of the bath engineering protocol by plotting the work output and efficiency vs
$\tau$ for constant values of $\Delta^{*}$. In this case also, the results
show that the engine performance can be enhanced by choosing appropriate
constant values of $\Delta^{*}$ as shown in Figure 5, thus highlighting the
effectiveness of the proposed protocol in practical scenarios.

Figure 5: $\mathcal{|W|}$ is plotted as a function of $\tau$ using constant
value of $\Delta^{*}$ for all $\tau$. Inset: $\eta$ is plotted as a function
of $\tau$. Here, $\gamma(\Delta^{*}=2.1)=9$ and $\gamma(\Delta^{*}=0.3)=62$.
Other parameters are same as in Figure 3.

We point out that one may be able to further simplify the control protocol by
implementing bath engineering in only one of the non-unitary strokes (single-
stroke BEQE). In Figure 6 it can be seen that even a single-stroke BEQE
performed better than the finite- and infinite-time engines. Therefore, this
simplified protocol can be helpful as long as the overall work output can be
increased, which one can calculate following the mechanism discussed above,
even though there might be scenarios in which this simplified protocol may not
suffice.

Figure 6: $\mathcal{|W|}$ is plotted as a function of $\tau$ using a constant
value of $\Delta^{*}$ for all $\tau$ with bath engineering only in the
$\bf{C}\to\bf{D}$ non-unitary stroke (single-stroke BEQE) and with bath
engineering in both the non-unitary strokes (BEQE). Inset: $\eta$ is plotted
as a function of $\tau$. Here $\Delta^{*}=0.3,\gamma=62$. Other parameters are
same as in Figure 3.

Even though we have set $h_{2}=1$ for Figures 3, 5, and 6, the improvement
shown by BEQE persists when one crosses the quantum critical point during the
unitary strokes. This is shown in Figure 7 in which we have used (see Equation
(21))

$\displaystyle\Delta^{*}=C_{2}(1-h_{2})+C_{3}$ (35)

to improve the output work and efficiency for $h_{2}<1$.

We note that as we decrease $\tau$, more $k$ modes get excited and become
detrimental to the performance of the engine, thereby resulting in a
diminishing work output in all cases, in the small $\tau$ regime. This, in
turn, is also reflected in the decreasing power output for small values of
$\tau$ , as shown in Figures 8 and 9 in Appendix B.

Figure 7: $\mathcal{|W|}$ is plotted as a function of $\tau$ for (a)
$h_{2}=0.5,C_{3}=0.02$ and (b) $h_{2}=0.8,C_{3}=0.08$, with $\Delta^{*}$ given
by Equation (21). Here
$L=1000,h_{1}=10,T_{H}=20,T_{C}=1,\tau_{1}=\tau_{2}=\tau,C_{2}=2$ (see
Equation (27)).

##  V Conclusions

We have proposed the bath-engineered quantum engine, in which, through the
appropriate choice of bath spectral functions, one can dramatically boost the
performance of quantum critical engines based on free fermionic WMs. The
operation of the BEQE inherently depends on the Kibble–Zurek mechanism;
consequently, knowledge about the universality class of the WM and the strokes
of the Otto cycle suffice to implement this method. This is in stark contrast
to conventional methods such as shortcuts to adiabaticity, in which one may
need detailed knowledge about the eigenspectrum for their application. We
emphasize that the improvement in performance in the case of the BEQE is also
accompanied by the simplicity of formulating the control protocol, as compared
to more conventional techniques, such as STA, in which the unitary strokes may
involve non-trivial many-body terms in quantum systems driven through quantum
critical points del Campo  _et al._ (2012); Kolodrubetz  _et al._ (2017); Sels
and Polkovnikov (2017). We note that the costs of applying STA and the BEQE
also raise important questions. The cost of STA has been widely studied in the
literature, for example, in Campbell and Deffner (2017); Abah and Paternostro
(2019); Torrontegui  _et al._ (2017); Guéry-Odelin  _et al._ (2019). These
costs depend on the control Hamiltonian, which in the case of critical systems
may involve long-range interactions del Campo  _et al._ (2012) and can be
expected to be strongly dependent on the details of the setup in question. In
the case of the BEQE, the control approach involves introducing lower cut-offs
in the bath spectral functions (see Equations (III) and (21)). As discussed in
Naseem  _et al._ (2020); Gelbwaser-Klimovsky  _et al._ (2013), one can
engineer these cut-offs by introducing filters in the form of harmonic
oscillators with appropriate frequencies. As for STA, we expect the cost of
the BEQE to depend on the details of the implementation. However, a detailed
discussion regarding the comparison of costs of BEQE and STA goes beyond the
scope of the current paper. Furthermore, interestingly, in spite of the
simplicity of the proposed control protocol, our analysis with TIM WM shows
that the BEQE can outperform quantum engines assisted through STA, and even
infinite-time quantum engines, thus highlighting the significant benefits
offered by this control method. We also exhibit the robustness of the BEQE
control protocol by considering constant values of $\Delta^{*}$.

The control protocol proposed here can be expected to be most relevant in
quantum critical engines, owing to the high probability of excitations in
systems driven through quantum phase transitions, and the presence of many-
body interaction terms in the corresponding STA Hamiltonians. However, one can
expect the BEQE, which depends on introducing lower cut-offs into bath
spectral functions, to be applicable even away from criticality, for multi-
level WMs involving non-equispaced energy levels. On the other hand, as one
can infer from Equations (III) and (21), this protocol becomes invalid in the
presence of equispaced energy levels.

Several existing setups can be suitable for the experimental realization of
the BEQE, such as trapped ions Roßnagel  _et al._ (2016); Ulm  _et al._
(2013); Pyka  _et al._ (2013); Maslennikov  _et al._ (2019); von Lindenfels
_et al._ (2019), optical lattices Schreiber  _et al._ (2015), superconducting
qubits, nitrogen vacancy centers in diamond Klatzow  _et al._ (2019), NMR
qubit systems Peterson  _et al._ (2019), etc. For example, quantum simulators
based on trapped ions have already been used to study the Kibble–Zurek
mechanism in momentum space Cui  _et al._ (2020).

Finally, we note that although this technique appears to be highly successful
in the case of free fermionic WMs, open questions remain in the case of its
application with non-integrable WMs, where such non-interacting $k$ modes may
not exist. For example, one can choose the WM to be the antiferromagnetic
transverse Ising model with a longitudinal field (LTIM), described by the
Hamiltonian

$H=J\sum_{i}\sigma_{i}^{z}\sigma_{i+1}^{z}-B_{x}(t)\sum_{i}\sigma_{i}^{x}-B_{z}\sum_{i}\sigma_{i}^{z}.$
(36)

Here $J$ is the strength of antiferromagnetic interaction, $B_{z}$ is a
longitudinal field, and $B_{x}$ denotes a time-dependent transverse field. The
competition between $J$ and $B_{z}$ leads to a quantum phase transition from
the antiferromagnetic state to the paramagnetic state at a critical value of
$B_{x}^{c}$ for a fixed value of $B_{z}$ Sharma  _et al._ (2015); de Alcantara
Bonfim  _et al._ (2019). One can model an Otto cycle using LTIM WM and
implement the BEQE as described above (see Appendix C). However, preliminary
studies suggest that, unlike the case of the integrable model, there is no
improvement in the output of the BEQE in this case (see Figure 10). This can
be attributed to the absence of non-interacting momentum modes, as obtained
for free fermionic systems. However, additional rigorous studies are needed to
acquire a deeper understanding of the possibility of the application of BEQE
to quantum engines based on more generic non-integrable WMs.

###### Acknowledgements.

R.B.S. and U.D. acknowledge CHANDRA super cluster of IIT Palakkad, in which
all numerical simulations were performed. V.M. acknowledges support from
Science and Engineering Research Board (SERB) through MATRICS (Project No.
MTR/2021/000055) and Seed Grant from IISER Berhampur.

##  Appendix A Adiabatic Evolution of TIM

The energies at the end of stroke $i$ are calculated using the expresssion

${E}_{i}=\sum_{k}\text{Tr}[H_{k}^{i}\rho_{k}^{i}].$ (37)

where $i={A},{B},{C},{D}$.

  * (i)

At B: The density matrix is given by Equation (25) and the Hamiltonian in the
diagonal basis takes the form

$H_{k}(h_{1})=\begin{bmatrix}-\epsilon_{k}(h_{1})&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&\epsilon_{k}(h_{1})\end{bmatrix}.$ (38)

The energy ${E}_{B}$ can be calculated as follows

$\text{Tr}[H_{k}(h_{1})\rho_{k}^{B}]=\frac{\epsilon_{k}(h_{1})}{Z(h_{1})}\left(e^{-\beta_{H}\epsilon_{k}(h_{1})}-e^{\beta_{H}\epsilon_{k}(h_{1})}\right)$
(39)

or

$\displaystyle{E}_{B}$ $\displaystyle=$
$\displaystyle\sum_{k}\text{Tr}[H_{k}(h_{1})\rho_{k}^{B}]$ (40)
$\displaystyle=$
$\displaystyle\sum_{k}\frac{\epsilon_{k}(h_{1})}{Z(h_{1})}\left(e^{-\beta_{H}\epsilon_{k}(h_{1})}-e^{\beta_{H}\epsilon_{k}(h_{1})}\right)$
(41)

  * (ii)

At C: If the evolution is purely adiabatic, the populations in the eigenenergy
levels do not change, resulting in

$\rho_{k}^{C,adia}=\rho_{k}^{B}=\begin{bmatrix}\frac{e^{\beta_{H}\epsilon_{k}(h_{1})}}{Z(h_{1})}&0&0&0\\\
0&\frac{1}{Z(h_{1})}&0&0\\\ 0&0&\frac{1}{Z(h_{1})}&0\\\
0&0&0&\frac{e^{-\beta_{H}\epsilon_{k}(h_{1})}}{Z(h_{2})}\end{bmatrix}$ (42)

and

$H_{k}(h_{2})=\begin{bmatrix}-\epsilon_{k}(h_{2})&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&\epsilon_{k}(h_{2})\end{bmatrix}$ (43)

with $\epsilon_{k}(h_{2})=2\sqrt{(h_{2}-\cos k)^{2}+{\sin k}^{2}}$.
Therefore,

${E}_{C}^{adia}=\sum_{k}\frac{\epsilon_{k}(h_{2})}{Z(h_{1})}\left(e^{-\beta_{H}\epsilon_{k}(h_{1})}-e^{\beta_{H}\epsilon_{k}(h_{1})}\right)$
(44)

  * (iii)

At D: This energy can be calculated similarly to that at ${B}$ so that

${E}_{D}=\sum_{k}\frac{\epsilon_{k}(h_{2})}{Z(h_{2})}\left(e^{-\beta_{C}\epsilon_{k}(h_{2})}-e^{\beta_{C}\epsilon_{k}(h_{2})}\right)$
(45)

  * (iv)

At A: Following the same procedure used to calculate the energy
${E}_{C}^{adia}$ in order to find the energy ${E}_{A}^{adia}$, we obtain

${E}_{A}^{adia}=\sum_{k}\frac{\epsilon_{k}(h_{1})}{Z(h_{2})}\left(e^{-\beta_{C}\epsilon_{k}(h_{2})}-e^{\beta_{C}\epsilon_{k}(h_{2})}\right)$
(46)

The input heat energy absorbed by the WM in the non-unitary stroke ${A}\to{B}$
can be easily calculated, and is given by

$\displaystyle\mathcal{Q}_{in}^{adia}$ $\displaystyle=$
$\displaystyle{E}_{B}-{E}_{A}^{adia}$ (47) $\displaystyle=$
$\displaystyle\sum_{k}\epsilon_{k}(h_{1})\\{\frac{\left(e^{-\beta_{H}\epsilon_{k}(h_{1})}-e^{\beta_{H}\epsilon_{k}(h_{1})}\right)}{Z(h_{1})}$
$\displaystyle-$
$\displaystyle\frac{\left(e^{-\beta_{C}\epsilon_{k}(h_{2})}-e^{\beta_{C}\epsilon_{k}(h_{2})}\right)}{Z(h_{2})}\\}$

Similarly,

$\displaystyle\mathcal{Q}_{out}^{adia}$ $\displaystyle=$
$\displaystyle{E}_{D}-{E}_{C}^{adia}$ (48) $\displaystyle=$
$\displaystyle\sum_{k}\epsilon_{k}(h_{2})\\{\frac{\left(e^{-\beta_{C}\epsilon_{k}(h_{2})}-e^{\beta_{C}\epsilon_{k}(h_{2})}\right)}{Z(h_{2})}$
$\displaystyle-$
$\displaystyle\frac{\left(e^{-\beta_{H}\epsilon_{k}(h_{1})}-e^{\beta_{H}\epsilon_{k}(h_{1})}\right)}{Z(h_{1})}\\}.$

We can now calculate the output work of the engine in the adiabatic limit,
which is given by

$\displaystyle\mathcal{|W|}_{adia}$
$\displaystyle=\sum_{k}\left(\epsilon_{k}(h_{1})-\epsilon_{k}(h_{2})\right)$
(49) $\displaystyle\\{$
$\displaystyle\frac{\left(e^{-\beta_{H}\epsilon_{k}(h_{1})}-e^{\beta_{H}\epsilon_{k}(h_{1})}\right)}{Z(h_{1})}$
$\displaystyle-$
$\displaystyle\frac{\left(e^{-\beta_{C}\epsilon_{k}(h_{2})}-e^{\beta_{C}\epsilon_{k}(h_{2})}\right)}{Z(h_{2})}\\}$

resulting in

$\displaystyle\eta_{adia}=1-\frac{\sum_{k}\epsilon_{k}(h_{2})\\{\frac{\left(e^{-\beta_{H}\epsilon_{k}(h_{1})}-e^{\beta_{H}\epsilon_{k}(h_{1})}\right)}{Z(h_{1})}-\frac{\left(e^{-\beta_{C}\epsilon_{k}(h_{2})}-e^{\beta_{C}\epsilon_{k}(h_{2})}\right)}{Z(h_{2})}\\}}{\sum_{k}\epsilon_{k}(h_{1})\\{\frac{\left(e^{-\beta_{H}\epsilon_{k}(h_{1})}-e^{\beta_{H}\epsilon_{k}(h_{1})}\right)}{Z(h_{1})}-\frac{\left(e^{-\beta_{C}\epsilon_{k}(h_{2})}-e^{\beta_{C}\epsilon_{k}(h_{2})}\right)}{Z(h_{2})}\\}}.$
(50)

##  Appendix B Power Output for BEQE

The power output for an engine is defined as

$\mathcal{P}=\frac{-\mathcal{W}}{\tau_{\text{total}}}$ (51)

where $\tau_{\text{total}}=\tau_{1}+\tau_{2}+\tau_{H}+\tau_{C}$, with
$\tau_{H}(\tau_{C})$ being the time duration of the A $\rightarrow$ B (C
$\rightarrow$ D) non-unitary stroke. As shown in Figures 8 and 9, the BEQE
outperforms finite-time engines without controls and with STA for a wide range
of $\tau$ values.

Figure 8: $\mathcal{|P|}$ is plotted as a function of $\tau$ using different
values of $\Delta^{*}$. All parameters are same as in Figure 5. We have set
$\tau_{H}+\tau_{C}=4$.  Figure 9: $\mathcal{|P|}$ is plotted as a function of
$\tau$ for $h_{2}<1$. Here,
$L=1000,h_{1}=10,T_{H}=20,T_{C}=1,\tau_{1}=\tau_{2},\Delta^{*}(h_{2}=0.5)=1.02,\Delta^{*}(h_{2}=0.8)=0.48$.
We have set $\tau_{H}+\tau_{C}=4$.

##  Appendix C BEQE Using LTIM WM

With LTIM as the WM of the quantum Otto cycle, the transverse field $B_{x}$ is
changed from an $h_{1}$ to an $h_{2}$ value during the unitary strokes.
Following the technique presented here, bath engineering can be implemented by
evaluating the corresponding $\Delta^{*}$ using the Kibble–Zurek mechanism and
then choosing an appropriate cut-off for bath spectral functions such that
energy gaps which are less than $\Delta^{*}$ do not participate in the
dynamics. LTIM falls under the same universality class as that of TIM.

Let us discuss the implementation of single-stroke bath engineering in LTIM.
When performing bath engineering, some levels will not be allowed to
thermalize, depending upon the energy gap. There are $2^{L}$ energy levels for
a system size $L$ and thus $2^{L}-1$ energy gaps. Those energy levels having
gaps less than the threshold value of $\Delta^{*}=\left(\frac{\nu
z(h_{1}-h_{2})}{\tau}\right)^{\frac{\nu z}{1+\nu z}}$ will not thermalize. On
the other hand, those with energy gaps greater than $\left(\frac{\nu
z(h_{1}-h_{2})}{\tau}\right)^{\frac{\nu z}{1+\nu z}}$ will thermalize
according to the equation

$p_{i}=p_{i-1}e^{-(E_{i}-E_{i-1})/T}$ (52)

where $p_{i}$ and $p_{i-1}$ are the populations in the $i^{th}$ and
$(i-1)^{th}$ energy levels.

  * •

In the C $\rightarrow$ D non-unitary stroke, the Hamiltonian is with a
transverse field $h_{2}$. The energy gaps are compared with the $\Delta^{*}$
value. Those energy levels having gaps greater than $\Delta^{*}$ are allowed
to interact with the decaying bath in the C $\rightarrow$ D stroke and thus
thermalize according to

$\frac{p_{i+1}^{D}}{p_{i}^{D}}=e^{-(E_{i+1}-E_{i})/T_{c}}$ (53)

where $E_{i+1}-E_{i}>\Delta^{*}$.

In order to apply bath engineering in the case of gaps that are less than
$\Delta^{*}$, i.e, when $E_{i+1}-E_{i}<\Delta^{*}$, we have two possibilities.

    * (i)

If $E_{i}-E_{i-1}>\Delta^{*}$, the populations are determined by the condition

$p_{i}^{D}+p_{i-1}^{D}=p_{i}^{C}+p_{i-1}^{C}.$ (54)

    * (ii)

If $E_{i}-E_{i-1}<\Delta^{*}$, the $i^{th}$ level does not interact with any
other level, leading to $p_{i}^{D}=p_{i}^{C}$.

Solving these system of equations, along with the condition that
$\sum_{i}p_{i}^{D}=1$, we obtain the populations of all the other energy
levels at D. Thus, we have $\rho^{\prime D}$, which is the state reached after
carrying out bath engineering.

  * •

From D to A, $h_{2}$ is changed back to $h_{1}$ from $\rho^{\prime D}$ using
the evolution equation

$\frac{d\rho}{dt}=-i[H,\rho]$ (55)

which gives the new density matrix at A, $\rho^{\prime A}$.

  * •

At A, the system with the Hamiltonian $H(h_{1})$ is connected to the
energizing bath. All energy levels interact with each other, resulting in the
steady state at B.

  * •

B $\rightarrow$ C stroke, $h_{1}$ changed to $h_{2}$ by evolving the system
from $\rho^{B}$ to obtain $\rho^{C}$.

Now, for the bath-engineered engine,

$\displaystyle\mathcal{Q}_{in}^{\prime}$ $\displaystyle=$
$\displaystyle{E}_{B}-{E}_{A}^{\prime}$ (56)
$\displaystyle\mathcal{Q}_{out}^{\prime}$ $\displaystyle=$
$\displaystyle{E}_{D}^{\prime}-{E}_{C}$ (57)
$\displaystyle\mathcal{W}^{\prime}$ $\displaystyle=$
$\displaystyle-(\mathcal{Q}_{in}^{\prime}+\mathcal{Q}_{out}^{\prime})$ (58)
$\displaystyle\eta^{\prime}$ $\displaystyle=$
$\displaystyle-\mathcal{W}^{\prime}/\mathcal{Q}_{in}^{\prime}$ (59)

Numerical analysis suggests that bath engineering failed to improve the
performance of the engine in this case (see Figure 10). We also investigated
the implementation of bath engineering in both energizing and decaying
strokes, and again failed to improve the performance of the engine.

Figure 10: $\eta$ is plotted as a function of $\tau$ for the BEQE using LTIM
as the WM compared with a finite-time engine. Inset: $\mathcal{|W|}$ as a
function of $\tau$. The parameters used are
$L=6,h_{1}=10,h_{2}=0.75,T_{H}=500,T_{C}=0.1,\tau_{1}=\tau_{2}=\tau$.

## References

  * Alicki and Kosloff (2018) Robert Alicki and Ronnie Kosloff, “Introduction to quantum thermodynamics: History and prospects”, in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by Felix Binder, Luis A. Correa, Christian Gogolin, Janet Anders, and Gerardo Adesso (Springer International Publishing, Cham, 2018) pp. 1–33.
  * Vinjanampathy and Anders (2016) Sai Vinjanampathy and Janet Anders, “Quantum thermodynamics”, Contemporary Physics 57, 545–579 (2016).
  * Gemmer  _et al._ (2009) Jochen Gemmer, Mathias Michel, and Günter Mahler,  _Quantum thermodynamics: Emergence of thermodynamic behavior within composite quantum systems_, Vol. 784 (Springer: Berlin/Heidelberg, Germany, 2009).
  * Kosloff (2013) Ronnie Kosloff, “Quantum thermodynamics: A dynamical viewpoint”, Entropy 15, 2100–2128 (2013).
  * Binder  _et al._ (2018) Felix Binder, Luis A Correa, Christian Gogolin, Janet Anders, and Gerardo Adesso, “Thermodynamics in the quantum regime”, Fundamental Theories of Physics 195, 1–2 (2018).
  * Bhattacharjee and Dutta (2021) Sourav Bhattacharjee and Amit Dutta, “Quantum thermal machines and batteries”, The European Physical Journal B 94, 239 (2021).
  * Myers  _et al._ (2022) Nathan M. Myers, Obinna Abah, and Sebastian Deffner, “Quantum thermodynamic devices: From theoretical proposals to experimental reality”, AVS Quantum Science 4, 027101 (2022).
  * Erdman and Noé (2022) Paolo A. Erdman and Frank Noé, “Identifying optimal cycles in quantum thermal machines with reinforcement-learning”, npj Quantum Information 8, 1 (2022).
  * Khait  _et al._ (2022) Ilia Khait, Juan Carrasquilla, and Dvira Segal, “Optimal control of quantum thermal machines using machine learning”, Phys. Rev. Research 4, L012029 (2022).
  * Quan  _et al._ (2007) H. T. Quan, Yu-xi Liu, C. P. Sun, and Franco Nori, “Quantum thermodynamic cycles and quantum heat engines”, Phys. Rev. E 76, 031105 (2007).
  * Mukherjee and Divakaran (2021) Victor Mukherjee and Uma Divakaran, “Many-body quantum thermal machines”, Journal of Physics: Condensed Matter 33, 454001 (2021).
  * Gelbwaser-Klimovsky  _et al._ (2015) David Gelbwaser-Klimovsky, Wolfgang Niedenzu, and Gershon Kurizki, “Chapter twelve - thermodynamics of quantum systems under dynamical control”, (Academic Press: Cambridge, MA, USA, 2015) pp. 329–407.
  * Mukherjee  _et al._ (2019) V. Mukherjee, A. Zwick, A. Ghosh, X. Chen, and G. Kurizki, “Enhanced precision of low-temperature quantum thermometry via dynamical control”, Commun Phys 2, 162 (2019).
  * del Campo  _et al._ (2014) A. del Campo, J. Goold, and M. Paternostro, “More bang for your buck: Super-adiabatic quantum engines”, Sci. Rep. 4, 6208 (2014).
  * Hartmann  _et al._ (2020) Andreas Hartmann, Victor Mukherjee, Wolfgang Niedenzu, and Wolfgang Lechner, “Many-body quantum heat engines with shortcuts to adiabaticity”, Phys. Rev. Research 2, 023145 (2020).
  * Deng  _et al._ (2013) Jiawen Deng, Qing-hai Wang, Zhihao Liu, Peter Hänggi, and Jiangbin Gong, “Boosting work characteristics and overall heat-engine performance via shortcuts to adiabaticity: Quantum and classical systems”, Phys. Rev. E 88, 062122 (2013).
  * Beau  _et al._ (2016) Mathieu Beau, Juan Jaramillo, and Adolfo del Campo, “Scaling-up quantum heat engines efficiently via shortcuts to adiabaticity”, Entropy 18, 168 (2016).
  * del Campo  _et al._ (2018) Adolfo del Campo, Aurélia Chenu, Shujin Deng, and Haibin Wu, “Friction-free quantum machines”, in  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, edited by Felix Binder, Luis A. Correa, Christian Gogolin, Janet Anders, and Gerardo Adesso (Springer International Publishing, Cham, 2018) pp. 127–148.
  * Sels and Polkovnikov (2017) Dries Sels and Anatoli Polkovnikov, “Minimizing irreversible losses in quantum systems by local counterdiabatic driving”, Proceedings of the National Academy of Sciences 114, E3909–E3916 (2017).
  * del Campo  _et al._ (2012) Adolfo del Campo, Marek M. Rams, and Wojciech H. Zurek, “Assisted finite-rate adiabatic passage across a quantum critical point: Exact solution for the quantum ising model”, Phys. Rev. Lett. 109, 115703 (2012).
  * Kolodrubetz  _et al._ (2017) Michael Kolodrubetz, Dries Sels, Pankaj Mehta, and Anatoli Polkovnikov, “Geometry and non-adiabatic response in quantum and classical systems”, Physics Reports 697, 1–87 (2017), geometry and non-adiabatic response in quantum and classical systems.
  * Sachdev (2011) Subir Sachdev,  _Quantum Phase Transitions_, 2nd ed. (Cambridge University Press, 2011).
  * Campisi and Fazio (2016) Michele Campisi and Rosario Fazio, “The power of a critical heat engine”, Nature communications 7, 1–5 (2016).
  * Chen  _et al._ (2019) Yang-Yang Chen, Gentaro Watanabe, Yi-Cong Yu, Xi-Wen Guan, and Adolfo del Campo, “An interaction-driven many-particle quantum heat engine and its universal behavior”, npj Quantum Information 5, 88 (2019).
  * Ma  _et al._ (2017) Yu-Han Ma, Shan-He Su, and Chang-Pu Sun, “Quantum thermodynamic cycle with quantum phase transition”, Phys. Rev. E 96, 022143 (2017).
  * Piccitto  _et al._ (2022) Giulia Piccitto, Michele Campisi, and Davide Rossini, “The ising critical quantum otto engine”, arXiv preprint arXiv:2205.09528 (2022).
  * Fogarty and Busch (2020) Thomás Fogarty and Thomas Busch, “A many-body heat engine at criticality”, Quantum Science and Technology 6, 015003 (2020).
  * Revathy  _et al._ (2020) B. S Revathy, Victor Mukherjee, Uma Divakaran, and Adolfo del Campo, “Universal finite-time thermodynamics of many-body quantum machines from kibble-zurek scaling”, Phys. Rev. Research 2, 043247 (2020).
  * Dutta  _et al._ (2015) A. Dutta, G. Aeppli, B. K. Chakrabarti, U. Divakaran, T. F. Rosenbaum, and D. Sen,  _Quantum phase transitions in transverse field spin models: from statistical physics to quantum information_ (Cambridge University Press, Cambridge, 2015).
  * Polkovnikov  _et al._ (2011) Anatoli Polkovnikov, Krishnendu Sengupta, Alessandro Silva, and Mukund Vengalattore, “Colloquium: Nonequilibrium dynamics of closed interacting quantum systems”, Rev. Mod. Phys. 83, 863–883 (2011).
  * Dziarmaga (2010) Jacek Dziarmaga, “Dynamics of a quantum phase transition and relaxation to a steady state”, Advances in Physics 59, 1063–1189 (2010).
  * Mukherjee  _et al._ (2016) V. Mukherjee, W. Niedenzu, A. G. Kofman, and G. Kurizki, “Speed and efficiency limits of multilevel incoherent heat engines”, Phys. Rev. E 94, 062109 (2016).
  * Keck  _et al._ (2017) M. Keck, S. Montangero, G. E. Santoro, R. Fazio, and D. Rossini, “Dissipation in adiabatic quantum computers: lessons from an exactly solvable model”, New Journal of Physics 19, 113029 (2017).
  * Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Deng  _et al._ (2011) Shusa Deng, Gerardo Ortiz, and Lorenza Viola, “Dynamical critical scaling and effective thermalization in quantum quenches: Role of the initial state”, Phys. Rev. B 83, 094304 (2011).
  * Roßnagel  _et al._ (2014) J. Roßnagel, O. Abah, F. Schmidt-Kaler, K. Singer, and E. Lutz, “Nanoscale heat engine beyond the carnot limit”, Phys. Rev. Lett. 112, 030602 (2014).
  * Niedenzu  _et al._ (2018) Wolfgang Niedenzu, Victor Mukherjee, Arnab Ghosh, Abraham G Kofman, and Gershon Kurizki, “Quantum engine efficiency bound beyond the second law of thermodynamics”, Nature communications 9, 1–13 (2018).
  * Zurek  _et al._ (2005) Wojciech H. Zurek, Uwe Dorner, and Peter Zoller, “Dynamics of a quantum phase transition”, Phys. Rev. Lett. 95, 105701 (2005).
  * Polkovnikov (2005) Anatoli Polkovnikov, “Universal adiabatic dynamics in the vicinity of a quantum critical point”, Phys. Rev. B 72, 161201 (2005).
  * Cui  _et al._ (2020) Jin-Ming Cui, Fernando Javier Gómez-Ruiz, Yun-Feng Huang, Chuan-Feng Li, Guang-Can Guo, and Adolfo del Campo, “Experimentally testing quantum critical dynamics beyond the kibble–zurek mechanism”, Communications Physics 3, 44 (2020).
  * Bando  _et al._ (2020) Yuki Bando, Yuki Susa, Hiroki Oshiyama, Naokazu Shibata, Masayuki Ohzeki, Fernando Javier Gómez-Ruiz, Daniel A. Lidar, Sei Suzuki, Adolfo del Campo, and Hidetoshi Nishimori, “Probing the universality of topological defect formation in a quantum annealer: Kibble-zurek mechanism and beyond”, Phys. Rev. Research 2, 033369 (2020).
  * Damski and Zurek (2006) Bogdan Damski and Wojciech H. Zurek, “Adiabatic-impulse approximation for avoided level crossings: From phase-transition dynamics to landau-zener evolutions and back again”, Phys. Rev. A 73, 063405 (2006).
  * Lieb  _et al._ (1961) Elliott Lieb, Theodore Schultz, and Daniel Mattis, “Two soluble models of an antiferromagnetic chain”, Annals of Physics 16, 407 – 466 (1961).
  * Pfeuty (1970) Pierre Pfeuty, “The one-dimensional ising model with a transverse field”, Annals of Physics 57, 79 – 90 (1970).
  * Bunder and McKenzie (1999) J. E. Bunder and Ross H. McKenzie, “Effect of disorder on quantum phase transitions in anisotropic xy spin chains in a transverse field”, Phys. Rev. B 60, 344–358 (1999).
  * Bandyopadhyay  _et al._ (2018) Souvik Bandyopadhyay, Sudarshana Laha, Utso Bhattacharya, and Amit Dutta, “Exploring the possibilities of dynamical quantum phase transitions in the presence of a markovian bath”, Scientific Reports 8, 11921 (2018).
  * Campbell and Deffner (2017) Steve Campbell and Sebastian Deffner, “Trade-off between speed and cost in shortcuts to adiabaticity”, Phys. Rev. Lett. 118, 100601 (2017).
  * Abah and Paternostro (2019) Obinna Abah and Mauro Paternostro, “Shortcut-to-adiabaticity otto engine: A twist to finite-time thermodynamics”, Phys. Rev. E 99, 022110 (2019).
  * Torrontegui  _et al._ (2017) E. Torrontegui, I. Lizuain, S. González-Resines, A. Tobalina, A. Ruschhaupt, R. Kosloff, and J. G. Muga, “Energy consumption for shortcuts to adiabaticity”, Phys. Rev. A 96, 022133 (2017).
  * Guéry-Odelin  _et al._ (2019) D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, “Shortcuts to adiabaticity: Concepts, methods, and applications”, Rev. Mod. Phys. 91, 045001 (2019).
  * Naseem  _et al._ (2020) M. Tahir Naseem, Avijit Misra, Özgür E. Müstecaplioğlu, and Gershon Kurizki, “Minimal quantum heat manager boosted by bath spectral filtering”, Phys. Rev. Research 2, 033285 (2020).
  * Gelbwaser-Klimovsky  _et al._ (2013) D. Gelbwaser-Klimovsky, R. Alicki, and G. Kurizki, “Minimal universal quantum heat machine”, Phys. Rev. E 87, 012140 (2013).
  * Roßnagel  _et al._ (2016) Johannes Roßnagel, Samuel T. Dawkins, Karl N. Tolazzi, Obinna Abah, Eric Lutz, Ferdinand Schmidt-Kaler, and Kilian Singer, “A single-atom heat engine”, Science 352, 325–329 (2016).
  * Ulm  _et al._ (2013) S Ulm, J Roßnagel, G Jacob, C Degünther, ST Dawkins, UG Poschinger, R Nigmatullin, A Retzker, MB Plenio, F Schmidt-Kaler, _et al._ , “Observation of the kibble–zurek scaling law for defect formation in ion crystals”, Nature communications 4, 1–7 (2013).
  * Pyka  _et al._ (2013) K. Pyka, J. Keller, H. L. Partner, R. Nigmatullin, T. Burgermeister, D. M. Meier, K. Kuhlmann, A. Retzker, M. B. Plenio, W. H. Zurek, A. del Campo, and T. E. Mehlstäubler, “Topological defect formation and spontaneous symmetry breaking in ion coulomb crystals”, Nat. Comm. 4, 2291 (2013).
  * Maslennikov  _et al._ (2019) Gleb Maslennikov, Shiqian Ding, Roland Hablützel, Jaren Gan, Alexandre Roulet, Stefan Nimmrichter, Jibo Dai, Valerio Scarani, and Dzmitry Matsukevich, “Quantum absorption refrigerator with trapped ions”, Nature Communications 10, 202 (2019).
  * von Lindenfels  _et al._ (2019) D. von Lindenfels, O. Gräb, C. T. Schmiegelow, V. Kaushal, J. Schulz, Mark T. Mitchison, John Goold, F. Schmidt-Kaler, and U. G. Poschinger, “Spin heat engine coupled to a harmonic-oscillator flywheel”, Phys. Rev. Lett. 123, 080602 (2019).
  * Schreiber  _et al._ (2015) Michael Schreiber, Sean S. Hodgman, Pranjal Bordia, Henrik P. Lüschen, Mark H. Fischer, Ronen Vosk, Ehud Altman, Ulrich Schneider, and Immanuel Bloch, “Observation of many-body localization of interacting fermions in a quasirandom optical lattice”, Science 349, 842–845 (2015).
  * Klatzow  _et al._ (2019) James Klatzow, Jonas N. Becker, Patrick M. Ledingham, Christian Weinzetl, Krzysztof T. Kaczmarek, Dylan J. Saunders, Joshua Nunn, Ian A. Walmsley, Raam Uzdin, and Eilon Poem, “Experimental demonstration of quantum effects in the operation of microscopic heat engines”, Phys. Rev. Lett. 122, 110601 (2019).
  * Peterson  _et al._ (2019) John P. S. Peterson, Tiago B. Batalhão, Marcela Herrera, Alexandre M. Souza, Roberto S. Sarthour, Ivan S. Oliveira, and Roberto M. Serra, “Experimental characterization of a spin quantum heat engine”, Phys. Rev. Lett. 123, 240601 (2019).
  * Sharma  _et al._ (2015) Shraddha Sharma, Sei Suzuki, and Amit Dutta, “Quenches and dynamical phase transitions in a nonintegrable quantum ising model”, Phys. Rev. B 92, 104306 (2015).
  * de Alcantara Bonfim  _et al._ (2019) O. F. de Alcantara Bonfim, B Boechat, and J Florencio, “Ground-state properties of the one-dimensional transverse ising model in a longitudinal magnetic field”, Physical Review E 99, 012122 (2019).


# Finite element discretization of a biological network formation system: a
preliminary study

Clarissa Astuto  King Abdullah University of Science and Technology (KAUST),
4700, Thuwal, Saudi Arabia  Daniele Boffi  King Abdullah University of Science
and Technology (KAUST), 4700, Thuwal, Saudi Arabia  Department of Mathematics
”F. Casorati”, University of Pavia, Pavia, Italy  Fabio Credali  King Abdullah
University of Science and Technology (KAUST), 4700, Thuwal, Saudi Arabia
Department of Mathematics ”F. Casorati”, University of Pavia, Pavia, Italy

###### Abstract

A finite element discretization is developed for the Cai-Hu model, describing
the formation of biological networks. The model consists of a non linear
elliptic equation for the pressure $p$ and a non linear reaction-diffusion
equation for the conductivity tensor $\mathbb{C}$. The problem requires high
resolution due to the presence of multiple scales, the stiffness in all its
components and the non linearities. We propose a low order finite element
discretization in space coupled with a semi-implicit time advancing scheme.
The code is verified with several numerical tests performed with various
choices for the parameters involved in the system. In absence of the exact
solution, we apply Richardson extrapolation technique to estimate the order of
the method.

##  1 Introduction

The branching structure of biological transport networks, such as leaf
venations and blood flow, has been the object of extensive research due to
their widespread existence in living organisms and natural phenomena, such as,
leaf venation in plants [1], vascular pattern formation [2], mammalian
circulatory systems and neural networks transporting electric charge [3, 4].

It has been suggested that the optimization of energy consumption, as a result
of natural selection, governs the branching structure. In 2013, Hu and Cai
[18] proposed a global energy functional approach that takes into
consideration a metabolic cost to sustain the network. This consumption is
also influenced by the the local gradient of pressure that can modify the
diameter of the vessels as an adaptive response to the stress [5].

The considered equations have been recently analyzed, from the numerical point
of view, making use of finite difference schemes [6, 7]. Finite element
methods (FEM) (see, e.g., [10]) are appealing for their flexibility and are
widely used by the mathematical community. In this work we present a low order
finite element method for solving the biological network PDE system in tensor
form and capture the steady state solution. The approximation of biological
network formation has been already studied by means of mixed finite elements
[8, 15] in the case of the corresponding vector model. A comparison between
the two different versions, tensor and vector models, has been showed in [7].

Our code is based on the finite elements library FreeFEM [12, 16], which
proposes a large variety of triangular finite elements to solve partial
differential equations. FreeFEM has been already used for numerical
computations of stationary states of fast-rotating Bose–Einstein condensates
[11], and for the minimization of energy functionals associated with Ginzburg-
Landau and Schrödinger models [21, 22].

The paper is organized as follows: in Section 2, we introduce the tensor PDE
model and we show the derivation of the weak formulation. In Section 3, we
describe the discrete version of the problem by coupling finite elements in
space with a semi-implicit time advancing scheme based on Backward Euler. In
Section 4, we report and discuss accuracy tests and results of the numerical
simulations. Finally, in Section 5, we summarize the results and draw some
conclusions.

##  2 Mathematical model

We consider an elliptic-parabolic system of partial differential equations for
the conductivity tensor, $\mathbb{C}\in\mathbb{R}^{2\times 2}$, and a pressure
scalar force, $p\in\mathbb{R}$. The model, introduced in [18], reads

$\displaystyle-\nabla\cdot\left((r\mathbb{I}+\mathbb{C})\nabla p\right)$
$\displaystyle=S\qquad{\rm in}\,\Omega$ (1)
$\displaystyle\frac{\partial\mathbb{C}}{\partial
t}-D^{2}\Delta\mathbb{C}-c^{2}\nabla p\otimes\nabla
p+\alpha\|\mathbb{C}\|^{\gamma-2}\mathbb{C}$ $\displaystyle=0\qquad{\rm
in}\,\Omega.$ (2)

The term $S=S(\vec{x})$ denotes the distribution of sources and sinks, which
is a known datum. The function $r:\Omega\to\mathbb{R}^{+}$, with $r(x)\geq
r_{0}>0$, describes the isotropic background permeability of the medium. In
Eq. (2) the diffusion coefficient ${D>0}$ controls the random effects in the
transportation medium, while the activation parameter $c>0$ describes the
tendency of the network to align with the pressure gradient. The reaction term
$\alpha\|\mathbb{C}\|^{\gamma-2}\mathbb{C}$, where $\|\cdot\|$ denotes the
Frobenius norm, models the metabolic cost of maintaining the network
structure, with metabolic coefficient $\alpha>0$ and metabolic exponent
$\gamma>0$. For instance, for leaf venation in plants, $1/2\leq\gamma\leq 1$
[18, 17]. For this range of values, there is no proof of convergence to a
unique steady state, while it has been proved for $\gamma>1$ in [13]. Both
cases have been numerically investigated in [7].

We consider Eqs.(1–2) in a bounded domain $\Omega\subset\mathbb{R}^{2}$, with
smooth boundary $\partial\Omega$. We choose homogeneous Neumann boundary
conditions for $p$ and $\mathbb{C}$

$\nabla
p(t,\vec{x})\cdot\nu=0,\quad\nabla\mathbb{C}(t,\vec{x})\cdot\nu=0,\quad\text{for
}\vec{x}\in\partial\Omega,\,t\geq 0,$ (3)

where $\nu$ denotes the outer normal vector to $\partial\Omega$. With the
homogeneous Neumann boundary condition for the pressure $p$, we impose the
global mass balance

$\int_{\Omega}S(\vec{x})\,d\vec{x}=0$ (4)

to ensure solvability of Eq. (1). We choose Neumann boundary conditions for
$\mathbb{C}$ because in [7] we observed that homogeneous Dirichlet boundary
conditions lead to the formation of boundary layers in the solutions. On the
contrary, homogeneous Neumann boundary conditions seem able to bypass this
problem. Let us also notice that the same type of conditions was selected for
the numerical simulations carried out in [19]. Homogeneous Neumann boundary
conditions for the pressure $p$ are the most natural choice since we impose
zero flux of the transported material through $\partial\Omega$, where the
material flux is given by ${\vec{q}=(r\mathbb{I}+\mathbb{C})\nabla p}$.

At the end, the system is closed by defining an initial condition
$\mathbb{C}_{0}\geq 0$ for the conductivity $\mathbb{C}$,

$\mathbb{C}(t=0,\vec{x})=\mathbb{C}_{0}(\vec{x})\quad\text{ in }\Omega.$ (5)

The system is derived as gradient flow of an energy functional, consisting of
a diffusive term, an activation term and a metabolic cost term. The expression
for the functional $\mathcal{E}[\mathbb{C}]$ reads

$\mathcal{E}[\mathbb{C}]:=\int_{\Omega}\left(\frac{D^{2}}{2}|\nabla\mathbb{C}|^{2}+c^{2}\nabla
p[\mathbb{C}]\cdot\mathbb{P}[\mathbb{C}]\nabla
p[\mathbb{C}]+\frac{\alpha}{\gamma}|\mathbb{C}|^{\gamma}\right)\mathrm{d}\vec{x},$
(6)

where $\mathbb{P}[\mathbb{C}]=r\mathbb{I}+\mathbb{C}$, and, in Section 4, we
show that our numerical method provides energy decay, in agreement with the
mathematical model.

Before presenting the variational formulation of Eqs.(1–2), we discuss the
functional analysis setting in which we are going to work.

The variable $p$ is the solution of the diffusion equation with Neumann
boundary conditions, consequently we impose a null-mean condition that
guarantees uniqueness of the solution. Therefore, we introduce the space $Q$
defined as

$Q=\biggl{\\{}q\in H^{1}(\Omega):\int_{\Omega}q\,{\rm d}\Omega=0\biggr{\\}}.$
(7)

For the conductivity tensor $\mathbb{C}$, we introduce the following tensor-
valued space, which takes into account the symmetry

$\mathbf{V}=\biggl{\\{}\mathbb{B}\in[H^{1}(\Omega)]^{2\times
2}:\mathbb{B}=\mathbb{B}^{\top}\biggr{\\}}.$ (8)

Multiplying Eqs.(1–2) by test functions $q\in Q$ and
$\mathbb{B}\in\mathbf{V}$, respectively, and integrating, we obtain

$\displaystyle-\int_{\Omega}\nabla\cdot\left((r\mathbb{I}+\mathbb{C})\nabla
p\right)\,q\,{\rm d}\Omega=\int_{\Omega}S\,q\,{\rm d}\Omega$
$\displaystyle\forall q\in Q$ (9)
$\displaystyle\int_{\Omega}\frac{\partial\mathbb{C}}{\partial
t}\,\mathbb{B}\,{\rm
d}\Omega-D^{2}\int_{\Omega}(\Delta\mathbb{C})\,\mathbb{B}\,{\rm d}\Omega$ (10)
$\displaystyle\hskip 42.67912pt-c^{2}\int_{\Omega}(\nabla p\otimes\nabla
p)\,\mathbb{B}\,{\rm
d}\Omega+\alpha\int_{\Omega}\|\mathbb{C}\|^{\gamma-2}\mathbb{C}\,\mathbb{B}\,{\rm
d}\Omega=0$ $\displaystyle\forall\mathbb{B}\in\mathbf{V}.$

Integrating by parts and taking into account the boundary conditions, we have

$\displaystyle\int_{\Omega}(r\mathbb{I}+\mathbb{C})\nabla p\cdot\nabla q\,{\rm
d}\Omega=\int_{\Omega}S\,q\,{\rm d}\Omega$ $\displaystyle\forall q\in Q$ (11)
$\displaystyle\int_{\Omega}\frac{\partial\mathbb{C}}{\partial
t}\,\mathbb{B}\,{\rm
d}\Omega+D^{2}\int_{\Omega}\nabla\mathbb{C}:\nabla\mathbb{B}\,{\rm d}\Omega$
(12) $\displaystyle\hskip 42.67912pt-c^{2}\int_{\Omega}(\nabla p\otimes\nabla
p)\,\mathbb{B}\,{\rm
d}\Omega+\alpha\int_{\Omega}\|\mathbb{C}\|^{\gamma-2}\mathbb{C}\,\mathbb{B}\,{\rm
d}\Omega=0$ $\displaystyle\forall\mathbb{B}\in\mathbf{V}.$

Finally, denoting by $(\cdot,\cdot)$ the scalar product in $L^{2}(\Omega)$,
our problem in variational formulation reads as follows.

######  Problem 1.

Given $S\in L^{2}(\Omega)$ and $\mathbb{C}_{0}\in\mathbf{V}$, such that
$\nabla\mathbb{C}^{0}(t,\vec{x})\cdot\nu=0$, find $p(t)\in Q$ and
$\mathbb{C}(t)\in\mathbf{V}$ such that, for almost every $t\in(0,T)$, it holds

$\displaystyle\big{(}(r\mathbb{I}+\mathbb{C})\nabla p,\nabla
q\big{)}=\big{(}S,q\big{)}$ $\displaystyle\forall q\in Q$ (13)
$\displaystyle\bigg{(}\frac{\partial\mathbb{C}(t)}{\partial
t},\mathbb{B}\bigg{)}+D^{2}\big{(}\nabla\mathbb{C}(t),\nabla\mathbb{B}\big{)}$
(14) $\displaystyle\hskip 42.67912pt-c^{2}\big{(}\nabla p(t)\otimes\nabla
p(t),\mathbb{B}\big{)}+\alpha\big{(}\|\mathbb{C}(t)\|^{\gamma-2}\mathbb{C}(t),\mathbb{B}\big{)}=0$
$\displaystyle\forall\mathbb{B}\in\mathbf{V}$ $\displaystyle\nabla
p(t,\vec{x})\cdot\nu=0$ $\displaystyle\text{on }\partial\Omega$ (15)
$\displaystyle\nabla\mathbb{C}(t,\vec{x})\cdot\nu=0$ $\displaystyle\text{on
}\partial\Omega$ (16)
$\displaystyle\mathbb{C}(0,\vec{x})=\mathbb{C}_{0}(\vec{x})$
$\displaystyle\text{in }\Omega.$ (17)

##  3 Discretization

We discretize Problem 1 in space by finite elements. Let us consider a mesh
$\mathcal{T}_{h}$ of $\Omega$ with mesh size $h$, so that we can introduce two
finite dimensional subspaces $Q_{h}\subset Q$ and
$\mathbf{V}_{h}\subset\mathbf{V}$. The discrete problems is the following one.

######  Problem 2.

Given $S\in L^{2}(\Omega)$ and $\mathbb{C}_{0,h}\in\mathbf{V}_{h}$, such that
$\nabla\mathbb{C}_{h}^{0}(t,\vec{x})\cdot\nu=0$, find $p_{h}(t)\in Q_{h}$ and
$\mathbb{C}_{h}(t)\in\mathbf{V}_{h}$ such that, for almost every $t\in(0,T)$,
it holds

$\displaystyle\big{(}(r\mathbb{I}+\mathbb{C}_{h})\nabla p_{h},\nabla
q_{h}\big{)}=\big{(}S,q_{h}\big{)}$ $\displaystyle\forall q_{h}\in Q_{h}$ (18)
$\displaystyle\bigg{(}\frac{\partial\mathbb{C}_{h}(t)}{\partial
t},\mathbb{B}_{h}\bigg{)}+D^{2}\big{(}\nabla\mathbb{C}_{h}(t),\nabla\mathbb{B}_{h}\big{)}-c^{2}\big{(}\nabla
p_{h}(t)\otimes\nabla p_{h}(t),\mathbb{B}_{h}\big{)}$ (19)
$\displaystyle\hskip
128.0374pt+\alpha\big{(}\|\mathbb{C}_{h}(t)\|^{\gamma-2}\mathbb{C}_{h}(t),\mathbb{B}_{h}\big{)}=0$
$\displaystyle\forall\mathbb{B}_{h}\in\mathbf{V}_{h}$ $\displaystyle\nabla
p_{h}(t,\vec{x})\cdot\nu=0$ $\displaystyle\text{on }\partial\Omega$ (20)
$\displaystyle\nabla\mathbb{C}_{h}(t,\vec{x})\cdot\nu=0$
$\displaystyle\text{on }\partial\Omega$ (21)
$\displaystyle\mathbb{C}_{h}(0,\vec{x})=\mathbb{C}_{0,h}(\vec{x})$
$\displaystyle\text{in }\Omega.$ (22)

We now discuss the discretization in time. We consider a uniform partition of
the time interval $[0,T]$ in $N\in\mathbb{N}$ parts of length $\Delta t$
denoting the nodes by $t^{n}=n\Delta t$ and
$\mathbb{C}_{h}^{n}\approx\mathbb{C}_{h}(t^{n})$. A first order semi-implicit
approximation of the time derivative of $\mathbb{C}_{h}$ is obtained with
Implicit Euler scheme

$\frac{\partial\mathbb{C}_{h}(t_{n+1})}{\partial
t}\approx\frac{\mathbb{C}_{h}^{n+1}-\mathbb{C}_{h}^{n}}{\Delta t}$

so that Eq. (19) can be rewritten as

$\displaystyle\left(\,\mathbb{C}_{h}^{n+1},\mathbb{B}_{h}\,\right)+\Delta
t\,D^{2}\left(\,\nabla\mathbb{C}_{h}^{n+1},\nabla\mathbb{B}_{h}\,\right)$
$\displaystyle-\Delta t\,c^{2}\left(\,\nabla p_{h}^{n}\otimes\nabla
p_{h}^{n},\mathbb{B}_{h}\,\right)$ $\displaystyle+\Delta
t\,\alpha\left(\,\left(\|\mathbb{C}_{h}^{n}\|+\varepsilon\right)^{\gamma-2}\mathbb{C}_{h}^{n+1},\mathbb{B}_{h}\right)=\left(\,\mathbb{C}_{h}^{n},\mathbb{B}_{h}\,\right),$

where we add a regularization parameter $\varepsilon>0$, to prevent the
instability coming from the division by zero in the metabolic term when
$\gamma<1$ and $\mathbb{C}$ is close to zero in the region without
ramifications (see [7, 6] for more details). In order to solve Problem 2, we
adopt a partitioned scheme; in particular, the nonlinear term $\left(\,\nabla
p_{h}^{n}\otimes\nabla p_{h}^{n},\mathbb{B}_{h}\,\right)$ is computed by means
of $p_{h}^{n}$, which is obtained by solving

$\left(\,\left(r\mathbb{I}+\mathbb{C}_{h}^{n}\right)\nabla p_{h}^{n},\nabla
q_{h}\,\right)=\left(\,S,q_{h}\right)$

at the previous time step. Moreover, we remark that the nonlinear metabolic
term

$\left(\,\left(\|\mathbb{C}_{h}^{n}\|+\varepsilon\right)^{\gamma-2}\mathbb{C}_{h}^{n+1},\mathbb{B}_{h}\right)$

is treated explicitly.

In this setting, our time advancing scheme is of order one. Anyway, a second
order scheme cannot be easily obtained because of the nonlinearity of the
system. Fully second order accuracy could be recovered by adopting, for
example, a predicted pressure at time $t^{n+1/2}$, as suggested in [20], or by
treating the metabolic term by means of a semi-implicit IMEX approach (see,
e.g., [9]).

##  4 Results

In this section we perform several simulations with different choices of the
parameters involved in the discrete model. In particular, we see that the
results agree with the previous work obtained with finite difference schemes
[6, 7].

The domain we consider is $\Omega=[0,1]\times[0,1]$ and the initial conditions
and source term are defined in Eqs.(23–24). In Table 1 we collect the values
of all the parameters related to our tests.

| $\alpha$ | $c^{2}$ | $D^{2}$ | $\varepsilon$ | $\gamma$ | $r$ | $T$
---|---|---|---|---|---|---|---
Test-A1: | 0.75 | 1 | $10^{-1}$ | $10^{-1}$ | 1.25 | $10^{-1}$ | 1
Test-A2: | 0.75 | 1 | $10^{-1}$ | $10^{-1}$ | 0.75 | $10^{-1}$ | 1
Test-A3: | 0.75 | 1 | $10^{-4}$ | $10^{-1}$ | 0.75 | $10^{-1}$ | 1
Test-Al1: | 1.5 | 25 | $10^{-3}$ | $10^{-3}$ | 0.75 | $10^{-3}$ | 10
Test-Al2: | 1.5 | 25 | $5\cdot 10^{-4}$ | $10^{-3}$ | 0.75 | $10^{-3}$ | 10
Test-DD1: | 0.75 | 25 | $2\cdot 10^{-3}$ | $10^{-3}$ | 0.75 | $10^{-3}$ | 10
Test-DD2: | 0.75 | 25 | $10^{-3}$ | $10^{-3}$ | 0.75 | $10^{-3}$ | 10
Test-DD3: | 0.75 | 25 | $5\cdot 10^{-4}$ | $10^{-3}$ | 0.75 | $10^{-3}$ | 10
Test-r1: | 0.75 | 25 | $10^{-3}$ | $10^{-3}$ | 0.75 | $10^{-1}$ | 10
Test-r2: | 0.75 | 25 | $10^{-3}$ | $10^{-3}$ | 0.75 | $10^{-2}$ | 10
Table 1: Parameters involved in the tests shown in Section 4. The first three
rows are related to the accuracy tests: their results are summarized in Table
2.

We define the initial condition for each component of the tensor
$\mathbb{C}_{0}(\vec{x})$ and the source function $S(\vec{x})$ as follows:

$\displaystyle\mathbb{C}_{0}^{(1,1)}(\vec{x})=1,\quad\mathbb{C}_{0}^{(1,2)}(\vec{x})=0,\quad\mathbb{C}_{0}^{(2,2)}(\vec{x})=1,$
(23) $\displaystyle S(\vec{x})=E(\vec{x})-\overline{E},\quad
E(\vec{x})=\exp(-\sigma(\vec{x}-\vec{x}_{0})^{2}),$ (24)
$\displaystyle\sigma=500,\quad\vec{x}_{0}=(0.25,0.25),$ (25)

where $\overline{E}$ is the average of $E$ over the domain $\Omega$. We remark
that the initial conditions and sink-source function $S$ are identical to
those already used in [6]. Moreover, Test-DD3 is reproduced from the same
paper.

The discretization of the problem is performed by considering continuous
piecewise linear finite elements for both $p$ and $\mathbb{C}$ on uniform
triangulations $\mathcal{T}_{h}$ of $\Omega$; in particular, we have

$\displaystyle Q_{h}=\\{q\in
H^{1}(\Omega):\,q_{|K}\in\mathcal{P}_{1}(K)\,\,\forall K\in\mathcal{T}_{h}\\}$
(26) $\displaystyle\mathbf{V}_{h}=\\{\mathbb{B}\in[H^{1}(\Omega)]^{2\times
2}:\,\mathbb{B}_{|K}\in[\mathcal{P}_{1}(K)]^{2\times 2}\,\,\forall
K\in\mathcal{T}_{h}\\}.$

###  4.1 Accuracy tests

First, we check the accuracy in time of the adopted scheme. In Table 2, we see
the error for the norm of the conductivity tensor, calculated with Richardson
extrapolation (see, e.g., [23]) by refining the mesh size $h$. We remark that
the time step $dt$ depends on $h$, since we have $dt=h$. The chosen parameters
are collected in Test-A1–Test-A3.

$h$ | error | order
---|---|---
1/20 | – | –
1/40 | 0.0117878 | –
1/80 | 0.00612581 | 0.94
1/160 | 0.00312414 | 0.97
1/320 | 0.00157781 | 0.99

$h$ | error | order
---|---|---
1/20 | – | –
1/40 | 0.0136681 | –
1/80 | 0.00721071 | 0.92
1/160 | 0.00370715 | 0.96
1/320 | 0.00188007 | 0.98

$h$ | error | order
---|---|---
1/20 | – | –
1/40 | 0.0136681 | –
1/80 | 0.00721071 | 0.92
1/160 | 0.00370715 | 0.96
1/320 | 0.00188007 | 0.98

Table 2: Accuracy study by computing the relative error in $L^{2}-$norm
between two solutions with two different mesh sizes. In order to compute the
error for $\|\mathbb{C}\|$ on $\mathcal{T}_{h}$, we consider as reference
solution $\|\mathbb{C}\|$ computed in the previous test on
$\mathcal{T}_{h/2}$. The problem has been solved with parameters defined in
Test-A1–Test-A3.

###  4.2 Numerical tests

In this subsection, we show and analyze the results of our tests in terms of
$\|\mathbb{C}\|$: in particular, we analyze its snapshots when the steady
state is reached. When this happens, the energy functional remains constant in
time. Another possible indicator may be the numerical time derivative of the
tensor, that should be approximately close to zero.

Again, the parameters of the tests are collected in Table 1, the initial
conditions and source term in Eqs.(23–24). Each test has been performed by
considering a triangulation with size $h=1/600$ and $dt=0.01$.

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures_FEM/r2.pdf}
\put(0.5,0.1){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\rightarrow\,x$}}
\put(0.01,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\uparrow$}}
\put(0.01,7.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$y$}}
\put(9.0,43.0){{Test-}r1: $\|\mathbb{C}\|$} \end{overpic}

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures_FEM/r1.pdf}
\put(0.5,0.1){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\rightarrow\,x$}}
\put(0.01,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\uparrow$}}
\put(0.01,7.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$y$}}
\put(9.0,43.0){{Test-}r2: $\|\mathbb{C}\|$} \end{overpic}

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures_FEM/DD1.pdf}
\put(0.5,0.1){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\rightarrow\,x$}}
\put(0.01,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\uparrow$}}
\put(0.01,7.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$y$}}
\put(7.0,43.0){{Test-DD1}: $\|\mathbb{C}\|$} \end{overpic}

Figure 1: In this figure we show three results considering different values of
the background permeability $r$. From left to right we have the results of
Test-r1 with $r=10^{-1}$, Test-r2 with $r=10^{-2}$ and Test-DD1 with
$r=10^{-3}$. We plot two different quantities for each test: in the first row,
the norm of the solution $\mathbb{C}$ at final time; in the second row, the
time dependent plot of the energy $\mathcal{E}[\mathbb{C}]$.

In Fig. 1 we show the results with different values of the background
permeability $r$. The tests we consider are Test-r1, Test-r2 and Test-DD1. In
terms of $\|\mathbb{C}\|$, the most interesting case is Test-DD1, with
$r=10^{-3}$, because a more detailed pattern appears, in agreement with the
results in [6, 7]. When $r\to 0$, the energy plots show some oscillations
while approaching the steady state. Moreover, as showed in [6], the iteration
matrix for the pressure becomes ill conditioned when we decrease the parameter
$r$, and this was causing the deterioration of the symmetry of the solution.
In this case, not only we see the lack of symmetry, but also the non monotonic
behaviour of the energy in time.

In Fig. 2, we show the results with different values of the diffusivity
$D^{2}$. The tests we consider are Test-DD1 with $D^{2}=2\cdot 10^{-3}$, Test-
DD2 with $D^{2}=10^{-3}$ and Test-DD3 $D^{2}=5\cdot 10^{-4}$. From the plot of
$\|\mathbb{C}\|$, we can see that the ramifications become thinner and thinner
when decreasing the diffusivity, while the plots of the energy show some
oscillations when approaching the steady state, since we are considering
$r=10^{-3}$. In particular, the pattern emerged from Test-DD3 is in line with
the results presented in [6]. On the other hand, the energy shows some
oscillations, probably caused by the considered time step.

In Fig. 3, we show the results with different values of the parameters
$\alpha$ and $D^{2}$. From left to right we have the results of Test-Al1 with
$\alpha=0.75$ and $D^{2}=10^{-3}$, Test-Al2 with $\alpha=1.5$ and
$D^{2}=5\cdot 10^{-4}$ and Test-DD3 with $\alpha=0.75$ and $D^{2}=5\cdot
10^{-4}$. We see more ramifications for $\alpha=1.5$, but the respective
energy plots show a more irregular behavior in time since the oscillations are
higher.

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures_FEM/DD1.pdf}
\put(0.5,0.1){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\rightarrow\,x$}}
\put(0.01,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\uparrow$}}
\put(0.01,7.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$y$}}
\put(7.0,43.0){{Test-DD1}: $\|\mathbb{C}\|$} \end{overpic}

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures_FEM/DD2.pdf}
\put(0.5,0.1){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\rightarrow\,x$}}
\put(0.01,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\uparrow$}}
\put(0.01,7.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$y$}}
\put(7.0,43.0){{Test-DD2}: $\|\mathbb{C}\|$} \end{overpic}

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures_FEM/DD3.pdf}
\put(0.5,0.1){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\rightarrow\,x$}}
\put(0.01,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\uparrow$}}
\put(0.01,7.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$y$}}
\put(7.0,43.0){{Test-DD3}: $\|\mathbb{C}\|$} \end{overpic}

Figure 2: In this figure we show three results considering different values of
the squared diffusivity $D^{2}$. From left to right we have the results of
Test-DD1 with $D^{2}=2\cdot 10^{-3}$, Test-DD2 with $D^{2}=10^{-3}$ and Test-
DD3 with $D^{2}=5\cdot 10^{-4}$. Same format as Fig. 1.

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures_FEM/AL1.pdf}
\put(0.5,0.1){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\rightarrow\,x$}}
\put(0.01,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\uparrow$}}
\put(0.01,7.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$y$}}
\put(7.0,43.0){{Test-AL1}: $\|\mathbb{C}\|$} \end{overpic}

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures_FEM/AL4.pdf}
\put(0.5,0.1){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\rightarrow\,x$}}
\put(0.01,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\uparrow$}}
\put(0.01,7.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$y$}}
\put(7.0,43.0){{Test-AL2}: $\|\mathbb{C}\|$} \end{overpic}

\begin{overpic}[abs,width=390.25534pt,unit=1mm,scale={.25}]{Figures_FEM/DD3.pdf}
\put(0.5,0.1){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\rightarrow\,x$}}
\put(0.01,2.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$\uparrow$}}
\put(0.01,7.0){{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}$y$}}
\put(7.0,43.0){{Test-DD3}: $\|\mathbb{C}\|$} \end{overpic}

Figure 3: In this figure we show three results considering different values
for $\alpha$ and $D^{2}$. From left to right we have the results of Test-Al1
with $\alpha=0.75$ and $D^{2}=10^{-3}$, Test-Al2 with $\alpha=1.5$ and
$D^{2}=5\cdot 10^{-4}$ and Test-DD3 with $\alpha=0.75$ and $D^{2}=5\cdot
10^{-4}$. Same format as Fig. 1.

##  5 Conclusions

In this paper we studied the Cai-Hu model, an elliptic-parabolic system to
describe the formation of biological networks, and, in particular, of leaf
venation networks.

We used a finite element method in space associated with a semi-implicit
discretization in time, due to the nonlinearities of the model. The use of the
library FreeFEM makes the presented space discretization very easy to
implement. The convergence rate in time and space is computed numerically and
confirms the first order accuracy.

We studied the evolution of the energy, that decays in time, but it is non-
monotonic when the parameter $r\to 0$. This depends on the non-solvability of
the Poisson equation when $r=0$.

We remark that a big variety of parameters is involved in the system and each
of them seems to influence the solution in different ways. A deeper
investigation should focus on their roles in the equations.

Our future research will investigate how to improve the numerical method by
computing the solutions in a monolithic way: for instance, this can be done by
considering a mixed formulation. In this case, in order to treat the nonlinear
terms in the system, an iterative method can be adopted, such as a Newton or a
fixed point iterator.

The ADI scheme, already applied in [6, 7] in its symmetric version, would
improve the time accuracy of the scheme. High order in time can be assured
also with the use of implicit-explicit (IMEX) Runge-Kutta schemes. These
methods are effective also in presence of stiffness in the system; in
particular, the family of Stiffly Accurate (SA) schemes have been created for
this purpose.

Moreover, the use of mesh adaptivity could accurately track the conductivity
flow by refining the mesh in those regions of interest where ramifications
start appearing, deeply examining the multiple scales hidden in the problem.
These features are already under investigation, together with the design of a
parallel solver able to accelerate the solution of the linear system at each
time step. In this paper we were able to solve the equations for the three
components of the tensor, plus one for the pressure, choosing a mesh size
equal to 1/600. This resolution is already computationally expensive for a
serial algorithm.

## Acknowledgments

F. C. is partially supported by IMATI/CNR.

## References

  * [1] Malinowski, R. (2013). _Understanding of leaf development – the science of complexity_. _Plants_ , 2(3), 396–415. Publisher: Multidisciplinary Digital Publishing Institute.
  * [2] Sedmera, D. (2011). _Function and form in the developing cardiovascular system_. _Cardiovascular research_ , 91(2), 252–259. Publisher: Oxford University Press.
  * [3] Eichmann, A., Le Noble, F., Autiero, M., & Carmeliet, P. (2005). _Guidance of vascular and neural network formation_. _Current opinion in neurobiology_ , 15(1), 108–115. Publisher: Elsevier.
  * [4] Michel, O., & Biondi, J. (1995). _Morphogenesis of neural networks_. _Neural Processing Letters_ , 2(1), 9–12. Publisher: Springer.
  * [5] Hu, D., Cai, D., & Rangan, A. V. (2012). _Blood vessel adaptation with fluctuations in capillary flow distribution_. _PLoS ONE_. Publisher: Public Library of Science San Francisco, USA.
  * [6] C. Astuto, D. Boffi, J. Haskovec, P. Markowich and G. Russo, Asymmetry and condition number of an elliptic-parabolic system for biological network formation. 2023\. doi: 10.1007/s42967-023-00297-3.
  * [7] C. Astuto, D. Boffi, J. Haskovec, P. Markowich and G. Russo, Comparison of Two Aspects of a PDE Model for Biological Network Formation. In: Mathematical and Computational Applications 27.5 (2022). issn: 2297-8747. doi: 10.3390/mca27050087, https://www.mdpi.com/2297-8747/27/5/87.
  * [8] Nicola Bellomo, Pierre Degond, and Eitan Tadmor. Active Particles, Volume 1: Advances in Theory, Models, and Applications. Birkhäuser, 2017.
  * [9] S. Boscarino, F. Filbet, and G. Russo. High Order Semi-implicit Schemes for Time Dependent Partial Differential Equations. In: Journal of Scientific Computing (2016).
  * [10] Philippe G. Ciarlet. The finite element method for elliptic problems. SIAM, 2002.
  * [11] I. Danaila and F. Hecht. A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose–Einstein condensates. In: Journal of Computational Physics 229.19 (2010), pp. 6946–6960. issn: 0021-9991. doi: 10.1016/j.jcp.2010.05.032.
  * [12] K. Ohtsuka, F. Hecht, O. Pironneau. FreeFem++ Manual.
  * [13] J. Haskovec, P. Markowich, and B. Perthame. Mathematical Analysis of a PDE System for Biological Network Formation. In: Communications in Partial Differential Equations 40.5 (2015), pp. 918–956. doi: 10.1080/03605302.2014.968792.
  * [14] J. Haskovec, P. Markowich, and S. Portaro. Emergence of biological transportation networks as a self-regulated process. In: Discrete Contin. Dyn. Syst. 43.3-4 (2023), pp. 1499–1515. issn: 1078-0947. doi: 10.3934/dcds.2022159.
  * [15] J. Haskovec, P. Markowich, B. Perthame and M. Schlottbom Notes on a PDE system for biological network formation. In: Nonlinear Analysis 138 (2016), pp. 127–155.
  * [16] F. Hecht. New development in FreeFem++. In: J. Numer. Math. 20.3-4 (2012), pp. 251–265. issn: 1570-2820.
  * [17] Dan Hu. Optimization, adaptation, and initialization of biological transport networks. In: Notes from lecture 1 (2013), pp. 3–1.
  * [18] Dan Hu and David Cai. Adaptation and optimization of biological transport networks. In: Physical review letters 111.13 (2013), p. 138701.
  * [19] Dan Hu and David Cai. An optimization principle for initiation and adaptation of biological transport networks. In: Communications in Mathematical Sciences 17.5 (2019), pp. 1427–1436.
  * [20] A. Raudino, A. Grassi, G. Lombardo, G. Russo, C. Astuto and M. Corti, Anomalous Sorption Kinetics of Self-Interacting Particles by a Spherical Trap. In: Communications in Computational Physics 31.3 (2022), pp. 707–738. issn: 1991-7120. doi: 10.4208/cicp.OA-2021-0101.
  * [21] N. Raza, S. Sial, and S. Siddiqi, Approximating time evolution related to Ginzburg–Landau functionals via Sobolev gradient methods in a finite-element setting. In: Journal of Computational Physics 229.5 (2010), pp. 1621–1625.
  * [22] N. Raza, S. Sial, and S. Siddiqi, Shahid S and T. Lookman, Energy minimization related to the nonlinear Schrödinger equation. In: Journal of Computational Physics 228.7 (2009), pp. 2572–2577.
  * [23] Lewis Fry Richardson. IX. The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. In: Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 210.459-470 (1911), pp. 307–357.