# Imposing equilibrium on measured 3-D stress fields using Hodge decomposition
and FFT-based optimization

Hao Zhou  Ricardo A. Lebensohn  Péter Reischig  Wolfgang Ludwig  Kaushik
Bhattacharya

###### Abstract

We present a methodology to impose micromechanical constraints, i.e. stress
equilibrium at grain and sub-grain scale, to an arbitrary (non-equilibrated)
voxelized stress field obtained, for example, by means of synchrotron X-ray
diffraction techniques. The method consists in finding the equilibrated stress
field closest (in $L^{2}$-norm sense) to the measured non-equilibrated stress
field, via the solution of an optimization problem. The extraction of the
divergence-free (equilibrated) part of a general (non-equilibrated) field is
performed using the Hodge decomposition of a symmetric matrix field, which is
the generalization of the Helmholtz decomposition of a vector field into the
sum of an irrotational field and a solenoidal field. The combination of: a)
the Euler-Lagrange equations that solve the optimization problem, and b) the
Hodge decomposition, gives a differential expression that contains the bi-
harmonic operator and two times the curl operator acting on the measured
stress field. These high-order derivatives can be efficiently performed in
Fourier space. The method is applied to filter the non-equilibrated parts of a
synthetic piecewise constant stress fields with a known ground truth, and
stress fields in Gum Metal, a beta-Ti-based alloy measured in-situ using
Diffraction Contrast Tomography (DCT). In both cases, the largest corrections
were obtained near grain boundaries.

###### keywords:

X-ray diffraction , experimental mechanics , stress equilibrium , Hodge
decomposition , fast Fourier transform , micromechanics

††journal: Mechanics of Materials

[1]Division of Engineering and Applied Science, California Institute of
Technology, Pasadena CA 91125, USA [2]Theoretical Division, Los Alamos
National Laboratory, Los Alamos, NM 87544, USA [3]European Synchrotron
Radiation Facility, Grenoble, France [4]InnoCryst Ltd, Daresbury, Warrington,
WA4 4AD, UK [5]Laboratoire Materiaux, Ingenierie et Science, INSA Lyon, Lyon,
France

Analytical formulation to impose stress equilibrium is obtained through a
variational optimization method.

A FFT-based numerical solver is developed and further modified for interface
discontinuities.

The algorithm is first applied to synthetic stress fields, which provide a
known ground truth.

The algorithm is applied to a 3-D X-ray diffraction experiment of a beta-Ti
alloy, resulting in large corrections near grain boundaries.

##  1 Introduction

Emerging non-destructive X-ray diffraction (XRD) techniques for 3-D material
characterization performed at light sources, collectively known as 3DXRD (e.g.
[1, 2, 3]), provide previously inaccessible in-situ microstructural and
micromechanical information on polycrystalline materials at grain and sub-
grain scales. In-situ 3-D measurements are enabled by the collection of
multiple diffraction patterns as the sample is rotated and deformed. These
patterns are then used in combination with forward-modelling of the
diffraction experiment to optimize the agreement between those measured
detector images and simulated patterns corresponding to every possible
orientation and micromechanical field configuration. Near-field (nf)
techniques [4, 5, 6] are used to obtain crystal orientation fields, resulting
in space-resolved voxelized microstructural images with intragranular
resolution. Further, far-field (ff) techniques [7, 8, 9] provide information
on local micromechanical fields. Due to a trade-off between direct and
reciprocal space resolutions, the original ff data inversion methods were able
to deliver average stresses/elastic strains in single crystal grains and the
volume and location of the centers of mass of those grains, but not voxelized
intragranular fields.

Recently, 3DXRD data inversion methods were improved to provide not only
voxelized crystal orientation but also stress fields with intragranular
resolution using a variety of experimental procedures, including a microbeam
(1D), line beam (2D) or box beam (3D) illumination, each demanding a tailored
optimization-based data processing approach [10, 11, 12, 13, 14]. However,
these types of adopted optimization methods have been exclusively based on
imposing diffraction constraints, by minimizing the difference between
experimental and simulated detector images obtained by forward-modelling of
the diffraction experiment, and retaining as measured fields the orientation
and stress fields that produced that minimum. In either case, using the
original data inversion based on grain averages, or the more advanced methods
that deliver voxelized information, the resulting stress fields are in general
unbalanced (not divergence-free, violating the stress equilibrium condition).

Since the micro-mechanical field that directly affects diffraction patterns—by
changing the local lattice spacing—is elastic strain, accounting for stress
equilibrium at the sub-grain level requires a fairly accurate knowledge of the
single crystal properties, such as zero-stress lattice parameters, and single
crystal elastic constants. Thus, failing to fulfill mechanical equilibrium in
the solution may be unavoidable. In other instances, ignoring the latter
allows for a simplified and faster solution process, where large-scale
optimization is already challenging and computationally intensive. However,
when the material properties are known, taking into account mechanical
equilibrium clearly provides a more constrained solution, and may
significantly improve the stress/strain fields determination, and
consequently, the spatial resolution, reliability, and application range of
3DXRD methods. The elastic deformation solvers are typically iterative, and
thus mechanical constraints could be employed in one of many ways: to correct
the final stress field solution, in every iteration step by enforcing
equilibrium in the latest solution, or as an additional component in an
optimization target function.

Given the aforementioned problem, an approach to solve the lack of equilibrium
of the measured stress field is to combine micromechanical modelling and
simulations with experimental observations. In different contexts, this
combination has been recently done by Pokharel and Lebensohn [15] using a fast
Fourier transform (FFT)-based elasto-viscoplastic (EVPFFT) model, and by
Chatterjee et al. [16] adopting a field dislocation mechanics finite element
model, to complement X-ray diffraction based experimental information and
characterize residual and internal stresses. Also, Pagan and Beaudoin [17]
used lattice orientation and crystal plasticity kinematics to recover the
geometry, and further calculated the stress through finite element simulation.

While the above authors solve physics-based field equations, McNelis, Dawson
and Miller [18] proposed a new approach where they seek to impose equilibrium
to match the lattice stresses inferred from XRD and elasticity. In particular,
they were interested in large mechanical parts (at the order of meters) where
the detection spots are limited. Therefore, they proposed a two-scale method,
where the continuum scale stress field imposes equilibrium to match the
lattice scale stresses. The approach was extended to three dimensions by Demir
et al. [19], and has found successful applications [20, 21]. However, this
approach does not impose equilibrium at the grain and sub-grain scale, which
is the focus of our work.

This paper presents a novel methodology to impose micromechanical constraints,
i.e. stress equilibrium at the subgrain and grain scale, to an non-
equilibrated voxelized stress field obtained, for example, by means of
synchrotron X-ray diffraction techniques. The main idea is to find the
equilibrated stress field closest (in $L^{2}$-norm sense) to a measured (and
possibly) non-equilibrated stress field, via the solution of an optimization
problem based on the Hodge decomposition of a symmetric 2nd rank tensorial
field. The Hodge decomposition is a generalization of the classical Helmholtz
decomposition that states that any sufficiently smooth vector field can be
decomposed into the sum of an irrotational vector field and a solenoidal
vector field. We use the version of the Hodge decomposition for symmetric
tensorial fields developed by Geymonat and Krasucki [22]. We present an
efficient numerical implementation of the Hodge decomposition of symmetric 2nd
rank tensorial fields using FFTs, and apply it to the analysis/determination
of equilibrated stress fields in polycrystalline materials, including 3-D
space-resolved stress measured by X-ray diffraction. As mentioned before,
these stress fields are in general not divergence-free, due to piecewise
approximation based on grain averages and/or absent the consideration of any
micromechanical constraint in the data inversion procedure. Using the proposed
decomposition and formulating an optimization problem, unbalanced stress
fields can be filtered to extract their divergence-free part. The proposed
methodology has been further modified/extended to deal with interfaces and
discontinuities occurring at grain boundaries. We begin by demonstrating that
the method is capable of good recovery of synthetic data obtained in the
elastic regime with an FFT-based micromechanical model that provides an
objective ground-truth (a validation methodology originally used in [13]), and
then apply it to a stress field actually obtained from diffraction data.

The outline of the paper is as follows. In Section 2 we present a summary of
the 3D orientation and strain mapping approach developed by Reischig and
Ludwig [10], and its application to Gum Metal, a beta-Ti-based alloy, deformed
and measured in-situ within the elastic regime. In Section 3 we recall the
fundamentals of the Hodge decomposition and formulate the optimization problem
that allow us to extract the divergence-free part of a general non-
equilibrated stress field. In Section 4, we present the FFT-based method that
enables an efficient numerical resolution of the previously formulated
optimization problem when applied to a voxelized field, the type of data
structure that naturally results from data reduction of 3DXRD experiments.
Section 5 shows applications of the proposed method to synthetic piecewise
constant stress fields with a known ground truth, and stress fields in Gum
Metal measured by 3DXRD. In both cases, we show that the largest corrections
are obtained near grain boundaries. In Section 6, we draw conclusions and give
perspectives of the adoption of the proposed method in 3DXRD data inversion
packages.

##  2 Experimental characterization of elastic strain tensors

###  2.1 Data acquisition

The experimental dataset used in this work to demonstrate the proposed stress
filtering method was recorded in a regular Diffraction Contrast Tomography
(DCT) scan [2, 23]. The method uses a parallel monochromatic synchrotron X-ray
beam and records diffraction spots from individual grains of the
polycrystalline sample on a near-field area detector while the sample is
rotated continuously over $360^{\circ}$. The image stack is subject to a
number of pre-processing steps, at the end of which the diffraction spots are
segmented and indexed according to their grain of origin and their $(hkl)$
Miller indices. The indexing is based on crystallographic principles and the
knowledge of the crystal structure and approximate undeformed lattice cell
parameters of the one (or possible several) phases present in the polycrystal.
Each diffraction spot is stored as a $(u,v,w)$ volume array representing the
diffracted intensity distribution as a 3D scalar field: $u$ and $v$ being the
horizontal and vertical image coordinates, and $w$ is the rotation angle. The
use of monochromatic X-rays and a high-resolution detector provides
sensitivity to the local orientation and local unit cell parameters within the
crystallites. Local variations in lattice parameters, i.e. a change of the
shape and dimensions of the unit cell, in this case is interpreted as
mechanical elastic strain, and any effects from a possible local change of
chemical composition is neglected.

The results presented in this study were obtained from a polycrystalline Gum
Metal sample with a composition of Ti-36Nb-2Ta-3Zr-0.3O $wt\%$ and mean grain
diameter of $61\mu$m, from the ”low load” scan presented in [12]. Gum Metal is
able to sustain elastic strains up to $2\%$ elongation and above, thus its
nickname. A $500\mu$m tall gauge volume of a tensile specimen with an
approximately $600\mu$m wide rectangular cross-section was illuminated in its
entirety with a box beam during the $360^{\circ}$ rotation, using 40 keV beam
energy, 1.5 sec exposure time per image, $0.05^{\circ}$ angular step size,
resulting in 7200 images in total. The detector pixel size was $1.4\mu$m, and
the rotation axis to detector distance was 7 mm. The sample was mounted in a
miniature tensile rig [24] and a DCT scan at a uniaxial external tensile
stress of 34 MPa was acquired. On average 29 diffraction spots from the first
three families of reflections could be recorded for each of the $\sim 1430$
grains in the illuminated sample volume, which served as the input data for
the reconstruction.

###  2.2 Reconstruction of the sub-grain elastic deformation field

The first stage of the data analysis follows the now standard DCT processing
route, with a calibration procedure based on the dataset itself [10, 12]. The
centroid positions and mean orientations of the grains were determined based
on their diffraction spot centroid metadata. The initial 3D shapes of the
grains were reconstructed iteratively from the diffraction spots using the
SIRT algorithm, and assuming a constant orientation and elastic strain
distribution throughout the grain volume.

The deformation analysis applied here corresponds to an earlier version of the
novel Iterative Tensor Field (ITF) Reconstruction method described in detail
in [12], the main difference being that here no Tikhonov regularization term
was applied in the elastic deformation optimization target function (Equation
(39) in [12]). The ITF method was used to retrieve the complete local elastic
strain tensor (6 parameters), the local lattice misorientations from a
reference grain orientation (3 parameters) and refine the grain boundaries
over the entire 3D voxelated gauge volume in the specimen. ITF is based on a
vector representation of the 3D distribution of X-ray diffracting powers of
the crystallites $\mathbf{p}$ (one scalar per voxel) and a 3D intragranular
deformation field $\mathbf{d}$ (9 deformation components per voxel) of their
crystal lattice, and aims to reconstruct this
$\mathscr{F}(x,y,z):\mathbb{R}^{3}\mapsto\mathbb{R}^{1+9}$ tensor field on a
grain-by-grain basis in a static deformed state, measured in a single scan.
Mosaicity, i.e. a variation of crystal lattice orientations without associated
stress, is allowed in the solution, hence the compatibility equations are not
enforced. A diffraction spot is effectively treated as a projection of the 3D
grain shape into the 3D detector domain $(u,v,w)$, where the sub-grain
deformations (local misorientations and elastic strains) determine the
projection geometry that is unknown. ITF builds on the principles of
kinematical diffraction and ray tracing, and the solution for the elastic
deformation field is found by minimizing the differences between the measured
pixel intensities of the detector $\mathbf{q}^{mes}$ (ordered as a column
vector) and the simulated intensities $\mathbf{q}$ in the diffraction spots of
a grain in an iterative large scale optimization. This is a strongly non-
linear and usually under-determined problem but in each iteration step a local
linear problem at the latest solution can be formulated in matrix form:

$\displaystyle\Delta\mathbf{q}=\mathbf{q}^{mes}-\mathbf{q}=\boldsymbol{A}_{d}\,\Delta\mathbf{d}$
(1)

where $\Delta\mathbf{d}$ (a column vector) contains the corrections sought to
be applied to each elastic deformation component of each active grain voxel in
the current iteration step, either positive or negative. The
$\boldsymbol{A}_{d}$ matrix contains the intensity contributions to each pixel
(one row per pixel) from each elastic deformation component correction of the
active grain voxels. The $\boldsymbol{A}_{d}$ matrix is non-sparse, and it
needs to be recomputed frequently due to the strong non-linearity of the
model. Translated into a least-squares problem, the solver minimizes the
objective function $\Gamma$ which in this case is the square of the L2 norm
(denoted by $\left\|.\right\|_{2}$) of the residual vector:

$\displaystyle\Gamma=\left\|\Delta\mathbf{q}-\boldsymbol{A}_{d}\Delta\mathbf{d}\right\|_{2}^{2}$
(2)
$\displaystyle\Delta\mathbf{d}^{*}=\underset{\Delta\mathbf{d}}{\operatorname*{arg\,min}}(\Gamma)$
(3)

The elastic deformation corrections $\Delta\mathbf{d}^{*}$ found in the latest
iteration step are applied to the current elastic deformation field of the
grain. A 3D smoothing operation is then performed on the elastic deformation
field in each step, which effectively acts as regularization in the
optimization and results in a more realistic field. This smoothness suppresses
but does not eliminate the deviations from local mechanical equilibrium across
a grain volume.

The corrections to the diffracting powers, which determine the grain shapes,
are recomputed in a separate step in a similar way but less frequently, using
the latest elastic deformation field (i.e. the latest projection geometry).

Errors in the elastic deformation fields result from a number of sources in
the acquisition and data processing. The strain sensitivity here is expected
to be in the high single digits of $10^{-4}$. The overall deformation
sensitivity of the setup was limited by the angular step size, pixel size and
detector distance. The relatively small number of diffraction spots meant the
problem was more ill-posed. Some systematic deviations in intensity were
caused by the simplifications in the model of the image formation process.
Random noise from the photon counting statistics, fluorescence and scattering
processes was often significant. The grain-by-grain reconstruction approach is
inherently prone to higher errors near grain boundaries. These aspects and
potential future improvements of the acquisition and data processing,
including the handling of all grains simultaneously and ways to include the
local mechanical equilibrium, are described in more detail in [12]. The fit
quality in the ITF deformation solver can be monitored through the residual.
In this study, the fitting is considered to have failed in several percent of
the grains, mainly at the top and bottom of the gauge volume and at the free
surface. Although the dataset does not constitute the state-of-art measurement
capabilities, it is perfectly suitable to demonstrate the proposed stress
filtering method.

##  3 Proposed Method

###  3.1 Hodge decomposition

We use the following notation. Let $\Omega\subset{\mathbb{R}}^{3}$ denote a
rectangular domain and $H^{1}_{per}(\Omega,{\mathbb{R}})$ be the set of all
functions $u\in H^{1}_{loc}({\mathbb{R}}^{3},{\mathbb{R}})$ that are periodic
with period $\Omega$. The Helmholtz decomposition is a fundamental theorem in
vector calculus that has applications in many fields including elasticity,
incompressible flows and electromagnetism (see for example, [25]).

###### Theorem 1

Given any vector field $v\in L^{2}(\Omega,{\mathbb{R}}^{3})$, there exists a
scalar field $\varphi\in H^{1}_{per}(\Omega,{\mathbb{R}})$ and a vector field
$w\in H^{1}_{per}(\Omega,{\mathbb{R}}^{3})$, such that,

$v=\nabla\times w+\nabla\varphi,\quad\nabla\cdot w=0.$

Further $\varphi,w$ are unique up to a constant.

The Hodge decomposition is a generalization to tensor fields, and we use a
version due to Geymonat and Krasucki [22] for symmetric fields. Let
${\mathbb{M}}^{3}$ denotes the linear space of all second-order matrices and
${\mathbb{M}}_{sym}^{3}$ the linear space of all second-order symmetric
matrices.

###### Theorem 2

Given a symmetric matrix field $A\in L^{2}(\Omega,{\mathbb{M}}^{3}_{sym})$,
there exist a vector field $y\in H^{1}_{per}(\Omega,{\mathbb{R}}^{3})$, a
symmetric matrix field $H\in H^{2}(\Omega,{\mathbb{M}}^{3}_{sym})$ and a
constant matrix field $c_{0}$ such that,

$A=\nabla\times H\times\nabla+\frac{1}{2}(\nabla y+\nabla
y^{T})+c_{0},\quad\nabla\cdot H=0$

where $(\nabla\times
H\times\nabla)_{kn}=\epsilon_{ijk}\epsilon_{lmn}\nabla_{i}\nabla_{l}H_{jm}$ in
indicial notation. Further $y,H$ up to a constant, and $c_{0}$ are unique.

###  3.2 Problem formulation

Given a symmetric matrix field $S_{exp}\in
L^{2}(\Omega,{\mathbb{M}}^{3}_{sym})$, we want to find the closest field
$S_{equil}\in H^{1}(\Omega,{\mathbb{M}}^{3}_{sym})$ in $L^{2}$ norm, such that
$\nabla\cdot S_{equil}=0$. In other words,

$S_{equil}=\operatorname*{arg\,min}\left\\{||S-S_{exp}||_{2}:S\in
H^{1}_{per}(\Omega,{\mathbb{M}}^{3}_{sym}),\nabla\cdot S=0\right\\}.$

It is easy to verify using the divergence theorem that $\langle(\nabla
y+\nabla y^{T}),S\rangle_{H^{1}}=0$ for any $y\in
H^{1}_{per}(\Omega,{\mathbb{R}}^{3})$ and $S\in
H^{1}_{per}(\Omega,{\mathbb{M}})$ that satisfies $\nabla\cdot S=0$. Since
$\nabla\cdot S_{equil}=0$, it follows that the Hodge decomposition of
$S_{equil}$ has $y=0$, or

$S_{equil}=\nabla\times H\times\nabla+c_{0},\quad\nabla\cdot H=0$ (4)

for some $H\in H^{2}_{per}(\Omega,{\mathbb{M}}^{3}_{sym})$ and
$c_{0}\in{\mathbb{M}}^{3}_{sym}$. It follows that

$S_{exp}=S_{equil}+\frac{1}{2}(\nabla y+\nabla y^{T})=\nabla\times
H\times\nabla+c_{0}+\frac{1}{2}(\nabla y+\nabla y^{T}).$ (5)

Now, taking the average over $\Omega$ and again using the divergence theorem,
$c_{0}=\overline{S}_{exp}$ where we use the overhead bar to denote spatial
average $\overline{A}=\frac{1}{|\Omega|}\int_{\Omega}AdX$. Putting all this
together, the original minimization problem can then be written as

$\begin{split}H=\operatorname*{arg\,min}\left\\{{\mathscr{L}}(A):A\in
H^{2}_{per}(\Omega,{\mathbb{M}}^{3}_{sym}),\nabla\cdot A=0\right\\}\\\
\quad\quad\quad\text{where}\quad{\mathscr{L}}(A):=||\nabla\times
A\times\nabla-\widetilde{S}_{exp}||_{2}.\end{split}$ (6)

where tilde denotes the demeaned part
$\widetilde{S}_{exp}=S_{exp}-\overline{S}_{exp}$. We use this variational
problem to find $S_{equil}$.

The first variation

$\delta{\mathscr{L}}=\int_{\Omega}(\nabla\times
H\times\nabla-\widetilde{S}_{exp})\cdot(\nabla\times\delta H\times\nabla)\
dX=0\quad\forall\ \delta H\in{\mathscr{A}}$

where ${\mathscr{A}}=\\{A\in
H^{2}_{per}(\Omega,{\mathbb{M}}^{3}_{sym}):\nabla\cdot A=0\\}$. Notice that
for any $A,B\in H^{1}_{per}(\Omega,{\mathbb{M}}^{3}_{sym})$,

$\int_{\Omega}(\nabla\times A)\cdot B-A\cdot(\nabla\times B)\
dX=\int_{\Omega}\nabla\cdot(A\times B)\ dX=0.$

where the second equality follows from the divergence theorem. Thus,

$\int_{\Omega}(\nabla\times A)\cdot B\ dX=\int_{\Omega}A\cdot(\nabla\times B)\
dX.$

Applying this equality twice to the functional, we obtain

$\delta{\mathscr{L}}=\int_{\Omega}(\nabla\times\nabla\times
H\times\nabla\times\nabla-\nabla\times\widetilde{S}_{exp}\times\nabla)\cdot\delta
H\ dX=0\quad\forall\ \delta H\in{\mathscr{A}}.$

We obtain the Euler-Lagrange equation,

$\nabla\times\nabla\times
H\times\nabla\times\nabla=\nabla\times\widetilde{S}_{exp}\times\nabla.$

Since $\nabla\times(\nabla\times H)=\nabla(\nabla\cdot H)-\nabla^{2}H$ and
$\nabla\cdot H=0$, the Euler-Lagrange equation simplifies to

$\nabla^{4}H=\nabla\times\widetilde{S}_{exp}\times\nabla.$ (7)

We solve this equation to find $H$ and then $S_{equil}$. Recall that $H$ is
unique up to the constant, but the constant does not affect $S_{equil}$: so we
choose the constant to make $\overline{H}=0$.

###  3.3 Solution strategy

Given $S_{exp}$, we average over the volume to find $c_{0}$ and
$\widetilde{S}_{exp}$. We then solve Eq. (7) to find $H$. Finally, we use Eq.
(4) to obtain $S_{equil}$. The solution of Eq. (7) poses a challenge since it
involves the fourth order biharmonic equation for $H\in H^{2}$, which requires
continuous functions and derivatives. This in turn requires either higher
order elements using finite-element discretization, or higher order
differences in finite-difference methods. In either situation, it gives rise
to stiff numerical problems [26, 27]. Another approach is to break it up into
two harmonic problems and to solve them iteratively [28]. Here, we solve Eq.
(7) using Fourier transforms. This poses two challenges. The first is that the
experimental stress-field 4 is not periodic. To overcome this, a buffer is
filled around the voxelized representation of the material to make the domain
a periodic box. The properties of the buffer near the lateral surfaces are
chosen to give zero stress. We refer to a detailed discussion in section 2.
Second, the stress may suffer discontinuities across grain boundaries,
interfaces and across the boundaries of the periodic unit cell. A Fourier
series does not converge to discontinuous functions in the $L^{\infty}$ norm,
and this in turn leads to ringing artifacts in practice. We address this in
the next section.

##  4 Different operators in Fourier space

Let $\Omega=(0,L)^{3}$ to be a cube and let us discretize it with a $N\times
N\times N$ uniform cubic grid for $N$ even. The corresponding domain in the
Fourier space is $(-\pi/L,\pi/L)^{3}$ and again discretized uniformly with a
$N\times N\times N$ grid. For any $f:\Omega\to{\mathbb{R}}^{3}$, the discrete
Fourier transform (DFT) $\hat{f}$ satisfies

$f(i,j,k)=\sum_{l=1}^{N}\sum_{m=1}^{N}\sum_{n=1}^{N}\hat{f}(l,m,n)exp(\frac{2\pi\mathrm{i}}{L}((i-1)(l-1)+(j-1)(m-1)+(l-1)(k-1)),$
(8)
$\hat{f}(l,m,n)=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{N}f(i,j,k)exp(\frac{-2\pi\mathrm{i}}{L}((i-1)(l-1)+(j-1)(m-1)+(l-1)(k-1)),$
(9)

where $\mathrm{i}=\sqrt{-1}$ (distinct from index $i$).

###  4.1 Continuous differential operator

The continuous differential operator (CDO) of $f$ in Fourier domain is
obtained point-wise multiplication. For example, the partial derivative of f
with respect to x is

$\frac{\partial f(i,j,k)}{\partial
x}=\sum_{l=1}^{N}\sum_{m=1}^{N}\sum_{n=1}^{N}\mathrm{i}\xi_{x}\hat{f}(l,m,n)exp(\frac{2\pi\mathrm{i}}{L}((i-1)(l-1)+(j-1)(m-1)+(l-1)(k-1))$
(10)

where

$\xi_{x}=\begin{cases}&\frac{2\pi}{L}l,\quad\text{if }1\leq l<\frac{N}{2}\\\
&0,\quad\text{if }l=\frac{N}{2}\\\ &\frac{2\pi}{L}(l-N),\quad\text{if
}\frac{N}{2}<l\leq N\end{cases}$ (11)

Higher order derivatives are a direct composition of first derivatives, with
modification on the coefficients of the highest frequency term (see for
example [29]).

###  4.2 Discrete differential operator

We are interested in studying situations with grain boundaries and interfaces
across which the stress may suffer discontinuities. A Fourier series does not
converge to such a function in the $L^{\inf}$ norm even as $N\to\infty$. In
practice, we encounter ringing artifacts or spurious oscillations at finite
$N$. A solution to this problem is to construct the differential operator by
taking finite differences and then using DFT [30, 31]. The first order central
difference approximation reads,

$\frac{\partial f(i,j,k)}{\partial x}=\frac{f(i+1,j,k)-f(i-1,j,k)}{2\Delta h}$
(12)

etc., and therefore the discrete differential operator is given by Eq. (10)
with the substitution

$i\xi_{x}\leftarrow\frac{i}{\Delta h}\sin(\frac{2\pi(m-1)}{N}),$

and analogously for the other partial derivatives. For example,

$(\nabla\times
H\times\nabla)_{kn}=\epsilon_{ijk}\epsilon_{lmn}\nabla_{i}\nabla_{l}H_{jm},$

we substitute $\nabla_{x}$ with $\frac{i}{\Delta h}\sin(\frac{2\pi(m-1)}{N})$
for the discrete approximation. Note that this first order operator
approximates the continuous differential operator for small $\xi_{x}$, but
then decays thereby acting as a low-pass filter. This resolves the
overshooting and ringing at the interfaces. However, the convergence with
respect to mesh size is reduced to second order from the exponential
convergence of the CDO.

The biharmonic operator adopts the following 25 point stencil that is
consistent,

$\displaystyle\nabla^{4}f(i,j,k)$
$\displaystyle=42f(i,j,k)-12(f(i-2,j,k)+f(i+2,j,k)+f(i,j-2,k)+f(i,j+2,k)$
$\displaystyle+f(i,j,k-2)+f(i,j,k+2))+f(i-4,j,k)+f(i+4,j,k)+f(i,j-4,k)$
$\displaystyle+f(i,j+4,k)+f(i,j,k-4)+f(i,j,k+4)+2(f(i-2,j-2,k)$
$\displaystyle+f(i-2,j+2,k)+f(i+2,j-2,k)+f(i+2,j+2,k)+f(i,j-2,k-2)$
$\displaystyle+f(i,j-2,k+2)+f(i,j+2,k-2)+f(i,j+2,k+2)+f(i-2,j,k-2)$
$\displaystyle+f(i-2,j,k+2)+f(i+2,j,k-2)+f(i+2,j,k+2)).$

Notice that this is based on the second neighbor to make it consistent with
the first order difference , Eq. (12). We obtain the corresponding operator in
Fourier space from the usual Fourier space biharmonic operator with the
substitution

$\displaystyle\xi_{x}^{4}+\xi_{y}^{4}+\xi_{z}^{4}$
$\displaystyle\leftarrow(42-24(cos(4\pi(m-1)/N)+cos(4\pi(n-1)/N))+cos(4\pi(l-1)/N))$
$\displaystyle+2(cos(8\pi(m-1)/N)+cos(8\pi(n-1)/N)+cos(8\pi(l-1)/N))$
$\displaystyle+8(cos(4\pi(m-1)/N)cos(4\pi(n-1)/N)+cos(4\pi(n-1)/N)cos(4\pi(l-1)/N)$
$\displaystyle+cos(4\pi(m-1)/N)cos(4\pi(l-1)/N)))\frac{1}{16\Delta h^{4}}.$

##  5 Numerical Examples

We study two examples, one using synthetic data with a known ground truth and
another with experimental data.

###  5.1 Synthetic data with known ground truth

Figure 1: Synthetic data: (a) lateral stress field $\sigma_{11}$ obtained with
FFT-based model; (b) grain averaged $\sigma_{11}$; (c) divergence based on
CDO; (d) divergence based on DDO.

We begin with synthetic data from the results of a simulation using the FFT-
based method [32] to obtain an equilibrated stress field on a periodic unit
cell representing a Cu polycrystal deformed in tension in the elastic regime.
This polycrystalline unit cell is discretized by a $64^{3}$ grid on a
$[0,64]\mu m^{3}$ box, and contains 100 Cu single crystal grains generated by
periodic Voronoi tesselation. The Cu single crystal elastic constants reflect
their cubic anisotropy, with $C_{11}=168.4$ GPa, $C_{12}=121.4$ GPa and
$C_{44}=75.4$ GPa, which determines heterogeneous, piece-wise constant elastic
properties associated with the different crystal orientations, and, therefore,
a non-uniform stress field, when the polycrystal is subjected to load in the
elastic regime. For this analysis, the unit cell was loaded to a strain of
$10^{-5}$ in uniaxial tension along axis $x_{3}$, corresponding to a
longitudinal stress $\Sigma_{33}=6.32$ MPa, and zero lateral stresses:
$\Sigma_{11}=\Sigma_{22}=0$, applied to the unit cell. Figure 1(a) shows the
$\sigma_{11}$ component of the local stress field (in MPa), which fluctuates
with respect to the macroscopic value $\Sigma_{11}=0$. Since this synthetic
data was generated using a physically meaningful model that uses CDO to solve
the governing equations of micro-mechanics (stress equilibrium and strain
compatibility), the resulting stress field is indeed equilibrated (and thus
divergence-free with respect to the CDO).

We start with this equilibrated stress field and obtain a grain-averaged
stress field by averaging the stress over the grains – Figure 1(b) shows the
$\sigma_{11}$ component of the grain-averaged stress. This piece-wise constant
field is not divergence-free, and the divergence is concentrated at the grain
boundaries. We compute the divergence of the grain-averaged stress using both
CDO and DDO – the results are shown in Figure 1(c) and Figure 1(d)
respectively. We observe that the divergence is generally concentrated at the
grain boundaries, but CDO smears it into the grain interior. This is a result
of the spurious oscillation and overshooting that is characteristic of CDO on
discontinuous functions.

Figure 2: Recovery from grain averaged data, $\sigma_{11}$. (a-c) Results of recovery with CDO: (a) recovered stress field, (b) difference between ground truth an recovered stress field and (c) difference between grain averaged and recovered stress field. (d-f) Results of recovery with DDO: (a) recovered stress field, (b) difference between ground truth an recovered stress field and (c) difference between grain averaged and recovered stress field. (g) $\sigma_{11}$ along the line $0.5\times 0.5\times[0,1]$ | Diff with ground truth | Diff with average
---|---|---
CDO ($L^{2}$/$L^{\infty}$) | 0.0024/4.6041 | 0.0018/2.6945
DDO ($L^{2}$/$L^{\infty}$) | 0.0025/4.6007 | 0.0018/2.6222
Table 1: Error in the recovery from grain averaged data

We now study if it is possible to obtain the original stress field from the
grain-averaged data using our method. We apply the proposed algorithm to the
grain-averaged stress field using both CDO and DDO. The results are shown in
Figure 2. We see that the recovery is imperfect, and this is confirmed in
Table 1 which lists the $L^{2}$ and $L^{\infty}$ norms of the error
($|f(x)|_{L^{\infty}}=sup\\{f(x)\text{$\forall$ $x$ in $\Omega$}\\}$).
However, both CDO and DDO yield similar errors. In fact, Figure 2(g) shows the
$\sigma_{11}$ component of the stress along the line $0.5\times
0.5\times[0,1])$. Both operators recover similar fields though the CDO has
more oscillations within the grains and an overshoot at grain boundaries (see
for example the grain boundary around $z=0.9$).

To understand the error in recovery, recall that our algorithm only filters
out the symmetrized gradient (curl-free) part and keeps all the divergence-
free part of the given stress-field. So, if the difference between the given
data and the ground truth deviate by both curl-free and divergence-free
fields, our algorithm will filter out the symmetrized gradient portion of the
deviation but retain the divergence-free part of that deviation. In our
example, the difference between grain-averaged stress field and ground truth
contains both curl-free and divergence-free components. Our algorithm filtered
out the former but not the latter and this is the error in both Figure 2 and
Table 1. We have verified this by taking the difference between grain-averaged
stress field and ground truth, and applying our algorithm to it. The resulting
residual equals (close to machine precision) to the recovery error (difference
between the ground truth and recovered) in both CDO and DDO.

Figure 3: (a,b) Recovery from a perturbation that is a symmetrized gradient
(curl-free) (a) perturbed and (b) recovered fields. Compared to the ground
truth in Figure 1(a), the recovery is very good – $9.016\times 10^{-6}$ in
$L^{2}$ norm and $0.014$ in $L^{inf}$ norm. (c,d) Recovery from a perturbation
that is divergence-free (c) perturbed and (d) recovered fields. The recovered
field agrees with the perturbed field.

We further verify this in Figure 3. We perturb the ground truth by adding the
symmetrized gradient $\nabla a+(\nabla a)^{T}$ of some vector field $a$ – see
Figure 3(a). For a comparison of the scale, $|\sigma_{11}|_{L^{\infty}}$ of
the original field is 6.02 MPa. $|\sigma_{11}|_{L^{\infty}}$ of the
perturbation is 1.00 MPa. Our algorithm is able to recover the ground truth
from the perturbation – see Figure 3(b) with minimal error: $9.016\times
10^{-6}$ in $L^{2}$ norm and $0.014$ in $L^{inf}$ norm. On the other hand,
when we perturb the ground truth by adding a divergence free field (Figure
3(c)), our algorithm returns the perturbed field and not the ground truth
(Figure 3(d)). For a comparison of the scale, $|\sigma_{11}|_{L^{\infty}}$ of
the original field is 6.02 MPa. $|\sigma_{11}|_{L^{\infty}}$ of the
perturbation is 4.81 MPa. We note that whether the perturbation remains does
not depend on the scale of perturbation.

In conclusion, our approach finds the best projection to divergence-free
fields. Therefore, it is able to filter out artifacts due to fields that are
symmetrized gradients/curl-free, but unable to filter out artifacts that
happen to be divergence-free. Indeed, it is not possible to filter out
divergence-free artifacts without additional knowledge in terms of material
behavior. For this reason, describing the average stress in each grain only
provides limited information about the actual state of stress in the material.

###  5.2 Experimental Sample

In this section, we apply the proposed filtering methodology to the DCT
experimental data on Gum Metal processed with the ITF method to obtain
voxelized stress fields as described in Section 2. We focus on the measured
stress field from a sample subjected to an average normal stress of 34 MPa
along the $z$-axis while the other two lateral surfaces are stress free. In
order to use our algorithm that assumes periodicity, we introduce a buffer
region around the specimen. The buffer region on the side (along $x$-axis and
$y$-axis) has zero stress to accommodate the free surface. We also note that
we implicitly make the top and bottom to be periodically connected. This
spurious connection between top and bottom in the longitudinal direction
evidently deteriorates the stress balance near those regions. However, we
should note that the experimental measurement is also less reliable in those
locations and in grains on the free surface of the sample that have high
mosaicity as a result of the spark cutting process. This treatment still
preserves the information in the bulk of the sample. The resulting unit cell
of $[0,640\mu m]^{3}$, resolved at $256^{3}$ voxels is shown in Figure 4(a,b)
along with the $\sigma_{11}$ and $\sigma_{33}$ components of the measured
stress. We note that the state of stress is quite heterogeneous. Figure 4(c)
shows the magnitude of the divergence of the measured stress-field calculated
using DDO, and we observe high values of divergence at the grain boundaries
and triple junctions. This shows the the difficulty of accurate measurement of
lattice strains and thus stresses where lattice structure is not preserved.

We present the filtering results with DDO (Fig. 5). After filtering, the peak
value is reduced. Most of the profile patterns are preserved. Comparing the
difference, we notice that most of the region undergoes minor modification.
Significant differences mostly occur around the grain boundaries. We also
observe a significant correction of stress in the upper left grain, likely
because the deformation solver failed for this grain and some other grains
close to the surface. The main trends in Figure 5(g) are well preserved though
and most of the smaller features within the grains are retained, even if the
stress magnitudes changed. This is an encouraging result, and suggests that
the elastic deformation solver based on diffraction constraints returns
realistic trends in the stress fields. The difference in the actual magnitudes
of this measurement are somewhat less concerning, as these were less reliable
due to the deformation resolution of the experimental setup not being fully
adequate. The deformation sensitivity of the setup will be improved in the
future.

There are a number of outlier grains where the strain/stress state is extreme
and significantly deviates from their neighbours. The ITF deformation fitting
has potentially failed in these grains. A notable result of the filtering is
that most of those grains are retained as outliers. Such high local deviations
do not seem realistic, although their validity cannot be ruled out. Thus, it
is not clear whether the filtered field is indeed more accurate in those
regions. The measurement errors might be so high in those grains that the
filtering cannot be expected to correct for them.

A detailed comparison of the two operators is conducted by looking closely at
the stress distribution along line $320\mu m\times 320\mu m\times[0,640]\mu
m)$ (Fig. 5)(g). Both methods yield similar results. In some grains, the
stress is not corrected much. But in the grain near $z\approx 300\mu m$, there
is a significant correction. The main difference between these two operators
lie in grain boundaries. CDO and DDO both have some oscillations around the
boundaries (e.g. $z\approx 250\mu m$). But CDO provides a higher estimation at
the boundaries. In contrast, DDO has a smoother transition but it is
underestimating the values at the interface.

Figure 4: Stress field as measured by DCT-ITF: (a) $\sigma_{11}$ (transparent
box is the computational domain); (b) $\sigma_{33}$. (c) Magnitude of the
divergence based on DDO.

Figure 5: Stress distribution: (a-c) $\sigma_{11}$; (d-f) $\sigma_{33}$; (g-i)
$\sigma_{11}$, $\sigma_{22}$, $\sigma_{33}$ along the line $320\mu m\times
320\mu m\times[0,640]\mu m$

##  6 Conclusions

We developed a FFT-based Hodge decomposition method for symmetric matrix
fields. To solve the problem of oscillations at interfaces and preserve the
orthogonality of operators, we further developed a modified discrete
differential operator. We applied the Hodge decomposition to recover a
equilibrated stress field from experimental stress measurements. Since the
method is based on a projection, it is able to filter out artifacts due to
fields that are symmetrized gradients (curl-free), but unable to filter out
artifacts that happen to be divergence-free. Indeed, it is not possible to
filter out divergence-free artifacts without additional knowledge in terms of
material behavior. We show that the algorithm can be effectively applied to
experimentally obtained stress fields. We have also investigated the
difference between two discrete differential operators, CDO and DDO. They
provide quantitatively similar results overall. However, DDO performs better
at interfaces damping out spurious oscillations.

The proposed methodology requires knowledge of the single crystal elastic
constants. If they are not precisely known, the method could be extended to
refine their values as part of the optimization procedure. 3DXRD
reconstruction where grain centroids and average elastic strains/stresses are
determined is now a fairly standard post-experimental analysis, and the
proposed stress filtering can be directly combined with those reconstructions
techniques as post-optimization correction. On the other hand, sub-grain
methods are novel and still under development, and the proposed methodology
could be combined with them in different ways. The sequential methodology
demonstrated here, based on the application of the FFT-based micromechanical
filter to stress fields resulting from diffraction-based optimization, could
be extended to fulfill the diffraction and micromechanical constraints
concurrently during data inversion and refinement. In both cases, but
specially in the latter (concurrent optimization), the extension and
application of the method has a great potential to increase both deformation
and spatial accuracy, as they are interlinked, including the quantitative
estimation of errors of the diffraction measurements.

In terms of numerical performance of the proposed method, since it is based on
the efficient FFT computations, a single application of the filter to a
$256^{3}$ voxel image takes 473 seconds in an Intel Core i7-8700 CPU machine.
This high efficiency of a single application of the filter bodes well for
standing up the more numerically involved (but still realizable) concurrent
optimization, mentioned above.

As a final comment, while the proposed method is based on imposing one type of
micro-mechanical constraint on the measured stress field, i.e. stress
equilibrium, the dual constraint is strain compatibility. Such condition can
also be expressed in terms of differential operators acting, in this case, on
the total strain field. Total strains are not directly accessible by the
existing diffraction-based data reduction methods, since the non-elastic part
of the deformation does not directly affect lattice spacing and, consequently,
diffraction patterns. However, if new methods for determining total strain
fields using 3DXRD experiments are further developed, e.g. using evolving 3D
microstructural images to perform DCT-based Digital Volume Correlation (DVC),
in what would be an extension of similar techniques based on Computed
Tomography (CT), strain compatibility imposed via FFT-based methods should be
certainly considered as part of the data inversion procedures.

## Acknowledgement

The core part of this work was conducted during HZ’s visit to Los Alamos
National Laboratory (LANL) to work with RAL under LANL’s LDRD program, whose
support is acknowledged. RAL’s work is also supported by LANL’s Science
Campaign 2 program. KB and HZ were also partially supported by the Air Force
Office of Scientific Research through the MURI grant FA9550-16-1-0566.

## References

  * [1] Henning Friis Poulsen.  Three-dimensional X-ray diffraction microscopy: mapping polycrystals and their dynamics, volume 205.  Springer Science & Business Media, 2004.
  * [2] Wolfgang Ludwig, P Reischig, Andrew King, Michael Herbig, EM Lauridsen, G Johnson, TJ Marrow, and Jean-Yves Buffiere.  Three-dimensional grain mapping by x-ray diffraction contrast tomography and the use of friedel pairs in diffraction data analysis.  Review of scientific instruments, 80(3):033905, 2009.
  * [3] U Lienert, SF Li, CM Hefferan, J Lind, RM Suter, JV Bernier, NR Barton, MC Brandes, MJ Mills, MP Miller, et al.  High-energy diffraction microscopy at the advanced photon source.  Jom, 63(7):70–77, 2011.
  * [4] RM Suter, D Hennessy, C Xiao, and U Lienert.  Forward modeling method for microstructure reconstruction using x-ray diffraction microscopy: Single-crystal verification.  Review of scientific instruments, 77(12):123905, 2006.
  * [5] SF Li and RM Suter.  Adaptive reconstruction method for three-dimensional orientation imaging.  Journal of Applied Crystallography, 46(2):512–524, 2013.
  * [6] Nicola Viganò, Alexandre Tanguy, Simon Hallais, Alexandre Dimanov, Michel Bornert, Kees Joost Batenburg, and Wolfgang Ludwig.  Three-dimensional full-field x-ray orientation microscopy.  Scientific reports, 6(1):1–9, 2016.
  * [7] Henning Friis Poulsen, Søren Fæster Nielsen, Erik M Lauridsen, Søren Schmidt, RM Suter, U Lienert, L Margulies, T Lorentzen, and D Juul Jensen.  Three-dimensional maps of grain boundaries and the stress state of individual grains in polycrystals and powders.  Journal of applied crystallography, 34(6):751–756, 2001.
  * [8] Joel Vincent Bernier, Nathan Rhodes Barton, Ulrich Lienert, and Matthew Peter Miller.  Far-field high-energy diffraction microscopy: a tool for intergranular orientation and strain analysis.  The Journal of Strain Analysis for Engineering Design, 46(7):527–547, 2011.
  * [9] Jette Oddershede, Søren Schmidt, Henning Friis Poulsen, Henning Osholm Sørensen, Jonathan Wright, and Walter Reimers.  Determining grain resolved stresses in polycrystalline materials using three-dimensional x-ray diffraction.  Journal of Applied Crystallography, 43(3):539–549, 2010.
  * [10] Péter Reischig and Wolfgang Ludwig.  Three-dimensional sub-grain mapping of elastic strain state and orientations in the bulk of polycrystals.  In 40th Risoe International Symposium on Metal Microstructures in 2D, 3D and 4D, 2019.
  * [11] Yujiro Hayashi, Daigo Setoyama, Yoshiharu Hirose, Tomoyuki Yoshida, and Hidehiko Kimura.  Intragranular three-dimensional stress tensor fields in plastically deformed polycrystals.  Science (New York, N.Y.), 366(6472):1492–1496, 2019.
  * [12] Péter Reischig and Wolfgang Ludwig.  Three-dimensional reconstruction of intragranular strain and orientation in polycrystals by near-field x-ray diffraction.  Current Opinion in Solid State and Materials Science, page 100851, 2020.
  * [13] Yu-Feng Shen, He Liu, and Robert M Suter.  Voxel-based strain tensors from near-field high energy diffraction microscopy.  Current Opinion in Solid State and Materials Science, 24(4):100852, 2020.
  * [14] N. Axel Henningsson, Stephen A. Hall, Jonathan P. Wright, and Johan Hektor.  Reconstructing intragranular strain fields in polycrystalline materials from scanning 3DXRD data.  Journal of Applied Crystallography, 53(2):314–325, apr 2020.
  * [15] Reeju Pokharel and Ricardo A Lebensohn.  Instantiation of crystal plasticity simulations for micromechanical modelling with direct input from microstructural data collected at light sources.  Scripta Materialia, 132:73–77, 2017.
  * [16] K Chatterjee, JYP Ko, JT Weiss, HT Philipp, J Becker, P Purohit, SM Gruner, and Armand Joseph Beaudoin.  Study of residual stresses in ti-7al using theory and experiments.  Journal of the Mechanics and Physics of Solids, 109:95–116, 2017\.
  * [17] Darren C Pagan and Armand J Beaudoin.  Utilizing a novel lattice orientation based stress characterization method to study stress fields of shear bands.  Journal of the Mechanics and Physics of Solids, 128:105–116, 2019\.
  * [18] Kevin P McNelis, Paul R Dawson, and Matthew P Miller.  A two-scale methodology for determining the residual stresses in polycrystalline solids using high energy x-ray diffraction data.  Journal of the Mechanics and Physics of Solids, 61(2):428–449, 2013\.
  * [19] E Demir, J-S Park, MP Miller, and PR Dawson.  A computational framework for evaluating residual stress distributions from diffraction-based lattice strain data.  Computer Methods in Applied Mechanics and Engineering, 265:120–135, 2013.
  * [20] J-S Park, Ulrich Lienert, Paul R Dawson, and Matthew P Miller.  Quantifying three-dimensional residual stress distributions using spatially-resolved diffraction measurements and finite element based data reduction.  Experimental Mechanics, 53(9):1491–1507, 2013.
  * [21] J-S Park, K Yildizli, Eralp Demir, Paul Richard Dawson, and Matthew Peter Miller.  Non-destructive characterization of subsurface residual stress fields and correlation with microstructural conditions in a shot-peened inconel component.  Experimental Mechanics, 58(9):1389–1406, 2018.
  * [22] Giuseppe Geymonat and Françoise Krasucki.  Hodge decomposition for symmetric matrix fields and the elasticity complex in lipschitz domains.  Communications on Pure & Applied Analysis, pages 295–309, 2009\.
  * [23] Péter Reischig, Andrew King, Laura Nervo, Nicola Viganó, Yoann Guilhem, Willem Jan Palenstijn, K. Joost Batenburg, Michael Preuss, and Wolfgang Ludwig.  Advances in X-ray diffraction contrast tomography: Flexibility in the setup geometry and application to multiphase materials.  Journal of Applied Crystallography, 46(2):297–311, 2013.
  * [24] N Gueninchault, H Proudhon, and W Ludwig.  Nanox: a miniature mechanical stress rig designed for near-field X-ray diffraction imaging techniques.  Journal of Synchrotron Radiation, 23(6):1474–1483, nov 2016.
  * [25] Jason Cantarella, Dennis DeTurck, and Herman Gluck.  Vector calculus and the topology of domains in 3-space.  The American mathematical monthly, 109(5):409–442, 2002.
  * [26] John H Argyris and Padraic C Dunne.  The finite element method applied to fluid mechanics.  In: Computational methods and problems in aeronautical fluid dynamics.(A77-12551 02-02) London and New York, pages 158–197, 1976.
  * [27] Murli M Gupta and Ram P Manohar.  Direct solution of the biharmonic equation using noncoupled approach.  Journal of Computational Physics, 33(2):236–248, 1979.
  * [28] Xiao-liang Cheng, Weimin Han, and Hong-ci Huang.  Some mixed finite element methods for biharmonic equation.  Journal of computational and applied mathematics, 126(1-2):91–109, 2000.
  * [29] Steven G Johnson.  Notes on fft-based differentiation.  MIT Applied Mathematics, 2011.
  * [30] Stéphane Berbenni, Vincent Taupin, Komlan Sénam Djaka, and Claude Fressengeas.  A numerical spectral approach for solving elasto-static field dislocation and g-disclination mechanics.  International Journal of Solids and Structures, 51(23-24):4157–4175, 2014.
  * [31] Ananthan Vidyasagar, WL Tan, and Dennis M Kochmann.  Predicting the effective response of bulk polycrystalline ferroelectric ceramics via improved spectral phase field methods.  Journal of the Mechanics and Physics of Solids, 106:133–151, 2017\.
  * [32] Ricardo A Lebensohn, Anand K Kanjarla, and Philip Eisenlohr.  An elasto-viscoplastic formulation based on fast fourier transforms for the prediction of micromechanical fields in polycrystalline materials.  International Journal of Plasticity, 32:59–69, 2012.


#  Macroscopic effects of localised measurements
in jammed states of quantum spin chains

Kemal Bidzhiev  Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France
Maurizio Fagotti<EMAIL_ADDRESS>Université
Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France  Lenart Zadnik  Université
Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France

###### Abstract

A quantum jammed state can be seen as a state where the phase space available
to particles shrinks to zero, an interpretation quite accurate in integrable
systems, where stable quasiparticles scatter elastically. We consider the
integrable dual folded XXZ model, which is equivalent to the XXZ model in the
limit of large anisotropy. We perform a jamming-breaking localised measurement
in a jammed state. We find that jamming is locally restored, but local
observables exhibit nontrivial time evolution on macroscopic, ballistic
scales, without ever relaxing back to their initial values.

### Introduction.

Recent remarkable progress in quantum information science owes much to the
advancement of cold-atom and trapped-ion setups [1, 2, 3, 4]. These have
become so refined [5, 6, 7, 8] to allow for the design of so-called “quantum
simulators” [9]. Understanding the effects of perturbations and measurements
has always been crucial in this regard, be it due to the role of perturbations
in the decoherence process, or the one of measurements in simulation protocols
[10]. Integrable [11, 12, 13, 14] and generic [15, 16, 17, 18, 19, 20] quantum
circuits are potentially useful in this endeavour: they allow for some degree
of exact treatment [21, 22, 23, 24, 25] and befit experimental realisation
[26, 27].

Localised measurements can be viewed as a type of “quantum quench”, i.e., the
non-equilibrium dynamics induced by a sudden perturbation, studied in the last
decades especially in the context of relaxation in quantum many-body systems
(see [28] and articles therein). Such perturbations break the homogeneity of
the system, making its study challenging. In integrable models the most
effective large-scale description of the dynamics in the presence of
inhomogeneities is arguably the so-called “generalised hydrodynamics” (GHD)
[29, 30, 31]. Although GHD correctly predicts the large-scale dynamics of
local observables in numerous quench protocols (e.g., in the two-temperatures
scenario), the information provided by the theory is sometimes incomplete. The
first example of this kind was exhibited in Ref. [32], considering the massive
Heisenberg model: the ingredients of GHD are blind to observables that are odd
under spin flip, entailing the inclusion of an additional independent
continuity equation. An even more striking example was considered in Refs [33,
34, 35, 36, 37], in which GHD keeps a symmetry that is instead broken in the
thermodynamic limit: observables not respecting that symmetry are affected by
a class of localised perturbations at arbitrarily long times and large
distances from the inhomogeneity.

In this Letter we study the effect of a localised projective measurement in a
quantum jammed state. To this end we consider a modification of the so-called
“wing-flap protocol” (see Fig. 1) introduced in Ref. [38] to provide insight
into the quantum information scrambling. The system is prepared in a low-
entangled stationary state, which, we assume, is also an eigenstate of the
operator measured by Alice. We then go back in time (a trivial step, since the
state is stationary) considering an alternative history in which, at some
unknown but fixed ancient time, Bob had performed an unknown projective
quantum measurement at an unknown but fixed large distance from Alice.
Returning to an alternative present we wonder, following Ref. [39], whether
Alice can still recover the information she had without Bob’s intervention.

$e^{iHt}$$e^{-iHt}$$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$$W$$O$$O$BobAliceAlice$(1)$$(2)$$(3)$$(4)$$(5)$
Figure 1: The “wing-flap” protocol: (1) Alice measures an observable $O$; (2)
the system evolves for a time $-t$; (3) Bob applies a “wing-flap” perturbation
$W$; (4) the system evolves for a time $t$; (5) Alice measures $O$ again. In
Ref. [39] $W$ is a “blind” local projective quantum measurement at the
position of $O$ (no information is kept), and step (5) is upgraded to a
quantum state tomography. There it was shown that, if $H$ describes a generic
system, the local state can be recovered with a limited amount of effort. We
simplify the protocol by trivialising the part highlighted in grey.

Generally, in a shift invariant quantum spin chain the distant Bob’s
measurement is not expected to have any effect on Alice’s subsystem. In
particular, Ref. [40] showed that, in noninteracting spin chains, relaxation
to a Generalised Gibbs ensemble (GGE) is not compromised by a localised
perturbation. In our setting this implies that, at late times, the effect of
Bob’s measurement on a finite subsystem fades away. We are not aware of any
physical argument against the generalisation of this result to interacting
integrable systems. Indeed, this conclusion is supported by numerical
investigations, which generally show the irrelevance of localised
perturbations on the state after long-enough time. Consider, for example, the
numerical tests of the generalised hydrodynamic predictions of time evolution
after two different states are joined together [41]. The tests always differ
in the way the states are joined, nevertheless the asymptotic behaviours at
large times match. In generic systems scrambling is even more pronounced, so a
very distant quantum measurement in a low-entangled stationary state is not
expected to have visible effects however large the time is. Incidentally, even
if the state were not stationary and Bob measured an observable in Alice’s
subsystem, Alice would be able to recover the original local state [39].

Here we apply this protocol to the dual folded XXZ spin-$1/2$ chain, which
belongs to a class of effective models emerging in strong coupling limits of
spin chains [42]. It corresponds to a special point of the two-component
Bariev model [43] and is described by the Hamiltonian

$H=J\sum_{\ell}(\sigma_{\ell}^{x}\sigma_{\ell+2}^{x}+\sigma_{\ell}^{y}\sigma_{\ell+2}^{y})\frac{\mathbbm{1}-\sigma_{\ell+1}^{z}}{2}\,.$
(1)

The model can be solved exactly by introducing two species of particles
associated with spins up on either even or odd sites [42] (see also [44]). The
solution is based on the observation that the unique nontrivial effect of $H$
corresponds to moving a spin up by two sites when it is adjacent to two spins
down:

In a configuration in which all spins down are isolated no hopping process can
occur: the aforementioned particles cannot move because they are stacked
together. There are exponentially many such configurations with clustering
properties, and we refer to them as jammed states.

We report a striking exception to the empirical rule that a distant localised
perturbation in a low-entangled state of a quantum many-body system described
by a translationally invariant Hamiltonian does not have visible effects at
infinitely large times and distances. To the best of our knowledge, only one
exception has been reported so far and it concerns systems prepared in the
ground state when a discrete symmetry is spontaneously broken [33, 34, 35, 36,
37]. Similarly, we consider an initial state that breaks the conserved charges
that completely characterise a basis of energy eigenstates, but, differently,
our initial state belongs to an _exponentially_ large degenerate sector. We
will show that, in the limit of large time, the jammed sector is
asymptotically stable under the wing-flap protocol. Quite exceptionally, we
also demonstrate that the measurement results in a macroscopic change of the
spin profiles on ballistic scales, namely its effect does not fade away at
large times, the expectation values of local observables instead approaching
nontrivial functions of the ratio between distance and time.

### Quantum jamming.

Our discussion is specialised to the sector spanned by states $\ket{\Phi}$
satisfying the jamming condition

$\mathcal{P}_{\downarrow\downarrow}(\ell)\equiv\bra{\Phi}\frac{\mathbbm{1}-\sigma_{\ell}^{z}}{2}\frac{\mathbbm{1}-\sigma_{\ell+1}^{z}}{2}\ket{\Phi}=0,\,\,\,\forall\ell\,.$
(2)

Due to (Introduction.), such states are clearly eigenstates of $H$ with zero
energy and are jammed (see also Refs [42, 45, 46]). Jammed states are special
since they can break the two-site shift invariance of the complete set of
charges exhibited in Ref. [42]. This is possible because the spectrum has huge
degeneracies associated with hidden symmetries, which have been shown, for
example, to play a key role in the quench dynamics within the non-interacting
sectors of the model [47].

The simplest basis of the jammed sector consists of product states that are
eigenstates of $\\{\sigma_{\ell}^{z}\\}$, but one can easily construct states
with any entanglement entropy density up to $\frac{1}{2}\log 2$. For example,
if $o$ denotes the sublattice of odd sites and $e$ the one of even sites, the
state $\ket{\uparrow\ldots\uparrow}_{o}\otimes\ket{\psi}_{e}$ is jammed and
the entanglement entropy of a spin block is half of that in the state
$\ket{\psi}$. Particularly interesting are $(2n)$-site shift invariant jammed
product states that are not eigenstates of $\\{\sigma_{\ell}^{z}\\}$, e.g.,

$\ket{\mathcal{L}_{n}}=e^{i\frac{\pi}{2n}\sum_{j}j\sigma_{j}^{z}}\ket{\uparrow\leftarrow\dots\uparrow\leftarrow}\,.$
(3)

This family of stationary states breaks the $U(1)$ symmetry generated by
$\sum_{\ell}\sigma_{\ell}^{z}$ and, for $n>1$, also the 2-site shift
invariance of the model’s charges. Such states belong to the fully interacting
sector of the model, which is characterised by the presence of spins up on
both even and odd positions [42].

The jammed sector is invariant under measurements of operators that commute
with all elements of the set $\\{\sigma_{\ell}^{z}\\}$. Note instead that
measuring other observables generally results in leaving the sector.

### Locally quasi-jammed states.

Criterion (2) can be extended to inhomogeneous states that are jammed
asymptotically in some scaling limit, for example,

$\lim_{\ell,t\to\infty\atop\zeta=\ell/t\text{
fixed}}\braket{\Phi_{t}}{\frac{\mathbbm{1}-\sigma_{\ell}^{z}}{2}\frac{\mathbbm{1}-\sigma_{\ell+1}^{z}}{2}}{\Phi_{t}}=0,\qquad\forall\zeta\,.$
(4)

We call them locally quasi-jammed states (LQJS).

We will focus on initial states belonging to the family
$\ket{\mathcal{L}_{n}}$ and, for the sake of simplicity, assume Alice to
measure
$\cos(\frac{j\pi}{n})\sigma_{j}^{x}-\sin(\frac{j\pi}{n})\sigma^{y}_{j}$, with
$j$ even, so that $\ket{\mathcal{L}_{n}}$ is not affected by the measurement.
The state after Alice’s measurement is then of the form
$\ket{\Phi}=\ket{\uparrow\ldots\uparrow}_{o}\otimes\ket{\varphi_{0}}_{e}$. Let
then Bob be in an odd position $2\ell^{\prime}-1$ and perform a blind
measurement of the spin. Since the spin is up before the measurement, the
density matrix after the measurement, $\rho(0)$, will commute with
$\sigma^{z}_{2\ell^{\prime}-1}$; in particular it reads

$\rho(0)=\frac{2}{3}\ket{\Phi}\bra{\Phi}+\frac{1}{3}\sigma_{2\ell^{\prime}-1}^{x}\ket{\Phi}\bra{\Phi}\sigma_{2\ell^{\prime}-1}^{x}\,.$
(5)

The coefficients $2/3$ and $1/3$ come from considering a (blind) projective
measurement – see the Supplementary material (SM) [48]. Nevertheless, the
structure of the density matrix would remain the same also under more general
measurement protocols. A local measurement can break condition (2) only
locally and, because the defect is moved over the jammed state as a
quasiparticle would be moved under the effect of a hopping Hamiltonian (in
some nontrivial background – see also Ref. [49]), it is reasonable to expect
the validity of (2) in the limit of large time. We can then foresee an LQJS to
emerge, i.e., a state where (4) holds. This expectation is supported by our
numerical investigation – see Fig. 2.

Figure 2: Profile of the rescaled jamming criterion
$t\,\mathcal{P}_{\downarrow\downarrow}(\ell/t)$ at times $t=10,50$ after the
flip of the central spin up (in an odd position) of the state
$\ket{\mathcal{L}_{1}}=\ket{\uparrow\leftarrow\uparrow\leftarrow\ldots}$. The
data inside the light cone are in excellent agreement with conjecture (18).
The jamming condition is locally restored as
$\mathcal{P}_{\downarrow\downarrow}(\ell/t)\sim 1/t$.

Notwithstanding the initial state being shift invariant by a finite number of
sites, Bob’s measurement locally breaks that symmetry. We find that, while
shift invariance is locally restored, it remains globally broken, as in the
systems that can be described by generalised hydrodynamics. Contrary to the
latter case, however, the two-site shift invariance, which is a symmetry of a
complete set of local conservation laws, is generally not restored, even
locally. Figure 3, for instance, shows that the $xy$-plane spin profiles
remain staggered on the sublattice of even sites, however long the time after
the initial perturbation to the state $\ket{\mathcal{L}_{2}}$ is: only
$4$-site shift invariance is locally restored. This property is shared with
non-abelian integrable systems like the quantum XY model, where the late time
behaviour is generally not captured by a macro-state characterised by the
complete set of one-site shift invariant charges [47].

Figure 3: Spin profiles $\braket{\sigma^{x}}_{\zeta}$ (red) and
$\braket{\sigma^{y}}_{\zeta}$ (blue) at times $t=10,50$ after having flipped
the central spin up (in an odd position) of the state
$\ket{\mathcal{L}_{2}}=\ket{\uparrow\leftarrow\uparrow\rightarrow\ldots}$. The
profile of $\sigma^{z}$ asymptotically approaches zero. Nontrivial profiles
emerge on _even_ sites, whereas spins on odd positions relax to their original
values along the $z$ direction, and are not shown. The data for
$\langle\sigma^{y}\rangle_{\zeta}$ are in excellent agreement with conjecture
(18).

### Dynamics after a local spin flip.

Since $\ket{\Phi}$ is an eigenstate, the time evolution of (5) can be
immediately deduced from that of

$\ket{\Psi(0)}=\sigma_{2\ell^{\prime}-1}^{x}\ket{\Phi}\equiv(\sigma^{x}_{0}\ket{\uparrow\ldots\uparrow})_{o}\otimes\ket{\varphi_{0}}_{e}\,,$
(6)

where, without loss of generality, we have set $\ell^{\prime}=0$. The state at
time $t$ can then be represented as follows

$\ket{\Psi(t)}=\sum_{n}(\sigma^{x}_{n}\ket{\uparrow\ldots\uparrow})_{o}\otimes\ket{\varphi_{n}(t)}_{e}\,,$
(7)

where the unnormalized wave functions $\ket{\varphi_{n}(t)}_{e}$ of spins on
the sublattice of even sites satisfy

$i\partial_{Jt}\ket{\varphi_{n}(t)}=K_{n-1,n}\ket{\varphi_{n}(t)}+\\\
+(\mathbbm{1}-\sigma^{z}_{n-1})\ket{\varphi_{n-1}(t)}+(\mathbbm{1}-\sigma^{z}_{n})\ket{\varphi_{n+1}(t)}\,$
(10)

with $\ket{\varphi_{n}(0)}=\delta_{n,0}\ket{\varphi_{0}}$. Here
$K_{n-1,n}=\sigma_{n-1}^{x}\sigma_{n}^{x}+\sigma_{n-1}^{y}\sigma_{n}^{y}$.
Note that this is a system of $\log_{2}\sqrt{D}$ equations for
$\ket{\varphi_{n}(t)}$, which are states belonging to subspaces of size
$\sqrt{D}$, where $D$ denotes the size of the Hilbert space.

The dynamical equation (10) is block-diagonal in the Fourier space defined by

$\ket{\tilde{\varphi}_{P}(t)}=\sum_{n}e^{inP}\Pi^{-n}\ket{\varphi_{n}(t)}\,,$
(11)

where $\Pi$ is the $1$-site shift operator on the sublattice, such that
$\Pi^{-1}=\Pi^{\dagger}$ and
$\Pi(O_{\ell})_{o(e)}\Pi^{\dagger}=(O_{\ell+1})_{o(e)}$. Specifically, we find

$\ket{\tilde{\varphi}_{P}(t)}=e^{-i\tilde{H}(P)t}\ket{\varphi_{0}}\,,$ (12)

where the independent time evolutions labelled by $P$ are generated by

$\tilde{H}(P)=J\left[K_{-1,0}+(e^{-iP}\Pi(\mathbbm{1}-\sigma^{z}_{-1})+h.c.)\right]\,.$
(13)

Since $\Pi$ is a sublattice shift, $P$ is the eigenvalue of the momentum
generating the two-site shifts on the full lattice. Indeed, denoting by $U$
the map $\ket{\varphi_{n}(t)}\mapsto\Pi\ket{\varphi_{n-1}(t)}$ we find
$U\ket{\tilde{\varphi}_{P}(t)}=e^{iP}\ket{\tilde{\varphi}_{P}(t)}$.

If the sites $-1$ and $0$ of the sublattice are both occupied by a spin up,
the state is destroyed by $\tilde{H}(P)$: a quarter of the Hilbert space is a
nullspace of $\tilde{H}(P)$. The nontrivial action of $\tilde{H}(P)$ in the
remaining space corresponds to moving a configuration of $N$ spins up on a
lattice of $L=\log_{2}\sqrt{D}$ sites with periodic boundary conditions
through a two-site defect at $-1$ and $0$. If both spins at $-1$ and $0$ are
down, $\tilde{H}(P)$ acts as a right or left global shift of the spins up by
one site. If instead there is a single spin up at $-1$ or $0$, either that
spin is moved to the other site or all spins up are globally shifted in the
opposite direction. Importantly, once all spins up have passed across the
defect, the relative distances between them are the same as in the initial
configuration. This key observation allows us to represent each configuration
as an effective particle propagating by a hopping Hamiltonian (with localised
defects that can be seen as a deformation of the space), on an extended
lattice of $L+N-1$ sites. The mapping is described in the SM [48]. We find the
energies to be parametrised as $E_{P}(k)=4J\cos(k+P)$, where the momentum $k$
of the effective particle and $P$ satisfy

$e^{i(L+N-1)k+iNP}=1,\qquad e^{iLP}=1\,.$ (14)

This mapping can also be exploited to compute the overlaps between the states
and the matrix elements of the spin operators [48], providing therefore all
the ingredients for the exact computation of the scaling profiles, such as the
ones depicted in Figs 2 and 3.

For example, for odd $x$ (for the general expression see SM [48]) the jamming
condition (2) after the local spin flip can be written as (after the
projective measurement there is instead an additional overall factor $1/3$)

$\displaystyle\begin{gathered}{\cal
P}_{\downarrow\downarrow}(x,t)=\frac{1}{L^{2}}\sum_{N}\sum_{\chi^{(1)}_{N},\chi^{(2)}_{N}}\braket{\varphi_{0}}{\chi_{N}^{(1)}}\braket{\chi_{N}^{(2)}}{\varphi_{0}}\times\\\
\braket{\chi_{N}^{(1)}}{\downarrow}\\!\braket{\downarrow}{{}_{0}\,\chi_{N}^{(2)}}e^{i\lceil\frac{x}{2}\rceil(P_{1}-P_{2})+it(E_{P_{1}}(k_{1})-E_{P_{2}}(k_{2}))}\,.\end{gathered}$
(17)

Here, $\chi_{N}=(k;P,\ldots)$ denotes the eigenstates and $\ldots$ stands for
additional quantum numbers characterising the exponentially large sectors at
fixed $k$ and $P$. The sum excludes the jammed states, since the latter are
destroyed by the observable. It turns out that the sum over the additional
quantum numbers can be carried out analytically, and one can reduce the entire
expression to a finite number of sums, which will be discussed in a more
technical work still in preparation [50]. We only anticipate that a thorough
analysis of (17) shows that, first, the terms with $N/L\sim 1/2$ contribute
the most and, second, the asymptotic behaviour is determined by some singular
points of the averaged matrix elements of the observable (the projector on two
neighbouring spins down).

### Numerical simulations.

Our effective sublattice description of time evolution also has some numerical
advantages. First, it almost doubles the system size accessible to exact
diagonalisation techniques. In particular, systems with $26$ spins can be
simulated in a reasonable time with an ordinary laptop. As a matter of fact,
it is also possible to perform semi-analytical calculations after having
reduced expressions like (17) to finite numbers of sums. In that way it is
possible to reach total system sizes of up to $\sim 60$ sites ($L\sim 30$) in
few days of single-processor computational time (the numerical effort is
expected to scale as $L^{5.5}-L^{6}$). Our main numerical checks have been
however based on DMRG algorithms, which are much more efficient.

### Profiles.

We observe that, at large times, the numerical data are in excellent agreement
with two ballistic-scale conjectures (i.e., $f(x,t)\to f(\zeta=\tfrac{x}{t})$)

$t\,\mathcal{P}_{\downarrow\downarrow}(\zeta)\approx\frac{a}{[1-(\tfrac{\zeta}{v_{\pm}})^{2}]^{\frac{1}{2}}},\
\ \langle\sigma^{y}\rangle_{\zeta}\approx
b[1-(\tfrac{\zeta}{v_{\pm}})^{2}]^{\frac{1}{2}}\,,$ (18)

with $v_{\pm}=\pm 16J/3$, and whose assessment of validity is still in
progress [50]. The constant prefactors in the scaling functions are well
approximated by $a\approx 0.06$ and $b\approx 0.24$. These conjectures are
consistent with the prediction for the maximal velocity of quasiparticle
excitations on top of a jammed state $v_{\pm}[\ket{\Phi}]=\pm
8J/\left(1+\tfrac{2}{L}\braket{\Phi}{S^{z}}{\Phi}\right)$, which was derived
in Ref. [51].

### Can Alice recover the local state?

After a time $t$, the reduced density matrix of Alice’s site, assumed to be at
an odd distance $x+1$ from Bob’s past blind measurement at site $-1$, reads

$\rho_{\rm Alice}=(\mathbbm{1}+\vec{p}_{x,t}\cdot\vec{\sigma})/2\,,$ (19)

with $\vec{p}_{x,t}=\tfrac{1}{3}(\braket{\sigma^{x}}_{\zeta}-2\cos\tfrac{\pi
x}{n},\braket{\sigma^{y}}_{\zeta}+2\sin\tfrac{\pi x}{n},0)$. Here $\zeta=x/t$,
and the expectation values refer to the dynamics after the local spin flip –
see Fig. 3 and SM [48] for details. For $\zeta$ in the light-cone $\vec{p}$ is
slightly tilted from the initial (pre-measurement) direction, the tilt itself
depending on Bob’s position. As far as we can see, this makes it impossible
for Alice to fully recover the information she had without Bob’s measurement.

### What is general.

We argue that a key role in the phenomenon is played by the fact that the
movement of excitations over jammed states, i.e., the jamming-breaking
impurities, is associated with a few-sites shift of a _string_ of jammed
particles (see also Ref. [49]). In this way the state can store memory of its
past.

The effect seems to remain stable even under generic Hamiltonian perturbations
preserving the jammed states. We have checked it for a perturbation of the
form
$V=g\sum_{\ell}\frac{\mathbbm{1}-\sigma_{\ell}^{z}}{2}\frac{\mathbbm{1}-\sigma_{\ell+1}^{z}}{2}K_{\ell+2,\ell+3}$.
Flipping a single spin again results in macroscopic reorganisation of the spin
profiles, albeit it now takes place in diffusive-scale coordinates
$\ell\sim\sqrt{t}$ (see SM [48] for numerical evidence). Thus, the phenomenon
could be observed also in non-integrable systems, in the appropriate scaling
limit.

### Summary.

We addressed the question of how a local measurement in a jammed state of an
interacting quantum system affects the late-time dynamics of distant local
observables. The measurement triggers a ballistic dynamics in which jamming is
locally restored but the profile of observables remains irremediably affected
by the perturbation. Arguably, this prevents the full recovery of locally
damaged information in the wing-flap protocol, which is instead expected quite
generally [39].

Our work gives rise to several questions. Firstly, the precise role of
symmetry breaking for the observed phenomenon is yet to be clarified. Our
initial states always break some symmetries (e.g., $U(1)$, $\mathbbm{Z}_{2}$,
or two-site shift invariance – see also Ref. [49]), however, a similar
behaviour has been observed (after the completion of this paper) in a
different setting where no symmetry seems to be broken [52].

The second question is that of repeated projective measurements, recently
identified as the cause of dynamical phase transition in certain quantum many-
body systems [53, 54, 55, 56]. Their effects on the observables, as well as on
the restoration of jamming pose an intriguing open problem that could be
investigated in the dual folded XXZ model even analytically.

### Acknowledgements.

We thank Viktor Eisler for useful discussions. This work was supported by the
European Research Council under the Starting Grant No. 805252 LoCoMacro.

During the peer-review process of this manuscript, Ref. [52] provided a new
example of a setting where a localised perturbation remains relevant at late
times.

## References

  * Jaksch and Zoller [2005] D. Jaksch and P. Zoller, The cold atom hubbard toolbox, Annals of Physics 315, 52 (2005), special Issue.
  * Bloch  _et al._ [2012] I. Bloch, J. Dalibard, and S. Nascimbène, Quantum simulations with ultracold quantum gases, Nature Physics 8, 267 (2012).
  * Barreiro  _et al._ [2011] J. T. Barreiro _et al._ , An open-system quantum simulator with trapped ions, Nature 470, 486 (2011).
  * Tan  _et al._ [2015] T. R. Tan _et al._ , Multi-element logic gates for trapped-ion qubits, Nature 528, 380 (2015).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Probing many-body dynamics on a 51-atom quantum simulator, Nature 551, 579 (2017).
  * Zhang  _et al._ [2017] J. Zhang _et al._ , Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator, Nature 551, 601 (2017).
  * Arute  _et al._ [2019] F. Arute _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator, Nature 595, 227 (2021).
  * Buluta and Nori [2009] I. Buluta and F. Nori, Quantum simulators, Science 326, 108 (2009).
  * Schlosshauer [2019] M. Schlosshauer, Quantum decoherence, Physics Reports 831, 1 (2019), quantum decoherence.
  * Gritsev and Polkovnikov [2017] V. Gritsev and A. Polkovnikov, Integrable Floquet dynamics, SciPost Phys. 2, 021 (2017).
  * Vanicat  _et al._ [2018] M. Vanicat, L. Zadnik, and T. Prosen, Integrable trotterization: Local conservation laws and boundary driving, Phys. Rev. Lett. 121, 030606 (2018).
  * Ljubotina  _et al._ [2019] M. Ljubotina, L. Zadnik, and T. Prosen, Ballistic spin transport in a periodically driven integrable quantum system, Phys. Rev. Lett. 122, 150605 (2019).
  * Sá  _et al._ [2021] L. Sá, P. Ribeiro, and T. Prosen, Integrable nonunitary open quantum circuits, Phys. Rev. B 103, 115132 (2021).
  * Bertini  _et al._ [2018] B. Bertini, P. Kos, and T. Prosen, Exact spectral form factor in a minimal model of many-body quantum chaos, Phys. Rev. Lett. 121, 264101 (2018).
  * Nahum  _et al._ [2017] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Quantum entanglement growth under random unitary dynamics, Phys. Rev. X 7, 031016 (2017).
  * Chan  _et al._ [2018] A. Chan, A. De Luca, and J. T. Chalker, Solution of a minimal model for many-body quantum chaos, Phys. Rev. X 8, 041019 (2018).
  * von Keyserlingk  _et al._ [2018] C. W. von Keyserlingk, T. Rakovszky, F. Pollmann, and S. L. Sondhi, Operator hydrodynamics, otocs, and entanglement growth in systems without conservation laws, Phys. Rev. X 8, 021013 (2018).
  * Khemani  _et al._ [2018] V. Khemani, A. Vishwanath, and D. A. Huse, Operator spreading and the emergence of dissipative hydrodynamics under unitary evolution with conservation laws, Phys. Rev. X 8, 031057 (2018).
  * Bensa and Žnidarič [2021] J. Bensa and M. Žnidarič, Fastest local entanglement scrambler, multistage thermalization, and a non-hermitian phantom, Phys. Rev. X 11, 031019 (2021).
  * Aleiner [2021] I. L. Aleiner, Bethe ansatz solutions for certain periodic quantum circuits (2021), arXiv:2107.05715 [cond-mat.mes-hall] .
  * Claeys  _et al._ [2021] P. W. Claeys, J. Herzog-Arbeitman, and A. Lamacraft, Correlations and commuting transfer matrices in integrable unitary circuits (2021), arXiv:2106.00640 [quant-ph] .
  * Gopalakrishnan [2018] S. Gopalakrishnan, Operator growth and eigenstate entanglement in an interacting integrable floquet system, Phys. Rev. B 98, 060302 (2018).
  * Alba  _et al._ [2019] V. Alba, J. Dubail, and M. Medenjak, Operator entanglement in interacting integrable quantum systems: The case of the rule 54 chain, Phys. Rev. Lett. 122, 250603 (2019).
  * Klobas  _et al._ [2021] K. Klobas, B. Bertini, and L. Piroli, Exact thermalization dynamics in the “rule 54” quantum cellular automaton, Phys. Rev. Lett. 126, 160602 (2021).
  * Salathé  _et al._ [2015] Y. Salathé _et al._ , Digital quantum simulation of spin models with circuit quantum electrodynamics, Phys. Rev. X 5, 021027 (2015).
  * Neill  _et al._ [2021] C. Neill _et al._ , Accurately computing the electronic properties of a quantum ring, Nature 594, 508 (2021).
  * Calabrese  _et al._ [2016] P. Calabrese, F. H. L. Essler, and G. Mussardo, Special issue on quantum integrability in out of equilibrium systems, Journal of Statistical Mechanics: Theory and Experiment 2016, 064001 (2016).
  * Bertini and Fagotti [2016] B. Bertini and M. Fagotti, Determination of the nonequilibrium steady state emerging from a defect, Phys. Rev. Lett 117, 130402 (2016).
  * Castro-Alvaredo  _et al._ [2016] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Emergent hydrodynamics in integrable quantum systems out of equilibrium, Phys. Rev. X 6, 041065 (2016).
  * Bertini  _et al._ [2016] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Transport in out-of-equilibrium XXZ chains: Exact profiles of charges and currents, Phys. Rev. Lett. 117, 207201 (2016).
  * Piroli  _et al._ [2017] L. Piroli, J. De Nardis, M. Collura, B. Bertini, and M. Fagotti, Transport in out-of-equilibrium xxz chains: Nonballistic behavior and correlation functions, Phys. Rev. B 96, 115124 (2017).
  * Eisler and Maislinger [2020] V. Eisler and F. Maislinger, Front dynamics in the XY chain after local excitations, SciPost Phys. 8, 037 (2020).
  * Gruber and Eisler [2021] M. Gruber and V. Eisler, Entanglement spreading after local fermionic excitations in the XXZ chain, SciPost Phys. 10, 005 (2021).
  * Zauner  _et al._ [2015] V. Zauner, M. Ganahl, H. Evertz, and T. Nishino, Time evolution within a comoving window: Scaling of signal fronts and magnetization plateaus after a local quench in quantum spin chains, J. Phys.: Condens. Matter 27, 425602 (2015).
  * Eisler  _et al._ [2016] V. Eisler, F. Maislinger, and H. G. Evertz, Universal front propagation in the quantum Ising chain with domain-wall initial states, SciPost Phys. 1, 014 (2016).
  * Eisler and Maislinger [2018] V. Eisler and F. Maislinger, Hydrodynamical phase transition for domain-wall melting in the XY chain, Phys. Rev. B 98, 161117(R) (2018).
  * Campisi and Goold [2017] M. Campisi and J. Goold, Thermodynamics of quantum information scrambling, Phys. Rev. E 95, 062127 (2017).
  * Yan and Sinitsyn [2020] B. Yan and N. A. Sinitsyn, Recovery of damaged information and the out-of-time-ordered correlators, Phys. Rev. Lett. 125, 040605 (2020).
  * Gluza  _et al._ [2019] M. Gluza, J. Eisert, and T. Farrelly, Equilibration towards generalized gibbs ensembles in non-interacting theories, SciPost Phys. 7, 038 (2019).
  * Alba  _et al._ [2021] V. Alba, B. Bertini, M. Fagotti, L. Piroli, and P. Ruggiero, Generalized-Hydrodynamic approach to Inhomogeneous Quenches: Correlations, Entanglement and Quantum Effects, arXiv e-prints , arXiv:2104.00656 (2021), arXiv:2104.00656 [cond-mat.stat-mech] .
  * Zadnik and Fagotti [2021] L. Zadnik and M. Fagotti, The Folded Spin-1/2 XXZ Model: I. Diagonalisation, Jamming, and Ground State Properties, SciPost Phys. Core 4, 10 (2021).
  * Bariev [1991] R. Z. Bariev, Integrable spin chain with two- and three-particle interactions, Journal of Physics A: Mathematical and General 24, L549 (1991).
  * Pozsgay  _et al._ [2021] B. Pozsgay, T. Gombor, A. Hutsalyuk, Y. Jiang, L. Pristyák, and E. Vernier, An integrable spin chain with hilbert space fragmentation and solvable real time dynamics (2021), arXiv:2105.02252 .
  * Menon  _et al._ [1997] G. I. Menon, M. Barma, and D. Dhar, Conservation laws and integrability of a one-dimensional model of diffusing dimers, Journal of Statistical Physics 86, 1237 (1997).
  * Yang  _et al._ [2020] Z.-C. Yang, F. Liu, A. V. Gorshkov, and T. Iadecola, Hilbert-space fragmentation from strict confinement, Phys. Rev. Lett. 124, 207602 (2020).
  * Fagotti [2014] M. Fagotti, On conservation laws, relaxation and pre-relaxation after a quantum quench, Journal of Statistical Mechanics: Theory and Experiment 2014, P03016 (2014).
  * [48] Supplementary material.
  * Zadnik  _et al._ [2021a] L. Zadnik, S. Bocini, K. Bidzhiev, and M. Fagotti, Measurement catastrophe and ballistic spread of charge density with vanishing current (2021a), arXiv:2111.06325 [quant-ph] .
  * [50] L. Zadnik and M. Fagotti, In preparation.
  * Zadnik  _et al._ [2021b] L. Zadnik, K. Bidzhiev, and M. Fagotti, The Folded Spin-1/2 XXZ Model: II. Thermodynamics and Hydrodynamics with a Minimal Set of Charges, SciPost Phys. 10, 99 (2021b).
  * Fagotti [2021] M. Fagotti, Global quenches after localised perturbations (2021), arXiv:2110.11322 [cond-mat.str-el] .
  * Li  _et al._ [2018] Y. Li, X. Chen, and M. P. A. Fisher, Quantum zeno effect and the many-body entanglement transition, Phys. Rev. B 98, 205136 (2018).
  * Chan  _et al._ [2019] A. Chan, R. M. Nandkishore, M. Pretko, and G. Smith, Unitary-projective entanglement dynamics, Phys. Rev. B 99, 224307 (2019).
  * Skinner  _et al._ [2019] B. Skinner, J. Ruhman, and A. Nahum, Measurement-induced phase transitions in the dynamics of entanglement, Phys. Rev. X 9, 031009 (2019).
  * Minato  _et al._ [2021] T. Minato, K. Sugimoto, T. Kuwahara, and K. Saito, Fate of measurement-induced phase transition in long-range interactions (2021), arXiv:2104.09118 [quant-ph] .

“Macroscopic effects of localised measurements in jammed states of quantum
spin chains”

Kemal Bidzhiev1, Maurizio Fagotti1, Lenart Zadnik

1Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France

##  .1 The effective dynamics generated by $\tilde{H}(P)$

Here we describe the idea behind the diagonalisation of the effective
Hamiltonian $\tilde{H}(P)$ from Eq. (12) of the main text. Detailed method
will be presented in a separate publication [50]. Consider a configuration of
$N$ spins up; the effective Hamiltonian

$\begin{gathered}\tilde{H}(P)=JK_{-1,0}+Je^{-iP}\Pi(\mathbbm{1}-\sigma^{z}_{-1})+\\\
+Je^{iP}\Pi^{-1}(\mathbbm{1}-\sigma^{z}_{0})\end{gathered}$ (S.1)

maps it into a superposition of identical configurations shifted by one site
to the right (left), with a phase $e^{-iP}$ (resp. $e^{iP}$). Repeated action
of $\tilde{H}(P)$ propagates configuration in this manner until a spin up
comes onto position $-1$ (or $0$). The term $JK_{-1,0}$ then kicks in: the
spin up jumps to position $0$ (resp. $-1$) without acquiring a phase, while
the rest of spins up retain their positions. Here is an example of nontrivial
microscopic dynamics considering only the part of a superposition that moves
in one direction:

$t$-$5$-$4$-$3$-$2$-$1$$0$$1$$2$$3$$4$$5$black$\uparrow$black$\uparrow$black$\uparrow$black$\uparrow$black$\uparrow$lightred$\uparrow$black$\uparrow$black$\uparrow$lightred$\uparrow$black$\uparrow$lightred$\uparrow$black$\uparrow$black$\uparrow$lightred$\uparrow$black$\uparrow$black$\uparrow$black$\uparrow$black$\uparrow$lightred$\uparrow$black$\uparrow$black$\uparrow$lightred$\uparrow$black$\uparrow$black$\uparrow$$2J$$2Je^{iP}$

For the coordinate of the effective particle (in blue) we have chosen a
shifted centre-of-mass position

$X(\underline{\ell})=\frac{1}{N}\sum_{j=1}^{N}\ell_{j}+\frac{N-1}{2N}(N-2N_{-}).$
(S.4)

Here $\ell_{j}$ are positions of spins up, $N$ their total number, and $N_{-}$
denotes the number of cases $\ell_{j}<0$. The second term in (S.4) ensures
that $X(\underline{\ell})$ changes by $1$ even if only one of the spins up
jumps between the sites $-1$ and $0$, while the rest of spins up retain their
positions (the corresponding hopping of the effective particle is represented
by the red curved arrows in the above diagram). The system is thus described
by a hopping Hamiltonian on the lattice
$\Lambda=X[\underline{\ell}]+\mathbbm{Z}$ with defects accounting for the
change in the coupling constants whenever a spin up jumps between the sites
$-1$ and $0$ (in red).

The defects can be removed via a unitary transformation, which leads to a
hopping model with single-particle energy levels

$E_{P}(k)=4J\cos(k+P)\,,$ (S.5)

where $k$ denotes the momentum of the effective particle. The eigenstates of
the effective Hamiltonian (S.1) are now parametrised as
$\ket{k;P,\underline{d}}$, where $\underline{d}$ is a collection of distances
between subsequent spins up, defined so as to be preserved by the dynamics:

$d_{j}=\ell_{j+1}-\ell_{j}-1+\theta(\ell_{j+1}<0)-\theta(\ell_{j}<0).$ (S.6)

For a finite sublattice size $L=\log_{2}\sqrt{D}$ one needs to consider
periodic boundary conditions, where $k$ and $P$ are instead quantised
according to

$e^{iLP}=1,\qquad e^{i(L+N-1)k+iNP}=1.$ (S.7)

Quantisation of $P$ follows from the momentum $P$ being associated with a
shift operator $\Pi$ on the sublattice. The factor $e^{iNP}$ in the
quantisation condition for $k$ is a result of the unitary transformation that
removes the $N$ defects and effectively adds a phase to each hopping process
caused by the defect. Factor $e^{i(L+N-1)k}$ comes from the lag acquired by
the configuration w.r.t. to the position of the effective particle, when the
latter traverses the system. Each time a spin up jumps from $0$ to $-1$, the
right-most spin up remains fixed, except if it itself jumps between these two
sites. The total acquired lag is therefore $N-1$ sites, so the configuration
needs $L+N-1$ steps before returning to its initial position.

As a proof of concept, Fig. S1(a) shows the comparison between the DMRG-based
algorithm and exact-diagonalisation of Hamiltonians $\tilde{H}(P)$.

Figure S1: Panel (a) shows the profile
$t\mathcal{P}_{\downarrow\downarrow}(x,t)$ at $t=0.8$, for $L=13$ (a chain of
$26$ spins), computed with a DMRG-based algorithm and exact diagonalisation of
the effective Hamiltonians (S.1). Panel (b) shows contributions to
$\mathcal{P}_{\downarrow\downarrow}$ coming from summing the elements with
given $N$ (and hence $\braket{S^{z}}$) at $x=0$ and $t=1/2$ for the sublattice
sizes $5\leq L\leq 13$. Data are fitted with Gaussian curves centered around
the sublattice magnetisation $2\langle S^{z}\rangle=2N-L=-1$. The full width
at half maximum scales sublinearly with the sublattice size. The term with the
highest contribution is $N=(L-1)/2$.

##  .2 Jamming condition $\mathcal{P}_{\downarrow\downarrow}(x,t)$

In this section we describe the ingredients for the computation of the local
jamming condition ${\cal
P}_{\downarrow\downarrow}(x,t)=\tfrac{1}{4}\braket{\Psi(t)}{(\mathbbm{1}-\sigma^{z}_{x})(\mathbbm{1}-\sigma^{z}_{x+1})}{\Psi(t)}$
in Eq. (14) of the main text. Representation

$\ket{\Psi(t)}=\sum_{n}(\sigma^{x}_{n}\ket{\uparrow\ldots\uparrow})_{o}\otimes\ket{\varphi_{n}(t)}_{e}$
(S.8)

of the wave function is useful to reduce the computation to the sublattice of
even sites, where the condition becomes

${\cal
P}_{\downarrow\downarrow}=\frac{1}{2}\begin{cases}\braket{\varphi_{\frac{x}{2}+1}(t)}{\mathbbm{1}-\sigma^{z}_{\frac{x}{2}}}{\varphi_{\frac{x}{2}+1}(t)}&\text{$x$
even}\,,\\\
\braket{\varphi_{\frac{x+1}{2}}(t)}{\mathbbm{1}-\sigma^{z}_{\frac{x+1}{2}}}{\varphi_{\frac{x+1}{2}}(t)}&\text{$x$
odd}\,.\end{cases}$ (S.9)

The operators on the right-hand side act on the state of the sublattice –
their indices have been changed accordingly. Using Eqs (10), (11) of the main
text, we then obtain

$\displaystyle\begin{gathered}{\cal
P}_{\downarrow\downarrow}=\frac{1}{L^{2}}\sum_{N}\sum_{\underline{d}_{N}^{(1)},\underline{d}_{N}^{(2)}}\sum_{P_{1},P_{2}\atop
e^{iLP}=1}\sum_{k_{1},k_{2}\atop
e^{i(L+N-1)k+iNP}=1}e^{i\big{(}\frac{x}{2}+\frac{3+(-1)^{x}}{4}\big{)}(P_{1}-P_{2})}e^{i4Jt(\cos(k_{1}+P_{1})-\cos(k_{2}+P_{2}))}\times\\\
\times\braket{\varphi_{0}}{k_{1};P_{1},\underline{d}_{N}^{(1)}}\braket{k_{1};P_{1},\underline{d}_{N}^{(1)}}{\frac{\mathbbm{1}-\sigma^{z}_{-[1+(-1)^{x}]/2}}{2}}{k_{2};P_{2},\underline{d}_{N}^{(2)}}\braket{k_{2};P_{2},\underline{d}_{N}^{(2)}}{\varphi_{0}}\,.\end{gathered}$
(S.12)

The matrix element and the overlaps between the wave functions can be computed
exactly by employing the mapping between a configuration of spins up and an
effective particle [50]. Here we only report the resulting formulas for the
case $\ket{\varphi_{0}}=\ket{\leftarrow\ldots\leftarrow}$ (i.e. after a local
measurement in the state $\ket{\mathcal{L}_{1}}$). For $k_{1}\neq k_{2}$ they
read

$\displaystyle\begin{aligned}
\braket{\leftarrow\ldots\leftarrow}{k;P,\underline{d}_{N}}&=\frac{1-e^{iP}}{e^{ik}-1}\sum_{\ell=1}^{N}\frac{e^{i(2\ell-1-2N+\sum_{n=N+1-\ell}^{N-1}d_{n})k}e^{i(\ell-1)P}}{\sqrt{2^{L}(L+N-1)}}\,,\\\
\braket{k_{1};P_{1},\underline{d}_{N}}{\tfrac{\mathbbm{1}-\sigma_{0}^{z}}{2}}{k_{2};P_{2},\underline{d}_{N}}&=\frac{1-e^{i(P_{2}-P_{1})}e^{i(k_{2}-k_{1})}}{e^{i(k_{2}-k_{1})}-1}\sum_{\ell=0}^{N-1}\frac{e^{i(2\ell+1-2N+\sum_{n=N-\ell}^{N-1}d_{n})(k_{2}-k_{1})+i\ell(P_{2}-P_{1})}}{L+N-1}\,,\\\
\braket{k_{1};P_{1},\underline{d}_{N}}{\tfrac{\mathbbm{1}-\sigma_{-1}^{z}}{2}}{k_{2};P_{2},\underline{d}_{N}}&=e^{i(k_{1}-k_{2})}\braket{k_{1};P_{1}\underline{d}_{N}}{\tfrac{\mathbbm{1}-\sigma_{0}^{z}}{2}}{k_{2};P_{2},\underline{d}_{N}}\,,\end{aligned}$
(S.13)

whereas the limit $k_{1}=k_{2}$ has to be performed carefully, taking the
quantisation conditions (S.7) into account. We point out that the numerical
computation of the fixed-$N$ terms in (S.12) for small system sizes (via exact
diagonalisation) suggests that the main contribution to the jamming condition
${\cal P}_{\downarrow\downarrow}(x,t)$ at fixed $x$ and $t$ comes from the
terms with $N\sim L/2$ – see Fig. S1(b). This numerical observation can be
proved analytically and holds for any sublattice size $L$ [50].

##  .3 Reduced density matrices

The purpose of this section is to explain Alice’s and Bob’s reduced density
matrices. Bob performs a blind projective measurement on one of the spins up
in a factorised state. Suppose that the axis of the measurement is fixed. It
can be obtained by rotating the $z$-axis by an angle $\phi\in[0,\pi)$ around
the $xy$-plane unit vector $\hat{n}_{\varphi}=(\cos\varphi,\sin\varphi,0)$,
parametrised by $\varphi\in[0,2\pi)$. Had Bob read off the result of the
measurement, the spin would have collapsed into the state

$\rho(\Phi,\varphi)=e^{i\tfrac{\phi}{2}\hat{n}_{\varphi}\cdot\vec{\sigma}}\ket{\uparrow}\bra{\uparrow}e^{-i\tfrac{\phi}{2}\hat{n}_{\varphi}\cdot\vec{\sigma}}\,.$
(S.14)

In our protocol, described in the main text, the measurement is instead blind:
Bob does not read off its result. Hence, due to a classical uncertainty, the
state after the measurement is computed as an average w.r.t. the Haar measure:

$\displaystyle\begin{gathered}\rho_{\rm Bob}=\int_{0}^{2\pi}\frac{{\rm
d}\varphi}{2\pi}\int_{0}^{\pi}\frac{{\rm
d}\phi}{\pi}\rho(\phi,\varphi)\left(\frac{\pi}{2}\cos\frac{\phi}{2}\right)=\\\
=\frac{2}{3}\ket{\uparrow}\bra{\uparrow}+\frac{1}{3}\ket{\downarrow}\bra{\downarrow}\,.\end{gathered}$
(S.17)

The weight $\tfrac{\pi}{2}\cos\tfrac{\phi}{2}$ in the integral is the
probability for the nonzero projection of the spin up onto the measurement
axis.

At time $t$ after Bob’s blind projective measurement at site $-1$, the density
matrix of the system reads

$\displaystyle\rho(t)=\frac{2}{3}\ket{\Phi}\bra{\Phi}+\frac{1}{3}\ket{\Psi(t)}\bra{\Psi(t)},$
(S.18)

where $\ket{\Psi(0)}=\sigma^{x}_{-1}\ket{\Phi}$, and
$\ket{\Phi}=\ket{\mathcal{L}_{n}}$ is a jammed state given in Eq. (4) of the
main text. Alice’s reduced density matrix at site $x$ (at an odd distance
$x+1$ from Bob’s measurement) in general reads

$\displaystyle\rho_{\text{Alice}}=(\mathbbm{1}+\vec{p}_{x,t}\cdot\vec{\sigma})/2\,,$
(S.19)

where

$\displaystyle\vec{p}_{x,t}=({\rm Tr}[\rho(t)\sigma^{x}_{x}],{\rm
Tr}[\rho(t)\sigma^{y}_{x}],{\rm Tr}[\rho(t)\sigma^{z}_{x}])\,.$ (S.20)

In the jammed state $\ket{\Phi}=\ket{\mathcal{L}_{n}}$ Alice finds (note that
$x$ is even)

$\displaystyle\begin{aligned}
\bra{\Phi}\sigma_{x}^{x}\ket{\Phi}&=-\cos\tfrac{\pi x}{n}\,,\\\
\bra{\Phi}\sigma_{x}^{y}\ket{\Phi}&=\sin\tfrac{\pi x}{n}\,,\\\
\bra{\Phi}\sigma_{x}^{z}\ket{\Phi}&=0\,.\end{aligned}$ (S.21)

On the other hand, in the time-evolved part of the state, $\ket{\Psi(t)}$, the
expectation values of Pauli matrices asymptotically depend only on the ray
$\zeta=x/t$ (we assume large distance $x$ and time $t$), i.e.,
$\bra{\Psi(t)}\sigma_{x}^{\alpha}\ket{\Psi(t)}=\braket{\sigma^{\alpha}}_{\zeta}$,
$\alpha\in\\{x,y,z\\}$. In particular, the $z$-component of spin at an even
site far from Bob’s measurement is asymptotically zero, i.e.,
$\braket{\sigma^{z}}_{\zeta}=0$ (see Fig. 3 of the main text). States
$\ket{\Phi}$ and $\ket{\Psi(t)}$ enter the classical mixture (S.18) with
probabilities $2/3$ and $1/3$, respectively, whence we finally obtain

$\displaystyle\vec{p}_{x,t}=\tfrac{1}{3}\big{(}\braket{\sigma^{x}}_{\zeta}-2\cos\tfrac{\pi
x}{n},\braket{\sigma^{y}}_{\zeta}+2\sin\tfrac{\pi x}{n},0\big{)}\,.$ (S.22)

##  .4 Perturbation preserving the jammed states

In this section we report the effects of a perturbation that breaks
integrability while preserving the jammed sector. Specifically, we consider
the Hamiltonian

$\displaystyle
H^{\prime}=H+g\sum_{\ell}\frac{\mathbbm{1}-\sigma_{\ell}^{z}}{2}\frac{\mathbbm{1}-\sigma_{\ell+1}^{z}}{2}K_{\ell+2,\ell+3}\,.$
(S.23)

Fig. S2 shows $\braket{\sigma^{x}}$ after that a spin in a jammed state of the
perturbed Hamiltonian $H^{\prime}$ is flipped. The expectation value exhibits
diffusive scaling. This allows for a description of the local observables by
macrostates depending on the diffusive-scale coordinate $\ell/\sqrt{t}$.

Figure S2: The profile $\braket{\sigma_{\ell}^{x}}_{t}$ in the diffusive-scale
coordinates $\ell/\sqrt{t}$, at times $t=4$ (blue) and $8$ (red) after
flipping up the central spin (in an odd position) of the state
$\ket{\mathcal{L}_{4}}=\ket{\uparrow\leftarrow\uparrow\rightarrow\cdots}$. The
initial state evolves under $H^{\prime}$ with $g=\sqrt{2}$. Nontrivial
profiles emerge on _even_ sites, while spins on odd positions (not shown)
relax to their original values along the $z$-axis. The inset shows the jamming
condition $\sqrt{t}\mathcal{P}_{\downarrow\downarrow}(\ell,t)$ for different
times in the diffusive-scale coordinates:
$\mathcal{P}_{\downarrow\downarrow}(\ell,t)$ vanishes with time, implying the
emergence of LQJS on diffusive scales. The decay of the jamming condition (up
to DMRG-accessible times) seems to point at the scaling $\sim 1/\sqrt{t}$, but
longer times would be required for a conclusive answer.