# Tensorial tomographic Fourier Ptychography with applications to muscle
tissue imaging

Shiqi Xu  Duke University, Durham, NC, USA 27708  Xiang Dai  Duke University,
Durham, NC, USA 27708  UC San Diego, La Jolla, CA, USA 92093  Paul Ritter
Friedrich-Alexander University, Erlangen, Germany 91052  Kyung Chul Lee  Duke
University, Durham, NC, USA 27708  Yonsei University, Seoul, South Korea 03722
Xi Yang  Duke University, Durham, NC, USA 27708  Lucas Kreiss  Duke
University, Durham, NC, USA 27708  Friedrich-Alexander University, Erlangen,
Germany 91052  Kevin C. Zhou  Duke University, Durham, NC, USA 27708  UC
Berkeley, Berkeley, CA, USA 94720  Kanghyun Kim  Duke University, Durham, NC,
USA 27708  Amey Chaware  Duke University, Durham, NC, USA 27708  Jadee Neff
Duke Hospital, Durham, NC, USA 27710  Carolyn Glass  Duke Hospital, Durham,
NC, USA 27710  Seung Ah Lee  Yonsei University, Seoul, South Korea 03722
Oliver Friedrich  Friedrich-Alexander University, Erlangen, Germany 91052
Roarke Horstmeyer  Duke University, Durham, NC, USA 27708  Corresponding
author<EMAIL_ADDRESS>
###### Abstract

We report Tensorial tomographic Fourier Ptychography (_T 2oFu_), a new non-
scanning label-free tomographic microscopy method for simultaneous imaging of
quantitative phase and anisotropic specimen information in 3D. Built upon
Fourier Ptychography, a quantitative phase imaging technique, _T 2oFu_
additionally highlights the vectorial nature of light. The imaging setup
consists of a standard microscope equipped with an LED matrix, a polarization
generator, and a polarization-sensitive camera. Permittivity tensors of
anisotropic samples are computationally recovered from polarized intensity
measurements across three dimensions. We demonstrate _T 2oFu_’s efficiency
through volumetric reconstructions of refractive index, birefringence, and
orientation for various validation samples, as well as tissue samples from
muscle fibers and diseased heart tissue. Our reconstructions of muscle fibers
resolve their 3D fine-filament structure and yield consistent morphological
measurements compared to gold-standard second harmonic generation scanning
confocal microscope images found in the literature. Additionally, we
demonstrate reconstructions of a heart tissue sample that carries important
polarization information for detecting cardiac amyloidosis.

_K_ eywords Computational imaging, 3D imaging, phase retrieval microscopy,
polarization-sensitive imaging, label-free imaging

##  1 Introduction

Quantitative phase imaging (QPI) is a well-known label-free microscopy
approach that can detect phase delay introduced by semi-transparent cells and
tissue [1]. Due in part to its ability to provide quantitative information
about primarily transparent biological specimens with low phototoxicity, QPI
has become an invaluable tool in scientific and clinical studies, including
for monitoring neuronal firing [2] and cancer cell line detection [3], to name
a few. Besides scalar phase contrast, transparent specimens also exhibit
alternative and important endogenous optical contrast mechanisms, including
anisotropic properties such as material birefringence and orientation. Indeed,
the orientation of molecular arrangements in lipid membranes can now be
quantitatively monitored with polarization-sensitive microscopes for studying
multi-organelle interactive activities [4, 5]. There are naturally a variety
of polarization-sensitive microscope arrangements, including early analog
designs [6, 7], differential interference contrast (DIC) methods [8], and
contemporary digital approaches that reconstruct quantitative specimen
retardance and orientation [9, 10, 11]. These methods have been applied to
study small model organisms [12, 13] and to assist clinical diagnosis [14, 15,
16, 17].

Recently, there has been increased interest in measuring polarization-
sensitive phase information from specimens. In general, polarization-sensitive
quantitative phase imaging (PS-QPI) methods can be divided into two
categories: those that rely on interferometric detection, and those that
utilize computational phase retrieval methods. Interferometric methods (e.g.,
off-axis holography) can create polarization-sensitive phase images with as
few as a single measurement [18, 19, 20, 21], and can be extended to 3D with
diffraction tomography approaches [22, 23, 24, 25]. While often impressive,
these methods usually require complex arrangements of coherent laser
illumination and careful system alignment, which sets additional design
requirements for use in clinical applications. Computational phase retrieval
methods, on the other hand, rely on image reconstruction algorithms that
convert multiple intensity measurements into phase-sensitive outputs [26, 27,
28, 29, 30, 31], and can be implemented with less expensive hardware. Due to
their simple instrumentation, computational polarization microscopes have been
increasingly applied to image biological samples, for instance, to study white
matter tracts within whole brain slices [30] and to diagnose malaria from
blood smears [28]. Furthermore, these methods have been extended to image 3D
samples, such as axons and cardiac tissue with axial scanning [32, 33]. Thick
samples such as organoids and tissue slices also naturally have intriguing 3D
structures that can include anisotropic material. Creating high-resolution
volumetric representations of the polarization properties of these samples is
essential to studying biology and pathology. While tomographic imaging methods
such as confocal-based approaches have been developed in the past to image
muscle tissue and neural organoids, for instance [34, 35], there remain
relatively few microscopic techniques to jointly capture quantitative phase
and anisotropy across a large three-dimensional volume at high resolution. One
very recent study attempted to create 3D anisotropy maps using off-axis LED
illumination but did not provide tomographic permittivity matrix
reconstructions [36].

Here, we propose a non-scanning polarization-sensitive tomography method,
termed Tensorial tomographic Fourier Ptychography (_T 2oFu_), to create
quantitative volumetric permittivity matrix reconstructions without any moving
parts. Our method is an extension of recently developed intensity optical
diffraction tomography principles [37, 38, 39, 40, 41, 42] and has the
potential to be extended to a video-rate system in the future [37]. We image a
variety of calibration targets, as well as 3D tissue and muscle fibers
samples, through which we demonstrate the ability to resolve the fine filament
structure and are consistent with one of the gold standards, i.e., second
harmonic images, found in the literature [34]. Additionally, we demonstrate
reconstructions of a heart tissue sample that carries important information
for detecting cardiac amyloidosis.

##  2 Methods

Figure 1: Workflow of the proposed method. (A) An illustration of _T 2oFu_’s
experimental setup. Green light from an addressable LED array is circularly
polarized with a left circular polarizer (LCP) to illuminate the sample.
Sequential illuminations from various angles are used to scan the object in
the spatial frequency domain. The sample is then imaged with an infinity-
corrected optical system. The polarized light intensity at
$0\degree,45\degree,90\degree$, and $135\degree$ are recorded with a
polarization-sensitive CMOS camera as images illustrated in (B). (C) Those
intensity measurements are then fused to form volumetric sample permittivity
tensor reconstructions by solving the corresponding inverse problem.

###  2.1 Experimental setup

The imaging setup of _T 2oFu_ is illustrated in Fig. 1(A). The illumination
consists of an array of $25\times 25$ addressable LEDs (µInventions Inc;
Portugal) and a left-circular polarizer (CP42HE; Edmund Co., USA).
$521\,\text{nm}$ wavelength light from the LED array is circularly polarized
and illuminates the 3D sample from different angles. A microcontroller (ARM
Cortex-M3) and a voltage level shifter (SN74AHCT) are used to turn on each
small LED pixel (WS2812b-2020) sequentially. The optical field then passes
through an optical system consisting of an infinity-corrected objective (0.25
NA or 0.4NA, Plan Achromatic; Olympus, JP) and a tube lens (Achromatic
Doublets, 180mm focal length; Thorlabs, Inc, USA). The largest illumination NA
($NA_{illum}$) is chosen to match the native NA of the objective ($NA_{obj}$).
Hence, the theoretical lateral and axial resolutions are
$\delta_{x,y}=\nicefrac{{\lambda}}{{2(NA_{obj}+NA_{illum})}}$ and
$\delta_{z}=\nicefrac{{\lambda}}{{(2-\sqrt{1-{NA}_{obj}^{2}}-\sqrt{1-{NA}_{illum}^{2}})}}$,
respectively [43]. For the 0.25 NA system, $\delta_{x,y}=5.2\,\mu\mathrm{m}$
and $\delta_{z}=8.2\,\mu\mathrm{m}$. For the 0.4 NA system,
$\delta_{x,y}=3.2\,\mu\mathrm{m}$ and $\delta_{z}=3.1\,\mu\mathrm{m}$. The
intensity image is captured with a polarization-sensitive CMOS camera
(BFS-U3-51S5P; Teledyne FLIR LLC; USA). The polarization camera analyzes the
light orientated at four different polarizations
($0\degree,45\degree,90\degree,$ and $135\degree$), achieved by placing
$2\times 2$ wire-grid polarizers between the pixel array and the microlens
array, as illustrated in the m-dashed box in Fig. 1(A). Exemplary intensity
images under illumination from different angles are shown in the box in Fig.
1(B).

###  2.2 Principles of Tensorial tomographic Fourier Ptychography

####  2.2.1 Notation

Here we introduce the notation we use in the rest of this article. First, we
use $\vec{\cdot}$ and $\bar{\bar{\cdot}}$ symbols to denote vectors and
matrices for variables, respectively. By default, all the vectors in this
article are column vectors. Second, we utilize the Mathematical Script font
enclosed in curly braces to denote spatial operators. As an example, one
frequently used operator is $\mathcal{F}\\{\cdot\\}$, which indicates the
Fourier transform in space. Next, bold letters in lowercase represents support
vectors in either frequency ($\mathbf{u}$) or space ($\mathbf{r}$). Finally,
we utilize the $\tilde{\cdot}$ symbol to denote the frequency-domain
counterparts of variables previously defined in the space domain, such as
$\tilde{I}(\mathbf{u})=\mathcal{F}\\{I(\mathbf{r})\\}$.

####  2.2.2 Vectorial light propagation

The optical properties of a 3D sample can be described by its $3\times 3$
permittivity matrix [44]

$\bar{\bar{\epsilon}}=\begin{bmatrix}\epsilon_{xx}(\mathbf{r})&\epsilon_{xy}(\mathbf{r})&\epsilon_{xz}(\mathbf{r})\\\
\epsilon_{yx}(\mathbf{r})&\epsilon_{yy}(\mathbf{r})&\epsilon_{yz}(\mathbf{r})\\\
\epsilon_{zx}(\mathbf{r})&\epsilon_{zy}(\mathbf{r})&\epsilon_{zz}(\mathbf{r})\end{bmatrix},$
(1)

where $\mathbf{r}=(x,y,z)$ is the voxel position in space. In the scope of
this work, we adopt the first Born approximation [32, 23], and the scattered
vectorial electric field $\vec{E}^{s}$ and illumination $\vec{E}^{0}$ is
related as

$\vec{E}^{s}(\mathbf{r})=\iiint\bar{\bar{G}}(\mathbf{r}-\mathbf{r}^{\prime})\bar{\bar{V}}(\mathbf{r}^{\prime})\vec{E}^{0}(\mathbf{r}^{\prime})d\mathbf{r}^{\prime},$
(2)

where
$\bar{\bar{V}}(\mathbf{r}^{\prime})=\bar{\bar{\epsilon}}-\bar{\bar{\epsilon}}_{0}$
is the sample scattering potential tensor with $\bar{\bar{\epsilon}}_{0}$ as
the permitivity tensor of the background medium. $\bar{\bar{G}}(\mathbf{r})$
is the dyadic Green’s tensor [45]. In addition, we use illuminations with
relatively small incident angles, which have weak polarization along the
z-axis. Therefore, we mainly consider the transverse polarization of the
electrical field, $\vec{E}_{\perp}=(E_{x},E_{y})^{T}$, which will be discussed
in more detail in the next subsection. The intensities of the vectorial
electric field $\vec{E}_{\perp}$ at different orientations are then analyzed
and recorded by a polarization-sensitive optical imaging system, which we
model as a $2\times 2$ pupil matrix at each spatial frequency $\mathbf{u}$
[26]. Therefore, the measured intensity under plane wave illumination with
angle $\mathbf{u}^{\prime}$ analyzed by a polarizer with Jones vector
$\vec{a}$ is

$I^{l}(\mathbf{r},\mathbf{u}^{\prime})=\Big{\lvert}{\vec{a}_{l}}^{T}\mathcal{F}_{2d}^{-1}\Big{\\{}\bar{\bar{P}}(\mathbf{u})\mathcal{F}_{2d}\big{\\{}\vec{E}_{\perp}^{s}(\mathbf{r},\mathbf{u}^{\prime})+\vec{E}_{\perp}^{0}(\mathbf{r},\mathbf{u}^{\prime})\big{\\}}\Big{\\}}\Big{\rvert}^{2}.$
(3)

For a linear polarization analyzer oriented at $\alpha$,
$\vec{a}=[\cos{\alpha},\sin{\alpha}]^{T}$ [46, 26]. The intuition behind this
is that light can only oscillate in the same direction as the linear polarizer
[26].

####  2.2.3 Forward model and inverse problem

For our initial demonstration of _T 2oFu_, we modify the model discussed in
Section 2.2.2 by making certain approximations. The objective of this is to
craft a meaningful forward model that accurately describes our experimental
measurements, whose inverse problem is less ill-posed. First, following Saba
et al. and other closely relevant works in literature [23, 27, 28, 26], we
begin with a paraxial approximation, which assumes weak polarization along the
optical axis of the illumination and negligible interaction between traverse
and axial polarization from the sample. With this, we can simplify the
$3\times 3$ permittivity matrix to a $2\times 2$ matrix, denoted as

$\bar{\bar{\epsilon}}=\begin{bmatrix}\epsilon_{xx}(\mathbf{r})&\epsilon_{xy}(\mathbf{r})\\\
\epsilon_{yx}(\mathbf{r})&\epsilon_{yy}(\mathbf{r})\end{bmatrix}.$ (4)

While this approximation may not always be accurate for every anisotropic
sample when illuminated at high angles, it is accurate under up to a
$25\degree$ oblique illumination, based on a study using finite element
analysis [23]. Additionally, we assume that the background media is isotropic
and uniform (i.e., not spatially dependent) with a diagonal permittivity
tensor $\bar{\bar{\epsilon}}_{0}=\epsilon_{0}\mathbb{I}$. This simplifies
Green’s tensor to a diagonal matrix with the same component for each
polarization [23],

$\bar{\bar{G}}(\mathbf{r},\mathbf{r}^{\prime})=\begin{bmatrix}G(\mathbf{r},\mathbf{r}^{\prime})&0\\\
0&G(\mathbf{r},\mathbf{r}^{\prime})\end{bmatrix}.$ (5)

$G(\mathbf{r},\mathbf{r}^{\prime})=G(r)$, where
$r=\lvert\mathbf{r}-\mathbf{r}^{\prime}\rvert$, is the scalar Green’s function
that has a Weyl expansion [44, 47]

$\frac{e^{jk_{0}r}}{r}=\frac{1}{j2\pi}\int
d\mathbf{u}\frac{e^{-j(\mathbf{u}\cdot\mathbf{x}+\eta\lvert z\rvert)}}{\eta},$
(6)

with wavenumber vectors in lateral ($\mathbf{u}=(k_{x},k_{y})$) and axial
($\eta=\sqrt{k_{0}^{2}-\lvert\mathbf{u}\rvert^{2}}$) directions. $k_{0}$ is
the wavenumber of the isotropic background medium. This expansion is generally
easier to work with when we prefer to represent the object in both lateral
frequency and axial space domain.

In addition, we make the assumption that the sample being imaged is
homogeneous, which is usually assumed for many types of crystals and
biological samples [48, 10]. The permittivity matrix is then symmetric and can
be decomposed into [48]:

$\bar{\bar{\epsilon}}=\begin{bmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{bmatrix}\begin{bmatrix}\epsilon_{e}&0\\\
0&\epsilon_{o}\end{bmatrix}\begin{bmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{bmatrix}.$ (7)

The variables $\epsilon_{o}$ and $\epsilon_{e}$ represent the permittivity
values along the ordinary and extraordinary axes, respectively. The parameter
$\theta$ is the angle between the principle axis and the extraordinary axis
(also known as the slow axis). Following the convention introduced by an early
literature [49], we rename elements in Eq.4 as

$\displaystyle\epsilon_{1}$ $\displaystyle=\epsilon_{xx}$ (8a)
$\displaystyle\epsilon_{2}$ $\displaystyle=\epsilon_{yy}$ (8b)
$\displaystyle\epsilon_{3}$ $\displaystyle=\epsilon_{xy}=\epsilon_{yx}.$ (8c)

Since we use left circularly polarized illumination, the transverse scattering
potential becomes

$\begin{bmatrix}V_{1}&V_{3}\\\ V_{3}&V_{2}\end{bmatrix}\begin{bmatrix}1\\\
j\end{bmatrix}=\begin{bmatrix}V_{1}+jV_{3}\\\ V_{3}+jV_{2}\end{bmatrix},$ (9)

where

$\displaystyle V_{1}$
$\displaystyle=4\pi{k}_{0}^{2}(\epsilon_{1}-\epsilon_{0})$ (10a)
$\displaystyle V_{2}$
$\displaystyle=4\pi{k}_{0}^{2}(\epsilon_{2}-\epsilon_{0})$ (10b)
$\displaystyle V_{3}$ $\displaystyle=4\pi{k}_{0}^{2}\epsilon_{3}.$ (10c)

Finally, we refer to Ling et. al. [50] for two further approximations of the
illumination and scattering processes. Our first assumption is that the
illumination from each LED at the sample plane is a plane wave, which is
commonly used in Fourier ptychographic computational microscopy [51]. The
second approximation involves utilizing a weak object assumption that ignores
the second-order scattering term [37, 52, 53]. Further, we disregard pupil
aberration and model it as a low-pass filter $P(\mathbf{u})$ with a cutoff
frequency based on the numerical aperture of the objective lens [37]. Jones
matrix has been demonstrated feasible in previous literature for the
correction of anisotropic aberration by jointly reconstructing the pupil [26,
31], which is planned for future research. Based on these approximations, the
forward model can be expressed as [37]

$\displaystyle\begin{split}\tilde{I}^{l,m}(\mathbf{u},z=0)\approx\tilde{I}_{0}^{l,m}(\mathbf{u},z=0)+\int\Bigr{[}{H}_{Re}^{m}&(\mathbf{u},z)\cdot\tilde{V}_{Re}^{l}(\mathbf{u},z)\\\
&+{H}_{Im}^{m}(\mathbf{u},z)\cdot\tilde{V}_{Im}^{l}(\mathbf{u},z)\Bigl{]}\,dz,\end{split}$
(11)

where $\tilde{I}^{l,m}(\mathbf{u},z)$ and $\tilde{I}_{0}^{l,m}(\mathbf{u},z)$
are 2D Fourier transform of the measurement and DC term from $m^{th}$ LED
illumination analyzed by $l^{th}$ polarizer, respectively.

${\begin{split}H_{Re}^{m}(\mathbf{u},z)=\frac{jk^{2}}{2}S(\mathbf{u}_{m})\biggl{\\{}P^{*}(-\mathbf{u}_{m})\frac{e^{-j[\eta_{i}+\eta(\mathbf{u}-\mathbf{u}_{m})]z}}{\eta(\mathbf{u}-\mathbf{u}_{m})}&P(\mathbf{u}-\mathbf{u}_{m})\\\
&-P(-\mathbf{u}_{m})\frac{e^{j[\eta_{i}+\eta(\mathbf{u}+\mathbf{u}_{m})]z}}{\eta(\mathbf{u}+\mathbf{u}_{m})}P(-\mathbf{u}-\mathbf{u}_{m})\biggr{\\}}\end{split}}$
(12)

and

${\begin{split}H_{Im}^{m}(\mathbf{u},z)=\frac{-k^{2}}{2}S(\mathbf{u}_{m})\biggl{\\{}P^{*}(-\mathbf{u}_{m})\frac{e^{-j[\eta_{i}+\eta(\mathbf{u}-\mathbf{u}_{m})]z}}{\eta(\mathbf{u}-\mathbf{u}_{m})}&P(\mathbf{u}-\mathbf{u}_{m})\\\
&+P(-\mathbf{u}_{m})\frac{e^{j[\eta_{i}+\eta(\mathbf{u}+\mathbf{u}_{m})]z}}{\eta(\mathbf{u}+\mathbf{u}_{m})}P(-\mathbf{u}-\mathbf{u}_{m})\biggr{\\}}\end{split}}$
(13)

are the diffractive transfer functions in frequency and space for the real and
imaginary part of the scattering potential under $m^{th}$ LED illumination
with shape ${S}^{m}(\mathbf{u}^{\prime})$ [52, 54]. Since each LED has a very
small die area ($<170$ µm in diameter), we represent it as a delta function.
$l\in\\{0\degree,45\degree,90\degree,135\degree\\}$.
$\tilde{V}_{Re}^{l}(\mathbf{u},z)$ and $\tilde{V}_{Im}^{l}(\mathbf{u},z)$ are
2D Fourier transforms of ${V}_{Re}^{l}(\mathbf{r})$ and
${V}_{Im}^{l}(\mathbf{r})$ along lateral directions at depth $z$,
respectively, and are related to Eq. 10 via

$\displaystyle{V}_{Re}^{0\degree}(\mathbf{r})+j{V}_{Im}^{0\degree}(\mathbf{r})=V_{1}(\mathbf{r})+jV_{3}(\mathbf{r})$
(14a)
$\displaystyle\begin{split}&{V}_{Re}^{45\degree}(\mathbf{r})+j{V}_{Im}^{45\degree}(\mathbf{r})=\frac{\sqrt{2}}{2}(V_{1}(\mathbf{r})+V_{3}(\mathbf{r}))\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+j\frac{\sqrt{2}}{2}(V_{2}(\mathbf{r})+V_{3}(\mathbf{r}))\end{split}$
(14b)
$\displaystyle{V}_{Re}^{90\degree}(\mathbf{r})+j{V}_{Im}^{90\degree}(\mathbf{r})=V_{3}(\mathbf{r})+jV_{2}(\mathbf{r})$
(14c)
$\displaystyle\begin{split}&{V}_{Re}^{135\degree}(\mathbf{r})+j{V}_{Im}^{135\degree}(\mathbf{r})=\frac{\sqrt{2}}{2}(-V_{1}(\mathbf{r})+V_{3}(\mathbf{r}))\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+j\frac{\sqrt{2}}{2}(V_{2}(\mathbf{r})-V_{3}(\mathbf{r})),\end{split}$
(14d)

which are the scattering potential components corresponding to each analyzer
angle. We want to point out that $V_{1,2,3}(\mathbf{r})$ are complex variables
that may have imaginary parts; hence, the above equation does not imply
$\tilde{V}_{Re}(\mathbf{r})=V_{1}(\mathbf{r})$, for example.

For concise expression, we define a new variable
$\mathbf{v}\in\mathbb{C}^{N\times M\times T\times 3}$ representing all the
potentials $V_{1}(\mathbf{r}),V_{2}(\mathbf{r})$, and $V_{3}(\mathbf{r})$,
where $N,M,T$ are width, height, and depth of the 3D sample, respectively.
Further, we introduce the operator $\mathcal{A}^{l,m}(\cdot)$ as the forward
model for the sample illuminated by the $m^{th}$ LED and analyzed by the
$l^{th}$ linear polarizer. To reconstruct the permittivity matrix, we
formulate the inverse problem as

$\mathbf{v}=\operatorname*{arg\,min}_{\mathbf{v}}\mathcal{L}(\mathbf{v}),$
(15)

with the loss function

$\mathcal{L}(\mathbf{v})=\sum_{l}\sum_{m}\lVert\mathcal{A}^{l,m}(\mathbf{v})-{\tilde{I}}^{l,m}(\mathbf{u},z)\rVert_{2}^{2}+\gamma\text{tv}(\mathbf{v}).$
(16)

tv$(\cdot)$ is the isotropic total variation operator. $\gamma$ is a
regularization parameter empirically set to be $1\times 10^{-6}$ for all
experiments. The forward model is implemented in Pytorch and the loss function
is optimized using a stochastic gradient descent method with Nesterov momentum
acceleration [55, 56].

Subsequently, we extract polarization properties of interest, such as
orientation and birefringence from _T 2oFu_ reconstructions
$\epsilon_{1,2,3}$. The permittivity of ordinary and extraordinary axes can be
computed as

$\displaystyle\epsilon_{0}$
$\displaystyle=\bar{\epsilon}-\nicefrac{{1}}{{2}}\Delta\epsilon$ (17a)
$\displaystyle\epsilon_{e}$
$\displaystyle=\bar{\epsilon}+\nicefrac{{1}}{{2}}\Delta\epsilon,$ (17b)

and the refractive index along the ordinary and extraordinary axes
$n_{o,e}=\sqrt{\epsilon_{o,e}}$ can be further derived, along with the
averaged refractive index ($\bar{n}=n_{o}+n_{e}$) and birefringence
($\Delta{n}=n_{e}-n_{o}$). Moreover, the orientation to the slow axis is
computed as

$\theta=\begin{cases}\frac{1}{2}\arctan\nicefrac{{2\epsilon_{3}}}{{\epsilon_{1}-\epsilon_{2}}},\;&\text{if}\;\epsilon_{1}-\epsilon_{2}>0\\\
\frac{1}{2}\arctan\nicefrac{{2\epsilon_{3}}}{{\epsilon_{1}-\epsilon_{2}}}+\frac{\pi}{2},\;&\text{otherwise.}\end{cases}$
(18)

##  3 Experimental Results

Figure 2: Reconstruction results for polystyrene microspheres. (A-B) Intensity
images of the sample illuminated with center LED. (A) The data is captured
when focused at z=0 µm. (B) Images were captured when focused at -5 µm and 7
µm by mechanically moving the sample. These images serve as a reference to be
compared with the reconstruction. (C)-(E) plots the reconstructed refractive
index and birefringence. All columns in (C) and (D) share the same scale bar,
all the reflective index reconstructions (C), as well as all birefringence
reconstructions (D) and all cross-sections (E) share the same colorbar. (E)
displays cross-sections of the reconstruction in places color labeled in (A).
(F) plots the profile of reconstruction averaged over 10 microspheres. Figure
3: Reconstruction results of potato starch grain. (A) Polarized intensity
measurements from the center LED. (B) Images taken when focused at -9 µm and 7
µm serve as a reference for validating the reconstruction. (C) Tomographic
reconstruction of the sample. The orientation is coded in color, while the
birefringence is represented as brightness. All the images in (C) share the
same colorbar.

To validate the performance of the proposed method, we first show results from
a variety of calibration targets. Due to the lack of commercially available
high-resolution 3D polarization phantom, we follow previous works to validate
different aspects separately [32, 28, 57, 26]. We then show reconstructions of
a single fixed muscle fiber. Finally, we show reconstructions of a tissue
section sample from a human heart biopsy that is predictive for lethal cardiac
amyloidosis. All the images presented here are captured and then reconstructed
using a $20\times$, 0.4 NA system, except for the cardiac tissue, which is
recorded and processed with a $10\times$, 0.25 NA system.

###  3.1 Tomographic reconstruction

In this subsection, we demonstrate tomographic reconstructions of averaged
refractive index $\bar{n}$ and birefringence $\Delta n$ from isotropic and
anisotropy calibration samples. These results are acquired with a $20\times$
objective ($0.4$NA, Olympus Corporation, Japan.) imaging system. Figure 2 (A)
shows the image captured with center LED illumination when the imaging system
is focused at the middle of a polystyrene microsphere sample. The sample is
made of two layers of $3$µm-diameter microspheres immersed in $n=1.575$ oil.
(B) shows images taken with center LED illumination when focused at different
depths. The arrows in (B) highlight the microspheres that are in focus,
suggesting two layers of microspheres are placed at $-7$µm and $5$µm planes.
Note that since the refractive index of polystyrene ($n=1.60$ @ $520$nm) is
very close to the background medium, the contrast of in-focus microspheres is
very low. (C)-(D) shows the reconstructed refractive index and birefringence
at two different depths. The accuracy of the tomographic depth reconstruction
is validated with images displayed in (B). Since polystyrene is an isotropic
material, the reconstruction shows no anisotropy properties except on the
edge, which agrees with the well-recognized edge birefringence phenomena[58]
reported in prior literature[32, 26]. (E) shows the cross-sections of the
reconstructed reflective index and birefringence. (F) plots the profile of the
reconstructed birefringence and reflective index averaged across 10
microspheres.

In addition, we show reconstructions of an anisotropic potato starch sample in
Fig. 3. The sample consists of two potato starch grains immersed in $n=1.515$
oil at different depths. Fig. 3(A) depicts captured intensity images
illuminated with the center LED at four different polarizations. We can see
the spiral patterns imaged with different polarization differs from each other
noticeably. Similarly, Fig. 3(B) shows images when the system is focused at
different depths of the sample. The black arrows point at grains that are in
focus, suggesting the two potato starch grains are suspended at $-9\,$µm and
$6\,$µm. Fig. 3(C) shows the reconstructed orientation and birefringence. To
best visualize the results, we follow the convention [21, 28, 26] to display
this multidimensional data using an HSV colormap, where saturation is set to
one, value is associated with birefringence, and orientation is coded in hue.
The reconstructed structures agree with starch grain reconstructions reported
in previous holography-based literatures [23, 25].

Figure 4: Reconstructions of monosodium urate (MSU) crystals. (A) shows the
reconstructed polarization orientation of the crystals at two different
depths, labeled with their respective structural directions. (B) renders the
birefringence reconstruction of the sample in 3D.

###  3.2 Orientation measurement

To verify the accuracy of the reconstructed orientation, we test our method on
a sample made with monosodium urate (MSU). MSU are needle-shaped crystals
precipitated from uric acid that could trigger robust inflammation such as
acute arthritis and other immune activations that cause severe pain in
patients [59]. Figure 4(A) displays the reconstructed birefringence and
orientation at two slightly different depths. These results were also obtained
with a $20\times$ objective imaging system and suggest that the reconstructed
orientation values of line-shaped MSU crystals follow the structural direction
(labeled next to each MSU crystal). This is in agreement with reconstruction
results from LC-PolScope-based methods found in the literature [28, 57].
Additionally, Figure 4(B) presents a 3D rendering of the reconstructed
birefringence.

###  3.3 Muscle fiber assessment

Figure 5: Reconstructions of a muscle fiber. (A) shows the image of a muscle
fiber with the center LED illumination. The imaging system is focused in the
middle of the muscle fiber. (B) plots the reconstructed birefringence. The
zoom-in region resembles reconstruction results reported in [34] using a
second harmonic generation imaging approach. (c) shows the image of the same
muscle fiber at a different region, where a non-muscle fiber with a ninety-
degree bend is placed below the muscle fiber (see -10µm in D). The imaging
system is focused between this and the muscle fiber. (D) shows the
reconstructed orientation at different depths, with a zoom-in showing the fine
sarcomere structure of muscle tissue.

High contrast and high-resolution structural imaging of intrinsic signals in
muscle skeletal fibers is an important task for the rapid detection of changes
in myofibrillar organization that can lead to skeletal myopathies. Currently,
3D muscle tissue is typically imaged by complex and expensive systems, such as
second harmonic generation (SHG) microscopy. SHG exploits the contrast of
polarization properties from the non-linear susceptibility in molecules like
myosin, by using high-power, ultra-short pulsed lasers in a point scanning
configuration [34]. Here we show 3D  _T 2oFu_ reconstruction results of an
isolated healthy muscle fiber using our inexpensive, LED-based and scanning-
free system. Figure 5(A) shows an image of a muscle fiber captured with center
LED illumination. Figure 5(B) shows a volumetric reconstruction of the muscle
fiber. The cross sections of regions highlighted in Fig. 5(A) are depicted.
The zoom-in region of the image in x-y resembles reconstructions reported in
[34]. Figure 5(C) shows an image of the same muscle fiber from a different
field of view, where a non-muscle fiber with a ninety-degree bend is placed
below the muscle fiber. Figure 5(D) depicts reconstructed orientation and
birefringence at two different depths, showing the regular pattern of a
healthy muscle fiber. The zoom-in region highlights the muscle grains with
consistent orientations, in agreement with the results reported by Both et al.
[34]. The reconstruction shows a change in orientation at the bend of the non-
muscle fiber (D1), while the orientation of the muscle fiber remains constant.
Furthermore, the spatial resolution of the reconstruction is sufficient to
clearly resolve the muscle filaments (D3).

###  3.4 Imaging cardiac amyloidosis

Figure 6: Images of a heart tissue sample with cardiac amyloidosis. (A)
displays a brightfield image. (B) shows a cross-polarized image taken with a
color microscope. (C)-(D) present reconstructed refractive index and
birefringence, along with zoom-ins of the boxed region at different depths
depicted in (E)-(F). (C) and (E) share the same color bar, while (D) and (F)
share another common color bar.

Finally, we apply our method to image a heart tissue sample that has cardiac
amyloidosis. Cardiac amyloidosis is a lethal disease that affects more than
12,000 patients in the US alone with $<5\%$ 10-year survival rate [60]. In
current practice, the biopsied tissue is first frozen and thinly sliced, then
stained with a congo red-colored dye and inspected under a cross-polarized
microscope. Figure 6(A) shows the brightfield image. (B) shows the same
region, imaged with a cross-polarized color microscope. The vibrant apple
green color suggests mostly likely amyloid protein has built up inside the
tissue sample. (C)-(D) depict the reconstructed refractive index and
birefringence, while zoom-ins of the boxed region at different depths are
shown in (E)-(F). Since the sample is thinly sliced, we do not observe
noticeable structural changes in different layers. However, the structure of
the birefringence reconstruction is correlated with the color-stained cross-
polarized image, which could potentially be useful for rapid on-site
inspections in the future.

##  4 Discussion and Conclusion

In this article, we introduce _T 2oFu_, a new non-scanning microscopy method
that reconstructs volumetric permittivity metrics of samples based on
computational illumination strategies to record polarized measurements and
retrieve phase information. Using relatively low-NA objectives, we demonstrate
that _T 2oFu_ can provide polarization-sensitive tomographic reconstruction
for various calibration samples and biological specimens that are potentially
useful for future scientific and clinical studies.

To ensure successful clinical translation, there are a few improvements that
can be made in the next step. First, to increase the frame rate of _T 2oFu_,
high flash-rate LEDs need to be deployed [61]. Novel sensitive camera sensors
[62, 63, 64] will also play important roles in maintaining a good signal-to-
noise ratio. In terms of the reconstruction algorithm, as a first
demonstration, we use a gradient-based method implemented with auto-
differentiation to optimize the loss function in Eq. 16 [55]. As the data
fidelity term has a closed-form solution, variable-splitting methods can be
very effective and allow incorporating advanced regularization that does not
have implicit forms [65, 66]. Moreover, due to spectral scanning, so far, _T
2oFu_ has a lower frame rate compared to most holographic methods in the
literature. To improve _T 2oFu_’s throughput, in addition to using higher-
speed LEDs and a more sensitive camera [37], adopting a multi-aperture
approach appears to be another promising direction [67, 68, 69, 70, 71].
Finally, concerning modeling, to relieve the ill-posedness of the problem, we
have approximated the permittivity tensor with its lateral components in this
work. Albeit widely used [23, 27, 28, 26, 72, 73, 74], this simplification
disregards out-of-plane anisotropy. Recent works have shown under uni-axial
approximation, tomography of 3D polarization orientation can be retrieved
[32]. Further investigation in this direction to extend the current method to
extract valuable out-of-plane information is planned [32, 75].

## Ethical approval

Animal tissue sample studies were approved by the Friedrich-Alexander-
University Erlangen-Nürnberg. All mice were maintained in the animal facility
in a 12-h light–dark cycle with access to food and water. Human tissue sample
studies were conducted with approval from the Duke University Health System
Institutional Review Board (DUHS IRB). The DUHS IRB determined that the
following protocol meets the criteria for a declaration of exemption from
further IRB review as described in 45 CFR 46.101(b), 45 CFR 46.102 (f), or 45
CFR 46.102 (d), satisfies the Privacy Rule as described in 45 CFR 164.512(i),
and satisfies Food and Drug Administration regulations as described in 21 CFR
56.104, where applicable.

## Funding

The authors would also like to thank to a Duke-Coulter Translational
Partnership and funding from a 3M Nontenured Faculty Award.

## Conflict of interest

S.X., X.D., L.K., J.N, C.G, and R.H. have submitted a patent application
related to this work, assigned to Duke University.

## References

  * [1] Y. Park, C. Depeursinge, and G. Popescu, “Quantitative phase imaging in biomedicine,” _Nature photonics_ , vol. 12, no. 10, pp. 578–589, 2018.
  * [2] T. Ling, K. C. Boyle, V. Zuckerman, T. Flores, C. Ramakrishnan, K. Deisseroth, and D. Palanker, “High-speed interferometric imaging reveals dynamics of neuronal deformation during the action potential,” _Proceedings of the National Academy of Sciences_ , vol. 117, no. 19, pp. 10 278–10 285, 2020.
  * [3] Z. El-Schich, A. Leida Mölder, and A. Gjörloff Wingren, “Quantitative phase imaging for label-free analysis of cancer cells—focus on digital holographic microscopy,” _Applied Sciences_ , vol. 8, no. 7, p. 1027, 2018\.
  * [4] K. Zhanghao, W. Liu, M. Li, Z. Wu, X. Wang, X. Chen, C. Shan, H. Wang, X. Chen, Q. Dai, P. Xi, and D. Jin, “High-dimensional super-resolution imaging reveals heterogeneity and dynamics of subcellular lipid membranes,” _Nature communications_ , vol. 11, no. 1, pp. 1–10, 2020.
  * [5] J. Lu, H. Mazidi, T. Ding, O. Zhang, and M. D. Lew, “Single-molecule 3d orientation imaging reveals nanoscale compositional heterogeneity in lipid membranes,” _Angewandte Chemie International Edition_ , vol. 59, no. 40, pp. 17 572–17 579, 2020.
  * [6] W. J. Schmidt, _Die Bausteine des Tierkörpers in polarisiertem Lichte_. F. Cohen, 1924.
  * [7] S. Inoue, “Polarization optical studies of the mitotic spindle,” _Chromosoma_ , vol. 5, no. 1, pp. 487–500, 1953.
  * [8] G. Nomarski, “Differential microinterferometer with polarized waves,” _J. Phys. Radium Paris_ , vol. 16, p. 9S, 1955.
  * [9] R. Oldenbourg, “Polarized light microscopy: principles and practice,” _Cold Spring Harbor Protocols_ , vol. 2013, no. 11, pp. pdb–top078 600, 2013\.
  * [10] S. B. Mehta, M. Shribak, and R. Oldenbourg, “Polarized light imaging of birefringence and diattenuation at high resolution and high sensitivity,” _Journal of Optics_ , vol. 15, no. 9, p. 094007, 2013.
  * [11] E. M. Spiesz, W. Kaminsky, and P. K. Zysset, “A quantitative collagen fibers orientation assessment using birefringence measurements: calibration and application to human osteons,” _Journal of structural biology_ , vol. 176, no. 3, pp. 302–306, 2011.
  * [12] M. Tadayon, S. Amini, A. Masic, and A. Miserez, “The mantis shrimp saddle: a biological spring combining stiffness and flexibility,” _Advanced Functional Materials_ , vol. 25, no. 41, pp. 6437–6447, 2015.
  * [13] A. Le Gratiet, M. d’Amora, M. Duocastella, R. Marongiu, A. Bendandi, S. Giordani, P. Bianchini, and A. Diaspro, “Zebrafish structural development in mueller-matrix scanning microscopy,” _Scientific reports_ , vol. 9, no. 1, pp. 1–10, 2019.
  * [14] L. C. U. Junqueira, G. Bignolas, and R. R. Brentani, “Picrosirius staining plus polarization microscopy, a specific method for collagen detection in tissue sections,” _The Histochemical journal_ , vol. 11, pp. 447–455, 1979\.
  * [15] H. V. Desai, W. S. Aronow, S. J. Peterson, and W. H. Frishman, “Cardiac amyloidosis: approaches to diagnosis and management,” _Cardiology in review_ , vol. 18, no. 1, pp. 1–11, 2010.
  * [16] C. W. Pirnstill and G. L. Coté, “Malaria diagnosis using a mobile phone polarized microscope,” _Scientific reports_ , vol. 5, no. 1, pp. 1–13, 2015\.
  * [17] T. Liu, M. Lu, B. Chen, Q. Zhong, J. Li, H. He, H. Mao, and H. Ma, “Distinguishing structural features between crohn’s disease and gastrointestinal luminal tuberculosis using mueller matrix derived parameters,” _Journal of Biophotonics_ , vol. 12, no. 12, p. e201900151, 2019\.
  * [18] Y. Jiao, M. E. Kandel, X. Liu, W. Lu, and G. Popescu, “Real-time jones phase microscopy for studying transparent and birefringent specimens,” _Optics Express_ , vol. 28, no. 23, pp. 34 190–34 200, 2020.
  * [19] S. Shin, K. Lee, Z. Yaqoob, P. T. So, and Y. Park, “Reference-free polarization-sensitive quantitative phase imaging using single-point optical phase conjugation,” _Optics express_ , vol. 26, no. 21, pp. 26 858–26 865, 2018.
  * [20] B. Ge, Q. Zhang, R. Zhang, J.-T. Lin, P.-H. Tseng, C.-W. Chang, C.-Y. Dong, R. Zhou, Z. Yaqoob, I. Bischofberger, and P. So, “Single-shot quantitative polarization imaging of complex birefringent structure dynamics,” _ACS Photonics_ , vol. 8, no. 12, pp. 3440–3447, 2021.
  * [21] T. Liu, K. de Haan, B. Bai, Y. Rivenson, Y. Luo, H. Wang, D. Karalli, H. Fu, Y. Zhang, J. FitzGerald, and A. Ozcan, “Deep learning-based holographic polarization microscopy,” _ACS photonics_ , vol. 7, no. 11, pp. 3023–3034, 2020.
  * [22] J. Van Rooij and J. Kalkman, “Polarization contrast optical diffraction tomography,” _Biomedical Optics Express_ , vol. 11, no. 4, pp. 2109–2121, 2020.
  * [23] A. Saba, J. Lim, A. B. Ayoub, E. E. Antoine, and D. Psaltis, “Polarization-sensitive optical diffraction tomography,” _Optica_ , vol. 8, no. 3, pp. 402–408, 2021.
  * [24] S. Shin, J. Eun, S. S. Lee, C. Lee, H. Hugonnet, D. K. Yoon, S.-H. Kim, J. Jeong, and Y. Park, “Tomographic measurement of dielectric tensors at optical frequency,” _Nature Materials_ , vol. 21, no. 3, pp. 317–324, 2022\.
  * [25] A. M. Taddese, M. Lo, N. Verrier, M. Debailleul, and O. Haeberlé, “Jones tomographic diffractive microscopy with a polarized array sensor,” _Optics Express_ , vol. 31, no. 5, pp. 9034–9051, 2023.
  * [26] X. Dai, S. Xu, X. Yang, K. C. Zhou, C. Glass, P. C. Konda, and R. Horstmeyer, “Quantitative jones matrix imaging using vectorial fourier ptychography,” _Biomedical Optics Express_ , vol. 13, no. 3, pp. 1457–1470, 2022.
  * [27] Q. Song, A. Baroni, R. Sawant, P. Ni, V. Brandli, S. Chenot, S. Vézian, B. Damilano, P. de Mierry, S. Khadir, P. Ferrand, and P. Genevet, “Ptychography retrieval of fully polarized holograms from geometric-phase metasurfaces,” _Nature communications_ , vol. 11, no. 1, pp. 1–8, 2020.
  * [28] S. Song, J. Kim, S. Hur, J. Song, and C. Joo, “Large-area, high-resolution birefringence imaging with polarization-sensitive fourier ptychographic microscopy,” _ACS Photonics_ , vol. 8, no. 1, pp. 158–165, 2021.
  * [29] S. Hur, S. Song, S. Kim, and C. Joo, “Polarization-sensitive differential phase-contrast microscopy,” _Optics Letters_ , vol. 46, no. 2, pp. 392–395, 2021.
  * [30] S.-M. Guo, L.-H. Yeh, J. Folkesson, I. E. Ivanov, A. P. Krishnan, M. G. Keefe, E. Hashemi, D. Shin, B. B. Chhun, N. H. Cho, M. D. Leonetti, M. H. Han, T. J. Nowakowski, and S. B. Mehta, “Revealing architectural order with quantitative label-free imaging and deep learning,” _elife_ , vol. 9, p. e55502, 2020.
  * [31] A. Baroni, M. Allain, P. Li, V. Chamard, and P. Ferrand, “Joint estimation of object and probes in vectorial ptychography,” _Optics Express_ , vol. 27, no. 6, pp. 8143–8152, 2019.
  * [32] L.-H. Yeh, I. E. Ivanov, S.-M. Guo, B. B. Chhun, E. Hashemi, M. H. Han, and S. B. Mehta, “upti: uniaxial permittivity tensor imaging of intrinsic density and anisotropy,” in _Novel Techniques in Microscopy_. Optical Society of America, 2021, pp. NM3C–4.
  * [33] S. Xu, X. Dai, X. Yang, K. C. Zhou, K. Kim, V. Pathak, C. Glass, and R. Horstmeyer, “Tensorial tomographic differential phase-contrast microscopy,” in _2022 IEEE International Conference on Computational Photography (ICCP)_. IEEE, 2022, pp. 1–11.
  * [34] M. Both, M. Vogel, O. Friedrich, F. v. Wegner, T. Ku¨ nsting, R. H. Fink, and D. Uttenweiler, “Second harmonic imaging of intrinsic signals in muscle fibers in situ,” _Journal of biomedical optics_ , vol. 9, no. 5, pp. 882–892, 2004.
  * [35] K. Seo, S. Cho, J.-H. Lee, J. H. Kim, B. Lee, H. Jang, Y. Kim, H. M. Cho, S. Lee, Y. Park _et al._ , “Symmetry breaking of hpscs in micropattern generates a polarized spinal cord-like organoid (psco) with dorsoventral organization,” _BioRxiv_ , pp. 2021–09, 2021.
  * [36] C. Joo, S. Song, J. Kim, T. Moon, B. Seong, W. Kim, C.-H. Yoo, and J.-K. Choi, “Polarization sensitive intensity diffraction tomography,” _researchsquare preprint_ , 2022.
  * [37] J. Li, A. Matlock, Y. Li, Q. Chen, C. Zuo, and L. Tian, “High-speed in vitro intensity diffraction tomography,” _Advanced Photonics_ , vol. 1, no. 6, pp. 066 004–066 004, 2019.
  * [38] R. Horstmeyer, J. Chung, X. Ou, G. Zheng, and C. Yang, “Diffraction tomography with fourier ptychography,” _Optica_ , vol. 3, no. 8, pp. 827–835, 2016\.
  * [39] T.-A. Pham, E. Soubies, A. Goy, J. Lim, F. Soulez, D. Psaltis, and M. Unser, “Versatile reconstruction framework for diffraction tomography with intensity measurements and multiple scattering,” _Optics express_ , vol. 26, no. 3, pp. 2749–2763, 2018.
  * [40] M. A. Lodhi, Y. Ma, H. Mansour, P. T. Boufounos, and D. Liu, “Inverse multiple scattering with phaseless measurements,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 1519–1523.
  * [41] J. Li, N. Zhou, J. Sun, S. Zhou, Z. Bai, L. Lu, Q. Chen, and C. Zuo, “Transport of intensity diffraction tomography with non-interferometric synthetic aperture for three-dimensional label-free microscopy,” _Light: Science & Applications_, vol. 11, no. 1, p. 154, 2022.
  * [42] S. Chowdhury, M. Chen, R. Eckert, D. Ren, F. Wu, N. Repina, and L. Waller, “High-resolution 3d refractive index microscopy of multiple-scattering samples from intensity images,” _Optica_ , vol. 6, no. 9, pp. 1211–1219, 2019.
  * [43] L. Tian and L. Waller, “3d intensity and phase imaging from light field measurements in an led array microscope,” _optica_ , vol. 2, no. 2, pp. 104–111, 2015.
  * [44] M. Born and E. Wolf, _Principles of optics: electromagnetic theory of propagation, interference and diffraction of light_. Elsevier, 2013.
  * [45] A. D. Yaghjian, _A direct approach to the derivation of electric dyadic Green’s functions_. Department of Commerce, National Bureau of Standards, Institute for Basic …, 1978, vol. 13.
  * [46] P. Ferrand, M. Allain, and V. Chamard, “Ptychography in anisotropic media,” _Optics letters_ , vol. 40, no. 22, pp. 5144–5147, 2015.
  * [47] H. Weyl, “A new extension of the theory of relativity,” _Ann. Physik_ , vol. 59, p. 101, 1919.
  * [48] R. A. Chipman, W. S. T. Lam, and G. Young, _Polarized light and optical systems_. CRC press, 2018.
  * [49] R. C. Jones, “A new calculus for the treatment of optical systemsi. description and discussion of the calculus,” _Josa_ , vol. 31, no. 7, pp. 488–493, 1941.
  * [50] R. Ling, W. Tahir, H.-Y. Lin, H. Lee, and L. Tian, “High-throughput intensity diffraction tomography with a computational microscope,” _Biomedical optics express_ , vol. 9, no. 5, pp. 2130–2141, 2018.
  * [51] P. C. Konda, L. Loetgering, K. C. Zhou, S. Xu, A. R. Harvey, and R. Horstmeyer, “Fourier ptychography: current applications and future promises,” _Optics express_ , vol. 28, no. 7, pp. 9603–9630, 2020.
  * [52] M. Chen, L. Tian, and L. Waller, “3d differential phase contrast microscopy,” _Biomedical optics express_ , vol. 7, no. 10, pp. 3940–3950, 2016.
  * [53] A. B. Ayoub, J. Lim, E. E. Antoine, and D. Psaltis, “3d reconstruction of weakly scattering objects from 2d intensity-only measurements using the wolf transform,” _Optics Express_ , vol. 29, no. 3, pp. 3976–3984, 2021.
  * [54] N. Streibl, “Three-dimensional imaging by a microscope,” _JOSA A_ , vol. 2, no. 2, pp. 121–127, 1985.
  * [55] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [56] Y. E. Nesterov, “A method of solving a convex programming problem with convergence rate o$\backslash$bigl(k^2$\backslash$bigr),” in _Doklady Akademii Nauk_ , vol. 269, no. 3. Russian Academy of Sciences, 1983, pp. 543–547.
  * [57] B. Bai, H. Wang, T. Liu, Y. Rivenson, J. FitzGerald, and A. Ozcan, “Pathological crystal imaging with single-shot computational polarized light microscopy,” _Journal of Biophotonics_ , vol. 13, no. 1, p. e201960036, 2020\.
  * [58] R. Oldenbourg, “Analysis of edge birefringence,” _Biophysical journal_ , vol. 60, no. 3, pp. 629–641, 1991.
  * [59] Y. Shi, A. D. Mucsi, and G. Ng, “Monosodium urate crystals in inflammation and immunity,” _Immunological reviews_ , vol. 233, no. 1, pp. 203–217, 2010\.
  * [60] T. P. Quock, T. Yan, E. Chang, S. Guthrie, and M. S. Broder, “Epidemiology of al amyloidosis: a real-world study using us claims data,” _Blood advances_ , vol. 2, no. 10, pp. 1046–1053, 2018.
  * [61] M. B. Reiser and M. H. Dickinson, “A modular display system for insect behavioral neuroscience,” _Journal of neuroscience methods_ , vol. 167, no. 2, pp. 127–139, 2008.
  * [62] L. Tian, Z. Liu, L.-H. Yeh, M. Chen, J. Zhong, and L. Waller, “Computational illumination for high-speed in vitro fourier ptychographic microscopy,” _Optica_ , vol. 2, no. 10, pp. 904–911, 2015.
  * [63] T. Aidukas, P. C. Konda, A. R. Harvey, M. J. Padgett, and P.-A. Moreau, “Phase and amplitude imaging with quantum correlations through fourier ptychography,” _Scientific reports_ , vol. 9, no. 1, p. 10445, 2019.
  * [64] X. Yang, P. C. Konda, S. Xu, L. Bian, and R. Horstmeyer, “Quantized fourier ptychography with binary images from spad cameras,” _Photonics Research_ , vol. 9, no. 10, pp. 1958–1969, 2021.
  * [65] Y. Sun, S. Xu, Y. Li, L. Tian, B. Wohlberg, and U. S. Kamilov, “Regularized fourier ptychography using an online plug-and-play algorithm,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 7665–7669.
  * [66] K. C. Zhou and R. Horstmeyer, “Diffraction tomography with a deep image prior,” _Optics express_ , vol. 28, no. 9, pp. 12 872–12 896, 2020.
  * [67] D. L. Wakefield, R. Graham, K. Wong, S. Wang, C. Hale, and C.-C. Yu, “Cellular analysis using label-free parallel array microscopy with fourier ptychography,” _Biomedical Optics Express_ , vol. 13, no. 3, pp. 1312–1327, 2022.
  * [68] M. Harfouche, K. Kim, K. C. Zhou, P. C. Konda, S. Sharma, E. E. Thomson, C. Cooke, S. Xu, L. Kreiss, A. Chaware _et al._ , “Imaging across multiple spatial scales with the multi-camera array microscope,” _Optica_ , vol. 10, no. 4, pp. 471–480, 2023.
  * [69] E. Thomson, M. Harfouche, K. Kim, P. Konda, C. W. Seitz, C. Cooke, S. Xu, W. S. Jacobs, R. Blazing, Y. Chen _et al._ , “Gigapixel imaging with a novel multi-camera array microscope,” _Elife_ , vol. 11, p. e74988, 2022.
  * [70] X. Yang, M. Harfouche, K. C. Zhou, L. Kreiss, S. Xu, P. C. Konda, K. Kim, and R. Horstmeyer, “Multi-modal imaging using a cascaded microscope design,” _Optics Letters_ , vol. 48, no. 7, pp. 1658–1661, 2023.
  * [71] X. Yao, V. Pathak, H. Xi, A. Chaware, C. Cooke, K. Kim, S. Xu, Y. Li, T. Dunn, P. C. Konda _et al._ , “Increasing a microscope’s effective field of view via overlapped imaging and machine learning,” _Optics Express_ , vol. 30, no. 2, pp. 1745–1761, 2022.
  * [72] X. Dai, P. C. Konda, S. Xu, and R. Horstmeyer, “Polarization and phase imaging using an led array microscope,” in _Polarized Light and Optical Angular Momentum for Biomedical Diagnostics_ , vol. 11646. SPIE, 2021, p. 116460U.
  * [73] X. Dai, P. C. Konda, S. Xu, and R. W. Horstmeyer, “Towards a vectorial treatment of fourier ptychographic microscopy,” in _Computational Optical Sensing and Imaging_. Optical Society of America, 2020, pp. CF2C–3.
  * [74] S. Xu, X. Dai, X. Yang, K. Zhou, P. Konda, and R. Horstmeyer, “Imaging anisotropy with vectorial fourier ptychography,” in _Computational Optics 2021_ , vol. 11875. SPIE, 2021, p. 118750F.
  * [75] B. Yang, N.-J. Jan, B. Brazile, A. Voorhees, K. L. Lathrop, and I. A. Sigal, “Polarized light microscopy for 3-dimensional mapping of collagen fiber architecture in ocular tissues,” _Journal of biophotonics_ , vol. 11, no. 8, p. e201700356, 2018.


# Fully kinetic simulations of strong steady-state collisional planar plasma
shocks

S. E. Anderson<EMAIL_ADDRESS>L. Chacón  W. T. Taitano  A. N. Simakov  B.
D. Keenan  Los Alamos National Laboratory, Los Alamos, NM 87545, USA

(August 27, 2024)

###### Abstract

We report on the first steady-state simulations of strong plasma shocks with
fully kinetic ions _and_ electrons, independently confirmed by two fully
kinetic codes (an Eulerian continuum and a Lagrangian particle-in-cell). While
kinetic electrons do not fundamentally change the shock structure as compared
with fluid electrons, we find an appreciable rearrangement of the preheat
layer, associated with nonlocal electron heat transport effects. The electron
heat flux profile qualitatively agrees between kinetic and fluid electron
models, suggesting a certain level of “stiffness”, though substantial
nonlocality is observed in the kinetic heat flux. We also find good agreement
with nonlocal electron heat-flux closures proposed in the literature. Finally,
in contrast to the classical hydrodynamic picture, we find a significant
collapse in the ‘precursor’ electric-field shock at the preheat layer edge,
which correlates with the electron-temperature gradient relaxation.

PACS numbers

52.65.-y, 52.25.Dg, 52.25.-b, 52.65.Ff, 52.35.Tc

††preprint: APS/123-QED

## Introduction:

Figure 1: Left: Spatial profiles of $T_{i}$ (solid lines) and $T_{e}$ (dashed
lines) normalized to the downstream limit for hybrid iFP (fluid electrons and
kinetic ions, Hyb.), fully kinetic iFP (Kin.), and VPIC. Center: Spatial
profiles of the temperature difference $T_{e}-T_{i}$ clearly showing a kinetic
suppression of the temperature separation in the preheat layer. Right: Spatial
profiles of $Kn_{e,\nabla T}$ for hybrid (Hyb.) and fully kinetic (Kin.) iFP
simulations.

Strong shocks are present in a variety of high-energy-density (HED)
environments, including inertial confinement fusion (ICF) capsule implosions
(Rosenberg  _et al._ , 2015, 2014; Rinderknecht  _et al._ , 2014). Given the
role of shock propagation in ICF compression and yield, it is important to
understand the impact that kinetic effects may have in the shock structure,
and its imprint on the imploding capsule (e.g., Ref. (Taitano  _et al._ ,
2018)). This is relevant, as the state-of-the-art for simulating ICF capsule
implosions is radiation-hydrodynamics (rad-hydro), which is only strictly
valid in systems with small Knudsen numbers and therefore only accurate for
weak shocks with $M\sim 1$ (Simakov  _et al._ , 2017). Nevertheless, rad-hydro
has been used extensively to investigate strong plasma shocks (Jukes, 1957;
Shafranov, 1957; Jaffrin and Probstein, 1964; Grewal, 1973; Masser, Wohlbier,
and Lowrie, 2011; Zel’dovich and Raizer, 1967). The classical structure of a
strong collisional plasma shock exhibits, from upstream to downstream (see
e.g., (Zel’dovich and Raizer, 1967)): (1) a prominent preheat layer wherein
the electron temperature, $T_{e}$, exceeds that of the ions, $T_{i}$, (2) an
embedded ion compression shock wherein $T_{i}$ increases rapidly and surpasses
that of the electrons, and (3) a region of ion-electron temperature
equilibration. The preheat and relaxation regions [(1) and (3)] are
approximately of width $\approx
Z_{i}^{2}\left(m_{i}/m_{e}\right)^{1/2}\lambda_{ii}^{DS}$ (where $m_{i}$,
$m_{e}$, $Z_{i}$, and $\lambda_{ii}^{DS}$ are the ion and electron masses, the
ion charge, and the ion-ion downstream mean-free-path, respectively), while
the embedded compression shock is of order several $\lambda_{ii}^{DS}$.

To date, kinetic studies of strong plasma shocks have focused on ions, and
have thus been performed almost exclusively with ‘hybrid’ kinetic codes
employing a fluid-electron model coupled with a kinetic-ion Vlasov-Fokker-
Planck (VFP) description (Casanova, Larroche, and Matte, 1991; Larroche, 1993;
Vidal  _et al._ , 1993; Keenan  _et al._ , 2017, 2018). In these codes, the
fluid electrons are quasi-neutral and ambipolar, with $T_{e}$ determined from
the electron energy equation and with the electron heat flux modeled as the
Braginskii/Spitzer-Härm expression (i.e.,
$\mathbf{\mathscr{Q}}_{Brag.,e}=-\kappa_{e}\nabla T_{e}$, $\kappa_{e}\propto
T_{e}^{5/2}$), usually with some variety of heat flux limiter:
$\left|\mathbf{\mathcal{\mathscr{Q}}}_{e}\right|=\min(f_{lim}\mathcal{\mathscr{Q}}_{FS,e},\left|\mathbf{\mathcal{\mathscr{Q}}}_{Brag.,e}\right|)$,
where $f_{lim}$ is effectively a ‘tuning knob’ to prevent faster-than-
streaming diffusion, and $\mathcal{\mathscr{Q}}_{FS,e}=n_{e}T_{e}v_{th,e}$ is
the free-streaming thermal flux. However, it is known that local fluid models
of the electron heat flux are only valid for sufficiently small electron
Knudsen numbers, (defined as the particle mean-free-path over a characteristic
gradient length scale, $Kn=\lambda/L$) $Kn_{e}\lesssim 10^{-3}$. For
$Kn_{e}\gtrsim 10^{-3}$, we expect to see deviations from the fluid results.

We note that sophisticated models for nonlocal electron heat flux have been
explored (such as the Luciani-Mora-Virmont (LMV) model (Luciani, Mora, and
Virmont, 1983; Bendib, Luciani, and Matte, 1988; Schurtz, Nicolaï, and
Busquet, 2000; Cao, Moses, and Delettrez, 2015)). Several of these models were
investigated in Ref. (Vidal  _et al._ , 1995) with strong planar plasma
shocks, using the kinetic ion code FPION (Casanova, Larroche, and Matte, 1991;
Larroche, 1993; Vidal  _et al._ , 1993). However, while comparison was made to
a kinetic electron code (Matte and Virmont, 1982), it was not self-
consistently coupled to the kinetic ions and utilized an expansion in
spherical harmonics, neglecting anisotropy in the collision operator. Thus, to
the authors’ knowledge, a comparison of these nonlocal electron heat flux
models with fully kinetic self-consistent VFP simulations is not yet
available.

Recently, Zhang _et al._ published results of fully kinetic particle-in-cell
simulations for an $M\approx 3$ piston-driven planar plasma shock (Zhang  _et
al._ , 2021). Unfortunately, there are some significant limitations in their
work. Firstly, the shock has not detached from the piston; secondly, the total
simulation time is less than the ion-electron thermal equilibration time
(Huba, 2013); thirdly, the simulation domain is roughly the shock width. Thus,
their solution is transient and has not yet relaxed to steady state. This is
confirmed by the fact that $T_{e}$ is non-monotonic and concave-upward in the
preheat layer [see their Fig. 4(b)], in significant contrast to what is
expected.

In this Letter, we present, for the first time, self-consistent simulations of
strong _steady-state_ collisional planar plasma shocks with fully kinetic
electrons and ions. We have performed these simulations with two different
fully kinetic codes (an Eulerian continuum and a Lagrangian particle-in-cell),
finding excellent agreement. The simulations are performed with a realistic
electron-ion mass ratio. We find that, while kinetic electrons do not
fundamentally change the structure of the shock as compared with fluid
electrons, they do modify appreciably the structure of the preheat layer. In
particular, the sharp electron-temperature gradient at the front of the
preheat layer disappears, and average ion-electron temperature separation
there decreases by about 20%. The rearrangement of the electron temperature
profile in the preheat layer can be traced to nonlocal electron heat transport
effects due to the kinetic nature of electrons. Moreover, we find that the
electron heat-flux profile displays “stiffness” in that it qualitatively
agrees remarkably well between fluid and kinetic electron models, demanding
adjustments in the corresponding electron temperature profile. We also verify
the LMV nonlocal electron heat-flux closures proposed in the literature
against self-consistent full VFP simulations (to the authors’ knowledge, the
first time this has been done), and find them to provide a reasonable
agreement with the high-fidelity simulations.

## Problem setup:

To study this problem, we employ the Eulerian VFP code iFP (Anderson  _et al._
, 2020; Taitano, Chacón, and Simakov, 2018; Taitano  _et al._ , 2021a, b), and
the Lagrangian code VPIC (Bowers  _et al._ , 2008a, b, 2009). The iFP code
solves the coupled VFP equations for each of the plasma species in a 1D planar
electrostatic approximation and the electric field is given by the 1D Ampère
equation:

$\displaystyle\partial_{t}\left(f_{\alpha}\right)+\partial_{x}\left(v_{\parallel}f_{\alpha}\right)+\frac{Z_{\alpha}e}{m_{\alpha}}E_{\parallel}\partial_{v_{\parallel}}\left(f_{\alpha}\right)$
$\displaystyle=\sum_{\beta=1}^{N_{s}}C_{\alpha\beta},$ (1)
$\displaystyle\epsilon_{0}\partial_{t}\left(E_{\parallel}\right)+\sum_{\alpha=1}^{N_{s}}Z_{\alpha}enu_{\parallel,\alpha}$
$\displaystyle=0,$ (2)

where $\alpha,\beta$ denote the species index, $f_{\alpha}$ is the particle
distribution function for each species; the symbols $Z_{\alpha}$ and
$m_{\alpha}$ indicate each species’ particle charge and mass, and $e$ is the
proton charge. The electric field and $\alpha$-species bulk velocity in the
shock propagation direction are denoted by $E_{\parallel}$ and
$u_{\parallel,\alpha}$, respectively. We denote vector components along the
spatial direction of the shock propagation ($x$) with the subscript ∥, while
perpendicular components are denoted with the subscript ⊥. The VPIC code is a
3D electromagnetic particle-in-cell code, though for the simulations presented
here, it has been limited to 1D in configuration-space.

We consider planar hydrogen plasma shocks with protons _p_ _+_ and electrons
_e_ _-_ , with the realistic mass ratio $m_{p}/m_{e}=1836$. The computations
are performed in the frame of the shock. We initialize the computations with
the Rankine-Hugoniot relations as upstream and downstream boundary conditions,
with a hyperbolic tangent transition in between. In order to allow room for
the downstream equilibration zone, the compression-shock position is off-
centered, at $x/L=0.4$, where $L$ is the total domain length. Unless otherwise
specified, the total simulation domain is 250 downstream ion-ion mean-free-
paths, i.e. $L=250\lambda_{ii}^{DS}$. In the regime considered here, the Debye
length $\lambda_{D}$ is much smaller than collisional mean-free-paths,
$\lambda_{D}\ll\lambda_{ii}$, everywhere.

The VPIC simulations were initialized by sampling from the steady-state iFP
ion and electron distribution functions over a truncated spatial domain
($L=130\lambda_{ii}^{DS})$ and run to steady state. To mitigate the noise
inherent in the particle-in-cell simulations, we have applied a Savitzky-Golay
filter (Savitzky, A.; Golay, 1964) to smooth the resulting VPIC moment
profiles. In all the figures, the vertical dotted lines indicate (from left to
right) the approximate locations of the upstream edge of the preheat layer,
the compression shock, and the downstream edge of the relaxation layer. The
spatial profiles presented here are truncated to a span of
$130\lambda_{ii}^{DS}$ to match the VPIC simulation domain.

## Results:

We consider a strong, $M=6$ shock. In Fig. 1-left and -center, we see the
spatial profiles of temperature and temperature difference ($T_{e}-T_{i}$) for
ions and electrons comparing hybrid (kinetic ions and fluid electrons) to
fully kinetic simulations. The profiles of the bulk velocity and density
moments do not change appreciably between hybrid and fully kinetic
simulations. In these moment profiles, the main difference between hybrid and
fully kinetic results is the suppression of the $\delta T=T_{e}-T_{i}$ in the
preheat layer ($x\sim 15$ – $55$) by approximately $\sim 20\%$, with $T_{e}$
relaxing towards $T_{i}$ (see. Fig. 1-center). The temperature profiles
through the compression shock and into the downstream equilibration layer are
essentially unchanged between the fluid and kinetic electron treatments.

In Fig. 1-right we compare the electron Knudsen number based upon the
temperature-gradient length scale [$Kn_{e,\nabla
T}\equiv\lambda_{e}\left|\nabla\left(\ln T_{e}\right)\right|$] for fully
kinetic and hybrid iFP simulations. We see that $Kn_{e}$ is approximately
$5\times 10^{-2}$ throughout the preheat layer (well above the $Kn_{e}\sim
10^{-3}$ threshold), which corresponds exactly to where the electron-
temperature profile has been adjusted in Fig. 1, and is much smaller
elsewhere. Notably the sharp peak in fluid-electron $Kn_{e}$ profile at the
preheat layer upstream edge (where $\partial_{x}T_{e}$ is largest) is smoothed
significantly in the kinetic-electron case.

Figure 2: Top: Spatial profiles within the preheat region for the
$\partial_{x}\mathscr{Q}_{\parallel,e}/T_{e}$ term in Eq. (3) for hybrid
(Hyb.) and fully-kinetic (Kin.) simulations. Bottom: Spatial profiles of the
inverse gradient length scales based on $T_{e}$ (solid) and $n_{e}$ (dashed)
for hybrid (Hyb.) and fully-kinetic (Kin.) simulations.  Figure 3: Left:
Spatial profiles of the electron heat flux comparing
$\mathscr{Q}_{\parallel,e}^{hyb.}$ (Hyb.), $\mathscr{Q}_{\parallel,e}^{kin.}$
(Kin.), and $\mathscr{Q}_{\parallel,e,Brag.}^{kin.}$ (Brag.) and
$\mathscr{Q}_{\parallel,e,LMV}^{kin.}$ (LMV) heat fluxes from the fully
kinetic electron moments. Note the similarity of
$\mathscr{Q}_{\parallel,e}^{hyb.}$ and $\mathscr{Q}_{\parallel,e}^{kin.}$.
Center: Spatial profiles of $\mathscr{Q}_{\parallel,i}^{kin.}$ (Kin.) and
$\mathscr{Q}_{\parallel,e,Brag.}^{kin.}$ (Brag.). Note
$\mathscr{Q}_{\parallel,i}^{kin.}$ is comparable in magnitude to
$\mathscr{Q}_{\parallel,e}^{kin.}$. Right: Phase-space plot of the electron
heat flux comparing $\mathscr{Q}_{\parallel,e}^{kin.}$ (Kin.) to the
$\mathscr{Q}_{\parallel,e}^{hyb.}$ (Hyb.) and the
$\mathscr{Q}_{\parallel,e,LMV}^{kin.}$ (LMV). The heat flux curves for
$\mathscr{Q}_{\parallel,e}^{kin.}$ and $\mathscr{Q}_{\parallel,e,LMV}^{kin.}$
have symbols colored by local Knudsen numbers.

To look deeper into differences between hybrid and fully kinetic simulations,
we examine the balance of terms in the electron energy equation (see Eq. (A3)
of Ref. (Simakov and Molvig, 2016)). Our simulations are performed in the
frame of the shock and allowed to reach a steady state, so the temporal term
is neglected. Further, our analysis showed the contribution of the electron
viscosity, friction, and thermal relaxation terms to be negligible. Upon
rearranging, we find the energy equation may be expressed as

$\frac{\partial_{x}\mathscr{Q}_{\parallel,e}}{T_{e}}\approx
n_{e}u_{\parallel,e}\left[\partial_{x}\left(\ln
n_{e}\right)-\frac{3}{2}\partial_{x}\left(\ln T_{e}\right)\right].$ (3)

As the particle flux density $n_{e}u_{\parallel,e}$ is constant through the
shock, Eq. (3) clearly expresses that the ratio of the divergence of the
electron heat flux to the electron temperature is dependent only on the two
inverse length scales, $L_{T_{e}}^{-1}\equiv\left|\partial_{x}\left(\ln
T_{e}\right)\right|$ and $L_{n_{e}}^{-1}\equiv\left|\partial_{x}\left(\ln
n_{e}\right)\right|$. This is demonstrated in Fig. 2, which shows the left-
hand side of Eq. (3) (Top) and the inverse gradient length scales for the
electron-temperature and electron-density profiles (Bottom). From this, we can
clearly conclude that the primary effect of the kinetic electrons is to smooth
the steep electron-temperature gradient at the preheat layer edge where the
Knudsen number (based on $T_{e}$) is the largest.

It follows from Fig. 1-right and Fig. 2 that, in the preheat layer,
$T_{e}^{Kin.}\sim T_{e}^{Hyb.}$ and $L_{T_{e}}^{Kin.}\sim 2L_{T_{e}}^{Hyb.}$.
Then, Eq. (3) predicts in this region
$\mathscr{Q}_{\parallel,e}^{Kin.}\approx\mathscr{Q}_{\parallel,e}^{Hyb.}$.
Figure 3-left demonstrates this quite clearly, where we see spatial profiles
of the kinetic electron parallel heat flux,
$\mathscr{Q}_{\parallel,e}^{kin.}$, compared to the Braginskii heat flux in
the hybrid simulation, $\mathscr{Q}_{\parallel,e}^{hyb.}$. Here, we compute
the kinetic heat flux for a species $\alpha$ from its distribution function
$f_{\alpha}$ using

$\mathscr{Q}_{\parallel,\alpha}^{kin.}=\int_{\bm{v}}\frac{1}{2}m_{\alpha}\left(\bm{v}-\bm{u}_{\alpha}\right)^{2}\left(v_{\parallel}-u_{\parallel,\alpha}\right)f_{\alpha}d\bm{v},$
(4)

where $\bm{u}_{\alpha}$ is the $\alpha$-species bulk velocity. For comparison,
we also compute from the kinetic-electron moments a Braginskii heat flux,
$\mathscr{Q}_{\parallel,e,Brag.}^{kin.}$, and a heat flux based from the
kinetic closure LMV model (with electric-field correction (Bendib, Luciani,
and Matte, 1988; Schurtz, Nicolaï, and Busquet, 2000)),
$\mathscr{Q}_{\parallel,e,LMV}^{kin.}$. The LMV model is a nonlocal closure of
the form

$\mathscr{Q}_{\parallel,LMV}=\int_{-\infty}^{+\infty}W\left(x,x^{\prime}\right)\mathscr{Q}_{\parallel,Brag.}\left(x^{\prime}\right)\frac{dx^{\prime}}{a\lambda_{e}\left(x^{\prime}\right)},$
(5)

where $W\left(x,x^{\prime}\right)$ is a phenomenologically-chosen
delocalization kernel for the Braginskii heat flux,
$\lambda_{e}\left(x^{\prime}\right)$ is an appropriate delocalization length
scale (related to the mean-free-path and potentially containing an electric
field correction), and $a$ is an adjustable parameter (Ref. (Schurtz, Nicolaï,
and Busquet, 2000) recommends $a=32$). We note that, as we expected from the
analysis of Eq. (3) and Fig. 2, the kinetic and hybrid heat fluxes
$\mathscr{Q}_{\parallel,e}^{kin.}$ and $\mathscr{Q}_{\parallel,e}^{hyb.}$ are
very similar in overall magnitude and shape. The difference in the profile is
essentially in a ‘smoothing’ of the gradients at the upstream edge of the
preheat layer (essentially the same modification as seen in the electron-
temperature profile). For the post-processed heat-flux models based on the
kinetic electron-moment profiles, we see that the LMV heat flux,
$\mathscr{Q}_{\parallel,e,LMV}^{kin.}$, does very well, while
$\mathscr{Q}_{\parallel,e,Brag.}^{kin.}$ is quite different. For reference,
the kinetic ion heat flux, $\mathscr{Q}_{\parallel,i}^{kin.}$, and the
corresponding Braginskii heat flux, $\mathscr{Q}_{\parallel,i,Brag.}^{kin.}$
are also included in Fig. 3-center. There, we observe that
$\mathscr{Q}_{\parallel,i}^{kin.}$ is of the same order of magnitude as
$\mathscr{Q}_{\parallel,e}^{kin.}$ and is substantially larger than
$\mathscr{Q}_{\parallel,i,Brag.}^{kin.}$.

Despite the similarity between hybrid and kinetic-electron heat fluxes, a
detailed analysis reveals deeper physical differences. In particular, we find
clear evidence of nonlocal transport effects in the kinetic-electron
simulation, as expected from weakly collisional conditions. This is evidenced
in Fig. 3-right, where we show a phase-space plot of the heat flux vs.
$\partial_{x}T_{e}$, with symbols colored according to the local Knudsen
number. We see that the kinetic heat flux demonstrates a significant departure
from the Braginskii model, with the nonlocal kinetic heat flux possessing a
distinct multi-valued dependence on the temperature gradient. Notably, the LMV
heat-flux model also demonstrates a __ very similar nonlocal behavior, though
it appears to over-predict slightly the nonlocality relative to the true
kinetic heat flux. We note the LMV without electric-field correction (i.e.,
the original model of Luciani _et al._ (Luciani, Mora, and Virmont, 1983), not
shown here) over-predicts the heat flux by as much as a factor of 2 (the peak
in Fig. 3-right would be nearer to
$-\mathscr{Q}_{\parallel,e}/T_{e}^{5/2}=30$), though with a similar
qualitative shape.

Figure 4: Top: Spatial profiles of $E_{\parallel}$ from hybrid (Hyb.) and
fully kinetic iFP (Kin.) simulations. Also included is the estimated electric
field from the semi-analytic solution of Ref. (Jaffrin and Probstein, 1964)
(J&P). Bottom: Spatial profiles of the charge separation quantity
$\delta_{E}$, again comparing hybrid and fully kinetic iFP and the semi-
analytic estimate of Ref. (Jaffrin and Probstein, 1964)

The electric field, $E_{\parallel}$, shows further differences between hybrid
and fully kinetic simulations. In the fully kinetic simulation,
$E_{\parallel}$ is solved consistently from Eq. 2, while in the hybrid
simulation it is obtained from the electron momentum equation by neglecting
electron inertia, viscosity, and the electron-ion friction term,

$E_{\parallel}=\frac{1}{q_{e}n_{e}}\partial_{x}\left[n_{e}T_{e}\right].$ (6)

In Fig. 4-top, we see the spatial profiles of $E_{\parallel}$ for hybrid and
fully kinetic iFP simulations. As expected, the most noticeable difference is
in the preheat layer, where the electron-temperature gradient has been
significantly reduced, suppressing the $E_{\parallel}$ ‘spike’ there. This
preheat layer spike corresponds to the ’precursor’ electric shock layer
described by Jaffrin and Probstein (Jaffrin and Probstein, 1964). Included in
Fig. 4-top is a profile based on their semi-analytic shock solution. We note
that the derivation of Ref. (Jaffrin and Probstein, 1964) is based on a
Navier-Stokes model. To make a more fair quantitative comparison, we have
altered the solution by making the transformations
$\epsilon\rightarrow\epsilon/2$ in their Eq. (3.4a) and $\epsilon\rightarrow
2\epsilon$ in their Eq. (3.5a) (to bring the ratio of electron and ion
transport coefficients more into line with the estimates of Braginskii
(Braginskii, 1965)). The result shows remarkably good agreement with the
hybrid iFP solution in the precursor shock layer at the upstream preheat layer
edge, but agrees only qualitatively elsewhere. We include an analysis of the
charge separation, $\delta_{E}=\epsilon_{0}\partial_{x}E_{\parallel}$, in Fig.
4-bottom. Given that $\lambda_{D}/\lambda_{ii}\ll 1$, we see that,
unsurprisingly, the charge separation for this problem is quite small. We also
again see the precursor shock collapse in the fully kinetic case. However, the
charge separation predicted by the semi-analytic solution in the preheat layer
is within an order of magnitude of that shown by the hybrid simulation.
Reference (Jaffrin and Probstein, 1964) estimates the precursor shock
thickness to be $l\sim M\lambda_{ii}^{US}/\sqrt{m_{e}/m_{i}}$ (here $l\sim
6$), which is about a factor of two smaller than the width estimated from the
hybrid and fully kinetic simulations.

## Conclusions:

We have performed the first converged fully kinetic simulations of strong
steady-state planar plasma shocks. We have found that the differences between
hybrid and fully kinetic simulations are limited to the electron preheat
region of the shock, and focused at the upstream edge, where the temperature
gradient length scale is smallest and the electron temperature Knudsen number
is largest. Due to its strong connection to the electron-temperature profile
[through Eq. (3)], we find the electron heat flux exhibits stiffness across
various models, with only small adjustments to the heat flux slope near the
upstream edge of the preheat layer. However, while its spatial profile is
qualitatively unchanged, the heat flux exhibits significant nonlocality in the
fully kinetic case. We find that a heat flux computed using the LMV model
(Luciani, Mora, and Virmont, 1983; Bendib, Luciani, and Matte, 1988; Schurtz,
Nicolaï, and Busquet, 2000) using kinetic-electron moments recaptures the
kinetic-electron heat flux reasonably well including the nonlocal phase-space
behavior, while the corresponding Braginskii heat flux shows significant
differences. The electric field, $E_{\parallel}$, at the preheat layer leading
edge is found to be very similar to that predicted by a semi-analytic fluid
solution (Jaffrin and Probstein, 1964), but with minimal charge separation due
to the gradient smoothing effect of the kinetic electrons.

###### Acknowledgements.

This research was performed under the auspices of the National Nuclear
Security Administration of the U.S. Department of Energy at Los Alamos
National Laboratory (LANL) under Contract No. 89233218CNA000001. And was
supported by funding from the Advanced Simulation and Computing Program (ASC)
Thermonuclear Burn Initiative (TBI) and Laboratory Directed Research and
Development (LDRD), with the use of LANL Institutional Computing (IC)
resources. The authors thank Dr. Bill Daughton for assistance with setting up
the VPIC particle-in-cell simulations.

## References

  * Rosenberg  _et al._ (2015) M. J. Rosenberg, F. H. Séguin, P. A. Amendt, S. Atzeni, H. G. Rinderknecht, N. M. Hoffman, A. B. Zylstra, C. K. Li, H. Sio, M. Gatu Johnson, J. A. Frenje, R. D. Petrasso, V. Y. Glebov, C. Stoeckl, W. Seka, F. J. Marshall, J. A. Delettrez, T. C. Sangster, R. Betti, S. C. Wilks, J. Pino, G. Kagan, K. Molvig, and A. Nikroo, Physics of Plasmas 22, 062702 (2015).
  * Rosenberg  _et al._ (2014) M. J. Rosenberg, H. G. Rinderknecht, N. M. Hoffman, P. A. Amendt, S. Atzeni, A. B. Zylstra, C. K. Li, F. H. Séguin, H. Sio, M. G. Johnson, J. A. Frenje, R. D. Petrasso, V. Y. Glebov, C. Stoeckl, W. Seka, F. J. Marshall, J. A. Delettrez, T. C. Sangster, R. Betti, V. N. Goncharov, D. D. Meyerhofer, S. Skupsky, C. Bellei, J. Pino, S. C. Wilks, G. Kagan, K. Molvig, and A. Nikroo, Physical Review Letters 112, 185001 (2014).
  * Rinderknecht  _et al._ (2014) H. G. Rinderknecht, H. Sio, C. K. Li, A. B. Zylstra, M. J. Rosenberg, P. Amendt, J. Delettrez, C. Bellei, J. A. Frenje, M. Gatu Johnson, F. H. Séguin, R. D. Petrasso, R. Betti, V. Y. Glebov, D. D. Meyerhofer, T. C. Sangster, C. Stoeckl, O. Landen, V. A. Smalyuk, S. Wilks, A. Greenwood, and A. Nikroo, Physical Review Letters 112, 135001 (2014).
  * Rinderknecht  _et al._ (2018a) H. G. Rinderknecht, H. S. Park, J. S. Ross, P. A. Amendt, S. C. Wilks, J. Katz, N. M. Hoffman, G. Kagan, E. L. Vold, B. D. Keenan, A. N. Simakov, and L. Chacón, Physics of Plasmas 25, 056312 (2018a).
  * Rinderknecht  _et al._ (2018b) H. G. Rinderknecht, H. S. Park, J. S. Ross, P. A. Amendt, D. P. Higginson, S. C. Wilks, D. Haberberger, J. Katz, D. H. Froula, N. M. Hoffman, G. Kagan, B. D. Keenan, and E. L. Vold, Physical Review Letters 120, 095001 (2018b).
  * Taitano  _et al._ (2018) W. T. Taitano, A. N. Simakov, L. Chacón, and B. D. Keenan, Physics of Plasmas 25 (2018), 10.1063/1.5024402.
  * Simakov  _et al._ (2017) A. N. Simakov, B. D. Keenan, W. T. Taitano, and L. Chacón, Physics of Plasmas 24, 092702 (2017).
  * Jukes (1957) J. D. Jukes, Journal of Fluid Mechanics 3, 275 (1957).
  * Shafranov (1957) V. D. Shafranov, Soviet Physics JETP 5, 1183 (1957).
  * Jaffrin and Probstein (1964) M. Y. Jaffrin and R. F. Probstein, Physics of Fluids 7, 1658 (1964).
  * Grewal (1973) M. S. Grewal, Physics of Fluids 16, 561 (1973).
  * Masser, Wohlbier, and Lowrie (2011) T. O. Masser, J. G. Wohlbier, and R. B. Lowrie, Shock Waves 21, 367 (2011).
  * Zel’dovich and Raizer (1967) Y. B. Zel’dovich and Y. P. Raizer,  _Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena_ , 2nd ed., edited by W. D. Hayes and R. F. Probstein, Vol. II (Academic Press, 1967).
  * Braginskii (1965) S. I. Braginskii, in  _Reviews of Plasma Physics_ , Vol. 1 (Consultants Bureau, 1965) pp. 205–311.
  * Spitzer and Härm (1953) L. Spitzer and R. Härm, Physical Review 89, 977 (1953).
  * Casanova, Larroche, and Matte (1991) M. Casanova, O. Larroche, and J. P. Matte, Physical Review Letters 67, 2143 (1991).
  * Larroche (1993) O. Larroche, Physics of Fluids B 5, 2816 (1993).
  * Vidal  _et al._ (1993) F. Vidal, J. P. Matte, M. Casanova, and O. Larroche, Physics of Fluids B: Plasma Physics 5, 3182 (1993).
  * Keenan  _et al._ (2017) B. D. Keenan, A. N. Simakov, L. Chacón, and W. T. Taitano, Physical Review E 96, 053203 (2017).
  * Keenan  _et al._ (2018) B. D. Keenan, A. N. Simakov, W. T. Taitano, and L. Chacón, Physics of Plasmas 25, 032103 (2018).
  * Luciani, Mora, and Virmont (1983) J. F. Luciani, P. Mora, and J. Virmont, Physical Review Letters 51, 1664 (1983).
  * Bendib, Luciani, and Matte (1988) A. Bendib, J. F. Luciani, and J. P. Matte, Physics of Fluids 31, 711 (1988).
  * Schurtz, Nicolaï, and Busquet (2000) G. P. Schurtz, P. D. Nicolaï, and M. Busquet, Physics of Plasmas 7, 4238 (2000).
  * Cao, Moses, and Delettrez (2015) D. Cao, G. Moses, and J. Delettrez, Physics of Plasmas 22, 082308 (2015).
  * Vidal  _et al._ (1995) F. Vidal, J. P. Matte, M. Casanova, and O. Larroche, Physics of Plasmas 2, 1412 (1995).
  * Matte and Virmont (1982) J. P. Matte and J. Virmont, Physical Review Letters 49, 1936 (1982).
  * Zhang  _et al._ (2021) W. S. Zhang, H. B. Cai, B. Du, D. G. Kang, S. Y. Zou, and S. P. Zhu, Physical Review E 103, 023213 (2021).
  * Huba (2013) J. D. Huba,  _Office of Naval Research_ (2013) pp. 1–71.
  * Anderson  _et al._ (2020) S. E. Anderson, W. T. Taitano, L. Chacón, and A. N. Simakov, Journal of Computational Physics 419, 109686 (2020).
  * Taitano, Chacón, and Simakov (2018) W. T. Taitano, L. Chacón, and A. N. Simakov, Journal of Computational Physics 365, 173 (2018).
  * Taitano  _et al._ (2021a) W. T. Taitano, B. D. Keenan, L. Chacón, S. E. Anderson, H. R. Hammer, and A. N. Simakov, Computer Physics Communications 263, 107861 (2021a).
  * Taitano  _et al._ (2021b) W. T. Taitano, L. Chacón, A. N. Simakov, and S. E. Anderson, Computer Physics Communications 258, 107547 (2021b).
  * Bowers  _et al._ (2008a) K. J. Bowers, B. J. Albright, L. Yin, B. Bergen, and T. J. Kwan, Physics of Plasmas 15, 055703 (2008a).
  * Bowers  _et al._ (2008b) K. J. Bowers, B. J. Albright, B. Bergen, L. Yin, K. J. Barker, and D. J. Kerbyson, in  _SC ’08: Proceedings of the 2008 ACM/IEEE conference on Supercomputing_, November 2008 (2008) pp. 1–11.
  * Bowers  _et al._ (2009) K. J. Bowers, B. J. Albright, L. Yin, W. Daughton, V. Roytershteyn, B. Bergen, and T. Kwan, Journal of Physics: Conference Series 180, 012055 (2009).
  * Savitzky, A.; Golay (1964) M. J. E. Savitzky, A.; Golay, Anal. Chem 36, 1627 (1964).
  * Simakov and Molvig (2016) A. N. Simakov and K. Molvig, Physics of Plasmas 23, 032115 (2016).