# Optimising hadronic collider simulations using amplitude neural networks

Ryan Moodie  Institute for Particle Physics Phenomenology, Ogden Centre for
Fundamental Physics, Department of Physics, University of Durham, South Road,
Durham, DH1 3LE, UK<EMAIL_ADDRESS>
###### Abstract

Precision phenomenological studies of high-multiplicity scattering processes
at collider experiments present a substantial theoretical challenge and are
vitally important ingredients in experimental measurements. Machine learning
technology has the potential to dramatically optimise simulations for
complicated final states. We investigate the use of neural networks to
approximate matrix elements, studying the case of loop-induced diphoton-plus-
jets production through gluon fusion. We train neural network models on one-
loop amplitudes from the NJet C++ library and interface them with the Sherpa
Monte Carlo event generator to provide the matrix element within a realistic
hadronic collider simulation. Computing some standard observables with the
models and comparing to conventional techniques, we find excellent agreement
in the distributions and a reduced total simulation time by a factor of
thirty.

IR

    infrared
ME

    matrix element
ML

    machine learning
MC

    Monte Carlo
NN

    neural network
SM

    Standard Model
LO

    leading order
NLO

    next-to-leading order
NNLO

    next-to-next-to-leading order
QCD

    quantum chromodynamics
PDF

    parton distribution function

##  1 Introduction

With the increasing Large Hadron Collider dataset driving ever more precise
experimental measurements, Standard Model (SM) predictions for high-
multiplicity scattering at hadronic colliders form a vital part of precision
phenomenology studies. Currently, these calculations mainly rely on automated
numerical codes [1] to calculate matrix elements up to high multiplicities,
including tree-level real corrections at next-to-leading order (NLO) and
double-real corrections at next-to-next-to-leading order (NNLO), and one-loop
real-virtual corrections at NNLO. These codes have been a theoretical
revolution, particularly at one-loop, but the evaluation time is relatively
high. Due to the high dimensionality of the phase space, these real-type
corrections are often the computational bottleneck in higher-order
calculations.

Following recent advances in precision quantum chromodynamics (QCD), there has
been a flurry of activity around NNLO QCD corrections to diphoton-plus-jet
production, including full-colour two-loop amplitudes [2] and leading-colour
NNLO distributions [3]. In the loop-induced gluon fusion channel, the full-
colour two-loop amplitudes were computed [4], leading to full-colour NLO
distributions [5]; it was found that the opening of this new channel enhances
the mixed quark- and gluon-initiated NNLO prediction [3] in certain regions,
disrupting the usual perturbative convergence. Also for diphoton production
through gluon fusion, the three-loop amplitudes were calculated [6], making
available the final remaining piece for its NNLO corrections. Therefore, we
study the loop-induced class of processes with two photons and many gluons as
they are extremely relevant for current phenomenology.

ML technology has found a wealth of application in high energy physics [7]. We
employ the ensemble neural network (NN) model of Ref. [8], which studied
$e^{+}e^{-}$ annihilation to jets, to emulate the gluon-initiated diphoton-
plus-jets MEs within a full Monte Carlo (MC) event generator simulation. This
tests the methodology against the additional complexity of hadronic collider
simulations, including parton distribution function (PDF) convolution and
variable centre-of-mass scales, complex phase space cuts and jet clustering,
and phase space sampling optimisation methods of integrators.

This contribution is organised as follows. We first discuss the gluon-
initiated diphoton-plus-jets amplitudes and their implementations within the
C++ NJet library [9, 10] which were used for training. We then describe the
phase space partitioning used to handle infrared (IR) divergent regions. Next,
we present the architecture of the NNs used. Then, we discuss the simulation
pipeline and interface of the NN model to the Sherpa MC event generator [11].
Finally, we study the performance of the model compared to the original
amplitude library for $gg\to\gamma\gamma gg$ and present some distributions
before concluding.

This contribution is based on Ref. [12]. Our code is publicly available [13].

##  2 Amplitudes

Figure 1:  Diagram of $gg\to\gamma\gamma gg$ ($N=6$) at LO. The photons couple
to an internal quark loop.

Figure 2:  Single-threaded CPU evaluation time of the ME for a single phase
space point. The value is the mean of 100 evaluations over a random phase
space. Results are shown for available implementations at various
multiplicities, including numerical and analytical evaluations using NJet and
inference on the NN model.

As there is no vertex coupling gluons to photons in the SM, diphoton-plus-jets
production through gluon fusion (Fig. 2) is loop induced. The leading order
(LO) process is $\mathcal{O}({\alpha_{s}}^{N-2})$ for multiplicity $N$,
appearing at NNLO in the perturbative expansion of the combined quark- and
gluon-initiated process. We study the channels with only gluons and photons in
the external particles, $gg\to\gamma\gamma+n\times g$. These proceed through a
quark loop at LO.

Conventional event generator simulations optimise virtual corrections in NLO
calculations by learning the phase space of the LO process and using this to
sample the virtual contribution. This technique fails for loop-induced
processes, where the expensive one-loop amplitude has no tree-level process to
optimise the phase space on. Therefore, new methods are required to improve
the efficiency of integrating these channels at high multiplicity.

We use the one-loop-squared MEs from the NJet library as the targets for our
NN emulation. These include two classes of amplitudes: an automated numerical
setup for arbitrary multiplicity; and hard-coded analytical expressions for
$N\in\\{4,5\\}$. The numerical implementation obtains the diphoton-plus-jets
amplitudes by summing permutations of pure-gluon primitive amplitudes [14],
which are themselves based on generalised unitary [15] and integrand reduction
[16]. While completely automated, evaluation time and numerical stability are
increasingly difficult to control. The hard-coded implementations offer
compact analytical expressions with extremely fast and stable evaluation,
although they are unavailable for higher multiplicity. The $N=5$ result is
obtained through a finite field reconstruction [17]. The evaluation timings of
these methods are compared to the NN model in Fig. 2.

##  3 Phase space partitioning

Training a single NN over the entire phase space results in a poor fit,
especially at higher multiplicity [8]. This is caused by regions where the ME
becomes IR divergent, which arise from soft ($s_{i}$) and collinear ($c_{ij}$)
emissions. These singularities are regulated with cuts, but the ME is rapidly
varying in local regions, which causes problems for the global fit. Therefore,
we train a separate NN on each of the IR structures of the phase space.

We first partition the phase space into a non-divergent region,
$\mathcal{R}_{\textrm{non-div}}$, and a divergent region,
$\mathcal{R}_{\textrm{div}}$. Points that pass a cut,
$\text{min}(\\{s_{ij}/s_{12}:i,j\in\\{1,\ldots,N\\}\\})<y$, are included in
$\mathcal{R}_{\textrm{div}}$. The threshold $y$ must be tuned to discriminate
points of similar scales into each region, while having sufficient points in
$\mathcal{R}_{\textrm{div}}$ to train on.

We then sub-divide $\mathcal{R}_{\textrm{div}}$ according to the decomposition
of the FKS subtraction scheme [18]. We define a set of FKS pairs,
$\mathcal{P}_{\text{FKS}}=\left\\{(i,j):s_{i}\lor s_{j}\lor c_{ij}\right\\}$,
corresponding to the $\binom{N}{2}-1$ singular configurations, which includes
redundancy (App. B of Ref. [8]). Each pair is assigned a partition function,
$\mathcal{S}_{ij}=1/\left(s_{ij}\sum_{j,k\in\mathcal{P}_{\text{FKS}}}1/s_{jk}\right),$
which smoothly isolates that divergence on multiplication with the ME.

We train a NN on $\left|\mathcal{A(\boldsymbol{p})}\right|^{2}$ for
$\boldsymbol{p}\in\mathcal{R}_{\textrm{non-div}}$, and a NN on each of the
partition-function-weighted MEs,
$\left\\{\mathcal{S}_{ij}\left|\mathcal{A(\boldsymbol{p})}\right|^{2}\,:\,i,j\in\mathcal{P}_{\text{FKS}}\,;\,\boldsymbol{p}\in\mathcal{R}_{\textrm{div}}\right\\}.$
We reconstruct the complete ME in $\mathcal{R}_{\textrm{div}}$ by summing the
weighted MEs,
$\left|\mathcal{A}\right|^{2}=\sum_{i,j\in\mathcal{P}_{\text{FKS}}}\mathcal{S}_{ij}\left|\mathcal{A}\right|^{2}.$
This ensemble of NNs, referred to as the model, can be used to accurately
infer the ME over the complete phase space, $\mathcal{R}_{\textrm{non-
div}}\cup\mathcal{R}_{\textrm{div}}$.

Note that increasing the cut $y$, which increases the proportion of points in
$\mathcal{R}_{\textrm{div}}$, incurs a performance penalty due to the higher
cost of inferring over several NNs in $\mathcal{R}_{\textrm{div}}$ compared to
the single NN in $\mathcal{R}_{\textrm{non-div}}$.

##  4 Model architecture

Although using fine-tuned architectures for each configuration (process, cuts,
etc.) would provide optimal performance, this would be prohibitively
expensive. We use a general setup as the gains of specialised NN optimisation
are beyond the scope of this pioneering work, performing hyperparameter
optimisation on the $gg\to\gamma\gamma g$ process.

Each NN uses a fully-connected architecture, a standard choice for a
regression problem, parameterised using the Keras Python interface [19] to the
TensorFlow machine learning (ML) library [20]. There are $4\times N$ input
nodes: one for each component of each momentum in the phase space point. The
three hidden layers are comprised of 20, 40, and 20 nodes respectively, all
with hyperbolic-tangent activation functions. There is a single output node
with a linear activation function, which returns the approximation of the ME.
We find that this is a sufficient number of layers and nodes to learn the MEs,
while remaining economical for computational performance.

We train with a mean-squared-error loss function, using Adam-optimised
stochastic gradient descent [21]. We expect that the model will learn the mean
of the target distribution using this loss function (App. A of Ref. [8]). The
number of training epochs is determined by Early Stopping regularisation, with
a patience of 100 epochs to mitigate the effects of the limited size of
$\mathcal{R}_{\text{div}}$ that may appear in the validation set. We use
32-bit floating-point numbers throughout.

##  5 Pipeline

Our ML pipeline for the $gg\to\gamma\gamma gg$ results presented is: generate
the training and validation datasets by running Sherpa with NJet while
uniformly sampling phase space; train the model; infer on the model to
estimate the MEs during event generation with Sherpa, using the original
integration grid.

Input data consists of a list of phase space points,
$p_{i}^{\mu}\in\mathbb{R}^{4N}$, and the corresponding colour- and helicity-
summed one-loop-squared ME, $\left|\mathcal{A}\right|^{2}\in\mathbb{R}$. Phase
space sampling is determined by the integrator, meaning the training is
optimal only for a specific integrator. The data is extracted from a run of
the integrator, generating 100k points which are split 4:1 into training and
validation datasets. A 3M point testing dataset is produced by a second run of
the integrator with a different random number seed and used to evaluate model
performance.

We perform inference on an ensemble of twenty models, each of which have
different random weight initialisation and shuffled training and validation
datasets. We take as the result the mean of the ensemble, with the standard
error providing the precision/optimality error [8].

While training was performed using Python, event generators are generally
written in C++. To use the model within a simulation, we wrote a C++ inference
code and a bespoke C++ interface for Sherpa. The weights of the trained models
are written to file and read by the inference code at runtime; the library
Eigen [22] is used to perform efficient linear algebra on the CPU. The
interface can also be used to call C++ amplitude libraries directly; we use
this to interface NJet to Sherpa to generate the datasets, which is performed
with 64-bit floats.

PDF are provided by LHAPDF [23] using the NNPDF3.1 set NNPDF31_nlo_as_0118
[24]. Cuts are adapted from those in Ref. [25]. Analysis is performed using
Rivet [26] with an adapted reference analysis script [27].

##  6 Results

Figure 3:  PDF-weighted and unweighted histograms by region of the logarithm
of the ratio between the ME returned by the model and NJet for each point in a
1M subset of the training data for $gg\to\gamma\gamma gg$. The region cut is
$y=10^{-3}$ and $\mathcal{R}_{\textrm{div}}$ contains 2.4% of the points.

Figure 4:  Histogram of the root mean squared error (RMSE) of the model
compared to NJet for a univariate slice of phase space in $x_{2}$, the
momentum fraction of the second incoming gluon, for $gg\to\gamma\gamma g$.
Also shown are the points in the training dataset, binned in $x_{2}$.

Comparing the output of the trained $gg\to\gamma\gamma gg$ model to the
amplitude library value by unweighted point-by-point ratio in Fig. 4, we see a
peaked and approximately symmetric error distribution with a shifted mean in
both regions. Both region histograms have a similar mean, indicating
comparable accuracy, with $\mathcal{R}_{\textrm{non-div}}$ performing slightly
better. The distributions are fairly broad. $\mathcal{R}_{\textrm{non-div}}$
shows a slight tail on the right, which arises from points near the cutoff
$y$.

Despite the per-point agreement being somewhat poor, the total cross section
is found to be in agreement, with $\sigma_{\text{\acs{NN}}}=(45\pm 6)\times
10^{-7}\thinspace\mathrm{pb}$ (precision/optimality error) and
$\sigma_{\mathrm{NJet}}=(49\pm 5)\times 10^{-7}\thinspace\mathrm{pb}$ (MC
error). Figure 4 shows that the regions that are sampled the most due to the
shape of the gluon PDF are those that have the lowest error. In addition, the
PDF-weighted histograms in Fig. 4 are more narrowly peaked and closer to being
unit-centred than the unweighted histograms. Thus, the agreement in the total
cross section is much better than for point-by-point comparison because poorly
performing points fall in PDF-suppressed regions. This indicates that the
accuracy of distributions inferred with the model is dependent on the choice
of process, cuts, and observable.

Figure 5:  Differential normalised cross sections for $gg\to\gamma\gamma gg$,
comparing NJet (MC error) to the model (precision/optimality error), in
R-separation between the hardest jet and photon (left) and the transverse
momentum of the hardest jet (right). Refer to Ref. [12] for definitions of
observables and cuts, and further distributions.

Fig. 5 shows excellent agreement between the distributions obtained from the
model and NJet for two differential slices of phase space. There are some
fluctuations in the tails although they appear statistical rather than
systematic and the model predictions mostly remain within the NJet MC error
bands. Normalised NN uncertainties are negligible compared to the MC error.

In Ref. [12], we also demonstrates how agreement can be improved in
$\mathcal{R}_{\textrm{div}}$ by reweighting event weights by the ratio of the
emulated and true MEs at known points from the training data, as well as
showing good performance for $gg\to\gamma\gamma g$ when relaxing cuts at
inference compared to training.

Subsequent to this work, the authors of Ref. [28] achieve improved per-point
agreement at tree-level by exploiting the factorisation properties of MEs.

##  7 Conclusion

We extend previous work which pioneered the emulation of MEs with NNs,
studying these techniques for the first time within a full hadronic collider
simulation. We focus on loop-induced diphoton-plus-jets production via gluon
fusion. The difficulties introduced by IR behaviour are tamed by partitioning
the phase space as prescribed by FKS subtraction. We provide a general
interface for trained models to Sherpa.

We find that amplitude NN models provide an efficient and general framework
for optimising high-multiplicity observables at hadronic colliders. Agreement
in differential distributions is excellent. As the cost of inference is
negligible compared to the amplitude library call in training, the speed up in
total simulation time (including training) compared to conventional methods is
given by the ratio of the number of points used for inference and training,
$N_{\mathrm{infer}}/N_{\mathrm{train}}$. For this study, this gave a factor of
thirty, although for studies with higher statistics or coverage of multiple
cut configurations, the factor would be much greater.

I would like to thank Joseph Aylett-Bullock and Henry Truong for useful
discussions, and Simon Badger for comments on the draft manuscript. I am
supported by UKRI-STFC ST/S505365/1 and ST/P001246/1.

## References

  * [1] Degrande C, Hirschi V and Mattelaer O 2018 Ann. Rev. Nucl. Part. Sci. 68 291–312
  * [2] Agarwal B, Buccioni F, von Manteuffel A and Tancredi L 2021 Phys. Rev. Lett. 127 262001 (Preprint 2105.04585)
  * [3] Chawdhry H A, Czakon M, Mitov A and Poncelet R 2021 JHEP 09 093 (Preprint 2105.06940)
  * [4] Badger S et al. 2021 JHEP 2021 (Preprint 2106.08664)
  * [5] Badger S, Gehrmann T, Marcoli M and Moodie R 2022 Phys. Lett. B 824 136802 (Preprint 2109.12003)
  * [6] Bargiela P, Caola F, von Manteuffel A and Tancredi L 2022 JHEP 02 153 (Preprint 2111.13595)
  * [7] Feickert M and Nachman B 2021 (Preprint 2102.02770)
  * [8] Badger S and Bullock J 2020 JHEP 06 114 (Preprint 2002.07516)
  * [9] Badger S, Biedermann B, Uwer P and Yundin V 2013 Comput. Phys. Commun. 184 1981–1998 (Preprint 1209.0100)
  * [10] Badger S, Biedermann B, Uwer P and Yundin V 2014 J. Phys. Conf. Ser. 523 012057 (Preprint 1312.7140)
  * [11] Bothmann E et al. (Sherpa) 2019 SciPost Phys. 7 034 (Preprint 1905.09127)
  * [12] Aylett-Bullock J, Badger S and Moodie R 2021 JHEP 2021 (Preprint 2106.09474)
  * [13] Aylett-Bullock J and Moodie R 2021 n3jet_diphoton v1 https://gitlab.com/JosephPB/n3jet_diphoton
  * [14] de Florian D and Kunszt Z 1999 Phys. Lett. B 460 184–188 (Preprint hep-ph/9905283)
  * [15] Badger S D 2009 JHEP 01 049 (Preprint 0806.4600)
  * [16] Ossola G, Papadopoulos C G and Pittau R 2007 Nucl. Phys. B 763 147–169 (Preprint hep-ph/0609007)
  * [17] Peraro T 2019 JHEP 07 031 (Preprint 1905.08019)
  * [18] Frederix R, Frixione S, Maltoni F and Stelzer T 2009 JHEP 10 003 (Preprint 0908.4272)
  * [19] Chollet F et al. 2015 Keras https://github.com/fchollet/keras
  * [20] Abadi M et al. 2015 TensorFlow https://www.tensorflow.org/
  * [21] Kingma D P and Ba J 2015 3rd International Conference for Learning Representations (Preprint 1412.6980)
  * [22] Guennebaud G, Jacob B et al. 2010 Eigen v3 https://eigen.tuxfamily.org
  * [23] Buckley A et al. 2015 Eur. Phys. J. C 75 132 (Preprint 1412.7420)
  * [24] Ball R D et al. (NNPDF) 2017 Eur. Phys. J. C 77 663 (Preprint 1706.00428)
  * [25] Badger S, Guffanti A and Yundin V 2014 JHEP 03 122 (Preprint 1312.5927)
  * [26] Bierlich C et al. 2020 SciPost Phys. 8 026 (Preprint 1912.05451)
  * [27] Aaboud M et al. (ATLAS) 2017 Phys. Rev. D 95 112005 (Preprint 1704.03839)
  * [28] Maître D and Truong H 2021 JHEP 11 066 (Preprint 2107.06625)

  *[SM]: Standard Model
  *[NLO]: next-to-leading order
  *[NNLO]: next-to-next-to-leading order
  *[QCD]: quantum chromodynamics
  *[NN]: neural network
  *[MEs]: matrix element
  *[MC]: Monte Carlo
  *[PDF]: parton distribution function
  *[IR]: infrared
  *[NNs]: neural network
  *[LO]: leading order
  *[ME]: matrix element
  *[ML]: machine learning


# Super-resolving 2D stress tensor field conserving equilibrium constraints
using physics informed U-Net

Kazuo Yonekura
Department of Systems Innovation
The University of Tokyo
Tokyo, JAPAN 113-8656
<EMAIL_ADDRESS>
Kento Maruoka
Department of Systems Innovation
The University of Tokyo
Tokyo, JAPAN 113-8656
Kyoku Tyou
Department of Systems Innovation
The University of Tokyo
Tokyo, JAPAN 113-8656
Katsuyuki Suzuki
Department of Systems Innovation
The University of Tokyo
Tokyo, JAPAN 113-8656

###### Abstract

In a finite element analysis, using a large number of grids is important to
obtain accurate results, but is a resource-consuming task. Aiming to real-time
simulation and optimization, it is desired to obtain fine grid analysis
results within a limited resource. This paper proposes a super-resolution
method that predicts a stress tensor field in a high-resolution from low-
resolution contour plots by utilizing a U-Net-based neural network which is
called PI-UNet. In addition, the proposed model minimizes the residual of the
equilibrium constraints so that it outputs a physically reasonable solution.
The proposed network is trained with FEM results of simple shapes, and is
validated with a complicated realistic shape to evaluate generalization
capability. Although ESRGAN is a standard model for image super-resolution,
the proposed U-Net based model outperforms ESRGAN model in the stress tensor
prediction task.

_Keywords_ Super resolution $\cdot$ Deep neural networks $\cdot$ Physics
guided machine learning $\cdot$ U-Net

##  1 Introduction

With huge advances in computational mechanics, finite element analysis is
widely used for designing and optimizing structures. Recently, to obtain
accurate computation results, extensive grids are being used for finite
element method (FEM) analysis [10, 13, 20] as well as structural optimization
[1, 12] because the accuracy of FEM is better with fine mesh than with coarse
mesh. However, the larger the FEM model becomes, the more computation time and
memory are required.

Reduction of computation time in FEM is desired from the points of view of
several applications. Real-time FEM simulation has recently become available
owing to advances in computation resources and software [11]. Additionally,
FEM-based real-time simulation software is commercially available [2].
However, such applications do not allow an extensive number of grids because
only a little computation time is allowed. Hence, the accuracy of the FEM
analysis is relatively low. Therefore, it is desirable to obtain accurate
results with coarse mesh .

Short computation time is also desired in shape and topology optimization. In
this regard, the FEM analysis is iteratively conducted to calculate the
sensitivity from the FEM results. FEM computation time directly affects the
total computation time; hence, short computation time is desired. Recently,
machine learning models have been utilized to realize a shorter FEM
computation time ; [16] proposed a surrogate model-based method to approximate
the FEM results.

In computational mechanics research, the superconvergent patch recovery (SPR)
method [8, 23] is proposed to obtain accurate results from coarse mesh. The
SPR method uses a shape function to interpolate stress fields. The accuracy of
the SPR method is guaranteed, but the computation time is not small because it
needs matrix products of shape functions.

Additionally, a deep neural network (DNN) model is also utilized to realize
the aforementioned aim. The advantage of DNN is that the computation time for
inference is short, e.g. less than a few seconds. However, the accuracy of the
result is not guaranteed. One approach using DNN directly predicts contour
plots from boundary conditions. [5] proposed a StressGAN that predicts two-
dimensional stress contour. Another approach is to conduct a super-resolution
of the contour plots from coarse FEM results. [16] proposed a convolution
neural network model to predict fine-mesh results from coarse-mesh results.
However, the training and test datasets are similar and generalization ability
is not well studied. In those existing literature models output only scalar
fields, such as von Mises stress or strain fields. However, in FEM analysis,
it is required to predict stress tensor components. In addition, the results
should satisfy the physical equations, which is not considered in the existing
models do not satisfy them. In particular, the generated contour image does
not satisfy equilibrium equations. In this study, a DNN-based super-resolution
model for a two-dimensional stress tensor field considering equilibrium
constraints is proposed.

Usually, in a super-resolution task, the pixel resolution is improved.
However, in the FEM super-resolution task, the pixel resolution need not need
be improved. The contour image is shown on a FEM mesh, and the resolution of
the contour should be improved.

In image processing research, several super-resolution models have been
proposed. The most successful models are SRGAN [9] and ESRGAN [17], which are
based on the generative adversarial network (GAN) [3] framework. One of the
issues in SRGAN and ESRGAN is that training the models is generally unstable
[18]. Therefore, in this study, a non-GAN-based model is utilized to avoid
training instability. The U-Net [14] is a fully convolution network model, and
U-Net++ [22] is a modified model of U-Net that outperforms U-Net. The U-Net
was first proposed for semantic segmentation of medical images. However, in
this study, it is utilized for the super-resolution task. The U-Net-based
models are non-GAN-based and can be trained more stably than GAN-based models.
In the semantic segmentation task, each pixel is classified for its
corresponding classes, e.g., sky, road, or human. Hence, the output is a
binary value on each pixel. By modifying the U-Net, the model outputs
continuous values corresponding to stress tensor components. Thus, it the
model can be utilized for super-resolution of stress tensor components. In
this study, U-Net based models are proposed and compared with ESRGAN-based
models.

The super-resolved images should satisfy the equilibrium equations. To obtain
such results, the physical constraints are considered using the framework of a
physics-informed neural network (PINN) [7, 19]. One of the methods to consider
physical equations is to add the residual of the physical equation to a loss
function. For example, [21] used the residual of Newton’s second law to
predict seismic response, and [4] used the residual of energy conservation to
predict lake temperature. The stress tensor follows equilibrium equations. In
this study, the residual of the equilibrium equations is employed as the
physical loss function. By employing the physical loss, the obtained super-
resolved image satisfies the equilibrium constraint.

This study is organized as follows. The overview of the proposed model and
dataset is explained in section 2. The proposed physics-informed models are
introduced in section 3. Numerical experiments are presented in section 4.
Finally, conclusions follow in section 5.

##  2 Proposed model overview and dataset

###  2.1 Problem definition and model overview

Herein, the task is to obtain reasonable smooth contour images from coarse FEM
results. The super-resolution model inputs stress tensor data from a coarse
mesh and outputs super-resolved data. In this study, only two-dimensional
images were used. The stress tensor $\sigma$ in the two-dimensional analysis
is written as

$\displaystyle\sigma=\begin{pmatrix}\sigma_{x}&\tau_{xy}\\\
\tau_{xy}&\sigma_{y}\end{pmatrix}.$ (1)

An overview of the proposed method is shown in Figure 1. The neural network
(NN) model inputs coarse mesh images of three stress tensor components and
outputs super-resolved images. The loss function is calculated between super-
resolved images and ground truth images. To handle non-uniform shapes, such as
rectangles and L-shapes, data are given as pixel images. The input coarse
images and ground truth fine images are obtained using the FEM analysis. It is
noted that even in the fine mesh result images, the mesh size is generally
larger than the pixel size. It is desirable if the mesh resolution of the
output super-resolved images is finer than that of ground truth.

Figure 1: Overview of the proposed model.

###  2.2 Dataset preparation

FEM analysis was conducted using commercial FEM software, and gray-scale
contour images are obtained. The contour range was defined when generating the
contour images. The intensity of pixels was calculated by the post-processing
module based on the stress components and contour map. We used the linear
contour map to calculate stress components from the intensity of images; the
stress components were calculated by multiplying a coefficient $C_{k}$ with
intensity $I^{1,2,3}_{k}(i,j)$ at each pixel $(i,j)$;

$\displaystyle\sigma_{x}(i,j)|_{k}=C_{k}\left(I^{1}_{k}(i,j)+s_{k}\right),$
(2)
$\displaystyle\sigma_{y}(i,j)|_{k}=C_{k}\left(I^{2}_{k}(i,j)+s_{k}\right),$
(3) $\displaystyle\tau_{xy}(i,j)|_{k}=C_{k}\left(I^{3}_{k}(i,j)+s_{k}\right),$
(4)

where $i$ and $j$ are the locations of the pixel; $k$ is the suffix of the
image; and $s_{k}$ is a constant for adjusting contour level. The superscripts
$1$, $2$, and $3$ on $I$ correspond to $\sigma_{x}$, $\sigma_{y}$, and
$\tau_{xy}$, respectively. Consequently, $I^{1}_{k}(i,j)$ is the intensity at
the pixel $(i,j)$ of $k$-th image corresponding to $\sigma_{x}$. $C_{k}$ and
$s_{k}$ are defined for each analysis case $k$, which implies that the contour
maps are the same for $\sigma_{x}$, $\sigma_{y}$, and $\tau_{xy}$ of one
analysis, but different between different analysis cases. The data consists of
a set of coarse mesh results and fine mesh results, i.e., three coarse contour
images and three fine contour images.

The dataset comprises training and validation datasets. From an application
point of view, it is important that a model trained with a small number of
data is valid for a wide range of analyses. Hence, in this study, we trained
the model using relatively simple data and validated it with a complicated
truss-like structure. The training dataset was composed of simple analysis
cases, i.e., cantilevers and L-shaped beams, whereas the validation data was a
truss-like cantilever; both datasets are shown in Figure 2. Linear
quadrilateral element was used. The mesh size $\ell$ was $\ell=H/10$ and
$\ell=d/5$ in the coarse mesh analysis, and $\ell=H/40$ and $\ell=d/20$ in the
fine mesh analysis; essentially, the fine mesh was $4$ times smaller than that
of the coarse mesh. Data were generated using the same shape on different
boundary conditions and loads. In the cantilever and L-shape, the constraint
condition was selected from two types, i.e., fixed boundary and sliding
boundary. The loads were selected from multiple cases listed in Table 1. The
magnitude of total loads was $1kN$ in all cases. The dimensions were set as
$W=2H$ and $L=3d$. After preparing the FEM results, data augmentation was
performed to increase the number of data. The image was inverted vertically
and horizontally. The contour image of the inside shape, i.e., only the
contour area excluding the white background area, was also inverted.
Essentially, black was inverted into white. The total number of datasets for
cantilever and L-shape was $504$, including augmented data. The dataset was
randomly split into $3:1$ for training and test data. A truss-like cantilever
is employed as the validation dataset, which has a similar shape to the
topology-optimized shape, to determine if the model could be applied to
topology optimization. The element type in the validation dataset was triangle
to determine if the model could work with different types of elements. For
data augmentation, left-right, top-bottom, and black-white flips were enabled.

Table 1: Loading conditions of FEM models

Cantilever

Load type | Direction | $\bar{y}$
---|---|---
Concentrated | $x$ | $id/5$
Concentrated | $y$ | $id/5$
Distributed | $x$ | —
Distributed | $y$ | —

($i=0,1,\dots,5$)

L-shape

Load type | Direction | $\bar{y}$
---|---|---
Concentrated | $x$ | $iH/10$
Concentrated | $y$ | $iH/10$
Distributed | $x$ | —
Distributed | $y$ | —

($i=0,1,\dots,5$)

(a) Cantilever data.

(b) L-shape data.

(c) Truss-like cantilever.

Figure 2: Close up view of coarse, fine, and super-resolved images.

##  3 Physics informed super-resolution models

###  3.1 Physical loss

The stress tensor $\sigma$ has to satisfy the equilibrium constraint;

$\displaystyle{\rm
div}\sigma+\boldsymbol{f}=\begin{pmatrix}\frac{\partial\sigma_{x}}{\partial
x}+\frac{\partial\tau_{xy}}{\partial y}+f_{x}\\\
\frac{\partial\sigma_{y}}{\partial y}+\frac{\partial\tau_{xy}}{\partial
x}+f_{y}\end{pmatrix}=\boldsymbol{0},$ (5)

where $\boldsymbol{f}=(f_{x},f_{y})^{\top}$ implies a body force. The physics
loss $\mathcal{L}_{\rm PL}$ is defined as a penalty for violating equation
(5).

$\displaystyle\mathcal{L}_{\rm PL}=\|{\rm div}\sigma+\boldsymbol{f}\|^{2}.$
(6)

To calculate derivatives in equation 5 from output data
$\boldsymbol{y}=(y_{1},y_{2},y_{3})^{\top}$, the central difference is applied
on pixels. The derivatives at pixel $(i,j)$ are calculated by

$\displaystyle\frac{\partial\sigma_{x}}{\partial
x}(i,j)=y_{1}(i+1,j)-y_{1}(i-1,j),$ (7)
$\displaystyle\frac{\partial\sigma_{y}}{\partial
y}(i,j)=y_{2}(i,j+1)-y_{2}(i,j-1),$ (8)
$\displaystyle\frac{\partial\tau_{xy}}{\partial
x}(i,j)=y_{3}(i+1,j)-y_{3}(i-1,j),$ (9)
$\displaystyle\frac{\partial\tau_{xy}}{\partial
y}(i,j)=y_{3}(i,j+1)-y_{3}(i,j-1).$ (10)

The loss $\mathcal{L}_{\rm PL}$ must be summed only the internal region of the
shapes. Hence, the mask function is defined to extract the internal region.

$\displaystyle m(i,j)=1{\rm~{}if~{}}y_{1}(i+1)>1-\varepsilon$ (11)

and $\hat{\mathcal{L}}_{\rm PL}=\sum_{i,j}m(i,j)\mathcal{L}_{\rm PL}(i,j)$. In
the dataset, the external loads $\boldsymbol{f}$ are applied only on the
boundary edges and points. Those edges and points are masked by the mask
function. Hence, $\boldsymbol{f}$ in (6) can be eliminated.

It is obvious that $\hat{\mathcal{L}}_{\rm PL}$ is minimum when
$\boldsymbol{y}$ is constant everywhere. However, in such a case, the MSE loss
becomes large.

###  3.2 Modified U-Net/U-Net++ and PI-UNet/PI-UNet++

U-Net [14] is an encoder-decoder model that is proposed for semantic
segmentation for medical images, which is not usually used for super-
resolution. In this study, the U-Net is slightly modified to allow its
applicability for super-resolution purpose. The architecture of U-Net is shown
in Figure 3(a). The input data is encoded to a low-level feature map using the
max-pooling layer and the feature map is decoded to the original resolution.
The feature map of each encoder level is connected to the feature map of the
same decoder level.

(a) U-Net architecture. (Cited from [14] Fig.1)

(b) U-Net++ architecture. (Cited from [22] Fig.1 (a))

Figure 3: Architecture of U-Net and U-Net++.

The U-Net NN inputs images and outputs the segmentation map. The segmentation
map consists of binary data on each pixel, i.e., true/false with respect to
the corresponding classes. The softmax cross-entropy is employed as the loss
function in the U-Net.

In this study, the aim is to predict the intensity of each pixel whose value
is in $[0,1]$ from original images. Hence, we removed the softmax function and
used the output from the sigmoid function at the last layer so that the model
outputs the value in $[0,1]$. In addition, the loss function is modified and
the mean squared error (MSE) is used;

$\displaystyle\mathcal{L}_{\rm
MSE}=\frac{1}{N}\left\|\boldsymbol{t}-\boldsymbol{o}\right\|^{2},$ (12)

where $\boldsymbol{t}$ and $\boldsymbol{o}$ are the ground truth and output
data, respectively; and $N$ is the number of pixels. In this study, the
modified model is called a modified U-Net. It is noted that the resolution of
both the input and output images is the same, which is a common situation in
U-Net.

In the physics-informed U-Net (PI-UNet), the architecture of the model is the
same as the modified U-Net, except that physical loss $\hat{\mathcal{L}}_{\rm
PL}$ is added to the loss function of the modified U-Net.

$\displaystyle\mathcal{L}_{\rm PI\mathchar 45\relax UNet}=\mathcal{L}_{\rm
MSE}+\hat{\mathcal{L}}_{\rm PL}.$

The pixel resolution of the output image is the same as that of the input
image, which is a normal condition in UNet++.

UNet++ [22] is a nested U-Net model and is shown in Figure 3(b). The UNet++
adds a decoder from every encoder level to the input level. Herein, the
U-Net++ outperformed U-Net in the medical image segmentation task.

###  3.3 ESRGAN and PI-ESRGAN

ESRGAN [17] is a GAN-based super-resolution model. In the GAN framework,
networks consist of a generator network and a discriminator network. The loss
function in the GAN is defined as follows.

$\displaystyle V(D,G)=E_{\boldsymbol{x}\sim p_{\rm data}}\left[\log
D(\boldsymbol{x})\right]+E_{\boldsymbol{z}\sim p_{z}}\left[1-\log
D(G(\boldsymbol{z}))\right]$ (13)

where G and D imply the generator and the discriminator, respectively. $G$
inputs noise $\boldsymbol{z}$. $\boldsymbol{z}$ is a noise vector,
$\boldsymbol{x}$ is a real data, $p_{\rm data}$ and $p_{z}$ are the
probability distribution of the real data and noise, respectively. The
generator minimizes $V(D,G)$, whereas the discriminator maximizes V.

In the ESRGAN model, the relativistic discriminator, proposed in the
relativistic GAN [6], is employed. The ordinal discriminator is described by
$D(\boldsymbol{x})={\rm sigmoid}\left(C\left(\boldsymbol{x}\right)\right)$,
where the function $C$ implies the output from the discriminator without the
activation function. Then, the relativistic discriminator, $D_{R_{a}}$, is
defined as

$\displaystyle
D_{R_{a}}\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)={\rm
sigmoid}\left(C(\boldsymbol{x})-E_{\boldsymbol{x}^{\prime}}\left[C\left(\boldsymbol{x}^{\prime}\right)\right]\right).$
(14)

Thus, the loss function of the discriminator network, $L^{\rm D}_{\rm
ESRGAN}$, is formulated using the relativistic discriminator as

$\displaystyle L^{\rm D}_{\rm ESRGAN}=-E_{\boldsymbol{x}\sim
p_{data}}\left[\log\left(D_{R_{a}}\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)\right)\right]-E_{\boldsymbol{x}^{\prime}\sim
p_{x^{\prime}}}\left[\log\left(1-D_{R_{a}}\left(\boldsymbol{x}^{\prime},\boldsymbol{x}\right)\right)\right].$
(15)

The loss function of the generator network in ESRGAN comprises the perceptual
loss, GAN loss, and $L_{1}$ loss;

$\displaystyle L^{\rm G}_{\rm ESRGAN}=L_{percep}+\lambda L^{R_{a}}_{G}+\eta
L_{1},$ (16)

where $\lambda$ and $\eta$ are the hyperparameters to balance each term;
$\lambda=0.005$ and $\eta=0.01$ are used in the numerical examples. The
perceptual loss is

$\displaystyle
L_{percep}=\frac{1}{W_{i,j}H_{i,j}}\sum_{x=1}^{W_{i,j}}\sum_{y=1}^{H_{i,j}}\left(\phi_{i,j}\left(I^{HR}\right)_{x,y}-\phi_{i,j}\left(G\left(I^{LR}\right)\right)_{x,y}\right)^{2},$
(17)

where $I^{HR}$ and $I^{LR}$ imply the pixel intensity of the high and low-
resolution image, respectively. $\phi_{i,j}$ represents the feature map
obtained by the $j$-th convolution before $i$-th max-pooling layer in the pre-
trained VGG19 network [15]. $H_{i,j}$ and $W_{i,j}$ are the dimensions of the
corresponding feature maps within the VGG19 network. The GAN loss,
$L^{R_{a}}_{G}$, is formulated as

$\displaystyle L^{R_{a}}_{G}=-E_{\boldsymbol{x}\sim
p_{data}}\left[\log\left(1-D_{R_{a}}\left(\boldsymbol{x},\boldsymbol{x}^{\prime}\right)\right)\right]-E_{\boldsymbol{x}^{\prime}\sim
p_{x^{\prime}}}\left[\log\left(D_{R_{a}}\left(\boldsymbol{x}^{\prime},\boldsymbol{x}\right)\right)\right],$
(18)

which has a symmetric form of discriminator loss, $L^{R_{a}}_{D}$, in (15).
The $L_{1}$ loss is

$\displaystyle L_{1}=E_{I^{LR}}\|G\left(I^{LR}\right)-\boldsymbol{t}\|_{1}$
(19)

In the original ESRGAN model, the resolution of the high-resolution training
data is different from that of the low-resolution data. However, as previously
discussed, both resolutions are set the same in this study.

In the physics informed ESRGAN (PI-ESRGAN), the loss function of the generator
is modified to include physics loss;

$\displaystyle\mathcal{L}_{\rm PI\mathchar 45\relax ESRGAN}^{\rm G}=L^{\rm
G}_{\rm ESRGAN}+\hat{\mathcal{L}}_{\rm PL}.$

The loss function of the discriminator in PI-ESRGAN is the same as in original
ESRGAN; $\mathcal{L}_{\rm PI\mathchar 45\relax ESRGAN}^{\rm D}=L^{\rm D}_{\rm
ESRGAN}$.

##  4 Numerical experiments

###  4.1 Training and testing on the similar data

The six models proposed in section 3 are trained using the same training
dataset. The results of cantilever and L-shape examples are shown in Figure
4.1 and Figure 4.1, respectively. The super-resolved images of ESRGAN and PI-
ESRGAN in Figure 4.1 and Figure 4.1, respectively, show that both models
failed to super-resolve images.

The loss values of U-Net-based models are listed in Table 2, where in the
U-Net and UNet++ models, physical loss is not included; however, it is
calculated to compare the results. U-Net exhibits the smallest MSE loss among
all models, whereas PI-UNet exhibits the smallest physical loss. The PI-UNet
reduced the physical loss at the cost of MSE loss, which led to a lower total
loss. The loss values of UNet++ and PI-UNet++ were worse than those of UNet
and PI-UNet, respectively. The network architecture of UNet++ is complicated
compared with UNet. Hence, the number of network parameters to be trained was
larger in UNet++ than in UNet. This difference causes a slight difference in
the loss values.

The close-up view of a PI-UNet result is shown in Figure 6. The mesh shape can
be seen in the coarse and fine mesh, but the image becomes smooth in the UNet
and PI-UNet results.

The difference between ESRGAN-based and UNet-based models arises from the
difference in loss functions; the UNet-based models contain MSE loss, whereas
the ESRGAN-based models do not. The ESRGAN-based models consider perceptual
loss, and the generator is trained in an adversarial manner with the
discriminator. Hence, the discriminator has to be well trained to obtain a
good generator. The generator does not know the ground truth; hence it is a
difficult task for the generator to learn. For example, the super-resolved
images of $\sigma_{x}$ in ESRGAN and $\sigma_{y}$ in PI-ESRGAN in Figure 4.1
are completely different from ground truth images. The generators output such
wrong images because they do not know the ground truth. This learning
instability is often observed in the GAN literature. In contrast, the loss
functions of UNet-based models contain MSE loss between output and ground
truth. Hence, the models can output at least similar data to the ground truth.
The super-resolved images in Figure 4.1 and Figure 4.1 qualitatively look the
same as ground truth images. Consequently, its loss function and simple
architecture are the reasons for the superiority of U-Net and PI-UNet.

Table 2: Loss values of UNet-based models | Training loss | Validation loss
---|---|---
| Total loss | MSE loss | physical loss | Total loss | MSE loss | physical loss
UNet | 14.58$\times 10^{-4}$ | 1.18$\times 10^{-4}$ | 13.40$\times 10^{-4}$ | 19.83$\times 10^{-4}$ | 6.53$\times 10^{-4}$ | 13.28$\times 10^{-4}$
PI-UNet | 8.63$\times 10^{-4}$ | 4.44$\times 10^{-4}$ | 4.19$\times 10^{-4}$ | 13.80$\times 10^{-4}$ | 9.58$\times 10^{-4}$ | 4.22$\times 10^{-4}$
UNet++ | 22.90$\times 10^{-4}$ | 10.45$\times 10^{-4}$ | 12.45$\times 10^{-4}$ | 21.79$\times 10^{-4}$ | 8.99$\times 10^{-4}$ | 12.80$\times 10^{-4}$
PI-UNet++ | 18.92$\times 10^{-4}$ | 14.55$\times 10^{-4}$ | 4.37$\times 10^{-4}$ | 19.57$\times 10^{-4}$ | 15.17$\times 10^{-4}$ | 4.40$\times 10^{-4}$
Model |  |  |
---|---|---|---
Figure 4: Cantilever example

Model |  |  |
---|---|---|---
Figure 5: L-shape example

(a) Coarse.

(b) Fine.

(c) U-Net.

(d) PI-UNet.

Figure 6: Close up view of coarse, fine, and super-resolved images.

###  4.2 Generalization ability to unseen tasks

The trained model is validated using the truss-like cantilever. The losses are
shown in Table 3, and the super-resolved images are shown in Figure 4.2. All
the super-resolved images reproduce the shape of the original truss-like
cantilever, although the models are not trained with the data. The model is
trained to draw contours only in the contour area of the input coarse mesh.
However, ESRGAN-based models failed to predict accurately. The MSE and
physical losses of both models are greater than those of coarse mesh.
Contrastingly , UNet-based models succeeded in super-resolving data. In
particular, the physical loss of PI-UNet is small. The PI-UNet can predict
unfamiliar data from these experiments, although similar data is excluded from
the training data. However, the prediction accuracy is worse than the
experiments in section 4.1. It is reasonable because, in section 4.1, the test
data and training data are similar. Consequently, the model can be applied to
unseen data; moreover, to improve accuracy, similar data should be added to
the training dataset.

Model |  |  |
---|---|---|---
Figure 7: Truss-like cantilever example Table 3: Loss values of an unseen task. | Total loss | MSE loss | physical loss
---|---|---|---
Coarse | 81.28$\times 10^{-4}$ | 76.74$\times 10^{-4}$ | 4.54$\times 10^{-4}$
ESRGAN | 133.09$\times 10^{-4}$ | 128.20$\times 10^{-4}$ | 4.89$\times 10^{-4}$
PI-ESRGAN | 117.30$\times 10^{-4}$ | 111.51$\times 10^{-4}$ | 5.79$\times 10^{-4}$
U-Net | 61.98$\times 10^{-4}$ | 57.75$\times 10^{-4}$ | 4.23$\times 10^{-4}$
PI-UNet | 60.14$\times 10^{-4}$ | 58.38$\times 10^{-4}$ | 1.76$\times 10^{-4}$
UNet++ | 62.52$\times 10^{-4}$ | 57.72$\times 10^{-4}$ | 4.81$\times 10^{-4}$
PI-UNet++ | 83.50$\times 10^{-4}$ | 81.17$\times 10^{-4}$ | 2.33$\times 10^{-4}$

##  5 Conclusion

In this study, the UNet-based super-resolution model was proposed for
predicting the stress tensor components from coarse mesh results. Because the
equilibrium equation had to be satisfied in the FEM solution, physical loss
was defined using divergence of the stress tensor. The PI-UNet model indicated
the best score while considering the MSE and physics loss. Although ESRGAN is
an outstanding model for super-resolution in image processing literature, the
UNet-based models outperformed the ESRGAN-based models for the present task.
The proposed model was validated using two types of data: test and validation
data. The test data was similar to the training data. The model showed good
agreement with the ground truth for the test data. Moreover, although the mesh
shape was visible in ground truth data, it was not apparent in the super-
resolved images, which is desirable from an application point of view.
Contrastingly, the validation data is completely different unseen shapes. The
model outputs a reasonable solution for such data. In particular, the physical
loss is as small as that of the test data.

The machine learning model should be trained with a sufficient amount of data.
However, in the FEM analysis application, it is not possible to collect all
the possible data. This study shows that once trained, the proposed model is
valid for data similar to the training data and outputs reasonable solutions
for unseen tasks.

## Acknowledgement

This work was partially supported by JSPS KAKENHI Grant number 21K14064.

## References

  * [1] Niels Aage, Erik Andreassen, Boyan Stefanov Lazarov, and Ole Sigmund.  Giga-voxel computational morphogenesis for structural design.  Nature, 550:84–86, 2017.
  * [2] Christopher Fleischmann, Irina Leher, Reinhold Hartwich, Marc Hainke, and Stefan Sesselmann.  A new approach to quickly edit geometries and estimate stresses and displacements of implants in real-time.  Current Directions in Biomedical Engineering, 5:553–556, 2019.
  * [3] Ian J. Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In Proceedings of the 27th International Conference on Neural Information Processing Systems - Volume 2, NIPS’14, pages 2672–2680, Cambridge, MA, USA, 2014. MIT Press.
  * [4] Xiaowei Jia, Jared Willard, Anuj Karpatne, Jordan S. Read, Jacob A. Zwart, Michael Steinbach, and Vipin Kumar.  Physics-guided machine learning for scientific discovery: An application in simulating lake temperature profiles.  ACM/IMS Transactions on Data Science, 2(3), may 2021.
  * [5] Haoliang Jiang, Zhenguo Nie, Roselyn Yeo, Amir Barati Farimani, and Levent Burak Kara.  StressGAN: A generative deep learning model for 2D stress distribution prediction.  Journal of Applied Mechanics, 88:1–11, 2021.
  * [6] Alexia Jolicoeur-Martineau.  The relativistic discriminator: a key element missing from standard GAN, 2019.
  * [7] George Em Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3:422–440, 2021.
  * [8] Mukesh Kumar, Trond Kvamsdal, and Kjetil André Johannessen.  Superconvergent patch recovery and a posteriori error estimation technique in adaptive isogeometric analysis.  Computer Methods in Applied Mechanics and Engineering, 316:1086–1156, 2017.
  * [9] Christian Ledig, Lucas Theis, Ferenc Huszar, Jose Caballero, Andrew Cunningham, Alejandro Acosta, Andrew Aitken, Alykhan Tejani, Johannes Totz, Zehan Wang, and Wenzhe Shi.  Photo-realistic single image super-resolution using a generative adversarial network.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), July 2017.
  * [10] Hongshuai Lei, Chuanlei Li, Jinxin Meng, Hao Zhou, Yabo Liu, Xiaoyu Zhang, Panding Wang, and Daining Fang.  Evaluation of compressive properties of SLM-fabricated multi-layer lattice structures by experimental test and $\mu$-CT-based finite element analysis.  Materials & Design, 169:107685, 2019.
  * [11] Dragan Marinkovic and Manfred Zehn.  Survey of finite element method-based real-time simulations.  Applied sciences, 9:2775, 2019.
  * [12] Sougata Mukherjee, Dongcheng Lu, Balaji Raghavan, Piotr Breitkopf, Subhrajit Dutta, Manyu Xiao, and Weihong Zhang.  Accelerating large-scale topology optimization: State-of-the-art and challenges.  Archives of Computational Methods in Engineering, 28:4549–4571, 2021.
  * [13] Shinji Nakazawa, Naoki Iwasaki, Ryuichi Matsuki, Gaku Hashimoto, Hiroshi Okuda, and Kazuya Goto.  High-resolution structural analysis of multilayer package substrate with open-source parallel FE software FrontISTR.  In 2019 IEEE CPMT Symposium Japan (ICSJ), pages 95–98, 2019.
  * [14] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-Net: Convolutional networks for biomedical image segmentation.  In Nassir Navab, Joachim Hornegger, William M. Wells, and Alejandro F. Frangi, editors, Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, pages 234–241, Cham, 2015. Springer International Publishing.
  * [15] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  In Yoshua Bengio and Yann LeCun, editors, 3rd International Conference on Learning Representations, ICLR 2015, 2015.
  * [16] Ren Kai Tan, Chao Qian, Michael Wang, and Wenjing Ye.  An efficient data generation method for ANN-based surrogate models.  Structural and Multidisciplinary Optimization, 65:90, 2022.
  * [17] Xintao Wang, Ke Yu, Shixiang Wu, Jinjin Gu, Yihao Liu, Chao Dong, Yu Qiao, and Chen Change Loy.  ESRGAN: Enhanced super-resolution generative adversarial networks.  In Proceedings of the European Conference on Computer Vision (ECCV) Workshops, September 2018.
  * [18] Zhang Wenlong, Liu Yihao, Chao Dong, and Yu Qiao.  RankSRGAN: Generative adversarial networks with ranker for image super-resolution.  IEEE Transactions on Pattern Analysis and Machine Intelligence, pages 1–1, 2021.
  * [19] Jared Willard, Xiaowei Jia, Shaoming Xu, Michael Steinbach, and Vipin Kumar.  Integrating physics-based modeling with machine learning: A survey, 2020\.
  * [20] Praveen Yadav and Krishnan Suresh.  Large scale finite element analysis via assembly-free deflated conjugate gradient.  Journal of Computing and Information Science in Engineering, 14(4):041008, 2014.
  * [21] Ruiyang Zhang, Yang Liu, and Hao Sun.  Physics-guided convolutional neural network (phycnn) for data-driven seismic response modeling.  Engineering Structures, 215:110704, 2020.
  * [22] Zongwei Zhou, Md Mahfuzur Rahman Siddiquee, Nima Tajbakhsh, and Jianming Liang.  UNet++: A nested U-Net architecture for medical image segmentation.  In Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support, DLMIA 2018, ML-CDS 2018. Lecture Notes in Computer Science, Vol. 11045, pages 3–11. Springer, Cham, 2018.
  * [23] O. C. Zienkiewicz and J. Z. Zhu.  The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique.  International Journal for Numerical Methods in Engineering, 33(7):1331–1364, 1992.


# Measurement of low energy nuclear recoil events with the Phonon-mediated
voltage-assisted Hybrid Detector for rare event searches

S. Maludze  W. Baker  M. Mirzakhani  M. Lee  C. Savage  H. Neog  R. Mahapatra
N. Mirabolfathi  M. Platt  A. Jastram  Department of Physics and Astronomy,
Texas A&M University, College Station, TX, 77843

###### Abstract

The phonon-mediated hybrid detector is made out of a monolithic silicon
crystal characterized by two interconnected regions linked through a narrow
neck. Operating solely on phonon signal measurements, the hybrid design
facilitates the differentiation between electron recoil and nuclear recoil
events, effectively discerning two types of interaction down to low energy
levels. With a newly implemented software triggering technique (SWT), low-
energy nuclear recoil events of approximately $500~{}eV_{ee}$ have been
measured. In addition to this, an electron recoil background reduction of up
to $95\%$ has been successfully demonstrated.

##  I Introduction

Numerous astronomical observations have presented evidence suggesting that the
predominant constituent of the Universe’s matter content is non-luminous and
non-baryonic in nature, commonly referred to as dark matter [1, 2]. At the
same time, Only $~{}5\%$ of the universe is composed of observable matter.
Among the leading contenders for the identity of dark matter is the Weakly
Interacting Massive Particle (WIMP), characterized by its interaction on the
weak scale [3, 4]. Various experiments [5, 6, 7, 8] employ direct detection
methodologies, wherein the interaction between dark matter and normal matter
is probed through sensitive detector technologies and background mitigation
techniques. To probe the parameter space at lower mass levels than the current
detection thresholds allow, it becomes imperative to possess detectors
featuring lower thresholds and substantial target masses.

In the conventional configuration of a Cryogenic Dark Matter Search (CDMS)
apparatus, the iZIP (interleaved Z-sensitive Ionization and Phonon) detector
is equipped with both phonon sensors and charge electrodes, enabling the
simultaneous measurement of phonon and ionization signals [9]. This two-signal
detection serves as a tool for identifying a type of particle interaction
within the detector, specifically distinguishing electron recoil (ER) and
nuclear recoil (NR) events. Given the greater ionization yield for ER events,
a fraction of two signals is used as a discriminator factor. In contrast, a
high-voltage cylindrical detector, characterized by a lower threshold compared
to iZIP, is equipped solely with phonon sensors and lacks this discrimination
capability [10].

A novel hybrid detector has been developed with the motivation to combine the
advantageous features of these two technologies [11]. This innovation employs
a monolithic silicon crystal characterized by two distinct regions: a low-
voltage domain with a substantial volume and a high-voltage domain with a
relatively smaller volume. Those two regions are linked by a narrow segment
facilitating the propagation of charge carriers to the high-voltage region.

Within this hybrid detector, the low-voltage region constitutes the fiducial
volume where particle interactions occur. Phonon sensors, positioned at the
base of the low-voltage region, capture phonon signals generated at the
interaction point. On the other hand, electrons produced during the
interaction are guided toward the high-voltage region due to an applied
electric field. While traversing the high-voltage domain, these electrons
generate more Neganov-Trofimov-Luke (NTL) phonons [12, 13] with a total signal
proportional to the applied voltage and the number of electron-hole pairs
generated as a result of an initial interaction. These NTL phonons are then
collected by sensors located at the top of the high-voltage region. As a
result of this setup, the hybrid detector, equipped solely with phonon
sensors, possesses the remarkable capability to capture and quantify both
phonon (lower region) and ionization (upper region) signals. The fraction of
these signals provides a robust discriminator for distinguishing between ER
and NR events, thereby enhancing the detector’s capacity for background
rejection on an event-by-event basis.

This paper presents the updated results of the hybrid performance with a newly
developed software triggering technique (SWT). With this technique, it has
been possible to demonstrate the hybrid detector’s remarkable capability to
distinguish between ER and NR down to low energies, which represents the
region of interest (ROI) for rare event searches.

##  II Working principle

The hybrid detector consists of a monolithic silicon crystal and has two
regions separated by a narrow neck. The lower side of the crystal has a
truncated conical geometry with dimensions of 1.4 cm in height, 7.6 cm in
bottom diameter, and 4.8 cm in top diameter. The conical shape pushes electric
field lines toward the high-voltage region inside the crystal. The upper side
takes the form of a small cylinder with a diameter of 2.5 cm and a height of 1
cm. A narrow 0.5 mm gap separates these two regions. The detector’s fiducial
mass (mass of the low-voltage side) is $230~{}g$, while the total mass is
$250~{}g$.

On the crystal’s top and bottom surfaces, phonon Transition-Edge-Sensors
(TESs) [14] are photolithographically patterned. The top surface has a single
readout channel, referred to as a top channel, while the bottom surface is
subdivided into four readout channels: three inner channels and one outer
channel, denoted as the bottom channels. To ensure its proper functioning, the
detector is placed in a copper housing that not only shields it from external
background radiation but also plays a pivotal role in electric field shaping
within the crystal (Figure 1). Proper electric field shaping is essential for
the operation of the detector.

Figure 1: Picture of the top side of the Hybrid detector in the experimental
setup.

When an interaction occurs within the low-voltage domain, which constitutes
the fiducial volume, primary phonons, and electron-hole pairs are generated in
close proximity to the interaction point. These primary phonons bounce back
and forth within the low-voltage domain until their eventual absorption by the
bottom channels. Simultaneously, the electron-hole pairs drift along the
electric field lines within the domain, generating NTL phonons in the low-
voltage region, which are also captured by the bottom channels. Notably, the
contribution of NTL phonons to the signal registered on the bottom channels is
relatively lower compared to primary phonons. At the same time, the electric
field lines direct electrons toward the high-voltage domain, where they also
generate NTL phonons. These phonons are collected by the top channel. Thus,
the signal observed on the top channel is directly proportional to the number
of electrons created as a result of the initial interaction, providing an
indirect measure of ionization.

A notable feature of the hybrid detector is its capability to estimate both
the ionization yield (the Lindhard factor for nuclear recoil events [15]) and
the recoil energy based on signals from the top and bottom channels,
considering that the total phonon signal has a contribution from both sides of
the crystal. Leakage factor parameters that account for this contribution, are
determined experimentally from the calibration dataset (Pair of formulas
presented in [11]).

To summarize the hybrid detector working principle, the signal acquired from
the bottom channels measures the contribution of primary phonons and a
comparably smaller amount of NTL phonons generated in the low-voltage region.
In contrast, the top channel detects ionization by capturing NTL phonons
generated in the high-voltage domain. A simple ratio of the signals from the
top and bottom channels denoted as $Y$, effectively carries information about
the ionization yield.

##  III Software trigger and data analysis pipeline

In this hybrid experiment, we introduced the software triggering technique
(SWT). In previous runs, event triggering relied on manual threshold settings
based on the real-time display of pulses via the digitizer software. However,
this approach came with significant drawbacks. Notably, it necessitated human
judgment to differentiate noise from genuine pulses, potentially leading to
the oversight of low-energy events that exceeded the detector threshold by a
considerable margin. Once data acquisition was completed, the trigger
threshold was fixed and unmodifiable. Under such circumstances, implementing a
new trigger threshold necessitated the collection of an entirely new dataset.
The hardware trigger method served its purpose well in calibration runs —
where calibration lines comfortably exceeded the threshold and pulse-noise
differentiation was evident.

Figure 2: Flaw chart illustrating steps for updated analysis pipeline with SWT
technique.

With the SWT method, the signal from the detector is continuously recorded.
After recording the data in a desired format, individual pulses are found by
the software and saved in a different directory. The first step in the process
is Cross-Correlation - a pre-selected pulse template is cross-correlated to
the entire readout signal. This step significantly reduces the noise level and
facilitates pulse finding. The next step is finding pulses in the data with
the peak finder algorithm. The algorithm selects pulses with a certain
prominence. Although this method has an extra step of pulse finding and takes
more space and computational resources, the trigger threshold becomes a
tunable parameter and could be adjusted as needed once the raw data has been
taken. After individual pulses are identified and saved, the Optimal Filter
(OF) [16, 17] algorithm is implemented for event reconstruction. The algorithm
takes noise and pulse templates generated from the data and makes a fit to the
pulses. The primary output of this step is the amplitude of a pulse which is
then converted into energy units. This is the efficient way to estimate event
energy considering the noise environment.

Figure 3: An example of a trace (readout signal) containing a pulse and its
cross-correlated version.

##  IV Experimental configuration and results

As described above, the hybrid detector requires bias voltage to be applied
across the crystal for its operation. During this experimental run, numerous
datasets were acquired to ascertain the optimal signal-to-noise ratios and
bias voltage settings. For more accurate electron recoil calibration, an
${}^{241}Am$ source [18] was strategically positioned at the midpoint of the
bottom side of the detector. For investigating the detector’s response to
nuclear recoil events, ${}^{252}Cf$ neutron source was placed $2~{}ft$ away
from the detector with an appropriate shielding for gamma radiation. Data was
acquired with a software-triggering scheme to record near-threshold events
more efficiently.

The hybrid detector characterization process follows several steps - first, a
dataset without any applied voltage across the detector is taken. With known
calibration source lines, calibration factors and phonon leakage parameters
are calculated. Subsequently, the detector is biased with its optimal
operation voltage to investigate its response to both electron and nuclear
recoil events for an extended period.

A detector operating under the applied voltage experiences a phenomenon known
as charge build-up, wherein charge carriers become trapped within crystal
impurities, gradually diminishing the overall effective electric field over
time. To alleviate the accumulation of trapped charge carriers and restore the
detector’s operational state, Light-Emitting Diodes (LEDs) are turned on
periodically for a brief moment of time (LED flashing). Following the
activation of an LED, a cooling period is required for the detector to return
to its equilibrium temperature, during which the data acquisition process is
stopped. Repeated steps of LED flashing, stopping, and resuming the data
acquisition process have been automated. In one cycle, after $20~{}mins$ of
continuous data taking, data collection is momentarily halted, and the LED is
activated for a duration of $1~{}s$. Subsequently, the data acquisition
resumes once the detector returns to its designated operational temperature,
which usually takes another $20~{}mins$.

Using this experimental configuration, a combined dataset of $\sim 50~{}hrs$
featuring ${}^{252}Cf$ source was collected. The source was positioned at an
approximate distance of $2~{}ft$ from the detector. It was shielded by a layer
of 4” lead bricks, blocking the excess of gamma radiation. The detector was
operated under the optimal potential difference of $16~{}V$ across the
detector.

Figure 4: A plot of Y against the signal on the low-voltage side from 0 to 10
keV, featuring the mean along with $1\sigma$ and $2\sigma$ of the fitted bands
for electron recoil (red) and nuclear recoil (blue) events.

Figure 4 showcases the discriminator (Y), a fraction of the signal on the top
and bottom channels, plotted against energy within the energy window of
$0-10~{}keV$. It is important to note that the energy is represented in units
of electron recoil equivalent ($keV_{ee}$). As anticipated, discernible bands
corresponding to electron recoil (ER) and nuclear recoil (NR) events are
evident in the figure. Additionally, the fitted mean, $1\sigma$, and $2\sigma$
bands are superimposed on the plot. Here the polynomial fit was made on energy
bins with the size of $0.1~{}keV$. The ER event band appears uniformly
distributed, while the NR band exhibits energy-dependent behavior consistent
with expectations derived from the Lindhard theory [15]. The primary interest
of this experimental run is to assess the detector’s performance and its
discriminating power on near-threshold events. In Figure 5, energy region up
to $3~{}keV$ is displayed (x-axis in the logarithmic scale). As we see, this
detector demonstrates the capacity to measure nuclear recoil events while
effectively rejecting electron recoil background down to energy levels of
approximately $\sim 500~{}eV_{ee}$.

Figure 5: A plot of Y against the signal on the low-voltage side from 0 to 3
keV with the same mean, $1\sigma$ and $2\sigma$ of the fitted bands. Nuclear
recoil events are collected down to the energy levels of approximately $\sim
500~{}eV_{ee}$. Red lines represent electron recoil events, while blue lines
correspond to nuclear recoil events. The blob below $500~{}eV_{ee}$
corresponds to noise events.

Additionally, an extended background dataset was obtained in the absence of
any source to quantify the detector’s efficiency in rejecting background
events.

This estimation accounts for a fraction of the events within the NR band
divided by the entire event population. In addition to actual NR events, it is
important to note that ER events may appear within the NR band due to the
incomplete charge collection. The effect of incomplete charge collection is
characteristic of events happening near the conical surface/edges of the
detector. The effect is caused by either not enough electric field at the
peripheries to direct electrons to the neck or it can be attributed to crystal
defects and impurities. With this extended source free $16~{}hrs$ run, on
average, approximately $\sim 5~{}\%$ of total population events appear within
the NR $1\sigma$ band, leading to a substantial $95\%$ reduction in the
background.

Figure 6: Background run of the hybrid detector conducted throughout
$16~{}hrs$ without any source. Approximately $3\%-10\%$ of the events manifest
within the nuclear recoil band. Scaled event rate observed within the nuclear
recoil band. The estimated background rate above $1~{}keV$ is
$<900~{}Counts/keV/kg/days$ in the test facility without any effective
shielding techniques.

##  V conclusion

Following the initial successful run that demonstrated the hybrid detector’s
performance, we achieved a further reduction in the energy threshold by
implementing a software triggering scheme within the data processing pipeline.
In addition, the detector operation under the bias voltage was automated with
periodic LED flashing and cooling down, preventing significant charge build-up
within the detector and enabling continuous data taking for an extended time.
Consequently, we attained an energy threshold of approximately $\sim
500~{}eV_{ee}$, successfully demonstrating nuclear recoil/electron recoil
(NR/ER) discrimination at lower energies. The ER band is uniform while the NR
band follows the shape described by the Lindhard theory. Notably, the
discrimination power improves at lower energies due to the low energy
threshold. As a result, the electron recoil background reduction of $3-10\%$
has been achieved. Future studies will be directed toward the elimination of
incomplete charge collection which would further increase the background
rejection efficiency and reduce the background level for rare event searches.
New designs with improved field-shaping capabilities are under consideration.

##  VI acknowledgments

This work was fully supported by DOE grants DE-SC0017859 and DE-SC0018981. We
acknowledge the contribution of the key cryogenic infrastructure (Bluefors
LD400) provided by NISER, India. We would like to further acknowledge the
support of DAE through project Research in Basic Sciences - Dark Matter and
SERB-DST through the J.C. Bose fellowship.

## References

  * Zwicky [1933] F. Zwicky, The redshift of extragalatic nebulae, Helv. Phys. Acta 6, 110 (1933).
  * Rubin and W. K. Ford [1970] V. C. Rubin and J. W. K. Ford, Rotation of the andromeda nebula from a spectroscopic survey of emission regions, Astrophys. J. 159, 379 (1970).
  * Steigman and Turner [1985] G. Steigman and M. S. Turner, Review of mathematics, numerical factors, and corrections for dark matter experiments based on elastic nuclear recoil., Nucl. Phys. B , 253 (1985).
  * Goodman and Witten. [1985] M. W. Goodman and E. Witten., Detectability of certain dark matter candidates., Phys. Rev. D 31, 3059–306 (1985).
  * Agnese, et al.(2017) [SuperCDMS Collaboration] R. Agnese, et al.(SuperCDMS Collaboration) (SuperCDMS Collaboration), Projected sensitivity of the supercdms snolab experiment, Phys. Rev. D 95, 082002 (2017).
  * Aguilar-Arevalo, et al. [2016] A. Aguilar-Arevalo, et al. (DAMIC Collaboration), Search for low-mass wimps in a 0.6 kg day exposure of the damic experiment at snolab, Phys. Rev. D 94, 082006 (2016).
  * Barak et al. [2020] L. Barak et al. (SENSEI Collaboration), Sensei: Direct-detection results on sub-gev dark matter from a new skipper ccd, Phys. Rev. Lett. 125, 171802 (2020).
  * Abdelhameed, et al. [2019] A. H. Abdelhameed, et al. (CRESST Collaboration), First results from the cresst-iii low-mass dark matter program, Phys. Rev. D 100, 102002 (2019).
  * Agnese, et al.(2013) [SuperCDMS Collaboration] R. Agnese, et al.(SuperCDMS Collaboration), Demonstration of surface electron rejection with interleaved germanium detectors for dark matter searches, Appl. Phys. Lett. 103, 164105 (2013).
  * et al. [2014] R. A. et al. (SuperCDMS Collaboration), Search for low-mass wimps using voltage-assisted calorimetric ionization detection in the supercdms experiment, Phys. Rev. Lett. 122 (2014).
  * Neog  _et al._ [2022] H. Neog, R. Mahapatra, N. Mirabolfathi, M. Platt, A. Jastram, G. Agnolet, H. Chen, B. Mohanty, and A. Kubik, Phonon-mediated high-voltage detector with background rejection for low-mass dark matter and reactor coherent neutrino scattering experiments, Nuclear Inst. and Methods in Physics Research, A 1033 (2022).
  * Neganov and Trofimov [1985] B. S. Neganov and V. N. Trofimov, Colorimetric method measuring ionizing radiation, Otkryt. Izobret. 146, 215 (1985).
  * Luke [1988] P. N. Luke, Voltage‐assisted calorimetric ionization detector, Journal of Applied Physics 64, 6858 (1988), https://doi.org/10.1063/1.341976 .
  * K. D. Nam and Young [1995] B. C. K. D. Nam, S. W. Cabrera and B. A. Young, A quasiparticle‐trap‐assisted transition‐edge sensor for phonon‐mediated particle detection, Review of Scientific Instruments 66, 5322 (1995), https://doi.org/10.1063/1.1146105 .
  * Lindhard  _et al._ [1963] J. Lindhard, V. Nielsen, M. Scharff, and P. Thomsen, Integral equations governing radiation effects, Kong.Dan.Vid.Sel.Mat.Fys.Med. 33 (1963).
  * Gatti and Manfredi [1986] E. Gatti and P. Manfredi, Processing the signals from solid-state detectors in elementary particle physics, Riv. Nuovo Cim. 9, 1 (1986).
  * Zadeh and Ragazzini [1952] L. A. Zadeh and J. R. Ragazzini, Optimum filters for the detection of signals in noise, Proceedings of the IRE 40, 1223 (1952).
  * Belanger  _et al._ [1979] R. Belanger, D. Buckley, and J. Swensen, Environmental assesment of ionization chamber smoke detectors containing am-241 (1979).


# Is a Quantum Gravity Era Necessary?

Bogdan Veklych

###### Abstract

We present the first published framework of the entirety of cosmological
history which is thoroughly classical (without any quantum-gravitational era
or singularities) and which passes all the known extensive consistency checks
on such a model, and discuss some of its possible cosmological implications,
such as its ability to account for the matter-antimatter asymmetry, dark flow,
and the Hubble tension, albeit at the cost of further assumptions.

Author affiliation: University of Wisconsin-Madison<EMAIL_ADDRESS><EMAIL_ADDRESS>
orcid.org/0000-0003-1897-7998

Keywords: cosmological evolution, past-eternal cosmology, dark matter, Hubble
tension

##  1 Introduction

As is well-known, the current established cosmological knowledge of the
history of the Universe extends only 13.8 billion years into the past to a hot
dense rapidly expanding phase, nicknamed the Big Bang, which in turn is widely
conjectured to have been produced by the decay of an inflationary field - but
further than that it is not known with any certainty what preceded it or how
it came about. While it is easy to come up with speculative quantum-
gravitational ideas for the ultimate past of the Universe (such as e.g. origin
of Universes by budding from "quantum foam" in an eternal empty progenitor
Universe, or, as another example, an empty 2+1-dimensional compactified Milne
past-eternal contracting phase, stable due to rigidity of 2+1-dimensional
vacuum, which then changes dimensionality at the bounce via some hypothetical
quantum-gravitational mechanism), apparently there has never been an example
in the literature of a thoroughly classical consistent cosmological model of
the full prehistory of the cosmic inflation, with no singularities or
speculative quantum-gravitational epochs or other handwaved unclear "gaps in
history". This is however not surprising as there are numerous very tight
consistency constraints and no-go results on such a model [1][2][3][4][5][6],
and the purpose of the present article is to fill this theoretical gap by
introducing a framework with such desirable features that satisfies all these
consistency constraints.

##  2 The General Framework and a Particular Toy Model Illustration

The basic description of the model is as follows. In the pre-inflationary
"burning fuse era", the Universe has two spatial dimensions compactified (i.e.
it is shaped like an "infinite corridor" with flat torus cross-sections), with
two waves approaching each other along the third dimension with the speed of
light from past-eternity, while the Universe is static, flat, and empty in the
shrinking middle part between these two waves. When these primordial waves
finally collide and react, part of their contents forms proto-inflationary
scalar field(s), and another part is matter with nonzero intrinsic angular
momentum (like e.g. ordinary photons with spin). The Universe may start
collapsing, but due to the matter it bounces back - here we, crucially, assume
validity of Einstein-Cartan theory - after which the true inflation starts.
(It is very plausible that this setup can be slightly altered to have a pure
Einstein-Cartan bounce with no inflation after the wave collision, but for
simplicity we won’t consider this here.)

We shall now give a particular proof-of-concept toy illustration of this
model, with some additional simplifying assumptions, and then discuss the
consistency of our past-eternal cosmology. Note that a long chain of highly ad
hoc assumptions, akin to an eternal Rube Goldberg machine, is fundamentally
inevitable in any realization of this basic idea (and this is the trade-off
for really avoiding quantum gravity); see the end of the next section for why
this is completely benign. Our starting point is the effect that a perfect
fluid with spin together with a scalar field have on the metric in Einstein-
Cartan theory (see the equation (77) in [7]):

$\displaystyle G_{kl}=\kappa\left(\rho-\frac{1}{8}\kappa
s_{bd}s^{bd}\right)u_{k}u_{l}+\kappa\left(p-\frac{1}{8}\kappa
s_{bd}s^{bd}\right)(g_{kl}+u_{k}u_{l})$ (1)
$\displaystyle[+\kappa\Theta_{kl}]-\frac{1}{2}\kappa(g^{bd}+u^{b}u^{d})\nabla_{d}(s_{bk}u_{l}+s_{bl}u_{k})$

where $G_{kl}$ is the Einstein tensor (in the ordinary sense, constructed from
the Levi-Civita connection, likewise for the covariant derivative),
$\kappa=8\pi G/c^{4}$ (henceforth we set $\kappa=c=1$), $g_{kl}$ is the metric
tensor, $\rho,p,u^{k},s_{kl}$ are the density, pressure, four-velocity, and
spin density of the fluid, correspondingly, and $\Theta_{kl}$ is the
(canonical and symmetric) stress-energy tensor of a scalar field $\phi$ (which
can be incorporated into $\rho$ and $p$, hence the brackets). From (1) we can
derive the effect on the metric of spinning null dust by a simple limit
procedure, as follows. Let the spin and velocity point along the
$x$-direction, i.e. let $s_{kl}=0$ for $k+l\neq 5$ and $u^{2}=u^{3}=0$ (note
that the expression $u^{b}u^{d}\nabla_{d}(s_{bk}u_{l}+s_{bl}u_{k})$ now
vanishes), and take the limit as $u^{0}\rightarrow\infty$ while keeping
$(u^{0})^{2}\rho,(u^{0})^{2}p,u^{0}s_{kl}$, and $u^{1}/u^{0}$ constant. In the
end, we obtain

$\displaystyle
G_{kl}=\Theta_{kl}+\left(\epsilon-\frac{1}{4}s_{bd}s^{bd}\right)n_{k}n_{l}-\frac{1}{2}g^{bd}\nabla_{d}(s_{bk}n_{l}+s_{bl}n_{k})$
(2)

where $n^{k}$ is a null vector and $\epsilon(n^{0})^{2}$ is the energy density
of the null dust. Further, in the case of a simple plane wave characterized by
$n^{0}=1,n^{1}=-1,s_{23}=-s_{32}=s(x+t),\epsilon=\epsilon(x+t),g_{kl}=\operatorname{diag}(-1,1,h(x+t)^{2},h(x+t)^{2}),\phi=\phi(x+t)$,
(2) reduces to

$-2\frac{h^{\prime\prime}}{h}=(\phi^{\prime})^{2}+\epsilon-\frac{s^{2}}{2h^{4}}$
(3)

With an increasing numerical density of identical particles, $s^{2}$ grows
faster than $\epsilon$, which means that for concentrated enough matter the
right hand side of (3) can become negative, making the space that these
particles reach expand (or "bounce") transversally. For maximum simplicity,
let’s hold the matter concentration at exactly "the brink of bounce",
furthermore taking $\epsilon$ to be negligibly small, reducing the metric to
flat ($h=\operatorname{const}$, say $h=1$) and (3) to
$(\phi^{\prime})^{2}=\frac{1}{2}s^{2}$, and similarly for the left wave. Let’s
also assume that once this right wave collides with the symmetric left wave at
$t=0$, the intermediate interaction region $t>0,|x|<t$ will be homogeneous
(with all quantities dependent at most on $t$) and flat (but anisotropic -
having a preferred direction, in which the spin is pointing), which, by (1),
is achievable if

$\frac{1}{4}s(t)^{2}=\frac{1}{2}(\phi(t)^{\prime})^{2}$ (4)

(as $\rho=p=\frac{1}{2}(\phi(t)^{\prime})^{2}$ for the scalar field). To this
end, let’s make the following further toy model assumptions: the spin density
in the waves is constant, and so is the kinetic energy of the scalar field
(i.e. the field is linearly changing with $t\pm x$); the spinning particles
immediately stop and become spinning dust when reaching the middle region (for
example, due to a reaction that amounts to an inelastic collision with some
spinless particles traveling in the opposite direction). Now observe that 1)
under these conditions spin becomes twice more concentrated in the middle
region than in the waves - a spinning particle in a wave which is at a
distance $2\delta$ away from the interaction front will move $\delta$ more
before being stopped - thus $\frac{1}{4}s(t)^{2}$ in the middle region is
twice larger than $\frac{1}{2}s^{2}$ in the waves; 2) if
$\phi(x+t)=f_{1}(x+t)$ in the right wave and $\phi(t)=f_{2}(t)$ in the middle
region, from the fact that they match at $t=x$ we conclude that
$f_{2}(t)=f_{1}(2t)$, thus $\frac{1}{2}(\phi(t)^{\prime})^{2}$ in the middle
region is twice larger than $(\phi^{\prime})^{2}$ in the waves. These two
observations imply that (4) is satisfied in the intermediate region if (3) was
satisfied in the waves, and so our assumptions are consistent. (As a reminder,
we assumed, for simplicity, that up to this point the metric was flat, which
simplified the equations (1)-(3) immensely.)

Now, to enable the beginning of inflation, we assume (yet further) that the
scalar field has a stepfunction potential, so that it is only effectively
massless as long as it is below a certain value. (A more natural, but less
computationally convenient, option is to take two scalar fields $\psi_{left}$
and $\psi_{right}$ with a combined polynomial or other smooth potential that
includes a factor of $\psi_{left}^{2}\psi_{right}^{2}$, so that either one is
effectively massless when the other one is zero, and have $\psi_{left}$ be
zero to the right of the left wave and $\psi_{right}$ be zero to the left of
the right wave.) Specifically, assume that when it hits that value at some
moment $t_{0}$ (which can be as large as needed) all of its kinetic energy
becomes potential energy (and we also assume that in the waves the growth of
$\phi$ is truncated at infinitesimally less than that value). From now on
we’re only interested in the region causally uninfluenced by anything outside
of the homogeneous middle part at the moment of this transition. By definition
this region is shrinking inward with the speed of light on its right and left
edges, and internally it is subject to FLRW-like equations (see [8]; in the
following $a(t)=g_{11}=g_{22}=g_{33}$ is the scale factor, depending only on
time, $a(t_{0})=1$, $a^{\prime}(t_{0})=0$, $V$ is the height of the potential
of the scalar field, and $s_{0}^{2}$ is constant)

$3\frac{(a^{\prime})^{2}}{a^{2}}+\frac{s_{0}^{2}}{4a^{6}}=V=2\frac{a^{\prime\prime}}{a}+\frac{(a^{\prime})^{2}}{a^{2}}-\frac{s_{0}^{2}}{4a^{6}}$
(5)

From the first equality in (5) we get $(a^{\prime}/a)^{2}\leq V/3$,
substituting this into the second equality we get $a^{\prime\prime}/a\geq
V/3$; in terms of $A(t)=\ln a(t)$ this reads
$(A^{\prime})^{2}+A^{\prime\prime}\geq V/3$, thus $A^{\prime}$ increases from
zero at first and cannot decrease if it is less than $\sqrt{V/3}$, so $A$ is
bounded below by a positively-sloped line, and therefore $a$ (eventually)
grows exponentially. From this it also follows that, if $t_{0}$ is large
enough, this region will expand exponentially in the $x$-direction as well
(that is, with its linearly inward-moving left and right boundaries never
overtaking its exponential expansion in the middle) - the (toy-model)
inflation has started.

##  3 Consistency

In this section we show how this general framework (with any secondary
details) circumvents the numerous tight consistency constraints mentioned
before.

Penrose-Hawking singularity theorems [1] do not apply to Einstein-Cartan
theory.

Compactification of all but one spatial dimensions automatically stabilizes
and homogenizes plane waves. For example, suppressing one spatial dimension
for convenience (and limiting ourselves only to special relativity), consider
an arbitrary "splash" of a massless scalar field $\psi$ in a flat 2+1
dimensional Universe with one compactified spatial dimension of (minimal
incontractible loop) length $2\pi$, in the appropriate units of length.
Applying the Klein-Gordon equation $\square\psi=0$ to the Fourier
decomposition of $\psi$ along the periodic coordinate,
$\psi(t,x,y)=\sum_{k=-\infty}^{\infty}\psi_{k}(t,x)e^{iky}$, we get the set of
equations $(\square_{t,x}-k^{2})\psi_{k}(t,x)=0,k\in\mathds{Z}$. In
particular, taking $k=0$, we see that
$\psi_{0}(t,x)=\frac{1}{2\pi}\int_{0}^{2\pi}\psi(t,x,y)dy$ satisfies the one-
dimensional wave equation for speed 1 and (after possibly splitting into left-
and right- moving parts) will be indefinitely stably propagating.

There is no quantum dispersion [2] because any relevant wavefunction is
relativistically "frozen" by motion with the speed of light. And even
Heisenberg uncertainty inequality does not apply in the case of a periodic
coordinate: for example, in the textbook "particle in a ring" scenario the
particle can have zero angular momentum, even though its angular position is
only finitely imprecise (having the maximum possible imprecision, $2\pi$).

This model doesn’t have $H>0$ [3] (nor $H<0$ [4] with its environment of
unlimited blueshift, where any fluctuation - vanishing in the limit of past-
infinity - becomes deadly magnified). In fact, in this model any geodesically
moving observer was always, before a certain moment, simply in static flat
empty space. (Note, importantly, that any proposals in which the Universe has
a "life cycle", e.g. originating in a black hole in an older Universe, which
did likewise, "and so on", do not pass [3].)

This model also gives an example of how it is possible for "the Universe [to
be] in a near-equilibrium state for the first "half" of eternity and then …
[to exit] this equilibrium" ([5], p.29); indeed, it can be explicitly shown
that there is no entropy change (in particular, no entropy decrease) before
the wave collision, as follows. Considering the effect of one of the waves, we
can replace distant causally not yet reached regions, including those with the
other wave, by more of empty space compactified in two dimensions; now, this
Universe looks exactly the same after some time as before - if we just shift
our vantage point by the distance the wave has traveled - so its entropy is
also the same. Thus we conclude that there is no entropy change in the pre-
collision era even under some yet-unknown possible future quantum-
gravitational understanding of entropy, in the spirit of [5].

Static flat empty space is not subject to the quantum instabilities of [6]
(confirmed by A. Vilenkin in personal correspondence) due to fundamental
problems with assigning consistent decay-governing wavefunctions to it. And,
finally, we also need to make the very reasonable conjecture that static flat
empty space does not spontaneously produce matter. (It would be quite a
challenge to disprove it - such a disproof would yield the simplest _perpetuum
mobile_ design of all, just a large empty container, inside of which matter is
spontaneously generated and collected!) Another way to put this point is that
quantum fluctuations in static flat vacuum do not increase entropy.

Note that this model is extremely conservative (in particular, the
characteristic peak density reached during an Einstein-Cartan bounce, called
the Cartan density, is about forty orders of magnitude less than the Planck
density, so quantum gravity is not relevant during the wave collision),
corresponding precisely to the classical picture argued for by natural
philosophers like Aristotle and Giordano Bruno: time, space, motion, and
matter have existed forever, with matter always coming only from the matter
that existed in the preceding moments. (As to its fine-tuning, it can be
assumed, for example, that our Universe is a part of a larger Multiverse with
all sorts of physical parameters in different places - a common suggestion,
also applicable to this model; this "anthropic principle" can cover any
remaining loose ends in the model that might seem strange or unlikely but are
not outright demonstrably excluded - which is why the focus of the above was
pure consistency. Any particular realization of the idea above will be
something that can be described impressionistically as messy, ad hoc, and
convoluted - but arguably the same can be said about the Standard Model of
particle physics, for example, so this is not an argument for its falsity.)

##  4 Possible cosmological implications?

If we additionally assume that inflation only existed briefly (rather than as
eternal chaotic inflation, which would dilute away any traces of its
prehistory), many intriguing possibilities open up if this framework is
correct - but, as always with this "Rube Goldberg cosmology", at the cost of
further convoluted assumptions:

1) Let’s assume that the lightest neutrino mass eigenstate is, surprisingly
enough, exactly zero, which in turn enables us to make another assumption -
that the lightest neutrino was present in the primordial waves, in greater
abundance than the corresponding antineutrino. Then the matter-antimatter
imbalance is trivially accounted for under our framework - there was simply
always more matter than antimatter, both before and after the collision of the
primordial waves!

2) Let’s also assume that dark matter, likewise, had its own lightspeed
precursor(s) in the primordial waves, which reacted upon collision and formed
(nonzero-mass) dark matter - and that neither this precursor nor dark matter
itself interacts with luminous matter other than gravitationally. Then there
is no reason to assume that the inertial frame in which the average momentum
of luminous matter post-collision was zero was the same as the corresponding
frame for dark matter. In other words, we can assume that (massive) dark and
luminous matter originally were moving with respect to each other in the
direction of the uncompactified direction, and could have kept some of this
relative motion even at the time when the cosmic microwave background
radiation was emitted. Later, by mutual local gravitational tugs, they would
naturally finally "stick together" and equalize their average velocities -
which would mean that the visible matter would be now slightly moving on
average with respect to the inertial system defined by the cosmic microwave
background radiation (which is the system reflecting its motion before the
slight residual "final adjustment" by dark matter). But this is indeed
precisely what we observe, as the "dark flow" phenomenon. In other words, our
framework provides an explanation for the mysterious dark flow - and if this
explanation is correct, then, in turn, the dark flow gives us the direction of
the uncompactified dimension.

3) Finally, let’s furthermore assume that dark matter is "fuzzy" (which has
been recently confirmed [9]) but not scalar (that is, its particles are
extremely light and have spin 1 or 2), and also (similarly to the speed
imbalance assumption in 2) above) that the spins of the dark matter lightspeed
precursors weren’t precisely balanced to zero average in the colliding
primordial waves. Then dark matter now has net nonzero intrinsic angular
momentum pointing in the direction of the uncompactified dimension, inherited
from its primordial wave lightspeed progenitors. If so, the Universe is
homogeneous but not isotropic, having a preferred direction (in which the
intrinsic angular momentum of fuzzy dark matter is pointing), and is described
by [8] instead of a FLRW metric - which can explain the Hubble tension neatly
[10], and this explanation is even potentially experimentally checkable, via
the observational consequences of the corresponding spacetime torsion!

That is, if dark matter is fuzzy and has intrinsic angular momentum (pointing
in the direction of the dark flow or opposite to it, by 2), under our chains
of assumptions), then the standard cosmological Lambda-CDM model can be
augmented, with the relevant FLRW-like equations becoming (1) in [8], and the
torsion produced by the corresponding nonzero value of $\mu=K_{0}^{2}$ (in the
notations of [8]) is a prediction that can (and should!) be tested
experimentally [11].

In fact, as the dark flow demonstrates, there might indeed be a preferred
direction in the Universe, and if the assumption of its isotropy is dropped,
then the above one-parameter extension of the Lambda-CDM model is by far the
simplest possible explanation of the Hubble tension even without the
corroborating background provided by the general classical framework proposed
in this article, and therefore it is a reasonable candidate for being given an
observational test. (And there might also be other independent corroborating
evidence explained by this model, such as hints of a macroscopic compactified
dimension [12].)

References

[1] S. W. Hawking, "Singularities in the Universe", _Phys. Rev. Lett. 17_ ,
444

[2] A. Aguirre, J. Kehayias, "Quantum instability of the emergent universe",
_Physical Review D, 88(10)_ , arxiv.org/abs/1306.3232

[3] A. Borde, A. H. Guth, A. Vilenkin, "Inflationary spacetimes are not past-
complete", _Phys. Rev. Lett. 90_ (2003) 151301, abs/gr-qc/0110012;

W. H. Kinney, N. K. Stein, "Cyclic Cosmology and Geodesic Completeness",
_JCAP_ (accepted for publication), abs/2110.15380

[4] T. Banks, W. Fischler, "Black Crunch", abs/hep-th/0212113

[5] A. C. Wall, "The Generalized Second Law implies a Quantum Singularity
Theorem", _Class. Quantum Grav. 30_ , 165003 (2013), abs/1010.5513

[6] A. T. Mithani, A. Vilenkin, "Collapse of simple harmonic universe",

abs/1110.4096;

A. T. Mithani, A. Vilenkin, "Did the universe have a beginning?",
abs/1204.4658

[7] K. Pasmatsiou, C. G. Tsagas, J. D. Barrow, "Kinematics of Einstein-Cartan
universes", _Phys. Rev. D, 95(10)_ :104007 (2017), abs/1611.07878

[8] W. Kopczynski, "An anisotropic universe with torsion", _Phys. Lett. A,
43(1)_ , 63-64 (1973)

[9] A. Amruth et al., "Anomalies in Gravitational-Lensed Images Revealing
Einstein Rings Modulated by Wavelike Dark Matter", _Nat. Astron. (2023)_ ,
abs/2304.09895

[10] F. Izaurieta, S. Lepe, O. Valdivia, "The Spin Tensor of Dark Matter and
the Hubble Parameter Tension", abs/2004.13163;

T. Liu et al., "Revisiting Friedmann-like cosmology with torsion: newest
constraints from high-redshift observations", _JCAP_ (accepted for
publication), abs/2304.06425

[11] K. Bolejko et al., "Cosmological signatures of torsion and how to
distinguish torsion from the dark sector", _Phys. Rev. D, 101_ (2020) 104046,
abs/2003.06528

[12] Y. Akrami et al., "The Search for the Topology of the Universe Has Just
Begun", abs/2210.11426