# Multi-component Predictions of Transient Solution Fields with Sequential
Deep Operator Network

Junyan He1 Shashank Kushwaha1 Jaewan Park1 Seid Koric1,2<EMAIL_ADDRESS>Diab Abueidda2<EMAIL_ADDRESS>Iwona Jasiuk1 1 Department of Mechanical
Science and Engineering, University of Illinois at Urbana-Champaign, Urbana,
IL, USA
2 National Center for Supercomputing Applications, University of Illinois at
Urbana-Champaign, Urbana, IL, USA

###### Abstract

The Deep Operator Network (DeepONet) structure has shown great potential in
approximating complex solution operators with low generalization errors.
Recently, a sequential DeepONet (S-DeepONet) was proposed to use sequential
learning models in the branch of DeepONet to predict final solutions given
time-dependent inputs. In this novel work, the S-DeepONet architecture is
further extended by modifying the information combination mechanism between
the branch and trunk networks to simultaneously predict vector solutions with
multiple components at multiple time steps of the evolution history. Two
example problems, one on transient fluid flow and the other on path-dependent
plastic loading were shown to demonstrate the capabilities of the model to
handle different physics problems. The use of a trained S-DeepONet model in
inverse parameter identification via the genetic algorithm is shown to
demonstrate the application of the model. In almost all cases, the trained
model achieved an $R^{2}$ value of above 0.99 and relative $L_{2}$ error of
less than 10% with only 3200 training data points, indicating superior
accuracy. The vector S-DeepONet model, having only 0.4% more parameters than a
scalar model, can predict two output components simultaneously at an accuracy
similar to the two independently trained scalar models with a 20.8% faster
training time. The S-DeepONet inference is at least three orders of magnitude
faster than direct numerical simulations, and inverse parameter
identifications using the trained model is highly efficient and accurate.

###### keywords:

Deep Operator Network (DeepONet) , Gated recurrent unit (GRU) , Computational
fluid dynamics (CFD) , Plastic Deformation , Genetic algorithm (GA) ,

††journal: Computational Methods in Applied Mechanics and Engineering

*[inlinelist,1]label=(),

##  1 Introduction

The recent convergence of high-performance computing and machine learning (ML)
techniques has led to significant strides in the neural network (NN)
applications in various fields, including structural optimization, flow
prediction, additive manufacturing, exploring structure-property relations,
and beyond [1, 2]. Various data-driven and physics-informed NN architectures
have been used to study nonlinear forward and inverse problems for various
applications, including additive manufacturing [3], materials design [4], and
constitutive modeling [5], and coupled problems [6]. Fuhg and Bouklas [7] and
He et al. [8] presented a deep energy method (DEM) utilizing potential energy
to resolve nonlinear material responses. Moreover, trained NN models combined
with genetic algorithms inspired by biological evolution opened a new era for
data mining, parameter tuning, optimizations, and many more [9]. The algorithm
creates offspring from the previous generation and inherits the elite genome,
having only competent ones survive. As an effective optimization method
itself, there are various applications of genetic algorithms to engineering
problems. Corriveau et al. [10] adopted the method in design optimization
combined with finite element analysis, while Abueidda et al. [11] combined
genetic algorithms with a neural network model to fine-tune the sparse data.
Wang et al. [12] proposed an inverse designing process of shell-based
mechanical metamaterials. A genetic algorithm is introduced as a solution of
inverse design, sorting out the best design among various random designs using
trained ANNs in the selection process.

Lately, Deep Operator Network has been introduced by Lu et al. [13]. Widely
known as the DeepONet, the neural network architecture learns an operator from
the governing differential equation of the system, which enables retraining
not necessary despite the input to the system, such as parametric functions,
boundary conditions, or even varying geometry. DeepONet consists of two
separate neural networks, named branch and trunk. The branch network takes the
parametric function as an input, while the trunk network takes domain
geometry. Both the branch and trunk were fully connected neural networks
(FNNs) in the very first form suggested by Lu et al. [13], encoding their
input information each, respectively. The output of the branch and trunk were
then underwent a dot product to map the parametric functions to the target,
the solution to the governing differential equation. Koric et al. [14]
successfully predicted stress field on a 2D domain with a plastic deformation.
Moreover, DeepONets lately have been frequently applied to solve science and
engineering problems such as inverse designing nanoscale heat transport system
[15], predicting elastic-plastic stress on convoluted topology-optimized
domain [16], earthquake localization by Haghighat et el. [17], heat conduction
problem by Koric and Abueidda [18], and nonequilibrium thermodynamics of
chemical mixtures combined with physics-informed neural networks by Li et al.
[19]. All these applications, while successfully converged to their targets,
did not provide full-field time-dependent output. FNN network in the DeepONet
has a limit to reflect the full nature of time-dependent input, abundantly
appearing in engineering problems.

The prediction of stress fields in solid mechanics and flow fields in fluid
mechanics has seen remarkable advancements with the adoption of NNs. Various
NN architectures, including the convolutional neural network named FR-CNN
established by Nie et al. [20] and conditional generative adversarial network
(GAN) tested by Hoq et al. [21], have been used to predict full-field stress
predictions. Similarly, in the field of fluid mechanics, CNNs and GANs have
been utilized to predict the velocity and pressure fields for various
geometries and flow conditions. [22, 23]. At the same time, the time-dependent
aspect of the input information can only be captured using recurrent neural
networks (RNNs) like long short-term memory (LSTM) [24] and gated recurrent
unit (GRU) [25]. While successfully avoiding the vanishing gradient issue of
simpler recurrent neural networks [26], both architectures can delete obsolete
memory intentionally, keeping only the vital information. Various works
underscore the effectiveness of these networks in time-dependent engineering
problems, such as plastic deformation [27] and flow prediction for the final
state [28, 29, 30]. RNNs used in these studies could accurately predict the
one-dimensional sequence of output data. The previous studies did not study
the effects of combining a recurrent neural network in the DeepONet structure
except the singular work on S-DeepONet [31], which discussed the architecture
combining RNN and DeepONet to predict outcomes for the last time step,
particularly for scalar output. Furthermore, there is an application using
ANNs in a genetic algorithm [12], yet no other work utilized DeepONets in
inverse problems with incorporated genetic algorithms.

Real-world problems like thermo-mechanical analysis, path-dependent
constitutive models, and lid-driven cavity flow often exhibit variable loading
conditions, resulting in a time-dependent response in the system. Using
classical techniques, including finite element analysis (FEA), topology
optimization (TO), and sensitivity analysis to obtain a time-dependent
response is computationally expensive and time-consuming. Hence, a scientific
need exists for more robust data-driven surrogate models for time-dependent
nonlinear problems capable of predicting full-field contours for multiple time
steps and vector components. However, as discussed earlier in the manuscript,
the classical DeepONet, although capable of predicting a full-field solution,
lacks the capability of causal relationships in input data with its FNN branch
network. Standalone recurrent neural networks like GRUs and LSTMs capture
time-dependent signals and predict a full-field output without considering
spatial information. Even when RNN is combined with DeepONet, only the final
stress state or temperature distribution was predicted by this architecture by
He et al. [31]. Therefore, using RNN models in the DeepONet branch-trunk
architecture to predict full-field contours over multiple time steps is a
first in the literature and is the most significant and novel scientific
contribution of this work. We have named the proposed sequential DeepONets as
S-DeepONets. In this work, we tested our approach with two example problems:
(1) a lid-driven cavity flow and (2) a path-dependent plasticity problem, both
involving significant time dependence and representing real-world engineering
use cases. We used the proposed innovative S-DeepONet formulation to predict
full-field solutions and compared performance with the classical DeepONet
formulation.

This paper is organized as follows: Section 2 introduces neural network
architectures and provides details on the data generation method. Section 3
presents and discusses the performance of the NN model. Section 4 summarizes
the outcomes and limitations and highlights future works.

##  2 Methods

###  2.1 S-DeepONet with multiple output dimensions

The key idea in the original S-DeepONet architecture as proposed by the
authors [31] is the separation of the temporal component (handled by a GRU
branch network) and the spatial component (handled by a FNN trunk network) of
the solution operator, and combining the components via a matrix-vector
product and a bias:

$\bm{G}_{n}=\sum_{h=1}^{HD}\bm{B}_{h}\bm{T}_{nh}+\beta,$ (1)

where $\bm{G}$, $\bm{B}$ and $\bm{T}$ denote the outputs of the S-DeepONet,
branch network, and trunk network, respectively. Dimension index $h$
represents the hidden dimension (HD) of branch and trunk networks, and index
$n$ represents the flattened spatial dimension, which contains $N$ nodes in
the simulation domain. Finally, $\beta$ is a bias added to the product. This
structure allows the prediction of a full-field solution discretized by $N$
nodes, where the time-dependent input load information is encoded in the
branch network and the spatial geometry information is encoded in the trunk
network.

Although including time-dependent information in the input load, the original
S-DeepONet only predicts a scalar solution field at the end of the load. To
extend the S-DeepONet architecture to predict vector solution fields at
different time steps, we further exploit the idea of separation of spatial and
temporal components into the trunk and branch networks. To this end, a novel
S-DeepONet structure is proposed, and its schematic is shown in Fig. 1.

Figure 1: Schematic of the improved S-DeepONet used in this work. $m_{V}$,
$x$, $y$ and $\hat{G}$ denote the load magnitude, X coordinates, Y
coordinates, and the approximated solution operator. Dimensions $HD$, $S$, $N$
and $C$ are the number of hidden dimensions, time steps, nodes in the domain
and output vector components, respectively. The black GRU blocks return a
sequence (2D outputs), while the green GRU block compresses the output into
1D.

In the proposed architecture, we leverage the GRU branch network to produce
encoded hidden outputs for all the output time steps as a tensor $\bm{B}$ of
shape $[HD\times S]$, where $S$ is the number of time steps. Similarly, when
provided an $[N\times 2]$ input vector containing 2D coordinates of all nodes
in the domain, the trunk network in the proposed S-DeepONet produces an
encoded output tensor $\bm{T}$ of shape $[N\times HD\times C]$, where $C$ is
the number of output vector components. $\bm{T}$ contains the encoded hidden
outputs for all nodes and all vector components. To account for the new output
dimensions, we combine the information from the branch and trunk via the
following product:

$\bm{G}_{nsc}=\sum_{h=1}^{HD}\bm{B}_{hs}\bm{T}_{nhc}+\beta,$ (2)

where the lower-case indices $s$ and $c$ correspond to the time step and
output component dimensions. The tensor product nature of this combination
allows for efficient simultaneous generation of full-field vector outputs at
multiple output time steps. This combination can be understood as a
simultaneous identification of a set of ”basis” shapes (from trunk network)
and corresponding weights (from branch network), where the final solution
contours are expressed as the weighted linear combination of those basis
contours. This idea of basis and weight identification of the S-DeepONet
architecture is illustrated further in Fig. 2.

Figure 2: Basis and weight interpretation of the S-DeepONet combination
operator. For each vector component $c$, S-DeepONet identifies $HD$ basis
shapes from the trunk network. Simultaneously, for each time step $s$, $HD$
weights from the branch network. The solution contour is expressed as the
linear combination of all basis shapes including the data mean.

By predicting $C$ output vector components together, we are eliminating the
need to train $C$ different single-component S-DeepONets, each specialized in
one output component, thus leading to tremendous computational time savings.

The branch network of the developed S-DeepONet consists of four GRU layers.
The encoding-decoding structure of the original S-DeepONet is adopted, with
the first two being the encoder and the last two being the decoder. All GRU
layers use a tanh activation function. Finally, a time-distributed dense layer
with linear activation is used to output the branch results with a hidden
dimension ($HD$) of 32 for all $S$ time steps. The trunk network is a FNN with
the following seven layers of neurons: $[2,101,101,101,101,101,HD\times C]$,
where $C$ takes different values in the numerical examples presented in this
work. All NNs were implemented in the DeepXDE framework [32] with a TensorFlow
backend [33]. All models were trained for 300000 epochs with a batch size of
64. The Adam optimizer [34] was used and the scaled mean squared error (MSE)
was used as the loss function, which is defined as:

${\rm{MSE}}=\frac{1}{N}\sum^{N}_{i=1}(f_{FE}-f_{Pred})^{2},$ (3)

where $N$, $f_{FE}$, and $f_{Pred}$ denote the number of data points, the FE-
simulated field value, and the NN-predicted field value, respectively.

###  2.2 Data generation

In this work, we demonstrate the application of the proposed improved
S-DeepONet in time-dependent computational fluid dynamics (CFD) and history-
dependent plastic deformation problems. Using the trained S-DeepONet, inverse
parameter identification is performed with genetic algorithm.

####  2.2.1 Lid-driven cavity flow

In the first example, we consider a classical CFD benchmark problem: the lid-
driven cavity flow. Consider a rectangular cavity of length 3 and height 1,
which is discretized into a $121\times 41$ uniform grid with 4961 nodes. A lid
is located at $y=1$ and moves at a time-dependent velocity profile
$\bar{u}(t)$ in the $X$ direction, driving the fluid motion in the cavity. The
radial basis interpolation (RBI) was used to generate the smooth lid velocity
profile, which are defined by six uniformly spaced control points. The
starting velocity (i.e., at $t=0$) is 0, and the velocity in all other control
points were sampled from the range $[-2,2]$. The fluid is assumed to be
incompressible, and the system is governed by the Navier-Stokes equation:

$\frac{\partial\bm{u}}{\partial
t}=\mu\nabla^{2}\bm{u}-(\bm{u}\cdot\nabla)\bm{u}-\frac{1}{\rho}\nabla P,$ (4)

where $\bm{u}$, $\mu$, $\rho$, and $P$ denote the velocity vector, viscosity,
mass density, and pressure, respectively. In this example, a hypothetical
fluid with $\rho=1$ and $\mu=0.1$ was used. In 2D, let the $X$ and $Y$
components of the velocity vector be denoted as $u$ and $v$. No-slip boundary
condition was used for all four boundaries:

$\displaystyle u=\bar{u}(t),\,v=0\;\;\rm{if}\;y=1,$ (5)
$\displaystyle\bm{u}=\bm{0},\;\;\rm{Otherwise}.$

For the pressure degree of freedom, the following boundary conditions were
used:

$\displaystyle P=0\;\;\rm{if}\;y=1,$ (6) $\displaystyle\frac{\partial
P}{\partial x}=0\;\;\rm{if}\;x=0,3,$ $\displaystyle\frac{\partial P}{\partial
y}=0\;\;\rm{if}\;y=0.$

The trivial initial conditions $\bm{u}=\bm{0}$ and $P=0$ were used. For this
simple geometry, the second-order central finite difference (FD) method was
used to discretize the governing equation, and the pressure-projection method
with explicit time integration was used to evolve the system for 10000 time
steps with a time step size of $2\times 10^{-4}$. An Python code was used to
solve this problem, which was adopted from the work of Barba et al. [35]. A
total of 4000 simulations were conducted with distinct lid velocity profiles.
The primary variables $P$, $u$, and $v$ were stored at 25 uniformly spaced
output time steps over the simulation period. For better NN training and
accuracy, it is best to scale the input and output data to suitable ranges. To
account for the change in solution scale as the fluid flow develops, a time-
dependent scale factor for each output variable was first calculated. For a
time step $i$, this scale is defined as the maximum absolute value of the
field values (over all cases and all points) at this time step. The data at
time step $i$ was then divided by this scale, leading to a new data scale of
$[-1,1]$.

####  2.2.2 Dog-bone axial loading

In the second example, we consider the plastic deformation of a dog-bone
specimen under time-dependent axial loads. The specimen has a length of 110 mm
and gauge width of 20 mm. A total of 4756 linear plane stress elements were
used, and the specimen is fixed on its left end with displacement applied on
the right end. A schematic of the specimen and boundary conditions are shown
in Fig. 3a.

(a)

(b)

Figure 3: LABEL:sub@geom Schematic of the dog-bone specimen and the mesh used
in FE simulations. The axial displacement is along the global X direction.
LABEL:sub@disp A typical applied displacement in the plastic deformation
problem.

In the absence of any body and inertial forces, the equilibrium equations and
boundary conditions are:

$\displaystyle\nabla\cdot\bm{\sigma}=\bm{0},\;\;\forall\bm{X}\in\Omega,$ (7)
$\displaystyle\bm{u}=\bar{\bm{u}},\;\;\forall\bm{X}\in\partial\Omega_{u},$

where $\bm{\sigma}$ and $\bar{\bm{u}}$ denote the Cauchy stress and prescribed
displacement, respectively. Small-strain assumption is applied, which leads to
the following definition of total strain tensor:

$\bm{\epsilon}=\frac{1}{2}(\nabla\bm{u}+\nabla\bm{u}^{T}),$ (8)

as well as its additive decomposition into elastic and plastic strain parts:

$\bm{\epsilon}=\bm{\epsilon}^{e}+\bm{\epsilon}^{p}.$ (9)

For linear elastic isotropic material in plane stress, stress can be found
from elastic strain via:

$\begin{bmatrix}\sigma_{11}\\\ \sigma_{22}\\\
\sigma_{12}\end{bmatrix}=\frac{E}{1-\nu^{2}}\begin{bmatrix}1&\nu&0\\\
\nu&1&0\\\ 0&0&1-\nu\end{bmatrix}\begin{bmatrix}\epsilon^{e}_{11}\\\
\epsilon^{e}_{22}\\\ \epsilon^{e}_{12}\end{bmatrix},$ (10)

where $E$ and $\nu$ are the Young’s modulus and Poisson’s ratio. Plasticity is
modeled with the simple $J_{2}$ plasticity with linear isotropic hardening
law:

$\sigma_{y}=\sigma_{y0}+H\bar{\epsilon}_{p},$ (11)

where $\sigma_{y}$, $\bar{\epsilon}_{p}$, $\sigma_{y0}$, and $H$ denote the
flow stress, equivalent plastic strain, initial yield stress, and the
hardening modulus, respectively. Relevant material properties are shown in
Table 1. $\bar{\epsilon}_{p}$ is an internal history variable in plasticity
that records the accumulation of plastic strain at different points in the
loading history. Due to its vital role in plasticity, this is considered one
of the output variables of the S-DeepONet, while the other output variable is
the von Mises stress:

$\bar{\sigma}=\sqrt{\sigma_{11}^{2}+\sigma_{22}^{2}+\sigma_{11}\sigma_{22}+3\sigma_{12}^{2}}.$
(12)

Table 1: Material properties of the elastic-plastic material model Property |  $E$ [MPa] |  $\nu$ [/] |  $\sigma_{y0}$ [MPa] | H [MPa]
---|---|---|---|---
Value |  2.09$\times 10^{5}$ | 0.3 | 235 | 800

The time-dependent displacement histories were generated similarly to Section
2.2.1 with six uniformly spaced control points. The applied displacement is 0
at $t=0$. The displacement magnitude at each control point was randomly
selected such that the nominal axial strain magnitude is below 5%. A typical
example of the applied displacement is shown in Fig. 3b. A total of 4000 FE
simulations were generated using Abaqus/Standard [36], and $\bar{\sigma}$ and
$\bar{\epsilon}_{p}$ were stored at 40 uniformly spaced time steps as the
ground truth labels in the NN training. Similar to Section 2.2.1, some data
scaling is needed for best model performance. For this example, the min-max
scaler in Scikit-Learn [37] was used. Two different scalers were used for the
two output components to account for the drastically different scales of
stress and plastic strain.

###  2.3 Inverse history identification with trained S-DeepONet

Once the improved S-DeepONet is trained, it can be used to efficiently infer
time-dependent full-field solution contours. Take the example of predicting
von Mises stress as described in Section 2.2.2. Further post-processing of the
predicted fields yields a plot of the mean von Mises stress (over the entire
dog bone specimen) as a function of load duration, which is a function of the
input load history. If given a known curve of the mean stress over time, it is
possible to infer the load history necessary to generate this stress history
using the trained S-DeepONet and an optimizer. The outputs of this inverse
identification are the five scalar displacement values at the control points,
since the smooth load curve in Section 2.2.2 can be uniquely characterized by
these values (the value at first control point is 0). For this work, the
genetic algorithm (GA) implementation with PyGAD [38] was used, which provides
a gradient-free optimization framework to leverage the trained S-DeepONet as a
black-box model. Using GA, 25 generations of optimization were performed with
a population size of 100 and the number of parents mating was set to 10. The
GA seeks to maximize a scalar fitness function value, and in the current case,
it is defined as the inverse of the mean absolute error (MAE) between the
predicted and known stress histories. The process of evaluating the fitness
value with a trained S-DeepONet is presented in the following algorithm:

Input: Weights and bias of the trained S-DeepONet, point sets $\bm{P}$ with
$N$ points for contour evaluation, five displacement values at control points,
desired output stress history $\bm{\sigma}_{ref}$

1

Output: Fitness value

2

/* Initialization */

3 Load trained weights and bias to the S-DeepONet model $\mathcal{M}$

/* Get current load history */

4 Compute full load curve $\bm{m}$ from control point values with RBI

/* Forward inference */

5 Form a single input set $(\bm{P},\bm{m})$

6 Infer with $\mathcal{M}$ to obtain output
$\bm{O}=\mathcal{M}(\bm{P},\bm{m})$

7 Compute predicted mean stress over time
$\sigma_{pred,s}=\frac{1}{N}\sum_{i=1}^{N}O_{is}$

8

/* Fitness calculation */

9 return fitness $F=MAE(\bm{\sigma}_{ref},\bm{\sigma}_{pred})^{-1}$

Algorithm 1 Evaluation of the fitness value

##  3 Results and discussion

All simulations were conducted with six high-end AMD EPYC 7763 Milan CPU
cores. All NN training and inference were conducted using a single Nvidia A100
GPU card on Delta, an HPC cluster hosted at the National Center for
Supercomputing Applications (NCSA).

To evaluate the model performance in the test set, three quantitative metrics
were used, they are the relative $L_{2}$ error, mean absolute error, and
$R^{2}$ value:

$\displaystyle{\rm{Relative\;L_{2}\;error}}=\frac{|f_{FE}-f_{Pred}|_{2}}{|f_{FE}|_{2}}\times
100\%,$ (13)
$\displaystyle{\rm{MAE}}=\frac{1}{N_{T}}\sum_{i=1}^{N_{T}}\left|f_{FE}-f_{Pred}\right|,$
$\displaystyle
R^{2}=1-\frac{\sum_{i=1}^{N_{T}}\left(f_{FE}-f_{Pred}\right)^{2}}{\sum_{i=1}^{N_{T}}\left(f_{FE}-\bar{f}_{FE}\right)^{2}},$

where $f_{FE}$, $f_{Pred}$, $N_{T}$, $\bar{f}_{FE}$ denote the finite element
(FE) simulated field value, NN-predicted field value, number of test cases,
and mean value of the FE-simulated field values, respectively.

###  3.1 Time-dependent fluid flow

To use the improved S-DeepONet to solve the lid-driven cavity problem, we set
$C=3$ to predict three vector components $P$, $u$ and $v$ and $S=25$ for 25
output time steps. Doing so yields a model with 800384 trainable parameters.
To test the repeatability of the model performance, we trained the S-DeepONet
three times with randomly divided training and testing data with a 80-20 data
split. The results show similar repeatable performance, with a mean (over all
repetitions and components) relative $L_{2}$ error of 5.891% and a standard
deviation of 0.498%. Detailed performance metrics for the best performing run
out of the three are shown in Table 2. The average training and inference (per
case) times over the three runs were 22278s and $8\times 10^{-3}$s,
respectively. While on average, running each FD simulation took 18s, making
the inference 2240 times faster than the direct numerical simulation.

Table 2: Performance metrics for the CFD example Component | Relative $L_{2}$ error [%] | Mean absolute error |  $R^{2}$ value
---|---|---|---
$P$ | 8.244 | 9.671$\times 10^{-3}$ | 0.998
$u$ | 3.465 | 4.060$\times 10^{-3}$ | 0.999
$v$ | 6.717 | 1.198$\times 10^{-3}$ | 0.998

The contour and quiver plots for pressure and velocity vector at different
time steps are shown in Fig. 4. They are ranked by the percentile of mean
relative $L_{2}$ error in the predictions, and the 0th (best case), 80th and
100th (worst case) percentiles are shown for representation. For all plots in
Fig. 4, the filled contour was colored by pressure, and the velocity vectors
are shown as arrows in the domain whose lengths are proportional to the
velocity magnitude.

(a) Step 1, best |  (b) Step 9, best |  (c) Step 17, best |  (d) Step 25, best
---|---|---|---
(e) Step 1, 80th pct |  (f) Step 9, 80th pct |  (g) Step 17, 80th pct |  (h) Step 25, 80th pct
(i) Step 1, worst |  (j) Step 9, worst |  (k) Step 17, worst |  (l) Step 25, worst
Figure 4: Contour and quiver plots for the flow field at different time and
percentiles of prediction accuracy. In each sub-figure, the top and bottom
rows show finite difference solution and S-DeepONet prediction, resp. The
background contour is colored by the pressure, and the flow velocity field is
rendered as arrows with length proportional to the velocity magnitude.

When predicting full-field solution fields at multiple time steps, two types
of errors are worth investigating. They are the time-averaged error (i.e., for
each prediction case, the average error over all prediction time steps) and
case-averaged error (i.e., for each time step, the average error over all
prediction cases). Time-averaged error provides a holistic measurement of
model prediction accuracy over all time steps, while the case-averaged error
provides insight on how the prediction accuracy changes with the magnitude of
the input function. For that, scatter plots of the relative $L_{2}$ prediction
versus the mean magnitude of the lid velocity are shown in Fig. 5.

(a)

(b)

(c)

Figure 5: Scatter plots and trend lines of case-averaged errors for:
LABEL:sub@lp3 Pressure ($P$), LABEL:sub@lp4 x velocity ($u$), LABEL:sub@lp5 y
velocity ($v$).

From the performance metrics in Table 2, we see that the S-DeepONet is able to
generate accurate predictions with relative $L_{2}$ errors less than 10% and
$R^{2}$ of above 0.99 for all components with only 3200 training data point.
The contour plots in Fig. 4 further confirm the accuracy of the model, as we
can see that the S-DeepONet predicted contours match closely with those
computed from direct numerical simulations. For the worst case scenario (last
row of Fig. 4), we see the pressure contour prediction is inaccurate at step
1, where the pressure magnitude is small. Even in this worst case (ranked by
time-averaged error), we see that the model is able to generate accurate
predictions for some time steps (e.g., step 17, third column), where the
magnitudes of the solution fields are reasonable. The relationships between
prediction error and mean lid velocity magnitude is further elucidated by Fig.
5, where we see a consistent decreasing trend for increasing lid velocity
magnitude. All three components in this case show relatively strong
correlations with the trend line $R^{2}\geq 0.7$, indicating that the
prediction error has a strong correlation with the lid velocity, a subject we
will further investigate and compare in Section 3.2.

###  3.2 History-dependent plastic deformation

The means and standard deviations of the metrics over three repetitions of the
S-DeepONet training are summarized in Table 3.

Table 3: Result repeatability, plastic deformation Component | Relative $L_{2}$ error [%] | Mean absolute error |  $R^{2}$ value
---|---|---|---
$\bar{\sigma}$ | 5.164 (8.691$\times 10^{-2}$) 111Data format: Mean (Standard deviation) | 2.544 (6.039$\times 10^{-2}$) MPa | 0.997 (9.764$\times 10^{-5}$)
$\bar{\epsilon}_{p}$ | 21.190 (0.190) | 8.376$\times 10^{-5}$ (9.084$\times 10^{-7}$) | 0.999 (1.279$\times 10^{-4}$)

To show the statistical distribution of prediction error among the 800 test
cases, stress contours that correspond to the 0th (best case), 80th, and 100th
(worst case) percentile prediction error are displayed in Fig. 6 for the
improved S-DeepONet model. For the ranking of prediction errors, the average
relative $L_{2}$ error over all time steps and all output components were used
to obtain a single scalar representation of prediction accuracy over all time
steps and output components. Due to the symmetry in the loading, only the top
half of the NN-predictions are shown in Fig. 6, with the bottom half reserved
for (the flipped top half of) the FE simulation, the two are separated by a
black dashed line.

(a) Step 1, best |  (b) Step 20, best |  (c) Step 40, best
---|---|---
(d) Step 1, 80th pct |  (e) Step 20, 80th pct |  (f) Step 40, 80th pct
(g) Step 1, worst |  (h) Step 20, worst |  (i) Step 40, worst
Figure 6: Contour plots for stress and plastic strain fields at different time
and percentiles of prediction accuracy. In each sub-figure, the top and bottom
rows show stress and plastic strain, resp. In each contour, only top half of
NN-prediction is shown due to symmetry, with FE solution shown in bottom half.

To highlight the efficiency of the current vector-output model, we trained two
scalar S-DeepONets (i.e., by setting $C=1$), one for $\bar{\sigma}$ and the
other for $\bar{\epsilon}_{p}$, and compare the computational efficiency and
accuracy of the models. For a fair comparison, training of the three models
were conducted using identical training and testing data points. The model
size, training, and inference time are compared in Table 4.

Table 4: Vector model vs. scalar models Model | Trainable parameters | Training time [s] | Inference time [s]
---|---|---|---
Vector | 800975 | 15547 | 4.025$\times 10^{-3}$
Scalar, $\bar{\sigma}$ | 797711 | 9424 | 3.675$\times 10^{-3}$
Scalar, $\bar{\epsilon}_{p}$ | 797711 | 9364 | 3.588$\times 10^{-3}$

The prediction quality metrics are summarized in Table 5.

Table 5: Vector model vs. scalar models, performance metrics Model | Rel. $L_{2}$,$\bar{\sigma}$ | Mean abs.,$\bar{\sigma}$ |  $R^{2}$,$\bar{\sigma}$ | Rel. $L_{2}$,$\bar{\epsilon}_{p}$ | Mean abs.,$\bar{\epsilon}_{p}$ |  $R^{2}$,$\bar{\epsilon}_{p}$
---|---|---|---|---|---|---
Vector | 5.164 | 2.544 | 0.997 | 21.190 | 8.376$\times 10^{-5}$ | 0.999
Scalar, $\bar{\sigma}$ | 5.185 | 2.648 | 0.996 | / | / | /
Scalar, $\bar{\epsilon}_{p}$ | / | / | / | 19.471 | 5.549$\times 10^{-5}$ | 0.999

The contour plots predicted by the vector S-DeepONet model and the two scalar
S-DeepONet models at different time steps are shown in Fig. 7 for model
comparison. The test case corresponding to the 50th (median case) percentile
prediction error is shown for representation.

(a) Step 1, 50th pct |  (b) Step 14, 50th pct |  (c) Step 27, 50th pct |  (d) Step 40, 50th pct
---|---|---|---
Figure 7: Contour plots for stress and plastic strain fields predicted by the
vector and scalar S-DeepONets. In each sub-figure, the top and bottom rows
show stress and plastic strain, resp. In each contour, only top half of the
contour is shown due to symmetry. The top half shows the prediction from the
vector S-DeepONet, and the bottom half shows that from the scalar S-DeepONets.

We compare the histograms of the time-averaged error for von Mises stress and
equivalent plastic strain for the vector and scalar S-DeepONets in Fig. 8. The
scatter plots and trend lines for the prediction error versus the mean
absolute magnitude of the applied displacement are shown in Fig. 9.

(a)

(b)

Figure 8: Histograms of time-averaged errors for: LABEL:sub@htr1 von Mises
stress, LABEL:sub@htr2 equivalent plastic strain.

(a)

(b)

Figure 9: Scatter plots and trend lines of case-averaged errors for:
LABEL:sub@lp1 von Mises stress, LABEL:sub@lp2 equivalent plastic strain.

From the repeatability results in Table 3, we see that the performance of the
proposed vector S-DeepONet is consistent and repeatable, with minimal
performance variation when trained with different data. With the von Mises
stress, the model achieves a relative $L_{2}$ error of about 5% and a $R^{2}$
value of 0.997, which is accurate considering that only 3200 data points were
used in training. The relative $L_{2}$ error for equivalent plastic strain is
as high as 21%, but this is expected as this metric is ill-defined before the
specimen yields (all elastic, no plastic strain). When inspecting the mean
absolute error in plastic strain, we see that the error is only about
$8.4\times 10^{-3}$%, and the $R^{2}$ value is 0.999, again indicating
accurate model predictions. It is also worth highlighting that the proposed
architecture solved two problems with completely different physics with only
minor changes in network structures to account for the different number of
time steps and output vector components, thereby demonstrating the
generalizability and versatility of the proposed framework. The contour plots
at different percentiles and time steps as shown in Fig. 6 provide a more
direct view of the prediction performance. In the best case (first row of Fig.
6), the S-DeepONet is able to accurately capture the hot spots in stress and
plastic strains at different time steps, showing close agreement with the
corresponding FE simulations. At the 80th percentile (second row of Fig. 6),
we see that the S-DeepONet was unable to predict the stress contour initially
when the stress is small and response is fully elastic, but is still able to
provide accurate predictions of both quantities at later time steps once the
magnitudes of the stress and plastic strain increase. A similar trend is
observed for the worst case (last row of Fig. 6), where the S-DeepONet is
unable to capture the stress distribution due to its small magnitude. Note
that the specimen remains completely elastic throughout this loading history
as evidenced from the reference FE simulation. The relation of prediction
error and solution magnitude will be further explored later in the discussion.

The key improvement of the current architecture is the ability to predict
multiple vector components at time steps with one model. For the case of two
output components, we see from Table 4 that the vector model only has 0.4%
more trainable parameters as compared to the scalar S-DeepONet model. Training
the vector model took 15547s, which is 20.8% faster than the combined training
time for the two scalar models, while the inference speed is about 80% faster.
On average, each FE simulation requires 48s of CPU time to solve, making the
S-DeepONet inference about 11900 times faster than running a FE simulation.
Moreover, as revealed by the performance metrics presented in Table 5, we see
that the performance of the vector S-DeepONet is highly similar to (albeit
slightly worse) than the individual scalar networks. From Fig. 7, we see that
the predicted contour plots of the von Mises stress and equivalent plastic
strain are highly similar at different time steps as well. Except at the
beginning of the loading history (first column of Fig. 7) where the specimen
is fully elastic, the vector and scalar models predicted different contours
for equivalent plastic strain (when they should be identically 0). Fig. 8
compares the histograms of the time-averaged errors from the three different
models. We see that the distribution of the prediction accuracy is similar
between the vector and scalar S-DeepONet models. However, considering the
significant time and model size savings, training the proposed vector
S-DeepONet model instead of training two scalar models individually is the
superior option from an efficiency perspective.

Lastly, for a multi-step prediction, it is worth studying how the prediction
error changes as a function of the mean magnitude of the applied displacement
(which changes with time). From Fig. 9, we see two slightly different
behavior. For the von Mises stress (Fig. 9a), the negative slope of the fitted
trend line indicates that the error is, in general, decreasing with increasing
load magnitude. However, the small $R^{2}$ value of the trend line indicates
that the correlation is minimal. For equivalent plastic strain (Fig. 9b), a
similar decreasing trend is observed with increasing load magnitude, but this
time with a much higher $R^{2}$ value. As a comparison, recall from Section
2.2.1 that all three components ($P,u$, and $v$) show a relative strong
correlation with trend line $R^{2}\geq 0.7$. This difference is reasonable
since plasticity is path-dependent; the current stress magnitude depends not
only on the current displacement magnitude but also on the past loading
history, hence the lower correlation with the current displacement magnitude.
On the other hand, the steady-state flow field can be uniquely defined by the
lid velocity. Hence, the time-dependent flow field magnitude shows a stronger
correlation with the current lid velocity. The limitation that the model
cannot accurately predict the contour when the magnitude of the field
variables is small is likely due to the use of the MSE loss function, where
those data points only account for a small portion to the total loss and
therefore not effectively minimized during training.

###  3.3 Application of the trained S-DeepONet model

Three additional load curves were randomly generated following the same
procedure outlined in Section 2.2.2. Those curves were not seen by the trained
S-DeepONet model during its training and testing. FE simulations were
performed on those load curves to obtain the reference stress history curves,
which were inputs to the inverse identification. To validate the results, FE
simulations were performed using the load paths identified by S-DeepONet. The
comparison of the load curves and stress histories are shown in Fig. 10.

(a)

(b)

(c)

Figure 10: Comparison of the load curves and stress histories identified by
S-DeepONet and GA.

With the help of efficient forward inference of the S-DeepONet model, the
three GA cases were completed in an average of 24s, which is fast considering
that a single forward FE simulation of the dog bone specimen takes around 48s.
From the results, we see that in all three cases, the S-DeepONet predicted
stress histories match closely with the given ground truths. Moreover, a
comparison of the S-DeepONet predictions and the corresponding FE validation
results (generated from the same load curve predicted by GA) shows that the
S-DeepONet predictions remain highly accurate for load curves outside of the
original training and testing data points, again showing the high accuracy of
the trained model. Despite the high similarity between the identified and
given stress histories, the load curves identified by GA and S-DeepONet do not
always match the known ground truths. We see two characteristic cases from
Fig. 10: (1) matching load curves (i.e., Fig. 10a), (2) load curves with equal
but opposite displacements (i.e., Fig. 10b). It is also interesting to see a
combination of both characteristic behaviors in Fig. 10c where the second half
of the predicted load curve matches closely with ground truth while the first
half shows opposite displacements. This is reasonable from a mechanics
perspective, since the von Mises stress is a positive quantity regardless of
tensile or compressive loading, and the material considered in this work does
not exhibit tension-compression asymmetry.

##  4 Conclusions, limitations, and future work

The sequential DeepONet (S-DeepONet) model previously proposed by the authors
[31] is a variant of the DeepONet architecture [13] that uses a gated
recurrent unit (GRU) in the branch network to capture the temporal information
in the time-dependent input functions. However, the original S-DeepONet
architecture only predicts the last time frame of a time-dependent evolution,
and the predicted field is a scalar. Realizing this limitation, we introduced
an improved version of the S-DeepONet capable of simultaneously predicting
multiple time steps for solution fields with multiple vector components. This
feature is made possible by the tensor-product structure of the combination
operation that combines the encoded temporal information from the branch
network and the encoded spatial information (for all vector components) from
the trunk network. We further elucidate that the tensor-product combination
process can be viewed as the simultaneous identification of a set of basis
shapes and weights, with which the final contour can be expressed as a
weighted linear combination. This architecture expands on the idea of
exploiting the powerful temporal encoding capability of the GRU and the
spatial encoding capability of the DeepONet architecture. To demonstrate the
application of the improved S-DeepONet, we showed an example of lid-driven
cavity flow and an example of plastic deformation, both subjected to time-
dependent input loads and having multiple output components at many time
steps. In both cases, the improved S-DeepONet was able to provide accurate
predictions for all output components using only 3200 data points in training.
For all components, the DeepONet predictions achieved a $R^{2}$ value of over
0.99 and a relative $L_{2}$ error of less than 10% (except for the equivalent
plastic strain in Section 3.2). We highlight the fact that minimum
architecture change is needed (except for change the output components and
time steps) for solving those two problems with totally different underlying
physics, which shows that the proposed S-DeepONet is highly versatile. Once
trained, the S-DeepONet can infer accurate full-field results at different
time steps at least three orders of magnitude faster than direct numerical
simulations. The trained model can also be used in conjunction with gradient-
free optimizers such as the genetic algorithm to perform accurate inverse
parameter identification. The example using the plasticity data shown that the
inverse identification is highly efficient (finished in about half the time of
one FE simulation) and accurate. Using the plasticity example, we also
demonstrated the effectiveness of predicting all output components at once
using a vector network instead of training a separate scalar network for each
of the output component. With only 0.4% more parameters, the vector S-DeepONet
trained 20.8% faster than the two scalar networks combined, while maintaining
a very similar level of prediction accuracy. Therefore, it is recommended to
use the vector S-DeepONet proposed in the current work, if multiple output
components are of interest at multiple time steps.

Through the analysis of the mean error magnitude as a function of the mean
input load magnitude, it was revealed that the current S-DeepONet suffers from
the limitation that it is unable to accurately predict the field contours when
the field magnitude is small, even with a time-dependent data scaling scheme.
This is likely due to the use of the MSE loss function, which only results in
a small loss contribution for those data points with small magnitudes.

With the capability to accurately and efficiently predict solution history at
multiple time steps for multiple different components, the improved S-DeepONet
architecture proposed in this work provides a versatile tool to the
engineering community to build surrogate models for complex nonlinear
numerical simulations. The versatility of this architecture is demonstrated in
this work through its application in both solid and fluid mechanics problems,
and can be further employed in different engineering applications. With the
weights and bias of the trained model, almost instant full-field forward
predictions can be evaluated even on low-end platforms like laptops, which
provides a novel tool for rapid preliminary designs, sensitivity analysis,
uncertainty quantification, online controls, and as a black-box surrogate
model for optimization, as we have demonstrated. In future work, we will
investigate the use of transformer models [39] and attention mechanisms [40]
in the DeepONet architecture to achieve higher prediction accuracy with lower
computational costs.

## Replication of results

The data and source code that support the findings of this study can be found
at https://github.com/Jasiuk-Research-Group. Note to editor and reviewers: the
link above will be made public upon the publication of this manuscript. During
the review period, the data and source code can be made available upon request
to the corresponding author.

## Conflict of interest

The authors declare that they have no conflict of interest.

## Acknowledgements

The authors would like to thank the National Center for Supercomputing
Applications (NCSA) at the University of Illinois, and particularly its
Research Consulting Directorate, the Industry Program, and the Center for
Artificial Intelligence Innovation (CAII) for their support and hardware
resources. This research is a part of the Delta research computing project,
which is supported by the National Science Foundation (award OCI 2005572) and
the State of Illinois, as well as the Illinois Computes program supported by
the University of Illinois Urbana-Champaign and the University of Illinois
System. Finally, the authors would like to thank Professors George
Karniadakis, Lu Lu, and the Crunch team at Brown, whose original work with
DeepONets inspired this research.

## CRediT author contributions

Junyan He: Methodology, Formal analysis, Investigation, Writing - Original
Draft. Shashank Kushwaha: Methodology, Formal analysis, Investigation, Writing
- Original Draft. Jaewan Park: Methodology, Investigation, Writing - Original
Draft. Seid Koric: Conceptualization, Methodology, Supervision, Resources,
Writing - Original Draft, Funding Acquisition. Diab Abueidda: Supervision,
Writing - Review & Editing. Iwona Jasiuk: Supervision, Writing - Review &
Editing.

## References

  * Kushwaha et al. [2023] Shashank Kushwaha, Junyan He, Diab Abueidda, and Iwona Jasiuk.  Designing impact-resistant bio-inspired low-porosity structures using neural networks.  _Journal of Materials Research and Technology_ , 2023.
  * He et al. [2023a] Junyan He, Shashank Kushwaha, Diab Abueidda, and Iwona Jasiuk.  Exploring the structure-property relations of thin-walled, 2d extruded lattices using neural networks.  _Computers & Structures_, 277:106940, 2023a.
  * Perumal et al. [2023] V Perumal, D Abueidda, S Koric, and A Kontsos.  Temporal convolutional networks for data-driven thermal modeling of directed energy deposition.  _Journal of Manufacturing Processes_ , 85:405–416, 2023.
  * Chen and Gu [2020] Chun-Teh Chen and Grace X Gu.  Generative deep neural networks for inverse materials design using backpropagation and active learning.  _Advanced Science_ , 7(5):1902607, 2020.
  * Fuhg et al. [2022] Jan N Fuhg, Michele Marino, and Nikolaos Bouklas.  Local approximate gaussian process regression for data-driven constitutive models: development and comparison with neural networks.  _Computer Methods in Applied Mechanics and Engineering_ , 388:114217, 2022.
  * Mandl et al. [2023] Luis Mandl, André Mielke, Seyed Morteza Seyedpour, and Tim Ricken.  Affine transformations accelerate the training of physics-informed neural networks of a one-dimensional consolidation problem.  _Scientific Reports_ , 13(1):15566, 2023.
  * Fuhg and Bouklas [2022] Jan N Fuhg and Nikolaos Bouklas.  The mixed deep energy method for resolving concentration features in finite strain hyperelasticity.  _Journal of Computational Physics_ , 451:110839, 2022.
  * He et al. [2023b] Junyan He, Diab Abueidda, Rashid Abu Al-Rub, Seid Koric, and Iwona Jasiuk.  A deep learning energy-based method for classical elastoplasticity.  _International Journal of Plasticity_ , 162:103531, 2023b.
  * Lambora et al. [2019] Annu Lambora, Kunal Gupta, and Kriti Chopra.  Genetic algorithm-a literature review.  In _2019 international conference on machine learning, big data, cloud and parallel computing (COMITCon)_ , pages 380–384. IEEE, 2019.
  * Corriveau et al. [2010] Guillaume Corriveau, Raynald Guilbault, and Antoine Tahan.  Genetic algorithms and finite element coupling for mechanical optimization.  _Advances in Engineering Software_ , 41(3):422–426, 2010.
  * Abueidda et al. [2019] Diab W Abueidda, Mohammad Almasri, Rami Ammourah, Umberto Ravaioli, Iwona M Jasiuk, and Nahil A Sobh.  Prediction and optimization of mechanical properties of composites using convolutional neural networks.  _Composite Structures_ , 227:111264, 2019.
  * Wang et al. [2022] Yongzhen Wang, Qinglei Zeng, Jizhen Wang, Ying Li, and Daining Fang.  Inverse design of shell-based mechanical metamaterial with customized loading curves based on machine learning and genetic algorithm.  _Computer Methods in Applied Mechanics and Engineering_ , 401:115571, 2022.
  * Lu et al. [2021a] Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  _Nature machine intelligence_ , 3(3):218–229, 2021a.
  * Koric et al. [2023] Seid Koric, Asha Viswantah, Diab W Abueidda, Nahil A Sobh, and Kamran Khan.  Deep learning operator network for plastic deformation with variable loads and material properties.  _Engineering with Computers_ , pages 1–13, 2023.
  * Lu et al. [2022] Lu Lu, Raphaël Pestourie, Steven G Johnson, and Giuseppe Romano.  Multifidelity deep neural operators for efficient learning of partial differential equations with application to fast inverse design of nanoscale heat transport.  _Physical Review Research_ , 4(2):023210, 2022.
  * He et al. [2023c] Junyan He, Seid Koric, Shashank Kushwaha, Jaewan Park, Diab Abueidda, and Iwona Jasiuk.  Novel deeponet architecture to predict stresses in elastoplastic structures with variable complex geometries and loads.  _Computer Methods in Applied Mechanics and Engineering_ , 415:116277, 2023c.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2023.116277.  URL https://www.sciencedirect.com/science/article/pii/S0045782523004012.
  * Haghighat et al. [2023] Ehsan Haghighat, Umair bin Waheed, and George Karniadakis.  A novel deeponet model for learning moving-solution operators with applications to earthquake hypocenter localization.  _arXiv preprint arXiv:2306.04096_ , 2023.
  * Koric and Abueidda [2023] Seid Koric and Diab W Abueidda.  Data-driven and physics-informed deep learning operators for solution of heat conduction equation with parametric heat source.  _International Journal of Heat and Mass Transfer_ , 203:123809, 2023.
  * Li et al. [2023] Wei Li, Martin Z Bazant, and Juner Zhu.  Phase-field deeponet: Physics-informed deep operator neural network for fast simulations of pattern formation governed by gradient flows of free-energy functionals.  _arXiv preprint arXiv:2302.13368_ , 2023.
  * Nie et al. [2020] Zhenguo Nie, Haoliang Jiang, and Levent Burak Kara.  Stress field prediction in cantilevered structures using convolutional neural networks.  _Journal of Computing and Information Science in Engineering_ , 20(1):011002, 2020.
  * Hoq et al. [2023] Enjamamul Hoq, Osama Aljarrah, Jun Li, Jing Bi, Alfa Heryudono, and Wenzhen Huang.  Data-driven methods for stress field predictions in random heterogeneous materials.  _Engineering Applications of Artificial Intelligence_ , 123:106267, 2023.
  * Bhatnagar et al. [2019] Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik.  Prediction of aerodynamic flow fields using convolutional neural networks.  _Computational Mechanics_ , 64:525–545, 2019.
  * Chen et al. [2020] Donglin Chen, Xiang Gao, Chuanfu Xu, Shizhao Chen, Jianbin Fang, Zhenghua Wang, and Zheng Wang.  Flowgan: A conditional generative adversarial network for flow prediction in various conditions.  In _2020 IEEE 32nd International Conference on Tools with Artificial Intelligence (ICTAI)_ , pages 315–322. IEEE, 2020.
  * Schmidhuber et al. [1997] Jürgen Schmidhuber, Sepp Hochreiter, et al.  Long short-term memory.  _Neural Comput_ , 9(8):1735–1780, 1997.
  * Cho et al. [2014] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio.  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  _arXiv preprint arXiv:1406.1078_ , 2014.
  * Hochreiter [1998] Sepp Hochreiter.  The vanishing gradient problem during learning recurrent neural nets and problem solutions.  _International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems_ , 6(02):107–116, 1998.
  * Abueidda et al. [2021] Diab W Abueidda, Seid Koric, Nahil A Sobh, and Huseyin Sehitoglu.  Deep learning for plasticity and thermo-viscoplasticity.  _International Journal of Plasticity_ , 136:102852, 2021.
  * Deng et al. [2019] Zhiwen Deng, Yujia Chen, Yingzheng Liu, and Kyung Chun Kim.  Time-resolved turbulent velocity field reconstruction using a long short-term memory (lstm)-based artificial intelligence framework.  _Physics of Fluids_ , 31(7), 2019.
  * Hassanian et al. [2022] Reza Hassanian, Ásdís Helgadóttir, and Morris Riedel.  Deep learning forecasts a strained turbulent flow velocity field in temporal lagrangian framework: Comparison of lstm and gru.  _Fluids_ , 7(11):344, 2022.
  * Jin et al. [2020] Xiaowei Jin, Shujin Laima, Wen-Li Chen, and Hui Li.  Time-resolved reconstruction of flow field around a circular cylinder by recurrent neural networks based on non-time-resolved particle image velocimetry measurements.  _Experiments in Fluids_ , 61:1–23, 2020.
  * He et al. [2024] Junyan He, Shashank Kushwaha, Jaewan Park, Seid Koric, Diab Abueidda, and Iwona Jasiuk.  Sequential deep operator networks (s-deeponet) for predicting full-field solutions under time-dependent loads.  _Engineering Applications of Artificial Intelligence_ , 127:107258, 2024.  ISSN 0952-1976.
  * Lu et al. [2021b] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis.  Deepxde: A deep learning library for solving differential equations.  _SIAM review_ , 63(1):208–228, 2021b.
  * Abadi et al. [2015] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015.  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Barba and Forsyth [2018] Lorena A Barba and Gilbert F Forsyth.  Cfd python: the 12 steps to navier-stokes equations.  _Journal of Open Source Education_ , 2(16):21, 2018.
  * SIMULIA [2020] SIMULIA.  Abaqus, 2020.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011.
  * Gad [2021] Ahmed Fawzy Gad.  Pygad: An intuitive genetic algorithm python library.  _arXiv preprint arXiv:2106.06158_ , 2021.
  * Jaderberg et al. [2015] Max Jaderberg, Karen Simonyan, Andrew Zisserman, et al.  Spatial transformer networks.  _Advances in neural information processing systems_ , 28, 2015.
  * Vaswani et al. [2017] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  _Advances in neural information processing systems_ , 30, 2017.


# Gravity with dynamical torsion

M. O. Katanaev
Steklov Mathematical Institute of Russian Academy of Sciences,
8 Gubkina St., Moscow 119991, Russia E-mail<EMAIL_ADDRESS>
###### Abstract

We propose four simple Lagrangians for gravity models with dynamical torsion
which are free from ghosts and tachyons. The torsion propagates two massive or
massless particles of spin $1^{\pm}$ and $0^{\pm}$ besides the massless
graviton $2^{+}$ propagated by metric.

##  1 Introduction

Introduction of dynamical torsion in general relativity is the natural
geometric generalization of Einstein’s gravity. It is the subject of research
for many years. Physical motivations and earlier references can be found in
[1, 2, 3, 4, 5]. In spite of much efforts, we still do not have a generally
acknowledged action for such model. It is usually assumed that gravity
interactions are described by metric propagating spin $2^{+}$ massless
gravitons, whereas torsion describes massive or massless particles of spin
$2^{\pm}$, $1^{\pm}$, or $0^{\pm}$. To answer the crucial question whether
torsion can yield a physically acceptable generalization, we need a model but
it is missing, rather we have too many proposals and cannot make a unique
choice. To my mind, the main problem here is the mathematical complexity.

Torsion has three tensor indices and can propagate several massive or massless
particles of spin not higher then 2. From physical standpoint, the model must
not include ghosts and tachyons. A number of such models were proposed forty
years ago [6, 7, 8, 9], but they are not analyzed until now if full detail,
mainly, because of technical difficulties. Moreover, these models were found
within a certain class of Lagrangians which do not exhaust all possibilities.
In the present paper, we consider another class of dynamical torsion models
without higher derivatives, considering metric and torsion as independent
variables, and found four Lagrangians without ghosts and tachyons. The final
answer seems to be much simpler and transparent then previous proposals.

##  2 General Lagrangian

We consider a space-time manifold ${\mathbb{M}}$ with local coordinates
$x^{\alpha}$, $\alpha=0,1,2,3$, and the Riemann–Cartan geometry. The geometry
can be described either by metric $g_{\alpha\beta}$ and torsion
$T_{\alpha\beta}{}^{\gamma}=-T_{\beta\alpha}{}^{\gamma}$ or, equivalently, by
Cartan variables: vierbein $e_{\alpha}{}^{a}$ and Lorentz connection
$\omega_{\alpha}{}^{ab}=-\omega_{\alpha}{}^{ba}$, $a,b=0,1,2,3$.

As usual, a vierbein defines unique metric

$g_{\alpha\beta}:=e_{\alpha}{}^{a}e_{\beta}{}^{b}\eta_{ab},\qquad\eta_{ab}:={\sf\,diag\,}(+---),$

and the Lorentz connection is defined through the affine connection
$\Gamma_{\alpha\beta}{}^{\gamma}$ by the vielbein postulate

$\nabla_{\alpha}e_{\beta}{}^{a}=\partial_{\alpha}e_{\beta}{}^{a}-\Gamma_{\alpha\beta}{}^{\gamma}e_{\gamma}{}^{a}+e_{\beta}{}^{b}\omega_{\alpha
b}{}^{a}=0.$ (1)

If the vierbein is given, then this equation provides the one-to-one
correspondence between Lorentz connection and the torsion
$T_{\alpha\beta}{}^{\gamma}:=\Gamma_{\alpha\beta}{}^{\gamma}-\Gamma_{\beta\alpha}{}^{\gamma}$.

The vierbein has 6 components more then the metric due to the extra symmetry
with respect to local Lorentz transformations acting on Latin indices
$a,b,\dotsc$. In what follows, transformation of Latin indices into Greek ones
and vice versa is always performed by the vierbein and its inverse
$e^{\alpha}{}_{a}$. Moreover we assume that coordinates are chosen in such a
way that the coordinate $x^{0}$ is timelike and all sections
$x^{0}={\sf\,const}$ are spacelike.

Torsion and curvature have the following expressions in terms of Cartan
variables

$\displaystyle T_{\alpha\beta}{}^{a}$
$\displaystyle=\partial_{\alpha}e_{\beta}{}^{a}-e_{\alpha}{}^{b}\omega_{\beta
b}{}^{a}-(\alpha\leftrightarrow\beta),$ (2) $\displaystyle R_{\alpha\beta
a}{}^{b}$ $\displaystyle=\partial_{\alpha}\omega_{\beta
a}{}^{b}-\omega_{\alpha a}{}^{c}\omega_{\beta
c}{}^{b}-(\alpha\leftrightarrow\beta).$ (3)

If we want vierbein and Lorentz connection to be really independent
propagating fields then they must satisfy differential equations of motion.
Let us require the field equations to be covariant and of second order. Then
the Lagrangian should be quadratic in curvature and torsion. The general form
of such Lagrangian is

$\displaystyle\frac{1}{\sqrt{|g|}}L=$ $\displaystyle\kappa
R-\frac{1}{4}T^{abc}\big{(}\beta_{1}T_{abc}+\beta_{2}T_{cab}+\beta_{3}\eta_{ac}T_{b}\big{)}-$
$\displaystyle-\frac{1}{4}R^{abcd}\big{(}\gamma_{1}R_{abcd}+\gamma_{2}R_{cdab}+\gamma_{3}R_{acbd}+\gamma_{4}\eta_{bd}R_{ac}+\gamma_{5}\eta_{bd}R_{ca}\big{)}-2\Lambda,$
(4)

where $T_{b}:=T_{ab}{}^{b}$ is the trace of torsion, $R_{ab}:=R_{acb}{}^{c}$
and $R:=R_{a}{}^{a}$ are the Ricci tensor and scalar curvature, respectively;
$\kappa$, $\beta_{1,2,3}$, and $\gamma_{1,2,3,4,5}$ are coupling constants,
and $\Lambda$ is a cosmological constant. We do not have to add the scalar
curvature $\widetilde{R}(e)$ built from Christoffel’s symbols for zero torsion
due to the identity

$R+\frac{1}{4}T_{\alpha\beta\gamma}T^{\alpha\beta\gamma}-\frac{1}{2}T_{\alpha\beta\gamma}T^{\gamma\alpha\beta}-T_{\alpha}T^{\alpha}-\frac{2}{\sqrt{|g|}}\partial_{\alpha}\big{(}\sqrt{|g|}T^{\alpha}\big{)}=\widetilde{R}.$
(5)

Here and in what follows the tilde sign marks geometric notions constructed
entirely from metric for zero torsion.

There are three independent invariants quadratic in torsion in Lagrangian (4)
because torsion has three irreducible components with respect to the Lorentz
group. Curvature tensor in Riemann–Cartan geometry has six irreducible
components but we added only five in Eq.(4) due to the Gauss–Bonnet type
identity

$R_{abcd}R^{cdab}-4R_{ab}R^{ba}+R^{2}\equiv\widetilde{\nabla}_{\alpha}X^{\alpha},$
(6)

where $X^{\alpha}$ are components of some vector.

Several models without ghosts and tachyons were found within the class of
Lagrangians (4) in [7, 8, 9], and the analysis is rather complicated. There is
the additional problem. Lagrangian (4) is not the most general one in
dynamical torsion models, because the existence or absence of higher
derivatives depend on the chosen independent variables. For example,
transformation of the Lorentz connection into torsion tensor by Eq.(1)
involves derivatives of the vierbein. Therefore Lagrangian (4) is a higher
derivative one for the vierbein-torsion variables.

In what follows, we need the decomposition of torsion tensor $T_{abc}$ into
irreducible pieces with respect to the Lorentz group. To this end we extract
the trace $T_{a}$ and totally antisymmetric part $T_{[abc]}$:

$T_{abc}=S_{abc}+\frac{1}{3}(\eta_{ac}T_{b}-\eta_{bc}T_{a})+T_{[abc]},$ (7)

where

$T_{b}:=T_{ab}{}^{a},\qquad T_{[abc]}:=\frac{1}{3}(T_{abc}+T_{cab}+T_{bca}).$

The irreducible component $S_{abc}$ is defined by the algebraic relations:

$S_{abc}=-S_{bac},\qquad S_{ab}{}^{a}=0,\qquad S_{abc}+S_{bca}+S_{cab}=0.$ (8)

In four dimensions, the totally antisymmetric component of torsion is
parameterized by the pseudovector field
$T^{*d}:=-\frac{1}{6}T_{abc}\varepsilon^{abcd}$, where $\varepsilon^{abcd}$ is
the totally antisymmetric tensor, $\varepsilon_{0123}=1$. Thus the torsion
tensor has three irreducible components.

We change variables: $e_{\alpha}{}^{a},\omega_{\alpha}{}^{ab}$ $\mapsto$
$e_{\alpha}{}^{a},T_{\alpha\beta\gamma}$. Now the quadratic torsion terms
become simply the mass terms for torsion which can be easily analyzed. The
transformation of variables restricts the possible choice of kinetic terms for
the vierbein, if we do not want it to be mixed with the torsion. The
corresponding invariant must have dimension $l^{-2}$, where $l$ is the
dimension of length. We have only five tensors of this dimensionality in our
disposal

$\widetilde{R}_{\alpha\beta\gamma\delta},\qquad
R_{\alpha\beta\gamma\delta},\qquad\widetilde{\nabla}_{\alpha}T_{\beta\gamma\delta},\qquad\nabla_{\alpha}T_{\beta\gamma\delta},\qquad
T_{\alpha\beta\gamma}T_{\delta\epsilon\zeta}.$

The invariants are built by different contractions of these tensors with the
metric. We see that the only possible invariants are $\widetilde{R}$ and
$\widetilde{\nabla}_{\alpha}T^{\alpha}$. The last term mixes kinetic term of
the vierbein with torsion components. To simplify matters we postulate that
this term must be absent in dynamical torsion theory. Thus we are left with
the unique possibility: the Hilbert–Einstein action for the vierbein.
Therefore, we further transform variables: $e_{\alpha}{}^{a}$,
$T_{\alpha\beta\gamma}$ $\mapsto$ $g_{\alpha\beta}$, $T_{\alpha\beta\gamma}$.

It is well known that metric in the Hilbert–Einstein action describes massless
gravitons $2^{+}$, and the full symmetry under general coordinate
transformations must be used for elimination of the unphysical degrees of
freedom of the metric. Therefore the problem is reduced to choosing the
Lagrangian for the torsion tensor, and, what is important, we cannot use
general coordinate transformations to eliminate its degrees of freedom.

##  3 (1+3)-decomposition of torsion

Torsion tensor has three irreducible components with respect to the Lorentz
group: $S_{abc}$, $T_{a}$, and $T^{*}_{a}$. To analyse positive definiteness
of quadratic forms of velocities and masses, we further decompose torsion
components with respect to the rotational subgroup
${\mathbb{S}}{\mathbb{O}}(3)\subset{\mathbb{S}}{\mathbb{O}}(1,3)$, because
invariant quadratic forms for the Lorentz group are not positive or negative
definite.

Positive definiteness of the quadratic form of torsion momenta or,
equivalently, velocities depends only on the metric signature. Hence, we put
$g_{\alpha\beta}\mapsto\eta_{ab}:={\sf\,diag\,}(+---)$, i.e. we shall work in
Minkowskian space. We perform the (3+1)-decomposition of all components:
$(a)=(0,i)$, $i:=1,2,3$. Space indices are denoted by Latin letters from the
middle of the alphabet $i,j,\dotsc:=1,2,3$, and summation over them is
performed by the Euclidean metric $\delta_{ij}$. For example,
$T^{a}T_{a}=T^{0}T_{0}-T^{i}T_{i}$. Covector components of torsion tensor have
simple decompositions $(T_{a}):=(T_{0},T_{i})$ and
$(T^{*}_{a}):=(T^{*}_{0},T^{*}_{i})$. In the massive case, components $T_{0}$
and $T^{*}_{0}$ describe scalar particles $0^{+}$ and $0^{-}$, respectively,
where indices denote parities. Covector components $(T_{i})$ and $(T^{*}_{i})$
correspond to vector particles of spin $1^{-}$ and $1^{+}$. Irreducible
component $S_{abc}$ describes more particles. Because of the antisymmetry in
the first two indices the following representations hold:

$S_{ij0}=\varepsilon_{ijl}H^{l}\quad\Leftrightarrow\quad
H^{l}:=\frac{1}{2}S_{ij0}\varepsilon^{ijl};\qquad
S_{ijk}=\varepsilon_{ijl}H^{l}{}_{k}\quad\Leftrightarrow\quad
H^{l}{}_{k}:=\frac{1}{2}S_{ijk}\varepsilon^{ijl},$ (9)

where $\varepsilon_{ijk}$ is the totally antisymmetric third rank tensor,
$\varepsilon_{123}=1$. The relation $S_{ijk}\varepsilon^{ijk}=0$ implies the
tracelessness condition $H^{i}{}_{i}=0$. Components $S_{0ij}$ are decomposed
into symmetric traceless, antisymmetric parts and the trace:

$S_{0ij}=S_{0\\{ij\\}}+S_{0[ij]}+\frac{1}{3}\delta_{ij}S,\qquad
S:=S_{0i}{}^{i}.$

Equality $S_{a0}{}^{a}=S_{00}{}^{0}+S_{i0}{}^{i}=0$ implies $S=0$. From the
relation $S_{ai}{}^{a}=S_{0i}{}^{0}+S_{ji}{}^{j}=0$ we deduce that

$S_{0i0}=-\varepsilon_{ijk}H^{jk}\quad\Leftrightarrow\quad
H^{[mn]}=-\frac{1}{2}S_{0i0}\varepsilon^{imn}.$ (10)

Finally, the equality $S_{0ij}+S_{ij0}+S_{j0i}=0$ yields the expression for
the antisymmetric part

$S_{0[ij]}=-\frac{1}{2}\varepsilon_{ijk}H^{k}\quad\Leftrightarrow\quad
H^{k}=-S_{0ij}\varepsilon^{ijk}.$ (11)

Thus, irreducible components of torsion with respect to the subgroup
${\mathbb{S}}{\mathbb{O}}(3)\subset{\mathbb{S}}{\mathbb{O}}(1,3)$ are:

$\displaystyle S_{0i0},$ $\displaystyle(1^{-})$ $\displaystyle-\quad\text{3
components},$ (12) $\displaystyle S_{0\\{ij\\}},\quad S_{0i}{}^{i}=0,$
$\displaystyle(2^{+})$ $\displaystyle-\quad\text{5 components},$
$\displaystyle H_{\\{ij\\}},\quad H_{i}{}^{i}=0,$ $\displaystyle(2^{-})$
$\displaystyle-\quad\text{5 components},$ $\displaystyle H_{i},$
$\displaystyle(1^{+})$ $\displaystyle-\quad\text{3 components},$
$\displaystyle T_{0},$ $\displaystyle(0^{+})$ $\displaystyle-\quad\text{1
component},$ $\displaystyle T_{i},$ $\displaystyle(1^{-})$
$\displaystyle-\quad\text{3 components},$ $\displaystyle T^{*}_{0},$
$\displaystyle(0^{-})$ $\displaystyle-\quad\text{1 component},$ $\displaystyle
T^{*}_{i},$ $\displaystyle(1^{+})$ $\displaystyle-\quad\text{3 components},$

the antisymmetric parts $S_{0[ij]}$ and $H_{[ij]}$ being expressed through
independent components by equalities (11) and (10).

In the massive case, all independent ${\mathbb{S}}{\mathbb{O}}(3)$-components
of torsion are propagated. If one or more components of torsion turn out to be
massless, then the model posses extra gauge symmetry and the number of
propagated components is reduced but positive definiteness of the velocity
quadratic form does not change. Therefore we focus our attention on the
massive case.

##  4 Torsion mass terms

Since $T_{abc}$ has three Lorentz-irreducible components, $S_{abc}$, $T_{a}$,
and $T^{*}_{a}$, then the most general mass term depends on three constants:

$L_{m}:=-\beta_{1}S^{abc}S_{abc}-\beta_{2}T^{a}T_{a}-\beta_{3}T^{*a}T^{*}_{a}.$
(13)

If there are no tachyons, then this quadratic form must be negative definite.
Simple calculations yield formulae:

$\begin{split}S^{abc}S_{abc}=&2S^{0\\{ij\\}}S_{0\\{ij\\}}-2H^{\\{ij\\}}H_{\\{ij\\}}-3S^{0i0}S_{0i0}+3H^{i}H_{i},\\\
T^{a}T_{a}=&(T_{0})^{2}-T^{i}T_{i},\\\
T^{*a}T^{*}_{a}=&(T^{*}_{0})^{2}-T^{*i}T^{*}_{i}.\end{split}$ (14)

We see that for any choice of $\beta_{1,2,3}$ the quadratic form (13) is not
negative definite. It means that those torsion components which produce
positive contribution must be nonpropagating and have no kinetic terms.

Suppose that $\beta_{1}>0$, then the negative contribution is given by
$S_{0\\{ij\\}}$ and $H_{i}$. Thus, if the kinetic term is correct, the
component $S_{abc}$ can describe only massive particles $2^{+}$ and $1^{+}$.
Components $H_{\\{ij\\}}$ and $S_{0i0}$ produce positive contribution,
correspond to tachyons, and must have no velocity squared terms in the
Lagrangian. Conversely, if $\beta_{1}<0$, then propagating modes of torsion
component $S_{abc}$ can be only $H_{\\{ij\\}}$ and $S_{0i0}$. In this case,
$S_{abc}$ describes only particles $2^{-}$ and $1^{-}$.

Vector $T_{a}$ and pseudovector $T^{*}_{a}$ components of torsion describe
particles of positive mass in the following cases. If $\beta_{2}>0$, then
propagating mode can be only $T_{0}$, which correspond to scalar particles
$0^{+}$. For $\beta_{2}<0$ the negative contribution is given by vector
component $T_{i}$ describing $1^{-}$ particles.

The similar situation happens for $T^{*}_{a}$ component. If $\beta_{3}>0$ or
$\beta_{3}<0$, then propagating particles can be only $T^{*}_{0}$ $(0^{-})$ or
$T^{*}_{i}$ $(1^{+})$, respectively.

##  5 Kinetic terms

The kinetic term has dimensionality $l^{-4}$ and is built from tensors [10]:

$\displaystyle\nabla_{\alpha}\nabla_{\beta}R_{\gamma\delta\epsilon\zeta},$
$\displaystyle\quad R_{\alpha\beta\gamma\delta}R_{\epsilon\zeta\eta\theta},$
$\displaystyle\quad\nabla_{\alpha}\nabla_{\beta}\nabla_{\gamma}T_{\delta\epsilon\zeta},$
$\displaystyle\quad\nabla_{\alpha}T_{\beta\gamma\delta}\nabla_{\epsilon}T_{\zeta\eta\theta},$
(15)
$\displaystyle\nabla_{\alpha}\nabla_{\beta}T_{\gamma\delta\epsilon}T_{\zeta\eta\theta},$
$\displaystyle\quad\nabla_{\alpha}T_{\beta\gamma\delta}T_{\epsilon\zeta\eta}T_{\theta\iota\kappa},$
$\displaystyle\quad
T_{\alpha\beta\gamma}T_{\delta\epsilon\zeta}T_{\eta\theta\iota}T_{\kappa\lambda\mu},$
$\displaystyle\quad\nabla_{\alpha}T_{\beta\gamma\delta}R_{\epsilon\zeta\eta\theta},$
$\displaystyle
T_{\alpha\beta\gamma}\nabla_{\delta}R_{\epsilon\zeta\eta\theta},$
$\displaystyle\quad
T_{\alpha\beta\gamma}T_{\delta\epsilon\zeta}R_{\eta\theta\iota\kappa},$

There are 151 invariants listed in [10] which are too many for analytical
analysis. Nevertheless, note the following. The velocities squared terms
$\dot{T}^{2}$, where the dot denote time derivatives $\partial/\partial x^{0}$
and tensor indices are dropped, are contained in invariants $R^{2}$, $\nabla
T\nabla T$, and $\nabla TR$, and also in invariants $\nabla\nabla T\,T$ and
$\nabla TR$ after integration by parts. In the linear approximation, all
invariants quadratic in velocities have the same form as invariants $\nabla
T\nabla T$. Therefore, it is sufficient to consider only invariants
constructed from the tensor
$\nabla_{\alpha}T_{\beta\gamma\delta}\nabla_{\epsilon}T_{\zeta\eta\theta}$ by
different contractions with the metric.

We start with the component $S_{abc}$. Because this component satisfies Eqs.
(8), there are only three independent invariants with respect to the full
Lorentz group including reflections which are quadratic in first derivatives.
We choose them in the form

$\nabla_{a}S_{bcd}\nabla^{a}S^{bcd},\qquad\nabla_{a}S_{b}{}^{a}{}_{c}\nabla_{d}S^{bdc}\qquad\nabla_{a}S_{b}{}^{a}{}_{c}\nabla_{d}S^{cdb}.$
(16)

The remaining invariants can be expressed through them. The proof of this
statement is by item-by-item examination using the following properties. Two
integration by parts imply relations:

$\begin{split}\nabla_{a}S_{b}{}^{a}{}_{c}\nabla_{d}S^{bdc}\overset{{\sf\,div\,}}{=}&\nabla_{a}S_{bcd}\nabla^{c}S^{bad},\\\
\nabla_{a}S_{b}{}^{a}{}_{c}\nabla_{d}S^{cdb}\overset{{\sf\,div\,}}{=}&\nabla_{a}S_{bcd}\nabla^{c}S^{dab},\\\
\nabla_{a}S_{bc}{}^{a}\nabla_{d}S^{bcd}\overset{{\sf\,div\,}}{=}&\nabla_{a}S_{bcd}\nabla^{d}S^{bca}.\end{split}$

Moreover, Eq. (8) implies identities:

$\begin{split}\nabla_{a}S_{bcd}\nabla^{a}S^{cdb}=&-\frac{1}{2}\nabla_{a}S_{bcd}\nabla^{a}S^{bcd},\\\
\nabla_{[a}S_{bc]d}\nabla^{[a}S^{bc]d}=&\frac{1}{3}\nabla_{a}S_{bcd}\big{(}\nabla^{a}S^{bcd}+2\nabla^{c}S^{abd}\big{)}.\end{split}$

Finally, there is the relation

$0=\nabla_{a}S_{bcd}\big{(}\nabla^{b}S^{cda}+\nabla^{c}S^{dba}+\nabla^{d}S^{bca}\big{)}=\nabla_{a}S_{bcd}\big{(}\nabla^{d}S^{bca}+2\nabla^{c}S^{dba}\big{)}=\\\
=\nabla_{a}S_{bcd}\big{(}\nabla^{d}S^{bca}-2\nabla^{c}S^{bad}+2\nabla^{c}S^{dab}\big{)},$

where Eqs. (8) are used. The last equality after integration by parts is used
to express the invariant $\nabla_{a}S_{bc}{}^{a}\nabla_{d}S^{bcd}$ through
invariants (16).

The absence of ghosts is defined by the velocity quadratic form. Therefore, we
express invariants (16) through $(1+3)$-components and, for simplicity, keep
only quadratic terms:

$\begin{split}\nabla_{a}S_{bcd}\nabla^{a}S^{bcd}\quad\ni\quad&2\dot{S}^{0\\{ij\\}}\dot{S}_{0\\{ij\\}}-2\dot{H}^{\\{ij\\}}\dot{H}_{\\{ij\\}}-3\dot{S}^{0i0}\dot{S}_{0i0}+3\dot{H}^{i}\dot{H}_{i},\\\
\nabla_{a}S_{b}{}^{a}{}_{c}\nabla_{d}S^{bdc}\quad\ni\quad&\dot{S}^{0\\{ij\\}}\dot{S}_{0\\{ij\\}}-\dot{S}^{0i0}\dot{S}_{0i0}+\frac{1}{2}\dot{H}^{i}\dot{H}_{i},\\\
\nabla_{a}S_{b}{}^{a}{}_{c}\nabla_{d}S^{cdb}\quad\ni\quad&\dot{S}^{0\\{ij\\}}\dot{S}_{0\\{ij\\}}-\frac{1}{2}\dot{H}^{i}\dot{H}_{i}.\end{split}$
(17)

The kinetic part of the Lagrangian for torsion component $S_{abc}$ has the
form

$L_{\textsc{k}}:=c_{1}\nabla_{a}S_{bcd}\nabla^{a}S^{bcd}+c_{2}\nabla_{a}S_{b}{}^{a}{}_{c}\nabla_{d}S^{bdc}+c_{3}\nabla_{a}S_{b}{}^{a}{}_{c}\nabla_{d}S^{cdb},$
(18)

where $c_{1,2,3}$ are coupling constants. In terms of $(1+3)$-decomposition
this Lagrangian yields the quadratic form

$L^{(2)}_{\textsc{k}}=\dot{S}^{0\\{ij\\}}\dot{S}_{0\\{ij\\}}(2c_{1}+c_{2}+c_{3})-2\dot{H}^{\\{ij\\}}\dot{H}_{\\{ij\\}}c_{1}-\\\
-\dot{S}^{0i0}\dot{S}_{0i0}(3c_{1}+c_{2})+\dot{H}^{i}\dot{H}_{i}\left(3c_{1}+\frac{1}{2}c_{2}-\frac{1}{2}c_{3}\right).$
(19)

The analysis of the mass term for $S_{abc}$ shows that propagating can be
either separate torsion ${\mathbb{S}}{\mathbb{O}}(3)$-components or pairs
($S_{0\\{ij\\}}$, $H_{i}$) or ($H_{\\{ij\\}}$, $S_{0i0}$).

Since we have only three invariants (16), and there are four independent
${\mathbb{S}}{\mathbb{O}}(3)$-components of torsion $S_{abc}$, then the
absence of three kinetic terms implies the absence of the fourth. Therefore
Lagrangian (18) cannot describe propagation of only one
${\mathbb{S}}{\mathbb{O}}(3)$-irreducible component of $S_{abc}$.

Now we consider the possibility of simultaneous propagation of two
${\mathbb{S}}{\mathbb{O}}(3)$-irreducible components of $S_{abc}$. If the only
propagating components are $S_{0\\{ij\\}}$ and $H_{i}$, then the following
equalities must hold

$c_{1}=0,\qquad c_{2}=0,$

the constant $c_{3}$ being arbitrary. Then the kinetic part of the Lagrangian
is

$L^{(2)}_{\textsc{k}}=c_{3}\dot{S}^{0\\{ij\\}}\dot{S}_{0\\{ij\\}}-\frac{1}{2}c_{3}\dot{H}^{i}\dot{H}_{i}.$

It is clear that for any $c_{3}\neq 0$ the model contains ghost. Consequently,
propagation of only $S_{0\\{ij\\}}$ and $H_{i}$ is impossible.

Similar situation happens with components $H_{\\{ij\\}}$, $S_{0i0}$. If they
are the only propagating ones, then the equalities

$c_{3}=2c_{1},\qquad c_{2}=-4c_{1}$

must hold, where $c_{1}$ is arbitrary. The corresponding kinetic part of the
Lagrangian is

$L^{(2)}_{\textsc{k}}=-2c_{1}\dot{H}^{\\{ij\\}}\dot{H}_{\\{ij\\}}+c_{1}\dot{S}^{0i0}\dot{S}_{0i0}.$

We see that the model contains ghost for any $c_{1}\neq 0$. Therefore
Lagrangian $L^{(2)}_{\textsc{k}}$ cannot describe propagation of only one pair
of components ($H_{\\{ij\\}}$, $S_{0i0}$).

Thus, there is no Lagrangian in the considered class which propagates the
component $S_{abc}$ without ghosts and tachyons.

Now we discuss the irreducible component $T_{a}$. It is the usual vector
field. There are two possibilities of models without ghosts and tachyons.
Lagrangians

$\displaystyle L_{1}:=$
$\displaystyle-\frac{1}{4}\gamma_{1}F^{ab}F_{ab}+\frac{1}{2}m^{2}_{1}T^{a}T_{a},$
(20) $\displaystyle L_{2}:=$
$\displaystyle~{}~{}\frac{1}{2}\gamma_{2}(\nabla_{a}T^{a})^{2}-\frac{1}{2}m^{2}_{2}T^{a}T_{a},$
(21)

where $F_{ab}:=\partial_{a}T_{b}-\partial_{b}T_{a}$, $\gamma_{1,2}>0$ and
$m_{1,2}>0$ are the only ones without ghosts and tachyons. Lagrangian (20)
describes propagation of vector particles $1^{-}$ of mass $m_{1}$. It is the
Proca field. The second Lagrangian (21) propagates scalar particles $0^{+}$ of
mass $m_{2}$.

The same happens for pseudovector torsion component $T^{*}_{a}$. There are
only two ghost and tachyon free Lagrangians:

$\displaystyle L_{3}:=$
$\displaystyle-\frac{1}{4}\gamma_{3}G^{ab}G_{ab}+\frac{1}{2}m^{2}_{3}T^{*a}T^{*}_{a},$
(22) $\displaystyle L_{4}:=$
$\displaystyle~{}~{}\frac{1}{2}\gamma_{4}(\nabla_{a}T^{*a})^{2}-\frac{1}{2}m^{2}_{4}T^{*a}T^{*}_{a},$
(23)

where $G_{ab}:=\partial_{a}T^{*}_{b}-\partial_{b}T^{*}_{a}$, $\gamma_{3,4}>0$
and $m_{3,4}>0$. Lagrangians (22) and (23) do not contain ghosts and tachyons
and propagate particles $1^{+}$ and $0^{-}$ of masses $m_{3}$ and $m_{4}$,
respectively.

Returning to the gravity models with dynamical models, we summarize. We have
found four Lagrangians

$\displaystyle\frac{1}{\sqrt{|g|}}L_{13}:=$
$\displaystyle\kappa\widetilde{R}(g)-\frac{1}{4}\gamma_{1}F^{\alpha\beta}F_{\alpha\beta}+\frac{1}{2}m^{2}_{1}T^{\alpha}T_{\alpha}-\frac{1}{4}\gamma_{3}G^{\alpha\beta}G_{\alpha\beta}+\frac{1}{2}m^{2}_{3}T^{*\alpha}T^{*}_{\alpha},$
(24) $\displaystyle\frac{1}{\sqrt{|g|}}L_{23}:=$
$\displaystyle\kappa\widetilde{R}(g)+\frac{1}{2}\gamma_{2}(\nabla_{\alpha}T^{\alpha})^{2}-\frac{1}{2}m^{2}_{2}T^{\alpha}T_{\alpha}-\frac{1}{4}\gamma_{3}G^{\alpha\beta}G_{\alpha\beta}+\frac{1}{2}m^{2}_{3}T^{*\alpha}T^{*}_{\alpha},$
(25) $\displaystyle\frac{1}{\sqrt{|g|}}L_{14}:=$
$\displaystyle\kappa\widetilde{R}(g)-\frac{1}{4}\gamma_{1}F^{\alpha\beta}F_{\alpha\beta}+\frac{1}{2}m^{2}_{1}T^{\alpha}T_{\alpha}+\frac{1}{2}\gamma_{4}(\nabla_{\alpha}T^{*\alpha})^{2}+\frac{1}{2}m^{2}_{4}T^{*\alpha}T^{*}_{\alpha},$
(26) $\displaystyle\frac{1}{\sqrt{|g|}}L_{24}:=$
$\displaystyle\kappa\widetilde{R}(g)+\frac{1}{2}\gamma_{2}(\nabla_{\alpha}T^{\alpha})^{2}-\frac{1}{2}m^{2}_{2}T^{\alpha}T_{\alpha}+\frac{1}{2}\gamma_{4}(\nabla_{\alpha}T^{*\alpha})^{2}+\frac{1}{2}m^{2}_{4}T^{*\alpha}T^{*}_{\alpha}.$
(27)

For $\kappa>0$, $\gamma_{1,2,3,4}>0$, $m_{1,2,3,4}>0$ these Lagrangians do not
contain ghosts and tachyons. In all four cases, the first term is the usual
Hilbert–Einstein Lagrangian describing the massless graviton $2^{+}$.
Lagrangians $L_{13}$, $L_{23}$, $L_{14}$, and $L_{24}$ propagate in addition
massive degrees of freedom: $(1^{-},1^{+})$, $(0^{+},1^{+})$, $(1^{-},0^{-})$,
and $(0^{+},0^{-})$, respectively.

Mass term in the obtained Lagrangians can be set to zero. Then torsion
describes massless particles of the same spins and parities. The extra gauge
symmetries appear in these cases, and the number of propagating degrees of
freedom reduces. For example, let us put $m_{1}=0$ in $L_{13}$ and $L_{14}$.
Then these Lagrangians are invariant with respect to gauge transformations

$T_{\alpha}\mapsto T_{\alpha}+\partial_{\alpha}\varphi,$

where $\varphi(x)$ is an arbitrary local transformation parameter. These
transformations coincide with the gauge transformations in electrodynamics.
Therefore, if we switch on the interaction with matter fields in the same way,
then the torsion trace $T_{a}$ can be identified with the electromagnetic
potential. Unfortunately, the geometric meaning of such identification is not
clear and requires further analysis.

If $m_{2}=0$ in Lagrangians $L_{23}$ and $L_{24}$, then they are invariant
with respect to gauge transformations

$T^{\alpha}\mapsto T^{\alpha}+\widetilde{\nabla}_{\beta}\omega^{\beta\alpha},$

where $\omega^{\beta\alpha}(x)=-\omega^{\alpha\beta}(x)$ are local parameters.
Note that the number of independent parameters is

$C^{2}_{4}-C^{3}_{4}+C^{4}_{4}=6-4+1=3,$

which is necessary for elimination of three unphysical degrees of freedom of
vector field $T^{\alpha}$.

The same happens with the pseudovector field $T^{*}_{a}$.

##  6 Conclusion

Thus, supposing that kinetic terms of vierbein $e_{\alpha}{}^{a}$ and torsion
$T_{\alpha\beta\gamma}$ do not mix, we proved that there are only four
Lagrangians $L_{13}$, $L_{23}$, $L_{14}$, and $L_{24}$ without ghosts and
tachyons. In all cases, the vierbein describes one massless graviton $2^{+}$.
In addition, torsion describes massive or massless particles $1^{-}$, $1^{+}$
($L_{13}$); $0^{+}$, $1^{+}$ ($L_{23}$); $1^{-}$, $0^{-}$ ($L_{14}$); or
$0^{+}$, $0^{-}$ ($L_{24}$). Moreover, we proved that any other Lagrangian for
torsion free of ghosts and tachyons must have the same linear approximation as
one of the listed.

We proved also that the irreducible component $S_{abc}$ of torsion cannot have
a Lagrangian without ghosts and tachyons. The same result was obtained in [11]
using different causal analysis. In addition, the causal analysis yields the
argument in favor of Lagrangian (24) of the present paper [12].

This work was performed at the Steklov International Mathematical Center and
supported by the Ministry of Science and Higher Education of the Russian
Federation (agreement no. 075-15-2019-1614).

## References

  * [1] F. W. Hehl, P. von der Heyde, G. D. Kerlick, and J. M. Nester.  General relativity with spin and torsion: foundations and prospects.  Rev. Mod. Phys., 48(3):393–416, 1976.
  * [2] F. W. Hehl, J. D. McCrea, E. W. Mielke, and Y. Ne’eman.  Metric-affine gauge theory of gravity: Field equations, Noether identities, world spinors, and breaking of dilaton invariance.  Phys. Rep., 258(1&2):1–171, 1995.
  * [3] I. L. Shapiro.  Physical aspects of the space-time torsion.  Phys. Rep., 357:113–213, 2002.  arXiv:1204.5505 [astro-ph.CO].
  * [4] M. Blagojević.  Gravitation and Gauge Symmetries.  IOP, Bristol and Philadelphia, 2002.
  * [5] S. Capozziello and M. De Laurentis.  Extended theories of gravity.  Phys. Rep., 509:167–321, 2011.  arXiv:1204.5505 [astro-ph.CO].
  * [6] D. E. Neville.  Gravity lagrangian with ghost-free curvature-squared terms.  Phys. Rev., D18(10):3535–3543, 1978.
  * [7] E. Sezgin and P. van Nieuwenhuizen.  New ghost-free gravity lagrangians with propagating torsion.  Phys. Rev., D21(12):3269–3280, 1980.
  * [8] E. Sezgin.  Class of ghost-free gravity lagrangians with propagating torsion.  Phys. Rev., D24(6):1677–1680, 1981.
  * [9] K. Hayashi and T. Shirafuji.  Gravity from Poincaré gauge theory of the fundamental particles. IV. Mass and energy of particle spectrum.  Prog. Theor. Phys., 64(6):2222–2241, 1980.
  * [10] S. M. Christensen.  Second- and fourth-order invariants on curved manifolds with torsion.  J. Phys., A13(9):3001–3009, 1980.
  * [11] L. Fabbri and M. Tecchiolli.  Restrictions on torsion-spinor field theory.  Mod. Phys. Lett. A, 34:1950311, 2019.  arXiv:1811.12060 [physics.gen-ph].
  * [12] L. Fabbri.  A discussion on the most general torsion-gravity with electrodynamics for Dirac spinor matter fields.  Int. J. Geom. Meth. Mod. Phys., 12:1550099, 2015.  arXiv:1409.2007 [gr-qc].