# The Bad Locus in the moduli of
Super Riemann Surfaces with Ramond punctures

Ron Donagi  Department of Mathematics, University of Pennsylvania,
Philadelphia, PA 19104<EMAIL_ADDRESS>and  Nadia Ott  Department of
Mathematics, University of Pennsylvania, Philadelphia, PA 19104
<EMAIL_ADDRESS>
###### Abstract.

The bad locus in the moduli of super Riemann surfaces with Ramond punctures
parametrizes those super Riemann surfaces that have more than the expected
number of independent closed holomorphic 1-forms. There is a super period map
that depends on certain discrete choices. For each such choice, the period map
blows up along a divisor that contains the bad locus. Our main result is that
away from the bad locus, at least one of these period maps remains finite. In
other words, we identify the bad locus as the intersection of the blowup
divisors. The proof abstracts the situation into a question in linear algebra,
which we then solve. We also give some bounds on the dimension of the bad
locus.

###### Contents

  1. 1 Introduction
  2. 2 Periods of super Riemann surfaces with Ramond punctures
  3. 3 Moduli spaces and the main theorem
  4. 4 Linear algebra
  5. 5 From Ramond punctures to linear algebra
  6. 6 Dimension counts
  7. A Proof of the Linear Algebra Result
  8. B Odd periods
    1. B.1  $\Omega_{1}=H^{0}(C,L)$
    2. B.2 Odd periods

##  1\. Introduction

The super vector space $\Omega$ of closed holomorphic differentials on a
general super Riemann surface of genus $g$ with $2r$ Ramond punctures has
dimension $g|r$. The super period map as defined by Witten in [W15] sends this
space to a $g|r$ dimensional maximal isotropic subspace $P(\Omega)$ of period
space, which is a $2g|2r$-dimensional super vector space with a
supersymplectic pairing. Unfortunately there is a bad locus $B$ in the moduli
space parametrizing super Riemann surfaces on which the odd dimension of
$\Omega$ is strictly greater than $r$, so the period map fails to be injective
there. This bad locus $B$ was studied by Witten in [W15], mostly in the case
$r=1$. The goal of this note is to study the bad locus for higher $r$.

When $r=1$, the bad locus has codimension 1 and it is the locus where Witten’s
super period map blows up. We will see in Theorem 8 that in general the bad
locus (or at least one of its components) has codimension $r$. The period map
depends on some discrete choices: a symplectic basis of the integral homology,
as in the bosonic case, plus an orientation $\psi$ of each of the Ramond
punctures. Each choice of $\psi$ determines a period matrix that blows up
along a certain divisor $Y^{\prime\prime}_{\sigma}$. Here
$\sigma=\sigma(\psi)$ is a sign choice, determined by $\psi$ but retaining
less information than $\psi$: there are $2^{2r}$ choices of $\psi$, but only
$2^{r}$ choices for $\sigma$. The blowup divisor $Y^{\prime\prime}_{\sigma}$
depends only on ${\sigma}$. Our main result, Theorem 2, is that the
intersection $Y$ of these divisors equals the bad locus $B$. The
$Y^{\prime\prime}_{\sigma}$ therefore give a more-or-less explicit set of
equations for the bad locus. (This story has a variant that depends
additionally on the choice of a pairing $\nu$ of the $2r$ punctures. It
follows from our result that this intersection of the
$Y^{\prime\prime}_{\sigma}$ is independent of the choice of $\nu$.)

We review the setup and state the main result in section 2. Our strategy for
the proof is to abstract the situation and translate it to a question in
linear algebra, Conjecture 3. A special case of this conjecture, sufficient
for our needs here, is proved in Appendix A. The full conjecture was proved in
[D22]. The reduction of the main theorem to this linear algebra result is
explained in section 5.

During the preparation of this work, Ron Donagi was supported in part by NSF
grant DMS 2001673 and by Simons HMS Collaboration grant #390287. Nadia Ott was
supported in part by Simons HMS Collaboration grant #390287.

##  2\. Periods of super Riemann surfaces with Ramond punctures

A genus g super Riemann surface $X$ with $2r$ Ramond punctures determines a
curve $C$ of the same genus $g$ with a divisor $D=\Sigma_{i=1}^{2r}x_{i}$
consisting of $2r$ distinct labeled points and a twisted spin structure, i.e.
a line bundle $L$ satisfying

(1) $L^{\otimes 2}=K_{C}(D).$

Riemann-Roch gives

(2) $\chi(L)=h^{0}(L)-h^{0}(L(-D))=r.$

For $r>0$ an easy general position argument shows that generically
$h^{0}(L)=r$. This means that the set of super Riemann surfaces with
$h^{0}(L)=r$ can be identified with an open subset of the moduli space

(3) $M:=\\{(C,D,L)\ |\ C\ \text{a curve of genus}\ g,D\in C^{2r},\
L\in{\text{Pic}}^{g+r-1}(C)\ {\text{satisfies \eqref{spin}}}\\}.$

Following Witten [W15], we refer to the complement

$B:=\\{(C,D,L)|h^{0}(L)\geq r+1\\}\subset M.$

of this open subset as the “bad locus”.

The super vector space $\Omega$ of closed 1-forms on $X$ has dimension
$g|h^{0}(L)$. (Our sign convention is that when $X$ is split, the forms pulled
back from the reduced curve $C$ are considered even, e.g, for local
coordinates $(z,\theta)$ on $X$, $dz$ is even and $d\theta$ is odd. Both this
and the opposite convention occur in the literature.) In the split case, the
even closed forms are pulled back from $H^{0}(C,K_{C})$ while the odd closed
forms are exact, of the form $da(z)\theta$ where $a=a(z)$ is a section of the
twisted spin bundle $L$, cf. [W15], section 5.1, or our appendix B.2. So
generically $\Omega$ is of dimension $g|r$.

We recall Witten’s superperiod map

$P:\Omega\to{\mathbb{C}}^{2g|2r}.$

This depends on some discrete choices. The $2g$ even periods, just as in the
classical theory, depend on the choice of a symplectic basis
$(a_{1},\dots,a_{g};b_{1},\dots,b_{g})$ of $H_{1}(X,\mathbb{Z})$. Since the
residue gives a natural trivialization of the fiber of $K_{C}(D)$ at each
puncture $x_{i}\in D$, equation (1) specifies a trivialization of the square
$(L|{x_{i}})^{2}$ of the fiber of $L$ at $x_{i}$. By an orientation $\psi_{i}$
of the puncture we mean a square root of this, i.e. a trivialization of
$L|{x_{i}}$ that squares to the above trivialization of $K_{C}(D)$. The $2r$
odd periods depend on a choice of orientation at each of the $2r$ punctures.

These discrete choices determine identifications

$\Lambda_{0}=H^{1}(C,{\mathbb{C}})$

and

$\Pi\Lambda_{1}=H^{0}(L/L(-D))=\bigoplus_{i=1}^{2r}L_{p_{i}}.$

The even part of the superperiod map is then restriction from $X$ to $C$
followed by integration on the symplectic basis
$(a_{1},\dots,a_{g};b_{1},\dots,b_{g})$. The odd part is just the restriction
from $C$ to $D$:

$H^{0}(C,L)\to H^{0}(D,L)=H^{0}(L/L(-D)),$

cf. Appendix B.2.

One bad thing about the “bad” locus $B$ is that this superperiod map $P$ is
not injective there: On $X\in B$ there are (odd) closed $1$-forms with all
periods equal to zero, as can be seen from the long exact sequence

(4) $0\to H^{0}(L(-D))\to H^{0}(L)\overset{P_{1}}{\to}H^{0}(L/L(-D))\to
H^{1}(L(-D))\to\dots$

where $h^{0}(C,L(-D))\geq 1$.

The target of the superperiod map is period superspace
$\Lambda\cong{\mathbb{C}}^{2g|2r}$. This comes equipped with a super
symplectic pairing $Q$ (symplectic on ${\mathbb{C}}^{2g}$, orthogonal on
${\mathbb{C}}^{0|2r}$) which we may think of as a super version of the usual
intersection pairing on $H^{1}(X,\mathbb{R})$. In the obvious notation,

$Q=\Sigma_{k=1}^{g}a_{k}\wedge b_{k}-\Sigma_{i=1}^{2r}(\psi_{i})^{2},$

where we recall that the orientations $\psi_{i}$ form a basis for
$\Lambda_{1}$. The period map can be described by a $(g|r)\times(2g|2r)$
matrix $\mathcal{P}$. In order to obtain a $(g|r)\times(g|r)$ period matrix,
we need to split period space as a sum of two complementary Q-Lagrangian
subspaces, $\Lambda=\Lambda_{A}\oplus\Lambda_{B}$. For the even part we use
the span of $(a_{1},\dots,a_{g})$ for $\Lambda_{A,0}$ and of
$(b_{1},\dots,b_{g})$ for $\Lambda_{B,0}$. For the odd part we fix a pairing
$\nu$ of the $2r$ punctures, e.g. we could pair $x_{2i-1}$ with $x_{2i}$. We
then take $\Lambda_{A,1}$ as the span of $x_{2i-1}+\sqrt{-1}x_{2i}$ and
$\Lambda_{B,1}$ as the span of $x_{2i-1}-\sqrt{-1}x_{2i}$. The period map then
decomposes: $P=(P_{A},P_{B})$, with corresponding decomposition of the matrix,
$\mathcal{P}=(\mathcal{P}_{A},\mathcal{P}_{B})$. The $(g|r)\times(g|r)$ period
matrix is then $\mathcal{P}_{A}^{-1}\mathcal{P}_{B}$ . This is defined only
where $\mathcal{P}_{A}$ is invertible. It depends on:

  * •

the data $X=(C,D,L)$ underlying the super Riemann surface with labeled Ramond
punctures,

  * •

the symplectic basis $(a_{1},\dots,a_{g};b_{1},\dots,b_{g})$ of
$H_{1}(X,\mathbb{Z})$, and

  * •

the orientations $\psi_{i}$.

Sometimes we may want to replace the specific pairing ($x_{2i-1}$ with
$x_{2i}$) by an arbitrary pairing $\nu$ and to emphasize the dependence on
this choice. A pairing $\nu$ is specified by relabeling the $2r$ points as:

$x_{ij},\ \ \ i=1,\dots,r,\ \ \ j\in\pm:=\\{+,-\\}.$

By gluing the points $(x_{i+},x_{i-})$ in each pair we obtain a nodal curve
$\overline{C}$ (of arithmetic genus $g+r$, with $r$ nodes). Specifying the
pairing of the points is equivalent to specifying the normalization map
$\nu:C\to\overline{C}$.

Having specified the pairing $\nu$, a sign choice $\sigma$ is a specification
of a spin bundle, aka theta characteristic, on the nodal curve $\overline{C}$,
_i.e_ , a line bundle $L_{\sigma}$ on the nodal curve $\overline{C}$
satisfying

(5) $L_{\sigma}^{\otimes 2}=K_{\overline{C}}.$

In particular, $L:=\nu^{*}L_{\sigma}$ is a twisted spin structure on $C$
satisfying (1). Equivalently, given $L$, $\sigma$ specifies for each
$i=1,\dots,r$ an isomorphism

(6) $\sigma_{i}:L_{x_{i+}}\to L_{x_{i-}}$

whose square:

$\sigma_{i}^{\otimes 2}:K_{C}(D)_{|{x_{i+}}}=L_{{x_{i+}}}^{\otimes 2}\to
L_{{x_{i-}}}^{\otimes 2}=K_{C}(D)_{|{x_{i-}}}$

is identified via residues with multiplication by $-1$,

$\sigma_{i}^{\otimes 2}=-1:{\mathbb{C}}\to{\mathbb{C}}.$

($L_{\sigma}$ clearly determines a $\sigma$ as in (6). Conversely, given
$\sigma$ as in (6), we recover $L_{\sigma}$ as the subsheaf of $L$ consisting
of sections $s$ of $L$ satisfying

(7) $s(x_{i-})=-\sqrt{-1}\sigma_{i}s(x_{i+})$

for all $1\leq i\leq r$. In particular, $h^{0}(\overline{C},L_{\sigma})>0$ if
$L$ has a _global_ section satisfying the gluing equations (7) for all $1\leq
i\leq r$. )

Recall that an orientation $\psi$ is a trivialization of $L|_{D}$ compatible
via residues with (1). In our current notation this amounts to trivializations
$\psi_{ij}$ of each $L_{x_{ij}}$. Such a $\psi$ combined with formula (6)
therefore determines a $\sigma=\sigma(\psi)$. Changing $\psi$ at both $x_{i+}$
and $x_{i-}$ leaves $\sigma$ unchanged. On a given $X$ there are $2^{2r}$
$\psi$’s but only $2^{r}$ $\sigma$’s.

Recall that the odd parts

$P_{1}:\Omega_{1}\to\Lambda_{1},\ P_{A,1}:\Omega_{1}\to\Lambda_{A,1}$

of the superperiod map depend on the orientation $\psi$. The kernels, however,
depend only on $\sigma=\sigma(\psi)$. Indeed, if we identify $\Omega_{1}$ with
$H^{0}(C,L)$, we find the identifications :

(8) $Ker(P_{1})=H^{0}(L(-D))\subset
Ker(P_{A,1})=H^{0}(\overline{C},L_{\sigma})\subset H^{0}(C,L).$

We will write $\mathcal{P}^{\nu,\sigma}$ when we want to emphasize that the
matrix $\mathcal{P}$ representing the period map depends on a choice of
pairing $\nu$ and sign $\sigma$.

##  3\. Moduli spaces and the main theorem

We are interested in describing the locus in moduli space where the period
matrix is not defined, i.e., the locus where $\mathcal{P}_{A}$ is singular for
all possible signs $\sigma$ and for some or all pairings $\nu$. To do this, we
need to consider some covers of the moduli space $M$ of (3). Let

$M^{\prime}:=\\{\nu:C\to\overline{C},D,L\\}$

be the (bosonic) moduli space parametrizing the data of a curve $C$ with
labeled divisor $D$, a pairing $\nu$ on the points of $D$, and a twisted spin
bundle $L$ satisfying (1). There is a forgetful map

$\pi:M^{\prime}\to M.$

There is a further $2^{r}$-sheeted cover

$\pi^{\prime}:M^{\prime\prime}\to M^{\prime}$

where

$M^{\prime\prime}:=\\{\nu:C\to\overline{C},D,L_{\sigma}\\}$

parametrizes curves $C$ with labeled divisor $D$, a pairing of the points of
$D$, and a choice of a spin bundle $L_{\sigma}$ satisfying (5).

Let $\mathcal{P}$ denote (the component of degree $g+r-1$ of) the relative
Picard of the universal curve over the Deligne-Mumford compactification of the
moduli space $M_{g+r}$ of curves of genus $g+r$, so

$\mathcal{P}:=\\{C^{\prime},L^{\prime}\ |\ C^{\prime}\ {\text{a stable curve
of genus}}\ g,L^{\prime}\in{\text{Pic}}^{g+r-1}(C^{\prime})\\}.$

In $\mathcal{P}$ there is a natural divisor $\mathcal{P}_{e}$ parametrizing
pairs $\\{C^{\prime},L^{\prime}\\}$ where $L^{\prime}$ is an effective line
bundle: $h^{0}(L^{\prime})>0$. There is a natural map of $M^{\prime\prime}$ to
$\mathcal{P}$ sending $\\{\nu:C\to\overline{C},D,L_{\sigma}\\}$ to
$\\{\overline{C},L_{\sigma}\\}.$ We let $Y^{\prime\prime}$ denote the inverse
image of the divisor $\mathcal{P}_{e}$:

$Y^{\prime\prime}=Y^{\prime\prime}_{\sigma}:=\\{(\nu:C\to\overline{C},D,L_{\sigma})\
|\ L_{\sigma}^{\otimes 2}=K_{\overline{C}},\ h^{0}(L_{\sigma})\geq 1\\}\subset
M^{\prime\prime}.$

This is a divisor in $M^{\prime\prime}$. If $(C,D,L,\nu,\sigma)\in
Y^{\prime\prime}$ then there exists a global section of $L$ satisfying (7) for
all $1\leq i\leq r$, i.e., there exists a global section $s$ whose odd
A-periods all vanish. So $P_{A,1}^{\nu,\sigma}$ has a non-trivial kernel and
the matrix $\mathcal{P}_{A,1}^{\nu,\sigma}$ is singular.

As shown in Appendix A of [W15], the period matrix (viewed as a map on
$M^{\prime\prime}$) blows up (has a pole) along this divisor. This divisor
$Y^{\prime\prime}_{\sigma}$ does not depend on a choice of sign $\sigma$; the
subscript $\sigma$ simply reminds us that it is a divisor in
$M^{\prime\prime}$, whose fiber over a point of $M^{\prime}$ is indexed by the
$\sigma$’s.) The reason that $\sigma$ occurs here rather than the full
orientation $\psi$ was explained at the end of the previous section:
$Y^{\prime\prime}$ parametrizes points where the (odd part of the) superperiod
map fails to be injective. According to (8), this depends only on $\sigma$
rather than $\psi$.

Descending to $M^{\prime}$, we have the closed locus

$Y^{\prime}:=\\{(\nu:C\to\overline{C},D,L)\ |\ L^{\otimes 2}=K_{C}(D),\
h^{0}(L_{\sigma})\geq 1\ \forall\sigma\\}=\\{y\in M^{\prime}\ |\
\pi^{\prime-1}(y)\subset Y^{\prime\prime}\\}\subset M^{\prime}$

parametrizing super Riemann surfaces with a given pairing of their $2r$ Ramond
punctures $D$ such that the period matrix blows up for every sign choice
$\sigma$, i.e., $(C,D,L,\nu)\in Y^{\prime}$ if $\mathcal{P}_{A}^{\nu,\sigma}$
is singular for all $\sigma$. So the complement of $Y^{\prime}$ is the open
subset where the period matrix, for at least one sign choice $\sigma$, makes
sense. Somewhat informally we write
$Y^{\prime}=\cap_{\sigma}Y^{\prime\prime}_{\sigma}$. (Informal because
$Y^{\prime\prime}\subset M^{\prime\prime}$ while $Y^{\prime}\subset
M^{\prime}$. )

Descending further to $M$, we have the closed locus

$Y:=\\{(C,D,L)|(\nu:C\to\overline{C},D,L)\in Y^{\prime}\ \forall\nu\\}=\\{y\in
M\ |\ \pi^{-1}(y)\subset Y^{\prime}\\}\subset M,$

where $(C,D,L)\in Y$ if $\mathcal{P}_{A}^{\nu,\sigma}$ is singular for all
pairings $\nu$ and signs $\sigma$, and also the a priori larger locus

$\tilde{Y}:=\\{(C,D,L)|(\nu:C\to\overline{C},D,L)\in Y^{\prime}\ {\text{for
some}}\ \nu\\}=\pi(Y^{\prime})\subset M.$

The distinction is that $(C,D,L)$ is in $\tilde{Y}$ if $(C,D,L,\nu)\in
Y^{\prime}$ for some pairing $\nu$, while in order to be in $Y$ this must hold
for all $\nu$ .

######  Lemma 1.

$B\subseteq Y\subseteq\tilde{Y}$

###### Proof.

If $(C,D,L)\in B$, then $h^{0}(L(-D))\geq 1$ and so there exists a global
section of $L$ satisfying the gluing equation (7) for all $1\leq i\leq r$ and
for all pairings $\nu$ and signs $\sigma$. ∎

Our main result can now be stated as follows:

######  Theorem 2.

$B=Y=\tilde{Y}$.

##  4\. Linear algebra

Consider a set of $2r$ abstract points

$p_{ij},i=1,\dots,r,\ j\in\pm:=\\{+1,-1\\}.$

By a section, or $n$-section, we mean a function

$\sigma:\\{1,\dots,r\\}\to\pm,$

or equivalently the set of $r$ points

$p_{i,\sigma(i)},\ (i=1,\dots,r).$

Our study of Ramond punctures leads to the following:

######  Conjecture 3.

Let $V$ be a vector space of dimension $\geq r$ containing $r$ pairs of points
$p_{ij},i=1,\dots,r,\ j\in\pm$ that span $V$. The following conditions (*),
(**) are equivalent:

(*) For each of the $2^{r}$ sections $\sigma$, the $r$ points
$p_{i,\sigma(i)},\ (i=1,\dots,r)$ are linearly dependent.

(**) There is a subset $I\subset\\{1,\dots,r\\}$ of some cardinality $k<r$
such that the $2k$ points $p_{ij},\ i\in I,\ j\in\pm$ span a subspace of
dimension $<k$.

In section A of the appendix we prove

######  Theorem 4.

The conjecture holds when $V$ is $r$ dimensional.

We will need only this case. The full Conjecture has been proved in [D22].

##  5\. From Ramond punctures to linear algebra

In this section we reduce Theorem 2 to the linear algebra result, Theorem 4.

We begin by proving that condition (*) in Conjecture 3 and Theorem 4
characterizes the locus $\tilde{Y}$. Let $V$ be the vector space

$V:=(H^{0}(L)/H^{0}(L-D))^{*}$

of linear functions on $H^{0}(L)$ that vanish on $H^{0}(L-D)$. (Equivalently,
it is the dual of the image of the odd superperiod map $P_{1}$.) By Riemann-
Roch, its dimension is $r$. We are interested mostly in the case that
$h^{0}(L)=r$, so $H^{0}(L-D)=0$ and $V=H^{0}(L)^{*}$. For each
$i\in\\{1,\dots,r\\}$ and $j\in\pm$, fix an orientation $\psi_{ij}$, i.e. a
trivialization of the fiber $L_{|x_{ij}}$ of $L$ at $x_{ij}$, compatible with
the residue trivialization of $L^{\otimes 2}=K_{C}(D)$. Evaluation at $x_{ij}$
then becomes a vector $q_{ij}\in V$. We let $p_{ij}:=q_{i+}+j\sqrt{-1}q_{i-}$,
or explicitly:

$p_{i+}:=q_{i+}+\sqrt{-1}q_{i-},\ \ \ p_{i-}:=q_{i+}-\sqrt{-1}q_{i-}.$

The collection of $2r$ vectors $q_{ij}$ clearly spans $V$, and therefore so do
the $p_{ij}$. A sign choice $\sigma$ picks, for each $i$, a vector
$p_{i\sigma(i)}$ which is one of $p_{i+},p_{i-}$. The condition for being in
$Y^{\prime\prime}_{\sigma}$ is that the $r$ points $p_{i\sigma(i)},\ \
i=1,\dots,r$ should be contained in a hyperplane in $V$. Condition (*) of
Theorem 4 therefore characterizes $\tilde{Y}$.

Next we describe a certain locus $W$ in moduli space which we will show is
equal to both $B$ (Lemma 7) and $\tilde{Y}$ (Lemma 6), and thereby we prove
Theorem 2. For every subset

$I\subset\\{1,\dots,r\\},$

the divisor

$D_{I}:=\Sigma_{i\in I}(x_{i+}+x_{i-})$

which is of degree $2k$ where $k:=\\#I$, determines the locus:

$W_{I}:=\\{(C,D,L)|h^{0}(L(-D_{I}))\geq r+1-k\\}\subset M.$

For example, the bad locus $B$ equals $W_{\emptyset}$. Now let

$W:=\cup_{I}W_{I}.$

######  Note 5.

These loci obey a duality. Let $k:=\\#I$ and let $I^{\prime}$ denote the
complement of $I$. We have

$h^{0}(L(-D_{I}))=\chi(L(-D_{I}))+h^{1}(L(-D_{I}))=r-2k+h^{0}(K\otimes
L^{-1}(D-D_{I^{\prime}}))=r-2k+h^{0}(L-D_{I^{\prime}}))$

so $W_{I^{\prime}}=W_{I}$.

######  Lemma 6.

Condition (**) of Conjecture 3 and Theorem 4 characterizes $W$. In particular,
$W=\tilde{Y}$.

###### Proof.

Let $k:=\\#I$. The condition that the $2k$ points $p_{ij},\ i\in I,\ j\in\pm$
are contained in a subspace of dimension $k-1$ implies that we are in $W_{I}$.
The converse fails in general, because $h^{0}(L)$ can be $>r$ if we are at a
point of $B$. However, if we are at a point of $\tilde{Y}$ that is not in
$B=W_{\emptyset}$, then $h^{0}(L)=r$, $H^{0}(L-D)=0$ and $V=H^{0}(L)^{*}$, and
then the condition that the $2k$ points $p_{ij},\ i\in I,\ j\in\pm$ are
contained in a subspace of dimension $k-1$ is equivalent to being in $W_{I}$.
So Condition (**) of Theorem 4 characterizes $W$. ∎

The proof of Theorem 2 is now reduced to the following result:

######  Lemma 7.

We have the inclusion and equality:

$B=W_{\emptyset}=\cup_{I}W_{I}=W.$

In particular, $B=Y=\tilde{Y}$.

###### Proof.

A point of $W_{I}$ is given by $(C,D,L)$ such that $h^{0}(L(-D_{I}))\geq
r+1-\\#I$. For elements of $H^{0}(L(-D_{I}))$, the conditions for being in
some $Y^{\prime\prime}_{\sigma}$ are clearly satisfied automatically at points
of $I$ since such an element vanishes along the divisor $D_{I}$, so we only
need to impose $r-\\#I$ conditions at the remaining points in $D-D_{I}$,
ending with $h^{0}(L_{\sigma})\geq 1$. Thus $W=\cup_{I}W_{I}\subset Y$.

On the other hand, each $W_{I}$ is contained in $B$. Indeed, say we have a
point of $W_{I}$ that is not in $B$. Then $h^{0}(L)=r$ and the divisor $D_{I}$
imposes at most $k-1$ conditions on sections of $L$, where we assume $1\leq
k:=\\#I\leq r-1$. But the equality $W_{I}=W_{I^{\prime}}$ implies that the
divisor $D_{I^{\prime}}$ imposes at most $r-k-1$ conditions on sections of
$L$. So $h^{0}(L-D)\geq r-(k-1)-(r-k-1)=2$, contradicting (2).

∎

##  6\. Dimension counts

Among general line bundles $L$ of degree $g-1+r$ on a curve $C$ of genus $g$,
for $0\leq r<g$, the locus where $h^{0}(L)>r$ has codimension $r+1$ by Brill-
Noether theory [BN, ACGH]. One might therefore naively expect that the bad
locus $B$ should have codimension $r+1$ in $M$. The residue theorem implies
that this expectation is wrong: it is off by at at least 1. In this section we
show that there are no further corrections:

######  Theorem 8.

  * •

In the unpunctured case $r=0$, the moduli space $M=M^{+}\cup M^{-}$ is
reducible. The bad locus $B$ consists of all of $M^{-}$ plus the irreducible
“vanishing thetanull” divisor in $M^{+}$.

  * •

For $1\leq r<g$, every component of $B$ has codimension $\leq r$ and there is
a component of codimension exactly $r$.

  * •

For $r\geq g$ the bad locus is empty.

###### Proof.

The case $r=0$ is well known: $M^{-}$ is the moduli of odd spin bundles $L$,
where $h^{0}(L)$ is odd and in particular $h^{0}(L)\geq 1$. $M^{+}$ is the
moduli of even spin bundles $L$, for which $h^{0}(L)=0$ generically and
$h^{0}(L)$ is even and $\geq 2$ along the irreducible divisor of vanishing
thetanulls [A71, M71].

For $r\geq g$, the degree of a twisted spin bundle $L$ is $\geq 2g-1$, so
$h^{1}(L)=0$ and $h^{0}(L)=r$, so the bad locus is empty.

Assume $1\leq r<g$. To see that the codimension is at most $r$ (and thus off
by 1 from the Brill-Noether expectation), it is convenient to switch to the
$2r$-sheeted cover

$\tilde{M}:=\\{(C,D,L,p)\ |\ (C,D,L)\in M\ {\text{and}}\ p\in D\\}$

and to $\tilde{B}$, the inverse image there of $B$. Consider the map

$f:\tilde{M}\to{\text{Pic}}^{g-r}(\mathcal{C}/M),$

sending

$(C,D,L,p)\mapsto N:=L(p-D).$

Here $\mathcal{C}$ is the universal curve over $M$. In the Picard, the
effective locus

$G^{0}_{g-r}:=\\{N\ |\ h^{0}(N)>0\\}$

has codimension $r$. We claim that

$\tilde{B}=f^{-1}(G^{0}_{g-r}),$

which implies that the codimension is at most $r$.

Indeed in one direction,

$L\in\tilde{B}\iff h^{0}(L)>r\iff h^{0}(L(-D))>0\implies h^{0}(N)>0.$

But conversely, a non-zero section $s\in H^{0}(N)$ gives a meromorphic 1-form

$s^{2}\in H^{0}(K(-(D-p))(p))$

with a possible pole only at $p$. By the residue theorem, it actually cannot
have a pole there, so $s$ must be a non-zero section in
$H^{0}(N-p)=H^{0}(L-D)$.

To conclude, we need to exhibit a component of codimension $r$. Let $A$ be a
generic odd theta characteristic:

$A^{\otimes 2}=K,h^{0}(A)=1.$

Let $a\in{\text{Sym}}^{g-1}(C)$ be the unique divisor of a non-zero section of
$A$, and let $d\in{\text{Sym}}^{r}(C)$ be the sum of any $r$ of these points.
(Recall we are assuming $r\leq g-1$.) Let $p\in d$ be one of these points.
Generically such $p$ is not in $a-d$, which we will assume. Take
$D:=2d,L:=A(d)$. We claim that $f(\tilde{M})$ is transversal to $G^{0}_{g-r}$
at $f(C,D,L)$, so the component of $\tilde{B}$ through $(C,D,L)$ has
codimension $r$.

This is seen by comparing the tangent spaces to $f(\tilde{M})$ and
$G^{0}_{g-r}$ at

$f(C,D,L)=N:=L(p-D)=A(d+p-2d)=A(-(d-p)).$

The tangent there to $G^{0}_{g-r}$ is spanned by the $g-r$ tangent lines to
the Abel-Jacobi image of $C$ in its Jacobian at the points of $a-d+p$, or
projectively by the corresponding $g-r$ points of the canonical image of $C$
in canonical space. The tangent to $f(\tilde{M})$ is likewise spanned by the
$2r$ tangent lines to the Abel-Jacobi (or: canonical) image of $C$ at the
points of $D=2d$. (Since our $D$ is non-reduced, this really means the $r$
osculating planes to $C$ at the points of $d$.) This follows from the relation

$N^{\otimes 2}=K(-(D-p)+p),$

which implies that the differential of $f$ at each of these points is
$\pm\frac{1}{2}$ times the corresponding differential of Abel-Jacobi, so the
spans are the same.

The transversality is therefore equivalent to saying that the divisor

$2d+(a-d+p)=a+d+p$

is not contained in any canonical divisor. From our assumption that
$h^{0}(A)=1$ it follows that the unique canonical divisor containing $a$
equals $2a$. This does not contain $a+d+p$ because of our assumption that
$p\notin a-d$.

∎

##  Appendix A Proof of the Linear Algebra Result

We rephrase Theorem 4 as follows:

######  Theorem 9.

Let $\\{P_{1},\dots,P_{r}\\},\ P_{i}\neq 0$ be a set of nonzero subspaces of
the $r$ dimensional vector space $V=(H^{0}(L)/H^{0}(L-D))^{*}$ such that
$\operatorname{dim}P_{i}\leq 2$. For each $P_{i}$, fix two spanning vectors
$\\{p_{i},p_{i}^{\prime}\\}$ (with one necessarily redundant if
$\operatorname{dim}P_{i}=1$) and define

$V_{r}=\\{\\{v_{1},\dots,v_{r}\\}\ |\ v_{i}\in\\{p_{i},p_{i}^{\prime}\\}\\}.$

If no element of $V_{r}$ spans $V$, then there exists a subset
$\\{P_{i_{1}},\dots,P_{i_{k}}\\}$ for $k\leq r$ such that

$\operatorname{dim}\operatorname{span}\\{P_{i_{1}},\dots,P_{i_{k}}\\}<k.$

Let $\rho=\operatorname{max}\\{\ \operatorname{dim}(\operatorname{span}v)\ |\
\forall v\in V_{r}\\}$ and let $V_{r}^{\rho}=\\{v\in V_{r}\
|\operatorname{dim}\operatorname{span}v=\rho\\}$. For each $v\in V_{r}$ let
the set $v^{\prime}\in V_{r}$ be defined by the rule that if $v_{i}\in v$ then
$v_{i}^{\prime}\in v^{\prime}$. Let

$v_{-j}=\\{v_{1},\dots,\hat{v}_{j},\dots,v_{r}\\},\ \
v_{+i^{\prime}}=\\{v_{1},\dots,v_{r},v_{i}^{\prime}\\}.$

For $A\subset(v\cup v^{\prime})\cap W$, the square brackets $[A]$ will denote
the indices of the elements in $A$. For example, if
$A=\\{v_{1},v_{1}^{\prime},v_{2},v_{3}\\}$, then $[A]$ is the multi-set
$\\{1,1,2,3\\}$ where the two copies of $1$ correspond to $v_{1}$ and
$v_{1}^{\prime}$. We say that $i\in[A]$ is a _double_ if $v_{i}\in A$ and
$v_{i}^{\prime}\in A$.

######  Definition 10.

Let $v\in V_{r}^{\rho}$, let $W=\operatorname{span}v$, and let
$W^{\prime}=\operatorname{span}v^{\prime}$. We say that $w\subset(v\cup
v^{\prime})\cap W$ is a _spanning set_ if it satisfies the following
conditions: (1) $\operatorname{span}w=W$, (2) $\operatorname{card}w=\rho$, and
(3) $[w]$ contains no doubles.

Note that a spanning set $w$ is just a certain choice of basis for $W$. We let
$\operatorname{max}(v)$ denote the set of all spanning sets in $(v\cup
v^{\prime})\cap W$ and let $w^{c}\subset v$ denote the complement of $w$.

######  Lemma 11.

Let $v\in V_{r}^{\rho}$ and let $w\in\operatorname{max}(v)$. If $i\in[w^{c}]$,
then $P_{i}\subset\operatorname{span}w$.

###### Proof.

Suppose $P_{i}\not\subset\operatorname{span}w$. Then
$v_{i}^{\prime}\not\in\operatorname{span}w$, or
$v_{i}\not\in\operatorname{span}w$. WLOG, let
$v_{i}^{\prime}\not\in\operatorname{span}w$ and note that, under this
assumption, we have that $\dim\operatorname{span}w_{+i^{\prime}}=\rho+1$.
Recall that we assumed that the maximal dimension of any element in
$V_{r}^{\rho}$ is equal to $\rho$. Since $w_{+i^{\prime}}\in V_{r}^{\rho}$,
the equality $\dim\operatorname{span}w_{+i}=\rho+1$ contradicts the maximality
of $\rho$ and thus, $P_{i}\subset\operatorname{span}w$. ∎

######  Definition 12.

Let $v\in V_{r}^{\rho}$ and let $P_{i}\subset W=\operatorname{span}v$. We say
that a nonzero subset $\iota(P_{i})\subset v$ is a _minimal set for $P_{i}$_
if $P_{i}$ is contained in $\operatorname{span}\iota(P_{i})$ but there exists
no proper subset of $\iota(P_{i})$ whose span contains $P_{i}$.

If we fix $w\in\operatorname{max}(v)$, then there is a unique minimal set
$\iota(P_{i})\subset w$ for
$P_{i}=\operatorname{span}\\{v_{i},v_{i}^{\prime}\\}$ consisting of the unique
set of generators in $w$ for $v_{i}$ and $v_{i}^{\prime}$, e.g, if

$v_{i}^{\prime}=a_{1}v_{i_{1}}+\cdots+a_{\rho}v_{i_{\rho}},$

is the linear combination for $v_{i}^{\prime}$ in terms of the basis
$w=\\{v_{i_{1}},\dots,v_{i_{\rho}}\\}$, then
$v_{i_{1}},\dots,v_{i_{\rho}}\in\iota(P_{i})$, and similarly for $v_{i}$.

######  Lemma 13.

Fix $w\in\operatorname{max}(v)$ and let $i\in[w^{c}]$. Let
$\iota(P_{i})\subset W$ be minimal for $P_{i}$, then for each
$j\in[\iota(P_{i})]$, we have that $P_{j}\subset W$. Furthermore, if
$\iota(P_{j})$ is a minimal set for $P_{j}$ and $l\in[\iota(P_{j})]$, then
$P_{l}\subset W$. More generally, any $n$-th descendant of $P_{i}$ is
contained in $W$.

###### Proof.

By definition, if $v_{j}\in\iota(P_{i})$, then
$P_{i}\not\subset\operatorname{span}\iota(P_{i})_{-j}$ and so either $v_{i}$
or $v_{i}^{\prime}$ is _not_ in $\operatorname{span}\iota(P_{i})_{-j}$. WLOG,
let $v_{i}^{\prime}\not\in\iota(P_{i})_{-j}$ and note that this assumption
implies that $v_{i}^{\prime}\not\in\operatorname{span}w_{-j}$. Since
$v_{i}^{\prime}\not\in\operatorname{span}w_{-j}$, the set $w_{-j,+i^{\prime}}$
is a spanning set and since $j\in[w_{-j,+i^{\prime}}^{c}]$ we find that
$P_{j}\subset\operatorname{span}w_{-j,+i^{\prime}}\subset W$ by Lemma 11. For
the same $w\in\operatorname{max}(v)$ as above, let $\iota(P_{j})\subset w$ be
the minimal set for $P_{j}$ and let $l\in[\iota(P_{j})]$. We can use the same
proof to prove that $P_{l}\subset W$ if we can find a minimal set $B$ for
$P_{j}$ such that $l\in[B]$ and $j\not\in[B]$. Since
$P_{j}\subset\operatorname{span}w_{-j,+i^{\prime}}$, there exists a unique
$B\subset w_{-j,+i^{\prime}}$ minimal for $P_{j}$. For $l\neq j,i^{\prime}$,
note that $l\in[\iota(P_{j})]$ if and only if $l\in[B]$. since we already know
that $P_{j}$ and $P_{i}$ are contained in $W$, we may assume that $l\neq
j,i^{\prime}$. Now since $j\not\in[B]$, we can use the same proof as above to
conclude that $P_{l}\subset W$.

∎

Fix $w\in\operatorname{max}(v)$ and define

$[\Upsilon_{n}]=\\{\ j\ |\ v_{j}\in[\iota(P_{i})],\ \text{such that}\
\iota(P_{i})\subset w\ \text{and}\ i\in\Upsilon_{n-1}\\}\cup\Upsilon_{n-1},$

where

$[\Upsilon_{0}]=\\{\ j\ |\ v_{j}\in[\iota(P_{i})],\ \text{such that}\
\iota(P_{i})\subset w\ \text{and}\ i\in[w^{c}],\\}$

is the “initial condition.” Note that we meant what we wrote when we defined
$[\Upsilon_{0}]$, i.e., we do _not_ want $\Upsilon_{0}=[w^{c}]$. Clearly, the
sets $\Upsilon_{1},\Upsilon_{2},\dots$ will eventually converge to some subset
$\Upsilon$ of $w$. That is, there exists some $n_{0}>>0$ such that
$\Upsilon_{n_{0}}=\Upsilon_{n_{0}+i_{0}},\forall i_{0}\geq 0$ and so
$\Upsilon=\Upsilon_{n_{0}}$.

######  Lemma 14.

If $j\in[\Upsilon]$, then $P_{j}\subset\operatorname{span}\Upsilon$.

###### Proof.

If $j\in[\Upsilon]$, then $j\in[\Upsilon_{n}]$ for some $n$. By construction,
$\iota(P_{j})\subset\Upsilon_{n+1}$ and so
$P_{j}\subset\operatorname{span}\Upsilon_{n+1}$. In particular,
$P_{j}\subset\operatorname{span}\Upsilon$ since
$\Upsilon_{n+1}\subset\Upsilon$.

∎

The proof of Theorem 9 is now trivial:

###### Proof of Theorem 9.

Let

$[\Upsilon]=\\{j_{1},\dots,j_{h}\\}$

and note that $\operatorname{dim}\operatorname{span}\Upsilon\leq h$ and that
$\operatorname{span}\Upsilon$ contains at least $h$ planes, i.e,
$P_{j_{1}},\dots,P_{j_{h}}\subset\operatorname{span}\Upsilon$ by Lemma 14.
Recall that by construction $j_{1},\dots,j_{h}\in[w]$, but
$P_{i}\subset\operatorname{span}\Upsilon$ for all $i\in[w^{c}]$ since we
defined $\bigcup_{i\in w^{c}}\iota(P_{i})=\Upsilon_{0}$. Since
$\operatorname{card}w^{c}\geq 1$, the set
$\\{P_{j_{1}},\dots,P_{j_{h}}\\}\cup\bigcup_{i\in w^{c}}P_{i}$ satisfies the
conclusion of Theorem 9. ∎

##  Appendix B Odd periods

In [W15], the odd periods are defined to be restrictions of elements of
$\Omega$ to the components of the Ramond divisor, see Section B.2. Throughout
the paper, we computed the odd periods as values of elements in $H^{0}(C,L)$
along the points of $D$. In this section of the appendix we prove (Lemma 17)
that these values agree with the values of the odd periods as defined by
Witten . Throughout this section we consider only split super Riemann
surfaces, e.g. ones defined over a point or more generally over a (bosonic)
base.

###  B.1.  $\Omega_{1}=H^{0}(C,L)$

We begin by recalling Witten’s argument that the odd elements of $\Omega$ can
be identified with global sections of $L$. Let
$\operatorname{Ber}_{X}^{\prime}$ denote the subsheaf of
$\operatorname{Ber}_{X}(R)$ whose sections have a $\theta$-independent residue
along $R$. One can check that

(9) $\operatorname{Ber}_{X}^{\prime}=\Omega_{C}^{1}\oplus\Pi L.$

In contrast to an ordinary curve, the dualizing sheaf $\operatorname{Ber}_{X}$
on the supercurve $X$ is not equal to the cotangent bundle $\Omega_{X}^{1}$:
this is of course obvious given that $\operatorname{Ber}_{X}$ has rank $0|1$
while $\Omega_{X}^{1}$ has rank $1|1$. However, there is a canonical
isomorphism $\alpha:H^{0}(X,\operatorname{Ber}_{X}^{\prime})\cong\Omega$ which
in local superconformal coordinates $(z,\theta)$ sends $s=(g(z)+\theta
f(z))[dz|d\theta]$ to

(10) $\alpha(s)=g(z)d\theta+(f(z)+\theta g^{\prime}(z))dz.$

Composing $\alpha^{-1}$ with (9) we find that
$\Omega=H^{0}(C,\Omega_{C}^{1})\oplus\Pi H^{0}(C,L)$ (we will write $=$
whenever an isomorphism is canonical), and so
$\operatorname{dim}\Omega=\operatorname{dim}H^{0}(X,\operatorname{Ber}_{X}^{\prime})=g|h^{0}(L)$.

###  B.2. Odd periods

One can check that around each $R_{i}\in R$, there exist superconformal
coordinates $(z,\theta)$ such that $R_{i}$ is defined by the equation $z=0$.
In these coordinates a general element of $\Omega$ is of the form
$f(z,\theta)dz+g(z,\theta)d\theta$ and its restriction to $R_{i}$ (to compute:
set $z=0$ and apply the condition $ds=0$) can be shown to be equal to

$g(0)\ d\theta,$

where the constant $g(0)\in{\mathbb{C}}$ is unique up to a superconformal
change of coordinates $\theta\mapsto-\theta$. Let
$\sigma=\\{\sigma_{i}\in\\{+,-\\}\\}$ denote a choice of $\theta$ or $-\theta$
for each $R_{i}\in R$.

######  Definition 15.

For each $R_{i}\in R$, the _odd period $\mathfrak{o}_{i}^{\sigma}(s)$ of
$s\in\Omega$ with respect $\sigma$ _ is the constant $\sigma_{i}\cdot g(0)$.

######  Lemma 16.

The odd (resp. even) periods of even (resp. odd) elements in $\Omega$ are
equal to zero.

###### Proof.

If $s$ is an even element in $\Omega$, then its restriction to $R_{i}$ must
also be even; however $g_{0}d\theta$ is odd for all non-zero $g_{0}$. Thus
$g_{0}=0$.

Now let $s\in\Omega_{1}$ be odd and let
$\\{a_{1},\dots,a_{g},b_{1},\dots,b_{g}\\}$ be A- and B-cycles on $C$. The
even (A- and B-) periods of $s$ are the integrals

$\int_{a_{i}}\tau^{*}(s),\ \ \int_{b_{i}}\tau^{*}(s)$

where $\tau:C\subset X$ is the canonical inclusion. In local superconformal
coordinates $(z,\theta)$, the odd closed one form $s$ can be written as
$s=f(z)\theta dz+g(z)d\theta$. The pullback $\tau^{*}(s)$ is equal to $0$
since $\tau^{*}(\theta)=0$ and thus both of the above integrals are equal to
zero. ∎

We will now show that odd periods can be computed as residues along $D$. We
need to explain what we mean by a residue of an element in $\Omega$. Using the
identification $\Omega=H^{0}(X,\operatorname{Ber}_{X}^{\prime})$ and the
splitting (9) we can write every element $s$ in $\Omega$ as a sum
$s_{0}+s_{1}$ of a holomorphic differential $s_{0}$ on $C$ and a global
section $s_{1}$ of the twisted spin structure $L$. For each $p_{i}\in D$, we
define

$\operatorname{res}_{i}(s)=\operatorname{res}_{p_{i}}(s_{0})+\operatorname{res}_{p_{i}}(s_{1})=\operatorname{res}_{p_{i}}(s_{1})$

where $\operatorname{res}_{p_{i}}(s_{0})=0$ because $s_{0}$ is holomorphic
near $p_{i}$. Since $s_{1}$ is not a differential, we need to explain what we
mean by $\operatorname{res}_{p_{i}}(s_{1})$: Locally near each $p_{i}\in D$,
there exists coordinates $U(z)$ such that $\Omega_{U}^{1}\cong\mathcal{O}_{U}$
and such that $p_{i}$ is defined by $z=0$. Let $dz$ be a generator for
$\Omega_{U}^{1}$. We then have two possible generators for $L|_{U}$, namely
$\theta:=\sqrt{dz/z}$ and $-\theta:=-\sqrt{dz/z}$. For a fixed generator
$\theta$, there is a canonical isomorphism
$L|_{U}=\Omega_{U}^{1}(D)=\Omega_{U}^{1}(p_{i})$ sending $\theta$ to $dz/z$.
Let $\sigma=\\{\sigma_{i}\in\\{+,-\\}\\}_{i=1}^{2r}$ denote the sign of the
generator for each $L|_{U_{i}}$. Then each $s\in\Omega_{1}$ is locally of the
the form $\sigma_{i}\cdot g(z)\theta$ and is sent to
$\sigma_{i}\cdot\frac{g(z)}{z}dz$ in $\Omega_{U}^{1}(p_{i})$: And,

(11)
$\operatorname{res}_{p_{i}}(s_{1})=\sigma_{i}\cdot\operatorname{res}_{p_{i}}\left(\frac{g(z)}{z}dz\right).$

######  Lemma 17.

Let $\sigma$ be fixed and let $s\in\Omega$. Then for each point $p_{i}\in D$,
we have

$(1/2\pi i)\mathfrak{o}_{p_{i}}(s)=\operatorname{res}_{p_{i}}(s).$

###### Proof.

The lemma is obvious for even elements $\Omega$ so let $s$ be an odd element
and let $U(z,\sigma\cdot\theta)$ be superconformal coordinates near $R_{i}$ so
that $R_{i}$ is defined by $z=0$. Then there exist unique $g$ such that

$s|_{U}=\sigma_{i}\cdot(g(z)d\theta+g^{\prime}(z)\theta dz)$

and $\alpha^{-1}(s)=(\sigma_{i}\cdot g(z))[dz|d\theta]$. The odd periods are
now easy to read off: namely, $\mathfrak{o}_{i}^{\sigma}(s)=\sigma_{i}\cdot
g(0)$.

The section $\alpha^{-1}(s)$ is identified with the section $\sigma_{i}\cdot
g(z)\theta$ of $L|_{U}$ by (9) where $\theta=\sqrt{\frac{dz}{z}}$. The section
$\sigma_{i}\cdot g(z)\theta$ is in turn identified with the section
$\sigma_{i}\cdot\frac{g(z)}{z}dz$ of $\Omega_{U}^{1}(p_{i})$ under the
identification $L|_{U}=\Omega_{U}^{1}(p_{i})$. It is now easy to see that

$2\pi i\operatorname{res}_{z=0}\left(\frac{g(z)}{z}dz\right)=\sigma_{i}\cdot
g(0)=\mathfrak{o}_{i}^{\sigma}(s)$

.

∎

## References

  * [W13] Witten, E., “Notes On Holomorphic String And Superstring Theory Measures Of Low Genus,” arXiv:1306.3621.
  * [W15] Witten, E., “The Super Period Matrix With Ramond Punctures,” arXiv:1501.02499.
  * [D22] Donagi, N., “Valid Widgets Contain Legal Subwidgets,” to appear.
  * [A71] Atiyah, Michael Francis (1971),“Riemann surfaces and spin structures”, Annales Scientifiques de l’École Normale Supérieure, Série 4, 4: 47–62, ISSN 0012-9593, MR 0286136.
  * [M71] Mumford, David (1971), “Theta characteristics of an algebraic curve”, Annales Scientifiques de l’École Normale Supérieure, Série 4, 4 (2): 181–192, MR 0292836.
  * [ACGH] Arbarello, Enrico; Cornalba, Maurizio; Griffiths, Philip A.; Harris, Joe (1985), “The Basic Results of the Brill–Noether Theory”. Geometry of Algebraic Curves. Grundlehren der Mathematischen Wissenschaften 267. Vol. I. pp. 203–224. doi:10.1007/978-1-4757-5323-3-5. ISBN 0-387-90997-4.
  * [BN] von Brill, Alexander; Noether, Max (1874). ”Ueber die algebraischen Functionen und ihre Anwendung in der Geometrie”. Mathematische Annalen. 7 (2): 269–316. doi:10.1007/BF02104804. JFM 06.0251.01.


# On the use of graph neural networks and shape-function-based gradient
computation in the deep energy method

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

A graph neural network (GCN) is employed in the deep energy method (DEM) model
to solve the momentum balance equation in 3D for the deformation of linear
elastic and hyperelastic materials due to its ability to handle irregular
domains over the traditional DEM method based on a multilayer perceptron (MLP)
network. Its accuracy and solution time are compared to the DEM model based on
a MLP network. We demonstrate that the GCN-based model delivers similar
accuracy while having a shorter run time through numerical examples. Two
different spatial gradient computation techniques, one based on automatic
differentiation (AD) and the other based on shape function (SF) gradients, are
also accessed. We provide a simple example to demonstrate the strain
localization instability associated with the AD-based gradient computation and
show that the instability exists in more general cases by four numerical
examples. The SF-based gradient computation is shown to be more robust and
delivers an accurate solution even at severe deformations. Therefore, the
combination of the GCN-based DEM model and SF-based gradient computation is
potentially a promising candidate for solving problems involving severe
material and geometric nonlinearities.

###### keywords:

Automatic differentiation , Elasticity , Hyperelasticity , Partial
differential equations , Physics-informed neural networks

††journal: arXiv

*[inlinelist,1]label=(),

##  1 Introduction

The deep energy method (DEM) is a physics-informed neural network (PINN) model
developed by Nguyen-Thanh et al. [1]. The method readily applies to
engineering systems governed by an energy functional, whose solution coincides
with the stationary point of the functional. The DEM model takes a series of
points in the simulation domain as inputs and predicts field variables like
displacements [2, 3, 4] (for mechanics), temperature (for heat transfer),
and/or electric potential (for piezoelectricity) [4] at the same nodal
locations. The values of the predicted field variables are then used to
numerically compute the energy functional, which is defined as the loss
function. Standard optimization packages can be applied to minimize the loss,
and the weights and biases of the DEM model are updated through
backpropagation.

Current DEM model implementations are mostly based on a multilayer perceptron
(MLP) network [5, 6], which consists of multiple layers of fully connected
neurons. While shown to be highly effective, it is of interest to study how
the underlying neural network (NN) structure affects the performance of the
resulting DEM model. Chadha et al. [7] studied the effects of changing network
architecture (number of neurons and hidden layers) using an MLP-based DEM
model. Zhuang et al. [8] investigated the use of autoencoders in the DEM
method. While outside of the DEM framework, graph convolutional networks (GCN)
have seen wide applications in the field of engineering mechanics, ranging
from mechanics of polycrystals [9, 10], the solution to PDEs [11, 12, 13],
fluid mechanics [14, 15, 16], and crack propagation [17]. GCNs are
particularly suitable for describing data with structure and a high sense of
locality [11, 9]. In DEM, the computational domain is discretized into a
series of structured or unstructured nodes, which renders it perfect to be
represented by a graph weighted by the Euclidean distance between the nodes.
In existing DEM implementations based on MLP networks, most utilized a
structured grid to discretized the domain. As shown in [11], PINN models that
are based on convolutional neural networks require a fine grid to resolve the
solution, and the convolution operator limited it to structured grid. While
the GCN-based frameworks can represent unstructured mesh in terms of graphs
and apply a graph convolution to it, greatly reducing the burden to generate a
structured discretization of the domain. Therefore, in this work, we
investigate the use of GCNs as the underlying NN for the DEM method, and study
how this affects solution accuracy and time.

Many PINNs, including the DEM model, contain differential operators in the
definition of their loss functions. For PINNs that are based on the strong
form of the governing differential equation, such as the deep collocation
method [18, 19, 20, 21], second-order (e.g., in elasticity, heat transfer and
the Navier-Stokes equations) or higher-order (e.g., plate bending) spatial
gradients of field variables are required. For PINNs that are based on a
variational formulation, such as the deep Ritz method [22, 23], deep Galerkin
method [11], and DEM, first-order spatial gradients of field variables are
required. Besides, researchers have used a mixed formulation, combining both
the energy method and the strong form, to capture the mechanical response with
high solution gradients and stress concentrations [24, 25, 26].

In almost all cases, these spatial gradients are obtained utilizing automatic
differentiation (AD) of the NN model [27, 28, 29]. This approach is widely
used since AD is already applied in the backpropagation step during NN model
training. Most importantly, it alleviates the need to form an element-based
discretization of the computational domain as is typically required for the
finite element method (FEM), yielding a meshless simulation technique [19].
However, the cost of AD can be quite expensive, especially when many nodes are
employed in the domain. Therefore, as an alternative, spatial gradients of the
field variables can be calculated on a mesh-based discretization of the domain
using Sobel filters (finite difference) [30, 31, 32, 32, 33]. In addition,
finite element shape functions (SFs) can be used [7, 34, 11, 35]. In this
case, the gradients are evaluated at the integration points of the finite
elements through the shape function gradients, identical to the treatment in
classical FEM [36, 37]. We note that a thorough comparison between the AD- and
SF-based gradient computations and their implications on the accuracy and
stability of the DEM model has been missing from the literature. Thus, the
second objective of this work is to perform such a comparison using case
studies. More importantly, we point out that the point-wise AD-based gradient
computation is inherently susceptible to instability within the DEM framework,
which is a novel finding in the literature.

This paper is organized as follows: Section 2 presents an overview of graph
convolutional networks, the deep energy method, and a discussion on
instability. Section 3 presents and discusses the results of four numerical
examples. Section 4 summarizes the outcomes and highlights possible future
works. In the following discussion, we denote the DEM model based on MLP
networks as MLP-DEM and that based on GCNs as GCN-DEM.

##  2 Methods

###  2.1 Graph convolutional network

A GCN contains a graph that is defined by a set of nodes and the edges between
them. For computational mechanics problems, where geometries are typically
discretized into meshes, the nodes of domain discretization is a natural
choice for constructing a graph. The pair-wise Euclidean distance between the
nodes becomes the weight of each edge in the graph. The graph’s locality, or
sparseness, is ensured by removing edges whose weights exceed a threshold
radius $r$. Therefore, only nodes within a distance $r$ from each other remain
connected in the graph. An example graph for a 1-by-1-by-1 cube discretized by
27 nodes, forming a 3-by-3-by-3 grid is shown in Fig. 1; a threshold radius
$r=1/3$ was used to generate the graph. In this work, the GCN consists of an
input graph of nodal coordinates $\bm{X}$, which is mapped to an output graph
of nodal displacement vectors $\bm{u}$. The threshold distance is set to
$r=L_{x}/(N_{x}-1)$, where $L_{x}$ and $N_{x}$ denote the length and number of
nodes along the X-axis.

Figure 1: A graph of a cube with 27 nodes and 54 edges. Note that only nodes
within a distance of $1/3$ to each other remain connected in the graph.

The graph convolution is a key operation in GCN, whose implementation in
PyTorch Geometric [38] is based on the Chebyshev spectral graph convolution
operator [39]. The message passing function is defined as [39]:

$\mathbf{X}^{i}=f_{\rm{act}}\left(\sum_{k=1}^{K}\mathbf{Z}^{k}\cdot\mathbf{\Theta}^{k}\right),$
(1)

where $i$ denotes the $i^{th}$ layer, K is the polynomial order, $\bm{\Theta}$
denotes the trainable parameters, and the $k^{th}$ basis vectors
$\bm{Z}^{(k)}$ is defined recursively as [39, 11]:

$\displaystyle\bm{Z}^{i-1,1}=\bm{X}^{i-1},$ (2)
$\displaystyle\bm{Z}^{i-1,2}=\bm{\hat{L}}\cdot\bm{X}^{i-1},$
$\displaystyle\bm{Z}^{i-1,k}=2\cdot\bm{\hat{L}}\cdot\bm{Z}^{i-1,k-1}-\mathbf{Z}^{i-1,k-2},$

and:

$\displaystyle\bm{\hat{L}}=\bm{L}-\bm{I},$ (3)
$\displaystyle\bm{L}=\bm{I}-\bm{D}^{-\frac{1}{2}}\bm{A}\bm{D}^{-\frac{1}{2}}.$

$\bm{I}$, $\bm{A}$, and $\bm{D}$ denote the identity, adjacency, and degree
matrix of the graph, respectively. The hyperbolic tangent (tanh) function is
chosen to be the nonlinear activation function for the Chebyshev convolution
layers, and the maximum polynomial degree K=1 is used in this work. Training
of the GCN refers to the iterative process where the parameters $\bm{\Theta}$
of each layer are updated by gradient descent [40] using the L-BFGS algorithm
[41].

###  2.2 The deep energy method based on graph convolution network

In this section, we describe the deep energy method based on graph convolution
network (GCN-DEM) in the context of the elasticity equation in structural
mechanics. This work considers two material models: linear elastic material in
small deformation and Neo-Hookean material in finite deformation. In the
absence of any body and inertial forces, the equilibrium equations and
boundary conditions (Dirichlet and Neumann) under small deformation can be
stated in terms of the Cauchy stress tensor $\bm{\sigma}$ as

$\displaystyle\nabla_{\bm{x}}\cdot\bm{\sigma}=\bm{0},\;\;\forall\bm{X}\in\Omega,$
(4) $\displaystyle\bm{u}=\bar{\bm{u}},\;\;\forall\bm{X}\in\partial\Omega_{u},$
$\displaystyle\bm{\sigma}\cdot\bm{n}=\bar{\bm{t}},\;\;\forall\bm{X}\in\partial\Omega_{t},$

where $\bm{n}$, $\bar{\bm{u}}$ and $\bar{\bm{t}}$ denote the outward boundary
normal, prescribed displacement, and prescribed traction, respectively.
$\nabla_{\bm{x}}$ denotes the gradient operator in the current configuration.
In the small deformation setting, the strain tensor is given by:

$\bm{\epsilon}=\frac{1}{2}(\nabla_{\bm{x}}\bm{u}+\nabla_{\bm{x}}\bm{u}^{T}).$
(5)

For linear elastic materials, the stress can be computed from the constitutive
law as:

$\bm{\sigma}=\frac{E}{1+\nu}\bm{\epsilon}+\frac{E\nu}{(1+\nu)(1-2\nu)}tr(\bm{\epsilon})\bm{I},$
(6)

where $E$ and $\nu$ are Young’s modulus and Poisson’s ratio. The strain energy
density $\Psi_{LE}$ in the linear elastic case is given by:

$\Psi_{LE}=\frac{1}{2}\bm{\sigma}:\bm{\epsilon}.$ (7)

In the finite deformation setting, the equilibrium equations and boundary
conditions (Dirichlet and Neumann) can be stated in terms of the first Piola-
Kirchhoff stress tensor $\bm{P}$ as:

$\displaystyle\nabla_{\bm{X}}\cdot\bm{P}=\bm{0},\;\;\forall\bm{X}\in\Omega,$
(8) $\displaystyle\bm{u}=\bar{\bm{u}},\;\;\forall\bm{X}\in\partial\Omega_{u},$
$\displaystyle\bm{P}\cdot\bm{N}=\bar{\bm{t}},\;\;\forall\bm{X}\in\partial\Omega_{t},$

where $\nabla_{\bm{X}}$ and $\bm{N}$ denote the gradient operator and outward
normal in the reference configuration, respectively. $\bm{P}$ is power
conjugate to the deformation gradient tensor:

$\bm{F}=\nabla_{\bm{X}}\bm{u}+\bm{I},$ (9)

and can be computed from the hyperelastic strain energy density $\Psi_{NH}$
as:

$\bm{P}=\frac{\partial\Psi_{NH}}{\partial\bm{F}}.$ (10)

For Neo-Hookean material, $\Psi_{NH}$ is given by [42]:

$\Psi_{NH}=C_{10}[tr(\bar{\bm{F}}\cdot\bar{\bm{F}}^{T})-3]+\frac{1}{D_{1}}[{\rm{det}}(\bm{F})-1]^{2},$
(11)

where:

$\bar{\bm{F}}={\rm{det}}(\bm{F})^{-\frac{1}{3}}\bm{F},$ (12)

is the deviatoric part of the deformation gradient, and $C_{10}$ and $D_{1}$
are material constants.

In either case, the DEM model seeks the solution to the equilibrium equations
via the principle of minimum potential energy (PMPE). For a body in static
equilibrium with no applied body forces, the potential energy of the system
reads:

$\psi(\bm{u})=\int_{\Omega}\Psi\,dV-\int_{\partial\Omega_{t}}\bar{\bm{t}}\cdot\bm{u}\,dA.$
(13)

The loss function ($\mathcal{L}$) in GCN-DEM is defined identically as the
potential of the system:

$\mathcal{L}(\bm{u})=\psi(\bm{u}).$ (14)

The Neumann boundary conditions are enforced by the boundary integral part of
Eq. (13), and Dirichlet boundary conditions are enforced directly as in the
work of He et al. [34] to avoid modification of the loss function. The
solution to the elasticity problem, as given by the GCN-DEM model, is defined
as the stationary point of the potential energy functional:

$\bm{u}^{*}={\rm{arg}}\min_{\bm{u}}\,\mathcal{L}(\bm{u}).$ (15)

###  2.3 Spatial gradient computation and instability

A central component in calculating the loss function value is the evaluation
of the spatial gradients of $\bm{u}$. In many previous studies [1, 4, 19, 2],
the spatial gradients of $\bm{u}$ are evaluated at discrete points inside the
domain through AD of the underlying NN model:

$\frac{\partial\bm{u}(\bm{X})}{\partial\bm{X}}=\frac{\partial\bm{u}}{\partial\bm{z}^{n}}\cdot\frac{\partial\bm{z}^{n}}{\partial\bm{z}^{n-1}}\dots\cdot\frac{\partial\bm{z}^{2}}{\partial\bm{z}^{1}}\cdot\frac{\partial\bm{z}^{1}}{\partial\bm{X}},$
(16)

where $\bm{z}^{i}$ denotes the output of the $i^{th}$ layer of the NN model.
The first and last layers are the input and output layer, respectively. After
calculating gradients at the nodes $\bm{X}$ in the domain, Simpson’s rule or
trapezoidal rule is typically used to perform the integration. The AD-based
approach is straightforward to implement, as AD is already used in the
training of the NN model during backpropagation. This point-based method also
alleviates the need to form a mesh of the computational domain like in FEM.

The studies by Chadha et al. [7] and He et al. [34] demonstrated gradient
computation and numerical integration through finite element SFs and Gauss
quadrature. In this case, the spatial gradients are given by:

$\frac{\partial\bm{u}(\bm{X})}{\partial\bm{X}}\approx\frac{\partial\bm{\phi}}{\partial\bm{\xi}}\bm{J}^{-1}\cdot\bm{u},$
(17)

where $\bm{\phi}$ , $\bm{\xi}$, and $\bm{u}$ denote the finite element SFs,
natural coordinates, and displacement vector (evaluated at discrete nodes),
respectively. $\bm{J}=\frac{\partial\bm{X}}{\partial\bm{\xi}}$ denotes the
Jacobian matrix of the isoparametric mapping of the finite element. Instead of
evaluating the gradients at the nodes, the gradients are evaluated at the
quadrature integration points of the ’elements’, reminiscent to classical FEM
procedures. Therefore, the SF-based approach requires the formation of
isoparametric elements from the nodes in the computational domain.

We shall demonstrate with a simple 1D example that, despite being
straightforward to implement, the AD-based spatial gradient computation can be
susceptible to instability that leads to divergence of the DEM method.
Consider a 1D bar whose dimensions and properties are all unity, subject to a
unit tensile force at its right end. The displacement is simply $u(x)=x$. Let
the bar be discretized into two nodes located at $x=0$ and $x=1$, and let
$u_{p}(x)$ be a perturbation to the displacement field given by:

$u_{p}(x)=\frac{\Delta u}{2}({\rm{tanh}}[20(x-0.5)]+1),$ (18)

where $\Delta u$ characterizes the magnitude of the perturbation111The factor
20 is chosen randomly and does not affect the argument we want to make.. This
perturbation can be thought of as a strain localization at the center of the
bar. The true displacement $u$ and the perturbed displacement
$\tilde{u}=u+u_{p}$ are plotted in Fig. 2a.

(a)

(b)

Figure 2: 1D bar example: LABEL:sub@theory_disp True and perturbed
displacements. $\Delta u$ is taken to be 0.5 in the plot.
LABEL:sub@theory_energy Comparison of system potential energy at different
$\Delta u$

First, we use AD (in this case, it corresponds to analytical differentiation)
and trapezoidal rule to evaluate Eq. (13) using the two nodes. Then, we use SF
to approximate the displacement gradient and use one-point Gauss quadrature to
evaluate Eq. (13). We repeat the process for multiple values of $\Delta u$ and
plot the results in Fig. 2b. At zero perturbation, both methods yield the true
potential energy of the system, which is -0.5. However, as the perturbation
increases, the integrated potential calculated from AD keeps decreasing,
thereby violating PMPE (the unperturbed equilibrium state has minimum
potential). While that calculated from the SF gradient remains larger than the
true minimum and does not violate PMPE. This behavior is expected, as the tanh
function has vanishing gradient away from its center. Therefore the point-wise
AD fails to detect strain localization that occurs in between the nodes, while
the boundary integral part of Eq. (13) keeps increasing with increased
perturbation. On the contrary, strain localization can always be detected by
an element-based SF gradient, as the gradient is computed from the difference
of the displacements of the two nodes. Therefore, we argue that the SF-based
gradient computation is more stable than its AD-based counterpart with respect
to strain localization.

The inability of AD to detect and penalize strain localization between nodes
can adversely affect the solution of energy-based methods, as the loss
function has a lower value than the true equilibrium solution when strain
localization occurs, thus driving the solution to the incorrect direction of
increasing strain localization. Although the argument is presented using a
simple example, it can be shown easily that the argument generalizes to
arbitrary geometries with arbitrary spacing between the nodes. In Section 3,
we demonstrate that this type of instability is indeed observed when AD-based
gradient computation is used, especially when a large load is applied.

##  3 Results and discussion

In this section, we critically compare GCN-DEM and MLP-DEM using four case
studies. In addition, we compare the solution stability of AD- and SF-based
gradient computations. In all cases except Section 3.3, the geometry
considered is a cantilever beam of dimensions 4-by-1-by-1 units, subject to
downward traction at its right surface. In the first two cases, 3700 nodes
were placed in the domain, forming a 37-by-10-by-10 grid. Although a
structured hexahedral grid was used in all the examples, we highlight that the
GCN-based DEM can be easily extended to unstructured tetrahedral (or
triangular for 2D) meshes as in the work of Fuhg et al. [24]. In the first
case, we consider a linear elastic material with $E=1000$ and $\nu=0.3$,
subjected to sequentially increasing loads. In the second case, we consider a
Neo-Hookean material with $C_{10}=192.31$ and $D_{1}=0.0024$, subjected to the
same loads as in case 1. Then, we demonstrate that instability can occur in
other AD-based DEM implementations, not limited to the ones presented in our
work. In the last case, we investigate how grid refinement can remedy the
instabilities caused by AD-based gradient computation. In all cases except
Section 3.3, we compared the displacement results with those obtained from FEM
using Abaqus/Standard [42]. The GCN-DEM and MLP-DEM models were implemented in
PyTorch (version 1.11.0) [43]. The MLP-DEM implementation was adopted from
Abueidda et al. [2]. All training of the NNs was done on an Intel i7-11800H
processor. For a fair and consistent comparison, the GNN in GCN-DEM and the
MLP model in MLP-DEM share the same network structure: they have 6 layers
(including input and output). The number of neurons in each layer is 3, 16,
32, 64, 32 , 16, and 3, respectively. The hyperbolic tangent function was used
as activation function for all layers except the output, which has linear
activation. The L-BFGS optimizer [41] with a fixed learning rate of 0.01 is
used to train the models. Training process is stopped when a maximum of 20
training iterations is reached, or when the relative change in loss function
value is less than $5\times 10^{-5}$. To test the robustness of both methods,
we applied the full magnitude of the external load in a single load step,
which is in contrast to FEM, where large loads are applied gradually
throughout several load steps for better convergence. This approach puts the
stability and robustness of the methods at severe deformation to test.

###  3.1 Linear elastic material

In this case, we gradually increased the magnitude of the applied traction in
six different simulations. The applied loads were: $t=$ -2.5, -5, -7.5, -10,
-15, and -25, respectively. To quantify the model accuracy, we compare the
GCN-/MLP-DEM solutions to FEM solutions generated using identical node
layouts. The relative difference in displacement is computed as:

$RD_{i}=\frac{|\bm{u}^{i}_{NN}-\bm{u}^{i}_{FE}|}{{\rm{max}}(|\bm{u}^{i}_{FE}|)}\times
100\%,\;\;\;i=x,y,z.$ (19)

The mean relative difference (averaged over all nodes and all three
displacement components), final loss value, and train time are presented in
Table 1 and plotted in Fig. 3. Specifically, we highlighted the cases that
failed to converge in red. To give a more direct visualization of the
distribution of the displacement error, contour plots of the displacement
error for different methods at the case $t=-7.5$ are presented in Fig. 4 and
Fig. 5. The deformed shapes at $t=-25$ are presented in Fig. 6 to highlight
the occurrence of instability for the AD-based gradient computation.

Table 1: Performance comparison of GCN-DEM and MLP-DEM, linear elastic model Method | t = -2.5 | t = -5 | t = -7.5 | t = -10 | t = -15 | t = -25
---|---|---|---|---|---|---
| Mean percent difference (%)
GCN-DEM (AD) | 1.17 | 1.82 | 1.95 | 964.60 | 772.20 | 721.00
MLP-DEM (AD) | 1.64 | 4.17 | 2.16 | 2.38 | 3.18 | 604.60
GCN-DEM (SF) | 2.99 | 4.65 | 1.77 | 3.05 | 1.58 | 2.43
MLP-DEM (SF) | 2.79 | 6.29 | 2.91 | 4.99 | 3.87 | 2.09
| Final loss function value
GCN-DEM (AD) | -0.83 | -3.33 | -7.51 | -592.05 | -1014.66 | -2146.11
MLP-DEM (AD) | -0.83 | -3.33 | -7.56 | -13.32 | -29.83 | -1981.89
GCN-DEM (SF) | -0.81 | -3.24 | -7.31 | -12.99 | -29.26 | -81.13
MLP-DEM (SF) | -0.81 | -3.25 | -7.29 | -12.96 | -29.17 | -80.80
| Train time [s]
GCN-DEM (AD) | 36.08 | 31.17 | 37.50 | 117.08 | 117.80 | 112.70
MLP-DEM (AD) | 76.06 | 58.81 | 70.98 | 79.70 | 53.60 | 115.20
GCN-DEM (SF) | 45.70 | 42.80 | 55.20 | 47.00 | 46.70 | 34.90
MLP-DEM (SF) | 79.10 | 83.40 | 74.00 | 68.26 | 62.00 | 51.20

(a)

(b)

(c)

Figure 3: Comparing GCN-DEM and MLP-DEM, AD-based and SF-based gradient computation: LABEL:sub@fig:LE_RD Relative difference in displacements. LABEL:sub@fig:LE_loss Final loss value. LABEL:sub@fig:LE_time Training time for simulation. In LABEL:sub@fig:LE_RD and LABEL:sub@fig:LE_loss, simulations that failed to converge are marked with a $\bm{\times}$. (a) GCN, RDx=0.42% |  |  (b) GCN, RDy=0.15% |  |  (c) GCN, RDz=5.30% |
---|---|---|---|---|---
(d) MLP, RDx=0.50% |  (e) MLP, RDy=0.49% |  (f) MLP, RDz=5.50%
Figure 4: Comparison of relative displacement difference between GCN-DEM and MLP-DEM, using AD for gradient computation, linear elastic material. The mean relative difference for each displacement component is reported in the caption. (a) GCN, RDx=0.99% |  |  (b) GCN, RDy=1.06% |  |  (c) GCN, RDz=3.27% |
---|---|---|---|---|---
(d) MLP, RDx=1.18% |  (e) MLP, RDy=1.25% |  (f) MLP, RDz=6.30%
Figure 5: Comparison of relative displacement difference between GCN-DEM and
MLP-DEM, using FE shape functions for gradient computation, linear elastic
material. The mean relative difference for each displacement component is
reported in the caption.

(a) GCN-DEM, t = -25

(b) MLP-DEM, t = -25

Figure 6: Comparing deformed shapes computed from GCN-DEM and MLP-DEM:
LABEL:sub@fig:LE_def_G Using GCN-DEM. LABEL:sub@fig:LE_def_D Using MLP-DEM.
The deformed shape computed by SF-based gradient computation is shown in
translucent purple with a scaling factor of 0.15. $\epsilon_{11}$ contour is
overlaid onto the deformed shape obtained from AD-based method with a scaling
factor of 0.015. Note that for the plot to look reasonable, the scaling factor
in this figure is _ten_ times smaller for the AD cases.

From Table 1, we immediately see that the AD-based gradient computation leads
to instability in several cases, especially at loads greater than -10. In
contrast, the SF-based gradient computation leads to stable results in all
cases. Inspecting the final loss value for the cases that failed to converge,
we see that they are much smaller than the converged value for the same
loading, which is in complete agreement with our analysis presented in Section
2.3, and serves as a direct proof that this is the cause of instability. Even
for the cases where both AD- and SF-based formulation converged, we notice
that the AD-based methods always yield a lower loss value than their SF-based
counterparts, hinting that there exists a threshold load magnitude beyond
which divergence may occur. Further, when comparing the first three columns of
Table 1, we see that at small loads, the AD-based formulation outperforms the
SF-based formulation, yielding a shorter solution time and higher accuracy.
When the gradient computation method is fixed, we see that GCN-DEM outperforms
MLP-DEM in 4 out of the 6 cases, again showing higher accuracy with a shorter
run time.

Fig. 4 and Fig. 5 provide much more information regarding the performances of
GCN-DEM and MLP-DEM. Since the loading is in the XY plane and the Z planes are
unconstrained except at the root, stresses in the Z direction are expected to
be small and only account for a small percentage of the overall system strain
energy $\Psi$. In this sense, it is reasonable to expect that the relative
error in the Z direction to be higher than those in the X- and Y-directions,
as it has a smaller weight in the loss function. Comparing the two rows in
Fig. 4, we see that GCN-DEM outperforms MLP-DEM in this case, yielding lower
relative error in all three displacement components. When inspecting the SF-
based cases in Fig. 5, we again notice the same trend, where GCN-DEM provides
better results in all displacement components. Comparing Fig. 4 and Fig. 5, it
is visually apparent that SF-based formulation leads to a higher error than
the AD-based formulation in this case. However, we remark that the errors
remain below 1.5% for the SF cases in the X- and Y-components.

These observations indicate that when the applied load is small, the GCN-DEM
method based on AD gradient computation tends to provide the best performance
in terms of accuracy and run time. When the load magnitude increases, the GCN-
DEM method based on SF gradient computation offers a much more robust solution
in general.

Fig. 6 provides a closer look into the strain localization instability that
occurs with AD-based gradient computations. From the grids that are overlaid
onto the deformed shape of the AD-based simulations, we see that severe
deformation happened at the root of the beam, while to the right of the strain
localization, negligible deformation is seen. From a continuum mechanics
standpoint, such severe localization must be accompanied by high local strain
values, which are grossly neglected by the point-based AD gradient
computation. The observed displacement field of the AD simulations also
resembles the 1D example presented in Section 2.3, where large displacement
jumps occur in between two nodes of the domain. This behavior further
strengthens our argument and shows that the inability to detect strain
localization between the nodes is the root cause of instability.

###  3.2 Neo-Hookean material

In the previous case, linear elastic material under a small strain formulation
provides neither geometric nor material nonlinearity to the system. Therefore,
we devote this section to study how GCN-DEM and MLP-DEM compare when both
sources of nonlinearity are present. Same loads were applied as in Section
3.1, and the Neo-Hookean model in Abaqus/Standard was used to generate FEM
comparisons. The mean relative difference, final loss value, and train time
are presented in Table 2 and plotted in Fig. 3 for graphical visualization.
Cases that failed to converge are highlighted in red. Contour plots of the
displacement error at the case $t=-10$ are presented in Fig. 8 and Fig. 9. The
deformed shapes at $t=-25$ are presented in Fig. 10.

Table 2: Performance comparison of GCN-DEM and MLP-DEM, Neo-Hookean model Method | t = -2.5 | t = -5 | t = -7.5 | t = -10 | t = -15 | t = -25
---|---|---|---|---|---|---
| Mean percent difference (%)
GCN-DEM (AD) | 8.06 | 6.03 | 6.26 | 1.52 | 254.30 | 292.60
MLP-DEM (AD) | 14.17 | 5.56 | 5.23 | 3.76 | 361.10 | 811.70
GCN-DEM (SF) | 2.72 | 2.96 | 3.57 | 1.42 | 2.76 | 1.33
MLP-DEM (SF) | 4.10 | 3.11 | 3.03 | 2.88 | 1.13 | 2.39
| Final loss function value
GCN-DEM (AD) | -0.79 | -3.15 | -6.80 | -11.57 | -186.68 | -343.66
MLP-DEM (AD) | -0.71 | -3.08 | -6.80 | -11.55 | -177.65 | -413.19
GCN-DEM (SF) | -0.80 | -3.11 | -6.70 | -11.29 | -22.71 | -50.79
MLP-DEM (SF) | -0.79 | -3.10 | -6.69 | -11.30 | -22.69 | -50.71
| Train time [s]
GCN-DEM (AD) | 20.85 | 22.62 | 28.40 | 48.37 | 130.10 | 127.80
MLP-DEM (AD) | 21.26 | 35.89 | 63.22 | 59.05 | 133.10 | 141.40
GCN-DEM (SF) | 92.73 | 85.43 | 101.30 | 72.35 | 74.19 | 62.31
MLP-DEM (SF) | 147.10 | 140.00 | 116.90 | 141.10 | 136.90 | 98.80

(a)

(b)

(c)

Figure 7: Comparing GCN-DEM and MLP-DEM, AD-based and SF-based gradient computation: LABEL:sub@fig:NH_RD Relative difference in displacements. LABEL:sub@fig:NH_loss Final loss value. LABEL:sub@fig:NH_time Training time for simulation. In LABEL:sub@fig:NH_RD and LABEL:sub@fig:NH_loss, simulations that failed to converge are marked with a $\bm{\times}$. (a) GCN, RDx=0.77% |  |  (b) GCN, RDy=0.54% |  |  (c) GCN, RDz=3.25% |
---|---|---|---|---|---
(d) MLP, RDx=0.56% |  (e) MLP, RDy=0.40% |  (f) MLP, RDz=10.33%
Figure 8: Comparison of relative displacement difference between GCN-DEM and MLP-DEM, using AD for gradient computation, Neo-Hookean material. The mean relative difference for each displacement component is reported in the caption. (a) GCN, RDx=0.53% |  |  (b) GCN, RDy=0.47% |  |  (c) GCN, RDz=3.28% |
---|---|---|---|---|---
(d) MLP, RDx=0.52% |  (e) MLP, RDy=0.42% |  (f) MLP, RDz=7.71%
Figure 9: Comparison of relative displacement difference between GCN-DEM and
MLP-DEM, using FE shape functions for gradient computation, Neo-Hookean
material. The mean relative difference for each displacement component is
reported in the caption.

(a) GCN-DEM, t = -25

(b) MLP-DEM, t = -25

Figure 10: Comparing deformed shapes computed from GCN-DEM and MLP-DEM:
LABEL:sub@fig:NH_def_G Using GCN-DEM. LABEL:sub@fig:NH_def_D Using MLP-DEM.
The deformed shape computed by SF-based gradient computation is shown in
translucent purple with a scaling factor of 0.3. $\epsilon_{11}$ contour is
overlaid onto the deformed shape obtained from AD-based method with a scaling
factor of 0.03. Note that the scaling factor in this figure is _ten_ times
smaller for the AD cases.

Table 2 shows very similar trends compared to Table 1. AD-based formulations
become unstable when load magnitude is greater than 10, while SF-based
formulations remain stable for all loads tested. Comparing the four cases
where AD-based formulations were stable, we see that in the hyperelastic case,
GCN-DEM and MLP-DEM provide very similar accuracy, each having higher accuracy
in 2 of the 4 cases. While the accuracy is similar, GCN-DEM still provides a
shorter run time than MLP-DEM. When using SF-based gradient computation, GCN-
DEM generally delivers a more accurate solution than MLP-DEM, with a run time
that is always shorter than MLP-DEM. Interestingly, AD-based formulation no
longer holds an advantage on accuracy against SF-based formulation in small
loads as in Section 3.1 and instead delivers worse performance than SF-based
formulation in all cases.

Comparing AD-based GCN-DEM and MLP-DEM in Fig. 8, we see that MLP-DEM
outperforms GCN-DEM in the X- and Y-components, which are the primary
deformation modes. For the SF-based formulation in Fig. 9, we see that the
performance of GCN-DEM and MLP-DEM is again very comparable. The AD-based and
SF-based formulations provide very similar levels of accuracy in the
hyperelastic case, which is distinct from the observations made in Section
3.1. Fig. 10 demonstrates similar strain localization at the root of the beam
for AD-based simulations, indicating that this phenomenon is due to the
gradient computation method and is independent of the material model. AD-based
gradient computation again failed to capture the severe strain localization at
the beam’s root.

To conclude, in the nonlinear hyperelastic case, the solution accuracy of GCN-
DEM and MLP-DEM is very comparable. However, GCN-DEM is more computationally
effective and trains in a shorter time. SF-based gradient computation still
holds a massive advantage over AD-based methods in terms of stability and
robustness in severe deformation.

###  3.3 Occurrence of instability in the MLP-DEM implementation by Nguyen-
Thanh et al.

In this example, we demonstrate that the strain localization instability is
not limited to our own implementation of GCN-DEM and MLP-DEM, but also exists
in other AD-based MLP-DEM implementations. For this purpose, we used the MLP-
DEM code developed in the work of Nguyen-Thanh et al. [1]. The model was
implemented in PyTorch, and AD was used to evaluate the displacement gradients
at the nodes. The trapezoidal rule was used to perform the domain integration.
For simplicity, we used the 2D Neo-Hookean cantilever beam example, which has
a dimension of 4-by-1 units. Materials properties are as defined in the
example, and no modifications were made. To match similar gird spacing as used
in our examples, we changed the grid dimension to 37-by-10. Three different
downward loads were applied at the right edge; they are -1, -2.5, and -5. A
total of 200 iterations were conducted during the training of the MLP-DEM
model. Contour plots of the Y component of the displacements are shown in Fig.
11.

(a) t = -1

(b) t = -2.5

(c) t = -5

Figure 11: MLP-DEM simulations of a Neo-Hookean cantilever beam subject to
increasing downward load. A scaling factor of 1 is used to plot the deformed
shape in LABEL:sub@fig:L1 and LABEL:sub@fig:L2, while a scaling factor of 0.05
is used in LABEL:sub@fig:L3. Note the occurrence of strain localization
instability in the largest load case.

From Fig. 11, we see the occurrence of strain localization instability when
downward traction of -5 was applied. Comparing Fig. 11c and Fig. 10, we notice
very similar strain localization near the root of the beam. This result
demonstrates that strain localization instability is not limited to our
implementation of the GCN-DEM and MLP-DEM models, but also exists in other
implementations that use AD for displacement gradient computations. However,
we do point out that both frameworks tested in this work were implemented in
PyTorch, and the existence of similar AD-induced instability is not checked in
models implemented in TensorFlow [44], another commonly used package for
building NNs.

###  3.4 Effect of grid refinement

The three previous sections established that instabilities might occur for AD-
based formulations. Therefore, it is interesting to study how this instability
can be remedied besides applying the load gradually in multiple steps. As
analyzed in Section 2.3, the instability roots from possible strain
localization that happens between two nodes. It is therefore intuitive to
refine the grid for NN model training. We doubled the node count from 3700 to
7436 to form a 44-by-13-by-13 grid to solve the hyperelastic beam problem
subjected to a load $t=-15$. Since the AD-based MLP-DEM method could not
converge for this grid size, we crated another grid with 21708 nodes, forming
a 67-by-18-by-18 layout. The results are summarized in Table 3. For all cases,
finite element meshes with the same grid layouts were created for the FEM
comparison. In addition to testing our implementations, the MLP-DEM model
implemented by Nguyen-Thanh et al. [1] was also tested to see if grid
refinement has similar effects on the results. Section 3.3 shows that a
downward traction of -5 induced instability on a 37-by-10 grid. We performed
uniform grid refinement twice and presented the contour plots of the Y
displacements in Fig. 12.

Table 3: Effects of grid refinement, Neo-Hookean model, $t=-15$ Method |  37$\times$10$\times$10 |  44$\times$13$\times$13 |  67$\times$18$\times$18
---|---|---|---
| Mean percent difference (%)
GCN-DEM (AD) | 254.30 | 2.28 | /
MLP-DEM (AD) | 361.10 | 212.90 | 4.26
GCN-DEM (SF) | 2.76 | 2.28 | /
MLP-DEM (SF) | 1.13 | 3.90 | /
| Final loss function value
GCN-DEM (AD) | -186.68 | -23.07 | /
MLP-DEM (AD) | -177.65 | -111.91 | -22.81
GCN-DEM (SF) | -22.71 | -22.72 | /
MLP-DEM (SF) | -22.69 | -22.64 | /
| Train time [s]
GCN-DEM (AD) | 130.10 | 69.01 | /
MLP-DEM (AD) | 133.10 | 221.50 | 302.50
GCN-DEM (SF) | 74.19 | 103.70 | /
MLP-DEM (SF) | 136.90 | 230.50 | /

(a) 37$\times$10 grid

(b) 74$\times$20 grid

(c) 111$\times$30 grid

Figure 12: Effect of grid refinement on strain localization instability.
Scaling factors of 0.05, 0.6 and 1 are used to plot the deformed shape in
LABEL:sub@fig:L3M1, LABEL:sub@fig:L3M2 and LABEL:sub@fig:L3M3, respectively.

From Table 3, we see that GCN-DEM recovered from instability upon doubling the
node count, achieving an accuracy similar to that obtained by the SF-based
GCN-DEM in the original coarse grid. MLP-DEM did not recover from instability
upon doubling the node count and instead required an even finer mesh to remain
stable. We also note that SF-based simulations remained stable in all cases
and achieved similar accuracy compared to their AD-based counterparts. A
similar trend is also observed in Fig. 12, where we see that the first uniform
refinement (Fig. 12b) reduced the amount of strain localization while the
second uniform refinement (Fig. 12c) produced reasonable displacement values.
From the results, we conclude that AD-based gradient computation should be
used in conjunction with a fine grid during training to avoid strain
localization instability. If a coarse grid is desired, then AD-based gradient
computation may not be suitable, and SF-based gradient computation is
recommended for increased robustness.

##  4 Conclusions and future work

In this work, we developed a DEM model based on GCN, which is different from
the typical DEM models that use an MLP network. In addition, we implemented
and tested two different spatial gradient computation techniques, one using AD
and the other one using SF. Numerical tests were performed using two material
models, linear elastic and hyperelastic, to test the model performance in
linear and nonlinear settings. Special attention is paid to studying the
stability of the algorithm changes with increasing load magnitude, and how
different spatial gradient calculation methods affect the stability.

The current work is novel in two main aspects: to the best of the authors’
knowledge, this is the first time in the literature that a graph convolutional
network is employed in the deep energy method, and we present a comparison of
accuracy and computational time with the traditional DEM models based on
multilayer perceptron networks. In addition, we performed a critical
evaluation of two spatial gradient computation methods to investigate their
effects on solution accuracy, solution time and most importantly, solution
stability. Through two 3D examples using linear elastic and hyperelastic
material models, we show that our GCN-DEM model shows higher computational
efficiency than the MLP-DEM model while maintaining similar, if not better,
solution accuracy. This comparison highlights the important role of network
type in the energy-based models, sheds light on the future deployment of
different types of NNs in the deep energy framework, and shows a promising
direction to pursue to improve the quality of deep energy-based methods to
solve PDEs.

Another key finding of this current study is recognizing that the AD-based
spatial gradient computation can lead to instability in the DEM framework due
to its inability to detect erroneous strain localization. This result means
that the AD-based DEM framework might fail to converge to the actual solution
if the applied load in a single load step is large or when a coarse mesh is
used. On the contrary, SF-based gradient computation demonstrated superior
stability, especially in large deformations, while delivering a solution
accuracy similar to that of AD-based DEM (when the AD-based one is stable).
This robustness is highlighted by the fact that the SF-based DEM remained
stable and gave accurate solutions in all examples presented in this work.
Although the AD-based DEM method recovered from instability upon mesh
refinement, we argue that instability would return for a larger load, and
using a fine mesh should _not_ be required for a method to remain stable. The
superior stability of the SF-based approach renders it suitable in simulations
involving large deformations, which is typical for rubber-like hyperelastic
materials. The examples in this work demonstrated the ability to apply the
full magnitude of an external load in a single load step while maintaining
high solution accuracy, something that FEM struggles to do. This feature
renders the method computationally efficient for very nonlinear materials,
where FEM typically requires breaking down the load into many smaller load
steps to achieve convergence.

The ability of the SF-based DEM model to solve hyperelastic material problems
in a single load step renders it very attractive in the topology optimization
of nonlinear hyperelastic materials. For classical FEM-based methods, multiple
Newton-Raphson iterations are needed in each design iteration to achieve
global force convergence, and the process is repeated for many topology
optimization iterations. Using the DEM framework to substitute nonlinear FE
simulations similar to our previous work [34] might speed up the overall
solution process and will be our future work. In addition, we also plan to
extend the SF-based GCN-DEM network to unstructured tetrahedral meshes to
perform simulations in irregular domains to harness the full power of the
graph convolution network.

## 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.

## CRediT author contributions

Junyan He: Conceptualization, Methodology, Software, Formal analysis,
Investigation, Data Curation, Writing - Original Draft. Diab Abueidda:
Conceptualization, Supervision, Writing - Review & Editing. Seid Koric:
Supervision, Writing - Review & Editing. Iwona Jasiuk: Supervision, Resources,
Writing - Review & Editing, Funding Acquisition.

## References

  * Nguyen-Thanh et al. [2020] Vien Minh Nguyen-Thanh, Xiaoying Zhuang, and Timon Rabczuk.  A deep energy method for finite deformation hyperelasticity.  _European Journal of Mechanics-A/Solids_ , 80:103874, 2020\.
  * Abueidda et al. [2022a] Diab W Abueidda, Seid Koric, Rashid Abu Al-Rub, Corey M Parrott, Kai A James, and Nahil A Sobh.  A deep learning energy method for hyperelasticity and viscoelasticity.  _European Journal of Mechanics-A/Solids_ , 95:104639, 2022a.
  * Nguyen-Thanh et al. [2021] Vien Minh Nguyen-Thanh, Cosmin Anitescu, Naif Alajlan, Timon Rabczuk, and Xiaoying Zhuang.  Parametric deep energy approach for elasticity accounting for strain gradient effects.  _Computer Methods in Applied Mechanics and Engineering_ , 386:114096, 2021.
  * Samaniego et al. [2020] Esteban Samaniego, Cosmin Anitescu, Somdatta Goswami, Vien Minh Nguyen-Thanh, Hongwei Guo, Khader Hamdia, X Zhuang, and T Rabczuk.  An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications.  _Computer Methods in Applied Mechanics and Engineering_ , 362:112790, 2020.
  * Noriega [2005] Leonardo Noriega.  Multilayer perceptron tutorial.  _School of Computing. Staffordshire University_ , 2005.
  * Hastie et al. [2009] Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman.  _The elements of statistical learning: data mining, inference, and prediction_ , volume 2.  Springer, 2009.
  * Chadha et al. [2022] Charul Chadha, Diab Abueidda, Seid Koric, Erman Guleryuz, and Iwona Jasiuk.  Optimizing hyperparameters and architecture of deep energy method.  2022\.
  * Zhuang et al. [2021] Xiaoying Zhuang, Hongwei Guo, Naif Alajlan, Hehua Zhu, and Timon Rabczuk.  Deep autoencoder based energy method for the bending, vibration, and buckling analysis of kirchhoff plates with transfer learning.  _European Journal of Mechanics-A/Solids_ , 87:104225, 2021\.
  * Vlassis et al. [2020] Nikolaos N Vlassis, Ran Ma, and WaiChing Sun.  Geometric deep learning for computational mechanics part i: anisotropic hyperelasticity.  _Computer Methods in Applied Mechanics and Engineering_ , 371:113299, 2020.
  * Frankel et al. [2022] Ari L Frankel, Cosmin Safta, Coleman Alleman, and Reese Jones.  Mesh-based graph convolutional neural networks for modeling materials with microstructure.  _Journal of Machine Learning for Modeling and Computing_ , 3(1), 2022.
  * Gao et al. [2022] Han Gao, Matthew J Zahr, and Jian-Xun Wang.  Physics-informed graph neural galerkin networks: A unified framework for solving pde-governed forward and inverse problems.  _Computer Methods in Applied Mechanics and Engineering_ , 390:114502, 2022.
  * Sanchez-Gonzalez et al. [2020] Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia.  Learning to simulate complex physics with graph networks.  In _International Conference on Machine Learning_ , pages 8459–8468. PMLR, 2020.
  * Hernández et al. [2022] Quercus Hernández, Alberto Badías, Francisco Chinesta, and Elías Cueto.  Thermodynamics-informed graph neural networks.  _arXiv preprint arXiv:2203.01874_ , 2022.
  * Ogoke et al. [2021] Francis Ogoke, Kazem Meidani, Amirreza Hashemi, and Amir Barati Farimani.  Graph convolutional networks applied to unstructured flow field data.  _Machine Learning: Science and Technology_ , 2(4):045020, 2021.
  * Chen et al. [2021] Junfeng Chen, Elie Hachem, and Jonathan Viquerat.  Graph neural networks for laminar flow prediction around random 2d shapes.  _arXiv preprint arXiv:2107.11529_ , 2021.
  * He et al. [2022a] Xiaodong He, Yinan Wang, and Juan Li.  Flow completion network: Inferring the fluid dynamics from incomplete flow information using graph neural networks.  _arXiv preprint arXiv:2205.04739_ , 2022a.
  * Schwarzer et al. [2019] Max Schwarzer, Bryce Rogan, Yadong Ruan, Zhengming Song, Diana Y Lee, Allon G Percus, Viet T Chau, Bryan A Moore, Esteban Rougier, Hari S Viswanathan, et al.  Learning to fail: Predicting fracture evolution in brittle material models using recurrent graph convolutional neural networks.  _Computational Materials Science_ , 162:322–332, 2019.
  * Raissi [2018] Maziar Raissi.  Deep hidden physics models: Deep learning of nonlinear partial differential equations.  _The Journal of Machine Learning Research_ , 19(1):932–955, 2018.
  * Abueidda et al. [2021] Diab W Abueidda, Qiyue Lu, and Seid Koric.  Meshless physics-informed deep learning method for three-dimensional solid mechanics.  _International Journal for Numerical Methods in Engineering_ , 122(23):7182–7201, 2021.
  * Guo et al. [2021] Hongwei Guo, Xiaoying Zhuang, and Timon Rabczuk.  A deep collocation method for the bending analysis of kirchhoff plate.  _arXiv preprint arXiv:2102.02617_ , 2021.
  * Haghighat et al. [2021] Ehsan Haghighat, Maziar Raissi, Adrian Moure, Hector Gomez, and Ruben Juanes.  A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics.  _Computer Methods in Applied Mechanics and Engineering_ , 379:113741, 2021.
  * Yu et al. [2018] Bing Yu et al.  The deep ritz method: a deep learning-based numerical algorithm for solving variational problems.  _Communications in Mathematics and Statistics_ , 6(1):1–12, 2018.
  * Liao and Ming [2019] Yulei Liao and Pingbing Ming.  Deep nitsche method: Deep ritz method with essential boundary conditions.  _arXiv preprint arXiv:1912.01309_ , 2019.
  * 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.
  * Abueidda et al. [2022b] Diab W Abueidda, Seid Koric, Erman Guleryuz, and Nahil A Sobh.  Enhanced physics-informed neural networks for hyperelasticity.  _arXiv preprint arXiv:2205.14148_ , 2022b.
  * Rezaei et al. [2022] Shahed Rezaei, Ali Harandi, Ahmad Moeineddin, Bai-Xiang Xu, and Stefanie Reese.  A mixed formulation for physics-informed neural networks as a potential solver for engineering problems in heterogeneous domains: comparison with finite element method.  _arXiv preprint arXiv:2206.13103_ , 2022.
  * Paszke et al. [2017] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  2017\.
  * Van Merriënboer et al. [2018] Bart Van Merriënboer, Olivier Breuleux, Arnaud Bergeron, and Pascal Lamblin.  Automatic differentiation in ml: Where we are and where we should be going.  _Advances in neural information processing systems_ , 31, 2018.
  * Ketkar and Moolayil [2021] Nikhil Ketkar and Jojo Moolayil.  Automatic differentiation in deep learning.  In _Deep Learning with Python_ , pages 133–145. Springer, 2021.
  * Zhu et al. [2019] Yinhao Zhu, Nicholas Zabaras, Phaedon-Stelios Koutsourelakis, and Paris Perdikaris.  Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data.  _Journal of Computational Physics_ , 394:56–81, 2019.
  * Geneva and Zabaras [2020] Nicholas Geneva and Nicholas Zabaras.  Modeling the dynamics of pde systems with physics-constrained deep auto-regressive networks.  _Journal of Computational Physics_ , 403:109056, 2020.
  * Zhang et al. [2020] Ruiyang Zhang, Yang Liu, and Hao Sun.  Physics-guided convolutional neural network (phycnn) for data-driven seismic response modeling.  _Engineering Structures_ , 215:110704, 2020.
  * Wandel et al. [2021] Nils Wandel, Michael Weinmann, and Reinhard Klein.  Teaching the incompressible navier–stokes equations to fast neural surrogate models in three dimensions.  _Physics of Fluids_ , 33(4):047117, 2021.
  * He et al. [2022b] Junyan He, Shashank Kushwaha, Charul Chadha, Seid Koric, Diab Abueidda, and Iwona Jasiuk.  Deep energy method in topology optimization applications, 2022b.  URL https://arxiv.org/abs/2207.03072.
  * Yao et al. [2020] Houpu Yao, Yi Gao, and Yongming Liu.  Fea-net: A physics-guided data-driven model for efficient mechanical response prediction.  _Computer Methods in Applied Mechanics and Engineering_ , 363:112892, 2020.
  * Logan [2016] Daryl L Logan.  _A first course in the finite element method_.  Cengage Learning, 2016.
  * Fish and Belytschko [2007] Jacob Fish and Ted Belytschko.  _A first course in finite elements_ , volume 1.  John Wiley & Sons New York, 2007.
  * Fey and Lenssen [2019] Matthias Fey and Jan E. Lenssen.  Fast graph representation learning with PyTorch Geometric.  In _ICLR Workshop on Representation Learning on Graphs and Manifolds_ , 2019.
  * Defferrard et al. [2016] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst.  Convolutional neural networks on graphs with fast localized spectral filtering.  _Advances in neural information processing systems_ , 29, 2016.
  * Pattanayak et al. [2017] Santanu Pattanayak, John S Pattanayak, and Suresh John.  _Pro deep learning with tensorflow_.  Springer, 2017.
  * Zhu et al. [1997] Ciyou Zhu, Richard H Byrd, Peihuang Lu, and Jorge Nocedal.  Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization.  _ACM Transactions on mathematical software (TOMS)_ , 23(4):550–560, 1997.
  * SIMULIA [2020] SIMULIA.  Abaqus, 2020.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * 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.