# A directional total variation minimization algorithm for isotropic
resolution in
digital breast tomosynthesis

Emil Y. Sidky, Xiangyi Wu, Xiaoyu Duan, Hailing Huang, Wei Zhao, Leo Y. Zhang,
John Paul Phillips, Zheng Zhang, Buxin Chen, Dan Xia, Ingrid S. Reiser,
Xiaochuan Pan E. Y. Sidky, J. P. Phillips, Z. Zhang, B. Chen, D. Xia, I. S.
Reiser, and X. Pan are with the Department of Radiology MC-2026, The
University of Chicago, 5841\. S. Maryland Ave., Chicago, IL, 60637.X. Wu, X.
Duan, H. Huang, and W. Zhao are with the Department of Radiology, Renaissance
School of Medicine, Stony Brook University, 101 Nicolls Road, Health Sciences
Center, Level 4, Stony Brook, NY, 11794.L. Y. Zhang is with the Department of
Mathematics, University of Washington, Box 354350, C-138 Padelford, Seattle,
WA 98195-4350.

###### Abstract

An optimization-based image reconstruction algorithm is developed for contrast
enhanced digital breast tomosynthesis (DBT) using dual-energy scanning. The
algorithm minimizes directional total variation (TV) with a data discrepancy
and non-negativity constraints. Iodinated contrast agent (ICA) imaging is
performed by reconstructing images from dual-energy DBT data followed by
weighted subtraction. Physical DBT data is acquired with a Siemens Mammomat
scanner of a structured breast phantom with ICA inserts. Results are shown for
both directional TV minimization and filtered back-projection for reference.
It is seen that directional TV is able to substantially reduce depth blur for
the ICA objects.

###### Index Terms:

Limited angular-range reconstruction, digital breast tomosynthesis, dual-
energy, contrast-agent imaging

##  I Introduction

Over the past two decades, X-ray based tomosynthesis has been an emerging
imaging modality that has been making great strides in clinical adoption; most
notably, digital breast tomosynthesis (DBT) is becoming the workhorse imaging
device for breast cancer screening. Although DBT yields partially tomographic
volumes, there is still a large gap between it and X-ray computed tomography
in terms of achieving isotropically high resolution and in the ability to
achieve quantitive imaging for shape and volume features of tissues within the
scanned subject. Tomosynthesis imaging is subject to significant depth blur,
and this blurring is proportional to the size of the imaged structures.
Because of this depth blurring it is not possible to outline the boundaries of
structures and estimate their volume. Moreover, the gray scale values in the
image do not directly correspond to a physical density.

In this work, an optimization-based algorithm is developed that may enable
quantitive imaging with limited angular range scanning such as what is
performed in a tomosynthesis acquisition. In order to accomplish this, we
exploit sparsity in the gradient magnitude of the scanned subject. For sparse-
view imaging over a full scanning angular range, the use of total variation
(TV) minimization has proven effective [1]. For limited angular range
scanning, however, directional TV constraints show promise [2]. The use of
directional TV constraints can be shown to be equivalent to directional TV
minimization where the gradient magnitudes in different directions have
different penalty weights [3]. The work presented in Ref. [3] showed only
simulation studies, and in this work we apply this algorithm to physical DBT
transmission data.

The application of interest is contrast enhanced (CE) DBT [4], where the goal
is to obtain a quantitative image of the Iodinated contrast agent (ICA)
distribution. If such a system is successfully developed it could enable the
use of CE-DBT in assessment of breast cancer therapy effectiveness. For this
preliminary work, DBT data is obtained for structured breast phantom with ICA
inserts. Volumes are reconstructed from scans at 30 and 49 kV, and the
weighted image subtraction method is used to isolate the ICA distribution.

##  II Methods

###  II-A Image reconstruction

In order to obtain reconstructed DBT volumes that have nearly isotropic
resolution for the ICA distribution, we exploit the fact that this
distribution is localized with a fairly simple structure and high contrast. In
the algorithm implementation directional gradient sparsity regularization is
employed and the resolution is reduced by a factor of eight from the native
resolution of the DBT system’s detector.

The image reconstruction algorithm is developed for the Siemens Mammomat
scanner that acquires 25 projections over a 50 degree scanning arc and the
detector resolution is 85 microns. The transmission data are processed with
the negative logarithm after dividing by a flood-field scan and the
relationship between this processed data and the DBT volume is assumed to be
standard linear X-ray projection as described by

$g=Xf,$

where $f$ is a discrete representation of the DBT volume, $X$ is a matrix
encoding X-ray projection, and $g$ is the processed projection data.

For sparsity regularization, image reconstruction is formulated using an
optimization model where weighted image directional total variation (TV) is
minimized subject to a data-fidelity constraint. Furthermore, an $\ell_{1}$
penalty term is included that encourages pixel sparsity and helps to confine
the image to the true object support. Finally, the data fidelity term compares
ramp-filtered data with ramp-filtered estimated data, which has two important
effects: this filtering deemphasizes low spatial-frequency discrepancy in the
sinogram estimation and it also performs preconditioning.

$\displaystyle
f=\operatorname*{arg\,min}_{f^{\prime}}\bigl{\\{}\alpha\|\nabla_{x}f^{\prime}\|_{1}+\alpha\|\nabla_{y}f^{\prime}\|_{1}+(2-\alpha)\|\nabla_{z}f^{\prime}\|\bigr{.}$
$\displaystyle\bigl{.}+\beta\|f^{\prime}\|_{1}\text{ such that
}\|R(g-Xf^{\prime})\|_{2}\leq\epsilon\text{ and }f^{\prime}\geq 0\bigr{\\}}.$

The $x$ and $y$ coordinates are in-plane with $y$ being along the direction of
travel for the X-ray source, and the $z$ direction is the depth direction;
i.e. perpendicular to the X-ray detector. The symbol $\nabla$ represents
finite differencing along the direction indicated by the subscript; $\alpha$
is a weighting parameter that should be chosen between 0 and 2; $\beta$ is the
weighting parameter for the $\ell_{1}$ penalty, $R$ represents ramp filtering
and $\epsilon$ is a data error tolerance parameter. The posed optimization
problem is convex and it can be solved with the primal-dual algorithm of
Chambolle and Pock [5].

Figure 1: CIRS phantom with solid ICA inserts. The background is hetergeneous
with 50% glandularity.

When $\alpha=1$, the gradient sparsity regularization term is equivalent to TV
and it has been shown that constrained TV minimization is effective for
accurate image reconstruction for sparse-view acquisition when the scanning
arc is complete. Recently, we have shown that the anisotropic weighting of the
directional TV terms can improve image reconstruction accuracy for limited
scanning angular ranges such as those used for DBT. This in combination with
use of a coarse volume grid allows for nearly isotropic resolution in the
reconstructed DBT images using directional TV minimization.

###  II-B Image post-processing

Even when the transmission data preprocessing accurately estimates the image
sinogram, Reconstructed DBT images are susceptible to low spatial-frequency
artifact due to the limited angular range of the scan. In order to facilitate
the dual-energy image subtraction the low spatial-frequency artifacts are
estimated and removed by normalization. Specifically, the low spatial-
frequency image component is estimated by robust polynomial fitting

$b=\operatorname*{arg\,min}_{c}\bigl{\|}m\cdot(s-Pc)\bigr{\|}_{1},$ (1)

where $s$ is an in-plane slice of the image $f$; $P$ is a matrix of low-order
polynomials; $c$ is a vector of coefficients; $m$ is an indicator function
which is one on the object support and zero otherwise, where the object
support is estimated by thresholding the image slice, $s$; and $b$ is the
estimated background for the slice $s$. The use of the $\ell_{1}$-norm for the
objective function encourages the difference image in its argument to be
sparse, and this design should preserve important small structures in the
image when the background drift is processed away.

Figure 2: Raw and processed in-plane images containing all ICA objects for the
directional TV minimization images for the 30 kV (Left) and 49 kV (Right) DBT
scans. Image processing is performed by polynomial fitting and normalization
described in Eqs. (1) and (2), respectively. The gray scale for the raw images
is selected to span the range of background variation, and the gray scale for
the normalized images is selected to be $\pm$ 30% of background.

The polynomial fitting is performed slice by slice with the result stacked
into the background image $f_{b}$. Likewise, the mask image $f_{m}$ is
assembled by stacking the slice masks. The algorithm for solving this robust
polynomial fitting is iteratively reweighted least-squares. The final
displayed image, $f_{d}$, is computed as follows

$f_{d}=f_{m}\cdot\frac{f}{f_{b}}.$ (2)

After the processed images are computed for low and high kV, a difference
image is computed

$f^{(ICA)}=\left(f^{(HE)}_{d}-1\right)-w\cdot\left(f^{(LE)}_{d}-1\right).$

For this work, the dual energy data is processed in the image domain by
performing a weighted subtraction of low and high kV reconstructed volumes,
where the weights are designed to highlight the ICA distribution. Because the
low spatial-frequency background variation is reduced by normalization, the
background of the processed images center on the value of 1; thus 1 is
subtracted from both images before performing the weighted subtraction. The
weighting parameter $w$ is selected by visual inspection.

Figure 3: In-plane images containing all ICA objects for directional TV
minimization (Left) and FBP (Right). The top and middle rows show the images
corresponding to the 49 and 30 kV acquisitions, respectively. The bottom row
shows the weighted subtraction designed to isolate the ICA distributions. The
gray scales on the left and right are matched and the images are masked to the
object support.

###  II-C DE-DBT transmission data acquisition

The DE-DBT scanning consists of two DBT acquisitions taken at source potential
of 30 and 49 kV using filtration with Al and Ti, respectively. DE-DBT scans of
a 4 cm thick BR3D phantom (Model 020, CIRS Inc. Norfolk, VA) with solid ICA
inserts, see Fig. 1, are used for testing the proposed image reconstruction
algorithm. The ICA inserts vary in size and concentration; the diameters are
2, 3, 5, and 8 mm and the Iodine concentrations are 1, 2, 3, and 5 mg/mL.

##  III Results

For the present implementation of directional TV minimization, the weighting
parameter $\alpha$ is set to 1.75, favoring the in-plane differentiation in
the objective function. The $\ell_{1}$ penalty parameter $\beta$ is set to
0.1, and the data error constraint $\epsilon$ is 0.01 in terms of root mean
square error. The imaging volume is discretized on a coarse grid of size
$188\times 417\times 175$ with cubic voxels 680 microns wide.

Raw slice images from applying directional TV minimization to the DBT phantom
data are shown in the top panels of Fig. 2, and the low spatial-frequency
artifacts are clearly visible. The slice shown cuts through all of the ICA
objects. These objects are easily seen in the 49 kV images, and not the 30kV
images, because the K-edge of Iodine is 33.2 keV. The spatial dependence of
the artifact is different for the low and high kV images, and therefore direct
image subtraction of the raw images will not yield good isolation of the ICA
distribution. For performing the post-processing, we employed an eighth degree
2D polynomial expansion (45 terms), and resulting images are shown in the
bottom row of Fig. 2, where it is seen that the low spatial-frequency
artifacts are greatly reduced.

The results of weighted subtraction for directional TV minimization filtered
back-projection (FBP) are shown in Fig. 3, where it is seen in the bottom row
that it effectively suppresses the background variation and isolates the ICA
objects in the image.

Figure 4: Transverse plane images cutting through the groups of ICA objects of
the same concentration with directional TV minimization and FBP images on the
left and right, respectively. Concentration of the ICA decreases going from
top to bottom. All gray scales are matched for quantitative comparison and the
images are masked to the object support.

Transverse slices, where $y$ is horizontal and $z$ is vertical, for the
weighted subtraction of directional TV minimization and FBP images are
displayed in Fig. 4. While the FBP images show the usual depth blurring of the
ICA objects, it is clear that directional TV minimization can effectively
mitigate this artifact. On the other hand, the directional TV difference
images have slightly stronger background clutter due to imperfect cancellation
of the background glandular structures, which also interferes with the imaging
of the 1mg/mL ICA inserts. This imperfect cancellation is likely due to the
non-linearity of the sparsity regularization.

##  IV Conclusion

In this preliminary study on image reconstruction by directional TV
minimization for CE-DBT, the physical phantom results show promise in
achieving an accurate ICA distribution estimate with nearly isotropic spatial
resolution; it appears that the depth blurring artifact common to DBT imaging
is well controlled with directional TV minimization. Because the dual-energy
processing, in this work, is image-based the gray values in the ICA images are
not quantitive and cannot be directly tied to a concentration value. Future
work will focus on spectral calibration of the DBT transmission data so that
both one-step and two-step [6] quantitative DE-DBT image reconstruction can be
developed. Adapting directional TV minimization to a quantitive DE-DBT is
expected to greatly reduce background clutter in the ICA images from glandular
tissue, and this should also provide accurate ICA concentrations.

## Acknowledgment

This research was supported in part by the National Institute of Biomedical
Imaging and Bioengineering of the National Institutes of Health under award
numbers R01EB023968 and R21CA263660. The CEDBT system and clinical research
was supported financially by Siemens Medical Solutions, USA.

## References

  * [1] E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” _Phys. Med. Biol._ , vol. 53, pp. 4777–4808, 2008.
  * [2] Z. Zhang, B. Chen, D. Xia, E. Y. Sidky, and X. Pan, “Directional-TV algorithm for image reconstruction from limited-angular-range data,” _Medical Image Analysis_ , vol. 70, p. 102030, 2021.
  * [3] L. Y. Zhang, E. Y. Sidky, J. P. Phillips, Z. Zhang, B. Chen, D. Xia, and X. Pan, “Limited scanning arc image reconstruction with weighted anisotropic TV, minimization,” in _17th International Meeting on Fully 3D Image Reconstruction in Radiology and Nuclear Medicine_ , 2023, pp. 491–494.
  * [4] H. Huang, D. A. Scaduto, C. Liu, J. Yang, C. Zhu, K. Rinaldi, J. Eisenberg, J. Liu, M. Hoernig, J. Wicklein, S. Vogt, T. Mertelmeier, P. R. Fisher, and W. Zhao, “Comparison of contrast-enhanced digital mammography and contrast-enhanced digital breast tomosynthesis for lesion assessment,” _J. Med. Imaging_ , vol. 6, p. 031407, 2019.
  * [5] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” _J. Math. Imag. Vis._ , vol. 40, pp. 120–145, 2011.
  * [6] B. Chen, Z. Zhang, D. Xia, E. Y. Sidky, T. Gilat-Schmidt, and X. Pan, “Accurate image reconstruction in dual-energy CT with limited-angular-range data using a two-step method,” _Bioengineering_ , vol. 9, no. 12, p. 775, 2022.


2022

[1,2]Nicasio Barrere

1]PDU Ciencias Físicas, Centro Universitario Regional Este, Universidad de la
República,Uruguay

2]Instituto de Física, Facultad de Ciencias, Universidad de la República,
Uruguay

# Lagrangian mixing of pulsatile flows in constricted tubes

<EMAIL_ADDRESS>Javier Brum  Maximiliano Anzibar  Felipe Rinderknecht
Gustavo Sarasúa  [  Cecilia Cabeza  [

###### Abstract

In this work several lagrangian methods were used to analyze the mixing and
transport processes in an experimental model of a constricted artery under a
pulsatile flow. Upstream Reynolds number $Re$ was changed between 1187 and
1999, while the pulsatile period $T$ was kept fixed at 0.96s. Velocity fields
were acquired using Digital Particle Image Velocimetry (DPIV) for a region of
interest (ROI) located downstream of the constriction. The flow is composed of
a central jet and a recirculation region near the wall where vortex forms and
sheds. To study the mixing processes, finite time Lyapunov exponents (FTLE)
fields and concentration maps were computed. Two lagrangian coherent
structures (LCS) responsible for mixing and transporting fluid were found from
FTLE ridges. A first LCS delimits the trailing edge of the vortex, separating
the flow that enters the ROI between successive periods. A second LCS delimits
the leading edge of the vortex. This LCS concentrates the highest particle
agglomeration, as verified by the concentration maps. Moreover, from the
particle residence time maps (RT) the probability for a fluid particle of
leaving the ROI before one cycle was measured. As $Re$ increases, the
probability of leaving the ROI increases from 0.6 to 0.95. Final position maps
$r{{}_{f}}$ were introduced to evaluate the flow mixing between different
subregions of the ROI. These maps allowed us to compute an exchange index
between subregions, $\overline{\mathrm{EI}}$, which shows the main region
responsible for the mixing increase with $Re$.

Finally by integrating the results of the different lagrangian methods (FTLE,
Concentration maps, RT and $r_{f}$ maps), a comprehensive description of the
mixing and transport of the flow was provided.

###### keywords:

Pulsatile flows, Blood flows, Coherent Structures, Vortex Dynamics

##  1 Introduction

Hemodynamics is a very important research topic for its implications in
cardiovascular diseases Vlachopoulos2011 ; naghavi . One of the most dangerous
diseases in arteries is atherosclerosis atlas_cardiovascular , characterized
by the aggregation of fat and cholesterol which leads to a partial obstruction
of the arterial lumen. An accurate description of blood flow under these
conditions has become key to understand the initiation and progression of such
diseases. In this context, mixing among regions within the artery is known to
play a major role because of its dependence to disturbances in blood flow
Liu2011 ; shadden_applied_lagrangian .

The mechanisms that govern the mixing processes are complex and diverse,
including advection, diffusion, and turbulence. Since the early works of
Winter and Nerem Winter1984 , several authors peacock ; trip ; Brindise2018 ;
Xu2020  have studied the transition to turbulence. For instance, the works of
Jain Jain2019 ; Jain2022  serve as a benchmark in studying the transition to
turbulence in an oscillatory flow through a rigid pipe with stenosis. These
works showed, among other things, that the pulsatile frequency variation does
not modify the critical Reynolds for transition. However, mixing occurs in a
non-turbulent regime, and a comprehensive description of it is still missing.

Lagrangian methods are based on particle tracking and they provide more
insight into the flow topology than other widely used Eulerian fields.
Therefore, they offer a framework to understand mixing processes and are
essential in quantifying mixing and other features of pulsatile flows
shadden_applied_lagrangian . Particularly, finite-time Lyapunov exponents
(FTLE) haller_2000 ; haller_2000_chaos ; shadden-tesis ; shadden_2005  are
used to identify Lagrangian coherent structures (LCS) that act as material
barriers organizing the flow.

FTLE have been successfully applied in hemodynamic models to study steady and
pulsatile flows arzani_2012 ; shadden_applied_lagrangian ; badass_2017 . Vétel
et al. vetel  studied the flow near the bifurcation of a carotid artery,
identifying the formation of structures, such as vortices, through FTLE
fields. Espa et al. espa  have applied these techniques to study the flow
within the left ventricle. Similarly, Arzani and Shadden arzani_2012
calculated FTLE fields from numerical data to study the transport of vortices,
determining the areas in which vortices are detached, and characterizing the
mixing in aortic aneurysms. However, the description of how mixing occurs
could be complemented by other tools such as concentration maps. The
combination of these lagrangian techniques would enable the measurement of
fluid from each structure taking part in the mixing.

Another widely-used Lagrangian method is the distribution of residence time
(RT) and RT maps residence_time_camassa_91 ; goulliart_2011 ; monsen_2002 ;
Rayz_2010 , which represent the time a given particle remains in a specific
region. Several authors have used RT to study blood-flow mixing in the left
ventricle Seo2013 ; Long2013 ; dilabbio_2018 ; badass_2017 . Di Labbio et al.
dilabbio_2018  studied an experimental model of a left ventricle with
regurgitation. The authors observed that this condition increases the RT of
the flow inside the ventricle, with the potential risk of suffering from
arrhythmia or other cardiovascular diseases. Other studies have reported that
RT maps can be used to analyze how long platelets remains in a specific region
Rayz_2010 ; fuster_2005 .

Particle residence time (PRT) is a path-dependent metric which quantifies the
time a fluid parcel spends in a region of interest jeronimo2019 ; jeronimo2020
. However, PRT differs from RT since it tracks fluid tracers experimentally
while RT is based on advecting particles from the eulerian velocity field. The
work of Jeronimo et al. jeronimo2019  studied PRT of a flow in a model of
stenosed artery. This allowed the identification of particles initial region
and subsequent trajectory, categorizing the flow in reverse flow, jet,
recirculating flow and others. However, a map showing the exact final
positions of particles would allow a more accurate description of the mixing.

In the present work, Lagrangian methods were used to study the mixing of a
pulsatile flow in an experimental model of an artery with stenosis. The study
focused on the region downstream of the constriction, analyzing the mixing
processes by means of FTLE fields and concentration maps. The FTLE fields
allowed the identification of the LCS that delimits two regions with well-
differentiated dynamics: the vortices against the wall and the central jet.
Overall, we explained the characteristics of the flow and the mechanism of
mixing using FTLE fields and Concentration maps. Quantitative analysis of
mixing is not easy since it usually depends on each specific application;
however, in this work, we present Lagrangian descriptors capable of measuring
mixing. The final state of the mixing was analyzed through RT maps and $r_{f}$
maps, which show the final position of the particles. Two relevant parameters
were computed from these maps: the probability of flow leaving the ROI and the
$\overline{\mathrm{EI}}$ index, which measures flow exchange between
subregions within the ROI. Finally, a detailed discussion is given on the
results obtained from the FTLE fields, $C$ maps, RT, and $r_{f}$ and how they
complement each other.

The work is organized as follows. In sec. 2 we describe the experimental setup
and the computation of aforementioned lagrangian methods. In sec. 3 we briefly
show the results obtained by each lagrangian method and in sec. 4 we discuss
the results altogether. Finally, in sec. 5 we present the conclusions.

##  2 Materials and methods

###  2.1 Experimental setup

The experimental set up is shown in fig. 1. Pulsatile flow generated by a
programmable pump (PP) was led into a constricted tube (CT). This constricted
tube consisted of a transparent acrylic tube of diameter $D=2.6\pm 0.1$ cm
containing a hollow, cylindrical annular constriction with an internal
diameter $d_{0}=1.6\pm 0.1$ cm, which reduced the area by a factor of 39%
Barrere2020 .

A flow development section (FDS) was connected upstream CT to ensure developed
flow Barrere2020 . The region of interest, ROI, was defined as -0.5$<r/D<$0.5,
0$<z/D<$1.5, see fig1(b). Digital Particle Image Velocimetry (DPIV) technique
was used to acquire the velocity fields piv1 . Distilled water was seeded with
neutrally buoyant 50$\mu$m particles. A 1 W Nd:YAG laser was used to
illuminate a 2 mm-thick section of the tube. Sixteen pulsatile cycles were
acquired at a frame rate of 180Hz using a CMOS camera (Pixelink, PL-B776F).
The velocity field was finally computed using PIVlab software Thielicke2019
with $32\times 32$ pixel2 windows and an overlap of 8 pixels in both
directions.

Upstream Reynolds number was defined as $Re=Dv_{u}/\nu$, where $v_{u}$ is the
upstream peak velocity at the centerline ($r=0$) and $\nu$ is the kinematic
viscosity of water ($\nu\mathrm{=1.0\times 10^{-6}m^{2}/s}$ ). Four
experiments were performed with Reynolds values $Re$=[1187, 1427, 1625, 1999]
and $T$=0.96s for all cases. Axial velocity at $z=0$ and $r=0$ as function of
time (fig. 1(c)) was used as time reference for all experiments. The
decelerating phase was defined as $0<t<0.5T$, whereas the accelerating phase
was defined as $0.5T<t<T$.

Figure 1:  (a) Experimental setup. FDS: Flow development section PP:
Programmable pump; R: Reservoir; LS: Laser sheet; CT: constricted tube; C:
camera. (b) Cross-sectional view of a tube of diameter $D=2.6\pm 0.1$ cm and
an annular constriction of diameter $d_{0}=1.6\pm 0.1$cm; (c) $V_{z}$ at the
inlet and $r=0$ for $Re$ =1187. Red circles indicates selected time locations.

###  2.2 Lagrangian Coherent Structures detection

LCS are defined as manifolds that act as flow separatrices, i.e. material
barriers haller_2000 ; haller_2000_chaos , and they correspond to the ridges
of the FTLE fields shadden_2005 ; shadden-exp . The FTLE fields measure the
exponential separation between initially adjacent particles during a finite
time interval ($t_{0},t_{0}+\tau$) haller_2000 ; haller_2002 . This means
that, if $\left\|\delta\mathbf{x}_{0}\right\|$ is the initial separation
between two fluid particles, the maximum separation will be given by

$\left\|\delta\mathbf{x}\right\|_{\textrm{max}}=e^{\left\|\tau\right\|\Lambda}\left\|\delta\mathbf{x}_{0}\right\|$
(1)

where $\Lambda$ is the FTLE field calculated in the $(t_{0},t_{0}+\tau)$ time
interval. To compute FTLE fields, particles trajectories $\mathbf{x}(t)$ were
computed by numerically solving
$\mathbf{\dot{x}}\left(t\right)=\mathbf{v}\left(\mathbf{x},t\right)$.

To this end, the numerical procedure sub-divides the grid of the two-
dimensional velocity field and seeds $3\times 3$ artificial particles in each
of its cells. For simplicity, artificial particles will be named particles
throughout this work. Particles were then advected following a 4th order
Runge-Kutta method and a cubic interpolation. Further refinement of the grid
produced negligible changes in the results. The FTLE fields are then
calculated based on the Cauchy-Green tensor ($CG$), given that
$\left\|\delta\mathbf{x}\right\|_{\textrm{max}}$ is aligned with the
eigenvector with maximum eigenvalue, $\lambda_{\textrm{max}}(CG)$, of the $CG$
tensor. This means that
$\left\|\delta\mathbf{x}\right\|_{\textrm{max}}=\sqrt{\lambda_{\textrm{max}}(CG)}\left\|\delta\mathbf{x}_{0}\right\|$.
The FTLE fields is calculated as

$\Lambda=\frac{1}{\left\|\tau\right\|}\mathrm{ln}(\sqrt{\lambda_{\textrm{max}}(CG)})$
(2)

Trajectories can be integrated forward or backward in time. If $\tau>0$, FTLE
fields measure separation forward in time whereas if $\tau<0$ FTLE fields
measure separation backward in time. These FTLE fields were named
$\Lambda^{+}$ and $\Lambda^{-}$ and they highlight the repelling and
attracting manifolds respectively haller2015 ; haller_2002 .

The integration time $\tau$ was determined based on experimental data such as
the time it takes for a particle located in the vortex boundary to go through
half of that boundary contour, yielding a value of $\tau$= 0.25$T$. To
identify the ridges associated to LCS, all FTLE values below a threshold of
50% of the maximum FTLE value were discarded.

###  2.3 Concentration, Residence time and final position maps

For this study, the Péclet number $Pe=v_{u}D/d$, which compares advection and
diffusion, is of the order of $Pe>$1$\times$106, where $d$ is the mass
diffusivity (10-9m2/s for distilled water). Consequently, the mixing and
transport are primarily due to advection. Therefore, we calculated
concentration maps $C$, residence time maps RT and final position maps of
particles $r_{f}$.

Instantaneous concentration maps $C$ were calculated from particle
trajectories computed in section 2.2. To this end, we get the exact particles
coordinates in each timestep and assign its location to its corresponding DPIV
cell. Finally, we count the number of particles in each DPIV cell cagney_2016
.

####  2.3.1 Residence Time maps and probability of a particle leaving the ROI

For the RT maps, particles were initially distributed uniformly throughout the
ROI with the same spatial resolution used to compute the FTLE fields. Then,
particles were advected from an initial time $t_{0}$ during a timespan of two
pulsatile cycles. Each pixel of the RT map corresponds to the fluid particle
initial position and the pixel value (i.e. color scale) is given by the time
spent by this particle within the ROI. With this computation already used in
badass_2017 , the RT map will depend on the initial time $t_{0}$. To make this
dependence explicit, RT maps will be denoted RT($t_{0}$).

We computed $\mathrm{\overline{RT}}$ by averaging over $t_{0}$ in the time
interval $(0,T)$. The $\mathrm{\overline{RT}}$ map gives the mean residence
time for any point of the ROI regardless of the initial time $t_{0}$. Finally,
the proportion of particles which leave the ROI before a pulsatile cycle,
$P\left(\mathrm{\overline{RT}}<T\right)$, is computed by counting the
particles with $\mathrm{\overline{RT}}$ smaller than a pulsatile period
goulliart_2011 ; Danckwerts1952 .

####  2.3.2 Final position maps and particle exchange index between regions

The first step to evaluate the mixing inside the ROI was to compute the final
position maps ($r_{f}$) where each pixel corresponds to the fluid particle
initial position and the color value will be given by its final position in
$r$. Since particles left the ROI through the boundary at $z/D=1.5$, we
considered its final positions only in $r$ coordinate. In this work we divided
the whole ROI in two subregions $R_{1}$ and $R_{2}$ in order to measure the
particle exchange between these subregions.

Thus, we evaluate the particles which leave through a specific subregion. For
instance, considering a time interval $(t_{0},t_{0}+t_{f})$, the total number
of particles leaving the ROI through $R_{1}$ is $N_{1}^{t_{f}}$ and can be
written as $N_{1}^{t_{f}}=N_{1}^{t_{o}\rightarrow t_{f}}+N_{2\rightarrow
1}^{t_{o}\rightarrow t_{f}}$. Where $N_{1}^{t_{o}\rightarrow t_{f}}$ is the
number of particles in $R_{1}$ at $t=t_{0}$ which left the ROI through $R_{1}$
and $N_{2\rightarrow 1}^{t_{o}\rightarrow t_{f}}$ is the number of particles
in $R_{2}$ at $t=t_{0}$ which left the ROI through $R_{1}$. Normalization by
$N_{1}^{t_{f}}$ yields $1=\mathrm{EI}_{1}+\frac{N_{2\rightarrow
1}^{t_{o}\rightarrow t_{f}}}{N_{1}^{t_{f}}}$, where

$\mathrm{EI}_{1}=\frac{N_{1}^{t_{o}\rightarrow t_{f}}}{N_{1}^{t_{f}}}$ (3)

is defined as an index which measures the proportion of particles initially in
$R_{1}$ to all particles leaving the ROI through $R_{1}$. Analogously, we get
$\mathrm{EI}_{2}=\frac{N_{2}^{t_{o}\rightarrow t_{f}}}{N_{2}^{t_{f}}}$ for
subregion $R_{2}$.

If $\mathrm{EI}_{1}<1$, it means that $\mathrm{R}_{1}$ received particles from
$\mathrm{R}_{2}$ while $\mathrm{EI}_{2}<1$ means that $R_{2}$ received
particles from $R_{1}$. Finally, in order to get a mean rate of exchange for
any time of the pulsatile cycle, we defined $\overline{\mathrm{EI}}_{j}$ as
the mean of $\mathrm{EI}_{j}$ over all the initial conditions $t_{0}$ for
$j=1,2$.

##  3 Results

Figure 2 shows streamlines for $Re=$1187 at $t=0.25T$, $t=0.5T$, $t=0.6T$ and
$t=0.9T$. This flow is axisymmetric, as reported in previous studies
Barrere2020 ; sherwin_2005 ; griffith ; usmani . We observed two different
structures: a high-velocity jet in the central region, and a recirculation
region between the jet and the walls. This recirculation region is where the
vortices form and shed. During the deceleration phase, the recirculation near
the constriction generates a new vortex griffith ; usmani , while the vortex
generated the previous period continues travelling until it finally
dissipates. This flow structure was consistent throughout all the experiments.
A complete description of the flow structure for these experiments was
previously reported in Barrere2020 .

According to previous studies sadan ; casanova ; varghese_2007 ; trip ; Xu2017
, there is no transition to turbulence for the parameters used in our work. In
such a comparison, Reynolds and Womersley
$\left(\alpha=\frac{D}{2}\sqrt{\frac{2\pi}{T\nu}}\right)$ dimensionless
numbers must be taken into consideration. For instance, the works of Cassanova
& Giddens casanova  and later confirmation made by Varghese varghese_2007  do
not find transition to turbulence for an axisymmetric constriction of 75%,
peak $Re$=2240, and $\alpha$=15.6 in the immediate downstream region. Although
$\alpha$=33.6 in our work, other studies Xu2017 ; trip ; sadan  indicate that
for $\alpha>12$, transition to turbulence is independent of $\alpha$.

Figure 2: Streamlines for $Re$ = 1187 at $t$=[0.25 0.5 0.6 0.9]$T$. Black
dashed line indicates the vortex centre.

###  3.1 Lagrangian Coherent Structures

Figure 3 shows vorticity field, $\Lambda^{+}$ and $\Lambda^{-}$ fields, for
$Re$=1187. Similar results were obtained for the other experiments (shown in
appendix A). Figure 3 (a) shows four snapshots of the vorticity field
$\omega(s^{-1})$ for 0.25$T$, 0.5$T$, 0.75$T$ and $T$. The vortex is already
formed at $t=$0.25$T$ and its center is located approximately at $z$=0.5$D$.
At $t$=0.75$T$ a second vortex starts to form near the constriction, while at
$t=T$ the first vortex is close to the edge of the ROI.

Figure 3: For $Re$=1187, at $t$=[0.25 0.5 0.75 1]$T$ (a) vorticity field
$\omega\,(s^{-1})$. Positive vorticity region(counterclockwise) against the
upper wall and negative vorticity region (clockwise) against the lower wall;
(b) positive FTLE fields $\Lambda^{+}\,(s^{-1})$; (c) negative FTLE fields
$\Lambda^{-}\,(s^{-1})$. Both FTLE fields were computed with $\tau$=0.25$T$.
Colorbars on the head of the figure.

Figures 3 (b) and (c) show the ridges of $\Lambda^{+}$ field and $\Lambda^{-}$
field, which represent the repelling and attracting LCS, respectively. In fig.
3 (b) a ridge of the $\Lambda^{+}$ field delimits the trailing edge of the
vortex whereas in fig. 3 (c) a ridge of the $\Lambda^{-}$ field delimits the
leading edge of the vortex. This set of ridges move together with the vortex.
Moreover, figs. 3 (b) and (c) show ridges parallel to the wall that separates
the vortex from the central jet. Analogue LCS representing the leading edge
and the trailing edge of the vortex were found for the other experiments (see
appendix A).

###  3.2 Concentration maps

Figure 4: Instantaneous normalized particle concentration maps $C$ for
$t=$0.25$T$, 0.5$T$, 0.75$T$ and $T$ for $Re$=1187.

Figure 4 shows the concentration maps of particles $C$ normalized by the
number of particles at each time $t$, for $Re=$1187\. During the pulsatile
cycle, the concentration in the leading edge of the vortex becomes larger. By
comparing figures 4 and 3(c), large concentration values correspond to the
$\Lambda^{-}$ ridge of the FTLE fields.

###  3.3 Residence time maps

Figure 5: Residence time maps, RT($t_{0}$=0), for (a) $Re$=1187, (b)
$Re$=1427, (c) $Re$=1625 and (d) $Re$=1999. Each pixel represents the fluid
particle initial position and the color value represents the time spent by the
particle within the ROI normalized by $T$. In (a) particles trajectories
corresponding to particles initial positions in the center of the vortex and
in the edge of the vortex are shown near the lower and upper wall
respectively. (e) Proportion of particles which leave the ROI before a cycle
is finished $P(\mathrm{\overline{RT}}<T)$ as function of $Re$.

Figure 5 shows the RT maps for all the experiments, where $t_{0}=0$ is the
initial time. The color scale of the RT map corresponds to the residence time
normalized by the pulsatile period. Figure 5 shows that the central jet
presents lower residence times compared to the recirculation region for all
experiments carried out in this work. For example, for $Re$=1999 (fig. 5(d)),
RT$<$0.3$T$ in the central jet while in the recirculation region RT$>$0.7$T$.
Moreover, figure 5(a) shows the trajectory of a particle located in the center
of the vortex (lower wall) and in the vortex periphery (upper wall),
evidencing two different dynamics within the vortex. Finally, fig. 5(e) shows
that the proportion of particles that leave the ROI before a period $T$
increases with $Re$. For $Re$ = 1187, 60% of particles left the region, while
95% of particles left for Re = 1999.

###  3.4 Final position maps

Figure 6: Final position maps, $r_{f}(t_{0}=0)$, for (a) $Re$=1187, (b)
$Re$=1427, (c) $Re$=1625 and (d) $Re$=1999. Each pixel represents the fluid
particle initial position and the color value represents its final position in
$r$. Subregions $R1$ and $R2$ are delimited by a dotted line that runs through
the center of the vortex. In (e), flow exchange between subregions is measured
by $\overline{\mathrm{EI}}_{1}$ (solid line) and $\overline{\mathrm{EI}}_{2}$
(dashed line) as a function of $Re$.

Figure 6 shows the $r_{f}$ maps for the initial time $t_{0}=$0\. Based on the
$r_{f}$ maps for all the initial conditions, we found that there is no flow
exchange between the two halves of the ROI delimited by $r=$0\. Therefore, the
number of particles in each region does not change substantially. The vortex
is the main structure transporting particles in the radial direction.
Consequently, we chose the two subregions $R_{1}$ and $R_{2}$ delimited by a
line that runs through the center of the vortex and parallel to the wall of
the tube, see Figs. 6 (a)-(d).

Figures 6 (a)-(d), show that most of the flow initially in either the jet or
the recirculation region, remains in the same region. However, the mixing
occurs within the vortex. For instance, particles initially in the centre of
the vortex have final positions close to the jet i.e. subregion $R_{2}$. While
particles initially close to the periphery of the vortex present layers of
different $r_{f}$ values, meaning that its final position is close to the
walls.

Figure 6(e) shows that as $Re$ increases, $\overline{\mathrm{EI}}_{1}$
decreases while $\overline{\mathrm{EI}}_{2}$ increases. Therefore, $R_{1}$
receives more particles from $R_{2}$, while $R_{2}$ receives less particles
from $R_{1}$. This constitutes R1 as the main region responsible for
increasing the mixing. These parameters can be used as an indicator to
quantify the mixing and complete the description made with the RT maps,
$r_{f}$ maps, concentration maps and $\Lambda^{-}$ fields.

##  4 Discussion

This work studies the mixing of a pulsatile flow in a tube with axisymmetric
constriction. To achieve this, we used Lagrangian-based methods to describe
mixing processes, the final state of the mixing, and also quantify the mixing.
Specifically, we calculated the FTLE fields, concentration maps $C$, and RT
maps, and in order to include a spatial description of the particle exchange,
we introduced the final position maps $r_{f}$.

We found that the ridges of $\Lambda^{+}$ and $\Lambda^{-}$ fields organize
the flow, as they delimit the vortices and divide the ROI into two distinct
regions shadden_2005 ; vetel ; espa ; cagney_2016 . In particular,
$\Lambda^{+}$ ridges, delimit the trailing edge while $\Lambda^{-}$ ridge
delimits the leading edge of the vortex. The flow entering the ROI does not
cross the $\Lambda^{+}$ ridge, separating the flow between successive periods.
Similar results were obtained by Shadden et al.shadden-exp  and Shadden and
Taylor Shadden2008  where LCS predicts the vortex boundaries which organize
the flow that enters and leaves the vortex.

By comparing figures 3(c) and 4 , we established the connection between the
$\Lambda^{-}$ ridges and the $C$ maps. High concentration values correspond
with the $\Lambda^{-}$ ridge because it is the attractive LCS. At $t=$0.25$T$
of fig. 4, particles at $z<0.7D$ correspond to particles upstream the leading
edge of the vortex (see $\Lambda^{-}$ ridge in figure 3(c)). As the flow
evolves, these particles move towards this edge. However, particles initially
located at $z>$0.7$D$ show different behaviors: particles close to the vortex
are attracted by its leading edge, while the others leave the ROI. Finally,
the large difference in concentration near the front suggests that an intense
mixing occurs locally cagney_2016 ; prat .

FTLE fields and $C$ maps give an insight into how the mixing occurs, however
they do not give the final state of the mixing or provide quantitative
information on mixing. To overcome this issue, we complete this information
with RT and $r_{f}$ maps.

A comparison of the RT maps (fig. 5) with the FTLE fields (fig. 3) shows that
regions with different dynamics also have different residence times. For
instance, the residence times in the centerline are approximately RT=[0.8,
0.8, 0.5, 0.3]$T$ for $Re$=[1187, 1427, 1625, 1999], respectively. We observed
that the mean RT in the recirculation region is approximately 2 times larger
than the RT in the centerline. Moreover, we observed that this ratio does not
change with $Re$. This leads to the accumulation of low velocity flow in the
recirculation region and close to the constriction. Jeronimo et al.
jeronimo2019 ; jeronimo2020  observed analogous qualitative behavior, although
RT values are not comparable due to differences between the $Re$ values and
the geometries used.

RT maps also highlighted regions with different dynamics within the vortex.
Figure 5(a) shows that in the center of the vortex, approximately at
$z$=0.25$D$ and $\mid r\mid=$0.4$D$, residence times are lower than in the
rest of the vortex. In order to study this difference in detail, we show the
trajectories of representative particles initially at the center of the vortex
(see fig. 5(a) against the lower wall) and from the edge of the vortex (see
fig. 5(a) against the upper wall). We observed that the particles initially in
the center of the vortex move together with the vortex and leave the ROI
before the particles initially in the edge, which move towards the walls.

In blood flows, the information given by RT maps can be related to flow
stagnation. Previous works zarins2017pathophysiology  show that flow
stagnation is given by high residence time. This can lead to the aggregation
of platelets and blood cells, and a low transport of inhibitors, leading to
the formation of thrombosis.

Previous works describe the final position of particles by lagrangian
techniques dilabbio_2018 ; hendabadi ; jeronimo2019 . Usually final position
maps are represented by color-coded subregions of the ROI. However, in this
work we used $r_{f}$ maps which show the exact particle final position in $r$
coordinate. These maps are more appropriate for our purposes since the $r$
coordinate is normal to the boundary between the main structures of the flow:
vortex and jet. Then, $r_{f}$ maps are a reliable representation of the final
state of the mixing between these structures.

These maps combined with RT maps, give the spatial and temporal final state of
the flow. For instance, we analyze spatio-temporal final state for particles
initially close to $z$=0.25D and $\mid r\mid=$0.4D, by comparing figs. 5 and
6. These results show that those particles which remain less time in the ROI,
have final positions either within the center of the vortex or the central jet
(approximately -0.2D$<r<$0.2D). On the other hand, the particles that were
initially at the edge of the vortex show the highest residence time (see fig.
5) and its final positions are close to the walls (see fig. 6). We conclude
that for all the values of $Re$, the $r_{f}$ and RT maps delimit the same
regions within the ROI. The regions with higher residence time have a final
position close to the walls, whereas the regions with lower residence time
have $r_{f}$ in the central jet.

By calculating the $P(\mathrm{\overline{RT}}<T)$ and $\mathrm{\overline{EI}}$
parameters, we were able to quantify the mixing as $Re$ increases. Figure 5(e)
showed that residence times decrease as $Re$ increases, which means that
mixing is increased with $Re$.

Although this is an expected result, the incidence of each structure in mixing
is not self-evident, considering that the flow is not turbulent. However,
further analysis of $\mathrm{\overline{EI}}$ gives an insight on the incidence
of the vortex in the mixing. For instance, fig. 6(e) shows that
$\mathrm{\overline{EI}}_{1}$ decreases and $\mathrm{\overline{EI}}_{2}$
increases with $Re$. This means that $R_{1}$ receives more particles from
$R_{2}$ as $Re$ increases. Since $R_{1}$ is within the vortex region, these
results shows that the vortex is the main structure responsible for the mixing
between regions.

Based on previous works gharib ; Barrere2020 , we suggest this mixing
behaviour is due to circulation implications on vortex shedding. In a flow
being pushed through a constriction, vortex shedding occurs after reaching a
specific value of circulation, while any excess of circulation goes to a trail
of vorticity. This mechanism explains that as $Re$ is increased, vortex area
is enlarged. Therefore, the increment of $Re$ increases the proportion of
particles that attach to the vortex finally leading to an increment of
particles in $R_{1}$. Noteworthy, the introduction of
$\mathrm{\overline{EI}}_{j}$ parameter provides a direct and quantitative
measurement of the effects caused by the aforementioned mechanism.

Overall, we used LCS and Concentration maps to explain the behavior of the
main flow structures -i.e., Vortex and central jet – and the mechanism
underlying mixing. However, these Lagrangian-based methods do not provide
quantitative information on mixing. Since mixing is an emergent spatiotemporal
phenomenon, it is difficult to quantify and depends on each application. We
consider that in this work, a step was made in giving a comprehensive
description and quantification of mixing in our specific system. Specifically,
the introduction of $r_{f}$ maps and its combination with RT maps and
$P(\mathrm{\overline{RT}}<T)$ and $\mathrm{\overline{EI}}_{j}$ parameters
provide a set of descriptors to study and measure mixing as a function of
Reynolds number.

Future studies should focus on applying these methods in more specific
hematological situations. For instance, a possible application of the
$P(\mathrm{\overline{RT}}<T)$ parameter is to measure the time of a process of
interest. Previous works specify the relevance of thrombi transport times
Rayz_2010 , or the residence time of a flow in a region in pathological
conditions dilabbio_2018 ; arzani_stagnation . Thus, RT maps, and particularly
$P(\mathrm{\overline{RT}}<T)$, can be used to measure the time of such
processes. Similarly, the combination of RT maps and $r_{f}$ can be useful in
applications involving, for example, the transport of blood particles or drug
delivery. These problems require identifying specific particles, finding
stagnation points, computing particle concentration in different regions, and
exchanging particles between specific regions.

##  5 Conclusions

In this work, we used FTLE fields, concentration maps, RT residence time maps,
and $r_{f}$ final position maps as tools to describe the mixing of a pulsatile
flow in a tube with an axisymmetric constriction. The evolution of how mixing
occurs was given by the evolution of FTLE fields and $C$ maps. The FTLE fields
identified a set of LCS that separates the vortex from the central jet and
separates the flow between successive periods and sates the vortex leading
edge as the structure where particles agglomerate the most.

We get an insight into the final state of the mixing through RT and $r_{f}$
maps. These maps are scalar fields computed from particle trajectories and
capable of giving information about the topology and mixing of the flow. Based
on them, we calculated the $P(\mathrm{\overline{RT}}<T)$ and
$\mathrm{\overline{EI}}_{j}$ parameters to quantify the mixing with $Re$.
These parameters confirmed that the mixing increases as $Re$ increases and
provided an interpretation of how this occurs. For instance, an analysis of
$\mathrm{\overline{EI}}_{j}$, RT, and $r_{f}$ maps showed that the vortex is
the main responsible in mixing process and gives the basis to suggest a
mechanism behind this behavior. Finally, we conclude that by cross-referencing
the information obtained by each of these tools, we accomplished a better
understanding of the mixing processes.

##  6 Acknowledgements

N.B. acknowledges the funds granted by ANII (Agencia Nacional de Investigación
e Innovación, Uruguay) as part of the Doctoral scholarship
POSNAC-2015-1-109843. This research was supported by CSIC (Comisión Sectorial
de Investigación Científica (CSIC), Uruguay), through the 2016 I+D project
“Estudio dinámico de un flujo pulsátil y sus implicaciones hemodinámicas
vasculares” and CSIC’s group grant “CSIC2018 - FID13 - grupo ID 722” and by
PEDECIBA, Uruguay.

##  Appendix A Additional results

Figures 7, 8 and 9 show the vorticity fields, $\Lambda^{+}$ and $\Lambda^{-}$
for $Re$=1427, $Re$=1625 and $Re$=1999, respectively. In all cases, we
observed structures analogous to those in fig. 3, with $\Lambda^{+}$ and
$\Lambda^{-}$ ridges which delimit the boundaries of the vortex. As $Re$
increases, the flow becomes disorganized, particularly in the deceleration
phase, and the propagation velocity of the flow increasesBarrere2020 . For
$Re=$1999 (fig.9), the phases of the vortex are represented by $t=$0.25$T$,
$t=$0.53$T$, $t=$0.85$T$, and $t=T$. At $t$=0.53$T$ we observed disorganized
patterns in the vorticity field and the $\Lambda^{-}$ field due to the vortex
tail that had previously left the ROI. Compared with lower $Re$, the ridges of
the $\Lambda^{-}$ field were more filamentous. For $t$=0.85$T$, on the other
hand, a new vortex that had not yet detached started to form, hence only the
boundary delimited by the ridges of the $\Lambda^{+}$ field appeared.

Figure 7: For $Re$=1427, at $t$=[0.25 0.5 0.75 1]$T$ (a) vorticity field
$\omega\,(s^{-1})$; (b) positive FTLE fields $\Lambda^{+}\,(s^{-1})$; (c)
negative FTLE fields $\Lambda^{-}\,(s^{-1})$. Both FTLE fields were computed
with $\tau$=0.25$T$.

Figure 8: For $Re$=1625, at $t$=[0.25 0.5 0.75 1]$T$ (a) vorticity field
$\omega\,(s^{-1})$; (b) positive FTLE fields $\Lambda^{+}\,(s^{-1})$; (c)
negative FTLE fields $\Lambda^{-}\,(s^{-1})$. Both FTLE fields were computed
with $\tau$=0.25$T$.

Figure 9: For $Re$=1999, at $t$=[0.25 0.53 0.85 1]$T$ (a) vorticity field
$\omega\,(s^{-1})$; (b) positive FTLE fields $\Lambda^{+}\,(s^{-1})$; (c)
negative FTLE fields $\Lambda^{-}\,(s^{-1})$. Both FTLE fields were computed
with $\tau$=0.25$T$.

## References

  *   * (1) Vlachopoulos, C., O’Rourke, M., Nichols, W.W.: McDonald’s Blood Flow in Arteries, Sixth Edition: Theoretical, Experimental and Clinical Principles. HODDER & STOUGHTON, ??? (2011). https://www.ebook.de/de/product/11728652/charalambos_vlachopoulos_michael_o_rourke_wilmer_w_nichols_mcdonald_s_blood_flow_in_arteries_sixth_edition_theoretical_experimental_and_clinical_principles.html
  * (2) Naghavi, M., Libby, P., Falk, E., Casscells, S.W., Litovsky, S., Rumberger, J., Badimon, J.J., Stefanadis, C., Moreno, P., Pasterkamp, G., Fayad, Z., Stone, P.H., Waxman, S., Raggi, P., Madjid, M., Zarrabi, A., Burke, A., Yuan, C., Fitzgerald, P.J., Siscovick, D.S., de Korte, C.L., Aikawa, M., Airaksinen, K.E.J., Assmann, G., Becker, C.R., Chesebro, J.H., A.Farb, Galis, Z.S., Jackson, C., Jang, I., Koenig, W., Lodder, R.A., March, K., Demirovic, J., Navab, M., Priori, S.G., Rekhter, M.D., Bahr, R., Grundy, S.M., Mehran, R., Colombo, A., Boerwinkle, E., Ballantyne, C., Insull, W., Schwartz, R.S., Vogel, R., Serruys, P.W., Hansson, G.K., Faxon, D.P., Kaul, S., Drexler, H., Greenland, P., Muller, J.E., Virmani, R., Ridker, P.M., Zipes, D.P., Shah, P.K., Willerson, J.T.: From vulnerable plaque to vulnerable patient: a call for new definitions and risk assessment strategies: Part i. Circulation 108(14), 1664–1672 (2003). https://doi.org/10.1161/01.cir.0000087480.94275.97
  * (3) Organization, W.H.: Global Atlas on Cardiovascular Disease Prevention and Control. WORLD HEALTH ORGN, ??? (2012). https://www.ebook.de/de/product/18306761/world_health_organization_global_atlas_on_cardiovascular_disease_prevention_and_control.html
  * (4) Liu, X., Fan, Y., Deng, X., Zhan, F.: Effect of non-newtonian and pulsatile blood flow on mass transport in the human aorta. Journal of Biomechanics 44(6), 1123–1131 (2011). https://doi.org/10.1016/j.jbiomech.2011.01.024
  * (5) Shadden, S.C., A.Arzani: Lagrangian postprocessing of computational hemodynamics. Ann Biomed Eng 43, 41–58 (2015)
  * (6) Winter, D.C., Nerem, R.M.: Turbulence in pulsatile flows. Annals of Biomedical Engineering 12(4), 357–369 (1984). https://doi.org/10.1007/bf02407780
  * (7) Peacock, J., Jones, T., Lutz, R.: The onset of turbulence in physiological pulsatile flow in a straight tube. Experiments in Fluids 24, 1–9 (1998)
  * (8) Trip, R., Kuik, D.J., Westerweel, J., Poelma, C.: An experimental study of transitional pulsatile pipe flow. Physics of fluids 24 (2012)
  * (9) Brindise, M.C., Vlachos, P.P.: Pulsatile pipe flow transition: Flow waveform effects. Physics of Fluids 30(1), 015111 (2018). https://doi.org/%****␣sn-article.bbl␣Line␣225␣****10.1063/1.5021472
  * (10) Xu, D., Varshney, A., Ma, X., Song, B., Riedl, M., Avila, M., Hof, B.: Nonlinear hydrodynamic instability and turbulence in pulsatile flow. Proceedings of the National Academy of Sciences 117(21), 11233–11239 (2020). https://doi.org/10.1073/pnas.1913716117
  * (11) Jain, K.: Transition to turbulence in an oscillatory flow through stenosis. Biomechanics and Modeling in Mechanobiology 19(1), 113–131 (2019). https://doi.org/10.1007/s10237-019-01199-1
  * (12) Jain, K.: The effect of varying degrees of stenosis on transition to turbulence in oscillatory flows. Biomechanics and Modeling in Mechanobiology 21(3), 1029–1041 (2022). https://doi.org/10.1007/s10237-022-01579-0
  * (13) Haller, G., Yuan, G.: Lagrangian coherent structures and mixing in two-dimendional turbulence. Physica D 147, 352–370 (2000)
  * (14) Haller, G.: Finding finite-time invariant manifolds in two-dimendional velocity ffields. Chaos 10, 1–10 (2000)
  * (15) Shadden, S.C.: A dynamical systems approach to unsteady systems. PhD thesis, California Institute of Technology (2006)
  * (16) Shadden, S.C., Leiken, F., Marsden, J.E.: Definition and properties of lagrangian coherent structures from finite-time lyapunov exponents in two-dimendional aperiodic flows. Physica D 212, 271–304 (2005)
  * (17) Arzani, A., Shadden, S.C.: Characterization of the transport topology in patient-specific abdonimal aortic aneurysm models. Physics of Fluids 24 (2012)
  * (18) Badas, M.G., Domenichini, F., Querzoli, G.: Quantification of the blood mixing in the left ventricle using finite time lyapunov exponents. Meccanica 52(3), 529–544 (2017)
  * (19) Vétel, J., Garon, A., Pelletier, D.: Lagrangian coherent structures in the human body carotid artery bifurcation. Exp. Fluids 46, 1067–1079 (2009)
  * (20) Espa, S., et al.: A lagrangian investigation of the flow inside the left ventricle. European Journal of Mechanics B 35, 9–19 (2012)
  * (21) Camassa, R., Wiggins, S.: Chaotic advection in a rayleigh-benard flow. Physical Review A 43(2), 774–797 (1991)
  * (22) Gouillart, E., Dauchot, O., Thiffeault, J.L.: Measures of mixing quality in open flows with chaotic advection. Physics of fluids 23, 013604 (2011)
  * (23) Monsen, N.E., Lucas, J.E., Monismith, L.V., Stephen, G.: A comment on the use of flushing time, residence time, and age as transport time scales. Limnology and Oceanography 47(5), 1545–1553 (2002)
  * (24) Rayz, V.L., Boussel, L., Ge, L., Leach, J.R., Martin, A.J., Lawton, M.T., McCulloch, C., Saloner, D.: Flow residence time and regions of intraluminal thrombus deposition in intracranial aneurysms. Annals of Biomedical Engineering 38(10), 3058–3069 (2010). https://doi.org/10.1007/s10439-010-0065-8
  * (25) Seo, J.H., Mittal, R.: Effect of diastolic flow patterns on the function of the left ventricle. Physics of Fluids 25(11), 110801 (2013). https://doi.org/10.1063/1.4819067
  * (26) Long, C.C., Esmaily-Moghadam, M., Marsden, A.L., Bazilevs, Y.: Computation of residence time in the simulation of pulsatile ventricular assist devices. Computational Mechanics 54(4), 911–919 (2013). https://doi.org/%****␣sn-article.bbl␣Line␣475␣****10.1007/s00466-013-0931-y
  * (27) Labbio, G.D., Vétel, J., Kadem, L.: material transport in left ventricle with aortic valve regurgitation. Physical Review fluids 3 (2018)
  * (28) Fuster, V., Fayad, A.Z., Moreno, P.R., Poon, M., Corti, R., Badimon, J.J.: Atherotrombosis and high-risk plaque. part ii: approaches by noninvasive computed tomographic/ magnetic resonance imaging. Journal of the American College of Cardiology 46(7) (2005)
  * (29) Jeronimo, M.D., Zhang, K., Rival, D.E.: Direct lagrangian measurements of particle residence time. Experiments in Fluids 60(4) (2019). https://doi.org/10.1007/s00348-019-2718-1
  * (30) Jeronimo, M.D., Rival, D.E.: Particle residence time in pulsatile post-stenotic flow. Physics of Fluids 32(4), 045110 (2020). https://doi.org/10.1063/1.5144388
  * (31) Barrere, N., Brum, J., L’her, A., Sarasúa, G.L., Cabeza, C.: Vortex dynamics under pulsatile flow in axisymmetric constricted tubes. Papers in Physics 12, 120002 (2020). https://doi.org/10.4279/pip.120002
  * (32) Westerweel, J.: Fundamentals of digital particle velocimetry. Measurement science and technology 8, 1379–1392 (1997)
  * (33) Thielicke, W., Eize J. Stamhuis: PIVlab - Time-Resolved Digital Particle Image Velocimetry Tool for MATLAB. figshare (2019). https://doi.org/%****␣sn-article.bbl␣Line␣575␣****10.6084/M9.FIGSHARE.1092508
  * (34) Shadden, S.C., Dabiri, J.O., Marsden, J.E.: Lagrangian analysis of fluid transport in empirical vortex ring flows. Physics of fluids (2006)
  * (35) Haller, G.: Lagrangian coherent structures from approximate velocity data. Physics of fluids 14, 1851–1861 (2002)
  * (36) Haller, G.: Lagrangian coherent structures. Annual Review of Fluid Mechanics 47(1), 137–162 (2015). https://doi.org/10.1146/annurev-fluid-010313-141322
  * (37) Cagney, N., Balabani, S.: Lagrangian strucures and mixing in the wake of a streamwise oscillating cylinder. Physics of Fluids 28 (2016)
  * (38) Danckwerts, P.V.: The definition and measurement of some characteristics of mixtures. Applied Scientific Research 3(4), 279–296 (1952). https://doi.org/10.1007/bf03184936
  * (39) Sherwin, S.J., Blackburn, H.M.: Three dimensional instabilities and transition of steady and pulsatile axisymmetric stenotic flows. J. Fluid Mech. 533, 297–327 (2005)
  * (40) Griffith, M.D., Leweke, T., Thompson, M.C., Hourigan, K.: Pulsatile flow in stenotic geometries: flow behaviour and stability. J.Fluid Mech. 622, 291–320 (2009)
  * (41) Usmani, A.Y., Muralidhar, K.: Pulsatile flow in a compliant stenosed asymmetric model. Exp. Fluids 57(12), 186 (2016)
  * (42) Stettler, J.C., Hussain, A.K.M.F.: On transition of pulsatile pipe flow. J. Fluid Mech. 170, 169–197 (1986)
  * (43) Cassanova, R.A., Giddens, D.P.: Disorder distal to modeled stenoses in steady and pulsatile flow. J. Biomechanics 11, 441–453 (1977)
  * (44) Varghese, S.S., Frankel, S., Fischer, P.: Direct numerical simulation of stenotic flows. part 2. pulsatile flow. J. Fluid Mech. 582, 281–318 (2007)
  * (45) Xu, D., Warnecke, S., Song, B., Ma, X., Hof, B.: Transition to turbulence in pulsating pipe flow. Journal of Fluid Mechanics 831, 418–432 (2017). https://doi.org/10.1017/jfm.2017.620
  * (46) Shadden, S.C., Taylor, C.A.: Characterization of coherent structures in the cardiovascular system. Annals of Biomedical Engineering 36(7), 1152–1162 (2008). https://doi.org/10.1007/s10439-008-9502-3
  * (47) Pratt, K.R., J.Meiss, Crimaldi, J.: Reaction enhancement of initially distant scalars by lagrangian coherent structures. Physics of fluids 27(3), 035106 (2015)
  * (48) Zarins, C.K., Xu, C.: Pathophysiology of human atherosclerosis. In: Vascular Surgery, pp. 37–60. CRC Press, ??? (2017)
  * (49) Hendabadi, S., et al.: Topology of blood transport in the human left ventricle by novel processing of doppler echocardiography. Annals of Biomedical Engineering 41, 2603–2616 (2013)
  * (50) Gharib, M., Rambod, E., Shariff, K.: A universal time scale for vortex ring formation. Journal of Fluid Mechanics 360, 121–140 (1998)
  * (51) Arzani, A., Gambaruoto, A.M., Chen, G., Shadden, S.C.: wall shear stress exposure time: a lagrangian measure of near-wall stagnation and concentration in cardiovascular flows. Biomechanics and Modeling in Mechanobiology 16, 787–803 (2017)


# Two-dimensional Lorentz-violating Casimir effect

1,2K. E. L. de Farias<EMAIL_ADDRESS>1M. A. Anacleto
<EMAIL_ADDRESS>1,2F. A. Brito<EMAIL_ADDRESS>1E. Passos
<EMAIL_ADDRESS>2 Herondy Mota<EMAIL_ADDRESS>1Departamento de
Física, Universidade Federal de Campina Grande,
Caixa Postal 10071, 58429-900, Campina Grande, Paraíba, Brazil.  2Departamento
de Física, Universidade Federal da Paraíba,
Caixa Postal 5008, João Pessoa, Paraíba, Brazil.

###### Abstract

In this study, we consider the four-dimensional Maxwell electrodynamics
extended with CPT-even Myers-Pospelov Lorentz-violating dimension-six
operators to investigate the associated two-dimensional properties in the
context of quantum vacuum fluctuation effects, namely, the Casimir effect.
Upon projecting out the 4D theory down to a 2D theory we obtain analogs of
these operators leading to a modified dispersion relation in a Lorentz
invariance violation (LIV) scalar model equivalent to the electromagnetic
theory. By making use of the modified dispersion relation, we derive exact
analytic expressions for the Casimir energy and force induced by imposing
Dirichlet boundary conditions on the scalar field. In the regime where the LIV
parameter becomes very small, we recover known results for the Casimir energy
and force plus correction terms due to the LIV.

###### pacs:

11.15.-q, 11.10.Kk

##  I Introduction

The research of Lorentz invariance violation (LIV) in physics has garnered
significant attention due to its potential to explain natural phenomena that
the standard model of particle physics cannot account for, particularly in
high-energy regimes near the Planck energy scale $(\sim 10^{19}\rm GeV)$. The
LIV can be explored in various ways, such as by the notorious extension of the
standard model [1], by the introduction of the aether field [2], by the
Hořava-Lifshitz anisotropic theory [3], and by the incorporation of high
derivative operators [4], which is the starting point of this study. Our focus
is to elaborate an investigation of low-dimensional physics associated with
LIV high derivative operators.

The study of several Lorentz-violating scenarious in the low-dimension field
theory has been considered in Refs. [5, 6, 7, 8, 9, 10, 11, 12]. Among these
works we can highlight as motivations for our investigation the one conducted
in Ref.[5], which considers the Lorentz-violating Chern-Simons model defined
in 4D and performs a dimensional reduction to 3D generating a new model that
mixes the gauge field, massless scalar field and an external vector which is
used to verify the LIV effects in planar systems. Also, in Ref.[6], the
authors have conducted a 2D LIV study attributed to the scalar effective
theory aiming to investigate the effects of LIV in the context of defect
structures (see also Ref.[9], which study a class of LIV travelling solitons
systems). In addition, we highlight the study of Ref.[11], which introduces
the 2D LIV higher derivative term via radiative corrections and four-
dimensional projection down to two dimensions. More specifically, this work
considers the Myers-Pospelov Lorentz-violating dimension-five operators to
investigate 4D-2D projection, and then, apply in the defect scenarios (see
also Ref.[12]). Here, in our investigation, we shall focus on 2D effects from
the 4D theories with high derivative Lorentz-violating operators such as the
4D effective theories with Lorentz-violating dimension-six operators.

We initially consider the gauge sector of the modified theory in order to
derive the physics in two dimensions with analogues of these Lorentz-violating
operators. As the main object of the dimensional reduction process, we shall
consider the electrodynamics higher derivative operators of dimension-six
given as
${\mathcal{L}}_{4D}\sim\beta^{(6)}/M_{\mathrm{QG}}^{2}F^{\mu\nu}u_{\mu}u^{\sigma}(u\cdot\partial)^{2}\tilde{F}_{\sigma\nu}$,
with $\beta^{(6)}$ being a dimensionless parameters that control the LIV
intensity over $M_{\mathrm{QG}}$ scale suppression, where
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, is the field-
strength tensor [13]. The main property of this term is that it modifies the
dispersion relation with a quadratic term to the energy content of the system
[14]. Thus, we take advantage of this procedure to obtain a 2D theory with
scalar fields coupled to dimensional-four operators. Our goal is to
investigate new LIV effects on physical observables associated with 2D
systems, as for instance, the Casimir energy and force.

The detection of a LIV event implies the existence of new physics to explore,
complementing the current standard model of particle physics. However,
experiments aimed at detecting and measuring LIV are inherently challenging.
The increasing sensitivity of experiments designed to detect Casimir forces
offers a novel approach to study this phenomenon. The Casimir effect is a
natural occurrence that arises when two parallel plates are placed close to
each other (typically within micrometers). This arrangement limits the
frequencies of vacuum fluctuations between the plates, resulting in a
difference in force between the interior and exterior of the plates. First
discovered by H. Casimir in 1948 [15], an attempt to detect it was initially
made by Sparnaay in 1958 [16]. However, recent advancements have allowed for
high-precision detection by various experiments [17, 18, 19], highlighting the
evolving sensitivity in detecting this force.

With the evolving sensitivity in detecting Casimir forces, the idea of a
Casimir theory that violates Lorentz symmetry becomes increasingly plausible
and extensively explored in the literature [20, 21, 22]. In such models, the
LIV parameter introduces corrections to the Casimir energy. However, these
proposals often result in expressions without exact analytic solutions,
necessitating the use of approximations to compute the Casimir energy and
force. In our approach, however, we provide exact analytic solutions without
resorting to approximations. The significance of obtaining a closed-form
solution lies in the ability to fully understand the behaviour of the system
and explore its limits, thereby elucidating the physical constraints of the
parameters involved.

The present work is structured as follows: In Section II, we develop the
dimensional projection of a 4D theory into a 2D one for the gauge sector of an
extended higher derivative electrodynamics. The Section III addresses the
derivation of the modified equation of motion and provides insights into the
modified dispersion relation. Moving on to Section IV, we apply Dirichlet
boundary conditions on the free solution of the modified equation of motion
and use the corresponding modified dispersion relation to obtain exact
analytic expressions for the 2D Casimir energy, including its approximation
for the regime of small values for the LIV parameter, demonstrating the
recovery of known results found in literature. Section V focuses on deriving
the analytic expression for the Casimir force, along with its approximation
for small values of the LIV parameter, again recovering the standard result.
Finally, in Section VI, we offer our conclusions based on the results
obtained. Throughout this paper, we utilize natural units with $\hslash=c=1$.

##  II 4D-2D DIMENSIONAL PROJECTION

In this initial point, we focus on the dimensional projection out of the four-
dimensional gauge field components of a $4\mathrm{D}$ theory into a 2D theory.
And, as previously mentioned, we consider the Maxwell electrodynamics plus a
CPT-even higher derivative six-operator [14]. The Lagrangian density for this
configuration of the system is then given by

$\displaystyle\mathcal{L}_{4D}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}-\frac{\bar{g}}{2}(u\cdot\partial)F^{\mu\nu}u_{\mu}u^{\sigma}(u\cdot\partial)\tilde{F}_{\sigma\nu},$
(1)

where $\bar{g}=\beta^{(6)}/M_{\mathrm{QG}}^{2}$ is the LIV parameter with
dimension of length squared. Notice that the other quantities in Eq.(1) have
already been introduced in Sec.I.

Our simple 4D-2D dimensional projection will follow the approach of Ref.[11],
which implies to say that $\partial_{1}=\partial_{2}=0$ when acting on any
field (therefore, $\left.F_{ab}=0\right)$ being $a,b=1,2$. In this way, the
Maxwell term assumes the form

$\displaystyle-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$
$\displaystyle=-\frac{1}{4}\left(F_{ij}F^{ij}+2\partial_{0}A_{a}\partial^{0}A^{a}+2\partial_{3}A_{a}\partial^{3}A^{a}\right)$
$\displaystyle=-\frac{1}{4}F_{ij}F^{ij}+\frac{1}{2}\partial_{i}\phi\partial^{i}\phi+\frac{1}{2}\partial_{i}\chi\partial^{i}\chi,$
(2)

where $i,j=0,3$. We will work with the metric signature $(+,-,-,-)$, where
$A_{1}=\phi$ and $A_{2}=\chi$. Let us now discuss the CPT-even extension.
Under the same conditions described above we obtain

$\displaystyle-\frac{\bar{g}}{2}(u\cdot\partial)F^{\mu\nu}u_{\mu}u^{\sigma}(u\cdot\partial)\tilde{F}_{\sigma\nu}=$
$\displaystyle-\frac{\bar{g}}{2}\left(F_{ij}u^{i}u_{k}\left(u_{l}\partial^{l}\right)^{2}F^{kj}+\left(\partial_{i}A_{a}\right)u^{i}u_{j}\left(u_{l}\partial^{l}\right)^{2}\left(\partial^{j}A^{a}\right)\right)$
$\displaystyle=$
$\displaystyle-\frac{\bar{g}}{2}F_{ij}u^{i}u_{k}\left(u_{l}\partial^{l}\right)^{2}F^{kj}+\frac{\bar{g}}{2}\left(\partial_{i}\phi\right)u^{i}u_{j}\left(u_{l}\partial^{l}\right)^{2}\left(\partial^{j}\phi\right)+$
$\displaystyle\frac{\bar{g}}{2}\left(\partial_{i}\chi\right)u^{i}u_{j}\left(u_{l}\partial^{l}\right)^{2}\left(\partial^{j}\chi\right).$
(3)

Therefore, from a simple method of 4D-2D dimensional projection, we can
construct other LIV effective models in two-dimensions. Among them, a LIV
scalar model equivalent to the electromagnetic theory given by Eq.(1) can be
similar to

$\mathcal{L}_{2D}=\frac{1}{2}\partial_{i}\phi\partial^{i}\phi+\frac{\bar{g}}{2}\left(\partial_{i}\phi\right)u^{i}u_{j}\left(u_{l}\partial^{l}\right)^{2}\left(\partial^{j}\phi\right)+\cdot\cdot\cdot,$
(4)

where $(\cdot\cdot\cdot)$ means terms associate to $\chi-$field. This
resulting model is made out of two-dimensional fields living in two
dimensions. In two dimensions, the new term corresponds to a four-dimensional
operator Lorentz-violating. We shall refer to it as the 2D CPT-even higher-
derivative scalar model.

##  III A Brief Analysis of the Model

Considering the result given by Eq.(4), we write the following $(1+1)D$
massless Lagrangian density:

$\mathcal{L}_{\phi}=\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi+\frac{\bar{g}}{2}\left(\partial_{\mu}\phi\right)u^{\mu}u_{\nu}\left(u_{\alpha}\partial^{\alpha}\right)^{2}\left(\partial^{\nu}\phi\right),$
(5)

where $\mu,\nu=0,1$ with the constant parameter,
$u_{\mu}=\left(u_{0},u_{1}\right)$, controlling the LIV effect in the theory.
Notice that we are working in $2\mathrm{D}$ such that the metric signature is
given as $(+,-)$. The higher-derivative equation of motion associated with the
Lagrangian, $\mathcal{L}_{\phi}$, can be derived from

$\frac{\partial\mathcal{L}_{\phi}}{\partial\phi}-\partial_{\lambda}\left(\frac{\partial\mathcal{L}_{\phi}}{\partial\left(\partial_{\lambda}\phi\right)}\right)+\partial_{\lambda}\partial_{\rho}\left(\frac{\partial\mathcal{L}_{\phi}}{\partial\left(\partial_{\lambda}\partial_{\rho}\phi\right)}\right)-\partial_{\lambda}\partial_{\rho}\partial_{\sigma}\left(\frac{\partial\mathcal{L}_{\phi}}{\partial\left(\partial_{\lambda}\partial_{\rho}\partial_{\sigma}\phi\right)}\right)+\cdots=0.$
(6)

Consequently, we obtain a modified equation of motion as follows

$\square\phi+\bar{g}u^{\alpha}u^{\beta}u^{\mu}u^{\nu}\left(\partial_{\alpha}\partial_{\beta}\partial_{\mu}\partial_{\nu}\phi\right)=0.$
(7)

This equation offers us the following set of equations:

$\ddot{\phi}-\phi^{\prime\prime}+\bar{g}\ddddot{\phi}=0,\quad{\rm time-
like\;approach}:u_{\mu}=\left(1,0\right),$ (8a)
$\ddot{\phi}-\phi^{\prime\prime}+\bar{g}\phi^{\prime\prime\prime\prime}=0,\quad{\rm
space-like\;approach}:u_{\mu}=\left(0,1\right),$ (8b)

where dots and primes stand for derivatives with respect to time $t$ and
position $x$, respectively.

The solution to Eqs. (8a) and (8b) is of the plane wave type, that is,
$\phi\propto e^{-iEt+ikx}$, with $k$ being the momentum in the $x$-direction.
Consequently, we can derive the modified dispersion relation associated with
this solution. Hence, by considering that, $E\sim k_{x}$ (for high energies),
the equation of motion associated with the time-like case, Eq.(8a), along with
the plane wave solution, provide the following energy-momentum equation:

$\displaystyle E^{2}-k^{2}_{x}+\bar{g}k^{4}_{x}=0.$ (9)

Solving we obtain the energy solution

$\displaystyle E(k_{x})=|k_{x}|\sqrt{1-\bar{g}k_{x}^{2}}.$ (10)

This solution is similar to the dispersion relation associated with 4D theory
(see Ref.[14]) and provides evidence for the LIV. Note also that the imaginary
part of Eq.(9), can lead to the loss of unitarity and instabilities (we plot
the energy solution, Eq.(10), in Fig.1).

Figure 1: The behaviour of the dispersion relation in the time-lime approach,
Eq. (10). The black dashed line represents the light cone and the blue solid
line shows the energy $E$.

On the other hand, the space-like case dispersion relation associated with the
motion of equation, Eq.(8b), is given as

$\displaystyle E^{2}-k^{2}_{x}-\bar{g}k^{4}_{x}=0.$ (11)

This leads us to the following energy solution

$\displaystyle E(k_{x})=|k_{x}|\sqrt{1+\bar{g}k_{x}^{2}}.$ (12)

Differently from the time-like case, this solution has no imaginary
contribution, and then, with no loss of stability for approximation of higher
momenta (see its behaviour in Fig.2).

Figure 2: The behaviour of the dispersion relation in the space-like approach,
Eq. (12). The black dashed line represents the light cone and the blue solid
line shows the energy $E$.

At this point, we present a simple analysis of group and phase velocities.
Considering only the dispersion relation, given by Eq.(12), leads to group and
phase velocities, respectively,

$\displaystyle
v_{g}(k_{x})=\frac{(1+2\bar{g}k_{x}^{2})}{\sqrt{1+\bar{g}k_{x}^{2}}},$ (13a)
$\displaystyle v_{p}(k_{x})=\sqrt{1+\bar{g}k_{x}^{2}}.$ (13b)

Note that the phase and group velocities above are related through Rayleigh’s
formula: $v_{p}/v_{g}=1-\big{(}E/v_{p}\big{)}\big{(}dv_{p}/dk_{x}\big{)}$.
Hence, from Eqs.(13a)–(13b), we find

$\displaystyle\frac{v_{p}-v_{g}}{v_{g}}$ $\displaystyle=$
$\displaystyle-\frac{\bar{g}k_{x}^{2}}{\sqrt{1+\bar{g}k_{x}^{2}}}$ (14)
$\displaystyle\approx$ $\displaystyle-\bar{g}k_{x}^{2}+{\cal
O}(\bar{g}k_{x}^{2}),$

for the condition $\bar{g}|\vec{k}|^{2}\ll 1$. In other words, Eq.(14) reveals
a subluminal propagation (for $v_{p}>v_{g}$, we have a normal dispersion
medium). Therefore, we can conclude that the model is truly isotropic, and
thus, being important for phenomenological analyzes.

##  IV Casimir Energy

As shown earlier, the current model presents two alternatives for expressing
energy: the time-like approach, given by (8a), as well as the space-like
approach, given by (8b). While both can be considered for analyzing the
system’s energy, the timelike approach introduces inconsistencies to the
energy when examining the dispersion relation, as explained in the last
section. Consequently, we exclusively use the spacelike approach to derive the
energy density.

By examining the dispersion relation obtained for the spacelike approach, we
can determine the vacuum energy fluctuation, i.e., the Casimir energy. In
order to do that we aim to consider a configuration constituting of two
parallel plates, similar to the original Casimir system, which leads us to an
energy density depending on the volume between the plates. However, due to the
dimensional reduction, the effect of the plates is in fact codified in two
points along the line in the $x$-direction. Hence, in our effective $(1+1)D$
system, by applying Dirichlet boundary conditions on the free solution of the
equation of motion discussed in the previous section, the dispersion relation
$E^{2}=k_{x}^{2}+\bar{g}k_{x}^{4}$ now with discretized momentum in the
$x$-direction becomes

$\displaystyle E=$
$\displaystyle\frac{1}{2}\sum_{n=0}^{\infty}\sqrt{\left(\frac{n\pi}{a}\right)^{2}+\bar{g}\left(\frac{n\pi}{a}\right)^{4}},$
(15)

where $a$ is the distance between the boundary condition points. It is clear
that the sum in (15) yields a divergent result. However, we can employ the
Abel-Plana formula to obtain the finite part of this sum, and discard the
infinite continuum contribution by means of the minimum subtraction scheme
renormalization.

For our case, the Abel-Plana formula is expressed as follows [23]:

$\sum_{n=0}^{\infty}f(n)=\int_{0}^{\infty}f(z)dz+\frac{1}{2}f(0)+i\int\frac{f(iz)-f(-iz)}{e^{2\pi
z}-1}dz,$ (16)

where, from Eq. (15), we have defined

$\displaystyle f(n)=$
$\displaystyle\frac{1}{2}\sqrt{\left(\frac{n\pi}{a}\right)^{2}+\bar{g}\left(\frac{n\pi}{a}\right)^{4}}.$
(17)

Upon using the above equation, the first integral in the r.h.s. of Eq. (16)
gives us the divergent free Minkowski contribution and, consequently, should
be discarded by means of the minimum subtraction scheme. Note also that
$f(0)=0$.

Hence, the finite part of the sum of (17) is provided by the third term in the
r.h.s. of the Abel-Plana formula (16), that is,

$\displaystyle E_{0}^{\rm ren}$
$\displaystyle=\frac{i}{2}\int_{0}^{\infty}\frac{f(iz)-f(-iz)}{e^{2\pi
z}-1}dz$
$\displaystyle=i\frac{a}{2\pi}\int_{0}^{\infty}\frac{du}{e^{2au}-1}\left[\sqrt{i^{2}u^{2}+\bar{g}u^{4}}-\sqrt{(-i)^{2}u^{2}+\bar{g}u^{4}}\right],$
(18)

where we have made the change of variable $u=\frac{z\pi}{a}$. In order to
solve the integral above one should note that for $u>\bar{g}^{-\frac{1}{2}}$
we obtain a null result since

$f(iz)-f(-iz)=0.$ (19)

However, for $u<\bar{g}^{-\frac{1}{2}}$, we have

$f(iz)-f(-iz)=2i\sqrt{u^{2}-\bar{g}u^{4}}.$ (20)

As a result, Eq. (18) becomes

$E_{0}^{\rm
ren}=-\frac{a}{\pi\bar{g}}\int_{0}^{1}\frac{\sqrt{y^{2}-y^{4}}}{e^{2xy}-1}dy,$
(21)

where we have made the new change of variable $y=\bar{g}^{\frac{1}{2}}u$ and
defined the dimensionless parameter $x=a\bar{g}^{-\frac{1}{2}}$. By using the
identity

$\frac{1}{e^{2xy}-1}=\sum_{k=1}^{\infty}e^{-2xyk},$ (22)

the integral (21) can be performed to provide [24]

$\displaystyle E_{0}^{\rm
ren}=-\frac{a}{\pi\bar{g}}\sum_{k=1}^{\infty}\left[\frac{2\mu_{k}-3\pi
I_{2}(\mu_{k})+3\pi L_{2}(\mu_{k})}{6\mu_{k}}\right],$ (23)

where $\mu_{k}=2xk$, $I_{n}(x)$ is the modified Bessel function of the first
kind and $L_{n}(x)$ is the modified Struve function. Alternatively, the
renormalized vacuum energy can be rewritten as

$E^{\rm ren}_{0}=\frac{x^{2}}{6\pi
a}-\frac{x}{4a}\sum_{k=1}^{\infty}\frac{1}{k}\left[L_{2}(2kx)-I_{2}(2kx)\right],$
(24)

where we have used the Riemann zeta function $\zeta(s)$ calculated at zero
argument, i.e., $\zeta(0)=-\frac{1}{2}$, to arrive at the first term of the
r.h.s. of (24). The expressions in Eqs. (21) and (24) are equally equivalent,
however, as we shall see below, the former is the best one for numerical
computations while the latter is more convenient to analyze the regime of
small values of the LIV parameter $\bar{g}$ or, equivalently,
$x=a\bar{g}^{-\frac{1}{2}}\gg 1$. This is an important regime to study since
it allows us to verify the consistence of our results by showing the recovery
of the standard result, i.e., the result in the absence of LIV.

The behaviour of the renormalized energy is described in Fig.3, which we have
plotted by using Eq. (21). The plot shows the curve traced by the energy, per
unit of $\frac{\pi}{24a}$, in terms of $x=a\bar{g}^{-\frac{1}{2}}$. We can see
that in the limit of small LIV parameter, i.e., $x=a\bar{g}^{-\frac{1}{2}}\gg
1$, the curve (red solid line) tends to the known result for the vacuum energy
in $(1+1)D$ [25], while in the opposite limit it tends to zero (blue dashed
line). However, from the experimental point of view, the limit $x\ll 1$ does
not make sense since it would indicate a strong signal of LIV. Thus, the part
of the curve that does make sense is indicated in red solid line.

Figure 3: The behaviour of Eqs. (21) and (24), expressed dimensionless as
$\frac{24\pi a}{\pi}E^{\rm ren}_{0}$, in terms of $x=a\bar{g}^{-\frac{1}{2}}$.

We can analytically investigate the limit where the Lorentz violating
parameter is small, $x\gg 1$, as previously mentioned. In order to obtain this
limit, we need to consider the approximation below, for large arguments, of
the modified Struve function [26]:

$\displaystyle
L_{\nu}(z)-I_{-\nu}(z)\simeq\frac{1}{\pi}\sum_{j=0}^{\infty}\frac{(-1)^{j+1}\Gamma(j+1/2)}{\Gamma(\nu+1/2-j)}\left(\frac{2}{z}\right)^{2j-\nu+1},$
(25)

where $\Gamma(s)$ is the Gamma function. Having in mind that
$I_{-n}(z)=I_{n}(z)$ ($n\in\mathbb{Z}$), the use of the above approximation in
Eq. (24) provides

$E^{\text{ren }}_{0}=\frac{x^{2}}{6\pi
a}-\frac{x}{4a\pi}\sum_{j=0}^{\infty}\frac{(-1)^{j+1}\Gamma(j+1/2)}{\Gamma(5/2-j)}\frac{\zeta(2j)}{x^{2j-1}}.$
(26)

In order to have a better understanding of what the above approximation gives,
let us consider it up to order $\bar{g}^{3}$, that is,

$\displaystyle E^{\text{ren
}}_{0}=-\frac{\pi}{24a}+\frac{\bar{g}\pi^{3}}{480a^{3}}+\frac{\bar{g}^{2}\pi^{5}}{4032a^{5}}+\frac{\bar{g}^{3}\pi^{7}}{7680a^{7}}+\mathcal{O}(\bar{g}^{4}).$
(27)

We can observe that the first term in the r.h.s. of the above expression is
the known result for the vacuum energy in $(1+1)D$, assuming Lorentz
invariance is preserved [25], as indicated in the plot of Fig.3. The other
terms are corrections to the latter due to the the high derivative operator,
which causes LIV, present in the Lagrangian (5). Therefore, the presence of
the LIV parameter allows us to access correction terms of high orders. These
corrections terms, as shown in Eq. (27) and confirmed by the plot in Fig.3,
makes the vacuum energy to increase. Again, the effect of the corrections
term, as it has happened also in the energy’s case, is to increases the vacuum
force making it less attractive.

##  V Casimir force

The Casimir force arising from the Dirichlet boundary condition is obtained by
taking the derivative of the energy with respect to the distance $a$. Hence,
starting from (21), which is an analytical solution without any approximation,
we obtain:

$\displaystyle F^{{\rm ren}}$ $\displaystyle=-\frac{\partial E^{\rm
ren}}{\partial a}$ $\displaystyle=\frac{x^{2}}{\pi
a^{2}}\int_{0}^{1}\frac{\sqrt{y^{2}-y^{4}}}{e^{2xy}-1}dy-\frac{x^{3}}{2\pi
a^{2}}\int_{0}^{1}\frac{y\sqrt{y^{2}-y^{4}}}{\text{sinh}^{2}(xy)}dy.$ (28)

Or equivalently, from Eq. (24),

$\displaystyle F^{{\rm ren}}=-\frac{x^{2}}{6\pi
a^{2}}+\frac{x^{2}}{4a^{2}}\sum_{k=1}^{\infty}\left\\{\left[L_{1}(2xk)-I_{1}(2xk)\right]+\left[L_{3}(2xk)-I_{3}(2xk)\right]\right\\}.$
(29)

Again, both expressions (28) and (29) are equivalent, however, Eq. (28) is
more convenient for numerical calculations while Eq. (29) for investigating
the analytical expression for the vacuum force in the regime of small LIV
parameter, i.e., $x\gg 1$. In Fig.4 we have plotted the dimensionless vacuum
force $\frac{24\pi a^{2}}{\pi}F^{\rm ren}$, in terms of $x$. One can see that
the force shares some very similar aspects with the energy, that is, in the
regime where the LIV parameter is very small $x\gg 1$, the force approaches
its standard contribution in the absence of Lorentz violation while in the
opposite regime $x\ll 1$ it goes to zero. However, just like in the energy
case, from the experimental point of view, it does no make sense to talk about
values of $x$ smaller than 1 $(\bar{g}>1)$, otherwise the LIV signal would be
very strong. Thus, the solid red line is the part of the curve which is
realistically acceptable. On the other hand, the blue dashed line is the part
of the curve for which $x<1$.

Figure 4: The behaviour of the vacuum force (28)-(29), expressed dimensionless
as $\frac{24\pi a^{2}}{\pi}F^{\rm ren}$, in terms of
$x=a\bar{g}^{-\frac{1}{2}}$.

Now, in order to obtain the limit $\bar{g}\to 0$ $(x\gg 1)$, we follow the
same procedure used for the Casimir energy, where we adopt the approximation
in Eq. (25).

$\displaystyle F^{{\rm ren}}=-\frac{x^{2}}{6\pi a^{2}}$
$\displaystyle+\frac{x^{2}}{4\pi
a^{2}}\left[\sum_{j=0}^{\infty}\frac{(-1)^{j+1}\Gamma(j+1/2)}{\Gamma(3/2-j)}\frac{\zeta(2j)}{x^{2j}}\right.$
$\displaystyle\left.+\sum_{j=0}^{\infty}\frac{(-1)^{j+1}\Gamma(j+1/2)}{\Gamma(7/2-j)}\frac{\zeta(2j-2)}{x^{2j-2}}\right].$
(30)

As a consequence, we are able to write the vacuum force up to order
$\bar{g}^{3}$ as

$\displaystyle F^{\text{ren
}}=-\frac{\pi}{24a^{2}}+\frac{\bar{g}\pi^{3}}{160a^{4}}+\frac{5\bar{g}^{2}\pi^{5}}{4032a^{6}}+\frac{7\bar{g}^{3}\pi^{7}}{7680a^{8}}+\mathcal{O}(\bar{g}^{4}).$
(31)

We can see that the first term in the r.h.s. of the above approximation is the
dominant standard contribution of the vacuum force in the absence of LIV [25]
while the other terms are corrections due to the presence of a high derivative
operator that induces the Lorentz invariance breaking.

##  VI Conclusion

In this work, we initiate by projecting the four-dimensional gauge field onto
a two-dimensional theory using a higher-derivative order operator. This
projection leads to modified equations of motion that rely on the nature of
the $n_{\mu}$ Lorentz invariance violation (LIV) parameter. Specifically, we
derive space-like and time-like equations of motion, as depicted in Eqs. (8a)
and (8b) respectively. Justifying the space-like approach, we demonstrate
subluminal behaviour through phase and velocity group analysis (Eq. (14)) and
energy behaviour, as depicted in Fig.2.

By using this dispersion relation, we obtain analytic expressions for the one-
dimensional Casimir energy (Eqs. (21) and (24)), which rely on the LIV
parameter $\bar{g}$. The behaviour of this energy is illustrated in Fig.3,
displaying convergence to a constant value. However, this closed expression
does not allow exploration of the limit where the LIV parameter is small
($\bar{g}\ll 1$). To address this, we approximate for large arguments,
obtaining an expression for the Casimir energy (Eq. (27)) that includes both
the usual term for $(1+1)D$ Casimir energy and a correction term. These
correction terms depend on the LIV parameter, and in the limit $\bar{g}\to 0$,
they vanish entirely, leaving only the conventional term for the one-
dimensional Casimir energy.

The Casimir force is also described by the system in Eqs. (28) and (29),
resulting in a closed expression similar to that obtained for the Casimir
energy. The behaviour of this analytic expression for the force is depicted in
Fig.4, illustrating the correct interval where the force (30) is valid.
However, like the Casimir energy (24), the analytic expression for the force
does not allow exploration of the limit $\bar{g}\ll 1$. Consequently, we
employ an approximation, leading to Eq. (30), which results in Eq. (31). This
equation reveals a repulsive term in the conventional one-dimensional Casimir
force, along with correction terms that can be adjusted according to the
desired precision. Regardless of the correction order, the force converges
negatively as $a$ tends to infinity and diverges positively as $a$ tends to
zero.

In conclusion, the combination of the high derivative order operator and
dimensional reduction yields analytic expressions for the Casimir energy and
Casimir force. Typically, obtaining closed-form expressions in the presence of
the LIV parameter requires approximations. However, our approach provides
analytic solutions, offering a more comprehensive understanding of the system.
Furthermore, the introduction of the LIV parameter allows us to manipulate
correction terms to achieve the best fit for both Casimir energy and force, as
previously demonstrated. Furthermore, by considering the temperature effect in
the present study, it is expected that additional phenomenological details
will emerge from the system.

###### Acknowledgements.

We would like to thank CNPq, CAPES and CNPq/PRONEX/FAPESQ-PB (Grant nos.
165/2018 and 015/2019), for partial financial support. K.E.L.F would like to
thank the Brazilian agency CNPq for financial support. M.A.A, F.A.B and E.P
acknowledge support from CNPq (Grant nos. 306398/2021-4, 309092/2022-1,
304290/2020-3). H.F.S.M is partially supported by CNPq under grant no.
311031/2020-0.

## References

  * [1] D. Colladay and V. A. Kostelecky, “Lorentz violating extension of the standard model,” Phys. Rev. D, vol. 58, no. 359, p. 116002, 1998.
  * [2] S. M. Carroll and H. Tam, “Aether Compactification,” Phys. Rev. D, vol. 78, p. 044047, 2008.
  * [3] P. Horava, “Quantum Gravity at a Lifshitz Point,” Phys. Rev. D, vol. 79, p. 084008, 2009.
  * [4] M. Pospelov and Y. Shang, “On Lorentz violation in Horava-Lifshitz type theories,” Phys. Rev. D, vol. 85, p. 105001, 2012.
  * [5] H. Belich, Jr., M. M. Ferreira, Jr., J. A. Helayel-Neto, and M. T. D. Orlando, “Dimensional reduction of a Lorentz and CPT violating Chern-Simons model,” Phys. Rev. D, vol. 67, p. 125011, 2003.  [Erratum: Phys.Rev.D 69, 109903 (2004)].
  * [6] D. Bazeia and R. Menezes, “Defect structures in Lorentz and CPT violating scenarios,” Phys. Rev. D, vol. 73, p. 065015, 2006.
  * [7] D. Bazeia, M. M. Ferreira, A. R. Gomes, and R. Menezes, “Lorentz-violating effects on topological defects generated by two real scalar fields,” Physica D, vol. 239, pp. 942–947, 2010.
  * [8] E. Passos and A. Y. Petrov, “Two-dimensional Lorentz-violating Chern-Simons-like action,” Phys. Lett. B, vol. 662, pp. 441–444, 2008\.  [Erratum: Phys.Lett.B 664, 318 (2008)].
  * [9] A. d. Souza Dutra and R. A. C. Correa, “Traveling solitons in Lorentz and CPT breaking systems,” Phys. Rev. D, vol. 83, p. 105007, 2011.
  * [10] R. Casana, E. S. Carvalho, and M. M. Ferreira, Jr, “Dimensional reduction of the CPT-even electromagnetic sector of the Standard Model Extension,” Phys. Rev. D, vol. 84, p. 045008, 2011.
  * [11] F. A. Brito, M. S. Guimaraes, E. Passos, P. Sampaio, and C. Wotzasek, “The 4D-2D projection of Lorentz-violating Myers-Pospelov QED,” Phys. Rev. D, vol. 86, p. 105036, 2012.
  * [12] E. Passos, C. A. G. Almeida, F. A. Brito, R. Menezes, J. C. Mota-Silva, and J. R. L. Santos, “Soliton solutions in two-dimensional Lorentz-violating higher derivative scalar theory,” Annals Phys., vol. 396, pp. 351–370, 2018.
  * [13] R. C. Myers and M. Pospelov, “Ultraviolet modifications of dispersion relations in effective field theory,” Phys. Rev. Lett., vol. 90, p. 211601, 2003.
  * [14] S. Liberati and L. Maccione, “Lorentz Violation: Motivation and new constraints,” Ann. Rev. Nucl. Part. Sci., vol. 59, pp. 245–267, 2009\.
  * [15] H. B. G. Casimir, “On the attraction between two perfectly conducting plates,” Indag. Math., vol. 10, no. 4, pp. 261–263, 1948.
  * [16] M. J. Sparnaay, “Measurements of attractive forces between flat plates,” Physica, vol. 24, pp. 751–764, 1958.
  * [17] S. K. Lamoreaux, “Demonstration of the Casimir force in the 0.6 to 6 micrometers range,” Phys. Rev. Lett., vol. 78, pp. 5–8, 1997.  [Erratum: Phys.Rev.Lett. 81, 5475–5476 (1998)].
  * [18] U. Mohideen and A. Roy, “Precision measurement of the Casimir force from 0.1 to 0.9 micrometers,” Phys. Rev. Lett., vol. 81, pp. 4549–4552, 1998.
  * [19] G. Bressi, G. Carugno, R. Onofrio, and G. Ruoso, “Measurement of the Casimir force between parallel metallic surfaces,” Phys. Rev. Lett., vol. 88, p. 041804, 2002.
  * [20] I. J. Morales Ulion, E. R. Bezerra de Mello, and A. Y. Petrov, “Casimir effect in Horava–Lifshitz-like theories,” Int. J. Mod. Phys. A, vol. 30, no. 36, p. 1550220, 2015.
  * [21] R. A. Dantas, H. F. S. Mota, and E. R. Bezerra de Mello, “Bosonic Casimir Effect in an Aether-like Lorentz-Violating Scenario with Higher Order Derivatives,” Universe, vol. 9, no. 5, p. 241, 2023.
  * [22] A. Erdas, “Magnetic corrections to the fermionic Casimir effect in Horava-Lifshitz theories,” Int. J. Mod. Phys. A, vol. 38, no. 22n23, p. 2350117, 2023.
  * [23] A. A. Saharian, “The generalized abel-plana formula with applications to bessel functions and casimir effect,” arXiv preprint arXiv:0708.1187, 2007\.
  * [24] I. S. Gradshteyn and I. M. Ryzhik, Table of integrals, series, and products.  Academic press, 2014.
  * [25] J. Ambjoern and S. Wolfram, “Properties of the vacuum. i. mechanical and thermodynamic,” Annals of Physics, vol. 147, no. 1, pp. 1–32, 1983.
  * [26] M. Abramowitz, I. A. Stegun, and R. H. Romer, “Handbook of mathematical functions with formulas, graphs, and mathematical tables,” 1988.