# Real-time super-resolution mapping of locally anisotropic grain orientations
for ultrasonic non-destructive evaluation of crystalline material

J. Singh1,∗, K. M. M. Tant1, A. Curtis2, A. J. Mulholland3

1Department of Mathematics and Statistics, University of Strathclyde
2School of Geosciences, University of Edinburgh
3Department of Engineering Mathematics, University of Bristol
∗Corresponding author<EMAIL_ADDRESS>
###### Abstract

Estimating the spatially varying microstructures of heterogeneous and locally
anisotropic media non-destructively is necessary for the accurate detection of
flaws and reliable monitoring of manufacturing processes. Conventional
algorithms used for solving this inverse problem come with significant
computational cost, particularly in the case of high dimensional non-linear
tomographic problems. In this paper, we propose a framework which uses deep
neural networks (DNNs) with full aperture, pitch-catch and pulse-echo
transducer configurations to reconstruct material maps of crystallographic
orientation. We also present the first ever application of generative
adversarial networks (GANs) to achieve super resolution of ultrasonic
tomographic images, providing a factor-four increase in image resolution and
up to a 50% increase in structural similarity. The importance of including
appropriate prior knowledge in the GAN training dataset to increase inversion
accuracy is highlighted; known information about the material’s structure
should be present in the training data. We show that after a computationally
expensive training process, the DNNs and GANs can be used in less that one
second (0.9 seconds on a standard desktop computer) to provide a high
resolution map of the material’s grain orientations.

##  1 Introduction

Ultrasonic non-destructive evaluation (NDE) is widely used across a number of
industries including aerospace, nuclear, and oil and gas. The technique
involves the generation, transmission and reception of high-frequency
mechanical waves through a component [10]. An image of the component’s
interior is then generated via post processing of this data to aid in the
detection of any internal defects [31]. Conventional ultrasonic imaging
algorithms within NDE typically assume that the material that is being
inspected is isotropic and homogeneous. However, metals can develop locally
anisotropic and heterogeneous microstructures, particularly when they are
subjected to extreme thermal cycles, such as those present in welding and
additive manufacturing processes [15, 37, 45]. Conventional ultrasonic imaging
algorithms which assume homogeneity or isotropy can fail to focus the energy
correctly in the image domain in such cases and are therefore unreliable [33,
41, 50]. Algorithms which incorporate a priori information about a material’s
spatially varying properties significantly improve the accuracy of defect
characterisation [41].

In recent years, much effort has been expended on generating material property
maps non-destructively using tomographic inversion, where material properties
such as wave speed, or microstructural descriptors such as grain orientation,
are estimated from the scattered wave field data recorded at the surface of an
object. A wide range of advanced tomographic algorithms are used across
geophysics [1, 2, 11, 29, 43, 51, 53], bio-medicine [16] and NDE [13, 27, 41,
42]. A common approach is to use iterative methods to improve the fit of the
measured data to forward modelled data which depend on an estimate of the
material map. They sample potential material maps from some multi-dimensional
parameter space, solve a forward problem for each new material property map,
and update the estimated map to improve the data fit [27]. In the case of
probabilistic sampling frameworks (for example, those built around Markov
chain Monte Carlo methods [42, 52]), there is the added benefit of extracting
uncertainty information on the parameter estimates, facilitating valuable
uncertainty quantification studies. Although these algorithms have
demonstrated impressive results in reconstructing wave speed and grain
orientation maps, they are computationally demanding, often requiring the
storage of large sample sets and compute times of several hours to several
weeks. This poses a problem for the NDE community, where there is an
increasing demand for the monitoring of dynamical processes employed during
manufacturing, for example in welding and additive manufacturing processes
[23, 24], and so it is desirable to carry out inspection in real-time.

Machine learning shows strong potential to solve material characterisation
inverse problems rapidly [17]. Specifically, we focus on the use of deep
neural networks (DNNs), which can approximate any non-linear relationship
between two parameter spaces, given a sufficiently large set of training data
(pairs of dependent and corresponding independent parameters [8]). The
training of a DNN is computationally expensive. However, the training process
is only performed once prior to using a DNN, and a trained network can be used
effectively in real time without the need for high-performance computing.

Inversion methods based on DNNs have become increasingly popular for
tomographic imaging of isotropic material properties, particularly in
geophysics [5, 7, 12, 17, 32] and bio-medicine [4, 48]. However, DNNs have not
yet been implemented for tomographic reconstruction of anisotropic material
properties. Although various deep learning algorithms have been used to solve
inverse problems in NDE, for example, to predict material fatigue behaviour
[3], to augment ultrasonic data [44], and for ultrasonic crack
characterisation [35] and crack detection using image recognition [14, 21],
the use of DNNs for tomography has yet to be explored in this context.

In addition to DNNs, generative adversarial networks (GANs) have more recently
been applied to various computer vision tasks, including achieving super-
resolution with upscaling by up to a factor of four [30], colourisation [19],
segmentation and labelling [22]. This family of algorithms has strong
potential to improve image resolution, and has been used increasingly in
remote sensing [25] and X-ray tomography [49], however there has been no
application of GANs in NDE to produce ultrasonic tomographic images.

In this paper, we present the first DNN framework for rapid, non-linear two-
dimensional tomography of heterogeneous and locally anisotropic materials. The
datasets used for the tomographic inversion are the arrival times of
ultrasonic waves which have been transmitted and received by an array of
sensors on the exterior of the component. The examples shown are inspired by
the NDE of polycrystalline materials but the methodology should naturally
extend to other domains, for example imaging anisotropic fibrous tissue [18,
20] or the Earth’s subsurface [54]. We compare the network’s performance for a
range of transducer configurations, model textures and different types of
simulated ultrasonic testing data (i.e. we move beyond inverse crime
scenarios). A novel GAN-based method for post-processing ultrasound
tomographic images to achieve super-resolution with a four-fold upscaling
factor is presented, achieving up to 50% improvement using structural
similarity metrics. We define the term super-resolution in the context of
image processing, as reconstructing images below the original lengthscale.
This is different to an alternative definition often used in physical
acoustics, which is to image below the wavelength in the data.

##  2 Method

We employ model-driven deep learning, where a large dataset of simulated
material maps and corresponding travel time measurements are used to train a
DNN and hence solve the tomographic inverse problem. The forward modelling
problem can be denoted as

$f(\mathbf{m},\mathbf{s})=\mathbf{T_{m}},$ (1)

where $f$ is a forward mechanical wave modeling operator, $\mathbf{m}$ is a
material model, $\mathbf{s}$ contains the locations of the elements in the
ultrasonic transducer array and $\mathbf{T_{m}}$ is the time-of-flight (ToF)
matrix between every pair of array elements. Within each database used for
network training, the transducer configuration $\mathbf{s}$ is fixed and
therefore $\mathbf{s}$ is omitted in the notation for the ToF matrix
$\mathbf{T_{m}}$. We use deep learning to obtain (or learn) an approximation
of $f^{-1}$, which maps the measured data $\mathbf{T_{m}}$ to a material map
$\mathbf{m}$ (i.e., $DNN\approx f^{-1}$). In this study, the training data
consists of two-dimensional material models with spatially varying crystal
orientations $\theta(x,y)$ and the travel time matrix $\mathbf{T_{m}}$
corresponing to each one. To generate models in such a way that the
distribution of orientations are randomly assigned but still exhibit some
structural correlation, an initial random Voronoi tessellation [39] with 30
seeds (a set of two-dimensional Cartesian coordinates lying within the domain
of interest) is computed and an orientation $\theta$ between $0^{\circ}$ and
$45^{\circ}$ is randomly assigned to each of the 30 resulting Voronoi regions
or cells (Fig.1a). We consider only in-plane crystal rotation, and therefore
the orientation $\theta$ relates to the orientation of a slowness curve in
each cell. This slowness curve plots the reciprocal of velocity in the crystal
over a range of incident wave directions [42]. The material models used in the
training data $\\{\mathbf{m}_{16},\mathbf{T_{m_{16}}}\\}$ are generated by
discretising the Voronoi tessellation into a regularly spaced $16\times 16$
grid and smoothing with a Gaussian kernel (Fig.1b, the subscript 16 denotes
the model resolution). The smoothing simplifies the inverse problem so that
only smooth models are inverted for. To demonstrate that this machine learning
approach can be generalised for any locally anisotropic media, the
longitudinal group slowness curve is obtained for an arbitrary anisotropic
material with a cubic stiffness tensor, where $c_{11}=256.45$ GPa,
$c_{12}=133.5$ GPa and $c_{44}=c_{12}$ and density $\rho=7874$ kgm-3. Three
configurations of ultrasonic transducer array locations $\mathbf{s}$ are
considered: a full aperture coverage of 16 elements (4 on each face as shown
in Fig. 1d), a two-sided aperture pitch-catch configuration with 16
transmitting elements at the top of the model, with the time-of-flights
measured at 16 receiving elements at the bottom of the model (Fig. 1e), and a
one-sided aperture pulse-echo configuration, where 16 elements are positioned
along the top face and the travel times of waves reflecting off the bottom
face and returning to the transducers on the top face are measured (Fig. 1f).
The measured data is the time-of-flight (ToF) of each propagating wave between
each pair of array elements, represented in a ToF matrix $\mathbf{T_{m_{16}}}$
shown in Figure 1(c).

Figure 1: Illustration of the procedure for generating training data: (a)
Randomly generated Voronoi tessellation with 30 seeds and random grain
orientations ranging between 0∘ and 45∘ (other angles are included due to
symmetry of the slowness curve). (b) Material map $\mathbf{m}_{16}$ generated
by discretising the Voronoi image (a) on a 16x16 grid and then smoothing with
a Gaussian kernel. (c) An example travel time matrix $T_{\mathbf{m}_{16}}$
populated by measurements of the time-of-flight between every pair of
transducer elements. In this paper, we position transducer elements in three
configurations: (d) full aperture coverage by 16 elements (4 on each face),
(e) two-sided aperture coverage with 16 transmitting elements at the top of
the model with the pitch-catch time-of-flight being measured at 16 receiving
elements at the bottom of the model, (f) one-sided aperture coverage, where 16
transducers are positions along the top face and the pulse-echo arrival time
of the wave reflecting from the bottom face is measured.

###  2.1 Forward Model Approaches

Following model parameterisation, an efficient forward model is required for
computing the time-of-flight matrix $\mathbf{T_{m}}$ (Fig.1d) corresponding to
a grain orientation model $\mathbf{m}_{16}$ for each source-receiver pair. We
take two approaches: a semi-analytic model using an anisotropic multi-stencil
fast marching method (AMSFMM) algorithm from [42], denoted as
$f_{\mbox{FMM}}$, and a finite element analysis (FEA) method, denoted as
$f_{\mbox{FEA}}$. The AMSFMM incorporates the effects of ray bending due to
variations in locally anisotropic grain orientations, and models the travel-
time field by solving the Eikonal equation using an upwind finite difference
scheme [36, 40, 42]. This allows the calculation of the shortest travel time
between transmitter and receiver locations and the matrix
$\mathbf{T_{m}}^{\mbox{FMM}}$ can be constructed (that is
$\mathbf{T_{m}}^{\mbox{FMM}}=f_{\mbox{FMM}}(\mathbf{m}_{16})$). As wave
reflections are not incorporated into the AMSFMM, a different approach is
required for the pulse-echo transducer array configuration. In this case, the
time of flight between the transmitter and receiver is calculated by the
summation of the time of flight between the transmitter to all points along
the back-wall and between the receiver and all points along the bottom face.
The output of this summation is an array of travel times corresponding to all
the reflection points along the bottom face, and the minimum value is taken to
be the time of flight for the pulse-echo transducer array configuration. The
FEA method incorporates more of the underlying physics in the model compared
to AMSFMM, as it models full wave propagation including multiple scattering
and diffraction. To measure the ToF of the received waves, an amplitude
threshold is selected and the time for the recorded wave amplitude to reach
this threshold is used as an element of the travel time matrix
$\mathbf{T_{m}}^{\mbox{FEA}}$ (that is
$\mathbf{T_{m}}^{\mbox{FEA}}=f_{\mbox{FEA}}(\mathbf{m}_{16})$). The FEA method
is significantly more computationally expensive than the AMSFMM. As a large
number of data-model pairs are required to train a deep neural network, the
more efficient AMSFMM method is used to generate travel time matrices
$\mathbf{T_{m}}^{\mbox{FMM}}$ for the training data. The more physically
realistic FEA generated data is then used to generate data to test the trained
networks’ performance (see FEA set-up in the Appendix). A total of 7500 models
are generated and the corresponding travel time matrices
$\mathbf{T_{m}}^{\mbox{FMM}}$ are computed using AMSFMM for the training data
set.

###  2.2 Deep neural network for orientation mapping

Deep neural networks (DNNs) are mathematical mappings that emulate the
relationship between two parameter spaces [17]. Here, we seek a map between
the grain orientation models $\mathbf{m}_{16}$ and the corresponding time of
flight data $\mathbf{T_{m}}$ (that is
$DNN(T_{16}^{\mbox{FMM}})=\mathbf{m}_{16}^{\mbox{pred}}$, where the pred
superscript denotes the DNN prediction). For each of the transducer
configurations $\mathbf{s}$, a different number of travel times are used as
input to the neural network. For a full aperture configuration (Fig. 1d), we
have $n$ source-receivers per side of our rectangular domain, and so there are
$6n^{2}$ unique travel times (accounting for source-receiver reciprocity and
excluding those between elements which lie on the same side). When $n=4$, a
set of 96 travel times is taken from each ToF matrix $\mathbf{T_{m}}$. For the
pitch-catch configuration configuration (Fig. 1e), all source-receiver paths
are unique, therefore the full ToF matrix is used and with $n=16$ there are
256 inputs to the neural network. Finally, for the pulse-echo configuration
(Fig. 1f), when $n$ is even, there are $n^{2}/2+n/2$ unique travel times
(accounting for source-receiver reciprocity), so when $n=16$, a total of 136
travel times are selected for the ToF matrix. For network training, both the
input travel times and the output orientations are scaled to have zero mean
and unit variance.

We configure three DNNs (corresponding to three transducer configurations),
each with five fully connected layers (illustrated in Fig. 2), using sigmoid
activation functions. The final output layer contains a single node
corresponding to the orientation of a single pixel in the imaging domain.
Therefore, following the approach of [17], a separate network is trained for
each pixel, so for a $16\times 16$ resolution image a total of 256 networks
are trained. The networks are trained using the Adam optimisation algorithm
[28]. A description of network hyper-parameters is provided in Appendix 6.2.
These hyper-parameters are selected using a stochastic optimisation library
[6] for each network architecture corresponding to different transducer
configurations. We use a mean-squared-error (MSE) loss function, given by:

$MSE=\frac{\sum^{N}_{i=1}(\mathbf{m}_{16}^{\mbox{true}}(i)-\mathbf{m}_{16}^{\mbox{pred}}(i))^{2}}{N},$
(2)

where $\mathbf{m}_{16}^{\mbox{true}}$ and $\mathbf{m}_{16}^{\mbox{pred}}$ are
the true and predicted grain orientation models, $i$ denotes the pixel index
and $N$ is the total number of pixels (for models $\mathbf{m}_{16}$, $N=256$).
A validation data set is created using 20% of the training data. To avoid
over-fitting the network to the training data, the cost function is
periodically evaluated over the validation data set, and we implement an early
stopping algorithm so that training stops once the validation loss stops
decreasing (with a patience of 10 iterations). The time to train 256 separate
networks sequentially using Google Colab [9] free GPUs is approximately 40
minutes, although training could be parallelised to reduce this time if
required. Once trained, the compute time of $\mathbf{m}_{16}^{\mbox{pred}}$ is
approximately 0.15 seconds per model inversion.

Figure 2: Schematic of the tomography algorithm using deep neural networks
(DNNs). First, $N$ travel times are selected from the travel time matrix
(where $N=96$, $N=256$ and $N=136$ for the full aperture, pitch-catch and
pulse-echo configurations respectively). Travel times are then used as input
into the DNN. The DNN consists of 3 hidden layers ($L1$, $L2$ and $L3$) and a
final output layer. The nodes are illustrated as circles and the number of
nodes in each layer is denoted in the bottom circles ($L1_{N}$, $L2_{N}$ and
$L3_{N}$). Each output corresponds to the crystal orientation of a single
pixel in the material map $\mathbf{m_{16}}$ so 256 separate networks are
trained in order to predict all pixel orientations.

###  2.3 Generative adversarial networks for super resolution

Figure 3: (a) Schematic of the Generative Adversarial Network (GAN) post-
processing algorithm for achieving super-resolution material maps. The 16x16
output of the DNN tomography algorithm $\mathbf{m}_{16}$ is input to the
generator network, which outputs a 64x64 image $\mathbf{m}_{64}^{G}$. The
known 64x64 image $\mathbf{m}_{64}^{T}$, which was used to generate the DNN
input data, as well as the output of the generator $\mathbf{m}_{64}^{G}$ are
input into the discriminator, which outputs a prediction of which image is
generated and which belongs to the training dataset. Three different model
data sets are used for GAN training: (b) a layer model with up to 5 layers
where positions of the interfaces are random, (c) a random Voronoi
tessellation with 6 seeds, and (c) a random Voronoi tessellation with 30
seeds.

Conditional GANs learn a mapping between two images [22] and so can be used
for post-processing of the DNN tomography output
($\mathbf{m}_{16}^{\mbox{pred}}$) to increase resolution and accuracy. The GAN
architecture, as illustrated in Figure 3(a), consists of two separate
trainable networks: a generator ($\mbox{GAN}_{G}$) and a discriminator
($\mbox{GAN}_{D}$).

Training a GAN for post-processing the output of the DNN tomography method
($\mathbf{m}_{16}^{\mbox{pred}}$) to achieve an increase in image resolution
(super-resolution) requires an additional training data set, where travel time
data is generated using higher resolution models (64x64). The GAN framework
assumes some prior knowledge of the structure of the material which is
incorporated into the GAN training data. For example, in layered structures
such as Carbon fiber reinforced polymers (CFRPs), the training data should
include models with locally anisotropic layers; or alternatively models
exhibiting crystalline grain structures should be used to train GAN’s for
cases such as welds, and knowledge on the average grain size could feed into
the complexity of the models included in the training data. We use three
separate training data sets of increasing complexity. The first high
resolution model $\mathbf{m}_{64}^{\mbox{true}}$ consists of up to 5
horizontal layers where the orientation and thickness of each layer is
randomly assigned (Fig. 3b). The second and third are generated by
discretising random Voronoi tessellations with 6 and 30 seed locations into a
64x64 grid, as shown in Figures 3(c) and (d), respectively. The travel time
matrices $\mathbf{T_{m}}_{64}^{FMM}$ are calculated using the AMSFMM algorithm
for 2000 models for each of the three data sets, which are input into the DNN
tomography algorithm described in the previous section, which outputs a 16x16
predicted model $\mathbf{m}_{16}^{\mbox{pred}}$. The generator is configured
to take the low resolution $\mathbf{m}_{16}^{\mbox{pred}}$ image as input and
to output a high resolution 64x64 image $\mathbf{m}_{64}^{G}$ (i.e.,
${\mbox{GAN}}_{G}(\mathbf{m}_{16}^{\mbox{pred}})=\mathbf{m}_{64}^{G}$). Here,
the generator is a modified U-net [38] based on fully convolutional layers
(see the Appendix for network architecture). The discriminator takes the
output of the generator $\mathbf{m}_{64}^{G}$, as well as the known 64x64 high
resolution image ($\mathbf{m}_{64}^{\mbox{true}}$) that was used to generate
the ToF data, and predicts which image is generated (fake) and which is part
of the training data (real). The accuracy of the discriminator prediction can
then be established. These competing networks are then trained against each
other; in each iteration of training, the accuracy of the discriminator is fed
into the loss function of the generator network. The generator seeks to create
images $\mathbf{m}_{64}^{G}$ that decrease the discriminator accuracy meaning
that $\mathbf{m}_{64}^{G}$ cannot be discriminated from the reference training
data $\mathbf{m}_{64}^{\mbox{true}}$. Following the training process, the
generator can be used to map from 16x16 images to 64x64 resolution images.

##  3 Results

###  3.1 DNN Results

Figure 4: (a) True material orientation maps $\mathbf{m}^{\mbox{true}}_{16}$
from the test data set and corresponding predicted models
$\mathbf{m}^{\mbox{pred}}_{16}$ using the DNN tomography algorithm for (b)
full aperture, (c) pitch-catch, and (d) pulse-echo transducer array
configurations. The associated mean absolute errors (MAE) and the structural
similarity index metrics (SSIM) between the predicted and true maps are
labelled above each reconstructed material map.

Following the training of the fully connected DNN, we predict material maps
$\mathbf{m}_{16}^{\mbox{pred}}$ using the three transducer array
configurations shown in Figure 1 following

$\mathbf{m}_{16}^{\mbox{pred}}=DNN(\mathbf{T_{m}}^{\mbox{FMM}}),$ (3)

where $\mathbf{T_{m}}^{\mbox{FMM}}$ is test data which has not been used in
the network training process. The test data are generated following the same
protocol as for the training data, using smoothed Voronoi models
$\mathbf{m}_{16}$ and the AMSFMM algorithm to generate a total of 200 test
models and data. Comparisons of the true models
$\mathbf{m}^{\mbox{true}}_{16}$ with the predicted models
$\mathbf{m}^{\mbox{pred}}_{16}$ using the DNN and with full aperture, pitch-
catch and pulse-echo transducer array configurations are shown in Figure 4. We
use two metrics for comparing predicted models with the true models: the mean
absolute error (MAE), which is a scalar value (MAE$\geq$ 0, where MAE=0
describes a perfect prediction), and the structural similarity index measure
(SSIM) [46] (-1$\leq$SSIM$\leq$1, where SSIM=1 describes a perfect
prediction). The SSIM incorporates the similarity of three independent
parameters: image luminescence, contrast and structure (see the Appendix).
These values are calculated with orientations that are scaled to have zero
mean and unit variance. Note that lower values of MAE indicate higher
similarity between the true and predicted models, whereas higher values of
SSIM indicate higher image similarity.

In all cases, the predicted material property maps resemble the true
orientation maps, predicting the magnitude and location of areas with similar
orientations. The DNN predictions with a full aperture experimental
configuration (Fig. 4b) perform the best (lower MAEs and higher SSIMs) and
predictions made using the pulse-echo configuration perform the worst (higher
MAEs and lower SSIMs). The histograms of MAE and SSIM values for the 200 test
models are shown in Figures 5(a) and (b). The distributions of the pixel mean
absolute error (averaged for each pixel across the 200 models) are shown for
each transducer array configuration in Figures 5(c-e), showing that
reconstruction accuracy generally decreases (increasing pixel MAE) in the
central region of the domain and with distance from the transmitting element
transducer array.

Figure 5: Histogram of (a) mean absolute error (MAE) and (b) structural
similarity index metric (SSIM) for DNN prediction on 200 test models using
different transducer array configurations. Lower values of MAE and higher
values of SSIM suggest a good tomographic reconstruction. The lower panels
show the mean absolute error (MAE) for each pixel averaged across the 200 test
models for (c) full aperture, (d) pitch-catch and (e) pulse-echo transducer
array configurations. Figure 6: (a) True material maps
($\mathbf{m}_{16}^{\mbox{true}}$). The predicted models using the AMSFMM
trained DNN tomography algorithm with (b)
$\mathbf{m_{16}}^{\mbox{pred}}=DNN(\mathbf{T_{m}}^{\mbox{FMM}}$) and (c)
$\mathbf{m_{16}}^{\mbox{pred}}=DNN(\mathbf{T_{m}}^{\mbox{FEA}}$).

So far, the same mathematical model has been used for both the training data
and the test data (a so called inverse crime [47]), and this is not a
sufficient challenge of the methodology [26]. We therefore now use a different
mathematical model to test the trained DNN. One further additional challenge
is to generate material maps using a different method from that used in the
training data; so not originating from Voronoi diagrams. The material maps in
Figure 6(a) show a range of structures including a homogeneous model, a
checkerboard structure, a layered structure and a single circular anomaly; all
of which are significantly dissimilar from the textures and structures found
within the training data. The FEA method is used to generate ToF data
$\mathbf{T_{m}}^{FEA}$ using a full aperture transducer array configuration,
which is then input into the DNN to predict the grain orientation map
$\mathbf{m}_{16}^{\mbox{pred}}$

$\mathbf{m}_{16}^{\mbox{pred}}=DNN(\mathbf{T_{m}}^{\mbox{FEA}}).$ (4)

The predicted material maps $\mathbf{m}_{16}^{\mbox{pred}}$ shown in Figures
6(b) and (c) show similar results using $\mathbf{T_{m}}^{\mbox{FEA}}$ and
$\mathbf{T_{m}}^{\mbox{FEA}}$ time of flight data. In the cases of the
homogeneous model and the single circular anomaly the results using
$\mathbf{T_{m}}^{\mbox{FEA}}$ are slightly improved (lower MAE). The
similarity of results between the two data types indicates that the DNN is
robust to changes in different data simulation methods and to the noise in the
FEA dataset associated with the identification of travel times. The presence
of this additional noise does not appear to have a significant effect on the
changes in measured travel time due to anisotropy, therefore the inversion
remains accurate. The accuracy of the predicted models is lower where the
material maps exhibit different textures to those used in the training data;
compare the MAE and SSIM values in Figure 6(c) with those in Figure 4(b). The
higher accuracy of the results in Figure 4(b) highlights that the texture of
the target application material for the DNN tomography algorithm should be
included as far as possible in the training dataset.

###  3.2 GAN Results

Three GANs are trained using the layered, 6-seed Voronoi and 30-seed Voronoi
models $\mathbf{m}_{64}^{\mbox{true}}$, and 200 additional models per GAN are
used for testing, of which 5 are shown in Figures 7(a), 8(a) and 9(a),
respectively. The AMSFMM method is used to compute travel time data
($\mathbf{T_{m}}^{\mbox{FMM}}$) using a full aperture transducer array
configuration, which are input into the trained DNN (as used for the
generation of DNN predictions in Fig. 4b). The DNN predicted outputs
$\mathbf{m}_{16}^{\mbox{pred}}$ are shown in Figures 7(b), 8(b) and 9(b) and
the GAN outputs $\mathbf{m}_{64}^{G}$ in Figures 7(c), 8(c) and 9(c) for the
layered, 6-seed Voronoi and 30-seed Voronoi models, respectively. In order for
image comparison with MAE and SSIM, the 16x16 resolution DNN outputs are
upscaled to 64x64 resolution using nearest neighbour interpolation. Histograms
of the changes in MAE ($\Delta MAE=MAE_{GAN}-MAE_{DNN}$) and SSIM ($\Delta
SSIM=SSIM_{GAN}-SSIM_{DNN}$) when using a GAN to post-process the DNN
tomography outputs are shown in Figure 10.

Figure 7: (a) True high resolution (64x64) grain orientation maps
$\mathbf{m}_{64}^{\mbox{true}}$ consisting of 5 horizontal layers, (b) 16x16
resolution DNN tomography output ($\mathbf{m}_{16}^{\mbox{pred}}$), and (c)
64x64 GAN output $\mathbf{m}_{64}^{G}$. For row (b), the MAE and SSIM are
calculated on an upscaled image to $64\times 64$ resolution using nearest
neighbour interpolation. Note the significant improvements in MAE and SSIM
using the GAN methodology and the clear improvements in reconstructing a
layered structure.  Figure 8: (a) True high resolution (64x64) grain
orientation maps $\mathbf{m}_{64}^{\mbox{true}}$ consisting of a Voronoi
Tesselation with 6 seeds, (b) 16x16 resolution DNN tomography output
($\mathbf{m}_{16}^{\mbox{pred}}$), and (c) 64x64 GAN output
$\mathbf{m}_{64}^{G}$. For row (b), the MAE and SSIM are calculated on an
upscaled image to $64\times 64$ resolution using nearest neighbour
interpolation. Note the significant improvements in MAE and SSIM using the GAN
methodology and the clear improvements in reconstructing piecewise constant
structures. Figure 9: (a) True high resolution (64x64) grain orientation maps
$\mathbf{m}_{64}^{\mbox{true}}$ consisting of a Voronoi Tesselation with 30
seeds, (b) 16x16 resolution DNN tomography output
($\mathbf{m}_{16}^{\mbox{pred}}$), and (c) 64x64 GAN output
$\mathbf{m}_{64}^{G}$. For row (b), the MAE and SSIM are calculated on an
upscaled image to $64\times 64$ resolution using nearest neighbour
interpolation. Note the marginal improvements in MAE and SSIM using the GAN
methodology. Figure 10: Histograms showing the change in (a) the mean absolute
error (MAE) and (b) the structural similarity index measure (SSIM) when using
a GAN to post-process 200 low resolution DNN tomography outputs ($\Delta
MAE=MAE_{GAN}-MAE_{DNN}$ and $\Delta SSIM=SSIM_{GAN}-SSIM_{DNN}$).

For the 5 layer models (Fig. 7), the GAN predictions are significantly more
accurate compared to DNN predictions, offering large improvements in MAE
(decrease up to $\Delta MAE=-0.85$) and SSIM (increase up to $\Delta
SSIM=0.5$). The GAN successfully learns to generate horizontal (layered)
structures, so very little horizontal variation exists in the GAN predictions.
The reconstructed grain orientation maps from the GAN exhibit discontinuous
grain boundaries and piecewise constant orientations for each layer, compared
to the smooth spatially varying DNN tomography outputs (Fig. 7b). The GAN also
performs well for the 6-seed Voronoi tessellation models (Fig. 8), where
reconstructed grain orientation maps from the GAN exhibit discontinuous,
piecewise constant orientations for each grain. The GAN improves MAE and SSIM
in all cases, however there is slight blurring across some grain boundaries.
The GAN results for the 30-seed Voronoi tessellation models (Fig. 9) exhibit
stronger blurring across grain boundaries. While the GAN prediction are
texturally more similar to the true models (piecewise constant and
discontinuous regions), the distributions of $\Delta MAE$ and $\Delta SSIM$ in
Figure 10 show the GAN offers only marginal improvements in reconstruction
accuracy, and in some cases the accuracy decreases when using the GAN ($\Delta
MAE>0$ and $\Delta SSIM<0$). The difference between the 6-seed and 30-seed
Voronoi models is in the model complexity due to smaller individual grains in
the 30-seed models. In these models, multiple grains can fit into a single
pixel of a low resolution DNN tomography image, resulting in a loss of spatial
information that the GAN cannot fully recover. These results show that a GAN
can be used for post-processing tomography results to improve reconstruction
accuracy and image resolution, particularly when prior information regarding
the spatial distribution of the material map is known (e.g., if the sample is
known to be layered, or similarly well-structured) and the spatial
distribution is simple.

##  4 Discussion

The framework presented includes several stages: (1) the generation of
training data using the AMSFMM method, (2) training of the DNN, and (3)
training of the GAN. However, each of these stages only need be performed
once. Thereafter, the DNN and GAN can be used in effectively real time ($<1$
second). Here, the time for generating 7500 ToF matrices
$T_{\mathbf{m}}^{\mbox{FMM}}$ was approximately 1 hour, for training the DNN
was approximately 40 minutes (until convergence), and for training the GAN was
approximately 8 hours (using Google Colab GPUs [9]). It is clear that when
repeated material map reconstructions are desired, as is the case for NDE
monitoring purposes, the deep learning framework excels in its ability to
provide real time results.

The benefits of real time inversions comes at the expense of a few limitations
that are yet to be overcome in the current work. Firstly, the DNN is trained
with a constant transducer configuration, so a trained DNN cannot be generally
extended to changes in relative transducer locations. This is not a problem
for many applications in NDE, as the transducer arrays are rigid and fixed,
and the test sample geometries do not change through time. However, limited
network flexibility may be problematic in cases where the configuration
changes, such as in-process monitoring of additive manufacturing: during the
building process, the shape of the sample changes therefore the distribution
of transducer elements also changes. One solution is to train many DNNs for
all the possible transducer configurations throughout the building process,
however this would require a significantly expensive training process. Another
solution, proposed in [17], is to train more flexible networks that account
for missing data by augmenting the training data set with additional input
samples taken from additional transducer locations. Travel times in the ToF
matrix can be set to zero to indicate a transducer is not used for a
particular transducer configuration, and then the trained network can invert
using multiple configurations.

The GAN is also limited in its applicability. This is highlighted when a
trained GAN is used to invert for textures that are dissimilar to those found
in the training data. This can be seen in Figure 11, where the GAN trained on
the 30-seed Voronoi models is applied to the DNN prediction of the
checkerboard, layered and circular inclusion material models as well as a
30-seed Voronoi model for reference. There are significant decreases in SSIM
and increases in MAE when using the GAN on the models with dissimilar textures
to the Voronoi models. This highlights the importance of the training data
used in network training and suggests the GAN should only be used if prior
knowledge of the material is known and the expected textures are present in
the training data. Training a GAN with a much broader training data set, for
example including all of the layered, 6-seed and 30-seed Voronoi models in the
same training data set would allow for more general application of the GAN
where less prior knowledge of the material is known. We leave this for future
work.

Figure 11: (a) True high resolution (64x64) grain orientation maps
$\mathbf{m}_{64}^{\mbox{true}}$, (b) 16x16 resolution DNN tomography output
using AMSFMM generated data ($\mathbf{m}_{16}^{\mbox{pred}}$), and (c) 64x64
GAN output $\mathbf{m}_{64}^{G}$. For row (b), the MAE and SSIM are calculated
on an upscaled image to $64\times 64$ resolution using nearest neighbour
interpolation. It can be seen that the similar MAE and SSIM values result for
the Voronoi diagram arise since this type of texture was used in the training
data of the DNN and GAN. However, the GAN performs significantly worse in the
cases where the material texture is not part of the training data in columns
2,3 and 4.

Where real time inversions are not required, more computationally expensive
tomography algorithms can be implemented. Algorithms such as the rj-MCMC [41]
offer more information including an estimate of the uncertainty of the
tomography results. A place for rapid deep learning-based tomography still
exists within this framework as it can provide a fast, coarse initial model
which can be used a a starting point for more sophisticated algorithms.
Additionally a GAN can be used in post-processing any tomographic image. Often
linearised image methods are often regularised and hence predict smoother
structures than are expected to exist in the true medium, therefore a GAN can
be trained to upscale resolution and sharpen these images. Even where the GAN
provides marginal improvements to the DNN tomography results, the GAN output
models exhibit discontinuous boundaries. It can be important that such
boundaries are present in tomography algorithms where entire waveforms are
modelled and matched to the recorded waveforms (that is, full waveform
inversion [43]). A GAN might also be extended to take the full waveform as an
input, though this would require expensive FEA modelling to generate the
training data, so that all internal reflections are modelled.

##  5 Conclusion

We present a deep learning based framework for the real time tomographic
reconstruction of spatially varying crystal orientations in locally
anisotropic media using ultrasonic array time-of-flight data. We train a
series of deep neural networks (DNNs) using 7500 models in a training data
set, to accurately reconstruct orientation maps using full aperture, pitch-
catch and pulse-echo transducer array configurations. We present the first
application of generative adversarial networks (GANs) on ultrasonic
tomographic data, where a series of GANs are trained with three sets of
training data exhibiting increasing levels of complexity in the model
textures. The GAN takes the low resolution DNN output and upscales the
resolution by a factor of four. We show that prior information used to create
the training data for both the DNN and the GAN are important factors in
providing accurate estimations of the orientation maps. Using the methods
presented unlocks a wide range of potential applications for ultrasonic
monitoring, allowing for faster and more accurate detection of flaws and in-
process inspection during manufacturing.

##  6 Appendix

###  6.1 Finite element analysis

We implement a finite element simulation of elastic wave propagation in
anisotropic media using OnScale [34]. We apply absorbing boundary conditions
on all sides of the domain so energy continues past boundaries with no
reflections. We use Ricker wavelets with central frequencies of 1 MHz as the
source-time function, and apply pressure loads following the full aperture
transducer array configuration as shown in Figure 1(d). The values for the
finite element node spacing ($\Delta x,\Delta y$) are selected to ensure
spatial stability conditions following $\Delta x,\Delta y=\frac{\lambda}{15}$,
where $\lambda$ is the shortest wavelength in the domain.

Following the simulation for each transmitting array element, the travel time
to each receiving transducer is automatically picked by selecting the time for
arriving energy to increase above a threshold. This threshold is taken to be
2% of the peak displacement in the recorded signal.

###  6.2 Network Architectures

The deep neural networks (DNNs) are trained using 5 layers, where each node
receives an input from every node in the previous layer and a sigmoidal
activation function. The number of nodes in each layer are shown in Table 1.

Array configuration | No. of Inputs Nodes | $L1_{N}$ | $L2_{N}$ | $L3_{N}$ | No. of Trainable Parameters
---|---|---|---|---|---
Full Aperture | 96 | 315 | 63 | 63 | 58,591
Pitch-Catch | 256 | 354 | 55 | 55 | 113,911
Pulse-Echo | 136 | 353 | 68 | 68 | 72,511

Table 1: Network configurations showing the number of nodes for each layer
including the three hidden layers (L1-L3) for the full aperture, pitch-catch
and pulse-echo transducer array configurations.

Figure 12: Network configurations for (a) the GAN generator U-net architecture
and (b) the GAN discriminator PatchGAN architecture. All convolutional layers
have a kernel size of 4

The GAN generator is a modified U-Net based on [22] consisting of an encoder-
decoder chain. Each block in the encoder is a convolution-batch normalisation-
leaky rectified linear unit (ReLu) activation sequence. Each block in the
decoder is a transposed convolution -batch normalisation-ReLu sequence with
skip connections between mirrored layers in the encoder and decoder stacks
[38] (as shown in Fig. 12a). All convolutional layers use a kernel size of 4.
The generator loss is the discriminator sigmoid cross entropy loss of the
generates image with an array of ones combined with the mean absolute error
between the generated and known target image.

The GAN discriminator (Fig 12b) follows a PatchGAN architecture [22], which
divides the image into smaller 30x30 patches and the discriminator tries to
classify each patch separately. This motivates the GAN to discriminate high
frequency structure. The discriminator receives the target and generated
images as well as the low resolution input. The discriminator loss it the
sigmoid cross entropy loss with the real image and an array of ones, combined
with the sigmoid cross entropy loss with generated image and an array of
zeros.

###  6.3 Structural Similarity Index Measure (SSIM)

We use the SSIM described by [46] for image comparison. The SSIM is defined as
a weighted combination of comparisons between image luminance $l(X,Y)$,
contrast $c(X,Y)$ and structure $s(X,Y)$, where $X$ and $Y$ describe an image
window in known and estimated images of size $N\times N$. The SSIM is
therefore

$SSIM(X,Y)=[l(X,Y)]^{\alpha}\cdot[c(X,Y)]^{\beta}\cdot[s(X,Y)]^{\gamma}$ (5)

where $\alpha$, $\beta$ and $\gamma$ are the weighting parameters. We use
$\alpha=\beta=\gamma=1$. Luminance, contrast and structure are calculated as

$l(X,Y)=\frac{2\mu_{X}\mu_{Y}+C_{1}}{\mu_{X}^{2}+\mu_{Y}^{2}+C_{1}},$ (6)
$c(X,Y)=\frac{2\sigma_{X}\sigma_{Y}+C_{2}}{\sigma_{X}^{2}+\sigma_{Y}^{2}+C_{2}},$
(7) $s(X,Y)=\frac{\sigma_{XY}+C_{3}}{\sigma_{X}\sigma_{Y}+C_{3}}$ (8)

where $\mu$ and $\sigma$ are the mean and variance of the windows X or Y and
$\sigma_{XY}$ is the covariance of X and Y. This is computed over a sliding
Gaussian window of $9\times 9$.

##  7 Acknowledgments

This work was funded by the Engineering and Physical Sciences Research Council
(UK): grant number EP/P005268/1.

##  8 Data Availability

The data and Python scripts required to reproduce these findings are available
at:

https://github.com/jonnyrsingh/DeepLearningAnisoTomo, which can be executed
within Google colabotory. This requires no additional software or downloads
for the user.

## References

  * Aki and Lee [1976] K. Aki and W. Lee.  Determination of three-dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes: 1. A homogeneous initial model.  _Journal of Geophysical research_ , 81(23):4381–4399, 1976.
  * Aki et al. [1977] K. Aki, A. Christoffersson, and E. S. Husebye.  Determination of the three-dimensional seismic structure of the lithosphere.  _Journal of Geophysical Research_ , 82(2):277–296, 1977.
  * Amiri et al. [2020] N. Amiri, G. Farrahi, K. R. Kashyzadeh, and M. Chizari.  Applications of ultrasonic testing and machine learning methods to predict the static & fatigue behavior of spot-welded joints.  _Journal of Manufacturing Processes_ , 52:26–34, 2020.
  * Antholzer et al. [2019] S. Antholzer, M. Haltmeier, and J. Schwab.  Deep learning for photoacoustic tomography from sparse data.  _Inverse problems in science and engineering_ , 27(7):987–1005, 2019.
  * Araya-Polo et al. [2018] M. Araya-Polo, J. Jennings, A. Adler, and T. Dahlke.  Deep-learning tomography.  _The Leading Edge_ , 37(1):58–66, 2018.
  * Bergstra et al. [2015] J. Bergstra, B. Komer, C. Eliasmith, D. Yamins, and D. D. Cox.  Hyperopt: a Python library for model selection and hyperparameter optimization.  _Computational Science & Discovery_, 8(1):014008, 2015.
  * Bianco and Gerstoft [2018] M. J. Bianco and P. Gerstoft.  Travel time tomography with adaptive dictionaries.  _IEEE Transactions on Computational Imaging_ , 4(4):499–511, 2018.
  * Bishop et al. [1995] C. M. Bishop et al.  _Neural networks for pattern recognition_.  Oxford university press, 1995.
  * Bisong [2019] E. Bisong.  Google colaboratory.  In _Building Machine Learning and Deep Learning Models on Google Cloud Platform_ , pages 59–64. Springer, 2019.
  * Blitz and Simpson [1995] J. Blitz and G. Simpson.  _Ultrasonic methods of non-destructive testing_ , volume 2.  Springer Science & Business Media, 1995.
  * Bodin and Sambridge [2009] T. Bodin and M. Sambridge.  Seismic tomography with the reversible jump algorithm.  _Geophysical Journal International_ , 178(3):1411–1436, 2009.
  * Cao et al. [2020] R. Cao, S. Earp, S. A. de Ridder, A. Curtis, and E. Galetti.  Near-real-time near-surface 3D seismic velocity and uncertainty models by wavefield gradiometry and neural network inversion of ambient seismic noise.  _Geophysics_ , 85(1):KS13–KS27, 2020.
  * Capineri et al. [1992] L. Capineri, H. Tattersall, M. Silk, and J. Temple.  Time-of-flight diffraction tomography for NDT applications.  _Ultrasonics_ , 30(5):275–288, 1992.
  * Cha et al. [2017] Y.-J. Cha, W. Choi, and O. Büyüköztürk.  Deep learning-based crack damage detection using convolutional neural networks.  _Computer-Aided Civil and Infrastructure Engineering_ , 32(5):361–378, 2017.
  * Chassignole et al. [2000] B. Chassignole, D. Villard, M. Dubuget, J. Baboux, and R. E. Guerjouma.  Characterization of austenitic stainless steel welds for ultrasonic NDT.  In _AIP Conference Proceedings_ , volume 509, pages 1325–1332. American Institute of Physics, 2000.
  * Duric et al. [2005] N. Duric, P. Littrup, A. Babkin, D. Chambers, S. Azevedo, A. Kalinin, R. Pevzner, M. Tokarev, E. Holsapple, O. Rama, et al.  Development of ultrasound tomography for breast imaging: Technical assessment.  _Medical Physics_ , 32(5):1375–1386, 2005.
  * Earp and Curtis [2020] S. Earp and A. Curtis.  Probabilistic neural network-based 2D travel-time tomography.  _Neural Computing and Applications_ , 32(22):17077–17095, 2020.
  * Eltony et al. [2020] A. M. Eltony, P. Shao, and S.-H. Yun.  Measuring mechanical anisotropy of the cornea with Brillouin microscopy.  _arXiv preprint arXiv:2003.04344_ , 2020.
  * Guadarrama et al. [2017] S. Guadarrama, R. Dahl, D. Bieber, M. Norouzi, J. Shlens, and K. Murphy.  Pixcolor: Pixel recursive colorization.  _arXiv preprint arXiv:1705.07208_ , 2017.
  * Hoffmeister et al. [1994] B. K. Hoffmeister, E. D. Verdonk, S. A. Wickline, and J. G. Miller.  Effect of collagen on the anisotropy of quasi-longitudinal mode ultrasonic velocity in fibrous soft tissues: A comparison of fixed tendon and fixed myocardium.  _The Journal of the Acoustical Society of America_ , 96(4):1957–1964, 1994.
  * Huang et al. [2018] H.-w. Huang, Q.-t. Li, and D.-m. Zhang.  Deep learning based image recognition for crack and leakage defects of metro shield tunnel.  _Tunnelling and underground space technology_ , 77:166–176, 2018.
  * Isola et al. [2017] P. Isola, J.-Y. Zhu, T. Zhou, and A. A. Efros.  Image-to-image translation with conditional adversarial networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 1125–1134, 2017.
  * Javadi et al. [2019] Y. Javadi, C. N. MacLeod, S. G. Pierce, A. Gachagan, D. Lines, C. Mineo, J. Ding, S. Williams, M. Vasilev, E. Mohseni, et al.  Ultrasonic phased array inspection of a Wire+ Arc Additive Manufactured (WAAM) sample with intentionally embedded defects.  _Additive Manufacturing_ , 29:100806, 2019.
  * Javadi et al. [2020] Y. Javadi, N. E. Sweeney, E. Mohseni, C. N. MacLeod, D. Lines, M. Vasilev, Z. Qiu, R. K. Vithanage, C. Mineo, T. Stratoudaki, et al.  In-process calibration of a non-destructive testing system used for in-process inspection of multi-pass welding.  _Materials & Design_, 195:108981, 2020.
  * Jiang et al. [2019] K. Jiang, Z. Wang, P. Yi, G. Wang, T. Lu, and J. Jiang.  Edge-enhanced GAN for remote sensing image superresolution.  _IEEE Transactions on Geoscience and Remote Sensing_ , 57(8):5799–5812, 2019.
  * Kelly et al. [2017] B. Kelly, T. P. Matthews, and M. A. Anastasio.  Deep learning-guided image reconstruction from incomplete data.  _arXiv preprint arXiv:1709.00584_ , 2017.
  * Khairi et al. [2019] M. T. M. Khairi, S. Ibrahim, M. A. M. Yunus, M. Faramarzi, G. P. Sean, J. Pusppanathan, and A. Abid.  Ultrasound computed tomography for material inspection: Principles, design and applications.  _Measurement_ , 146:490–523, 2019.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Lebedev and Van Der Hilst [2008] S. Lebedev and R. D. Van Der Hilst.  Global upper-mantle tomography with the automated multimode inversion of surface and S-wave forms.  _Geophysical Journal International_ , 173(2):505–518, 2008.
  * Ledig et al. [2017] C. Ledig, L. Theis, F. Huszár, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang, et al.  Photo-realistic single image super-resolution using a generative adversarial network.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 4681–4690, 2017.
  * Lines [1998] D. Lines.  Rapid inspection using integrated ultrasonic arrays.  _Insight: Non-Destructive Testing and Condition Monitoring_ , 40(8):573–577, 1998.
  * Moya and Irikura [2010] A. Moya and K. Irikura.  Inversion of a velocity model using artificial neural networks.  _Computers & geosciences_, 36(12):1474–1483, 2010.
  * Nageswaran et al. [2009] C. Nageswaran, C. Carpentier, and Y. Tse.  Microstructural quantification, modelling and array ultrasonics to improve the inspection of austenitic welds.  _Insight-Non-Destructive Testing and Condition Monitoring_ , 51(12):660–666, 2009.
  * [34] OnScale.  770 Marshall Street, Redwood City, CA 94063.
  * Pyle et al. [2020] R. J. Pyle, R. L. Bevan, R. R. Hughes, R. K. Rachev, A. A. S. Ali, and P. D. Wilcox.  Deep learning for ultrasonic crack characterization in NDE.  _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , 2020.
  * Rawlinson and Sambridge [2004] N. Rawlinson and M. Sambridge.  Wave front evolution in strongly heterogeneous layered media using the fast marching method.  _Geophysical Journal International_ , 156(3):631–647, 2004.
  * Rodrigues et al. [2019] T. A. Rodrigues, V. Duarte, J. A. Avila, T. G. Santos, R. Miranda, and J. Oliveira.  Wire and arc additive manufacturing of HSLA steel: Effect of thermal cycles on microstructure and mechanical properties.  _Additive Manufacturing_ , 27:440–450, 2019.
  * Ronneberger et al. [2015] O. Ronneberger, P. Fischer, and T. Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _International Conference on Medical image computing and computer-assisted intervention_ , pages 234–241. Springer, 2015.
  * Senechal [1993] M. Senechal.  Spatial tessellations: Concepts and applications of Voronoi diagrams.  _Science_ , 260(5111):1170–1173, 1993.
  * Sethian [1999] J. A. Sethian.  _Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science_ , volume 3.  Cambridge university press, 1999.
  * Tant et al. [2018] K. M. Tant, E. Galetti, A. Mulholland, A. Curtis, and A. Gachagan.  A transdimensional Bayesian approach to ultrasonic travel-time tomography for non-destructive testing.  _Inverse Problems_ , 34(9):095002, 2018.
  * Tant et al. [2020] K. M. M. Tant, E. Galetti, A. Mulholland, A. Curtis, and A. Gachagan.  Effective grain orientation mapping of complex and locally anisotropic media for improved imaging in ultrasonic non-destructive testing.  _Inverse Problems in Science and Engineering_ , pages 1–25, 2020\.
  * Virieux and Operto [2009] J. Virieux and S. Operto.  An overview of full-waveform inversion in exploration geophysics.  _Geophysics_ , 74(6):WCC1–WCC26, 2009.
  * Virkkunen et al. [2021] I. Virkkunen, T. Koskinen, O. Jessen-Juhler, and J. Rinta-Aho.  Augmented ultrasonic data for machine learning.  _Journal of Nondestructive Evaluation_ , 40(1):1–11, 2021.
  * Wang et al. [2019] L. Wang, J. Xue, and Q. Wang.  Correlation between arc mode, microstructure, and mechanical properties during wire arc additive manufacturing of 316L stainless steel.  _Materials Science and Engineering: A_ , 751:183–190, 2019\.
  * Wang et al. [2004] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli.  Image quality assessment: from error visibility to structural similarity.  _IEEE transactions on image processing_ , 13(4):600–612, 2004.
  * Wirgin [2004] A. Wirgin.  The inverse crime.  _arXiv preprint math-ph/0401050_ , 2004.
  * Yoo et al. [2019] J. Yoo, S. Sabir, D. Heo, K. H. Kim, A. Wahab, Y. Choi, S.-I. Lee, E. Y. Chae, H. H. Kim, Y. M. Bae, et al.  Deep learning diffuse optical tomography.  _IEEE transactions on medical imaging_ , 39(4):877–887, 2019.
  * You et al. [2019] C. You, G. Li, Y. Zhang, X. Zhang, H. Shan, M. Li, S. Ju, Z. Zhao, Z. Zhang, W. Cong, et al.  CT super-resolution GAN constrained by the identical, residual, and cycle learning ensemble (GAN-CIRCLE).  _IEEE transactions on medical imaging_ , 39(1):188–203, 2019.
  * Zhang et al. [2012] J. Zhang, A. Hunter, B. W. Drinkwater, and P. D. Wilcox.  Monte Carlo inversion of ultrasonic array data to map anisotropic weld properties.  _IEEE transactions on ultrasonics, ferroelectrics, and frequency control_ , 59(11):2487–2497, 2012.
  * Zhang and Curtis [2020] X. Zhang and A. Curtis.  Variational full-waveform inversion.  _Geophysical Journal International_ , 222(1):406–411, 2020.
  * Zhang et al. [2018] X. Zhang, A. Curtis, E. Galetti, and S. De Ridder.  3-D Monte Carlo surface wave tomography.  _Geophysical Journal International_ , 215(3):1644–1658, 2018.
  * Zhao et al. [2020] X. Zhao, A. Curtis, and X. Zhang.  Bayesian seismic tomography using normalizing flows.  _EarthArXiv preprint https://doi.org/10.31223/X53K6G_ , 2020.
  * Zhu et al. [2017] H. Zhu, D. Komatitsch, and J. Tromp.  Radial anisotropy of the north american upper mantle based on adjoint tomography with usarray.  _Geophysical Journal International_ , 211(1):349–377, 2017.


# Characterization of passive CMOS sensors with RD53A pixel modules

Franz Glessgen1, Malte Backhaus1, Florencia Canelli2, Yannick Manuel Dieter3,
Jochen Christian Dingfelder3, Tomasz Hemperek3, Fabian Huegging3, Arash
Jofrehei2, Weijie Jin2, Ben Kilminster2, Anna Macchiolo2, Daniel
Muenstermann4, David-Leon Pohl3, Branislav Ristic1, Rainer Wallny1, Tianyang
Wang3, Norbert Wermes3, Pascal Wolf3 1 Institute for Particle Physics and
Astrophysics, ETH Zürich, Zürich, Switzerland  2 Universität Zürich, Zürich,
Switzerland  3 University of Bonn, Physikalisches Institut, Bonn, Germany  4
RheinMain University of Applied Sciences, Wiesbaden, Germany<EMAIL_ADDRESS>
###### Abstract

Both the current upgrades to accelerator-based HEP detectors (e.g. ATLAS, CMS)
and also future projects (e.g. CEPC, FCC) feature large-area silicon-based
tracking detectors. We are investigating the feasibility of using CMOS
foundries to fabricate silicon radiation detectors, both for pixels and for
large-area strip sensors. A successful proof of concept would open the market
potential of CMOS foundries to the HEP community, which would be most
beneficial in terms of availability, throughput and cost. In addition, the
availability of multi-layer routing of signals will provide the freedom to
optimize the sensor geometry and the performance, with biasing structures
implemented in poly-silicon layers and MIM-capacitors allowing for AC
coupling. A prototyping production of strip test structures and RD53A
compatible pixel sensors was recently completed at LFoundry in a 150nm CMOS
process. This presentation will focus on the characterization of pixel
modules, studying the performance in terms of charge collection, position
resolution and hit efficiency with measurements performed in the laboratory
and with beam tests. We will report on the investigation of RD53A modules with
25x100 $\upmu$m2 cell geometry.

##  1 Introduction

The High-Luminosity upgrade of the Large Hadron Collider (HL-LHC) will bring
the instantaneous luminosity of the LHC to $\mathcal{L}=5-7.5\times 10^{34}$
cm-2 s-1 and the pileup to between 140 and 200 for a target integrated
luminosity of 3 to 4 ab-1 over more than 10 years of operation. The CMS
tracker will have to endure a fluence of $2\times 10^{16}$ MeV
n${}_{\text{eq}}$.cm-2 for its first layer (at a radius of $r\approx 3$ cm
from the interaction region) and a TID of $12$ MGy [1]. To cope with the
expected irradiation levels and the increased trigger rate, the CMS Inner
Tracker (IT) will go through a complete replacement, known as the Phase II
upgrade. Pixel sensor prototypes were developed in collaboration with LFoundry
[2] to address the challenges of operating in the Phase II Tracker system. The
samples are passive planar n-in-p sensors for hybrid detector modules built in
CMOS technology using the 150 nm production line of LFoundry. CMOS technology
is largely used in the semiconductor industry and large-scale production of
sensors could significantly lower the cost and increase the production speed
compared to current standard processes utilized in HEP. The considered samples
have a 150 $\upmu$m thickness with a pixel size of 25$\times$100 $\upmu$m2 and
are DC coupled to the bump bonds of the RD53A readout chip (ROC). The RD53A
ROC is a prototype for the developement of the final readout chip for the
Phase II upgrade of the IT [3]. The measurements were partly done using the
DESY testbeam (TB) infrastructures [4].

##  2 CMOS sensors and stitching

The wafers produced by LFoundry are produced using stitching, a technology
that allows the designer to fabricate an image sensor that is larger than the
field of view of the lithographic equipment. This is achieved by implementing
different building blocks of the design in a reticule cell. By a dedicated
programming of the lithographic tool, each individual block of the reticle can
be selected and can be transferred into the photo-resist on the wafer,
possibly with different orientations. It is then possible to “stitch” the
various blocks together on the wafer during the lithographic process. In this
way, sensor wafers of several square cm of area, as normally implemented in
HEP experiments, can be built. Additionally, CMOS technology allows for the
implementation of very small on-pixel structures which gives access to new
sensor features (for HEP). One of them is the addition of metal layers on top
of the sensor for signal redistribution. Further, AC-coupled sensors can be
produced which cancels the need for leakage current compensation circuits in
the preamplifier of the ROC. Also, low and high resistivity polysilicon layers
can be deposited on the sensor to be used as field plates and bias resistors
respectively.

##  3 Qualification measurements

The sensors used in the IT upgrade have to match performance criteria. Those
are summarized in Table 1. The following measurements were performed to
compare the performance of the LFoundry sensors with the needed requirements.
All measurements and requirements are quoted for a sensor at a temperature of
20 oC.

Table 1: Phase II requirements for IT planar sensors. V${}_{\text{dep}}$ denotes the full depletion voltage.  Parameter | Requirement
---|---
Breakdown voltage |  $>$ 350 V
Leakage current |  $<$ 0.75 $\upmu$A.cm-2 at V${}_{\text{dep}}$ \+ 50 V
Efficiency |  $>$ 99 % at V${}_{\text{dep}}$ \+ 50 V
|

###  3.1 IV

The first two requirements in Table 1 constrain the IV curve of the sensors.
Figure 1 shows IV measurements of sensors from two different wafers. The bare
sensors from the first wafer showed a leakage current much higher than the
Phase II requirements. This was tracked down to a low implant dose in the
sensors backside. Two wafers were re-implanted to a higher dose before post-
processing for Under Bump Metallization (UBM). The IV curve of a sensor with a
higher implant concentration after flip-chipping is also shown in Figure 1.
The breakdown voltage and the leakage current at $V_{\text{dep}}$ \+ 50 V
comply with the CMS Phase II requirements.

Figure 1: IV measurements of LFoundry sensors before and after the increase of
the backside implant concentration. The dotted lines correspond to the sensors
with low backside implant concentration, measured at wafer level. The dashed
blue line is from the sensor with increased backside doping concentration,
measured at module level after flip-chipping. The maximal current allowed for
the Phase II requirements is indicated in light grey. All the IV measurements
were done at a temperature of 20 oC.

###  3.2 Efficiency

The efficiency is defined as the probability of detecting a hit on the Device
Under Test (DUT) within 500 $\upmu$m of each telescope track hit. The
efficiency of a non-irradiated LFoundry passive sensor, interconnected to a
RD53A chip, was investigated in a beam test at DESY [4] using the EUDET
telescope, with 5.2 GeV electrons. The module was read-out with the BDAQ53
system [5]. Figure 3 shows that the Phase II efficiency requirement is matched
from 10 V.

Figure 2: The blue markers show the efficiency of the sensor w.r.t the bias
voltage. The blue line shows the 99% efficiency limit and the red line is
positioned at the depletion voltage of the sensor.

Figure 3:  The histogram shows the number of pixels per row that did not
detect all telescope hits. The red line shows the same for the two rows
adjacent to the stitching line.

The sensor was processed using two reticle blocks stiched together (see
Section 2). To verify that the stitching process does not reduce the
efficiency, pixels next to the stitched part of the sensor (stitching line)
are compared to the rest of the pixel matrix. The number of pixels per row of
the sensor that were not fully efficient (under 85% efficiency for a mean
number of hits of 42) is plotted on Figure 3. The number of not fully
efficient pixels per row should be seen as a part of a set of 384 pixels. Some
rows reach an efficiency of 1 for all pixels but most rows show around 10 not-
fully-efficient pixels. The plot also shows this value for the rows adjacent
to the stitching line. It is consistent with the values of the non-stitched
rows and thus a strong argument for showing that stitching does not noticeably
reduce the efficiency.

###  3.3 Charge collection

The charge collected from a minimum ionizing particle in a thin silicon sensor
is of approximately 70 electrons per $\upmu$m, using Bichsel’s formula [6].
The expected charge collected in a 150 $\upmu$m thick sensor is thus
approximately 10500 electrons. The distribution of the charge collected in a
cluster has been obtained using the Time-over-Threshold (ToT) distribution
measured by the RD53A during a testbeam run at perpendicular incidence of the
beam on the DUT. The electron beam had an energy of 5.2 GeV. The ToT was then
converted to a physical charge using an X-ray source. Metal targets (Ag, Cu,
Mo, Sn, Zn) were irradiated by this source and emitted transition radiation of
known energy towards the RD53A. The deposited charge was then measured and the
relation between ToT and charge obtained. The most probable deposited charge
value (MPV) is obtained by fitting the cluster charge distribution with a
convolution of a Gaussian and a Landau distribution. The fit must then be
deconvolved to obtain the MPV of the Landau distribution. Figure 5 shows the
convoluted distribution and their fit for different bias voltages. Figure 5
shows the increase of the MPV after deconvolution w.r.t the bias voltage.

Figure 4: Cluster charge distribution with respect to the bias voltage. The
blue distribution corresponds to a bias of 10 V, the orange one to a bias of
20 V and the green one to a bias of 50 V. The MPVs can be found in Figure 5.

Figure 5: MPV of the collected charge as a function of the bias voltage after
deconvolution of the fitted function. The error bars come from the charge
calibration uncertainties.

###  3.4 Spatial resolution

The hit position on the DUT is independently computed in the two directions,
along the short and long pixel cell sides, on the sensor surface. The
following results focus on the resolution along the 25 $\upmu$m side of the
sensor. The hit position on the DUT is computed using the charge-weighted
center of the cluster [7]. The resolution of the telescope has to be taken
into account and projected on the sensor plane. Thus, the resolution of the
DUT is extracted from the formula

$\sigma_{\Delta
x}=\sqrt{\sigma_{\text{DUT}}^{2}+\left(\frac{\sigma_{\text{telescope}}}{\cos(\theta)}\right)^{2}}$
(1)

where $\theta$ is the angle between the beam axis and the sensor and the
standard deviation $\sigma_{\Delta x}$ comes from the distribution of the
residues between the telescope tracks and the hits on the DUT. Figure 6 shows
the resolution of the sensor as a function of the beam incidence angle as well
as the average cluster size along the 25 $\upmu$m side.

Figure 6:  The red markers show the DUT resolution w.r.t the beam incidence
angle in the direction of the short side of the pixel cell. The blue ones show
the mean cluster size in the same direction.

The minimal resolution of 2.2 $\upmu$m is reached between 8.75 and 9.75
degrees. It is compatible with the theoretical value of
$\arctan(\frac{p}{T})\approx 9.5$ degrees which is the angle maximizing the
number of clusters of size 2 in the 25 $\upmu$m direction. Here, $p$ denotes
the pixel pitch (25 $\upmu$m) and $T$ the thickness of the sensor.

###  3.5 Inter-pixel crosstalk

Capacitive coupling between neighbouring pixels (known as crosstalk) leads to
spurious hits and can deteriorate the sensor resolution. An enhanced cross-
talk has been measured on $25\times 100$ $\upmu$m2 sensors of other producers
with respect to the $50\times 50$ $\upmu$m2 geometry [8]. This was due to the
overlay or close vicinity of the metal system of the UBM of one pixel with the
implant layer of the neighbouring pixel, see Figure 8. Crosstalk can be
computed by injecting a known charge through the RD53A and reading out the
neighbouring pixels. Figure 8 shows the proportion of detected hits in pixels
that have received a charge injection as well as those that did not but are
crosstalk-coupled to the injected ones. The former detect all hits after their
threshold charge $T_{1}$ is crossed while the latter do not reach the 50 %
detection efficiency point (which defines the second threshold $T_{2}$) even
at the highest charge that can be injected using the RD53A. The crosstalk
$\alpha$ is defined as

$\alpha=\frac{1}{1+\frac{T_{2}}{T_{1}}}<2.5\%$ (2)

and is calculated by inserting the values of the measured threshold $T_{1}$
and the highest charge $T_{2}$ injected through the ROC.

Figure 7:  Mean proportion of detected hits w.r.t the injected charge. The red
curve corresponds to detected hits in the pixels that have been injected with
a charge and the blue curve to crosstalk coupled, neighbouring pixels, that
have not been injected with a charge.

Figure 8:  LFoundry pixel design showing the proximity of the bump bond pads
(octogonal shapes) to the implants (in blue).

##  4 Conclusion

The first laboratory and TB measurements of an unirradiated LFoundry passive
CMOS sensor flip-chipped to a RD53A ROC show promising results towards the
application in the Phase II upgrade of the CMS IT. No noticeable loss of
efficiency is found due to the stitching process. The crosstalk is shown to be
smaller than 2.5% and the minimal resolution along the short side of the pixel
cell of the sensor reaches 2.2 $\upmu$m. RD53A modules assembled with LFoundry
sensors will be characterized after irradiation to the highest fluence level
foreseen in the CMS IT system.

##  5 References

## References

  * [1] 2017 The Phase-2 Upgrade of the CMS Tracker Tech. Rep. CERN-LHCC-2017-009, CMS-TDR-014 CMS Collaboration CERN, Geneva
  * [2] LFoundry GmbH http://www.lfoundry.com , accessed: 2021-06-15
  * [3] Garcia-Sciveres M (RD53 Collaboration) 2017 The RD53A Integrated Circuit Tech. rep. CERN Geneva URL https://cds.cern.ch/record/2287593
  * [4] Diener R, Dreyling-Eschweiler J, Ehrlichmann H, Gregor I M, Kötz U, Krämer U, Meyners N, Potylitsina-Kube N, Schütz A, Schütze P et al. 2019 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 922 265–286
  * [5] Daas M, Dieter Y, Dingfelder J, Frohne M, Giakoustidis G, Hemperek T, Hinterkeuser F, Hügging F, Janssen J, Krüger H et al. 2021 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 986 164721
  * [6] Bichsel H 2006 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 562 154–197
  * [7] Turchetta R 1993 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 335 44–58
  * [8] Steinbrück, Georg and Group, CMS Tracker 2020 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 978 164438

The measurements leading to these results have been performed at the Test Beam
Facility at DESY Hamburg (Germany), a member of the Helmholtz Association
(HGF).


# Inflating and Reheating the Universe
with an Independent Affine Connection

Alberto Salvio Physics Department, University of Rome and INFN Tor Vergata,
Italy

###### Abstract

It is shown that a component of the dynamical affine connection, which is
independent of the metric, can drive inflation in agreement with observations.
This provides a geometrical origin for the inflaton. It is also found that the
decays of this field, which has spin 0 and odd parity, into Higgs bosons can
reheat the universe up to a sufficiently high temperature.

††journal: Physical Review D

###### Contents

  1. 1 Introduction
  2. 2 The key idea and inflation
  3. 3 Reheating
  4. 4 Conclusions
  5. A  Analytic expressions for inflationary quantities

##  1 Introduction

Einstein’s general relativity (GR) explains gravity in geometrical terms, the
distances are measured through the metric and the gravitational force is
determined by the (affine) connection, which is the essential building block
of covariant derivatives. This beautiful construction accounts for all
gravitational observations performed so far, including today’s nearly-
exponential accelerated expansion of the universe if the cosmological constant
is present.

It is generically accepted that another, but much more rapid, nearly-
exponential expansion occurred during the early stages of the universe
(inflation). This can be driven by a spin-0 field, the inflaton, with an
appropriate potential, which guarantees that such an expansion not only
occurred, but also eventually came to an end. Indeed, a reheating must take
place after inflation in order to generate all particles we observe.

From the purely geometrical point of view the metric and the connection,
unlike in GR, can be completely independent objects and, moreover, can contain
extra degrees of freedom besides the spin-2 graviton. This generalized
scenario is known as metric-affine gravity (see [1] for a recent discussion
and further references).

The goal of this paper is to discover whether the role of the inflaton can be
played by an extra dynamical component of the connection. The main motivation
behind this goal is to provide a geometrical origin for the inflaton too,
linking it to one of the essential geometrical objects, the connection. In
order to achieve this goal not only an inflaton with an appropriate potential
should be identified among the components of the connection, but also the
current constraints given by cosmic microwave background (CMB) observations
should be satisfied. These are provided by Planck [2] and, more recently, the
BICEP and Keck collaborations [3]. Moreover, as discussed, the universe must
be appropriately reheated after inflation. This requires an efficient
production of known particles, such as electrons, quarks and Higgs bosons.

In the following sections we show that all this is possible and we work out
the predictions in a simple, yet well-motivated, model.

##  2 The key idea and inflation

When the connection ${\cal A}_{\mu\leavevmode\nobreak\
\sigma}^{\leavevmode\nobreak\ \,\rho}$ and the metric, $g_{\mu\nu}$, are
independent there are two rather than one invariant that are linear in the
curvature

${\cal R}_{\mu\nu\leavevmode\nobreak\ \leavevmode\nobreak\
\sigma}^{\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\rho}\equiv\partial_{\mu}{\cal A}_{\nu\leavevmode\nobreak\
\sigma}^{\leavevmode\nobreak\ \,\rho}+{\cal A}_{\mu\leavevmode\nobreak\
\lambda}^{\leavevmode\nobreak\ \,\rho}{\cal A}_{\nu\leavevmode\nobreak\
\sigma}^{\leavevmode\nobreak\ \,\lambda}-(\mu\leftrightarrow\nu).$ (1)

The first one is the usual Ricci-like scalar111Greek indices are lowered and
raised by $g_{\mu\nu}$. ${\cal R}\equiv{\cal R}_{\mu\nu}^{\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \mu\nu}$ and the second one is the
parity-odd Holst invariant ${\cal
R^{\prime}}\equiv\epsilon^{\mu\nu\rho\sigma}{\cal
R}_{\mu\nu\rho\sigma}/\sqrt{-g}$ [4, 5, 6], where
$\epsilon^{\mu\nu\rho\sigma}$ is the totally antisymmetric Levi-Civita symbol
with $\epsilon^{0123}=1$ and $g$ is the determinant of the metric. In the GR
case, where ${\cal A}_{\mu\leavevmode\nobreak\ \sigma}^{\leavevmode\nobreak\
\,\rho}$ equals the Levi-Civita connection, ${\cal R}$ coincides with the
Ricci scalar, $R$, but ${\cal R^{\prime}}$ vanishes. For this reason in
metric-affine gravity ${\cal R^{\prime}}$ can be understood as a component of
the connection.

The key idea here is to identify the inflaton with ${\cal R^{\prime}}$. To do
so ${\cal R^{\prime}}$ has to be a dynamical field, which is independent of
the metric, and the simplest inflationary action that realizes this is

$S_{I}=\int d^{4}x\sqrt{-g}\left(\alpha{\cal R}+\beta{\cal R^{\prime}}+c{\cal
R^{\prime}}^{2}\right).$ (2)

Indeed, for $c=0$ one can easily show, by solving the connection equations,
that $S_{I}$ is equivalent to the Einstein-Hilbert action for any $\beta$,
having identified $\alpha=M_{P}^{2}/2$, where $M_{P}$ is the reduced Planck
mass. While for $c\neq 0$, standard auxiliary field methods show that an extra
spin-0 parity odd dynamical field $\zeta^{\prime}$, which is introduced as an
auxiliary field, is present and precisely equals ${\cal R^{\prime}}$ on shell
[7, 8, 9]: we can equivalently write

$S_{I}=\int d^{4}x\sqrt{-g}\left[\alpha{\cal
R}+\left(\beta+2c\zeta^{\prime}\right){\cal R^{\prime}}-c\zeta^{\prime
2}\right],$ (3)

which coincides with the expression in (2) after using the $\zeta^{\prime}$
equation. Because of its symmetry properties we call $\zeta^{\prime}$ the
pseudoscalaron. The $\beta{\cal R^{\prime}}$ term, known as the Holst term, is
also necessary to obtain a suitable inflaton potential, as we will see; the
quantity $M_{P}^{2}/(4\beta)$ is called the Barbero-Immirzi parameter [10,
11]. After using the ${\cal A}_{\mu\leavevmode\nobreak\
\sigma}^{\leavevmode\nobreak\ \,\rho}$ equation a non-canonical kinetic term
of $\zeta^{\prime}$ appears. It is possible to canonically normalize the
pseudoscalaron by considering the field redefinition

$\zeta^{\prime}(\omega)=\frac{1}{2c}\left(\frac{M_{P}^{2}\tanh
X(\omega)}{4\sqrt{1-\tanh^{2}X(\omega)}}-\beta\right),$ (4)

where

$X(\omega)\equiv\sqrt{\frac{2}{3}}\frac{\omega}{M_{P}}+\tanh^{-1}\left(\frac{4\beta}{\sqrt{16\beta^{2}+M_{P}^{4}}}\right),$
(5)

such that, after using the connection equations, $S_{I}$ becomes a standard
scalar-tensor action [9] (we use the mostly plus convention for the metric)

$S_{I}=\int
d^{4}x\sqrt{-g}\left[\frac{M_{P}^{2}}{2}R-\frac{(\partial\omega)^{2}}{2}-U(\zeta^{\prime}(\omega))\right],$
(6)

where $U(\zeta^{\prime}(\omega))=c\zeta^{\prime}(\omega)^{2}$ and clearly
$c>0$ for stability reasons. As clear from (6) this model does not contain any
ghost (the Einstein-Hilbert term has the usual sign and the kinetic term of
$\omega$ contributes positively to the kinetic energy). Furthermore, for $c>0$
the mass of $\omega$ (defined as the mass of the fluctuations of this field
around a Lorentz invariant solution) is positive, namely $\omega$ is not a
tachyon. The potential $U(\zeta^{\prime}(\omega))$ is symmetric in the
exchange $\\{\omega,\beta\\}\to\\{-\omega,-\beta\\}$. Therefore, an arbitrary
value of $\beta$ and its opposite are physically equivalent.

The $c{\cal R^{\prime}}^{2}$ term, besides being the simplest one leading to
an extra spin-0 field, is also motivated by scale invariance and Weyl
invariance at high energies: by replacing $g_{\mu\nu}\to\Omega^{2}g_{\mu\nu}$
that term is invariant, not only when $\Omega$ is spacetime independent (scale
invariance), but also when it is spacetime depend (Weyl invariance). This is
because in metric-affine gravity the metric and the connection are independent
and a rescaling of the metric does not imply any change in the connection. The
extension of the present model to a fully scale invariant one is beyond the
scope of the present work because is not mandatory to assess the viability of
this scenario: mass scales can also be added by hand. We thus leave such an
extension as an interesting outlook for future work. This may be realized
perhaps along the lines of [12] (where the vacuum expectation value of a
scalar field generates the mass scales, in our case $\alpha$ and $\beta$) or
[13] (where the mass scales are induced through a gravitational version of the
Coleman-Weinberg mechanism [14]). It is also interesting to note that the same
inflationary predictions generically emerge if one substitutes $c{\cal
R^{\prime}}^{2}$ with a general quadratic function of both ${\cal R}$ and
${\cal R^{\prime}}$, which is still compatible with scale invariance. This is
because such more general function leads to the same potential, as recently
shown in [9]. The inflationary predictions that we will find are, therefore,
quite robust.

The slow-roll approximation can be used when

$\epsilon\equiv\frac{M_{P}^{2}}{2}\left(\frac{1}{U}\frac{dU}{d\omega}\right)^{2}\ll
1,\quad\eta\equiv\frac{M_{P}^{2}}{U}\frac{d^{2}U}{d\omega^{2}}\ll 1$ (7)

and in this case the number of e-folds $N_{e}$ as a function of the field
$\omega$ is given by

$N_{e}(\omega)=N(\omega)-N(\omega_{\rm end}),$ (8)

where

$N(\omega)=\frac{1}{M_{P}^{2}}\int^{\omega}d\omega^{\prime}\,U\left(\frac{dU}{d\omega^{\prime}}\right)^{-1}$
(9)

and $\omega_{\rm end}$ satisfies $\epsilon(\omega_{\rm end})=1$ (see details
below). The scalar spectral index $n_{s}$, the tensor-to-scalar ratio $r$ and
the curvature power spectrum $P_{R}$ (at horizon exit) are then given by

$n_{s}=1-6\epsilon+2\eta,\quad r=16\epsilon,\quad
P_{R}=\frac{U/\epsilon}{24\pi^{2}M_{P}^{4}}.$ (10)

One finds analytic expressions not only for $\epsilon$, $\eta$, $n_{s}$, $r$
and $P_{R}$, but also for the e-fold functions $N$ and $N_{e}$ (see A).
Indeed, the equation $\epsilon(\omega_{\rm end})=1$ can be solved for real
$\omega_{\rm end}$ whenever $192\beta^{2}\geq 4M_{P}^{2}$ and one finds two
solutions, which we call $\omega_{\pm}$ and whose analytic expression is given
in Appendix A. Note that, as always the case in slow-roll inflation,
$\epsilon$, $\eta$, $n_{s}$ and $r$ are independent of the overall constant in
the potential ($1/c$ in this case), while $P_{R}$ is proportional to it. So
the observed value of $P_{R}$ (namely $(2.10\pm 0.03)\times 10^{-9}$ [2]) can
always be obtained by choosing $c$ appropriately.

Figure 1: Upper plot: the pseudoscalaron mass and the corresponding value of
$c$ (in the inset) that gives $P_{R}=2.1\times 10^{-9}$ [2] $N_{e}$ e-folds
before the end of inflation as a function of $\beta$. Down plot: the
corresponding pseudoscalaron potential for $\beta=-80M_{P}^{2}$. For all
curves $P_{R}=2.1\times 10^{-9}$ by construction. Also, the black points
correspond to the values of the inflaton for which $N_{e}=49$ (upper curve)
and $N_{e}=60$ (lower curve) are realized; the corresponding predictions for
$n_{s}$ and $r$ (in good agreement with the Planck, BICEP and Keck
observations) are provided.

We want $\omega_{\rm end}$ to be the value of $\omega$ at the end of inflation
so, given the shape of $U(\zeta^{\prime}(\omega))$, we take $\omega_{\rm end}$
such that $|\omega_{\rm end}|=\min(|\omega_{+}|,|\omega_{-}|)$. In Fig. 1 we
show the potential of $\omega$ (down plot) and its mass
$m_{\omega}=m_{\zeta^{\prime}}$ (upper plot) by setting $c$ in a way that
$P_{R}=2.1\times 10^{-9}$ at $N_{e}$ e-folds before the end of inflation. In
the $\omega$ potential there is a plateau, which increases for larger
$|\beta|$ and disappears when $\beta=0$. This is the reason why the
$\beta{\cal R^{\prime}}$ term in $S_{I}$ is necessary. In the down plot of
Fig. 1 $|\beta|=80M_{P}^{2}$ is chosen and is enough to even have 60 e-folds.

Figure 2: The scalar spectral index and the tensor-to-scalar ratio as
functions of the canonically normalized pseudoscalaron. The pseudoscalaron
values corresponding to $N_{e}=49,60$ e-folds before the end of inflation are
explicitly indicated: they correspond to the green points below the numbers
$49$ and $60$, respectively. In the inset the slow-roll parameters are shown.
We have set $\beta=-300M_{P}^{2}$.

In Figs. 2 and 3 it is shown that slow-roll inflation not only occurs, but is
also remarkably compatible with the most recent CMB observations provided by
Planck and BICEP/Keck (BK18 henceforth) for large $|\beta|$ (i.e. small values
of the Barbero-Immirzi parameter) and for an appropriate number of e-folds
$N_{e}$ [2]. Fig. 2 shows that viable slow-roll inflation with an appropriate
$N_{e}$ occurs for $\omega$ slightly above the Planck scale222However, the
corresponding values of the energy density, $\sim U$, is well below the cutoff
that is around the Planck scale (see the down plot of Fig. 1).; in that figure
$\beta=-300M_{P}^{2}$. In Fig. 3 we compare the observations and the
theoretical predictions as functions of $\beta$ and show that viable slow-roll
inflation with $N_{e}\simeq 49$ occurs already for $|\beta|\gtrsim
20M_{P}^{2}$ and with $N_{e}\simeq 60$ for $|\beta|\gtrsim 60M_{P}^{2}$. These
values of the mass parameter $\sqrt{|\beta|}$ are above $M_{P}$, but not much
larger than $M_{P}$: for $N_{e}\simeq 49$ and $N_{e}\simeq 60$ we have
respectively $\sqrt{|\beta|}\gtrsim 4M_{P}$ and $\sqrt{|\beta|}\gtrsim
8M_{P}$. In that figure $r_{0.002}$ is the value of $r$ at the reference
momentum scale $0.002\leavevmode\nobreak\ {\rm Mpc}^{-1}$, used by Planck and
BK18. In Figs. 3 we also report the predictions of Starobinsky inflation333In
Starobinsky inflation the inflationary action $S_{I}$ also features a
quadratic-in-curvature term, $S_{I}=\int
d^{4}x\sqrt{-g}\left(M_{P}^{2}R/2+cR^{2}\right)$, but with connection equal to
the Levi-Civita one. So it is interesting to compare the predictions of
pseudoscalaron inflation with those of Starobinsky inflation. [15] for $n_{s}$
and $r$; the predictions of pseudoscalaron inflation approach (but do not
quite reach) those of Starobinsky inflation for $|\beta|\to\infty$, while for
a finite value of $\beta$ they differ significantly.

It is interesting to note that the predictions for $n_{s}$ and $r$ of
pseudoscalaron inflation are within the reach of the future space mission
LiteBIRD [16], which will be, therefore, able to test this scenario.

Figure 3: The slow-roll parameters, the scalar spectral index and the tensor-
to-scalar ratio as functions of $\beta$. The green dots in the bottom plot are
the predictions of Starobinsky inflation.

##  3 Reheating

Reheating the universe after inflation is mandatory for the viability of any
model and, to achieve this, couplings between the inflaton and the Standard
Model (SM) particles are needed. If $\omega$ decays into some SM particles
with width $\Gamma_{\omega}$ the reheating temperature $T_{\rm RH}$ is at
least

$T_{\rm
RH}\gtrsim\min\left(\left(\frac{45\Gamma^{2}_{\omega}M_{P}^{2}}{4\pi^{3}g_{*}}\right)^{1/4},\left(\frac{30\rho_{\rm
vac}}{\pi^{2}g_{*}}\right)^{1/4}\right),$ (11)

where $g_{*}$ is the effective number of relativistic species in thermal
equilibrium at temperature $T_{\rm RH}$ and $\rho_{\rm vac}$ is the vacuum
energy density due to $\omega$ (note that $\rho_{\rm vac}$ represents the full
energy budget of the system). This is the standard perturbative contribution
to reheating; it is important to keep in mind that there may also be non-
perturbative contributions to the particle production after inflation [17, 18,
19, 20]. However, we leave their detailed calculation for future work because
they are not crucial to assess the viability of the present scenario. This is
because, as we will see, the value of $T_{\rm RH}$ estimated through the
standard perturbative approach can be large enough.

Let us first consider a fermion $f$ represented by a Dirac spinor $\Psi$
minimally coupled to gravity and with mass $m_{f}$, i.e. with action

$S_{f}=\int\sqrt{-g}\,\frac{1}{2}\overline{\Psi}(i\not{\cal D}-m_{f})\Psi+{\rm
h.c.}\,,$ (12)

where $\not{\cal D}\Psi\equiv\gamma^{a}e_{a}^{\mu}{\cal D}_{\mu}\Psi$ (the
$e_{a}^{\mu}$ satisfy [21] $e_{a}^{\mu}e_{b}^{\nu}g_{\mu\nu}=\eta_{ab}$),
$\overline{\Psi}\equiv\Psi^{\dagger}\gamma^{0}$, the Dirac gamma matrices
$\gamma^{a}$ satisfy $\\{\gamma^{a},\gamma^{b}\\}=-2\eta^{ab}$, ${\cal
D}_{\mu}\Psi=\partial_{\mu}\Psi+{\cal
A}^{ab}_{\mu}[\gamma_{a},\gamma_{b}]\Psi/8$ and ${\cal
A}_{\mu\leavevmode\nobreak\ b}^{\leavevmode\nobreak\ \,a}=e^{a}_{\nu}{\cal
A}_{\mu\leavevmode\nobreak\ \lambda}^{\leavevmode\nobreak\
\,\nu}e^{\lambda}_{b}-e^{\lambda}_{b}\partial_{\mu}e^{a}_{\lambda}$. By using
the connection equations with the formalism of [9], one finds the following
effective pseudoscalaron-fermion-fermion interaction

$\mathscr{L}_{\omega ff}=\frac{c_{\omega
ff}}{M_{P}}\,\partial_{\mu}\omega\,\overline{\Psi}\gamma_{5}\gamma^{\mu}\Psi,$
(13)

where

$c_{\omega
ff}=\left[\frac{3M_{P}}{1+16B^{2}}\frac{dB}{d\omega}\right]_{\omega=0}=\sqrt{\frac{3M_{P}^{4}}{8(M_{P}^{4}+16\beta^{2})}},$
(14)

$\gamma_{5}=i\gamma_{0}\gamma_{1}\gamma_{2}\gamma_{3}$ and
$B=(\beta+2c\zeta^{\prime}(\omega))/M_{P}^{2}$. This effective interaction
leads to the decay $\omega\to ff$ with width

$\Gamma_{\omega\to ff}=|c_{\omega ff}|^{2}\frac{m_{\omega}m_{f}^{2}}{2\pi
M_{P}^{2}}\sqrt{1-\frac{4m^{2}_{f}}{m_{\omega}^{2}}}.$ (15)

This channel can efficiently reheat the universe up to a temperature above the
electroweak scale if $m_{f}$ is very large compared to that scale. Such a
fermion is not present in the SM. It is possible to engineer a model where
there is a very heavy fermion with sizable couplings to SM particles such that
this channel is sufficient. For example, this is the case in the well-
motivated model [22, 23], which was proposed to solve the strong CP problem.

However, in order to keep our analysis as model independent as possible, we
consider another channel: the decay of $\omega$ into two identical real scalar
particles, e.g. two Higgs bosons. This channel can be active when there is a
non-minimal coupling between the real (canonically normalized) scalar field
$\phi$ in question and ${\cal R}$ in the action:

$S_{\rm nm}=\int\sqrt{-g}\frac{\xi\phi^{2}}{2}{\cal R}.$ (16)

This term is known to be generated by quantum corrections and, therefore, it
is more natural to include it. If one solves the connection equation in the
presence of (16) (using the results in [9]) one finds the following effective
pseudoscalaron-scalar-scalar interaction

$\mathscr{L}_{\omega\phi\phi}=\frac{c_{\omega\phi\phi}}{M_{P}}\,\partial_{\mu}\omega\,\phi\partial^{\mu}\phi,$
(17)

where

$c_{\omega\phi\phi}=\left[\frac{48\xi
M_{P}B}{1+16B^{2}}\,\frac{dB}{d\omega}\right]_{\omega=0}=\frac{4\sqrt{6}\beta\xi}{\sqrt{M_{P}^{4}+16\beta^{2}}}.$
(18)

This effective parity-violating operator only arises in the presence of the
Holst term because $c_{\omega\phi\phi}\to 0$ as $\beta\to 0$. The effective
interaction $\mathscr{L}_{\omega\phi\phi}$ leads to the decay
$\omega\to\phi\phi$ with width

$\Gamma_{\omega\to\phi\phi}=|c_{\omega\phi\phi}|^{2}\frac{m_{\omega}^{3}}{16\pi
M_{P}^{2}}\sqrt{1-\frac{4m^{2}_{\phi}}{m_{\omega}^{2}}},$ (19)

where $m_{\phi}$ is the mass of $\phi$. The produced Higgs particles
subsequently decay into other SM particles, such as leptons and quarks. The
channel $\omega\to\phi\phi$ can efficiently and naturally reheat the universe
up to a temperature much above the electroweak scale, even if one identifies
$\phi$ with the SM Higgs, so per se it does not require any beyond-the-SM
physics. For example, taking $m_{\phi}\ll m_{\omega}$, $g_{*}\sim 10^{2}$ and
$\beta\gtrsim M_{P}^{2}$ one finds $T_{\rm RH}\gtrsim 10^{9}|\xi|$ GeV. This
reheating temperature is compatible with all numbers of e-folds considered in
Sec. 2 for natural values of $|\xi|$ of order 1 or smaller. Since
$c_{\omega\phi\phi}\to 0$ as $\beta\to 0$ this reheating channel occurs thanks
to the presence of an independent connection: the Holst term would be absent
if the full connection were exactly the Levi-Civita one.

##  4 Conclusions

It has been found that a pseudoscalar component of a dynamical connection,
which is independent of the metric, can drive inflation in agreement with
current data. This pseudoscalaron is identified with the parity odd Holst
invariant and inflationary predictions in excellent agreement with data have
been found for small values of the Barbero-Immirzi parameter, where the
inflaton potential develops a plateau. The predictions approach, but do not
quite reach, those of Starobinsky inflation as the Barbero-Immirzi parameter
goes to zero; for finite values, on the other hand, the predictions
significantly differ. Pseudoscalaron inflation can be tested by future CMB
observations, such as those of LiteBIRD.

Moreover, the decays of the pseudoscalaron into Higgs particles (which occur
thanks to the presence of an independent connection) can efficiently reheat
the universe after inflation up to a high enough temperature. This temperature
could be further increased by other channels, such as decays into very massive
fermions, which we have computed too.

As an outlook example, it would be interesting to calculate the non-
perturbative particle production after inflation (preheating). Moreover, it
would also be interesting to engineer a fully scale invariant version of this
model. Indeed, a crucial ingredient of the present construction is a
quadratic-in-curvature term, $c{\cal R^{\prime}}^{2}$, which is compatible
with scale (and even Weyl) invariance.

## Acknowledgments

I thank G. Pradisi for useful discussions. This work has been partially
supported by the grant DyConn from the University of Rome Tor Vergata.

##  Appendix A
Analytic expressions for inflationary quantities

The analytic expressions for $\epsilon$, $\eta$, $n_{s}$, $r$,$P_{R}$ and $N$
are

$\displaystyle\epsilon(\omega)$ $\displaystyle=$
$\displaystyle\frac{4M_{P}^{4}\cosh^{2}X(\omega)}{3\left(M_{P}^{2}\sinh
X(\omega)-4\beta\right)^{2}},$ $\displaystyle\eta(\omega)$ $\displaystyle=$
$\displaystyle\frac{4M_{P}^{2}\left(M_{P}^{2}\cosh\left(2X(\omega)\right)-4\beta\sinh
X(\omega)\right)}{3\left(M_{P}^{2}\sinh X(\omega)-4\beta\right)^{2}},$
$\displaystyle N(\omega)$ $\displaystyle=$
$\displaystyle\frac{3}{4}\log\left(\cosh
X(\omega)\right)-\frac{3\beta\arctan\left(\sinh X(\omega)\right)}{M_{P}^{2}},$
$\displaystyle n_{s}(\omega)$ $\displaystyle=$ $\displaystyle
1-\frac{8M_{P}^{4}\cosh^{2}X(\omega)}{\left(M_{P}^{2}\sinh
X(\omega)-4\beta\right)^{2}}$
$\displaystyle+\frac{8M_{P}^{2}\left(M_{P}^{2}\cosh\left(2X(\omega)\right)-4\beta\sinh
X(\omega)\right)}{3\left(M_{P}^{2}\sinh X(\omega)-4\beta\right)^{2}},$
$\displaystyle r(\omega)$ $\displaystyle=$
$\displaystyle\frac{64M_{P}^{4}\cosh^{2}X(\omega)}{3\left(M_{P}^{2}\sinh
X(\omega)-4\beta\right)^{2}},$ $\displaystyle P_{R}(\omega)$ $\displaystyle=$
$\displaystyle\frac{\left(\beta-\frac{M_{P}^{2}\sinh
X(\omega)}{4}\right)^{2}\left(M_{P}^{2}\sinh
X(\omega)-4\beta\right)^{2}\text{sech}^{2}X(\omega)}{128\pi^{2}cM_{P}^{8}}.$

Moreover, the analytic expressions of $\omega_{\pm}$ (the two solutions of
$\epsilon(\omega_{\rm end})=1$) are

$\displaystyle\omega_{\pm}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{3}{2}}M_{P}\left(\sinh^{-1}\left(\pm\sqrt{\frac{192\beta^{2}}{M_{P}^{4}}-4}-\frac{12\beta}{M_{P}^{2}}\right)\right.$
$\displaystyle\left.-\tanh^{-1}\left(\frac{4\beta}{\sqrt{16\beta^{2}+M_{P}^{4}}}\right)\right).$

## References

## References

  * [1] A. Baldazzi, O. Melichev and R. Percacci, “Metric-Affine Gravity as an effective field theory,” Annals Phys. 438 (2022), 168757 doi:10.1016/j.aop.2022.168757 [arXiv:2112.10193].
  * [2] P. A. R. Ade et al. [Planck Collaboration], “Planck 2015 results. XX. Constraints on inflation,” Astron. Astrophys. 594 (2016) A20 doi:10.1051/0004-6361/201525898 [arXiv:1502.02114]. Y. Akrami et al. [Planck Collaboration], “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020), A10 doi:10.1051/0004-6361/201833887 [arXiv:1807.06211].
  * [3] P. A. R. Ade et al. [BICEP and Keck], “Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season,” Phys. Rev. Lett. 127 (2021) no.15, 151301 doi:10.1103/PhysRevLett.127.151301 [arXiv:2110.00483].
  * [4] R. Hojman, C. Mukku and W. A. Sayed, “Parity violation in metric torsion theories of gravitation,” Phys. Rev. D 22 (1980), 1915-1921 doi:10.1103/PhysRevD.22.1915
  * [5] P. C. Nelson, “Gravity With Propagating Pseudoscalar Torsion,” Phys. Lett. A 79 (1980), 285 doi:10.1016/0375-9601(80)90348-5
  * [6] S. Holst, “Barbero’s Hamiltonian derived from a generalized Hilbert-Palatini action,” Phys. Rev. D 53 (1996), 5966-5969 doi:10.1103/PhysRevD.53.5966 [arXiv:gr-qc/9511026].
  * [7] R. D. Hecht, J. M. Nester and V. V. Zhytnikov, “Some Poincare gauge theory Lagrangians with well posed initial value problems,” Phys. Lett. A 222 (1996), 37-42 doi:10.1016/0375-9601(96)00622-6
  * [8] J. Beltrán Jiménez and F. J. Maldonado Torralba, “Revisiting the stability of quadratic Poincaré gauge gravity,” Eur. Phys. J. C 80 (2020) no.7, 611 doi:10.1140/epjc/s10052-020-8163-8 [arXiv:1910.07506].
  * [9] G. Pradisi and A. Salvio, “(In)equivalence of metric-affine and metric effective field theories,” Eur. Phys. J. C 82 (2022) no.9, 840 doi:10.1140/epjc/s10052-022-10825-9 [arXiv:2206.15041].
  * [10] G. Immirzi, “Real and complex connections for canonical gravity,” Class. Quant. Grav. 14 (1997), L177-L181 doi:10.1088/0264-9381/14/10/002 [arXiv:gr-qc/9612030].
  * [11] G. Immirzi, “Quantum gravity and Regge calculus,” Nucl. Phys. B Proc. Suppl. 57 (1997), 65-72 doi:10.1016/S0920-5632(97)00354-X [arXiv:gr-qc/9701052].
  * [12] M. Shaposhnikov and D. Zenhausern, “Quantum scale invariance, cosmological constant and hierarchy problem,” Phys. Lett. B 671 (2009), 162-166 doi:10.1016/j.physletb.2008.11.041 [arXiv:0809.3406].
  * [13] A. Salvio and A. Strumia, “Agravity,” JHEP 1406 (2014) 080 doi:10.1007/JHEP06(2014)080 [arXiv:1403.4226]. A. Salvio, “Dimensional Transmutation in Gravity and Cosmology,” Int. J. Mod. Phys. A 36 (2021) no.08n09, 2130006 doi:10.1142/S0217751X21300064 [arXiv:2012.11608].
  * [14] S. R. Coleman and E. J. Weinberg, “Radiative Corrections as the Origin of Spontaneous Symmetry Breaking,” Phys. Rev. D 7 (1973) 1888.
  * [15] A. A. Starobinsky, “A New Type of Isotropic Cosmological Models Without Singularity,” Phys. Lett. B 91 (1980), 99-102 doi:10.1016/0370-2693(80)90670-X.

  * [16] E. Allys et al. [LiteBIRD], “Probing Cosmic Inflation with the LiteBIRD Cosmic Microwave Background Polarization Survey,” [arXiv:2202.02773].
  * [17] A. D. Dolgov and D. P. Kirilova, “On particle creation by a time dependent scalar field,” Sov. J. Nucl. Phys. 51 (1990), 172-177 JINR-E2-89-321.
  * [18] J. H. Traschen and R. H. Brandenberger, “Particle Production During Out-of-equilibrium Phase Transitions,” Phys. Rev. D 42 (1990), 2491-2504 doi:10.1103/PhysRevD.42.2491
  * [19] L. Kofman, A. D. Linde and A. A. Starobinsky, “Reheating after inflation,” Phys. Rev. Lett. 73 (1994), 3195-3198 doi:10.1103/PhysRevLett.73.3195 [arXiv:hep-th/9405187].
  * [20] L. Kofman, A. D. Linde and A. A. Starobinsky, “Towards the theory of reheating after inflation,” Phys. Rev. D 56 (1997), 3258-3295 doi:10.1103/PhysRevD.56.3258 [arXiv:hep-ph/9704452].
  * [21] $\eta_{ab}$ is the Minkowski metric, which we use to lower and raise the indices $a,b$.
  * [22] J. E. Kim, “Weak Interaction Singlet and Strong CP Invariance,” Phys. Rev. Lett. 43 (1979), 103 doi:10.1103/PhysRevLett.43.103
  * [23] M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, “Can Confinement Ensure Natural CP Invariance of Strong Interactions?,” Nucl. Phys. B 166 (1980), 493-506 doi:10.1016/0550-3213(80)90209-6