# Hodge Decomposition of Wall Shear Stress Vector Fields characterizing
Biological Flows

Faniry H. Razafindrazaka  Freie Universität Berlin, Germany  Institute for
Imaging Science and Computational Modelling in Cardiovascular Medicine,
Charité-Universitätsmedizin Berlin, Germany  Pavlo Yevtushenko  Institute for
Imaging Science and Computational Modelling in Cardiovascular Medicine,
Charité-Universitätsmedizin Berlin, Germany  Konstantin Poelke  Freie
Universität Berlin, Germany  Konrad Polthier  Freie Universität Berlin,
Germany  Leonid Goubergrits  Institute for Imaging Science and Computational
Modelling in Cardiovascular Medicine, Charité-Universitätsmedizin Berlin,
Germany

###### Abstract

A discrete boundary-sensitive Hodge decomposition is proposed as a central
tool for the analysis of wall shear stress (WSS) vector fields in aortic blood
flows. The method is based on novel results for the smooth and discrete Hodge-
Morrey-Friedrichs decomposition on manifolds with boundary and subdivides the
WSS vector field into five components: gradient (curl-free), co-gradient
(divergence-free), and three harmonic fields induced from the boundary, which
are called the center, Neumann and Dirichlet fields. First, an analysis of WSS
in several simulated simplified phantom geometries (duct and idealized aorta)
was performed in order to understand the impact of the five components. It was
shown that the decomposition is able to distinguish harmonic blood flow
arising from the inlet from harmonic circulations induced by the interior
topology of the geometry. Finally, a comparative analysis of 11 patients with
coarctation of the aorta (CoA) before and after treatment as well as 10
controls patient was done.

The study shows a significant difference between the CoA patients and the
healthy controls before and after the treatment. This means a global
difference between aortic shapes of diseased and healthy subjects, thus
leading to a new type of WSS-based analysis and classification of pathological
and physiological blood flow.

##### Keywords:

Hodge decomposition, vector fields, wall shear stress, computational fluid
dynamics, coarctation of the aorta

##  1 Introduction

Biological flows or hemodynamics of the cardiovascular system play an
important role in the genesis, progress and treatment of cardiovascular
pathologies including congenital or acquired diseases of the heart, heart
valves and vessels. This is because wall remodeling including wall thickness
and wall constitution is triggered by hemodynamics. The major hemodynamic
parameter describing an interaction between hemodynamics and a vessel wall,
which is covered by endothelial cells, is the wall shear stress (WSS). The WSS
is an area-normalized tangential force component of the blood flow acting on
the wall and/or endothelial cells. In turn, endothelial cells trigger and
modulate adaptation, inflammation and remodeling of the vessel wall as well as
a respective remodeling of the vessel lumen [1, 2]. Consequently, abnormal WSS
is considered an important local risk factor for a set of diseases or
pathological processes. These include, for example, atherosclerosis of carotid
arteries [3] or coronary artery disease [4], rupture risk of cerebral
aneurysms [5, 6] or abdominal aortic aneurysms [7], aortic dilatation [8], and
thrombus formation [9]. Furthermore, the analysis of WSS is also of great
interest for the study of the hemodynamic impact of a treatment or a change of
the hemodynamics caused by a certain treatment device. These studies include,
for example, an analysis of post-treatment flow conditions after a treatment
of cerebral aneurysms with a flow diverter [10] or a change of flow conditions
after an aortic valve replacement [11]. The use of WSS as a reliable
biomedical marker characterizing disease, disease progress or initiation and
also characterizing hemodynamic outcome of a treatment procedure is
challenging. This is because WSS is a surface bounded vector field that means
that WSS is described by a magnitude and direction varying in space and time.
This allows for a definition of a set of parameters, which were proposed
during the last years as hemodynamic risk parameters for endothelial
dysfunction and related wall remodeling. A characterization of WSS magnitude,
direction, time and space gradients as well as topological features results in
a relatively large set of parameters, which are well summarized in [12] and
[13]. The majority of studies investigating WSS in biological flows are
numerical studies investigating hemodynamics by an image-based computational
fluid dynamics approach [14]. 4D VEC MRI based assessment of the WSS is also
proposed in [15]. The primary source of data for the WSS analysis, however, is
CFD, since an accurate WSS assessment requires a high spatial resolution as
shown by mesh independence studies for CFD solutions [16].

Vector fields modelling fluid flow often tend to exhibit a complicated
behaviour on various scales and are hard to understand. This poses a
particular problem for clinical applications where the behaviour of blood flow
in vessels serves as an indicator for potential abnormalities. The classical
Helmholtz decomposition was a first step to classify and analyze vector fields
by decomposing them into a divergence-free component and a component having a
potential. With the advent of Hodge theory, Helmholtz’ results generalize to
decomposition rules for differential forms on closed manifolds in arbitrary
dimensions. Since then a tremendous amount of research—both on the theoretical
and on the applied side—has been carried out to include manifolds with
boundary, differential forms of Sobolev class and various flavours of Hodge-
type decomposition statements, see e.g. [17] for an overview of Hodge-type
decompositions and the survey [18].

An important landmark in this evolution is the $L^{2}$-orthogonal
decomposition of $k$-forms on manifolds with boundary as

$\Omega^{k}=\mathrm{d}\Omega_{D}^{k-1}\oplus\delta\Omega_{N}^{k+1}\oplus\mathrm{d}\Omega^{k-1}\cap\delta\Omega^{k+1}\oplus({\cal{H}}^{k}_{N}+{\cal{H}}_{D}^{k})$

where the spaces ${\cal{H}}^{k}_{N}$ and ${\cal{H}}_{D}^{k}$ of harmonic
Neumann and Dirichlet fields, respectively, reflect the absolute and relative
cohomology of the manifold. Specifically for vector fields, the first two
spaces in this decomposition correspond to divergent and rotational
irregularities in the interior of the geometry, whereas the latter three
spaces represent steady flows through the domain, as each field in these
spaces is harmonic. A fairly recent result [19] provides a further orthogonal
decomposition of these spaces into subspaces

${\cal{H}}^{k}_{N}={\cal{H}}_{N,\mathrm{co}}^{k}\oplus{\cal{H}}_{N,\partial\mathrm{ex}}^{k}\hskip
10.00002pt\text{ and }\hskip
10.00002pt{\cal{H}}_{D}^{k}={\cal{H}}_{D,\mathrm{ex}}^{k}\oplus{\cal{H}}_{D,\partial\mathrm{co}}^{k}$

which permits a precise distinction between harmonic flows induced by boundary
components, represented by the subspaces ${\cal{H}}_{N,\mathrm{co}}^{k}$ and
${\cal{H}}_{D,\mathrm{ex}}^{k}$, from those induced by the interior topology
of the manifold, represented by ${\cal{H}}_{N,\partial\mathrm{ex}}^{k}$ and
${\cal{H}}_{D,\partial\mathrm{co}}^{k}$.

For the numerical treatment of vector fields it is therefore important to seek
for a discretization which on the one hand provides a good approximation with
predictable error, and on the other hand preserves the structural
decomposition results from the smooth theory.

In this work we focus on a discretization by piecewise constant vector fields
(PCVF) resulting from CFD-based analyses of the blood flow. PCVFs are a very
intuitive and simple to implement approximation while at the same time a
concise theoretical framework has been developed in recent years, which
includes the aspects of convergence and structural consistency. The recent
work [20, 21] establishes a consistent discretization for PCVFs of the smooth
refined decomposition results for vector fields on surfaces with boundary, now
including distinguished subspaces for effective boundary analysis and control.
Previous to that, a first strategy for the analysis of vector fields is
provided by the decomposition in [22], with a convergence analysis on closed
surfaces in [23], and a discrete connection for PCVFs is proposed in [24],
both without an effective boundary control.

The aim of our study presented here is a proof of concept for the novel Hodge-
type decomposition analysis of the WSS vector fields for blood flows in
general and specifically for the aortic flow. The paper is structured as
follows: first, a theoretical analysis of each vector field component with
respect to a WSS vector field is given. Second, a detailed description of the
data acquisition and blood flow simulation is exposed. Finally, a statistical
analysis of several patients will summarize the results.

###  1.1 Discrete Hodge-type Decomposition

The most important results on discrete Hodge-type decompositions on simplicial
meshes concerning our application can be summarized by two fundamental
theorems: the traditional Hodge-Helmholtz decomposition decomposes vector
fields on closed surfaces into three components. In contrast, on surfaces with
boundary a refined decomposition is provided by the so-called Hodge-Morrey-
Friedrichs decomposition. The main ingredients of the discretization and the
related spaces are given in the appendix. These decompositions constitute the
building block of all analysis in the present work. For the theoretical
foundations see [20, 21].

######  Theorem 1.1 (Hodge-Helmholtz decomposition)

The space of piecewise constant vector fields $\Lambda^{1}(M_{h})$ on a closed
simplicial surface $M_{h}$ decomposes into an $L^{2}$-orthogonal sum of the
spaces of gradient fields, co-gradient fields and harmonic fields:

$\displaystyle\Lambda^{1}(M_{h})$ $\displaystyle=\nabla S_{h}\oplus J\nabla
S_{h}^{\ast}\oplus\left(H:=\ker\operatorname{curl}_{h}^{\ast}\cap\ker\operatorname{div}_{h}\right)$
$\displaystyle\mathcal{X}$
$\displaystyle=\underbrace{\nabla\varphi}_{\mathrm{\operatorname{curl}}_{h}^{\ast}\nabla\varphi=0}\oplus\underbrace{J\nabla\psi}_{\mathrm{\operatorname{div}}_{h}J\nabla
g=0}\oplus\underbrace{\mathcal{Y}}_{\mathrm{\operatorname{curl}}_{h}^{\ast}\mathcal{Y}=\mathrm{\operatorname{div}}_{h}\mathcal{Y}=0}$

Figure 1: Example of a Hodge-Helmholtz decomposition of a PCVF on a torus into
gradient, co-gradient and harmonic field.

The fields belonging to $\nabla S_{h}$ are free of turbulence and contain only
flow induced by sources and sinks. $J\nabla\psi$ is divergence-free and
contains the rotational part of the field (see figure 1). Furthermore, if
$M_{h}$ is homeomorphic to a sphere with $m$ boundaries, then the harmonic
fields can be decomposed into three components:

######  Theorem 1.2 (Hodge-Morrey-Friedrichs decomposition HMF)

On a surface $M_{h}$ homeomorphic to a sphere with $m$ boundaries, the space
of harmonic fields can be decomposed into Neumann fields, center fields, and
Dirichlet fields:

$\displaystyle\begin{split}\Lambda^{1}(M_{h})&=\nabla S_{0}\oplus J\nabla
S_{0}^{\ast}\oplus\nabla S_{h}\cap J\nabla S_{h}^{\ast}\oplus
J\nabla\psi\oplus\mathcal{H}_{N}\oplus\mathcal{H}_{C}\oplus\mathcal{H}_{D}\end{split}$
(1)

One of the main studies of this paper is to understand the nature of these
harmonic spaces on simulated CFD WSS vector fields. Intuitively the space
$\nabla S_{h}\cap J\nabla S_{h}^{\ast}$ of center vector fields behaves
similarly to the space of smooth vector fields forming an $\approx 45^{\circ}$
angle with the boundaries, the Neumann vector fields are orthogonal to the
boundaries, and the Dirichlet are parallel to the boundaries. By the
Pythagorian theorem it is

$\|\mathcal{X}\|^{2}=\|\nabla\varphi\|^{2}+\|J\nabla\psi\|^{2}+\|\mathcal{H}_{N}\|^{2}+\|\mathcal{H}_{C}\|^{2}+\|\mathcal{H}_{D}\|^{2}$

which enables a full quantification of the input vector fields according to
their decomposition components. Figure 2 shows an example of a HMF-
decomposition on the WSS of a simple flow on a cylinder. Notice how the field
is dominated by $\mathcal{H}_{D}$.

Figure 2: A HMF-decomposition of a perturbed WSS vecor field on a cylinder
into five components: gradient, co-gradient, center, Neumann, and Dirichlet
vector field.

###  1.2 WSS Component Analysis

In this section, we study each component of the HMF-decomposition with respect
to the WSS of several phantom as well as real patient models obtained from
CFD. The phantom models are either hand-designed or real patient models with
mathematical deformation and boundary conditions. The observations are used to
emphasize on possible changes of WSS encoded in each HMF-components with
respect to anatomy/topology of the geometry, and parameters used for blood
flow simulation.

####  1.2.1 Perturbed WSS

Consider a smooth cylinder with a WSS of a laminar flow. We add a moderate
amount of rotational noise to the vector field within the interval
$(-\alpha,\alpha)$ where $\alpha$ bounds the frequency of the noise. High
values of $\alpha$ correspond to high overall frequencies while small values
alter slightly the global smoothness of the flow. The HMF-decomposition shows
that the Dirichlet field $\mathcal{H}_{D}$ recovers the original field in its
unperturbed state, behaving similarly to a vector field denoising. The
increase of $\alpha$ decreases $\mathcal{H}_{D}$ and increases the co-gradient
field. Figure 3 is a quantitative comparison of each decomposition where
$\alpha$ varies from $0^{\circ}$ to $90^{\circ}$ degree. The diagram shows
that $\mathcal{H}_{D}$ is a good reference to understand the global structure
of the WSS. In general, harmonic fields depend only on the topology of the
shape, not the field. In figure 5 (second row), for example, $\mathcal{H}_{D}$
stays invariant even though the input velocity profile is changed.
Quantitatively half of the WSS component is Dirichlet. One logic behind this
is reflected in the nature of fluids, being mostly dominated by a laminar
component in order to move only in one direction.

Figure 3: Perturbation of a laminar WSS on a cylinder starting from
$0^{\circ}$ to $90^{\circ}$. An increase in angle deviation decreases the
Dirichlet component and increases the co-gradient components.

####  1.2.2 Coarctation analysis

Aortic coarctation is a common congenital heart disease. It represents a local
narrowing of the aortic vessel causing abnormal blood flow and pressure and in
the cardiovascular system. Generally, the WSS vector field of a pre- and post-
operative patient does not provide enough information about amelioration in
the patient blood flow. The HMF-decomposition enables us in a theoretical
setting to identify important changes between the two states. We took a
segmented MRI scan of a patient before and after operation, deformed the
coarctation linearly from pre to post and analyzed the WSS evolution during
the diffusion process. The simulation is performed with a plug profile and
settings given in section 2.1. The results are shown in Figure 4. We notice a
significant increase in the Dirichlet field amortized with a reduction in co-
gradient field. The improvement in the Dirichlet field component corresponds
to the improvement of the overall blood flow as proven previously. Notice how
the gradient, Neumann, and center fields remain almost unchanged. The nature
of these components is explained in the next sections.

Figure 4: Linear deformation of a pre to post intervention of a patient
stenosis. A constant input velocity profile is used for the simulation. An
increase in the Dirichlet component and a reduction in the co-gradient field
is observed within 10 frames of the deformation.

####  1.2.3 Plug vs MRI profile

The boundary conditions used in CFD are, most of the time, either a constant
input velocity field (plug) or velocity information acquired from 4D MRI scans
using specialized software and sequences. MRI profiles are noisy, however, and
sometimes the resulting WSS field looks more perturbed than a WSS field
obtained by a plug profile. Using the HMF-decomposition, one can classify
which components of the WSS are more affected by the inlet velocity profile. .

Figure 5: First row: Plug vs MRI input profile encoded in the center and
Neumann components. Second row: invariance of the Dirichlet component under
the change of input profile.

Figure 5 presents a WSS analysis of the same patient with a different inlet
profile. The two vector fields are very different but by analyzing each HMF-
decomposition component, one can see perturbations in the center and Neumann
fields. The Dirichlet fields in both cases are topologically the same (similar
streamline and same number of singularities). +Figure 6 is a comparison of
both inlet boundary conditions for ten control patients. The statistical
analysis (paired t-student test) of normally distributed data (Komogorov-
Smirnov test) found no significant difference for the gradient component
(p=0.571). +However, significantly (p=0.038) smaller Dirichlet components for
MRI-measured inlet velocity profiles accompanied with significantly larger co-
gradient (p=0.007), center (p=0.002) and Neumann (p=0.002) components of the
HMF-decomposition have been observed. Notice that the $L^{2}$-norms of the
center and Neumann components in both cases are relatively small compared to
the other components.

Figure 6: Comparison of the WSS of 10 healthy patients with MRI vs Plug inlet
velocity profile. The noise produced from the MRI can be identified by a
significant increase of center and Neumann fields.

####  1.2.4 Number of Branches

The following study shows the effect of branches on an idealized aorta.
Starting with a curved cylinder with zero branch, artificial branches are
successively added and a blood flow is simulated on each geometry using a plug
inlet profile. The results show that for this ideal situation the gradient
field increases with the number of branches. Geometrically, branches induce a
high curvature and hence more divergence. There are still several parameters
not taken into account such as tapering or twisted cylinders. The proposed
setup with the correct geometry can be used to analyze these extra cases.

Figure 7: Analysis of artificial aorta models with a varying number of
outlets. The second row shows the flow of the gradient field. The diagram
shows that the number of branching outlets is closely related to the gradient
and Dirichlet field.

####  1.2.5 Unsteady flow

Finally, the effect of time-varying flow boundary-conditions is examined in
this section. For this purpose, an unsteady CFD simulation of a whole cardiac
cycle performed earlier for a MICCAI CFD Challenge [25] was used. For the
analysis however, only the systolic part of the cardiac cycle is considered,
since the diastolic part shows only little to zero flow and therefore
negligible WSS in the aorta. Twenty time-points have been evaluated in total
and are presented in figure 8, along with inlet and outlet flow-curves.
Additionally, decomposed WSS vector field plots are presented for five time
points with a more detailed picture of the WSS distribution. As can be seen
from the decomposition at the various time points, the respective components
of the HMF-decomposition change over time, with a significant increase in co-
gradient component and decrease in Dirichlet component. Furthermore, it
appears that the variations of the HMF-components do not only arise from
variations of the flow rate (i.e. Reynolds Number) but also from acceleration
and deceleration effects. This can be seen by comparing two timepoints with
equivalent flow rates, as for example timepoints 3 and 17, both of which show
a flow rate of about 160 ml/s. Despite that, timepoint 17, at which the flow
is being decelerated, shows a significantly higher co-gradient and lower
gradient component than timepoint 3, where the flow is being accelerated.
Neumann and center components however remain at almost zero throughout the
whole systole.

---
Figure 8: Evolution of the WSS HMF-decomposition in a CFD simulation with an
unsteady flow. The diagram shows 20 time points of the simulation. The close-
ups are five phases from the 20 time points of the unsteady flow simulation
(green dots). The colors are relative to the min-max magnitude of each input
vector field.

##  2 Method

A diagram summarizing the analysis pipeline is given in figure 9. The
implementation of the HMF-decomposition is done following the iterative
$L^{2}$-projection approach [20, 22] with the discretization given in the
appendix. The choice of basis functions for each harmonic field subspace
follows [20]. Our system takes as input a mesh with a vector field and return
the five components decomposition of the vector field, assuming that the
surface fulfills the topological requirement. Our system takes as input a mesh
with a vector field and returns the five-term decomposition eq. (1) of the
vector field, assuming that the surface fulfills the topological requirement.
Our implementation is done in Java using the JavaView (www.javaview.de)
geometry processing package. The line integral convolution (LIC) implemented
in ZIBAmira 2015.28 (Zuse Institute Berlin) is used for the field
visualization. Maximum magnitude is colored with red while close to zero
vectors are colored in violet. Most of the data used in this paper is from
real patient biological models..

Figure 9: Analysis pipeline of WSS vector fields extracted from a simulated
model and analyzed via Hodge decomposition.

###  2.1 Data Input

##### MRI:

The HMF-decomposition analysis was done for WSS vector fields of the aorta
from a MRI based CFD analysis of the aortic flow. These are subdivided in two
groups: controls and coarctation of the aorta (CoA) patients before and after
treatment.

The study was carried out according to the principles of the Declaration of
Helsinki and approved by the local ethics committee. Written informed consent
was obtained from the participants and/or their legal guardians.

MRI examinations used to set boundary conditions for the CFD analysis were
performed using a 1.5 Tesla Achieva R5.1.8 MRI scanner with a five-element
cardiac phased-array coil (Philips Medical Systems, Best, The Netherlands).
MRI protocols including a routine three-dimensional anatomical imaging in end-
diastole are used to reconstruct the geometry of the aorta (3D MRI). The
sequence parameters used were: acquired voxel size $0.66\times 0.66\times 3.2$
mm, reconstructed voxel size $0.66\times 0.66\times 1.6$ mm, repetition time 4
ms, echo time 2 ms, flip angle $90^{\circ}$, number of signal averages 3.
Four-dimensional velocity-encoded MRI (4D VEC MRI) was used to capture the
flow data of the ascending aorta and the thoracic aorta (acquired voxel size
$2.5\times 2.5\times 2.5$ mm, reconstructed voxel size $1.7\times 1.7\times
2.5$ mm, repetition time 3.5 ms, echo time 2.2 ms, flip angle $5^{\circ}$, 25
reconstructed cardiac phases, number of signal averages 1). High velocity
encoding (3-6 m/s) in all three directions was used in order to avoid phase
wraps in the presence of valve stenosis or secondary flow. All flow
measurements were completed with automatic correction of concomitant phase
errors. These data were used to set inflow and outflow boundary conditions.

Figure 10: Segmented aorta reconstructed from MRI images and used for the CFD
simulation.

##### CFD:

CFD requires geometries. Geometries of human aortas were segmented and
reconstructed using ZIBAmira 2015.28 (Zuse Institute Berlin, Berlin, Germany)
according to the previous description [26]. Briefly, intensity based image
segmentation was done semi-automatically with an intense manual interaction.
Rough surface geometries were then generated from segmentations with a
subvoxel accuracy and subsequently smoothed using Meshmixer (v. 3.3, Autodesk,
Inc., San Rafael, USA). These procedures were described in more detail earlier
[26]. Figure 10 shows all aorta models used for our analysis.

With the exception of the unsteady case, all simulations were performed as
steady-state simulations of the peak-systolic aortic flow using STAR-CCM+ (v.
12.06, Siemens PLM Software, Plano, USA). Vessel walls were assumed to be
rigid and a no-slip boundary condition was applied at all walls. To model
turbulence observed in systolic aortic hemodynamics, a $k-\omega$ SST
turbulence model with a turbulence intensity of 5 percent at the velocity
inlet was used. Blood was modelled as a non-Newtonian fluid with a constant
density of $1050\,\,kg/m^{3}$ and a Carreau-Yasuda viscosity model [27].
Patient-specific flow rates as measured with GTFlow (GyroTools LLC, Zurich,
Switzerland) from 4D VEC MRI data were set at the LVOT inlet and the
descending aorta outlet. Furthermore, patient-specific velocity profiles at
peak systolic flow rate were extracted using MEVISFlow (v. 10.3, Fraunhofer
MEVIS, Bremen, Germany) and set as inlet boundary conditions. The used CFD
pipeline was earlier validated by a comparison with 4D VEC MRI measured
velocity fields as well as clinically validated against catheter measured
pressure drops in cases of CoA [28]. Furthermore, to validate results of our
simulations we compare velocity fields calculated by CFD against velocity
fields measured by 4D flow MRI, both visualized by velocity magnitude color
coded path lines [29]. Calculated wall shear stress values are in the range of
published results [30].

###  2.2 Statistical analysis

Statistical analysis of the Hodge Decomposition results was done using the
software package IBM SPSS Statistic, version 25 (IBM, USA). Measured data are
presented as mean and standard deviation (SD) for normally distributed data or
as a median with IQR. All data were tested for normality using the Kolmogorov-
Smirnov-Test. Depending on the results of the normality test, the T-student
test or Mann-Whitney-U test were used for the group comparison. Paired tests
were used to compare pre- and post-treatment results. A $p$ value $<0.05$ was
considered significant.

##  3 Results

Results of the HMF-decomposition analysis of 11 CoA patients before and after
treatment as well as 10 controls are illustrated in figure 11. The T-student
test found significantly lower gradient and significantly higher Dirichlet in
CoA cases before treatment vs. controls: 0.3 (SD=0.083) vs. 0.46 (SD=0.065)
gradient, and 0.54 (SD=0.125) vs. 0.42 (SD=0.065) Dirichlet. The co-gradient
in the CoA group was higher as in controls with median 0.119 IQR [0.069-0.147]
vs. median 0.086 IQR [0.065-0.098], approaching significance (Mann-Whitney
test, p=0.061). Overall significant reduction (paired Wilcoxon test, p=0.041)
in co-gradient has been observed from pre (median 0.119 IQR [0.069-0.147]) to
post intervention (median 0.070 IQR [0.064-0.113]) as expected from the
theoretical experimentation exposed previously. However, no significant
changes in the major flow descriptors of gradient (p=0.174) and Dirichlet
(p=0.073) were found between pre and post treatment WSS vector fields (paired
T-Student test).

Figure 11: Comparison of the WSS of 11 patients before and after intervention,
and 10 healthy patients. An improvement in gradient and Dirichlet together
with a reduction in co-gradient is observed.

##  4 Discussion

Our first results on the application of a discrete HMF-decomposition analysis
of the aortic flow and especially an analysis of the WSS vector fields of the
coarctation of the aorta (congenital narrowing of the aorta) disease revealed
great potential for computational biofluid mechanics. Based on the results
shown in figure 11 we suppose that the HMF-decomposition analysis allows us to
find anatomical shapes forming pathological hemodynamics before disease
progress becomes symptomatic.

Our findings show an added value of the HMF-decomposition analysis if compared
with the usually used analysis of WSS vector fields by visualization or
quantification of time- and surface-averaged WSS values, areas with low WSS
values (e.g. WSS values below 0.5 Pa) or areas with high OSI as well as an
analysis of WSS critical points [6, 30, 12]. This approach, however, does not
allow, for example, a quantitative analysis of two different abnormal WSS
vector fields or a quantitative analysis of different impact factors (boundary
conditions) forming abnormal hemodynamics.

The results shown in figure 11 together with the theoretical analysis on ideal
models raise several open questions. Could a pathological anomaly such us
stenosis present in the aorta be identified by its amount of WSS co-gradient?
Control healthy patients have less co-gradient field. The pre vs post
operative patient also show a significant improvement in co-gradient field. An
objective classification has not been achieved with our current analysis
because of the limited number of patient models. Nevertheless the theoretical
deformation shown in figure 4 suggests that it should generally be the case.

The current analysis is focusing only on the studying the differences in WSS
vector fields shown by the HMF-decomposition due to treatment aiming at
restoring the stenosed region towards a physiological diameter. The
differences between diseased and control groups aiming to identify hemodynamic
reasons for the development of a pathological anatomy are emphasized. Future
research could be also focused on the impact or decomposition of hemodynamic
and/or morphometric boundary conditions on the resulting HMF WSS vector field
decomposition. This is, however, a challenging task since the hemodynamics
depend on a set of non-linear effects of all boundary conditions including the
flow rates distributions, vessel curvature, branching topology and others.

A perfect flow would have a pure WSS Dirichlet field, but due to branches and
taperings in the shape, gradient and co-gradient components are also present.
On the one hand, there are higher gradient than Dirichlet components in the
control group, on the other hand there are higher Dirichlet than gradient
components in the pre and post operative groups. The theoretical analysis on
the number of branches shows that the nature of the gradient field may change
and become dominant. Understanding the correlation between the gradient and
the Dirichlet field will be a good direction for future research.

We applied the HMF-decomposition first to analyse WSS vector fields, since WSS
is a known risk factor for the genesis and progress of pathological processes
associated with an interaction between blood flow and vessel wall. The
analysis allows for an integral characterization of the WSS distribution.
However, it does not replace an analysis of WSS magnitudes, which are also
associated with abnormal blood flow conditions: regions with low WSS promote
development of atherosclerosis and thrombus formations, whereas high WSS could
cause an injury of endothelial cells. As part of our study we investigated the
impact of side branches, degree of stenosis and/or treatment procedure, inlet
flow profile boundary conditions and the impact of laminar flow disturbances
on WSS vector fields as characterized by the HMF-decomposition.

The HMF-decomposition analysis of simulated WSS vector fields was based on the
Reynolds-averaged Navier-Stokes (RANS) solver using the $k-\omega$ SST
turbulence model. However, flow simulations of hemodynamics allowing
assessment of pressure and velocity fields and hence WSS are not limited to
the RANS CFD. The Lattice-Boltzman method (LBM), Large-Eddy simulations (LES)
or Smoothed-Particle Hydrodynamics (SPH) are possible CFD alternatives. For
example, LES is supposed to be better suited in order to simulate accurately
transition to turbulence and to assess turbulent structures [31]. Finally, the
choice of the CFD approach should be done based on validation studies
comparing simulation results vs. in vivo measurements [32]. The HMF-
decomposition analysis is, however, independent from the CFD approach.

Further possible and planned studies include, for example, an analysis of
pulsatile flows, analysis of flow differences due to different turbulence
models, the extension of an analysis to other parts of circulation (e.g.
coronary arteries, carotid bifurcations or cerebral vessels) and other
diseases (e.g. abdominal aortic aneurysms, cerebral aneurysms or coronary
artery disease).

Summarizing our results, the HMF-decomposition is able to support (1) basic
research of the flow mediated disease, (2) predictive computational modelling
of the treatment procedure as well as (3) quantitative analysis of the
hemodynamic treatment outcome. Altogether, it supports a clinical translation
of the computational modelling approach.

##  5 Conclusion

The novel discrete Hodge-Morrey-Friedrichs decomposition was for the first
time applied to analyze the WSS vector fields of simulated patient-specific
aortic blood flows. The approach seems to be a powerful tool to distinguish
between pathological and physiologic blood flows, +and to characterize the
impact of inflow boundary conditions as well as the impact of a treatment.

## Appendix A

In this appendix, we give a brief introduction to the calculus on discrete
surfaces. Only the most relevant notions necessary to understand the discrete
Hodge decomposition are given. A complete overview can be found in [22].

###  A.1 Simplicial Surfaces

A 2-dimensional simplicial surface $M_{h}$ is a set of triangles glued at
their edges with a manifold structure. In finite element analysis, this type
of discrete surface is called triangle mesh. For actual FEM computations on
such meshes one often uses the space $S_{h}$ of linear Lagrange functions, or
the space $S^{*}_{h}$ of Crouzeix-Raviart functions. They are defined by

$\displaystyle S_{h}$
$\displaystyle:=\left\\{\varphi:M_{h}\rightarrow{\mathbb{R}}\text{
}\left|\text{ }\varphi_{|T}\text{ is linear on each triangle }T\text{, and
globally continuous}\right.\right\\}$ $\displaystyle S_{h}^{\ast}$
$\displaystyle:=\left\\{\psi:M_{h}\rightarrow{\mathbb{R}}\text{ }\left|\text{
}\psi_{|T}\text{ is linear, and continuous at edge midpoints}\right.\right\\}$

Figure 12: Examples of a function $\varphi\in S_{h}$ together with
$\nabla\varphi$, $J\nabla\varphi$ (left) and $\psi\in S_{h}^{*}$ with
$\nabla\psi$, and $J\nabla\psi$ (right) defined over the triangle $T$.

A geometric realization of example functions on $S_{h}$ and $S^{*}_{h}$ is
shown in figure 12. Two additional subspaces $S_{0}\subset{S_{h}}$ and
$S_{0}^{\ast}\subset{S_{h}^{\ast}}$ for surfaces with boundary are given by

$\displaystyle S_{0}$ $\displaystyle:=\left\\{\varphi\in S_{h}\text{ }|\text{
}\varphi(v)=0\text{ for all boundary vertices }v\right\\}$ $\displaystyle
S_{0}^{\ast}$ $\displaystyle:=\left\\{\psi\in S_{h}^{\ast}\text{
}|\,\,\psi(m_{e})=0\text{ for all boundary edge mid-points }m_{e}\right\\}.$

The _gradient field_ $\nabla\varphi$ of a function $\varphi\in S_{h}$ or
$S_{h}^{\ast}$ is a constant tangent vector in each triangle. The _co-gradient
field_ $J\nabla\varphi$ is obtained by a rotation $J$ of the gradient
$\nabla\varphi$ by $\frac{\pi}{2}$ in each triangle (see figure 12). The idea
of having functional spaces is a common technique in finite element analyses
to solve complicated partial differential equations. For example a temperature
map $u$ which assigns a scalar value to each vertex of $M_{h}$ is an element
of $S_{h}$. It can be expressed with respect to the nodal basis functions
$(\varphi_{i})_{i}$ of $S_{h}$, i.e $u=\sum_{i}u_{i}\varphi_{i}$, where
$\varphi_{i}$ is the Kronecker delta, $\varphi_{i}(v_{j})=1$ if $i=j$ and $0$
otherwise, for a vertex $v_{j}$ of $M_{h}$. Then, the gradient of an arbitrary
function in $S_{h}$ can simply be expressed as a linear combination of the
$\nabla\varphi_{i}$’s.

###  A.2 Vector fields on Simplicial Surfaces

######  Definition A.1 (PCVF)

The space of _piecewise constant tangential vector fields_
$\Lambda^{1}(M_{h})$ on a $2$-dimensional simplicial surface
$M_{h}\subset{\mathbb{R}^{n}}$ is given by:

$\Lambda^{1}(M_{h}):=\left\\{\mathcal{X}:M_{h}\rightarrow TM_{h}\ \left|\
\mathcal{X}_{|\text{triangle }T}\text{ is a constant tangent vector in
}T\right.\right\\}.$

Here, $TM_{h}$ denotes the (piecewise) tangent bundle of $M_{h}$. The gradient
field $\nabla\varphi$ introduced previously is an example of a tangential
vector field.

######  Definition A.2 ($L^{2}$-product)

The $L^{2}$-product of two vector fields, $\mathcal{X}=(\mathcal{X}_{T})_{T\in
M_{h}}$ and $\mathcal{Y}=(\mathcal{Y}_{T})_{T\in M_{h}}$, where
$\mathcal{X}_{T},\mathcal{Y}_{T}$ are tangent vectors in the triangle $T$, is
defined by the area-weighted Euclidean sum

$\left<\mathcal{X},\mathcal{Y}\right>_{L^{2}}=\sum_{T\in
M_{h}}\left<\mathcal{X}_{T},\mathcal{Y}_{T}\right>\text{Area(T)}.$

In particular, two vector fields $\mathcal{X}$ and
$\mathcal{Y}\in\mathcal{X}_{h}$ are _$L^{2}$ -orthogonal_ if
$\left<\mathcal{X},\mathcal{Y}\right>_{L^{2}}=0$. A vector field subspace
$A\subseteq\Lambda^{1}(M_{h})$ is the _$L^{2}$ -orthogonal decomposition_ of
two subspaces $B,C\subseteq A$ (written $A=B\oplus C$) if every
$\mathcal{X}\in A$ can be written uniquely as a sum
$\mathcal{X}=\mathcal{Y}+\mathcal{Z}$ with $\mathcal{Y}\in B$ and
$\mathcal{Z}\in C$, and furthermore $\left<\mathcal{Y},\mathcal{Z}\right>=0$.
The sum of two vector fields is a new vector field obtained by the sum of the
components.

###  A.3 Discrete Calculus

######  Definition A.3 (Discrete Curl)

The discrete curl of a vector field $\mathcal{X}=(\mathcal{X}_{T})_{T\in
M_{h}}$ at a vertex $p$ and an edge midpoint $m_{e}$ of $M_{h}$ is computed by

$\displaystyle\operatorname{curl}_{h}\mathcal{X}(p)$
$\displaystyle:=\frac{1}{2}\oint\limits_{\partial\operatorname{star}p}\mathcal{X}=\frac{1}{2}\sum_{i=1}^{k}\left<\mathcal{X}_{|T_{i}},e_{i}\right>$
$\displaystyle\operatorname{curl}_{h}^{\ast}\mathcal{X}(m_{e})$
$\displaystyle:={\oint\limits_{\partial\operatorname{star}e}}\mathcal{X}=-\left<\mathcal{X}_{|T_{1}},e\right>+\left<\mathcal{X}_{|T_{2}},e\right>$

where the $e_{i}$’s are the edges of the oriented boundary of $\text{star }p$,
the $T_{i}$’s the triangles adjacent to $p$ and $e$ the edge with midpoint
$m_{e}$ (see figure 13).

######  Definition A.4 (Discrete Divergence)

The discrete divergence of a vector field $\mathcal{X}$ at a vertex $p$ and an
edge midpoint $m_{e}$ of $M_{h}$ is computed by

$\displaystyle\operatorname{div}_{h}\mathcal{X}(p)$
$\displaystyle:=\frac{1}{2}{\displaystyle\oint\limits_{\partial\operatorname{star}p}}\left\langle\mathcal{X},\mathbf{n}\right\rangle
ds=-\frac{1}{2}\sum_{i=1}^{k}\left\langle\mathcal{X}_{|T_{i}},Je_{i}\right\rangle$
$\displaystyle\operatorname{div}_{h}^{\ast}\mathcal{X}(m_{e})$
$\displaystyle:={\displaystyle\oint\limits_{\partial\operatorname{star}e}}\left\langle\mathcal{X},\mathbf{n}\right\rangle
ds=\left\langle\mathcal{X}_{|T_{1}},J_{|T_{1}}e\right\rangle+\left\langle\mathcal{X}_{|T_{2}},J_{|T_{2}}e\right\rangle$

where $\mathbf{n}$ is the outer unit normal along
$\partial\operatorname{star}p$ resp. $\partial\operatorname{star}e$. Discrete
rotation and divergence are related by
$\operatorname{curl}_{h}J\mathcal{X}=\operatorname{div}_{h}\mathcal{X}$ and
$\operatorname{curl}_{h}^{\ast}J\mathcal{X}=\operatorname{div}_{h}^{\ast}\mathcal{X}$,
compare Figure 13.

Figure 13: Computation of the curl and divergence of a vector field at a point
on $\nabla S_{h}$ and $J\nabla S_{h}^{*}$.

######  Definition A.5 (Dirichlet and Neumann field)

A vector field $\mathcal{X}$ is a Dirichlet field (resp. Neumann field) if
$\mathcal{X}_{|T\in\partial M_{h}}$ is orthogonal (resp. “almost” parallel) to
the boundary edge of $T$.

In the discrete case, the definition of Neumann fields is subtle due to
technical properties of the chosen function spaces $S_{h}$ and $S_{h}^{*}$.
They are not strictly parallel along the boundary as one might expect from the
smooth case, but can deviate slightly. However, they are overall mostly
parallel, so for simplicity one may imagine them as being just parallel, in
perfect duality to the definition of Dirichlet fields. For technical details
we refer the reader to [20, Sec. 3.1].

The harmonic Dirichlet field $\mathcal{H}_{D}$ is for example a divergence-
free and a curl-free field orthogonal to $\partial M_{h}$. Note that Dirichlet
fields and Neumann fields do not exists on a closed surface. One can
nevertheless define them using hard directional constraints on certain
features of the underlying surface, e.g sharp features.

##  2 Statement

##### Ethics

The MRI data of CoA patients and volunteers were acquired in frames of a
study, which was carried out according to the principles of the Declaration of
Helsinki and approved by the local Ethics Committee at Charité-
Universitätsmedizin Berlin. Written informed consent was obtained from the
participants and/or their guardians. The clinical study has been registered
with `ClinicalTrials.gov`

##### Data accessibility

There is not data deposition applicable for the paper.

##### Funding

This research was carried out in the framework of MATHEON supported by the
Einstein Foundation Berlin within the ECMath Project CH18.

##### Competing interests

The authors declare no competing interests.

##### Authors’ contributions

F.H. Razafindrazaka, L.Goubergrits, and K. Polthier developed the project from
conception to design. P. Yevtushenko was responsible in acquiring the patient
models, reconstructing the geometry, and doing the numerical simulations. K.
Poelke, and K. Polthier designed and implemented the Hodge decomposition
algorithm. F.H. Razafindrazaka analysed the data using the method. L.
Goubergrits, and P. Yevtushenko interpreted the results with . F.H.
Razafindrazaka, L. Goubergrits, P. Yevtushenko, and K. Poelke wrote the
manuscript. All authors gave final approval for publication.

## References

  * [1] Chatzizisis YS, Coskun A, Jonas M, Edelman E, Feldman C, Stone P.  Role of Endothelial Shear Stress in the Natural History of Coronary Atherosclerosis and Vascular Remodeling: Molecular, Cellular, and Vascular Behavior.  Journal of the American College of Cardiology. 2007;49(25):2379 – 2393\.
  * [2] Chiu J, Chien S.  Effects of Disturbed Flow on Vascular Endothelium: Pathophysiological Basis and Clinical Perspectives.  Physiological Reviews. 2011;91(1):327–387.
  * [3] Zhang B, Gu J, Qian M, Niu L, Ghista D.  Study of correlation between wall shear stress and elasticity in atherosclerotic carotid arteries.  BioMedical Engineering OnLine. 2018 Jan;17(1):5.
  * [4] Soulis J, Fytanidis D, Papaioannou V, Giannoglou G.  Wall shear stress on LDL accumulation in human RCAs.  Medical Engineering & Physics. 2010;32(8):867 – 877.
  * [5] Meng H, Wang Z, Hoi Y, Gao L, Metaxa E, Swartz D, et al.  Complex Hemodynamics at the Apex of an Arterial Bifurcation Induces Vascular Remodeling Resembling Cerebral Aneurysm Initiation. 2007;38(6):1924–1931.
  * [6] Boussel L, Rayz V, McCulloch C, Martin A, Acevedo-Bolton G, Lawton M, et al.  Aneurysm Growth Occurs at Region of Low Wall Shear Stress Patient-Specific Correlation of Hemodynamics and Growth in a Longitudinal Study. 2008 08;39:2997–3002.
  * [7] Stevens R, Grytsan A, Biasetti J, Roy J, Lindquist Liljeqvist M, Gasser C.  Biomechanical changes during abdominal aortic aneurysm growth.  PLOS ONE. 2017 11;12(11):1–16.
  * [8] Saikrishnan N, Mirabella L, Yoganathan A.  Bicuspid aortic valves are associated with increased wall and turbulence shear stress levels compared to trileaflet aortic valves.  Biomechanics and Modeling in Mechanobiology. 2015 Jun;14(3):577–588.
  * [9] Lozowy R, Kuhn D, Ducas A, Boyd A.  The Relationship Between Pulsatile Flow Impingement and Intraluminal Thrombus Deposition in Abdominal Aortic Aneurysms.  Cardiovascular Engineering and Technology. 2017 Mar;8(1):57–69.
  * [10] Jing L, Zhong J, Liu J, Yang X, Paliwal N, Meng H, et al.  Hemodynamic Effect of Flow Diverter and Coils in Treatment of Large and Giant Intracranial Aneurysms.  World Neurosurgery. 2016;89:199 – 207.
  * [11] von Knobelsdorff-Brenkenhoff F, Trauzeddel R, Barker A, Gruettner H, Markl M, Schulz-Menger J.  Blood flow characteristics in the ascending aorta after aortic valve replacement, a pilot study using 4D-flow MRI.  International Journal of Cardiology. 2014;170(3):426 – 433.
  * [12] Arzani A, Shadden S.  Wall shear stress fixed points in cardiovascular fluid mechanics.  Journal of Biomechanics. 2018;73:145 – 152.
  * [13] Goubergrits L, Schaller J, Kertzscher U, Woelken T, Ringelstein M, Spuler A.  Hemodynamic impact of cerebral aneurysm endovascular treatment devices: coils and flow diverters.  Expert Review of Medical Devices. 2014;11(4):361–373.
  * [14] Morris P, Narracott A, von Tengg-Kobligk H, Silva Soto D, Hsiao S, Lungu A, et al.  Computational fluid dynamics modelling in cardiovascular medicine.  Heart. 2015;.
  * [15] Rodriguez-Palomares J, Dux-Santoy L, Guala A, Kale R, Maldonado G, Teixidó-Turà G, et al.  Aortic flow patterns and wall shear stress maps by 4D-flow cardiovascular magnetic resonance in the assessment of aortic dilatation in bicuspid aortic valve disease.  Journal of Cardiovascular Magnetic Resonance. 2018 Apr;20(1):28.
  * [16] Prakash S, Ethier C.  Requirements for Mesh Resolution in 3D Computational Hemodynamics. 2001 05;123:134–44.
  * [17] Schwarz G.  Hodge decomposition: a method for solving boundary value problems.  Lecture notes in mathematics. Springer; 1995.
  * [18] Bhatia H, Norgard G, Pascucci V, Bremer P.  The Helmholtz-Hodge Decomposition - A Survey.  IEEE Transactions on Visualization and Computer Graphics. 2013;19(8):1386–1404.
  * [19] Shonkwiler C.  Poincaré duality angles and the Dirichlet-to-Neumann operator.  Inverse Problems. 2013;29(4).
  * [20] Poelke K, Polthier K.  Boundary-aware hodge decompositions for piecewise constant vector fields.  Computer-Aided Design. 2016;78:126 – 136.  {SPM} 2016.
  * [21] Poelke K.  Hodge-Type Decompositions for Piecewise Constant Vector Fields on Simplicial Surfaces and Solids with Boundary.  Freie Universität Berlin; 2017.
  * [22] Polthier K, Preuss E.  Identifying Vector Field Singularities using a Discrete Hodge Decomposition.  In: Hege HC, Polthier K, editors. Visualization and Mathematics III. Springer Verlag; 2003. p. 113–134.
  * [23] Wardetzky M.  Discrete Differential Operators on Polyhedral Surfaces - Convergence and Approximation.  Freie Universität Berlin; 2006.
  * [24] Azencot O, Ovsjanikov M, Chazal F, Ben-Chen M.  Discrete Derivatives of Vector Fields on Surfaces–An Operator Approach.  ACM Transactions on Graphics (TOG). 2015;34(3):29.
  * [25] Schaller J, Goubergrits L, Yevtushenko P, Kertzscher U, Riesenkampff E, Kuehne T.  Hemodynamic in Aortic Coarctation Using MRI-Based Inflow Condition.  Lecture Notes in Computer Science. 2014;8330:65 – 73.
  * [26] Hellmeier F, Nordmeyer S, Yevtushenko P, Bruening J, Berger F, Kuehne T, et al.  Hemodynamic Evaluation of a Biological and Mechanical Aortic Valve Prosthesis Using Patient-Specific MRI-Based CFD.  Artificial Organs. 2017;42(1):49–57.
  * [27] Safoora K, Mahsa D, Paritosh V, Mitra D, Bahram D, Payman J.  Effect of rheological models on the hemodynamics within human aorta: CFD study on CT image-based geometry.  Journal of Non-Newtonian Fluid Mechanics. 2014;207:42 – 52.
  * [28] Goubergrits L, Riesenkampff E, Yevtushenko P, Schaller J, Kertzscher U, Hennemuth A, et al.  MRI-based computational fluid dynamics for diagnosis and treatment prediction: Clinical validation study in patients with coarctation of aorta.  Journal of Magnetic Resonance Imaging. 2015;41(4):909–916.
  * [29] Goubergrits L, Mevert R, Yevtushenko P, Schaller J, Kertzscher U, Meier S, et al.  The impact of MRI-based inflow for the hemodynamic evaluation of aortic coarctation.  Annals of Biomedical Engineering. 2013;41:2575–2587.
  * [30] LaDisa JF, Figueroa CA, Vignon-Clementel IE, Kim HJ, Xiao N, Ellwein LM, et al.  Computational simulations for aortic coarctation: representative results from a sampling of patients.  Journal of Biomechanical Engineering. 2011;133:091008–1 – 9.
  * [31] Vergara C, LeVan D, Quadrio M, Formaggia L, Domanin M.  Large eddy simulations of blood dynamics in abdominal aortic aneurysms.  Medical Engineering and Physics. 2017;47:38 – 46.
  * [32] Miyazaki S, Itatani K, Furusawa T, Nishino T, Sugiyama M, Takehara Y, et al.  Validation of numerical simulation methods in aortic arch using 4D Flow MRI.  Heart Vessels. 2017;32:1032 – 1044.


# Nonlinear spectroscopy of bound states in perturbed Ising spin chains

GiBaik Sim  Department of Physics TQM, Technische Universität München, $\&$
James-Franck-Straße 1, D-85748 Garching, Germany  Munich Center for Quantum
Science and Technology (MCQST), 80799 Munich, Germany  Johannes Knolle
Department of Physics TQM, Technische Universität München, $\&$ James-Franck-
Straße 1, D-85748 Garching, Germany  Munich Center for Quantum Science and
Technology (MCQST), 80799 Munich, Germany  Blackett Laboratory, Imperial
College London, London SW7 2AZ, United Kingdom  Frank Pollmann  Department of
Physics TQM, Technische Universität München, $\&$ James-Franck-Straße 1,
D-85748 Garching, Germany  Munich Center for Quantum Science and Technology
(MCQST), 80799 Munich, Germany

###### Abstract

We study the nonlinear response of non-integrable 1D spin models using
infinite matrix-product state techniques. As a benchmark and demonstration of
the method, we first calculate the 2D coherent spectroscopy for the exactly
soluble ferromagnetic transverse field Ising model, where excitations are
freely moving domain-walls. We then investigate the distinct signatures of
confined bound states by introducing a longitudinal field and observe the
emergence of strong non-rephasing like signals. To interpret the observed
phenomena, we use a two-kink approximation to perturbatively compute the 2D
spectra. We find good agreement in comparison with the exact results of the
infinite matrix-product state method in the strongly confined regime. We
discuss the relevance of our results for quasi$-$1D Ising spin chain materials
such as $\mathrm{CoNb}_{2}\mathrm{O}_{6}$.

Introduction. Spectroscopic tools have played a crucial role in our
understanding of complex quantum systems [1]. However, we can only partially
measure their correlations with existing tools. Terahertz 2D coherent
spectroscopy, one of developing spectroscopic tools [2, 3, 4], stands out as a
technique for a deeper understanding of strongly-correlated condensed matter
systems. It probes the nonlinear optical response which has been used to
identify the ground state symmetry of magnetic [5] and superconducting [6, 7]
materials, geometric phase in topological materials [8, 9, 10], novel ground
states in correlated systems [11, 12], and quasiparticle decay processes in a
disordered system [13]. Besides, recent experimental advances with terahertz
sources put the technique in a proper energy range to study rotational
dynamics in molecules [14], spin waves in conventional magnets [14], and
exotic excitations in quantum magnets [15, 16, 17, 18, 19].

In contrast to more common 1D spectroscopy, the 2D extension unravels not only
the optical excitations but also their interplay [4, 3]. The advantage of this
experimental technique has been widely adopted by chemists to reveal the
structure of complex molecules with great success. However, such achievements
rely on powerful numerical methods that help to interpret complicated
experimental data starting from concrete microscopic models [20, 21, 22]. In
this regard, it is desirable to develop an efficient numerical platform for
future 2D spectroscopy experiments on quantum magnets similar to the
successful use of matrix product state (MPS) techniques for conventional 1D
spectroscopy [23, 24, 25, 26, 27, 28, 29]. However, the calculation of
nonlinear response is less explored and the need for multiple time evolutions
makes it much more challenging.

In this work, we propose an efficient numerical tool using infinite MPS (iMPS)
and study the nonlinear response of 1D spin model. To benchmark the method, we
first focus on the 1D transverse field Ising model (TFIM) whose nonlinear
response can be analytically calculated using Jordan-Wigner (JW)
transformations [15]. Motivated by the quasi-1D structure of
$\mathrm{CoNb}_{2}\mathrm{O}_{6}$ – one of the best material example of an
Ising chain magnet – (albeit with more complicated magnetic interactions [28,
29]) we then include longitudinal field terms which capture the effects of
inter-chain interactions and lead to the emergence of confined bound state
excitations [30, 31, 32]. As a consequence, new signals appear in 2D
spectroscopy which include strong non-rephasing and rephasing like peaks. Our
results from the iMPS method are furthermore corroborated by perturbative
calculations starting from the projected two-kink (TK) low energy subspace. We
find quantitative agreement in the strongly confined regime, which allows us
to understand the origin of sharp peaks in the 2D spectrum as transitions
between bound states.

Model. We first introduce the 1D TFIM

$\displaystyle
H_{0}=-J\sum_{n}\sigma_{n}^{z}\sigma_{n+1}^{z}-h^{x}\sum_{n}\sigma_{n}^{x}$
(1)

with $J,h_{x}>0$. For $h_{x}<h^{c}_{x}=J$, it stabilizes a doubly degenerate
ferromagnetic ground state polarized along the easy axis $\hat{z}$. When
$h_{x}>h^{c}_{x}$, the system has a unique paramagnetic ground state. In the
ferromagnetic regime, the experimental excitation, i.e., a local spin flip,
splits into two freely moving kinks (domain-walls) between two degenerate
states. In the context of non-linear spectroscopy these fractionalized
excitations have been shown to be manifest as sharp signatures in the third
order magnetic susceptibilities [15].

We now include a longitudinal field and focus on the Hamiltonian given by
$H=H_{0}-h^{z}\sum_{n}\sigma^{z}_{n}$ with $J,h_{x},h_{z}>0$. In the
ferromagnetic regime, the longitudinal field lifts the degeneracy and selects
one of the polarized ground states. Besides, it induces a linear confining
potential between the kinks leading to bound states. As a result, the broad
continuum of free kink excitations as probed in linear response fragments into
sharp peaks [33]. At a low transverse field, the splitting can be understood
via a Schrödinger equation for the relative kink separation with a linear
potential [34] which can also be generalized to include lattice effects [35,
36, 37, 38, 39, 40, 41]. The main questions of our work are: How to
efficiently simulate the nonlinear response of the TFIM with a longitudinal
field using MPS methods? What are the robust signatures of confined bound
states in nonlinear 2D spectroscopy?

2D spectroscopy. Here, we introduce a two-pulse protocol which following
previous work Ref. 15. In this setup, two Dirac-delta pulses $B_{0}$ and
$B_{\tau}$ which are polarized along $\hat{\alpha}$ and $\hat{\beta}$
directions, respectively, reach the sample at time $T=0$ and $T=\tau>0$
successively. These magnetic pulses couple to the local moments of the sample
and the induced magnetization along $\hat{\gamma}$ direction is recorded as
$M^{\gamma}_{\mathrm{0\tau}}(T)$ at time $T=\tau+t$ where $t>0$ is the time
interval between the second pulse $B_{\tau}$ and the measurement. To subtract
the signal from the linear response, two different experiments are repeated
but with pulse $B_{0}$ or $B_{\tau}$ alone to measure $M^{\gamma}_{0}(T)$ and
$M^{\gamma}_{\tau}(T)$. The nonlinear signal field emerging from the sample at
$T=\tau+t$ in the $\hat{\gamma}$ direction is defined as

$\displaystyle
M^{\gamma}_{NL}(T)\\!\equiv\\!M^{\gamma}_{\mathrm{0\tau}}(T)-M^{\gamma}_{0}(T)-M^{\gamma}_{\mathrm{\tau}}(T).$
(2)

The nonlinear signal depends only on the nonlinear responses and directly
measures the second and higher order magnetic susceptibilities [18]:

$\displaystyle M^{\gamma}_{\mathrm{NL}}(t,\tau)$ $\displaystyle=$
$\displaystyle B_{0}B_{\tau}\chi_{\gamma\beta\alpha}^{(2)}(t,\tau+t)$ (3)
$\displaystyle+$
$\displaystyle(B_{0})^{2}B_{\tau}\chi_{\gamma\beta\alpha\alpha}^{(3)}(t,\tau+t,\tau+t)$
$\displaystyle+$ $\displaystyle
B_{0}(B_{\tau})^{2}\chi_{\gamma\beta\beta\alpha}^{(3)}(t,t,\tau+t)+O(B^{4}).$

The 2D spectrum is the Fourier transform of $M^{\gamma}_{\mathrm{NL}}(t,\tau)$
over both time domains $t$ and $\tau$. In Eq. (3), the leading contribution to
the nonlinear response in the two-pulse setup, i.e., the second order
nonlinear susceptibility $\chi_{\gamma\beta\alpha}^{(2)}(t,\tau+t)$, is given
as

$\displaystyle\chi^{(2)}_{\gamma\beta\alpha}(t,\tau+t)=$ $\displaystyle-$
$\displaystyle\frac{\theta(t)\theta(\tau)}{4L}\mathrm{Re}\sum_{j,l,m}\big{[}S^{\gamma\beta\alpha}_{j,l,m}(\tau+t,\tau,0)$
(4) $\displaystyle-$ $\displaystyle
S^{\beta\gamma\alpha}_{j,l,m}(\tau,\tau+t,0)\big{]}$

with the three point spin correlation function in the ground state
$|\psi\rangle$,

$\displaystyle
S^{\gamma\beta\alpha}_{j,l,m}(T_{1},T_{2},0)=\langle\psi|\sigma^{\gamma}_{j}(T_{1})\sigma^{\beta}_{l}(T_{2})\sigma^{\alpha}_{m}(0)|\psi\rangle$
(5)

where $\sigma^{\gamma}_{j}(T)\equiv e^{iHT}\sigma^{\gamma}_{j}e^{-iHT}$. When
the Hamiltonian and $|\psi\rangle$ preserves the lattice translation symmetry,
the site index $m$ in Eq. (4) and Eq. (5) can be fixed, e.g. as $c\equiv L/2$
the central site of the system, which we use in the following. Then,
$\chi^{(2)}_{\gamma\beta\alpha}(t,\tau+t)$ is obtained by evaluating

$\displaystyle-\frac{\theta(t)\theta(\tau)}{4}\mathrm{Re}\sum_{j,l}\big{[}S^{\gamma\beta\alpha}_{j,l,c}(\tau+t,\tau,0)-S^{\beta\gamma\alpha}_{j,l,c}(\tau,\tau+t,0)\big{]}.$
(6)

Method. A promising tool for calculating
$S^{\gamma\beta\alpha}_{j,l,c}(\tau+t,\tau,0)$ in Eq.(6) for a whole range of
site indices $j$ and $l$ is to use the iMPS method. We only need to perform
two different real time evolution runs to calculate
$\sum_{j,l}S^{\gamma\beta\alpha}_{j,l,c}(\tau+t,\tau,0)\\!=\\!\sum_{j,l}e^{iE(\tau+t)}\langle\psi|\sigma^{\gamma}_{j}e^{-iHt}\sigma^{\beta}_{l}e^{-iH\tau}\sigma^{\alpha}_{c}|\psi\rangle$
where $E$ is the ground state energy. In addition, finite size effects are
avoided. Such effect originates from the bouncing of correlations following a
local quench at site $j=1$ or $L$, boundary sites of the system. Below, we
explain a procedure to obtain
$\langle\psi|\sigma^{\gamma}_{j}e^{-iHt}\sigma^{\beta}_{l}e^{-iH\tau}\sigma^{\alpha}_{c}|\psi\rangle$
(See Section I of Supplementary Information (SI) for
$S^{\beta\gamma\alpha}_{j,l,c}(\tau,\tau+t,0)$).

Figure 1: (color online) Exploiting infinite boundary conditions for
translation-invariant systems, two real time evolution runs are sufficient to
evaluate $\sum_{j,l}S^{\gamma\beta\alpha}_{j,l,c}(\tau+t,\tau,0)$. (a) iMPS
representation of a state $e^{iHt}\sigma^{\gamma}_{c}|\psi\rangle$ with L
sites in the unit cell. (b) Transfer matrix of two distinct iMPS “bra” and
“ket” which represent $e^{iHt}\sigma^{\gamma}_{c+2}|\psi\rangle$ and
$\sigma^{\beta}_{c-2}e^{-iH\tau}\sigma^{\alpha}_{c}|\psi\rangle$ respectively.
The color gradient illustrates the light cone spreading of correlations
following a local quench.

  1. 1.

Find a ground state and perform a time evolution following a local quench,
$\sigma_{c}^{\gamma}$ or $\sigma_{c}^{\alpha}$, using infinite time evolving
block decimation (iTEBD) method[42, 43, 44] to obtain an iMPS for
$e^{iHt}\sigma^{\gamma}_{c}|\psi\rangle$, which is shown in Fig. 1(a), or
$e^{-iH\tau}\sigma^{\alpha}_{c}|\psi\rangle$.

  2. 2.

Shift every $B$ tensor of an iMPS, which represents
$e^{iHt}\sigma^{\gamma}_{c}|\psi\rangle$, $(c-j)$ sites to the right within a
window of size $L$ to obtain a new iMPS “bra” associated to
$e^{iHt}\sigma^{\gamma}_{j}|\psi\rangle$.

  3. 3.

Apply a local operator $\sigma^{\beta}_{l}$ to an iMPS associated to
$e^{-iH\tau}\sigma^{\alpha}_{c}|\psi\rangle$ and get a new iMPS “ket” which
represents $\sigma^{\beta}_{l}e^{-iH\tau}\sigma^{\alpha}_{c}|\psi\rangle$.

  4. 4.

Evaluate an overlap of two iMPS “bra” and “ket” within the window by
calculating the dominant left and right eigenvector of the corresponding
transfer matrix, which is shown in Fig. 1(b), and obtain
$\langle\psi|\sigma^{\gamma}_{j}e^{-iHt}\sigma^{\beta}_{l}e^{-iH\tau}\sigma^{\alpha}_{c}|\psi\rangle$
[33].

Figure 2: (color online) Second order susceptibility in the ferromagnetic
phase of the TFIM with $h_{x}/J=0.2$. (a) $\chi^{(2)}_{xxx}(t,\tau+t)$ from
iMPS method with a window of size $L=80$. The data is rescaled such that the
maximal absolute value is 1. (b) $\chi^{(2)}_{xxx}(t,\tau+t)$ at $Jt=15$ from
iMPS method and JW formalism. For the latter case, we set a size of the system
$L=80$ with PBC. (c) Real part of Fourier transformed
$\chi^{(2)}_{xxx}(t,\tau+t)$. Inset: TK excitation continua of the TFIM.

Before investigating the 2D spectrum of the non-integrable TFIM with
longitudinal field, we first focus on the free TFIM and compare the result of
$\chi^{(2)}_{xxx}(t,\tau+t)$ using two different schemes, i.e., the numerical
iMPS method and analytic calculations via the JW transformation with periodic
boundary condition (PBC). Here and below, we set $h_{x}/J=0.2$ and all iMPS
simulations are done with a spatial window of size $L=80$ sites and over the
time range $Jt,J\tau=30$. Within such temporal range, the light cone spreading
of correlations, which follows a local quench at the center of a spatial
window, does not reach the boundary of the window. In Fig. 2(a), we plot the
result of $\chi^{(2)}_{xxx}(t,\tau+t)$ from the iMPS method. In order to check
the errors of our method, we tracked the truncation error, the truncated
weight of many-body wave function at each time step in iTEBD, which quantifies
an upper limit for the truncation effect on local observables ($\lesssim
10^{-8}$ for every result given in our study). Besides, we also followed the
dependence of $\chi^{(2)}_{xxx}(t,\tau+t)$ on the time step $\delta t$ and the
bond dimension $\chi$, fixing to $\delta t=0.03/J$ and $\chi_{\text{max}}=30$.
In Fig. 2(b), we compare $\chi^{(2)}_{xxx}(t,\tau+t)$ at $Jt=15$ from the iMPS
method with the one from the JW formalism which confirms exact agreement. In
Fig. 2(c), we plot Re$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$, the real
part of the Fourier transformed $\chi^{(2)}_{xxx}(t,\tau+t)$. It contains a
sharp vertical line of intensity centered at $\omega_{t}=0$. Regarding
$\omega_{t}$ and $\omega_{\tau}$ as the detecting and pumping frequencies, the
response is known as a rectification signal. It also contains a diffusive,
weak non-rephasing signal in the first frequency quadrant, mirroring the
energy range of the free TK continuum, which is shown in the inset of Fig.
2(c) [15].

Results. Next, we focus on the 2D spectrum of the 1D TFIM with longitudinal
field using the iMPS method. Fig. 3(a) shows the
Re$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$, which is calculated in weakly
confined regime with $h_{z}/J=0.03$ (See Section II of SI for the
Im$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$, which is related to the
Re$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$ by the dispersion relation[45,
46, 47]). This value is similar to the one used in Ref. 35 to describe
$\mathrm{CoNb}_{2}\mathrm{O}_{6}$. New spectroscopic signals are encoded in
$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$ in the presence of a longitudinal
field. First, it contains a dominant non-rephasing signal which appears as
diagonal peaks in the first quadrant. At the same time, a weakly diffusive
terahertz rectification signal is also detected as a streak along the
$\omega_{\tau}$ axis. In Fig. 4(a), we plot
Re$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$ in strongly confined regime
with $h_{z}/J=0.4$. Unlike the previous regime, it contains a non-rephasing
like signal which appears as strong cross (off-diagonal) peaks in the first
quadrant. Besides, a subdominant rephasing like signal appears as cross peaks
in the fourth quadrant.

Figure 3: (color online) The 2D spectra
$\chi^{(2)}_{xxx}(\omega_{t},\omega_{\tau})$ of TFIM with longitudinal field
for $h_{x}/J=0.2$ and $h_{z}/J=0.03$. For the spectra, first and fourth
quadrant are only shown. The other half are obtained by complex
conjugation.(a) Result of iMPS method. (b) Result of perturbative calculation
within projected TK subspace. Inset : The energy spectrum of TK bound states
in the weakly confined regime. Figure 4: (color online) 2D spectra
$\chi^{(2)}_{xxx}(\omega_{t},\omega_{\tau})$ of TFIM with longitudinal field
for $h_{x}/J=0.2$ and $h_{z}/J=0.4$ (a) Result of iMPS method. (b) Result of
perturbative method. Inset: The energy spectrum of TK bound states in the
strongly confined regime.

To interpret the 2D spectra, we use a projected TK model and can study the
excitations of Eq.(1) perturbatively [37, 40]. The idea is to project the full
Hilbert space down to the Hilbert space of TK states, where regions of
opposite magnetization are separated by the two different domain-walls. In
this model, each TK state is represented as
$|j,l\rangle\equiv|...\uparrow\uparrow\downarrow_{j}\downarrow...\downarrow\downarrow_{(j+l-1)}\uparrow\uparrow...\rangle$,
where $j$ is the starting site of down spins. The TK model has been adopted to
phenomenologically understand the confinement of excitations observed in the
dynamical (1D) neutron response of $\mathrm{CoNb}_{2}\mathrm{O}_{6}$ [35, 37].
The projected model is expected to capture the low-energy excitations of the
original model well as long as $h_{x}/J\ll h_{z}/J$. This can be understood by
looking into the energy gap between TK states and four-kink states [48]. The
TK Hamiltonian, $\mathcal{H}_{TK}\equiv\mathcal{P}H\mathcal{P}$ with the
projector $\mathcal{P}$, acts as follows:

$\displaystyle\mathcal{H}_{TK}|j,l\rangle$ $\displaystyle=$ $\displaystyle
4J|j,l\rangle-h_{x}\big{[}|j,l+1\rangle+|j,l-1\rangle$ $\displaystyle+$
$\displaystyle|j+1,l-1\rangle+|j-1,l+1\rangle\big{]}+2h_{z}l|j,l\rangle.$

For our translational invariant model, the total momentum $p$ of the bound
state is a good quantum number. In the momentum basis
$|p,l\rangle=\sum_{j}\exp(ipj)|j,l\rangle$, the Hamiltonian is diagonal in $p$
and acts on $|p,l\rangle$ as

$\displaystyle\mathcal{H}_{TK}|p,l\rangle$ $\displaystyle=$ $\displaystyle
4J|p,l\rangle-h_{x}\big{[}(1+e^{\text{i}p})|p,l+1\rangle$ (8) $\displaystyle+$
$\displaystyle(1+e^{-\text{i}p})|p,l-1\rangle\big{]}+2h_{z}l|p,l\rangle.$

By solving the eigen-equation in Eq. (8) for a given momenta $p$ and a band
index $n$, one get TK bound state $|\Phi_{n}(p)\rangle$ with the excitation
energy $E_{n}(p)$, which is shown in the inset of Fig. 3(b) and 4(b). (See
Section III of SI for details). Then, following Ref. 37, which calculates the
linear response, we find (after some algebraic manipulation)

$\displaystyle\chi_{xxx}^{(2)}(t,\tau+t)$ $\displaystyle=$
$\displaystyle\theta(t)\theta(\tau)(Q_{1}+Q_{2}+Q_{3}+Q_{4})$ (9)

with

$\displaystyle Q_{1}$ $\displaystyle=$
$\displaystyle\sum_{n,p}C_{n,n}(p)\cos\big{[}E_{n}(p)(t+\tau)\big{]},$
$\displaystyle Q_{2}$ $\displaystyle=$
$\displaystyle\sum_{n,p}C_{n,n}(p)\cos\big{[}E_{n}(p)\tau\big{]},$
$\displaystyle Q_{3}$ $\displaystyle=$
$\displaystyle\sum_{n,m,p}C_{n,m}(p)\cos\big{[}E_{n}(p)t+E_{m}(p)\tau\big{]},$
$\displaystyle Q_{4}$ $\displaystyle=$
$\displaystyle\sum_{n,m,p}C_{n,m}(p)\cos\big{[}\big{(}E_{n}(p)-E_{m}(p)\big{)}t-E_{m}(p)\tau\big{]}$

where $C_{n,m}(p)$ is the optical matrix element which depends on band indices
$n,m$ and momentum $p$ (See Section III of SI for details).

The interpretation of 2D spectra now becomes transparent: $Q_{1}$ gives rise
to diagonal non-rephasing peaks at $\omega_{t},\omega_{\tau}=E_{n}(p)$ and
$Q_{2}$ produces dominant terahertz rectification signals at $\omega_{t}=0$
and $\omega_{\tau}=\pm E_{m}(p)$ [Fig. 3(b) and Fig. 4(b)]. In the strongly
confined regime, $Q_{3}$ gives rise to dominant cross peaks in the first
frequency quadrant originating from the non-rephasing like process [Fig.
4(b)]. To be more precise, such peaks sharply appear at
$\omega_{t}\\!=\\!E_{1}(0)$ and $\omega_{\tau}\\!=\\!E_{2}(0)$ or
$\omega_{t}\\!=\\!E_{2}(0)$ and $\omega_{\tau}\\!=\\!E_{1}(0)$, indicating the
presence of multiple TK excited states. Such sharp peaks do not appear in the
free 1D TFIM, which can be mapped to an independent two-level systems with
each having a single excited state [15]. $Q_{4}$ contains terms which induce
sharp (non-)rephasing like signals in the first (fourth) frequency quadrant
which are visible in the strongly confined regime [Fig. 4(b)]. Such signals
also originate from the presence of multiple excited states and appear at
$\omega_{t}=E_{2}(0)-E_{1}(0)$, an energy gap between first and second excited
states, and $\omega_{\tau}=E_{2}(0)$ or $-E_{1}(0)$, see inset.

Conclusions. In the present work, we have developed an iMPS method for
calculating the nonlinear response of 1D spin systems. As a demonstration, we
calculated the second order susceptibility, which dominates the nonlinear
response, for the ferromagnetic 1D TFIM where a single spin flip is
fractionalized into two freely moving domain-walls. We benchmarked our
numerical results with exact analytical calculations. We then included a
longitudinal field, which induces a linear confining potential between kink
excitations. In the presence of a longitudial field, the second order
susceptibility contains new signals which give rise to strong non-rephasing
and rephasing like peaks. To understand the emergence of such signals, we
employ a simplified two-kink description, which describes the low-energy
excitations in the strongly confined regime, and calculate the second order
susceptibility perturbatively. The approximate method captures the nonlinear
response of the system in the strongly confined regime and allows for a simple
interpretation.

As a future direction, it would be interesting to apply our iMPS method near
the quantum critical point between the ferromagnetic and paramagnetic states
where a hidden E8 symmetry emerges [35, 33]. A crucial question regards then
the existence of robust signals in 2D spectrum which detect the emergent
symmetry. Regarding the microscopic description and 3D nature of
$\textrm{CoNb}_{2}\textrm{O}_{6}$ [35, 36]
($\textrm{BaCo}_{2}\textrm{V}_{2}\textrm{O}_{8}$ [49, 50]) it would be
interesting to study a more quantitative model [28, 29] and go beyond simple
chains, e.g. by extending these to coupled Ising ladders [36]. Terahertz 2D
coherent spectroscopy holds the promise of uncovering the nature of exotic
excitations in strongly correlated quantum materials. A challenging but very
worthwhile direction will be an extension of our method to quantum magnets
beyond one-dimension where the nature of fractionalized excitations and
confinement thereof remains poorly understood.

Note added. When finalizing the manuscript, related works appeared that
investigate nonlinear response from quasiparticle interactions [51, 52].

##  I Acknowledgments

We thank N. P. Armitage, R. Coldea, M. Drescher, H.-K. Jin, and W. Choi for
insightful discussions related to this work. G.B.S. is funded by the European
Research Council (ERC) under the European Unions Horizon 2020 research and
innovation program (grant agreement No. 771537). F.P. acknowledges the support
of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy EXC-2111-390814868. J. K. acknowledges support
from the Imperial-TUM flagship partnership. The research is part of the Munich
Quantum Valley, which is supported by the Bavarian state government with funds
from the Hightech Agenda Bayern Plus. Tensor network calculations were
performed using the TeNPy Library [53].

## References

  * Devereaux and Hackl [2007] T. P. Devereaux and R. Hackl, Reviews of modern physics 79, 175 (2007).
  * Shen [1984] Y.-R. Shen,  _Principles of nonlinear optics_ (Wiley-Interscience, New York, NY, USA, 1984).
  * Hamm and Zanni [2011] P. Hamm and M. Zanni,  _Concepts and methods of 2D infrared spectroscopy_ (Cambridge University Press, 2011).
  * Mukamel [1999] S. Mukamel,  _Principles of nonlinear optical spectroscopy_ , 6 (Oxford University Press on Demand, 1999).
  * Fiebig  _et al._ [2005] M. Fiebig, V. V. Pavlov, and R. V. Pisarev, JOSA B 22, 96 (2005).
  * Chu  _et al._ [2020] H. Chu, M.-J. Kim, K. Katsumi, S. Kovalev, R. D. Dawson, L. Schwarz, N. Yoshikawa, G. Kim, D. Putzky, Z. Z. Li, _et al._ , Nature communications 11, 1 (2020).
  * Schwarz  _et al._ [2020] L. Schwarz, B. Fauseweh, N. Tsuji, N. Cheng, N. Bittner, H. Krull, M. Berciu, G. Uhrig, A. Schnyder, S. Kaiser, _et al._ , Nature communications 11, 1 (2020).
  * Wu  _et al._ [2017] L. Wu, S. Patankar, T. Morimoto, N. L. Nair, E. Thewalt, A. Little, J. G. Analytis, J. E. Moore, and J. Orenstein, Nature Physics 13, 350 (2017).
  * Shao  _et al._ [2021] Y. Shao, R. Jing, S. H. Chae, C. Wang, Z. Sun, E. Emmanouilidou, S. Xu, D. Halbertal, B. Li, A. Rajendran, _et al._ , Proceedings of the National Academy of Sciences 118 (2021).
  * He  _et al._ [2021] P. He, H. Isobe, D. Zhu, C.-H. Hsu, L. Fu, and H. Yang, Nature communications 12, 1 (2021).
  * Zhao  _et al._ [2017] L. Zhao, C. Belvin, R. Liang, D. Bonn, W. Hardy, N. Armitage, and D. Hsieh, Nature Physics 13, 250 (2017).
  * Zhao  _et al._ [2016] L. Zhao, D. Torchinsky, H. Chu, V. Ivanov, R. Lifshitz, R. Flint, T. Qi, G. Cao, and D. Hsieh, Nature Physics 12, 32 (2016).
  * Mahmood  _et al._ [2021] F. Mahmood, D. Chaudhuri, S. Gopalakrishnan, R. Nandkishore, and N. Armitage, Nature Physics 17, 627 (2021).
  * Lu  _et al._ [2016] J. Lu, Y. Zhang, H. Y. Hwang, B. K. Ofori-Okai, S. Fleischer, and K. A. Nelson, Proceedings of the National Academy of Sciences 113, 11800 (2016).
  * Wan and Armitage [2019] Y. Wan and N. Armitage, Physical Review Letters 122, 257401 (2019).
  * Li  _et al._ [2021] Z.-L. Li, M. Oshikawa, and Y. Wan, Physical Review X 11, 031035 (2021).
  * Choi  _et al._ [2020] W. Choi, K. H. Lee, and Y. B. Kim, Physical Review Letters 124, 117205 (2020).
  * Nandkishore  _et al._ [2021] R. M. Nandkishore, W. Choi, and Y. B. Kim, Physical Review Research 3, 013254 (2021).
  * Fava  _et al._ [2021] M. Fava, S. Biswas, S. Gopalakrishnan, R. Vasseur, and S. Parameswaran, Proceedings of the National Academy of Sciences 118 (2021).
  * Woutersen  _et al._ [2002] S. Woutersen, R. Pfister, P. Hamm, Y. Mu, D. S. Kosov, and G. Stock, The Journal of chemical physics 117, 6833 (2002).
  * Cho [2008] M. Cho, Chemical reviews 108, 1331 (2008).
  * Terranova and Corcelli [2014] Z. Terranova and S. Corcelli, The Journal of Physical Chemistry B 118, 8264 (2014).
  * Paeckel  _et al._ [2019] S. Paeckel, T. Köhler, A. Swoboda, S. R. Manmana, U. Schollwöck, and C. Hubig, Annals of Physics 411, 167998 (2019).
  * Vanderstraeten  _et al._ [2015] L. Vanderstraeten, F. Verstraete, and J. Haegeman, Physical Review B 92, 125136 (2015).
  * Bera  _et al._ [2017] A. Bera, B. Lake, F. Essler, L. Vanderstraeten, C. Hubig, U. Schollwöck, A. Islam, A. Schneidewind, and D. Quintero-Castro, Physical Review B 96, 054423 (2017).
  * Vanderstraeten  _et al._ [2018] L. Vanderstraeten, M. Van Damme, H. P. Büchler, and F. Verstraete, Physical Review Letters 121, 090603 (2018).
  * Van Damme  _et al._ [2021] M. Van Damme, R. Vanhove, J. Haegeman, F. Verstraete, and L. Vanderstraeten, Physical Review B 104, 115142 (2021).
  * Fava  _et al._ [2020] M. Fava, R. Coldea, and S. Parameswaran, Proceedings of the National Academy of Sciences 117, 25219 (2020).
  * Morris  _et al._ [2021] C. Morris, N. Desai, J. Viirok, D. Hüvonen, U. Nagel, T. Room, J. Krizan, R. Cava, T. McQueen, S. Koohpayeh, _et al._ , Nature Physics , 1 (2021).
  * Lee  _et al._ [2010] S. Lee, R. K. Kaul, and L. Balents, Nature Physics 6, 702 (2010).
  * Kinross  _et al._ [2014] A. Kinross, M. Fu, T. Munsie, H. Dabkowska, G. Luke, S. Sachdev, and T. Imai, Physical Review X 4, 031008 (2014).
  * Xu  _et al._ [2022] Y. Xu, L. Wang, Y. Huang, J. Ni, C. Zhao, Y. Dai, B. Pan, X. Hong, P. Chauhan, S. Koohpayeh, _et al._ , Physical Review X 12, 021020 (2022).
  * Kjäll  _et al._ [2011] J. A. Kjäll, F. Pollmann, and J. E. Moore, Physical Review B 83, 020407 (2011).
  * McCoy and Wu [1978] B. M. McCoy and T. T. Wu, Physical Review D 18, 1259 (1978).
  * Coldea  _et al._ [2010] R. Coldea, D. Tennant, E. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl, and K. Kiefer, Science 327, 177 (2010).
  * Morris  _et al._ [2014] C. Morris, R. V. Aguilar, A. Ghosh, S. Koohpayeh, J. Krizan, R. Cava, O. Tchernyshyov, T. McQueen, and N. Armitage, Physical review letters 112, 137403 (2014).
  * Rutkevich [2010] S. Rutkevich, Journal of Statistical Mechanics: Theory and Experiment 2010, P07015 (2010).
  * Rutkevich [2008] S. Rutkevich, Journal of Statistical Physics 131, 917 (2008).
  * Shinkevich and Syljuåsen [2012] S. Shinkevich and O. F. Syljuåsen, Physical Review B 85, 104408 (2012).
  * Kormos  _et al._ [2017] M. Kormos, M. Collura, G. Takács, and P. Calabrese, Nature Physics 13, 246 (2017).
  * Liu  _et al._ [2019] F. Liu, R. Lundgren, P. Titum, G. Pagano, J. Zhang, C. Monroe, and A. V. Gorshkov, Physical review letters 122, 150601 (2019).
  * Vidal [2003] G. Vidal, Physical review letters 91, 147902 (2003).
  * Vidal [2004] G. Vidal, Physical review letters 93, 040502 (2004).
  * Vidal [2007] G. Vidal, Physical review letters 98, 070201 (2007).
  * Kogan [1963] S. M. Kogan, Sov. Phys. JETP 16, 217 (1963).
  * Caspers [1964] W. J. Caspers, Physical Review 133, A1249 (1964).
  * Bassani and Scandolo [1991] F. Bassani and S. Scandolo, Physical Review B 44, 8446 (1991).
  * Verdel  _et al._ [2020] R. Verdel, F. Liu, S. Whitsitt, A. V. Gorshkov, and M. Heyl, Physical Review B 102, 014308 (2020).
  * Zou  _et al._ [2021] H. Zou, Y. Cui, X. Wang, Z. Zhang, J. Yang, G. Xu, A. Okutani, M. Hagiwara, M. Matsuda, G. Wang, _et al._ , Physical review letters 127, 077201 (2021).
  * Faure  _et al._ [2018] Q. Faure, S. Takayoshi, S. Petit, V. Simonet, S. Raymond, L.-P. Regnault, M. Boehm, J. S. White, M. Månsson, C. Rüegg, _et al._ , Nature Physics 14, 716 (2018).
  * Fava  _et al._ [2022] M. Fava, S. Gopalakrishnan, R. Vasseur, F. H. Essler, and S. Parameswaran, arXiv preprint arXiv:2208.09490 (2022).
  * Hart and Nandkishore [2022] O. Hart and R. Nandkishore, arXiv preprint arXiv:2208.12817 (2022).
  * Hauschild and Pollmann [2018] J. Hauschild and F. Pollmann, SciPost Physics Lecture Notes , 005 (2018).
  * McCulloch [2008] I. P. McCulloch, arXiv preprint arXiv:0804.2509 (2008).
  * Schollwöck [2011] U. Schollwöck, Annals of physics 326, 96 (2011).
  * Phien  _et al._ [2012] H. N. Phien, G. Vidal, and I. P. McCulloch, Physical Review B 86, 245107 (2012).
  * Milsted  _et al._ [2013] A. Milsted, J. Haegeman, T. J. Osborne, and F. Verstraete, Physical Review B 88, 155116 (2013).
  * Binder and Barthel [2018] M. Binder and T. Barthel, Physical Review B 98, 235114 (2018).
  * Ejima  _et al._ [2021] S. Ejima, F. Lange, and H. Fehske, SciPost Physics 10, 077 (2021).

##  II Supplementary Material

###  II.1 A procedure to obtain $S^{\beta\gamma\alpha}_{j,l,c}(\tau,\tau+t,0)$
with iMPS method

In this section, we provide steps to calculate
$S^{\beta\gamma\alpha}_{j,l,c}(\tau,\tau+t,0)=e^{iE\tau}\langle\psi|\sigma^{\beta}_{j}e^{iHt}\sigma^{\gamma}_{l}e^{-iH(\tau+t)}\sigma^{\alpha}_{c}|\psi\rangle$
using iMPS.

  1. 1.

Find an iMPS approximation of the ground state $|\psi\rangle$ with energy $E$
[54, 55].

  2. 2.

Allow the tensors of the iMPS with a spatial window of size $L$ to vary in
time as in Refs. 56, 57, 58, 59.

  3. 3.

Apply a local operator $\sigma^{\beta}_{c}$ ($\sigma^{\alpha}_{c}$) at the
center of the window to get $\sigma^{\beta}_{c}|\psi\rangle$
($\sigma^{\alpha}_{c}|\psi\rangle$).

  4. 4.

Perform a real time evolution following the local quench $\sigma^{\beta}_{c}$
($\sigma^{\alpha}_{c}$) using iTEBD method [42, 43, 44] to obtain an iMPS
which represents $e^{-iHt}\sigma^{\beta}_{c}|\psi\rangle$
($e^{-iH(\tau+t)}\sigma^{\alpha}_{c}|\psi\rangle$).

  5. 5.

Shift every $B$ tensor of an iMPS, which represents
$e^{-iHt}\sigma^{\beta}_{c}|\psi\rangle$, $(c-j)$ sites to the right within
the window to obtain a new iMPS “bra” associated to
$e^{-iHt}\sigma^{\beta}_{j}|\psi\rangle$.

  6. 6.

Apply an operator $\sigma^{\gamma}_{l}$ to an iMPS, which approximates
$e^{-iH(\tau+t)}\sigma^{\alpha}_{c}|\psi\rangle$, and get a new iMPS “ket”
associated to
$\sigma^{\gamma}_{l}e^{-iH(\tau+t)}\sigma^{\alpha}_{c}|\psi\rangle$.

  7. 7.

Evaluate an overlap of two iMPS “bra” and “ket” within the window by
calculating the dominant eigenvalue of the corresponding transfer matrix to
obtain
$\langle\psi|\sigma^{\beta}_{j}e^{iHt}\sigma^{\gamma}_{l}e^{-iH(\tau+t)}\sigma^{\alpha}_{c}|\psi\rangle$.

  8. 8.

Multiply $e^{iE\tau}$ and
$\langle\psi|\sigma^{\beta}_{j}e^{iHt}\sigma^{\gamma}_{l}e^{-iH(\tau+t)}\sigma^{\alpha}_{c}|\psi\rangle$.

###  II.2 Imaginary part of $\chi^{(2)}_{xxx}(\omega_{t},\omega_{\tau})$

In this section, we focus on Im$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$ in
the ferromagnetic phase of the Ising model. In Fig. S1(a), we plot
Im$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$ for $h_{x}/J=0.2$ without
longitudinal field using iMPS method. To obtain the result, we directly
perform fast Fourier transformations of the real time correlation functions
$\chi^{(2)}_{xxx}(t,\tau+t)$. Fig. S1(b) shows
Im$\chi^{(2)}_{xxx}(\omega_{t},\omega_{\tau})$ in weakly confined regime with
$h_{z}/J=0.03$. One can clearly observe the emergence of strong non-rephasing
like signals in the first frequency quadrant. In Fig. S1(c), we plot
Im$\chi^{(2)}_{xxx}(\omega_{t},\omega_{\tau})$ in strongly confined regime
with $h_{z}/J=0.4$. Simliar to the real part, it contains new signals which
include dominant non-rephasing and strong rephasing like peaks.

(a) Im$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$ without longitudial field

(b) Im$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$ in weakly confined regime

(c) Im$\chi_{xxx}^{(2)}(\omega_{t},\omega_{\tau})$ in strongly confined regime

Figure S1: (color online)

###  II.3 Three point spin correlation functions within TK model

In this section, we follow Ref. 37, which calculates the two point spin
correlator, and formulate the three point correlator within the projected TK
low energy subspace. In this subspace, each TK state is represented as
$|j,l\rangle\equiv|...\uparrow\uparrow\downarrow_{j}\downarrow...\downarrow\downarrow_{(j+l-1)}\uparrow\uparrow...\rangle$,
where $j$ is the starting site of down spins. The TK Hamiltonian
$\mathcal{H}_{TK}$ acts as follows:

$\displaystyle\mathcal{H}_{TK}|j,l\rangle$ $\displaystyle=$ $\displaystyle
4J|j,l\rangle-
h_{x}\big{[}|j,l+1\rangle+|j,l-1\rangle+|j+1,l-1\rangle+|j-1,l+1\rangle\big{]}+2h_{z}l|j,l\rangle.$
(S1)

In the momentum basis $|p,l\rangle=\sum_{j}\exp(ipj)|j,l\rangle$, the
Hamiltonian is diagonal in $p$ and acts on $|p,l\rangle$ as

$\displaystyle\mathcal{H}_{TK}|p,l\rangle$ $\displaystyle=$ $\displaystyle
4J|p,l\rangle-
h_{x}\big{[}(1+e^{\text{i}p})|p,l+1\rangle+(1+e^{-\text{i}p})|p,l-1\rangle\big{]}+2h_{z}l|p,l\rangle.$
(S2)

Now, the eigen-equation with the excitation energy $E_{n}(p)$,

$\displaystyle\mathcal{H}_{TK}|\Phi_{n}(p)\rangle=E_{n}(p)|\Phi_{n}(p)\rangle$
(S3)

where
$|\Phi_{n}(p)\rangle\equiv\sum_{l=1}^{\infty}\exp(\frac{\text{i}pl}{2})\psi_{n}(l,p)|p,l\rangle/\sqrt{\sum_{l=1}^{\infty}|\psi_{n}(l,p)|^{2}}$
with the site integer $l$, momentum $p$, and discrete band index $n$ takes the
form as below.

$\displaystyle[4J+2h_{z}l-E_{n}(p)]\psi_{n}(l,p)-2h_{x}\cos(p/2)[\psi_{n}(l+1,p)+\psi_{n}(l-1,p)]=0$
(S4)

We rewrite Eq. (S4) with dimensionless parameters as follows

$\displaystyle\left(-\lambda_{n}+\mu\,l\right)\psi_{n}(l,p)-\frac{\psi_{n}(l+1,p)+\psi_{n}(l-1,p)}{2}=0$
(S5)

where

$\displaystyle\lambda_{n}\equiv\frac{E_{n}(p)-4J}{4h_{x}\cos(p/2)},\quad\mu\equiv\frac{h_{z}}{2h_{x}\cos(p/2)}$
(S6)

with the boundary conditions, $\lim_{l\to 0}\psi_{n}(l,p)=0$ and
$\lim_{l\to+\infty}\psi_{n}(l,p)=0$ [37]. The result for the eigenvalues
$\lambda_{n}$ reads as

$\displaystyle\lambda_{n}=-\mu\,\nu_{n}$ (S7)

where $\nu_{n}$ are the solutions of the equation

$\displaystyle J_{\nu_{n}}(1/\mu)=0.$ (S8)

Here, $J_{\nu}(x)$ is the Bessel function of order $\nu$. Then, the solution
of the Eq. (S4) reads as

$\displaystyle E_{n}(p)=4J-2h_{z}\,\nu_{n}.$ (S9)

As given in the main text, $\chi^{(2)}_{xxx}(t,\tau+t)$ is written as

$\displaystyle\chi^{(2)}_{xxx}(t,\tau+t)=$ $\displaystyle-$
$\displaystyle\frac{\theta(t)\theta(\tau)}{4L}\mathrm{Re}\sum_{j,k,l}\big{[}S^{xxx}_{j,k,l}(\tau+t,\tau,0)-S^{xxx}_{j,k,l}(\tau,\tau+t,0)\big{]}$
(S10)

where
$S^{xxx}_{j,k,l}(s_{1},s_{2},s_{3})\equiv\langle\sigma^{x}_{j}(s_{1})\sigma^{x}_{k}(s_{2})\sigma^{x}_{l}(s_{3})\rangle$
is the three point correlation function. The first term of Eq. (S10) can be
rewritten as follows :

$\displaystyle\sum_{j,k,l}S^{xxx}_{j,k,l}(\tau+t,\tau,0)$ $\displaystyle=$
$\displaystyle\sum_{j,k,l}\langle\psi|e^{\text{i}H(\tau+t)}\sigma^{x}_{j}e^{-\text{i}Ht}\sigma^{x}_{k}e^{-\text{i}H\tau}\sigma^{x}_{l}|\psi\rangle$
(S11) $\displaystyle=$
$\displaystyle\sum_{j,k,l}\sum_{n,m,p,q}e^{-\text{i}(E_{n}(p)t+E_{m}(q)\tau)}\langle\psi|\sigma^{x}_{j}|\Phi_{n}(p)\rangle\langle\Phi_{n}(p)|\sigma^{x}_{k}|\Phi_{m}(q)\rangle\langle\Phi_{m}(q)|\sigma^{x}_{l}|\psi\rangle$
$\displaystyle=$
$\displaystyle\sum_{j,k,l}\sum_{n,m,p}e^{-\text{i}(E_{n}(p)t+E_{m}(p)\tau)}\langle\psi|\sigma^{x}_{j}|\Phi_{n}(p)\rangle\langle\Phi_{n}(p)|\sigma^{x}_{k}|\Phi_{m}(p)\rangle\langle\Phi_{m}(p)|\sigma^{x}_{l}|\psi\rangle$

where $|\psi\rangle$ is the non-degenerate ferromagnetic ground state with no
kinks. In Eq. (S11), $\langle\Phi_{n}(p)|\sigma^{x}_{k}|\Phi_{m}(p)\rangle$
can be expanded as

$\displaystyle\sum_{k}\langle\Phi_{n}(p)|\sigma^{x}_{k}|\Phi_{m}(p)\rangle=\frac{\sum_{k}\sum_{l^{\prime},l}\langle
p,l^{\prime}|\sigma^{x}_{k}|p,l\rangle
e^{-\text{i}p(l^{\prime}-l)/2}\psi_{n}(l^{\prime},p)\psi_{m}(l,p)}{\sqrt{\sum_{l=1}^{\infty}|\psi_{n}(l,p)|^{2}}\sqrt{\sum_{l=1}^{\infty}|\psi_{m}(l,p)|^{2}}}$
(S12)

with

$\displaystyle\langle
p,l^{\prime}|\sigma^{x}_{k}|p,l\rangle=\sum_{j^{\prime},j}\exp(-\text{i}p(j^{\prime}-j))\langle
j^{\prime},l^{\prime}|\sigma^{x}_{k}|j,l\rangle.$ (S13)

Eq. (S13) clearly shows that Eq. (S12) is finite only when one of the
following four conditions is met.

$\displaystyle 1.~{}|j,l\rangle=|k+1,l\rangle,~{}~{}$
$\displaystyle|j^{\prime},l^{\prime}\rangle=|k,l+1\rangle$ $\displaystyle
2.~{}|j,l\rangle=|k-l,l\rangle,~{}~{}$
$\displaystyle|j^{\prime},l^{\prime}\rangle=|k-l,l+1\rangle,$ $\displaystyle
3.~{}|j,l\rangle=|k,l\rangle,~{}~{}$
$\displaystyle|j^{\prime},l^{\prime}\rangle=|k+1,l-1\rangle,$ $\displaystyle
4.~{}|j,l\rangle=|k-l+1,l\rangle,~{}~{}$
$\displaystyle|j^{\prime},l^{\prime}\rangle=|k-l+1,l-1\rangle$

Then, Eq. (S12) can be simplified as the following.

$\displaystyle\sum_{k}\langle\Phi_{n}(p)|\sigma^{x}_{k}|\Phi_{m}(p)\rangle=\frac{\sum_{l}\cos{\frac{p}{2}}\big{[}\psi_{n}(l,p)\psi_{m}(l+1,p)+\psi_{n}(l,p)\psi_{m}(l-1,p)\big{]}}{\sqrt{\sum_{l=1}^{\infty}|\psi_{n}(l,p)|^{2}}\sqrt{\sum_{l=1}^{\infty}|\psi_{m}(l,p)|^{2}}}$
(S14)

To formulate Eq. (S14), we first rewrite Eq. (S5) with two different band
indices, $n$ and $m$.

$\displaystyle(-\lambda_{n}+\mu\,l)\psi_{n}(l,p)=\frac{\psi_{n}(l+1,p)+\psi_{n}(l-1,p)}{2}$
(S15)
$\displaystyle(-\lambda_{m}+\mu\,l)\psi_{m}(l,p)=\frac{\psi_{m}(l+1,p)+\psi_{m}(l-1,p)}{2}$
(S16)

Now, we multiply $\psi_{m}(l+1,p)$ to Eq. (S15) and $\psi_{n}(l-1,p)$ to Eq.
(S16), subtract one from the other, and sum it over $n$. Then, we get an
expression for $\sum_{l}\psi_{n}(l,p)\psi_{m}(l+1,p)$ :

$\displaystyle\sum_{l}\psi_{n}(l,p)\psi_{m}(l+1,p)=\frac{\psi_{n}(0,p)\psi_{m}(2,p)-\psi_{n}(1,p)\psi_{m}(1,p)}{2(\lambda_{n}-\lambda_{m}-\mu)}=\frac{-2}{\lambda_{n}-\lambda_{m}-\mu}$
(S17)

where we set $\psi_{n}(1,p)=-2$ without loss of generality [37]. We can
proceed similarly and get

$\displaystyle\sum_{l}\psi_{n}(l,p)\psi_{m}(l-1,p)=\frac{-2}{\lambda_{m}-\lambda_{n}-\mu}$
(S18)

In the end, we get expressions for $\sum_{j,k,l}S^{xxx}_{jkl}(\tau+t,\tau,0)$
as the following.

$\displaystyle\sum_{j,k,l}S^{xxx}_{jkl}(\tau+t,\tau,0)=\sum_{n,m,p}C_{n,m}(p)e^{-\text{i}(E_{n}(p)t+E_{m}(p)\tau)}$
(S19)

with the optical matrix element

$\displaystyle
C_{n,m}(p)\equiv-2\cos{\frac{p}{2}}(\frac{1}{\lambda_{m}-\lambda_{n}-\mu}+\frac{1}{\lambda_{n}-\lambda_{m}-\mu})I_{n}(p)I_{m}(p).$
(S20)

Here, the relative intensity of the $n$-th mode is defined as

$\displaystyle
I_{n}(p)\equiv\frac{|\psi_{n}(1,p)|^{2}}{\sum_{l=1}^{\infty}|\psi_{n}(l,p)|^{2}}=\frac{4}{\sum_{l=1}^{\infty}|\psi_{n}(l,p)|^{2}}=2\mu\left\\{\frac{\partial}{\partial\nu}\left[\frac{J_{\nu}(1/\mu)}{J_{\nu+1}(1/\mu)}\right]\right\\}^{-1}\bigg{|}_{\nu\to\nu_{n}}$
(S21)

where $\nu_{n}$ is the $n$-th solution of Eq. (S8) [37]. Similiar formulation
gives

$\displaystyle\sum_{j,k,l}S^{xxx}_{jkl}(\tau,\tau+t,0)=\sum_{n,m,p}C_{n,m}(p)e^{-\text{i}E_{m}(p)\tau}e^{-\text{i}(E_{m}(p)-E_{n}(p))t}.$
(S22)