# General simulation method
for quantum-sensing systems

Felix Riexinger1,2∗, Mirco Kutas1,2, Björn Haase1,2,
Michael Bortz1, and Georg von Freymann1,2

1Fraunhofer Institute for Industrial Mathematics ITWM
Fraunhofer-Platz 1, Kaiserslautern, 67663, Germany
2Department of Physics and Research Center OPTIMAS, Technische Universität
Kaiserslautern
Kaiserslautern, 67663, Germany

∗To whom correspondence should be addressed; E-mail:
<EMAIL_ADDRESS>
> Quantum sensing encompasses highly promising techniques with diverse
> applications including noise-reduced imaging, super-resolution microscopy as
> well as imaging and spectroscopy in challenging spectral ranges. These
> detection schemes use biphoton correlations to surpass classical limits or
> transfer information to different spectral ranges. Theoretical analysis is
> mostly confined to idealized conditions. Therefore, theoretical predictions
> and experimental results for the performance of quantum-sensing systems
> often diverge. Here we present a general simulation method that includes
> experimental imperfections to bridge the gap between theory and experiment.
> We develop a theoretical approach and demonstrate the capabilities with the
> simulation of aligned and misaligned quantum-imaging experiments. The
> results recreate the characteristics of experimental data. We further use
> the simulation results to improve the obtained images in post-processing. As
> simulation method for general quantum-sensing systems, this work provides a
> first step towards powerful simulation tools for interactively exploring the
> design space and optimizing the experiment’s characteristics.

## Introduction

Quantum sensing has become a rapidly evolving field with developments fueled
by the outlook of applications that surpass classical limitations (?). These
range from noise-reduced imaging (?, ?, ?) over super-resolution microscopy
(?, ?, ?) to imaging and spectroscopy with undetected photons (?, ?, ?, ?, ?).
Especially one novel quantum-imaging technique – quantum imaging with
undetected photons (QIUP) – which has been recently experimentally
demonstrated provides many potential applications (?, ?, ?). Using QIUP, an
image of an object is obtained without detecting the light that interacted
with it. This concept has been adapted to applications in spectroscopy (?, ?)
and optical coherence tomography (?, ?, ?). The principle of these measuring
techniques is based on the quantum superposition of the biphoton states of two
spatially separated sources (?, ?). The light that interacts with the object
is never detected. Meanwhile the light that does not interact with the object
is used to form an image of that object. The registered image depends on the
state of the undetected light due to the superposition of the biphoton states.
This distinctive feature of QIUP allows for a probing wavelength at which
detection is difficult, while detecting the image at a convenient wavelength
(?).

Some properties of QIUP have been analyzed analytically and numerically (?, ?,
?, ?). However, each work investigated only isolated aspects of QIUP, while a
holistic description considering the interplay of various effects is still
missing to date. A detailed understanding is of great importance for the
assessment of limiting factors and the potential of this technique. A full
model of the system further allows for the exploration of trade-offs between
different imaging characteristics such as resolution and visibility. The
availability of detailed simulations is becoming more relevant since
applications are on the verge of a breakthrough, while their characteristic
properties still need to be evaluated experimentally. Using simulations, the
design process can be sped up as the system’s properties can be explored and
optimized virtually.

In this article, we present and demonstrate a method for simulating the
detected images of a QIUP setup. This method includes a full model of the
biphoton sources and the optical system. With this, a large range of
experimental setups including their misaligned states can be simulated. We
further present an approximative variant of the method which speeds up the
computation times substantially. For this variant, the model of biphoton
sources is simplified. The approximative variant is still applicable to a
large range of experimental setups but is limited to fewer types of
misalignments. We demonstrate our method by simulating detector images and
imaging properties of a real experimental setup. With this, we present a
simulation suite enabling the virtual analysis of quantum-sensing experiments
for the first time.

Figure 1: Layout of the experimental setup. The dichroic mirrors separate the
wavelengths such that the detected light never interacts with the object. Only
the components included in the model are displayed. NL: nonlinear crystal. L:
lens. (D)M: (dichroic) mirror. BS: beam splitter. $d_{\mathrm{t}}$ denotes a
transversal displacement of the lens. The blue lines denote the source, object
and detector planes. Length of optical paths are not to scale.

QIUP schemes with a Mach-Zehnder geometry (Fig. 1) are driven by two
coherently pumped identical nonlinear crystals (NL1 and NL2) which produce
biphoton states $\ket{\psi_{\mathrm{NL1/NL2}}}$ consisting of a signal (s) and
an idler (i) photon. Schemes with a Michelson geometry use one nonlinear
crystal through which the pump laser passes twice. We treat the first and
second pass as two sources and denote them with NL1 and NL2 to keep the
notation consistent. We assume the power of the laser to be small enough, such
that stimulated emission can be neglected (?). The signal parts of the
biphoton states emitted from each source (yellow paths) are superposed and
subsequently detected with a camera. The idler part (red paths) emitted by NL1
is propagated through NL2 and aligned with the corresponding beam emitted by
NL2. In an ideal setup the probability amplitudes and modes emitted by NL1 and
NL2 are identical. Therefore, the idler states are indistinguishable and no
information about the emitting source can be obtained in the detector plane.
Due to the entanglement of the photons s and i, the observed interference on
the camera depends on the idler path even though only signal photons are
detected (?, ?, ?). To measure an object’s properties, it is placed in the
path of the undetected idler beam between the two crystals. The transmission
coefficient and the phase introduced by the object can then be reconstructed
from the detected image (?).

## Results

### Theory of quantum imaging

Our simulation model consists of three parts: the source model, the
propagation model and the detection model. The source model describes the
quantum states created by the biphoton sources. The propagation model
calculates the propagation of these states through the optical system. The
detection part evaluates the interference of the two paths and the resulting
signal in the detector plane. In the following, we will limit ourselves to the
description of QIUP setups in a far-field configuration.

The source model is used for the calculation of the transition amplitude
$g(\mathbf{k}_{\mathrm{s}},\mathbf{k}_{\mathrm{i}})$ for each pair of photon
momenta $\mathbf{k}_{\mathrm{s}}$ (signal) and $\mathbf{k}_{\mathrm{i}}$
(idler). The quantum state of light generated by the sources NL1 and NL2 is
given as

$\ket{\psi_{\mathrm{NL1/NL2}}}=\int\mathrm{d}\mathbf{k}_{\mathrm{s}}\mathrm{d}\mathbf{k}_{\mathrm{i}}g_{1/2}(\mathbf{k}_{\mathrm{s}},\mathbf{k}_{\mathrm{i}})\ket{\mathbf{k}_{\mathrm{s}},\mathbf{k}_{\mathrm{i}}}$
(1)

where $g_{1/2}(\mathbf{k}_{\mathrm{s}},\mathbf{k}_{\mathrm{i}})$ are the
transition amplitudes for the nonlinear crystals NL1 and NL2, respectively. We
omit the polarization component of the biphoton state for simplicity. Many
biphoton sources used for QIUP emit mainly one combination of photon
polarizations.

Figure 2: Sketch illustrating the notation at different parts of the setup.
Panels (A) and (C) show the origin of the biphoton states. The relation of
$\mathbf{k}_{\mathrm{i1}}$ to a point $\rho^{\prime}$ in the object plane is
illustrated in panel (B). The set of $\textbf{k}_{\mathrm{s}}$-vectors in (A)
and (C) is chosen such that the positions $\rho$ match in the detector plane
as shown in panel (D). The k-vectors are deliberately not parallel before the
lenses to demonstrate the case of a misaligned system denoted by the distances
${\sim}f_{\mathrm{i}}$.

Next we consider the propagation, which encompasses the change of momenta and
positions, losses, and phase-shifts of the biphoton states. We model losses
due to propagation and from the imaged object as beam splitters with
transmittance $t$ and a phase shift $\phi$ such that they can be described by
the transformation $\hat{{a}}^{\dagger}\rightarrow
t\mathrm{e}^{i\phi}\hat{{a}}^{\dagger}+\sqrt{1-t^{2}}\hat{{a}}^{\dagger}_{\mathrm{loss}}$.
Here, $\hat{{a}}^{\dagger}$ are the photon creation operators and
$\hat{{a}}^{\dagger}_{\mathrm{loss}}$ represents all losses along a beam path
or from the object. In general, the transmittances and phase shifts depend on
the photon momenta. We omit this dependence in the following equations for the
sake of readability. The propagation of the states
$\ket{\psi_{\mathrm{NL1/NL2}}}$ through the beam paths from a source plane to
the detectors plane is modeled by the transformation of the photon momenta
$\mathbf{k}_{m}\rightarrow\mathbf{k}_{m}^{\mathrm{D}}=h_{m}(\mathbf{k}_{m})$
for $m=\mathrm{i1,i2,s1,s2}$. Here $h_{m}$ is a transformation function
modeling the changed photon momenta after the propagation from the source
planes P1 and P2 to the detector plane D. The definition of the planes is
shown in Fig. 2.

In an ideal setup, a plane wave emerging from one of the crystals is focused
onto a point on the detector (?). This point $\rho(\mathbf{k}_{\mathrm{s}})$
depends on the momentum vector $\mathbf{k}_{\mathrm{s}}$ which defines the
plane wave. This is illustrated in Fig. 2D. The same holds for the relation of
momentum vector $\mathbf{k}_{\mathrm{i}}$ and a point on the object
$\rho^{\prime}(\mathbf{k}_{\mathrm{i}})$ (Fig. 2B). For a more realistic
setup, these assumptions do not hold as the components are not perfect and
further might be misaligned. In order to still have a momentum-to-point
relation, we make the assumption that the spatial extent of the sources can be
neglected.

With this, we can establish a relation between the photons from the two
sources (Figs. 2A and C). Fig. 2 illustrates these relations for a Mach-
Zehnder setup as well as the employed notation. A biphoton state emitted from
NL1 with signal momentum $\mathbf{k}_{\mathrm{s1}}$ is detected at the point
$\rho(\mathbf{k}_{\mathrm{s1}})$ on the detector. For a biphoton state from
NL2 to be detected at the same point the signal part needs to have the
momentum $\mathbf{k}_{\mathrm{s2}}^{\rho}$ in the source plane NL2 which is
the momentum that fulfills
$\rho(h_{\mathrm{s2}}(\mathbf{k}_{\mathrm{s2}}^{\rho}))=\rho(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}})$.
We can establish an analogous relation for the idler photons. With these
relations from the propagation model we can formulate the detection model. The
quantum state at the detector plane can be expressed as

$\displaystyle\ket{\psi_{\mathrm{D}}}=\int\mathrm{d}\mathbf{k}_{\mathrm{s1}}$
$\displaystyle\mathrm{d}\mathbf{k}_{\mathrm{i1}}\left[g_{1}(\mathbf{k}_{\mathrm{s1}},\mathbf{k}_{\mathrm{i1}})t_{\mathrm{s1}}t_{\mathrm{i}}t_{\mathrm{o}}\mathrm{e}^{i\Delta\phi_{1}}\hat{{a}}^{\dagger}(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}})\hat{{a}}^{\dagger}(\mathbf{k}_{\mathrm{i1}}^{\mathrm{D}})\right.$
$\displaystyle+g_{1}(\mathbf{k}_{\mathrm{s1}},\mathbf{k}_{\mathrm{i1}})t_{\mathrm{s1}}\sqrt{1-(t_{\mathrm{i}}t_{\mathrm{o}})^{2}}\mathrm{e}^{i\Delta\phi_{2}}\hat{{a}}^{\dagger}(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}})\hat{{a}}^{\dagger}_{\mathrm{i,loss}}$
$\displaystyle+g_{1}(\mathbf{k}_{\mathrm{s1}},\mathbf{k}_{\mathrm{i1}})t_{\mathrm{i}}t_{\mathrm{o}}\sqrt{1-t_{\mathrm{s1}}^{2}}\mathrm{e}^{i\Delta\phi_{3}}\hat{{a}}^{\dagger}_{\mathrm{s,loss}}\hat{{a}}^{\dagger}(\mathbf{k}_{\mathrm{i1}}^{\mathrm{D}})$
$\displaystyle+\left.g_{2}(\mathbf{k}_{\mathrm{s2}}^{\rho},\mathbf{k}_{\mathrm{i2}}^{\rho})t_{\mathrm{s2}}\mathrm{e}^{i\Delta\phi_{4}}\hat{{a}}^{\dagger}(\mathbf{k}_{\mathrm{s2}}^{\mathrm{D}})\hat{{a}}^{\dagger}(\mathbf{k}_{\mathrm{i2}}^{\mathrm{D}})\right]\ket{0}$
(2)

where the $\Delta\phi_{n}$ denote the phase shifts from the propagation
through the various optical paths which are provided in detail in Eqs. 10 to
13, $t_{\mathrm{o}}$ is the transmittance of the imaged object,
$t_{\mathrm{s1/s2}}$ are the transmittances of the signal paths between source
and detector planes and $t_{\mathrm{i}}$ is the transmittance of the idler
path between the two source planes. The first three terms describe the
biphoton states emitted from crystal NL1. Starting with the case where both
photons are transmitted through the system and then the cases where only
signal or only idler photons are transmitted. The fourth term describes
biphotons from the second crystal. We omitted further terms which do not
contribute to the measured signal such as idler losses after NL2.

From the state at the detector we obtain the count rate for a point on the
detector plane as

$\displaystyle\Gamma^{\pm}(\rho(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}}))$
$\displaystyle=\bra{\psi_{\mathrm{D}}}\left[\hat{{a}}^{\dagger}(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}})+\hat{{a}}^{\dagger}(\mathbf{k}_{\mathrm{s2}}^{\mathrm{D}})\right]\left[\hat{a}(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}})+\hat{a}(\mathbf{k}_{\mathrm{s2}}^{\mathrm{D}})\right]\ket{\psi_{\mathrm{D}}}$
(3)
$\displaystyle=\int\mathrm{d}\mathbf{k}_{\mathrm{i1}}\left[g_{1}(\mathbf{k}_{\mathrm{s1}},\mathbf{k}_{\mathrm{i1}})^{2}t_{\mathrm{s1}}^{2}+g_{2}(\mathbf{k}_{\mathrm{s2}}^{\rho},\mathbf{k}_{\mathrm{i2}}^{\rho})^{2}t_{\mathrm{s2}}^{2}\right.$
$\displaystyle\left.\quad\quad\pm
2g_{1}(\mathbf{k}_{\mathrm{s1}},\mathbf{k}_{\mathrm{i1}})g_{2}(\mathbf{k}_{\mathrm{s2}}^{\rho},\mathbf{k}_{\mathrm{i2}}^{\rho})t_{\mathrm{s1}}t_{\mathrm{s2}}t_{\mathrm{i}}t_{\mathrm{o}}\cos(\Delta\phi)\right]$
(4)

where $\Delta\phi=\Delta\phi_{1}-\Delta\phi_{4}$ is the phase difference of
biphoton states created in NL1 and NL2. In general, the phase terms depend on
the photon momenta and therefore the cosine term cannot be extracted from the
integral. The different signs of the third term stem from the $\pi$ phase
shift the beam splitter introduces to one of its output ports. To obtain the
count rate for a detector pixel we integrate $\Gamma^{\pm}$ over the area
covered by this pixel and multiply with the detector efficiency.

From the measurements at the two outputs (see Fig. 1), we can calculate the
pointwise visibility which gives an approximation of the optical properties of
the object in the idler path:

$\displaystyle V(\rho(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}}))=$
$\displaystyle[\Gamma^{+}(\rho(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}}))-\Gamma^{-}(\rho(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}}))]/[\Gamma^{+}(\rho(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}}))+\Gamma^{-}(\rho(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}}))]$
(5) $\displaystyle=$
$\displaystyle\frac{\int\mathrm{d}\mathbf{k}_{\mathrm{i}}2g_{1}(\mathbf{k}_{\mathrm{s1}},\mathbf{k}_{\mathrm{i1}})g_{2}(\mathbf{k}_{\mathrm{s2}}^{\rho},\mathbf{k}_{\mathrm{i2}}^{\rho})t_{\mathrm{s1}}t_{\mathrm{s2}}t_{\mathrm{i}}t_{\mathrm{o}}\cos(\Delta\phi)}{\int\mathrm{d}\mathbf{k}_{\mathrm{i}}[g_{1}(\mathbf{k}_{\mathrm{s1}},\mathbf{k}_{\mathrm{i1}})^{2}t_{\mathrm{s1}}^{2}+g_{2}(\mathbf{k}_{\mathrm{s2}}^{\rho},\mathbf{k}_{\mathrm{i2}}^{\rho})^{2}t_{\mathrm{s2}}^{2}]}$
(6)

where we have assumed identical imaging systems for the two measurements. In
this general case, the relation of this quantity and the object properties is
not straightforward. For an idealized setup with perfect momentum correlation,
the idler momentum integrals vanish which yields $V\propto
t_{\mathrm{o}}\cos(\phi_{\mathrm{o}})$ (?). We can make some approximations to
simplify Eq. 6 while still retaining important properties of the system. For
this we assume the crystals to be identical ($g=g_{1}=g_{2}$), constant and
equal losses in the optical paths
($t_{\mathrm{s}}=t_{\mathrm{s1}}=t_{\mathrm{s2}}$), vanishing phases
($\Delta\phi=\phi_{\mathrm{o}}$), and an optical system that ensures
$\mathbf{k}_{\mathrm{s1}}=\mathbf{k}_{\mathrm{s2}}^{\rho}$ and
$\mathbf{k}_{\mathrm{i1}}=\mathbf{k}_{\mathrm{i2}}^{\rho}$. This allows us to
reduce Eq. 6 to:

$V(\rho(\mathbf{k}_{\mathrm{s1}}^{\mathrm{D}}))=\frac{\int\mathrm{d}\mathbf{k}_{\mathrm{i}}g(\mathbf{k}_{\mathrm{s1}},\mathbf{k}_{\mathrm{i1}})^{2}t_{\mathrm{s}}^{2}t_{\mathrm{i}}t_{\mathrm{o}}\cos(\phi_{\mathrm{o}})}{\int\mathrm{d}\mathbf{k}_{\mathrm{i}}g(\mathbf{k}_{\mathrm{s1}},\mathbf{k}_{\mathrm{i1}})^{2}t_{\mathrm{s}}^{2}}$
(7)

where only the transition probabilities and object properties $t_{\mathrm{o}}$
and $\phi_{\mathrm{o}}$ depend on the photon momenta.

We further simplify and consider biphoton sources where the normalized
transition probabilities
$\hat{g}=g(\mathbf{k}_{\mathrm{s}},\mathbf{k}_{\mathrm{i}})^{2}/\int\mathrm{d}\mathbf{k}_{\mathrm{i}}g(\mathbf{k}_{\mathrm{s}},\mathbf{k}_{\mathrm{i}})^{2}$
depend only on the difference in $\mathbf{k}_{\mathrm{s}}$ and
$\mathbf{k}_{\mathrm{i}}$, i.e.
$\hat{g}=\hat{g}(\mathbf{k}_{\mathrm{s}}-\mathbf{k}_{\mathrm{i}})$. In the
case of a purely absorptive object ($\phi_{\mathrm{o}}=0$) this reduces Eq. 7
to a convolution of the objects transmittance with the $\hat{g}$ kernel:

$V(\rho(\mathbf{k}_{\mathrm{s}}^{\mathrm{D}}))=t_{\mathrm{i}}\int\mathrm{d}\mathbf{k}_{\mathrm{i}}\hat{g}(\mathbf{k}_{\mathrm{s}}-\mathbf{k}_{\mathrm{i}})t_{\mathrm{o}}(\rho^{\prime}(\mathbf{k}_{\mathrm{i}}))$
(8)

where the dependence of the object transmittance is now explicitly given. Due
to the constant convolution kernel this approach yields a computationally fast
method for calculating the measurement for any object. From this, properties
of the imaging setup such as resolution limit and maximum visibility can be
inferred. The employed assumption can readily be applied if the region of
interest is the center of the spots in the detector plane or when the higher
relative noise at lower count rates can be neglected. The assumption of
vanishing phases can be relaxed giving a more general case where the
convolution kernel becomes shift variant as it also contains the cosine term.
Including the signal $\phi_{\mathrm{s}}$ and idler phases $\phi_{\mathrm{i}}$
allows for the emulation of misalignments, such as transversally shifted
lenses in the two optical arms. An improvement in computational efficiency is
retained as long as $\hat{g}$ needs to be evaluated only sparsely. Both the
quasi-Monte Carlo simulation and the convolution approximation can be extended
to a broader spectral range by integrating over $\mathbf{k}_{\mathrm{s1}}$ as
well.

The insight, that the measured image can be seen as a convolution of the
object combined with a way to calculate the convolution kernel, paves the way
for using established deconvolution and image reconstruction techniques. Using
these techniques in the post-processing of the detected images allows to
improve their quality without any improvements in the detection system.
Applied correctly, deconvolution techniques provide a more accurate
approximation of the inferred object properties. To demonstrate this, we use
the Richardson-Lucy algorithm (?, ?) to deconvolve the simulated measurement
results with the convolution kernel $\hat{g}$.

### Implementation

To demonstrate the capabilities of our simulation methods, we model a far-
field Mach-Zehnder interferometer setup presented in (?, ?). The setup’s
sources create biphoton states with central wavelengths
$\lambda_{\mathrm{s}}=810\,$nm and $\lambda_{\mathrm{i}}=1550\,$nm. A detailed
account of the model of the setup is given in the Materials and Methods
section.

To describe the biphoton states produced by the two spontaneous parametric
downconversion (SPDC) sources NL1 and NL2 we use a highly accurate and general
model for predicting absolute photon rates (?). This sophisticated model takes
many of the material properties into account, such as the wavelength
dependence of refractive indices and nonlinear susceptibility yielding an
accurate and complex SPDC model that needs to be evaluated numerically. Note
that the approximations that have been made for the derivation of the SPDC
model do not limit the generality of the presented method, as it can be
substituted with a different source model.

The magnification of the setup depends on the wavelengths of signal and idler
(?, ?), it can be approximated by

$M=\frac{f_{\mathrm{D}}\lambda_{\mathrm{s}}}{f_{\mathrm{i}}\lambda_{\mathrm{i}}}$
(9)

In general, the magnification is different along the crystal axis due to the
birefringence of many nonlinear crystals. For the setup presented in this
paper, these differences can be neglected.

The transformation functions $h_{m}$ are modeled using geometrical optics. The
momentum $\mathbf{k}$ and position $\rho$ in the plane of interest is
evaluated by propagating a ray defined by the photon momentum to the
corresponding plane. This is used to determine the position on the object and
detector, as well as the photon momenta in the plane of the second source
$\mathbf{k}_{\mathrm{s2}}^{\rho}$. The phase terms are determined from the
optical path lengths $l_{m}$ as $\phi_{m}=\omega_{m}l_{m}/c$. For the aligned
interferometer we assume that only the phase shift from the object is
relevant, i.e. $\Delta\phi=\phi_{\mathrm{o}}$. We further assume that no
losses occur in the optical paths ($t_{\mathrm{s}}=t_{\mathrm{i}}=1$), as
these are either constant factors which reduce the visibility. Or they are
complex functions which would diminish the clarity of the demonstration of our
simulation method. For a full simulation of the QIUP setup we evaluate Eq. 3
over the detector area using a quasi-Monte Carlo integration scheme. For the
convolution variant, we evaluate the normalized transition probabilities
$\hat{g}(\mathbf{k}_{\mathrm{s}}-\mathbf{k}_{\mathrm{i}})$ for the collinear
$\mathbf{k}_{\mathrm{s}}$ defined by the optical axis of the system and the
central signal wavelength of the setup $\lambda_{\mathrm{s}}$. This is used to
calculate the convolution of the object given in Eq. 8. The computational
expense of evaluating $\hat{g}$ is equivalent to the evaluation of one pixel
on the detector, reducing the theoretical computation time by six orders of
magnitude for this setup.

Figure 3: Transmission object used for the simulations. The object measures
4$\times 4$ mm2 and consists of horizontal and vertical opaque bars at
different distances.

### Numerical Results

Figure 4: Simulation results. Simulated measurements of the transmission
object shown in Fig. 3. (A) and (B): Detector images for the destructive
($\Gamma^{-})$ and constructive ($\Gamma^{+}$) configuration. The noise stems
from the finite sampling rate and is perfectly correlated between the two
images. (C): The resulting measured transmittance. For the left half Gaussian
noise was added to the detector images, the right half is free from artificial
noise. (D): Constructive detector image for the misaligned setup. The lens
misalignment leads to vertical stripes in the detector plane. The top panel
shows simulated and theoretical modulation from the phase shift at
$\mathrm{y}=50$.

We simulate the detector images with the absorption object shown in Fig. 3. It
measures $4\times 4$ mm2, has three horizontal, $200$ µm wide bars with the
same distance between them, and two vertical bars of the same width but
separated by $1000$ µm. The pump beam waist used for this simulation is $300$
µm. The results of the simulated detector images and measurement of the
transmittance of a test object are shown in Fig. 4. The results for the
detector count rates $\Gamma^{-}$ and $\Gamma^{+}$ are shown in Figs. 4A and
4B, respectively. The two simulated images show the expected behavior, where
the sum of both images gives the intensity profile of the signal radiation
without distinctive features, whereas the difference shows an image of the
object. These figures also show the spot size of the signal radiation. The
noise visible in the simulated detector images is a numerical artifact which
is caused by the finite number of points we evaluated. The noise is perfectly
correlated between the two images such that the resulting pointwise visibility
shows no artifacts from this noise. Fig. 4C shows the measured pointwise
visibility of the object which is defined by Eq. 5. To demonstrate more
realistic conditions, we add Gaussian noise with a standard deviation of $50$
counts to the left half of the detector images for the calculation of this
image. We clip values below $0$ and above $1$. The right half is free from
artificial noise to illustrate the influence of this added noise. The
resulting measured visibility shows a smeared out image of the object. The
visibility in the noise-free half of the image ranges from $0$ to $1$ as we
neglected misalignment and losses from the optical paths and components which
would reduce it. The visibility does not decrease towards the outer parts of
the features. In the noisy half of the image, the visibility ranges from $0.1$
to $1$ and the noise is more prominent in the corners where the count rate is
lower. The measured distance of the centers of the vertical bars is $1255$ µm.
The distance on the object is $1200$ µm (see Fig. 3) which gives a
magnification value of $1.046$, in agreement with the theoretical prediction
from Eq. 9 ($M=1.045$). The accuracy of the measured magnification is limited
by the detector pixel size.

In real experimental setups, noise is not the only issue. Even with included
noise, they will always fall short of this ideal setup. As components have
imperfections and cannot be aligned perfectly, our model needs to be able to
reflect these effects. Due to the flexible propagation model, this is
possible. To demonstrate our method’s capabilities for misaligned systems, we
shift the idler lens Li1 in front of the object by $d_{\mathrm{t}}=0.3$ mm
transversally respective to the optical axis as shown in Fig. 1. This
introduces a phase shift which depends on the emission angle of the idler
photons. Fig. 4D shows the introduced vertical stripes in the simulated
detector images as one would expect from misaligning an interferometer. The
width of the stripes depends on the magnification and the misalignment. An
analytical expression for the phase shift in the detector plane is derived in
Eq. 17. The analytical and simulated results are shown in the top panel of
Fig. 4D. Our simulated stripe pattern matches the theoretical predictions
closely. The pointwise visibility in this setup is decreased to the range from
$0$ to $0.64$ due to the phase variation in the idler states.

We further calculated the resolution limit for both setups numerically (see
supplementary materials) to compare with the measured results of Fuenzalida et
al. (?). Our results show the same characteristic curve, but are consistently
smaller than the measured results for both experimental setups. Our simulated
resolution limits are $\sim$ $17$% smaller. This is to be expected when
comparing an ideal with a real setup. Reasons for the differences are the
common difficulties in real experiments: imperfect overlap of the beam paths,
losses, spherical aberrations, noise and others. Our results present a
theoretical limit to the achievable resolution for these setups and show that
these experimental realizations might still be improved. A detailed analysis
of the experimental setups can be performed using our simulation method. This
allows the identification of sensitive parts and limiting factors, which can
be used to improve the setups.

Figure 5: Comparison of simulation methods. The figure shows the measured
transmittance from the full numerical simulation (blue) and the convolution
method (orange). The perfect image of the object is given as a reference
(black).

With these examples we have demonstrated the capability and versatility of the
presented simulation method. Under certain conditions we can formulate the
pointwise visibility as a convolution, which leads to a huge speedup of the
simulation times while still giving accurate results. To demonstrate this, we
change the beam waist to $200$ µm and use an absorption object that consists
of three vertical slits with transmittance $t=1$. The width and distance
between the slits is $128$ µm which is the resolution limit for the setup in
this configuration. A horizontal cut at the center ($y=500$) of the detector
is shown in Fig. 5. The figure shows the results for the full simulation and
the convolution method. The convolution method yields an accurate
approximation of the fully simulated image. The maximum difference in
visibility over the shown area is $0.007$, showing that this method provides a
fast and accurate approximation for the full simulation results. An extension
of this method to phase objects and misaligned setups and the corresponding
simulations are shown in the supplementary materials. A downside of this
method is that the limited spot size is not considered. For system
characteristics such as visibility and resolution limit the convolution method
provides a fast calculation method that enables the exploration of the
influence of different parameters and the optimization of setups.

Figure 6: Deconvolution results. Shown is the deconvolution of the simulated
visibility of Fig. 4C using the normalized transition probabilities $\hat{g}$
of the system. The image becomes sharper and the noise is notably reduced. The
root mean square error to the real object is improved by $37$% and $35$% in
the noisy and noise-free half.

The convolution kernel used for the convolution variant can also be used to
improve the information we obtain from the two measurements. For this we apply
the Richardson-Lucy deconvolution with the collinear convolution kernel
$\hat{g}$ to the simulated visibility shown in Fig. 4C. The result after 50
iterations, shown in Fig. 6, is a much sharper image with narrower edges and
an increased contrast, especially in the center. Outside of the beam spot the
large relative noise leads to large artifacts. We therefore limit our analysis
to a box defined by $x\in[100,900]$ and $y\in[300,700]$ which cuts off the
parts where noise is dominant. In this box, the deconvolution algorithm
introduces only minor artifacts. With a maximum amplitude of $0.05$, they are
barely visible in the figure. In this area, the noise leads to minor
distortions. The addition of noise does negligibly increase the amplitude of
the artifacts compared to the noise-free half. To compare the deconvolution to
the regular simulation results, we calculate the root mean square error (RMSE)
to the real object for both methods. The deconvolution method improves the
RMSE in the noisy and noise-free half by $37$% and $35$%, respectively. The
range for the pointwise visibility in the center of the image is improved from
[$0.1,\,0.8$] to [$0,\,1$], which is the maximum achievable contrast. In the
noise-free half, the maximum improvement for a single pixel is $0.44$, the
mean improvement on this side is $0.1$. To quantify the improvement in the
sharpness, we have calculated the resolution limits including the
deconvolution step. Compared to the simulations without deconvolution, we
found a consistent improvement of 13% across a large range of beam waists.

## Discussion

We have presented and demonstrated a general method for simulating quantum
imaging with undetected photons. For the first time, these simulations enable
the characterization of QIUP systems in the design phase. We presented a
method for full fledged simulations of QIUP systems which allows analyzing a
large range of setups in full detail, including the influence of many
misalignments. The convolution approximation provides a trade-off between
computation time and accuracy and also a fast way to evaluate characteristics
of a setup with slightly limited options for misalignments. The numerical
simulations enable the identification of limiting factors of QIUP systems. The
fast computation times even allow for an interactive exploration of the design
space and optimization of the imaging characteristics – an invaluable tool for
improving applications based on QIUP. Both methods yield similar results for
perfect imaging systems. A comparison with analytical and experimental results
shows that the simulation results show realistic benchmarks of experimental
setups.

We further demonstrated the improvement of imaging characteristics by
deconvolution techniques. In our simulation experiment, we showed that the
quality of an image can be improved greatly by applying this technique.

Our results demonstrate the principle of the methods. Using more accurate
optical models to reproduce the images from experimental setups is subject of
further research. The convolution technique can also be improved. By
incorporating multiple convolution kernels, it is possible to get a better
approximation for the variation in different detector regions, allowing for an
estimation of the absolute count rate as well. Further work is also possible
in applying the deconvolutional technique to real experimental setups and
data. This requires an accurate model of the setup and its imperfections.
Using a deconvolution algorithm for shift-variant kernels allows adding the
phase information to the deconvolution method which can improve the image
further, since in real setups phase errors always occur.

The presented approach can also be applied to a wide range of other
experimental setups. We have presented an approach specifically for collinear
spectroscopy experiments before (?). The application to a Michelson geometry
or a near-field configuration can be realized by adapting the transfer
functions (?) or by using position-space operators instead of momentum-space
operators (?), respectively. The simulation method we present can be further
extended to include the spatial component of the biphoton sources. However,
this requires a different source model and is beyond the scope of this proof-
of-concept paper. At the current stage, the presented model and methods
already provide the necessary basis for the simulation of experiments and
applications based on quantum imaging with undetected photons and enables
obtaining improved images.

## Materials and Methods

### Phase terms

The phase terms $\Delta\phi_{n}$ in Eq. Theory of quantum imaging are composed
of the phase shifts from the propagation of the biphoton states through the
optical system. These can be separated into individual components $\phi_{m}$.
Where the indices $m=\mathrm{p1,p2}$ denote the propagation of the pump laser
to the source planes, $m=\mathrm{s1,i1,s2,i2}$ denote the phase shift of
propagation from source to detector plane and ${\mathrm{o}}$ denotes the phase
shift introduced by the object. The phase terms are then defined as

$\displaystyle\Delta\phi_{1}$
$\displaystyle=\phi_{\mathrm{p1}}+\phi_{\mathrm{s1}}+\phi_{\mathrm{i1}}+\phi_{\mathrm{o}}$
(10) $\displaystyle\Delta\phi_{2}$
$\displaystyle=\phi_{\mathrm{s1}}+\xi_{\mathrm{i1}}+\xi_{\mathrm{o}}$ (11)
$\displaystyle\Delta\phi_{3}$
$\displaystyle=\xi_{\mathrm{s1}}+\phi_{\mathrm{i1}}+\phi_{\mathrm{o}}$ (12)
$\displaystyle\Delta\phi_{4}$
$\displaystyle=\phi_{\mathrm{p2}}+\phi_{\mathrm{s2}}+\phi_{\mathrm{i2}}$ (13)

where $\xi_{m}$ stand for phase shifts which depend on the position where the
photon gets absorbed. In our model the precise form of $\xi_{m}$ is not
further specified as these terms do not influence the detected signal. Using
the form of these terms we find the phase shift modulating the image

$\Delta\phi=\phi_{\mathrm{o}}+\phi_{\mathrm{p1}}+\phi_{\mathrm{s1}}+\phi_{\mathrm{i1}}-\phi_{\mathrm{p2}}-\phi_{\mathrm{s2}}-\phi_{\mathrm{i2}}$
(14)

which shows that we observe constructive interference when $\Delta\phi$ is an
even multiple of $\pi$ and destructive interference when it is an odd
multiple.

### Experiment models

The setup models we evaluate are based on the experimental setups presented by
Fuenzalida et al. (?). These setups have a Mach-Zehnder geometry and share the
same optical components. The crystals used are periodically poled potassium
titanyl phosphate (PPKTP) crystals of the dimensions $2\times 2\times 1$ mm3
(L$\times$W$\times$H) with a type-0 phase-matching scheme. The setups differ
in the crystal poling periods $\Lambda_{1/2}$ and pump lasers which results in
different signal and idler wavelengths. The first setup, which is used for the
results in the main text, contains PPKTP crystals with $\Lambda_{1}=9.675$ µm
which are pumped by a $532$ nm laser. The second setup contains crystals with
$\Lambda_{2}=5.33$ µm which are pumped by a 405 nm laser. The first setup’s
sources create biphoton states with central wavelengths
$\lambda_{\mathrm{s}}=810$ nm and $\lambda_{\mathrm{i}}=1550$ nm. The second
setup creates states with central wavelengths $\lambda_{\mathrm{s}}=842$ nm
and $\lambda_{\mathrm{i}}=780\,$nm. The temperatures of the crystals are
$85^{\circ}$C in the first setup and $75^{\circ}$C in the second setup. Both
lasers are modeled with a power of $50$ mW. We choose a detector pixel size of
$5\times 5$ µm2. All lenses are modeled as ideal lenses defined only by their
focal lengths. The focal lengths of the lenses are $f_{\mathrm{i}}=75$ mm for
Li1, Li2, Ls1, and Ls1. The lens LD in front of the detector has a focal
length of $f_{\mathrm{D}}=150$ mm. The theoretical value for the magnification
is $M=2.159$ for setup 2. As the focal lengths are the same in both setups,
the difference stem from the ratio of wavelengths. The measured distance of
the centers of the vertical bars is $2590$ µm, which gives a magnification of
$2.158$ in good agreement with the theoretical value.

### Phase shift in misaligned setup

The setup is misaligned by a transversal shift of lens Li1. To obtain an
analytical expression for the resulting phase shift in the detector plane we
calculate the path-length difference with ray optics. We start with the phase
shift from a transversal displacement of a perfect lens in the parabolical
approximation. The phase shift from a lens aligned on the optical axis is
given by

$\phi_{\mathrm{L}}(x_{\mathrm{L}})=\frac{\pi x_{\mathrm{L}}^{2}}{f\lambda}$
(15)

where $x_{\mathrm{L}}$ is the position where the ray passes through the lens,
$f$ is the focal length of the lens and $\lambda$ is the wavelength of the
light. Assuming that the imaging system is aligned before the lens
displacement, only the additional phase introduced by the lens shift is
relevant. This is given by

$\displaystyle\Delta\phi_{\mathrm{L}}(x_{\mathrm{L}})$
$\displaystyle=\phi_{\mathrm{L}}(x_{\mathrm{L}}+{\mathrm{t}})-\phi_{\mathrm{L}}(x_{\mathrm{L}})$
(16) $\displaystyle=\frac{\pi{\mathrm{t}}^{2}+2\pi
x_{\mathrm{L}}{\mathrm{t}}}{f\lambda}$ (17)

where ${\mathrm{t}}$ is the transverse lens shift. To obtain the phase shift
in dependence of the position in the detector plane, we use the relation
$x_{\mathrm{D}}=Mx_{\mathrm{L}}$ which is valid for the two lenses Li1 and Li2
in the idler path.

## References

  * 1. M. Gilaberte Basset, F. Setzpfandt, F. Steinlechner, E. Beckert, T. Pertsch, M. Gräfe, Perspectives for applications of quantum imaging, Laser & Photonics Rev 13, 1900097 (2019).
  * 2. G. Brida, M. Genovese, I. Ruo Berchera, Experimental realization of sub-shot-noise quantum imaging, Nat. Photonics 4, 227-230 (2010).
  * 3. N. Samantaray, I. Ruo-Berchera, A. Meda, M. Genovese, Realization of the first sub-shot-noise wide field microscope, Light-Sci. Appl. 6, e17005-e17005 (2017).
  * 4. J. Sabines-Chesterking, A. R. McMillan, P. A. Moreau, S. K. Joshi, S. Knauer, E. Johnston, J. G. Rarity, J. C. F. Matthews, Twin-beam sub-shot-noise raster-scanning microscope, Opt. Express 27, 30810–30818 (2019).
  * 5. A. Classen, J. von Zanthier, M. O. Scully, G. S. Agarwal, Superresolution via structured illumination quantum correlation microscopy, Optica 4, 580–587 (2017).
  * 6. M. Unternährer, B. Bessire, L. Gasparini, M. Perenzoni, A. Stefanov, Super-resolution quantum imaging at the heisenberg limit, Optica 5, 1150–1154 (2018).
  * 7. R. Tenne, U. Rossman, B. Rephael, Y. Israel, A. Krupinski-Ptaszek, R. Lapkiewicz, Y. Silberberg, D. Oron, Super-resolution enhancement by quantum image scanning microscopy, Nat. Photonics 13, 116-122 (2019).
  * 8. G. B. Lemos, V. Borish, G. D. Cole, S. Ramelow, R. Lapkiewicz, A. Zeilinger, Quantum imaging with undetected photons, Nature 512, 409–412 (2014).
  * 9. M. Gilaberte Basset, A. Hochrainer, S. Töpfer, F. Riexinger, P. Bickert, J. R. León-Torres, F. Steinlechner, M. Gräfe, Video-rate imaging with undetected photons, Laser & Photonics Rev. 15, 2000327 (2021).
  * 10. J. Fuenzalida, A. Hochrainer, G. B. Lemos, E. A. Ortega, R. Lapkiewicz, M. Lahiri, A. Zeilinger, Resolution of Quantum Imaging with Undetected Photons, Quantum 6, 646 (2022).
  * 11. C. Lindner, J. Kunz, S. J. Herr, S. Wolf, J. Kießling, F. Kühnemann, Nonlinear interferometer for fourier-transform mid-infrared gas spectroscopy using near-infrared detection, Opt. Express 29, 4035–4047 (2021).
  * 12. M. Kutas, B. Haase, J. Klier, D. Molter, G. von Freymann, Quantum-inspired terahertz spectroscopy with visible photons, Optica 8, 438–441 (2021).
  * 13. D. A. Kalashnikov, A. V. Paterova, S. P. Kulik, L. A. Krivitsky, Infrared spectroscopy with visible light, Nat. Photonics 10, 98–101 (2016).
  * 14. A. V. Paterova, H. Yang, C. An, D. A. Kalashnikov, L. A. Krivitsky, Tunable optical coherence tomography in the infrared range using visible photons, Quantum Sci. Technol. 3, 025008 (2018).
  * 15. A. Vanselow, P. Kaufmann, I. Zorin, B. Heise, H. M. Chrzanowski, S. Ramelow, Frequency-domain optical coherence tomography with undetected mid-infrared photons, Optica 7, 1729–1736 (2020).
  * 16. M. Kutas, B. Haase, P. Bickert, F. Riexinger, D. Molter, G. von Freymann, Terahertz quantum sensing, Sci. Adv. 6, eaaz8065 (2020).
  * 17. L. J. Wang, X. Y. Zou, L. Mandel, Induced coherence without induced emission, Phys. Rev. A 44, 4614–4622 (1991).
  * 18. X. Y. Zou, L. J. Wang, L. Mandel, Induced coherence and indistinguishability in optical interference, Phys. Rev. Lett. 67, 318–321 (1991).
  * 19. M. Kutas, B. E. Haase, F. Riexinger, J. Hennig, P. Bickert, T. Pfeiffer, M. Bortz, D. Molter, G. von Freymann, Quantum sensing with extreme light, Adv. Quantum Technol. 5, 2100164 (2022).
  * 20. M. Lahiri, R. Lapkiewicz, G. B. Lemos, A. Zeilinger, Theory of quantum imaging with undetected photons, Phys. Rev. A 92, 013832 (2015).
  * 21. M. Lahiri, A. Hochrainer, R. Lapkiewicz, G. B. Lemos, A. Zeilinger, Nonclassicality of induced coherence without induced emission, Phys. Rev. A 100, 053839 (2019).
  * 22. W. H. Richardson, Bayesian-based iterative method of image restoration, J. Opt. Soc. Am. 62, 55–59 (1972).
  * 23. L. B. Lucy, An iterative technique for the rectification of observed distributions, Astron. J. 79, 745 (1974).
  * 24. F. Riexinger, M. Kutas, B. Haase, P. Bickert, D. Molter, M. Bortz, G. von Freymann, General simulation method for spontaneous parametric down- and parametric up-conversion experiments, arXiv preprint arXiv:2112.07243 (2021).
  * 25. B. Viswanathan, G. B. Lemos, M. Lahiri, Position correlation enabled quantum imaging with undetected photons, Opt. Lett. 46, 3496–3499 (2021).

## Acknowledgments

#### Funding:

This project was funded by the Fraunhofer-Gesellschaft within the Fraunhofer
Lighthouse Project Quantum Methods for Advanced Imaging Solutions (QUILT).

#### Author contributions:

G.v.F. and M.B. initiated this research. F.R. performed the theoretical and
numerical analysis. All authors discussed the results and contributed to the
writing of the manuscript.

#### Competing interests:

The authors declare that they have no competing interests.

#### Data and materials availability:

All data needed to evaluate the conclusions in the paper are present in the
paper. Additional data related to this paper may be requested from the
authors.

## Supplementary materials

### Resolution limit

We use the definition of the resolution limit given in ref. (?). This
definition uses an object with three slits of width $d$, separated by the same
length. We define $R$ as the ratio of visibility values at the local minima
between slits and the maximum at the center (see Fig. 5) and the minimum
resolvable distance as the slit width $d_{\mathrm{limit}}$ for which $R$ is
equal to $0.81$. For the numerical calculation of the resolution limit, we
first define a routine to calculate the ratio $R$. For this, we evaluate the
collinear convolution kernel of a given setup. This is used to calculate the
pointwise visibility of slits with width $d$ using the convolution variant of
our method. From the resulting image we calculate the visibility ratio $R$.
The minimum resolvable distance $d_{\mathrm{limit}}$ is then computed by an
optimization algorithm.

Figure S1: Resolution limits. Experimental and theoretical resolution limits
for both setup 1 and setup 2. The experimentally measured values (?) are
larger than the simulation of the ideal setups. The deconvolution improves the
resolution limit by $13$% for both setups.

This process is repeated with various pump beam waists $w_{\mathrm{p}}$ for
both setup 1 and 2. In a similar fashion, we evaluate the resolution limit
when a deconvolution is applied to the image before calculating $R$. The
resulting resolution limits are shown in Fig. S1. The experimental results are
limited not only by measurement accuracy, but also because the used test
target only provides discrete slit widths. The simulated results provide
continuously calculated values.

### Extension of convolution approximation

The cosine term in Eq. 7 vanishes under the assumption of an ideal setup and
an absorption object. The convolution method can still be applied when these
assumptions are relaxed. For this, we model the phase shift from the object
and misalignments affecting the phase. An example of such a misalignment is
defined by Eq. 17. The phase shifts in this model depend on the locations and
momenta at certain stages in the system. The locations can be determined from
the momenta in the source planes. The resulting term for the visibility is
similar to Eq. 7, but with the phase term $\phi_{\mathrm{o}}$ replaced by

$\Delta\phi=\sum_{m}\phi_{m}(\mathbf{k}_{m},x_{m})$ (18)

where $\phi_{m}$ are the phase shifts introduced by the object and the optical
system, $\mathbf{k}_{m}$ are the momenta and $x_{m}$ are the transverse
positions at the corresponding locations. The resulting equation can be
interpreted as a convolution with a shift-variant kernel, i.e. for each point
$\rho(\mathbf{k}_{\mathrm{s}})$ in the detector plane a different kernel has
to be used. This slightly increases the complexity of the computations, but
the resulting images can still be calculated within seconds.

Figure S2: Comparison of simulation methods for misaligned system. The
simulated visibility for the convolution approximation (A) and the quasi-Monte
Carlo simulation (B).The transversely shifted lens leads to the vertical
stripes in the measured visibility. The maximum visibility is reduced to
$0.63$ due to the misalignment.

We apply this method to the misaligned setup used for Fig. 4D. The resulting
simulation of the measured visibility is shown in Fig. S2. The image closely
resembles the result from the quasi-Monte Carlo simulation. The difference of
the two simulations is characterized by a RMSE of $0.013$ and a maximum error
of $0.07$ over all pixels. The differences are larger than in the case without
phases. This is explained by the more dynamic image caused by the stripes. Due
to this, numerical inaccuracies such as in the magnification lead to larger
errors.


capbtabboxtable[][]

# X-ray phase-contrast micro tomography of soft tissues using a compact
laboratory system with two-directional sensitivity

Carlos Navarrete-León  Department of Medical Physics and Biomedical
Engineering. University College London, Malet Place, Gower Street, WC1E 6BT
London, United Kingdom  Adam Doherty  Department of Medical Physics and
Biomedical Engineering. University College London, Malet Place, Gower Street,
WC1E 6BT London, United Kingdom  Savvas Savvidis  Department of Medical
Physics and Biomedical Engineering. University College London, Malet Place,
Gower Street, WC1E 6BT London, United Kingdom  Mattia F. M. Gerli  UCL
Division of Surgery and Interventional Science, Royal Free Hospital, London,
United Kingdom  Stem Cell and Regenerative Medicine Section, Great Ormond
Street Institute of Child Health, University College London, WC1N 1EH, United
Kingdom  Giovanni Piredda  Research Center for Microtechnology, Vorarlberg
University of Applied Sciences, Hochschulstr. 1, 6850, Dornbirn, Austria
Alberto Astolfo  Department of Medical Physics and Biomedical Engineering.
University College London, Malet Place, Gower Street, WC1E 6BT London, United
Kingdom  David Bate  Nikon X-Tek Systems Ltd, Tring, Herts, HP23 4JX, United
Kingdom  Department of Medical Physics and Biomedical Engineering. University
College London, Malet Place, Gower Street, WC1E 6BT London, United Kingdom
Silvia Cipiccia  Department of Medical Physics and Biomedical Engineering.
University College London, Malet Place, Gower Street, WC1E 6BT London, United
Kingdom  Charlotte K. Hagen  Department of Medical Physics and Biomedical
Engineering. University College London, Malet Place, Gower Street, WC1E 6BT
London, United Kingdom  Alessandro Olivo  Department of Medical Physics and
Biomedical Engineering. University College London, Malet Place, Gower Street,
WC1E 6BT London, United Kingdom  Marco Endrizzi  Department of Medical Physics
and Biomedical Engineering. University College London, Malet Place, Gower
Street, WC1E 6BT London, United Kingdom

###### Abstract

X-ray micro tomography is a non-destructive, three-dimensional inspection
technique applied across a vast range of fields and disciplines, ranging from
research to industrial, encompassing engineering, biology and medical
research. Phase-contrast imaging extends the domain of application of X-ray
micro tomography to classes of samples that exhibit weak attenuation, thus
appear with poor contrast in standard X-ray imaging. Notable examples are low-
atomic-number materials, like carbon-fibre composites, soft matter and
biological soft tissues. We report on a compact and cost effective system for
X-ray phase-contrast micro tomography. The system features high sensitivity to
phase gradients and high resolution, requires a low-power sealed X-ray tube, a
single optical element, and fits in a small footprint. It is compatible with
standard X-ray detector technologies: single-photon-counting offers higher
sensitivity whereas flat-panels are preferred for a larger field of view. The
system is benchmarked against known-material phantoms and its potential for
soft-tissue three-dimensional imaging is demonstrated on small-animal organs:
a piglet oesophagus and a rat heart.

## Introduction

X-ray micro tomography has become an invaluable tool for the non-destructive
volumetric imaging of samples for a wide variety of fields, ranging from
medical sciences to metrology and manufacturing [1, 2, 3]. Despite continuous
progress in X-ray micro tomography technology, three dimensional imaging of
low-density samples, such as soft tissue specimens, remains challenging due to
the weak X-ray attenuation contrast which is linked to low atomic number
materials. X-ray phase-contrast micro tomography (XPC$\mu$T) can overcome this
limitation by exploiting also the phase-shift experienced by the wave when
traversing the object for generating image contrast [4], thus extending the
applicability of micro tomography to a broader range of samples and
disciplines.

The technology for X-ray phase-contrast imaging is more demanding in
comparison to what is needed for conventional X-ray imaging, specifically in
terms of the required spatial and temporal coherence of the radiation [4].
Whilst synchrotron radiation facilities can achieve high degrees of coherence,
that is not often the case for laboratory-based systems, where the power
density for X-ray generation is limited and a small footprint is required. In
turn, this means that one has to choose between having a small or a powerful
source, and working with monochromatic radiation often requires long exposure
times. Many solutions have been identified and pursued over the past twenty
five years to overcome these limitations and facilitate the translation of
X-ray phase-contrast imaging to a laboratory setting. Without compiling a
complete history, here we summarise some of the milestones: microfocus source
and crystal combination [5], free space propagation with polychromatic
radiation [6], Talbot-Lau interferometry [7, 8, 9], Edge Illumination [10,
11], Fourier analysis [12], liquid metal jet sources [13], universal Moiré
effect [14], speckle-based imaging also in combination with liquid metal jet
sources [15, 16, 17, 18] and implicit tracking [19]; we refer to some
excellent reviews on the topic for a more comprehensive list and discussion
[20, 21, 22, 23, 24, 25, 26]. Notable laboratory-based XPC$\mu$T applications
include: imaging of soft-tissue specimens like lung [27, 28], breast [29],
heart [30], oesophagus [31] and brain [32], multi-material phantoms [33, 34,
35] and for composites like carbon fibre-reinforced composite materials [36].

In this context, we have identified a solution for lab-XPC$\mu$T that has
minimum requirements on the experimental set-up. Our approach is based on a
low-power, fixed target and sealed X-ray tube, it requires a single optical
element (modulator) and it is compatible with flat-panel detectors as well as
with more sophisticated single-photon counters. All of these components are
readily available and require a minimal amount of maintenance. The modulator
is fabricated through laser ablation of a tungsten foil and does not need to
be positioned with high accuracy or repeatability. The total length of the
system is less than one meter, making it one of the most compact for lab-
XPC$\mu$T. Our approach shares similarities with a Shack-Hartmann wavefront
sensor and in our work it evolved from one-directional to two-directional
sensitivity of beam tracking [37, 38], it can also be viewed as a highly
parallelised version of a single-probe scanning system [39]. We call this
approach two-directional beam tracking (2DBT), this solution offers a compact,
robust and cost-effective way of producing multi-contrast three-dimensional
images of the inner structure of samples, without compromising on sensitivity.

## Materials and methods

Figure 1: System’s set-up and working principle. (a) A schematic of the
laboratory set-up, including a sealed microfocus X-ray tube, the modulator,
sample stage and X-ray detector. An experimental image of a PMMA calibration
sphere with an inset comparing beamlets without and with the sample in place
(b). The effects of the presence of the sample are barely detectable by eye
inspection, however the dedicated data analysis algorithm reliably extracts
refraction images (c, d). These are subsequently integrated leading to phase
images (e) that offer superior signal-to-noise ratio to the standard
attenuation-contrast image (f). The refraction images are retrieved with a
sub-pixel cross-correlation algorithm, we observe good agreement between the
phase shift expected from theory, obtained through the numerical simulation
and from experimental data (f). All scale bars are $500$ $\mu$m.

### Model

The core element of the imaging system (Figure 1a) is the structured
illumination, obtained by spatially modulating the amplitude of the radiation
field before it reaches the sample. Structuring the illumination in two
dimensions enables tracking small changes in the position of each probe, i.e.
measuring refraction along two orthogonal directions. It redefines sampling
and spatial resolution of the imaging system. Sampling is no longer governed
by the detector pixel pitch, it is instead defined by the distance between
apertures in the modulator. The spatial resolution is redefined as being equal
to or better than the width of the apertures in the modulator [40, 41]. The
main drawback of this approach is the sacrifice in flux available for imaging.
We note, however, that this sacrifice does not affect dose efficiency because
the radiation field is shaped before it reaches the sample. We also note that
sampling is not necessarily limited by the pitch of the modulator: finer
sampling can be achieved by recombining subsequent exposures, having moved
either the modulator or the sample by a fraction of the pitch. We increase
sampling by displacing the modulator by equally spaced sub-pitch increments.
This requires a reference image for each modulator position, however it
preserves the conventional cone beam geometry and allows using standard 3D
reconstruction algorithms, even for the data sets acquired with higher
sampling. Each probe is detected, and subsequently analysed, independently.
The X-ray intensity distribution of the probe is modelled through a series of
convolutions:

$I_{0}(x,y)=(S*M*PSF)(x,y)$ (1)

where $*$ denotes the two-dimensional convolution operator, $S$ is the source
intensity distribution, geometrically scaled by
$-1+(z_{sm}+z_{mo}+z_{od})/z_{sm}$, $M$ the transmission of the modulator
(unity over a circle representing the aperture and zero elsewhere)
geometrically scaled by $(z_{sm}+z_{mo}+z_{od})/z_{sm}$, and $PSF$ the
detector point spread function. Referring to Figure 1a, $z_{sm}$ is the
distance between the source and the modulator, $z_{mo}$ is the distance
between the modulator and the axis of rotation and $z_{od}$ is the distance
between the axis of rotation and the detector.

When a sample is in place, the intensity distribution of the probes is
attenuated and shifted by the local properties of the sample:

$I(x,y)=tI_{0}(x-\Delta x,y-\Delta y)$ (2)

where $t$ is the transmission through the sample (at the position of the
probe) and $\Delta x$ and $\Delta y$ are the probe displacements on the
detector plane due to the refraction induced by the sample. These shifts are
measured by comparing the probe intensity distributions with and without the
object (Figure 1b). The transmission is calculated as the ratio of the total
intensity in the probe:

$t=\frac{\sum I(x,y)}{\sum I_{0}(x,y)},$ (3)

where sum operates over all the pixels illuminated by a single probe (Figure
1b, zoom-ins). For the retrieval of the displacements $\Delta x$ and $\Delta
y$ we used a sub-pixel cross-correlation algorithm [42]. These displacements
are related to the refraction angle and to the object-to-detector distance,
$z_{od}$, by:

$\alpha_{x}(x,y)=\tan^{-1}\left(\frac{\Delta
x}{z_{od}}\right)\;,\;\alpha_{y}(x,y)=\tan^{-1}\left(\frac{\Delta
y}{z_{od}}\right).$ (4)

By linking the two refraction images to the orthogonal gradients of the phase
shift (Figure 1, panels c and d) we obtain:

$\alpha_{x}(x,y)=\frac{1}{k}\Delta\Phi_{x}(x,y)\;,\;\alpha_{y}(x,y)=\frac{1}{k}\Delta\Phi_{y}(x,y)$
(5)

where $k$ is the wave-number. The phase-shift of the wavefront
$\Delta\Phi(x,y)$ (Figure 1e) can be obtained by means of the Fourier
derivative theorem [43]:

$\Delta\Phi(x,y)=\mathcal{F}^{-1}\left[\frac{\mathcal{F}[\Delta\Phi_{x}+i\Delta\Phi_{y}](u,v)}{2\pi
i(u+iv)}\right](x,y)$ (6)

where $(u,v)$ represent the reciprocal space coordinates to $(x,y)$. A better
signal-to-noise ratio can be seen by direct visual comparison with the
attenuation-image contrast of the same phantom (Figure 1f). The retrieved
quantities $t$ and $\Delta\Phi$ are line integrals along the photon path of
two physical properties of the sample, the linear attenuation coefficient
($\mu$) and the real part of the refractive index ($\delta$):

$-\ln t(x,y)=\int_{o}\mu(x^{\prime},y^{\prime},z^{\prime})dz$ (7)

$-\frac{\Delta\Phi(x,y)}{k}=\int_{o}\delta(x^{\prime},y^{\prime},z^{\prime})dz$
(8)

these relationships allow one to obtain volumetric reconstructions of $\mu$
and $\delta$ from projections taken at different viewing angles using the
standard algorithms for tomography.

### Experiments

The X-ray system features a fixed W-target microfocus Hamamatsu X-ray source
(L12161-07) which we run at 40 kVp and 10 W. The modulators were manufactured
with a laser ablation procedure from readily available 100 $\mu$m-thick
Tungsten foils (Goodfellow). Two pitches, 50 and 100 $\mu$m, were tested. The
apertures have a conical shape with diameters of $\sim$15 $\mu$m and $\sim$30
$\mu$m in the front and back apertures, respectively, the narrower aperture
was facing the source.

Three detectors were tested: i) Pixirad-2/PIXIE-III photon-counter with a 650
$\mu$m CdTe sensor (98% detection efficiency up to 50 keV), 62 $\mu$m pixel
pitch and a $50\times 32$ mm2 active area; ii) MerlinX photon-counter (Quantum
Detectors) with a Medipix 3RX chip (256x256), a 500 $\mu$m Si sensor and 55
$\mu$m pixel pitch; and iii) Hamamatsu Flat Panel (C9732DK-11) with directly
deposited CsI scintillator, 50 $\mu$m pixel pitch and $12\times 12$ cm2 active
area. The geometry was defined as follows: $z_{sm}=(150,150,140)$ mm,
$z_{mo}=(30,30,25)$ mm and $z_{od}=(680,580,535)$ mm, for the Pixirad, MerlinX
and Hamamatsu detectors, respectively. Single-photon-counters are preferable
when sensitivity is the priority, while flat panels offer a much larger field
of view albeit with lower sensitivity.

To characterise the system and fine-tune the quantitative retrieval algorithm,
a test phantom consisting of a Poly(methyl methacrylate) (PMMA) calibration
sphere (Goodfellow, $\diameter=3.18\pm 0.05$mm) was used. For the experimental
assessment of spatial resolution we used a phantom composed of soda-lime glass
micro-spheres (Fischer Scientific, monodisperse, 50 $\mu$m diameter) deposited
on a kapton substrate. The detector was the Pixirad with the 100 $\mu$m-pitch
modulator and noise threshold set at 10 keV. A total of $16\times 16$ frames
were acquired for both phantoms, with a sampling step in both $x$ and $y$
directions of $6.25$ $\mu$m and 1 s exposure time per frame.

The angular sensitivity was quantified on planar images of a phantom composed
of soda-lime glass micro-spheres (Fischer Scientific, monodisperse, 50 $\mu$m
diameter) embedded in wax and polyethylene foam. Images were acquired at
$4\times 4$ positions, with $x$ and $y$ steps of $12.5$ $\mu$m. The detector
was the MerlinX with the 50 $\mu$m-pitch modulator and noise threshold set at
5 keV. For each step 256 frames of 0.25 s exposure time were acquired to study
the sensitivity as a function of exposure time. The angular sensitivity was
measured in an area without the sample by calculating the mean and standard
error of the standard deviation of the measured refraction angles [44], 9
different windows of $10\times 10$ pixels each were used.

The ability of the system to measure $\mu$ and $\delta$ in three dimensions
was assessed with a multi-material phantom composed of three spheres
(Goodfellow) in a plastic straw. The materials were Polypropylene (PP),
Polystyrene (PS) and PMMA with diameters: $\diameter_{PP,PMMA}=3.18\pm 0.05$
mm, $\diameter_{PS}=3.5\pm 0.1$ mm. The detector was the Pixirad with the 100
$\mu$m-pitch modulator and noise threshold set at 8 keV. We acquired 900
projections over 360∘ with 1 s exposure time per projection. A full rotation
was acquired for each modulator position, giving a total of $8\times 8$ scans
with sampling every 12.5 $\mu$m in $x$ and $y$.

The performance of the system for soft-tissue imaging was demonstrated with
two animal-derived small organs; a heart and an oesophagus, extracted from a
300 g Sprague-Dawley rat and a 3 kg piglet respectively. The heart sample was
obtained from the UCL BSU, whilst the oesophagus from piglets procured by JSR
genomics, a Home Office approved supplier. Animals were euthanised via
Schedule 1 methods and sample sizing has been implemented following the NC3R
principles. The organs were prepared for imaging through paraformaldehyde
(PFA) fixation, followed by critical point drying (as described in previous
work [31]). For the piglet oesophagus we acquired 900 equally spaced
projections over 360∘ with 1 s exposure time using the Pixirad detector. The
$100\,\mu$m-pitch modulator was stepped by 12.5 $\mu$m in $x$ and $y$ for a
total of $8\times 8$ scans and an exposure time of $900\times 8\times 8\times
1$ s $=16$ h. For rat heart we acquired 2000 equally spaced projections over
360∘ with 1.2 s exposure time using the Hamamatsu flat-panel. The
$50\,\mu$m-pitch modulator was stepped $4\times 4$ times on a grid with in
$12.5$ $\mu$m displacements in $x$ and $y$. The total exposure time was
$2000\times 4\times 4\times 1.2$ s $=10.7$ h.

All reconstructions were performed using the Feldkamp-Davis-Kress algorithm
implementation in CUDA of the Astra Toolbox [45]. The volume renderings
presented in Figure 5 were obtained using the Avizo software (Thermo Fischer
Scientific).

Figure 2: System characterisation and benchmarking. (a) Resolution and (c)
accuracy of phase retrieval, as a function of sampling, in images of the
monodisperse soda-lime glass micro-spheres on Kapton substrate and the PMMA
calibration sphere, respectively. (a) The ability to visualise small details
depends critically on a fine sampling and is mainly defined by the width
aperture in the modulator. The MTF (b) matches well a real-space full width
half maximum of $16$ $\mu$m. (c) Phase is retrieved accurately in all cases,
with higher noise linked to coarser sampling. This can also be seen in the
line profiles plot (d) where the phase shifts retrieved in the four
experimental cases are compared to the theoretical phase shift. (e) Images of
the angular sensitivity phantom: the system sensitivity depends on the
exposure time, and longer exposure times lead to better signal to noise
ratios. Exposures as short as 250 ms are sufficient for robust phase
integration. Longitudinal slices of (f) $\mu$ and (g) $\delta$ for the multi-
material phantom. (h) The ability to separate the three materials is greatly
improved by the simultaneous availability of the two contrast channels.

Exp time (s) |  $\sigma(\alpha_{x})$ $(\mu rad)$ |  $\sigma(\alpha_{y})$ $(\mu rad)$
---|---|---
0.25 | 2.14 $\pm$ 0.05 | 2.33 $\pm$ 0.07
0.5 | 1.62 $\pm$ 0.04 | 1.66 $\pm$ 0.03
1 | 1.10 $\pm$ 0.02 | 1.19 $\pm$ 0.03
2 | 0.82 $\pm$ 0.02 | 0.85 $\pm$ 0.02
4 | 0.61 $\pm$ 0.01 | 0.65 $\pm$ 0.01
8 | 0.43 $\pm$ 0.01 | 0.46 $\pm$ 0.01
16 | 0.35 $\pm$ 0.01 | 0.34 $\pm$ 0.01
32 | 0.27 $\pm$ 0.01 | 0.31 $\pm$ 0.01
64 | 0.22 $\pm$ 0.01 | 0.22 $\pm$ 0.01

Figure 3: System sensitivity for different exposure times, MerlinX detector
and 50 $\mu$m-pitch modulator.

Figure 4: System sensitivity trend compared to quantum-limited ideal system,
MerlinX detector and 50 $\mu$m-pitch modulator.

## Results and discussion

Figure 2a and 2c report the resolution phantom and the PMMA calibration sphere
at different sampling steps ranging from $50$ to $6.25$ $\mu$m. The
positioning of the modulator has to be accurate and repeatable in comparison
to these numbers, namely approximately half of the aperture width at best,
which is compatible with most micro-positioning technologies. Whilst the phase
shift of a large object is correctly retrieved in all cases, the ability to
represent finer details and distinguish between adjacent micro-spheres
critically depends on the sampling step. The modulation transfer function of
the system at $6.25$ $\mu$m sampling pitch is shown in 2b, with 50% at $20$
lp/mm and 10% at 38 lp/mm. It matches well with the expectation from a real-
space full-width-half-maximum of $16$ $\mu$m, comparable to the aperture size
in the modulator. Figure 2d shows four line plots extracted from the phase
images in panel 2d against the theoretically expected phase shift. There is
good agreement between the phase shifts measured with different sampling steps
and theory, with higher noise linked to lower sampling. The decrease in noise
that is linked to progressively finer sampling is also visible by visual
comparison of the four images in panel 2c.

The angular sensitivity depends on the exposure time and closely follows the
trend expected from the Poisson statistics (Table 4 and Figure ). This
indicates that the main limitation in the phase retrieval process is photon
statistics. The gain in sensitivity with exposure time is slightly worse than
in a perfect quantum-limited system, this is due to environmental
instabilities including temperature drifts and vibrations, but these effects
become noticeable only with tens of seconds of integration. The two
differential phase contrast images and the integrated phase of the sensitivity
phantom are presented in Figure 2e. Smaller details emerge from the background
noise as the exposure time is increased. We report sub-micro-radian
sensitivities starting at 2 s of exposure time and down to minimum of 220
$\pm$ 10 $n$rad with 64 s per frame. These values are obtained with a single
exposure image, thus with a resolution limited by the modulator pitch at 50
$\mu$m. With the methods described earlier, the resolution can be increased up
to approximately $16$ $\mu$m while maintaining the same angular sensitivity
and increasing the total exposure time. The sensitivity with the Hamamatsu
flat-panel was slightly worse, ($\sigma(\alpha_{x})$, $\sigma(\alpha_{y})$) =
(1.96 $\pm$ 0.05, 2.10 $\pm$ 0.05) $\mu$rad and ($\sigma(\alpha_{x})$,
$\sigma(\alpha_{y})$) = (1.31 $\pm$ 0.04, 1.41 $\pm$ 0.02) $\mu$rad for 1 s
and 2 s exposure time, respectively. However, the field of view was extended
to 28 mm $\times$ 28 mm from the 3.2 mm $\times$ 3.2 mm achievable with the
MerlinX. To position our system in the landscape of existing solutions for
lab-based XPC$\mu$T we report in Table 1 a summary of imaging system
parameters from the literature, including source power, spatial resolution,
system dimensions and sensitivity.

Reference Technique Resolution Sensitivity Scan time (s) X-ray spot size
Source power Dimension Diemoz el al [46] EI 12 or 66.8 $\mu$m 270 nrad 112 or
14s 70 $\mu$m 875W 200 cm Havariyoun el al [47] EI 10 or 79 $\mu$m 230 nrad
384 or 24s 70 $\mu$m 400W 100 cm Revol et al [48] GBI 104 $\mu$m 110 nrad
80.4s 1x1 mm2 250-1000W 140 cm Thuring et al [49] GBI 4.1-7.1 $\mu$m 250-550
nrad 128s 5-10 $\mu$m 4-12W 20-33 cm Birnbacher et al [50] GBI 170 $\mu$m 5
nrad 275s 132x226 $\mu$m2 2.8kW 230 cm Villa-Comamala et al [51] GBI 21.5
$\mu$m 45 nrad 250s 10 $\mu$m 60W 90.4 cm Zanette et al [15] SBI 84 $\mu$m 240
nrad 300s 7.8x8.7 $\mu$m2 30W 300 cm Quenot et al [19] SBI 36 $\mu$m 5
$\mu$rad 600s 4 $\mu$m - 53 cm Our set-up 2DBT 15 or 50 $\mu$m 345 nrad 256 or
16s 5-20 $\mu$m 10W 70-86 cm

Table 1: Comparison of laboratory-based XPC$\mu$T systems in literature

In comparison with EI and the first GBI systems, our approach achieves
comparable sensitivity with a low-power source. The highest sensitivity was
achieved with GBI [50] and the most powerful X-ray source, at the lowest
spatial resolution. A more recent lab-based implementation of GBI [51]
achieved high sensitivity at high resolution, and required fine-pitch gratings
with high aspect ratio. SBI [15] also reports high sensitivity and requires a
large footprint and specialised source technology.

Tomographic reconstructions of $\mu$ and $\delta$ of the multi-material
phantom are presented in Figure 2f and 2g. The lower noise of phase tomography
is apparent. This can be seen quantitatively in the histograms in panel 2h
where the values extracted from the boxes highlighted in colours are compared.
Phase contrast allows for a better separation of PP and PS with respect to
attenuation contrast, and the simultaneous availability of both contrast
channels improves material separation even further. Measured values are
compared against theoretically expected ones [52] in Table 2. Effective
energies vary slightly for $\delta$ and $\mu$ as expected [53], however the
range of values is very close to the mean energy of the X-ray spectrum ($15.6$
keV). The contrast-to-noise ration is calculated as:

$CNR=\frac{\mu_{m}-\mu_{air}}{\sigma_{air}}$ (9)

where $\mu_{m}$ is the linear attenuation coefficient of each material and
$\sigma_{air}$ is the standard deviation of the background. An analogous
formula is used for $\delta$. CNRs were found to be substantially higher for
phase tomography.

Material $\mu_{exp}$ (cm-1) $\mu_{theo}$ (cm-1) $E_{\mu}$ (keV) CNR of $\mu$
$\delta_{exp}$ (x10-6) $\delta_{theo}$ (x10-6) $E_{\delta}$ (keV) CNR of
$\delta$ PP 0.76 $\pm$ 0.06 0.76 13 4.1 1.14 $\pm$ 0.01 1.14 14 31.6 PMMA 1.23
$\pm$ 0.06 1.23 15 7.1 1.22 $\pm$ 0.01 1.22 15 34.0 PS 0.86 $\pm$ 0.06 0.85 14
4.7 1.03 $\pm$ 0.01 1.03 15 28.4

Table 2: Experimental and theoretical $\mu$ and $\delta$ of the phantom
materials, along with effective energies and CNRs measured with respect to the
background (yellow box in Figure 2g). Good agreement is observed on both
attenuation- and phase-contrast channels with expectations from theory. Phase
tomography offers consistently higher CNRs.

Lab-XPC$\mu$T of the animal-derived small organs are presented in Figure 5.
The volume renderings and coronal, longitudinal and axial slices of both
samples are shown. Access to the full volume provides structural information
of the specimens which are represented with isotropic spatial resolution of
approximately $16$ $\mu$m.

Figure 5: X-ray phase-contrast micro tomography of a piglet esophagus (a-c)
and a rat heart (d, e). Piglet esophagus: (a) volume rendering and (b) phase
contrast micro tomography in the axial, longitudinal and coronal slices (from
left to right) and same slices for (c) conventional attenuation-contrast micro
tomography. The numerical labels correspond to five different tissue layers
that were identified in the phase contrast channel: epithelium (1), lamina
propria (2), submucosa (3), inner circular muscular layer (4) and longitudinal
outer muscular layer (5). These layers are hardly if at all visible in the
images of panel (c). For the rat heart: (d) volume rendering of the contrast
micro tomography (e) and its axial, longitudinal and coronal slices (from left
to right). The following structures were identified: the left pulmonary vein
(LPV), the aorta, the aortic valve (AV), the aortic arch (AA), the mitral
valve (MV), the left ventricle (LV), the right ventricle (RV), the papillary
muscles (PM), the right atrium (RA), and the caudal vena cava (CVC). We note
the presence of a blood clot (BC) within the RA. The circumferentially
oriented fibres in the mid-miocardium and the longitudinally oriented fibres
closer to the outer epicardium wall can also be observed, as highlighted by
the two zoom-in insets in panel (e).

For the piglet oesophagus (Figure 5a-c), we identified five different tissue
layers in the phase reconstructions as labelled in the images: epithelium (1),
lamina propria (2), submucosa (3), inner circular muscular layer (4) and
longitudinal outer muscular layer (5). These tissue layers are poorly, if at
all, visible in the attenuation-contrast reconstructions. For the rat heart
(Figure 5d and 5e), lab-XPC$\mu$T allowed to visualise the orientation of
micro-fibres in the muscle bundles and to identify the following structures:
the left pulmonary vein (LPV), the aorta, the aortic valve (AV), the aortic
arch (AA), the mitral valve (MV), the left ventricle (LV), the right ventricle
(RV), the papillary muscles (PM), the right atrium (RA), and the caudal vena
cava (CVC). We note the presence of an iron-rich blood clot (BC) appearing
brighter within the RA. The insets in the axial slice show the
circumferentially oriented fibres in the mid-miocardium (line-like structures)
and the longitudinally oriented fibres closer to the outer epicardium wall
(point-like structures). The tomography datasets were acquired in a step-and-
shoot fashion which carried overheads. Without optimisation the total
acquisition times were in the range of approximately 64 h, this mismatch is
only due to a need for optimising the sequences of positioning and detector
exposure. Work is currently underway to implement fly-scans for bringing the
data acquisition time much closer to the actual exposure times.

## Conclusion

We presented a compact and cost-effective system for performing X-ray phase-
contrast micro tomography in a laboratory setting. The system is based on a
single modulator, a low-power sealed X-ray tube and is compatible with a range
of readily available detectors. We report sub-micron angular sensitivity
starting at 2 seconds of exposure time at $50$ $\mu$m spatial resolution and
32 seconds at $16$ $\mu$m spatial resolution, and a maximum angular
sensitivity of 220 nrad in $64$ s at $50$ $\mu$m spatial resolution. These
results are compatible with the state of the art and were achieved with a
simple and compact set-up. The proposed approach is quantitative, as validated
through a phantom composed of known materials, and flexible in terms of
spatial resolution and sampling. By introducing the amplitude modulator,
sampling and resolution are driven by the structure of the illumination and
can be tuned to the requirement of the sample by adapting the data acquisition
strategy. The potential of the system for soft tissue imaging was demonstrated
on two biological specimens prepared without staining: five different tissue
layers were identified in a piglet oesophagus and the orientation of micro-
fibers within the myocardium of a rat heart was visualised. We believe that
the concurrent simplicity, sensitivity, robustness, compactness and efficiency
of our approach will be instrumental in making X-ray phase-contrast micro
tomography more accessible and available to a wider community.

## Author contributions

Conceptualization ME; Methodology ME CN-L AD GP SC AA; Software CN-L AD AA
CKH; Investigation CN-L AD ME; Resources ME GP SC AA SS MFMG AO; Writing CN-L
ME; Review CN-L AD SS MFMG GP AA SC DB CKH AO ME; Supervision AO DB CKH SC ME;
Funding acquisition AO DB CKH ME.

## Acknowledgements

This work was supported by the EPSRC EP/T005408/1; Wellcome Trust
221367/Z/20/Z; National Research Facility for Lab X-ray CT (NXCT) through
EPSRC grants EP/T02593X/1 and EP/V035932/1. AO is supported by the Royal
Academy of Engineering under their Chair in Emerging Technologies scheme
(CiET1819/2/78). CKH is supported by the Royal Academy of Engineering under
the Research Fellowship Scheme. MFMG is supported by an H2020 Marie Skłowdoska
Curie Action Fellowship (AmnioticID, grant agreement: 843265). SS is a UKRI
EPSRC Doctoral Prize Fellow (EP/T517793/1).

## References

  * [1] P. J. Withers, C. Bouman, S. Carmignato, V. Cnudde, D. Grimaldi, C. K. Hagen, E. Maire, M. Manley, A. D. Plessis, and S. R. Stock, “X-ray computed tomography,” Nature Reviews Methods Primers, vol. 1, Feb. 2021.
  * [2] S. R. Stock, “Recent advances in x-ray microtomography applied to materials,” International Materials Reviews, vol. 53, pp. 129–181, May 2008.
  * [3] E. Maire and P. J. Withers, “Quantitative x-ray tomography,” International Materials Reviews, vol. 59, pp. 1–43, Dec. 2013.
  * [4] D. Paganin, Coherent X-Ray Optics.  Oxford University Press, Jan. 2006.
  * [5] T. J. Davis, D. Gao, T. E. Gureyev, A. W. Stevenson, and S. W. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard x-rays,” Nature, vol. 373, pp. 595–598, Feb. 1995.
  * [6] S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. W. Stevenson, “Phase-contrast imaging using polychromatic hard x-rays,” Nature, vol. 384, pp. 335–338, Nov. 1996.
  * [7] F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance x-ray sources,” Nature Physics, vol. 2, pp. 258–261, Mar. 2006.
  * [8] M. Bech, T. H. Jensen, R. Feidenhans, O. Bunk, C. David, and F. Pfeiffer, “Soft-tissue phase-contrast tomography with an x-ray tube source,” Physics in Medicine and Biology, vol. 54, pp. 2747–2753, Apr. 2009.
  * [9] J. Herzen, T. Donath, F. Pfeiffer, O. Bunk, C. Padeste, F. Beckmann, A. Schreyer, and C. David, “Quantitative phase-contrast tomography of a liquid phantom using a conventional x-ray tube source,” Optics Express, vol. 17, p. 10010, May 2009.
  * [10] A. Olivo and R. Speller, “A coded-aperture technique allowing x-ray phase contrast imaging with conventional sources,” Applied Physics Letters, vol. 91, p. 074106, Aug. 2007.
  * [11] C. K. Hagen, P. R. T. Munro, M. Endrizzi, P. C. Diemoz, and A. Olivo, “Low-dose phase contrast tomography with conventional x-ray sources,” Medical Physics, vol. 41, p. 070701, June 2014.
  * [12] H. H. Wen, E. E. Bennett, R. Kopace, A. F. Stein, and V. Pai, “Single-shot x-ray differential phase-contrast and diffraction imaging using two-dimensional transmission gratings,” Optics Letters, vol. 35, p. 1932, June 2010.
  * [13] O. Hemberg, M. Otendal, and H. M. Hertz, “Liquid-metal-jet anode electron-impact x-ray source,” Applied Physics Letters, vol. 83, pp. 1483–1485, Aug. 2003.
  * [14] H. Miao, A. Panna, A. A. Gomella, E. E. Bennett, S. Znati, L. Chen, and H. Wen, “A universal moiré effect and application in x-ray phase-contrast imaging,” Nature Physics, vol. 12, pp. 830–834, Apr. 2016\.
  * [15] I. Zanette, T. Zhou, A. Burvall, U. Lundström, D. Larsson, M. Zdora, P. Thibault, F. Pfeiffer, and H. Hertz, “Speckle-based x-ray phase-contrast and dark-field imaging with a laboratory source,” Physical Review Letters, vol. 112, June 2014.
  * [16] I. Zanette, M.-C. Zdora, T. Zhou, A. Burvall, D. H. Larsson, P. Thibault, H. M. Hertz, and F. Pfeiffer, “X-ray microtomography using correlation of near-field speckles for material characterization,” Proceedings of the National Academy of Sciences, vol. 112, pp. 12569–12573, Sept. 2015.
  * [17] T. Zhou, F. Yang, R. Kaufmann, and H. Wang, “Applications of laboratory-based phase-contrast imaging using speckle tracking technique towards high energy x-rays,” Journal of Imaging, vol. 4, p. 69, May 2018.
  * [18] H. Wang, Y. Kashyap, and K. Sawhney, “From synchrotron radiation to lab source: advanced speckle-based x-ray imaging using abrasive paper,” Scientific Reports, vol. 6, Feb. 2016.
  * [19] L. Quénot, H. Rougé-Labriet, S. Bohic, S. Berujon, and E. Brun, “Implicit tracking approach for x-ray phase-contrast imaging with a random mask and a conventional system,” Optica, vol. 8, p. 1412, Nov. 2021.
  * [20] S. W. Wilkins, Y. I. Nesterets, T. E. Gureyev, S. C. Mayo, A. Pogany, and A. W. Stevenson, “On the evolution and relative merits of hard x-ray phase-contrast imaging methods,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 372, p. 20130021, Mar. 2014.
  * [21] A. Bravin, P. Coan, and P. Suortti, “X-ray phase-contrast imaging: from pre-clinical applications towards clinics,” Physics in Medicine and Biology, vol. 58, pp. R1–R35, Dec. 2012.
  * [22] R. Hanke, T. Fuchs, M. Salamon, and S. Zabler, “X-ray microtomography for materials characterization,” in Materials Characterization Using Nondestructive Evaluation (NDE) Methods, pp. 45–79, Elsevier, 2016.
  * [23] A. Momose, “X-ray phase imaging reaching clinical uses,” Physica Medica, vol. 79, pp. 93–102, Nov. 2020.
  * [24] S. Tao, C. He, X. Hao, C. Kuang, and X. Liu, “Principles of different x-ray phase-contrast imaging: A review,” Applied Sciences, vol. 11, p. 2971, Mar. 2021.
  * [25] L. Birnbacher, E.-M. Braig, D. Pfeiffer, F. Pfeiffer, and J. Herzen, “Quantitative x-ray phase contrast computed tomography with grating interferometry,” European Journal of Nuclear Medicine and Molecular Imaging, vol. 48, pp. 4171–4188, Apr. 2021.
  * [26] B. Müller, S. R. Stock, G. Wang, and J. G. Brankov, “Special Section Guest Editorial: Hard X-Ray Tomography with Micrometer Resolution,” Journal of Medical Imaging, vol. 9, no. 3, p. 031501, 2022.
  * [27] A. Velroyen, A. Yaroshenko, D. Hahn, A. Fehringer, A. Tapfer, M. Müller, P. Noël, B. Pauwels, A. Sasov, A. Yildirim, O. Eickelberg, K. Hellbach, S. Auweter, F. Meinel, M. Reiser, M. Bech, and F. Pfeiffer, “Grating-based x-ray dark-field computed tomography of living mice,” EBioMedicine, vol. 2, pp. 1500–1506, Oct. 2015.
  * [28] R. P. Murrie, F. Werdiger, M. Donnelley, Y. wei Lin, R. P. Carnibella, C. R. Samarage, I. Pinar, M. Preissner, J. Wang, J. Li, K. S. Morgan, D. W. Parsons, S. Dubsky, and A. Fouras, “Real-time in vivo imaging of regional lung function in a mouse model of cystic fibrosis on a laboratory x-ray source,” Scientific Reports, vol. 10, Jan. 2020.
  * [29] L. Massimi, T. Suaris, C. K. Hagen, M. Endrizzi, P. R. T. Munro, G. Havariyoun, P. M. S. Hawker, B. Smit, A. Astolfo, O. J. Larkin, R. M. Waltham, Z. Shah, S. W. Duffy, R. L. Nelan, A. Peel, J. L. Jones, I. G. Haig, D. Bate, and A. Olivo, “Detection of involved margins in breast specimens with x-ray phase-contrast computed tomography,” Scientific Reports, vol. 11, p. 3663, Feb 2021.
  * [30] M. Reichardt, M. Töpperwien, A. Khan, F. Alves, and T. Salditt, “Fiber orientation in a whole mouse heart reconstructed by laboratory phase-contrast micro-CT,” Journal of Medical Imaging, vol. 7, p. 1, Mar. 2020.
  * [31] S. Savvidis, M. F. Gerli, M. Pellegrini, L. Massimi, C. K. Hagen, M. Endrizzi, A. Atzeni, O. K. Ogunbiyi, M. Turmaine, E. S. Smith, C. Fagiani, G. Selmin, L. Urbani, N. Durkin, S. Shibuya, P. D. Coppi, and A. Olivo, “Monitoring tissue engineered constructs and protocols with laboratory-based x-ray phase contrast tomography,” Acta Biomaterialia, vol. 141, pp. 290–299, Mar. 2022\.
  * [32] M. Töpperwien, M. Krenkel, D. Vincenz, F. Stöber, A. M. Oelschlegel, J. Goldschmidt, and T. Salditt, “Three-dimensional mouse brain cytoarchitecture revealed by laboratory-based x-ray phase-contrast tomography,” Scientific Reports, vol. 7, Feb. 2017.
  * [33] M. Willner, M. Viermetz, M. Marschner, K. Scherer, C. Braun, A. Fingerle, P. Noël, E. Rummeny, F. Pfeiffer, and J. Herzen, “Quantitative three-dimensional imaging of lipid, protein, and water contents via x-ray phase-contrast tomography,” PLOS ONE, vol. 11, p. e0151889, Mar. 2016\.
  * [34] T. Weber, G. Pelzer, J. Rieger, A. Ritter, and G. Anton, “Report of improved performance in talbot-lau phase-contrast computed tomography,” Medical Physics, vol. 42, pp. 2892–2896, May 2015.
  * [35] X. Ji, R. Zhang, K. Li, and G.-H. Chen, “Dual energy differential phase contrast CT (DE-DPC-CT) imaging,” IEEE Transactions on Medical Imaging, vol. 39, pp. 3278–3289, Nov. 2020.
  * [36] D. Shoukroun, L. Massimi, F. Iacoviello, M. Endrizzi, D. Bate, A. Olivo, and P. Fromme, “Enhanced composite plate impact damage detection and characterisation using x-ray refraction and scattering contrast combined with ultrasonic imaging,” Composites Part B: Engineering, vol. 181, p. 107579, Jan. 2020.
  * [37] F. A. Vittoria, G. K. N. Kallon, D. Basta, P. C. Diemoz, I. K. Robinson, A. Olivo, and M. Endrizzi, “Beam tracking approach for single–shot retrieval of absorption, refraction, and dark–field signals with laboratory x–ray sources,” Applied Physics Letters, vol. 106, no. 22, p. 224102, 2015.
  * [38] E. Dreier, A. Bergamaschi, G. K. Kallon, R. Brönnimann, U. L. Olsen, A. Olivo, and M. Endrizzi, “Tracking based, high-resolution single-shot multimodal x-ray imaging in the laboratory enabled by the sub-pixel resolution capabilities of the mönch detector,” Applied Physics Letters, vol. 117, no. 26, p. 264101, 2020.
  * [39] M. D. de Jonge, B. Hornberger, C. Holzner, D. Legnini, D. Paterson, I. McNulty, C. Jacobsen, and S. Vogt, “Quantitative phase imaging with a scanning transmission x-ray microscope,” Phys. Rev. Lett., vol. 100, p. 163902, Apr 2008.
  * [40] P. C. Diemoz, F. A. Vittoria, and A. Olivo, “Spatial resolution of edge illumination x-ray phase-contrast imaging,” Opt. Express, vol. 22, pp. 15514–15529, Jun 2014.
  * [41] C. K. Hagen, F. A. Vittoria, M. Endrizzi, and A. Olivo, “Theoretical framework for spatial resolution in edge-illumination x-ray tomography,” Phys. Rev. Applied, vol. 10, p. 054050, Nov 2018.
  * [42] J. R. Fienup, M. Guizar-Sicairos, and S. T. Thurman, “Efficient subpixel image registration algorithms,” Optics Letters, Vol. 33, Issue 2, pp. 156-158, vol. 33, pp. 156–158, jan 2008.
  * [43] C. Kottler, C. David, F. Pfeiffer, and O. Bunk, “A two-directional approach for grating based differential phase contrast imaging using hard x-rays,” Opt. Express, vol. 15, pp. 1175–1181, Feb 2007.
  * [44] P. Modregger, B. R. Pinzer, T. Thüring, S. Rutishauser, C. David, and M. Stampanoni, “Sensitivity of x-ray grating interferometry,” Opt. Express, vol. 19, pp. 18324–18338, Sep 2011.
  * [45] W. van Aarle, W. J. Palenstijn, J. De Beenhouwer, T. Altantzis, S. Bals, K. J. Batenburg, and J. Sijbers, “The astra toolbox: A platform for advanced algorithm development in electron tomography,” Ultramicroscopy, vol. 157, pp. 35–47, 2015.
  * [46] P. C. Diemoz, C. K. Hagen, M. Endrizzi, and A. Olivo, “Sensitivity of laboratory based implementations of edge illumination x-ray phase-contrast imaging,” Applied Physics Letters, vol. 103, no. 24, p. 244104, 2013.
  * [47] G. Havariyoun, F. A. Vittoria, C. K. Hagen, D. Basta, G. K. Kallon, M. Endrizzi, L. Massimi, P. Munro, S. Hawker, B. Smit, A. Astolfo, O. J. Larkin, R. M. Waltham, Z. Shah, S. W. Duffy, R. L. Nelan, A. Peel, T. Suaris, J. L. Jones, I. G. Haig, D. Bate, and A. Olivo, “A compact system for intraoperative specimen imaging based on edge illumination x-ray phase contrast,” Physics in Medicine & Biology, vol. 64, p. 235005, nov 2019\.
  * [48] V. Revol, C. Kottler, R. Kaufmann, U. Straumann, and C. Urban, “Noise analysis of grating-based x-ray differential phase contrast imaging,” Review of Scientific Instruments, vol. 81, no. 7, p. 073709, 2010.
  * [49] T. Thüring, S. Hämmerle, S. Weiss, J. Nüesch, J. Meiser, J. Mohr, C. David, and M. Stampanoni, “Compact hard X-ray grating interferometry for table top phase contrast micro CT,” in Medical Imaging 2013: Physics of Medical Imaging (R. M. Nishikawa, B. R. Whiting, and C. Hoeschen, eds.), vol. 8668, pp. 293 – 299, International Society for Optics and Photonics, SPIE, 2013.
  * [50] L. Birnbacher, M. Willner, A. Velroyen, M. Marschner, A. Hipp, J. Meiser, F. Koch, T. Schröter, D. Kunka, J. Mohr, F. Pfeiffer, and J. Herzen, “Experimental Realisation of High-sensitivity Laboratory X-ray Grating-based Phase-contrast Computed Tomography,” Scientific Reports 2016 6:1, vol. 6, pp. 1–8, apr 2016.
  * [51] J. Vila-Comamala, L. Romano, K. Jefimovs, H. Dejea, A. Bonnin, A. C. Cook, I. Planinc, M. Cikes, Z. Wang, and M. Stampanoni, “High sensitivity x-ray phase contrast imaging by laboratory grating-based interferometry at high talbot order geometry,” Optics Express, vol. 29, p. 2049, Jan. 2021.
  * [52] “X-ray interactions: Photoabsorption, scattering, transmission, and reflection at e = 50-30,000 ev, z = 1-92,” Atomic Data and Nuclear Data Tables, vol. 54, no. 2, pp. 181–342, 1993.
  * [53] P. R. Munro and A. Olivo, “X-ray phase-contrast imaging with polychromatic sources and the concept of effective energy,” Physical Review A, vol. 87, p. 053838, may 2013.