# A Deep Learning Approach to Nonconvex Energy Minimization for Martensitic
Phase Transitions

Xiaoli Chen<EMAIL_ADDRESS>Phoebus Rosakis<EMAIL_ADDRESS>Zhizhang Wu
<EMAIL_ADDRESS>Zhiwen Zhang<EMAIL_ADDRESS>Institute for Functional Intelligent
Materials & Department of Mathematics, National University of Singapore,
119077, Singapore  Department of Mathematics and Applied Mathematics,
University of Crete, Greece  Institute of Applied and Computational
Mathematics, Foundation for Research and Technology-Hellas, Heraklion, Greece
Department of Mathematics, The University of Hong Kong, Pokfulam Road, Hong
Kong SAR, China.

###### Abstract

We propose a mesh-free method to solve nonconvex energy minimization problems
for martensitic phase transitions and twinning in crystals, using the deep
learning approach. These problems pose multiple challenges to both analysis
and computation, as they involve multiwell gradient energies with large
numbers of local minima, each involving a topologically complex microstructure
of free boundaries with gradient jumps. We use the Deep Ritz method, whereby
candidates for minimizers are represented by parameter-dependent deep neural
networks, and the energy is minimized with respect to network parameters. The
new essential ingredient is a novel activation function proposed here, which
is a smoothened rectified linear unit we call SmReLU; this captures the
structure of minimizers where usual activation functions fail. The method is
mesh-free and thus can approximate free boundaries essential to this problem
without any special treatment, and is extremely simple to implement. We show
the results of many numerical computations demonstrating the success of our
method.

AMS subject classification: 35R05, 65N30, 68T99, 74A50, 74G65.

###### keywords:

Deep learning method; variational problems; mesh-free method; phase transition
problems; nonconvex energy minimization;

##  1 Introduction

Physics-informed deep learning methods [22, 19, 34, 13, 14, 23, 37, 38, 40]
have recently achieved great success in solving PDE problems arising in
different fields of STEM. This entails using deep neural networks to
approximate the solution of PDEs, often in many dimensions, leading to a mesh-
free numerical method in place of more established schemes, such as finite-
element or finite-difference methods.

In this approach [26, 38], the loss function to be minimized is the mean
square of the PDE residual. This method is rigorously known to converge for
linear elliptic and parabolic problems [40], and is successful with a number
of nonlinear PDEs [7, 32]. It is not clear how it would perform in situations
where there are weak solutions of the PDE that suffer from reduced smoothness.
Examples are hyperbolic conservation laws, as in nonlinear elastodynamics, and
solid-solid phase transitions described by nonconvex gradient energy functions
in nonlinear elastostatics [15, 1, 3]. In the latter case, the associated
Euler-Lagrange equations change type and lose ellipticity at some values of
the solution gradient [1]. In these situations, there are weak solutions with
discontinuous gradient across interfaces which are unknown a priori but are
part of the solution.

Another type of interface problem involves heterogeneous media whose
properties have a priori specified high-contrast or discontinuous features,
such as PDEs whose coefficients are discontinuous across a fixed interface,
where solutions suffer a loss of smoothness. A number of sophisticated finite
element and finite difference methods for interface problems have been
developed [2, 30, 17, 8, 9, 36, 28]. These methods are quite accurate, but
they often require special treatment of the interface and jump conditions
across it in creating the finite element mesh and/or associated basis
functions for elements that intersect it.

Recently, deep learning was used to solve interface problems involving linear
elliptic PDE systems with discontinuous coefficients and/or discontinuous
forcing [42], motivated from biomechanics [43]. The problems studied in [42]
are solved by energy minimization. The energy functional is strictly convex,
and a weak solution of the PDE (Euler-Lagrange equation) is its unique
minimizer. The interface is fixed at the location of the jump in the
coefficients of the PDE. The Deep Ritz method [14] used in [42] approximates
solutions by a deep neural network (DNN) and minimizes the energy as a
function of the DNN parameters (weights and biases). The DNN method is mesh-
free and does not require any special structure to capture the location of the
interface inside the domain spontaneously and accurately.

In phase transition problems, interfaces separating parts of the body in
different phases (solid and liquid, or austenite and martensite) are free
boundaries whose position or shape is not known a priori. Rather, they are
part of the solution and can evolve in time-dependent problems. Phase
boundaries (domain walls, twin boundaries, etc. in different physical
settings) are characterized by a jump discontinuity or narrow transition layer
of the unknown scalar or vector field (e.g., order parameter, deformation,
polarization) or its spatial derivatives. To identify the unknown phase
boundary, special augmented models have been developed, notably the phase
field and level set methods. The prototypical phase field models are described
by the Allen-Cahn and Cahn-Hilliard equations. An essential ingredient of most
continuum models of phase transitions is a nonconvex, multiwell energy density
function. The DNN approach has been used to solve free boundary problems of
the Stefan type [41] and to describe phase boundaries in the Allen-Cahn and
Cahn-Hilliard equations [33, 37]. The method uses standard DNNs without any
modification to identify phase boundaries.

Martensitic phase transitions [15, 3, 1] are often modeled using a free energy
density that is a nonconvex, multi-well potential depending on the gradient of
the unknown deformation field. In fact, the martensitic phase itself involves
two or more symmetry-related variants, known as twins, that correspond to
different energy wells.

(a)

(b)

(c)

Figure 1: Twinned Microstructures: (a) Experimental in CuAlNi single crystal
[10](reproduced here with permission of the authors). (b), (c): Comparison of
results of different numerical methods for a simplified 2D model: (b) Finite
Difference solution [20] (reproduced here with permission of the authors) and
(c) our DNN solution of the 2D regularized problem (10).

The fact that the energy depends on the gradient of the unknown function
(deformation) introduces special challenges. Local energy minimizers involve
phase boundaries separating parts of the body in different phases. Across
these interfaces, the deformation field remains continuous, but its gradient
jumps. Continuity of the deformation field imposes restrictions on the jump of
its gradient (Hadamard compatibility conditions [1, 3]). In many crystalline
alloys, the austenite and martensite energy wells violate these conditions,
which are nonetheless satisfied between martensitic twin variants. As a
result, there is no absolute energy minimum achievable [3]. Instead, the
energy tends to its infimum in the limit, along a sequence of deformations
with a laminated microstructure of parallel twin interfaces. Such minimizing
sequences involve a larger and larger number of twin boundaries, and converge
to a limit state compatible with austenite [3]. Such laminates are observed in
experiments, Fig. 1(a). This microstructure is called finely twinned
martensite; it can also occur because of incompatibility between twins of
different orientations and/or rigid (null Dirichlet) boundary conditions, Fig.
1.

Near an incompatible boundary, twin layers taper into needles, and often an
individual layer is observed splitting into two or more needles in
experiments, Fig. 1(a) [10, 5]. This is known as twin branching and has been
studied analytically [25, 11, 20, 39] and numerically [29, 21, 20, 12] as a
mechanism for attainment of an energy minimum in the presence of interfacial
energy. In order to capture this complex behaviour, it is possible to consider
some lower-dimensional problems. Many of these have an energy of the form

$\int_{\Omega}W(\nabla u(\bm{x}))d\bm{x},$ (1)

where $\Omega\subset\mathbb{R}^{n}$, $u:\Omega\to\mathbb{R}^{m}$ and $n=2$ or
$3$, while $m=n$ or $m=1$ for the vector and scalar case, respectively. In all
cases, $W:\mathbb{R}^{m\times n}\to\mathbb{R}$ is a nonconvex function with
multiple wells. Letting the derivative of $W$ be $S:\mathbb{R}^{m\times
n}\to\mathbb{R}^{m\times n}$, the Euler-Lagrange equation of (1) are

$\hbox{div}\,S(\nabla u)=0\quad\hbox{in $\Omega$.}$ (2)

Because of loss of rank one convexity [3, 1] for values of $\nabla u$ between
wells, this PDE system loses ellipticity. Because of ellipticity loss, (2) has
weak solutions with phase boundaries, across which $\nabla u$ jumps
discontinuously from phase to phase [1] and second derivatives fail to exist.
This precludes using the physics-informed neural network method (PINN), which
requires the strong form of (2) to be valid [37]. An added problem is the
extensive loss of uniqueness of weak solutions of (2) [15], many of which are
energetically unstable (not even local energy minima). The alternative is to
minimize the energy (1) [3, 25, 12, 39, 11]. We use the Deep Ritz method [13],
namely energy minimization in terms of DNNs. An additional advantage is that
piecewise smooth local minimizers of (1) are also weak solutions of (2) which
are at least metastable, in contrast to some solutions of (2).

Here we demonstrate that the Deep Ritz method spontaneously captures complex
microstructures, such as the finely twinned one, including gradient jumps
across multiple curved interfaces and splitting/tapering topological
transitions near boundaries [20, 25, 11, 12, 39, 21, 5]. Notably, the method
is mesh-free; this allows it to accurately describe curved interfaces of
arbitrary orientation, in contrast to FEM and other mesh dependent methods
[29]. No special structure of the DNN is required for this, except a specially
designed activation function whose form is suggested by the structure of
solutions to simple problems.

The outline of this paper is as follows: In Section 2, we introduce the Deep
Ritz method for energy minimization in the setting of nonlinear elasticity,
and apply it to various problems in one and two dimensions, involving
nonconvex potentials of the deformation gradient. A crucial role is played by
the activation function. Our first choice is guided by the observation that in
the simplest one-dimensional problem, an exact solution is provided by the
piecewise linear ReLU activation function. We then consider a regularization
of the problem, by including the second derivative of the deformation in the
energy [6]. Minimizers are smooth, with transition layers replacing derivative
jumps. This motivates our introduction of a smoothened version of the ReLU
activation function, which we call SmReLU. In two-dimensional problems, use of
SmReLU captures spatially complex microstructures with fine twinning and twin
branching, Fig. 1(c), whereas traditional activation functions fail. Section 3
reports on extensive numerical computations. A discussion of our results,
their significance, and comparison to previous work is in Section 4.

##  2 Methods

###  2.1 The Minimization Problem

The problem we study in this paper concerns the attempt to minimize a
nonconvex functional that is a low-dimensional model for an austenite-
martensite phase transition. The problem involves finding

$\min E\\{u\\},\quad E\\{u\\}=\int_{\Omega}W(\nabla{u}(\bm{x}))d\bm{x}$ (3)

over suitable functions $u:{\Omega}\to\mathbb{R}$, where the function
$W:\mathbb{R}^{2}\to\mathbb{R}$ is a two-well potential. Here we choose [25,
11]

$W(\nabla u)=\frac{1}{2}[u_{x}^{2}(u_{x}-1)^{2}+u_{y}^{2}],$

with $\nabla u=(u_{x},u_{y})$. The wells (minima) of $W$ are at $(0,0)$ and
$(1,0)$. The domain ${\Omega}=[0,L]\times[0,1]\subset\mathbb{R}^{2}$. Here $u$
is subject to the Dirichlet boundary conditions

$u(x,y)=\gamma x\quad\forall\;(x,y)\in\partial\Omega.$ (4)

Of particular interest is the case $\gamma=1/2$ as the linear function
$u(x,y)=x/2$ for all $(x,y)\in{\Omega}$ satisfying the boundary conditions has
$\nabla u=(1/2,0)$ which is a saddle point of $W$. The zero-energy minima
$u\equiv 0$ and $u=x$ on ${\Omega}$ are incompatible with the boundary
conditions, whereas there are sequences of functions approaching zero energy,
namely problem (3) does not have a minimizer but only minimizing sequences,
namely, one can construct sequences of functions $u_{n}$ with $E\\{u_{n}\\}\to
0$ as $n\to\infty$. This occurs due to geometric incompatibility of the
energy-minimal states $u=\text{const}$ and $u=x+\text{const}$ on ${\Omega}$
with the boundary conditions (4).

On the other hand, the mixed boundary conditions

$u(0,y)=0,\quad u(1,y)=\gamma,\quad 0<y<1,$ (5)

at the vertical sides $x=0$, $L$ with the horizontal sides $y=0$, $1$ free,
are compatible with the ansatz

$u(x,y)=u(x),\quad\;0<y<1,$

(using $u$ for both functions by notation abuse) which reduces (3) to the
following one-dimensional problem. Let $\Omega=[0,1]$, define:

$W(z)=z^{2}(1-z)^{2},\quad z\in\mathbb{R},$

so that $W$ is a double well potential with minima at 0,1. Suppose
$u:[0,1]\to\mathbb{R}$ satisfies $u(0)=0$ and $u(1)={\gamma}$ for a given
constant ${\gamma}\in\mathbb{R}$. Then, minimize

$\min\int_{0}^{1}W(u^{\prime}(x))dx,$ (6)

with $u^{\prime}$ the derivative of $u$. Letting

$f(x)={\frac{x+|x|}{2}},\quad x\in\mathbb{R},$ (7)

a solution of the one-dimensional problem (6) (with zero energy) is

$u(x)=f\left(x-1+{\gamma}\right).$ (8)

This is piecewise smooth with a derivative discontinuity at $x=1-{\gamma}$.

######  Remark 2.1.

Any partition of $[0,1]$ into intervals where the slope of $u$ alternates
between 0 and 1 (while $u$ remains continuous) is a minimizer (with zero
energy). This shows the massive loss of uniqueness of minimizers of the one
dimension (1D) problem (6), but also the two dimension (2D) problem with
boundary conditions (5).

######  Remark 2.2.

In order to make problem (6) well posed, it is common [6] to add a higher
gradient regularization, also known as capillarity, to the energy which
becomes

$\min E_{\varepsilon}\\{u\\},\quad
E_{\varepsilon}\\{u\\}=\int_{0}^{1}\left[W(u^{\prime}(x))dx+{\frac{\varepsilon^{2}}{2}}[u^{\prime\prime}(x)]^{2}\right]dx,$
(9)

where $\varepsilon>0$ is the higher-gradient coefficient, a small parameter.
It is known [6] that (9) has an essentially unique solution, where the
gradient discontinuity of (8) is replaced by a single, smooth transition layer
between values 0 and 1 of the derivative $u^{\prime}$. For small $\varepsilon$
the thickness of this layer is of order $\varepsilon$ and the gradient term
contributes an interfacial or surface energy to the interface, also of order
$\varepsilon$.

The analogous 2D regularized energy is

$E_{\varepsilon}\\{u\\}=\int_{\Omega}\left[W(\nabla{u}(x,y))+\frac{\varepsilon^{2}}{2}\Bigl{[}u_{xx}(x,y)\Bigr{]}^{2}\right]dxdy.$
(10)

###  2.2 Definition of deep neural networks (DNNs)

We briefly recall the definition and properties of a DNN. There are two
ingredients in defining a DNN. The first one is a (vector) linear function of
the form $T:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$, defined as $T(x)=Ax+b$,
where $A=(a_{ij})\in\mathbb{R}^{m\times n}$, $x\in\mathbb{R}^{n}$ and
$b\in\mathbb{R}^{m}$. The second one is a nonlinear activation function
$\sigma:\mathbb{R}\rightarrow\mathbb{R}$. A frequently used activation
function in the artificial neural network literature, the sigmoid function, is
defined as $\sigma(x)=(1+e^{-x})^{-1}$. Here, for reasons that will become
clear below, we will use another activation function, known as the rectified
linear unit (ReLU), defined by $\sigma(x)=\max(0,x)$ [27]. Observe that
$\sigma(x)=f(x)$ defined in (7), which provides a solution (8) of the one-
dimensional minimization problem (6). By applying the activation function to
each component $x_{j}$ ($j=1,\ldots n$) of the vector $x\in\mathbb{R}^{n}$,
one can define a vector activation function
$\sigma:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ with components
$\sigma(x_{j})$, $j=1,\ldots n$.

Equipped with those definitions, we are able to construct a continuous
function $F(x)$ as an alternating composition of linear transformations and
(vector) activation functions as follows:

$F(x)=T^{k}\circ\sigma\circ T^{k-1}\circ\sigma\cdot\cdot\cdot\circ
T^{1}\circ\sigma\circ T^{0}(x).$ (11)

Here $T^{i}(x)=W_{i}x+b_{i}$ is an affine transformation, with the weights
$W_{i}$ as of yet unspecified matrices and the biases $b_{i}$ unspecified
vectors. Also $\sigma(\cdot)$ is the component-wise defined activation
function. Dimensions of $W_{i}$ and $b_{i}$ are chosen to make (11)
meaningful. The parametrization (11) is called a $(k+1)$-layer DNN, which has
$k$ hidden layers. Denoting all the undetermined coefficients (namely $W_{i}$
and $b_{i}$) in (11) as $\theta\in\Theta$, where $\theta$ is a high-
dimensional vector and $\Theta$ is the space of $\theta$. The DNN
representation of a continuous function can be viewed as

$\displaystyle F=F(x;\theta).$ (12)

Let $\mathbb{F}=\\{F(\cdot,\theta)|\theta\in\Theta\\}$ denote the set of all
functions expressible by DNNs parametrized by $\theta\in\Theta$ as in (11).
Then $\mathbb{F}$ provides an efficient way to represent unknown continuous
functions [19], compared to a linear solution space used in classic numerical
methods, e.g., a trial space spaced by linear nodal basis functions in the
FEM.

######  Remark 2.3.

In the computation, we propose a Smoothened Rectified Linear Unit (SmReLU)
activation function to approximate $u$ for the regularized energy problems (9)
and (10). We refer readers to Section 3 for details.

###  2.3 The DNN representation of the nonconvex energy minimization

Given the energy functional $E\\{u\\}$ from (3), we can derive a DNN-based
numerical method to solve the minimization problem in order to compute the
solution $u$, in the spirit of [13]. Here we seek the solution of the energy
minimization problem

$u=\operatorname*{arg\,min}_{v\in\mathbb{H}^{1}(\Omega)}E\\{v\\},$ (13)

subject to the specified boundary condition $u(x)=g(x)$ on $\partial\Omega$.

From the perspective of scientific computing, the energy minimization problem
(13) can be solved using numerical methods, such as FDMs or FEMs. From the
perspective of machine learning however, the numerical solution of $u(x)$ is
interpreted as a function with $x\in\mathbb{R}^{2}$ as its input and
$u\in\mathbb{R}^{1}$ as its output, that can be approximated by a DNN $F(x)$
defined in (11). As a result, problem (13) is replaced by the problem

$\tilde{u}=\operatorname*{arg\,min}_{F\in\mathbb{F}_{b}}E\\{F\\},$ (14)

Here $\tilde{u}$ denotes the DNN solution of the problem of minimizing the
energy over the DNNs (11) in $\mathbb{F}$. Let

$J(\theta)=E\\{F(\cdot,\theta)\\}=\int_{\Omega}W(\nabla
F(x,\theta))dx,\quad\theta\in\Theta,$ (15)

where $F\in\mathbb{F}$ is as in (11) and (12). Then, the minimization problem
(14) reduces to the (finite-dimensional) optimization problem

$\tilde{\theta}=\operatorname*{arg\,min}_{\theta\in\Theta}J(\theta),\quad\tilde{u}(x)=F(x,\tilde{\theta}).$
(16)

This yields $\tilde{u}$, the DNN representation of the solution of (14).

Notice that the original energy minimization problem (13) is non-convex. The
variational problem (16) is non-convex as well. Clearly, the issue of local
minima and saddle points is nontrivial, which brings essential challenges to
many existing optimization methods. Since the parameter space $\Theta$ is
typically very large, one often uses the Stochastic Gradient Descent (SGD)
method [4] to solve (16).

How to impose boundary conditions in the DNN method is an important issue
[26]. For the Dirichlet boundary condition, we use the equivalent of the
Nitsche method in FEM for DNN, known as the Deep Nitsche Method [31]. We add a
boundary penalty term in the energy function to address this issue.
Specifically, one adds a soft constraint (a boundary integral term) to the
energy $J(\cdot)$ defined in (15) and obtains

$\tilde{u}_{\tau}=\operatorname*{arg\,min}_{\theta\in\Theta}\Big{[}J(\theta)+\tau\int_{\partial\Omega}\big{(}F(x,\theta)-g(x)\big{)}^{2}dx\Big{]},$
(17)

where $\tau>0$ is the penalty parameter to scale and balance the losses from
the energy part and the boundary part. The idea is that the term
$\int_{\partial\Omega}\big{(}F(x,\theta)-g(x)\big{)}^{2}dx$ will approach zero
when we increase the parameter $\tau$ in the calculation, whereas the
Dirichlet boundary condition $u=g$ on $\partial{\Omega}$ will be satisfied in
a weak $L^{2}$ sense, and one hopes that $\tilde{u}_{\tau}\to\tilde{u}$ in
$\mathbb{F}$ as $\tau\to\infty$.

It has been recently pointed out that the penalty parameter $\tau$ should not
be chosen arbitrarily but adaptively in a rational way stemming from the
specific problem [16]. Here we also point out that the $\tau$ can be given
physical meaning as the stiffness of the austenite phase that exists beyond
the boundary $y=1$, whereas the martensite phase occupies ${\Omega}$. The
exact satisfaction of boundary conditions such as (4) would correspond to
infinitely rigid austenite.

###  2.4 Implementation of the Deep Ritz DNN method

Since the number $\hbox{dim}\,\Theta$ of degrees of freedom in the
optimization problem (16) is quite large, we apply the SGD method on the
parameter space $\Theta$ to solve it. Notice that $\theta$ is a high-
dimensional vector and let $\theta_{k}$ be any component of $\theta$. We
approximate the corresponding component $\partial
J(\theta)/\partial\theta_{k}$ of the gradient of $J(\theta)$ as follows:

$\displaystyle\frac{\partial J(\theta)}{\partial\theta_{k}}$
$\displaystyle=\int_{\Omega}\frac{\partial\big{(}W(\nabla
F(x,\theta))\big{)}}{\partial\theta_{k}}dx+\tau\int_{\partial\Omega}\frac{\partial\big{(}F(x,\theta)-g(x)\big{)}^{2}}{\partial\theta_{k}}ds$
$\displaystyle\approx\frac{vol(\Omega)}{N}\sum_{i=i}^{N}\frac{\partial\big{(}W(\nabla
F(x_{i},\theta))\big{)}}{\partial\theta_{k}}+\frac{\tau
S(\partial\Omega)}{N_{b}}\frac{\partial\big{(}(F(y_{j},\theta)-g(y_{j}))^{2}\big{)}}{\partial\theta_{k}},$
(18)

where the collocation points $x_{i}\overset{i.i.d.}{\sim}Unif(\Omega)$ are
uniform random points that are sampled from the physical domain $\Omega$,
$y_{j}\overset{i.i.d.}{\sim}Unif(\partial\Omega)$ are uniform random points
that are sampled from the boundary $\partial\Omega$, $vol(\Omega)$ is the
volume of the domain, and $S(\partial\Omega)$ is the perimeter of the domain.
In the context of the deep learning method, $N$ and $N_{b}$ are called batch
numbers, which means the number of training samples utilized in one iteration.

After we compute the approximation of the gradient with respect to
$\theta_{k}$, we can update each component of $\theta$ as

$\theta_{k}^{n+1}=\theta_{k}^{n}-\eta\frac{\partial
J(F(\cdot,\theta))}{\partial\theta_{k}}|_{\theta_{k}=\theta_{k}^{n}},$ (19)

where $\eta$ is the learning rate. In our numerical computation, we use the
Adam optimizer [24], a modified version of the SGD method. Our results
indicate that the Adam optimizer is very effective in solving this non-convex
optimization problem.

######  Remark 2.4.

From the derivation of the DNN formulation, one can see that the proposed
method automatically deals with the phase boundary (or discontinuity in the
gradient of the solution) without knowing the locations of the discontinuity a
priori.

Figure 2: Schematic of a DNN.

Fig. 2 shows a possible network layout for approximating $u$ with input
$(x,y)$. For the output $u$, we can use automatic differentiation techniques
to compute derivatives of $u$ that appear in the energy in the minimization
problem.

##  3 Numerical Results

We use a fully connected DNN (Fig.2) to perform our computations. The weights
are initialized with truncated normal distributions with variances set as two
over the sum of input and output unit numbers. The biases are initialized as
zero. We use the Adam optimizer with a learning rate of $10^{-3}$. In the main
text, the training takes $100,000$ iterations for one-dimensional problems,
while it takes $300,000$ iterations for two-dimensional problems.

###  3.1 One-dimensional problem

We demonstrate how to use DNNs to minimize the energy in a one-dimensional
problem. We consider the following energy minimization problem:

$\min\int_{0}^{1}W(u^{\prime}(x))dx,$ (20)

where

$W(z)=z^{2}(1-z)^{2}\;\;\forall z\in\mathbb{R},\quad u(0)=0,\quad
u(1)=\gamma.$ (21)

The integrand $W$ has two wells at 0 and 1 and we choose typically choose
$\gamma=1/2$. The functions $u$ are approximated by DNNs (see Methods)

$u_{NN}(x)=F(x;\theta),$

where $x$ is the input of the DNN and $\theta$ is the vector of weights
$w_{i}$ and biases $b_{i}$ as in (11), which are the trainable parameters in
the DNN.

On one hand, the parameters $\theta$ of the DNN $u_{NN}(x)$ should minimize
the energy in (20) (Deep Ritz Method [13]). We construct the corresponding
loss function as follows:

$loss_{e}(\theta)=\frac{1}{N}\sum_{i=1}^{N}W\left(\frac{\partial F}{\partial
x}(x_{i},\theta)\right),$ (22)

where $\\{x_{i}\\}_{i=1}^{N}$ are the collocation points in $\Omega$. The sum
here approximates the integral in (20). The derivative
$u^{\prime}_{NN}=\partial F/\partial x$ is computed with automatic
differentiation.

On the other hand, $u_{NN}(x)$ should satisfy the boundary conditions, which
corresponds to minimizing the loss function

$loss_{b}(\theta)=(F(0,\theta)-0)^{2}+(F(1,\theta))-\gamma)^{2}.$ (23)

This corresponds to the Deep Nitsche method [31]. Combining the two loss
functions (22) and (23), the total loss function to be minimized is

$loss(\theta)=loss_{e}(\theta)+\tau loss_{b}(\theta),$ (24)

where $\tau$ is the penalty parameter to scale and balance the losses from the
energy part and the boundary part. During training, the shared parameters
$\theta$ are adjusted by back-propagating the error obtained by minimizing a
loss function (24) that is the weighted sum of the above two constraints.

####  3.1.1 Activation function

In this problem, an important role is played by the activation function. As is
observed in Methods, Section 2.1, an exact solution of (20) is provided by
(8). We observe that in (8), the function $f(x)$ defined by (7) is identical
with the well-known Rectified Linear Unit activation function

$\hbox{ReLU: }\;\sigma(x)=\max\\{0,x\\}={\frac{x+|x|}{2}},\quad
x\in\mathbb{R}.$ (25)

Thus by choosing one weight equal to 1 and one bias equal to ${\gamma}-1$, the
DNN can provide the exact solution (8) to this simple one-dimensional problem.
This strongly suggests that we choose ReLU as the activation function until
further notice. To test the suitability of our choice, we explore how the
activation function affects the learning results. The results are shown in
Fig. 3. It is clear that only ReLU captures an exact minimizer of the problem
(any piecewise smooth function with slope 0 or 1 almost everywhere that
satisfies the boundary conditions), whereas other choices including Tanh,
Sigmoid, and even Leaky ReLU perform poorly. It is remarkable here that the
Deep Ritz method with ReLU captures exact global minimizers in this simple
problem, due to the piecewise smooth character of this activation function.
The other activation functions struggle and are indeed unable to capture these
weak solutions with derivative jumps. Later on we will propose a modification
of ReLU to treat more complex problems.

Figure 3: Results of solving problem (20) with different activation function.
Shown is the plot of $u(x)$.

####  3.1.2 Effect of learning rate for 1D problem

We investigate how the learning rate affects the results for 3 values of the
parameter $\gamma=0.25,0.5,0.75$. We use the ReLU activation function. The
results are shown in Fig. 4. We can see if we use very large learning rate
$10^{-1}$ to train the loss function, the solution will go to a constant
$u(x)=\frac{\gamma}{2}$, which is a local minimum of the elastic energy, but
does not satisfy the boundary conditions. If we use very small learning rate,
the solution will get stuck in another local minimum $u(x)=\gamma x$ for the
case $\gamma=0.25$, which however is not a global minimum, although it
satisfies the boundary conditions. For intermediate values of the learning
rate, a global minimum is achieved. We therefore should exercise some care in
choosing the learning rate carefully.

(a)

(b)

(c)

Figure 4:  1d results with ReLU activation function and different learning
rates with $iteration=100,000$, NN: $3\times 128$. (a) $\gamma=0.25$; (b)
$\gamma=0.5$; (c) $\gamma=0.75$.

####  3.1.3 Effect of DNN size for 1D problem

Then we explore how the depth and width affect the results when
$\gamma=0.25,0.5,0.75$. We use the ReLU activation function. The results are
shown in Fig. 5. For small $\gamma$, we need to use a larger DNN to
approximate the solution $u$.

(a)

(b)

(c)

Figure 5:  1d results with different size DNN. Here $iteration=100,000$ and
learning rate $\eta=10^{-2}$. (a) $\gamma=0.25$; (b) $\gamma=0.5$; (c)
$\gamma=0.75$.

####  3.1.4 Regularized 1D problem

In order to cure the massive loss of uniqueness of minimizers for problem (20)
it is common to add higher gradients to the energy; see Remarks 2.1 and 2.2,
also [6]. Thus we consider the minimization problem

$\min E_{\varepsilon}\\{u\\},\quad
E_{\varepsilon}\\{u\\}=\int_{0}^{1}\left[W(u^{\prime}(x))dx+{\frac{\varepsilon^{2}}{2}}[u^{\prime\prime}(x)]^{2}\right]dx,$
(26)

where $\varepsilon>0$ is a small parameter, with boundary conditions $u(0)=0$
and $u(1)=\gamma$.

It is well known that the higher gradients in the energy (26) smoothen
derivative discontinuities, and replace them with a smooth transition whose
width is of order $\varepsilon$, e.g. [6]. Thus we expect that the DNN with a
ReLU activation function will not be able to capture this accurately due to
its lack of smoothness. This leads us to construct a new activation function
by smoothening ReLU. In (25), replace the absolute value $|x|=\sqrt{x^{2}}$ by
$\sqrt{x^{2}+\rho^{2}}$, where $\rho$ is a small parameter. We call the
resulting function a Smoothened Rectified Linear Unit or

$\hbox{SmReLU: }\;\sigma(x)=\frac{x+\sqrt{x^{2}+\rho^{2}}}{2},\quad
x\in\mathbb{R}.$ (27)

SmReLU is a smooth function for $\rho>0$; it reduces to ReLU for $\rho=0$. We
typically use $\rho=0.1$ in computations.

Next, we consider how the initial condition, the value of the higher gradient
coefficient $\varepsilon$, and the activation function (ReLU vs SmReLU) affect
our computations.

The results are shown in Fig. 6. We use a $3\times 128$ DNN with ReLU and
SmReLU activation function to approximate the solution $u$. The learning rate
is $\eta=10^{-2}$ and $\tau=500$.

If we use the ReLU activation function, the optimal $u$ is not differentiable
at some points, and moreover, more than one derivative jumps (kinks) occur
when the initial condition has oscillations; see Fig. 6. On the other hand,
problem (26) has a unique smooth minimizer (modulo reflection about $x=1/2$)
[6] with a single smoothened kink (transition layer between slopes 0 and 1).
This is well captured by the SmReLU DNN, which gives solutions with one smooth
kink. The only exception which has a local minimum with two kinks occurs in
Fig. 6 (a4). For larger $\varepsilon$, the DNN learns a more smooth solution
with a wider transition layer, in agreement with theory [6]. If $\varepsilon$
is large enough, the regularization term dominates the loss function. So to
minimize the second derivative energy term, $u$ is very close to a linear
function for this case. We conclude that the SmReLU activation function (27)
we proposed works well for this problem and captures the behavior expected
from theory.

(a1)

(a2)

(a3)

(a4)

(b1)

(b2)

(b3)

(b4)

(c1)

(c2)

(c3)

(c4)

Figure 6:  1d regularized problem results. Here $iteration=100,000$,
$\eta=10^{-2}$ and NN: $3\times 128$. (a) $\gamma=0.25$; (b) $\gamma=0.5$; (c)
$\gamma=0.75$. First column: Random initial condition and ReLU activation
function; Second column: $\gamma x+0.1\sin(4x)$ initial condition and ReLU
activation function; Third column: Random initial condition and SmReLU
activation function; Last column: $\gamma x+0.1\sin(4x)$ initial condition and
SmReLU activation function.

###  3.2 Two-dimensional problem

Our main goal is to consider the following two dimensional nonconvex problem:

$\min\int_{\Omega}W(\nabla u(x,y))dxdy,$ (28)

and its regularized counterpart (33). Here

$W(\nabla u)=\frac{1}{2}[u_{x}^{2}(1-u_{x})^{2}+u_{y}^{2}],$ (29)

and $\Omega=[0,L]\times[0,1]$ with $L=1$ or $2$. We construct a fully
connected DNN for $u$ with input $(x,y)$. The loss function is defined as
follows:

$loss=loss_{e}+\tau loss_{b},$ (30)

where

$\displaystyle loss_{e}$ $\displaystyle=\frac{1}{N}\sum_{i=1}^{N}W(\nabla
u_{NN}(x_{i},y_{i})),$ $\displaystyle loss_{b}$
$\displaystyle=\frac{1}{N_{b}}\sum_{j=1}^{N_{b}}(u_{NN}(x_{j}^{b},y_{j}^{b})-u_{true}(x_{j}^{b},y_{j}^{b}))^{2},\
$

$\\{(x_{i},y_{i})\\}_{i=1}^{N}$ are the collocation points in the domain
$\Omega$, $\\{(x_{j}^{b},y_{j}^{b})\\}_{i=1}^{N_{b}}$ are the collocation
points at the parts of the boundary where Dirichlet boundary conditions are
specified, while $\tau$ is the weight to scale and balance the losses from the
energy ($loss_{e}$) and the boundary conditions ($loss_{b}$).

####  3.2.1 Mixed Boundary Conditions

Consider the mixed boundary conditions:

$u(0,y)=0,~{}u(1,y)=\gamma,~{}~{}~{}0<\gamma<1,$ (31)

and the other two sides are free. Functions $u$ that are independent of $y$
minimize the energy provided that they are solutions of the one-dimensional
problem (20). We use both ReLU and SmReLU activation function, and $3\times
128$ DNNs to approximate $u$ when $\gamma=0.25,~{}0.5,~{}0.75$. The learned
solutions are shown in Fig. 7. In the top row of Fig. 7, we show results from
the ReLU DNN. The bottom row of Fig. 7 shows the SmReLU DNN minimizers. In the
latter case solutions only have one gradient jump, whereas ReLU solutions
typically have two.

Our numerical results match the exact $y$-independent solutions of the one-
dimensional problem.

(a)

(b)

(c)

Figure 7:  2d mixed boundary conditions. (a) $\gamma=0.25$; (b) $\gamma=0.5$;
(c) $\gamma=0.75$. Top: ReLU; bottom: SmReLU activation function.

####  3.2.2 Dirichlet Boundary Conditions

In the following, we will explore Dirichlet boundary condition

$u(x,y)=\gamma x,\quad\forall(x,y)\in\partial\Omega.$ (32)

A typical value of interest is ${\gamma}=1/2$. In order to get a zero energy
deformation $u_{x}$ should take values 0 or 1 with $u_{y}=0$, which is what
happens in Fig. 7, the mixed boundary condition case considered in subsection
3.2.1. On the other hand, here the boundary conditions (32) dictate that
$u_{x}$ be close to 1/2 near the boundary so there is an incompatibility
between the energy wells and the boundary conditions. In fact it is possible
to show that problem (28) subject to (32) does not have a minimizer [3, 25].
Instead, the energy can be effectively minimized in the limit by a sequence of
deformations that involve a laminated microstructure, consisting of vertical
parallel interfaces between domains with gradients in alternating wells
$(u_{x},u_{y})=(0,0)$ and $(1,0)$. Such a minimizing sequence involves a
larger and larger number of twin boundaries (jumps in $\nabla u$ along lines
$x=$const.), in effect converging to a mixture of phases, where in the limit
$u_{x}$ converges weakly to the average value $1/2$, so that it becomes
compatible with the boundary conditions $u=x/2$ on $\partial{\Omega}$.

We show results using the SmReLU activation function and $\gamma=0.5$ in Fig.
8 and 9. Laminated microstructures prevail with alternating bands having
$u_{x}$ jumping between values near 0 (blue) and 1 (yellow) to minimize the
energy density $W$, except near horizontal boundaries where these values are
incompatible with boundary conditions (32) that require $u_{x}=1/2$. This
causes each band to split into two or more near the boundary as is also
observed in experiments Fig. 1(a) [10]. In Fig. 8, we study the effect of DNN
depth (number of layers) in capturing local minima with more and more bands.
With one layer only, the solution is stuck at an unstable state (saddle point
of $W$) with $u_{x}\approx 0.5$, Fig. 8(a). The number of bands increases as
we increase the number of layers (depth) of the DNN, albeit with some
oscillations. The energy also decreases and then appears to tend to a
constant. In Fig. 9, we fix the depth as $5$ and increase the width of the
DNNs. The energy will decrease for wider DNNs. If we use larger DNNs (deeper
or wider), we can get more refined microstructures in our results, but we
recall that there is no limit to the number and fineness of bands in this
problem all the way to infinity [3, 25], as there is no global minimum, but
apparently a sequence of local minima with an increasing number of bands, like
the ones observed here, that approach the infimum of energy. In a wider
network, there will be more parameters (weights and biases) to be trained,
rendering optimization more costly.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 8:  2d Dirichlet boundary condition with SmReLU activation function.
Here $\gamma=0.5$, $iteration=300,000$, $\eta=10^{-3}$.(a) NN: $1\times 128$;
(b) NN: $3\times 128$; (c) NN: $5\times 128$; (d) NN: $7\times 128$; (e) NN:
$9\times 128$; (f) NN: $11\times 128$; (g) NN: $13\times 128$; (h) NN:
$15\times 128$.

(a)

(b)

(c)

(d)

Figure 9:  2d Dirichlet boundary condition with SmReLU activation function.
Here $\gamma=0.5$, $iteration=300,000$, $\eta=10^{-3}$. Top: the energy with
different (a) NN: $5\times 16$; (b) NN: $5\times 32$; (c) NN: $5\times 64$;
(d) NN: $5\times 256$.

####  3.2.3 Regularized Scalar Problem with Dirichlet Boundary Conditions

In order to study a problem that actually has a global minimum we use the
common approach of adding higher gradients and studying the regularized
problem

$\min\int_{\Omega}\bigg{[}W(\nabla
u(x,y))+\frac{\varepsilon^{2}}{2}\big{[}\frac{\partial^{2}u(x,y)}{\partial
x^{2}}\big{]}^{2}\bigg{]}dxdy,$ (33)

where $W$ is given by (29), $\Omega=[0,2]\times[0,1]$ and the boundary
condition is the linear Dirichlet one (32). The higher gradients render local
minima smooth and gradient discontinuities become smooth transition layers
with width of order $\varepsilon$ [6, 20, 25, 12, 11]. Also there is a global
minimizer, but a careful bifurcation analysis [20] has revealed multitudes of
stable local minima with multiple bands. A comparison of one of our solutions
with one of the equilibria found in [20] is shown in Fig. 1(b), (c).

We explore how the regularization parameter $\varepsilon$ affects the energy.
The results are shown in Fig. 10-12, where $\gamma=0.25,0.5,0.75$
respectively. The energy decreases when the regularization parameter
$\varepsilon$ is smaller, while solutions for smaller $\varepsilon$ have a
finer microstructure with more bands, e.g., Fig. 11, as is also observed in
[20]. Also if we use a larger DNN, the energy will be lower, since the greater
number of parameters allows capturing a local minimum with lower energy and
eventually hopefully a global minimum. In the following, we will explore how
the structure of DNNs affect the results.

(a1)

(a2)

(a3)

(a4)

(b1)

(b2)

(b3)

(b4)

(c1)

(c2)

(c3)

(c4)

Figure 10:  2d regularized problem with SmReLU activation function. Here
$\gamma=0.25$, $iteration=300,000$ and learning rate $\eta=10^{-3}$. Top: the
energy with different (a) NN: $3\times 128$; (b) NN: $5\times 128$; (c) NN:
$7\times 128$. First column: $\varepsilon=\frac{0.1}{4}$; Second column:
$\varepsilon=\frac{0.1}{8}$; Third column: $\varepsilon=\frac{0.1}{16}$; Last
column: $\varepsilon=\frac{0.1}{32}$.

(a1)

(a2)

(a3)

(a4)

(b1)

(b2)

(b3)

(b4)

(c1)

(c2)

(c3)

(c4)

Figure 11:  2d regularized problem with SmReLU activation function. Here
$\gamma=0.5$, $iteration=300,000$ and learning rate $\eta=10^{-3}$. Top: the
energy with different (a) NN: $3\times 128$; (b) NN: $5\times 128$; (c) NN:
$7\times 128$. First column: $\varepsilon=\frac{0.1}{4}$; Second column:
$\varepsilon=\frac{0.1}{8}$; Third column: $\varepsilon=\frac{0.1}{16}$; Last
column: $\varepsilon=\frac{0.1}{32}$.

(a1)

(a2)

(a3)

(a4)

(b1)

(b2)

(b3)

(b4)

(c1)

(c2)

(c3)

(c4)

Figure 12:  2d regularized problem with SmReLU activation function. Here
$\gamma=0.75$, $iteration=200,000$ and learning rate $\eta=10^{-3}$. Top: the
energy with different (a) NN: $3\times 128$; (b) NN: $5\times 128$; (c) NN:
$7\times 128$. First column: $\varepsilon=\frac{0.1}{4}$; Second column:
$\varepsilon=\frac{0.1}{8}$; Third column: $\varepsilon=\frac{0.1}{16}$; Last
column: $\varepsilon=\frac{0.1}{32}$.

####  3.2.4 Effect of the structure of the DNNs

Here, we fix the width of the DNN as 128 and $\varepsilon=\frac{1}{160}$. We
explore how the depth affects the results, which are shown in Fig. 13. Once
again, one layer is insufficient as it does not capture the banded
microstructure. The energy decreases with the increasing depth of the DNN
after some oscillations.

The number of bands more or less increases with increasing depth, the maximum
observed being 7 yellow bands. This is probably because the global minimum is
probably reached, although this is difficult to prove. More bands would
increase the interfacial energy of transition layers, which is proportional to
$\varepsilon$ times the total length of transition layers between the blue and
yellow phases.

We also see a refinement mechanism in Fig. 13(b4) where a yellow band near the
middle is being split almost completely into two bands. The energy eventually
plateaus out as the depth increases.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 13:  2d regularized problem with SmReLU activation function. Here
$\gamma=0.5$, $iteration=300,000$, $\eta=10^{-3}$ and
$\varepsilon=\frac{0.1}{16}$. Top: the energy with different (a) NN: $1\times
128$; (b) NN: $3\times 128$; (c) NN: $5\times 128$; (d) NN: $7\times 128$; (e)
NN: $9\times 128$; (f) NN: $11\times 128$; (g) NN: $13\times 128$; (h) NN:
$15\times 128$; (i) NN: $17\times 128$; (j) NN: $19\times 128$.

In order to explore how the width affects the results, we fix the depth of the
DNN as 5 and $\varepsilon=\frac{1}{160}$. The results are shown in Fig. 14. We
can see the energy and the number of the yellow bars will decrease with
increasing width of the DNN. For a fixed depth of 5 layers of the DNN,
increasing the width beyond $128$ does not result in more than 5 yellow bands
in our solutions.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 14:  2d regularized problem with SmReLU activation function. Here
$\gamma=0.5$, $iteration=300,000$, $\eta=10^{-3}$ and
$\varepsilon=\frac{0.1}{16}$. Top left: the energy with different (a) NN:
$5\times 8$; (b) NN: $5\times 16$; (c) NN: $5\times 32$; (d) NN: $5\times 64$;
(e) NN: $5\times 128$; (f) NN: $5\times 256$; (f) NN: $5\times 512$.

####  3.2.5 Collocation Point Placement

In computed minimizers we observe that the energy density near the boundaries
is higher, whereas we use a random choice of collocation points in the domain
$\Omega$. In order to minimize the energy contribution near the boundary more
effectively, we select adaptively a higher proportion of collocation points
near the top and bottom parts of domain. We divide the domain $\Omega$ into
three parts, $\Omega_{1}=(0,2)\times(0,0.15)$,
$\Omega_{2}=(0,2)\times(0.15,0.85)$ and $\Omega_{3}=(0,2)\times(0.85,1)$. And
in each domain, we randomly choose $N_{1}$, $N_{2}$ and $N_{3}$ points, and
let $N_{1}=N_{3}$. For details see Fig. 15. When we increase the number
$N_{1}=N_{3}$, keeping the sum of $N_{1}$, $N_{2}$ and $N_{3}$ fixed, the
energy decreases.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 15: Adapted collocation points 2d regularized problem with SmReLU
activation function, $\gamma=0.5$, $iteration=300,000$, $\eta=10^{-3}$ and
$\varepsilon=\frac{0.1}{16}$. First figure: the position of collocation point;
second figure: the energy with different (a) $N_{1}=N_{3}=1500$ and
$N_{2}=7000$; (b) $N_{1}=N_{3}=2000$ and $N_{2}=6000$; (c) $N_{1}=N_{3}=3000$
and $N_{2}=4000$; (d) $N_{1}=N_{3}=3500$ and $N_{2}=3000$; (e)
$N_{1}=N_{3}=4000$ and $N_{2}=2000$; (f) $N_{1}=N_{3}=4500$ and $N_{2}=1000$.

###  3.3 Curved Interfaces

We consider the following energy

$W(\nabla
u)=\frac{1}{2}[(\cos(\phi)u_{x}+\sin(\phi)u_{y})^{2}(1-(\cos(\phi)u_{x}+\sin(\phi)u_{y}))^{2}+(-\sin(\phi)u_{x}+\cos(\phi)u_{y})^{2}],$

with Dirichlet boundary condition

$u(x,y)=\gamma(\cos(\phi)x+\sin(\phi)y),\quad\forall(x,y)\in\partial\Omega.$
(34)

This corresponds to rotating one of the minima of $W$ from $(0,1)$ to
$(\cos\phi,\sin\phi)$, where $\phi$ is a fixed parameter. Incompatibility with
the boundary conditions now forces the twin boundaries (interfaces) to bend
and curve, demonstrating the advantage of the mesh-free nature of the method.
It is easy to see that previously we have considered the case with $\phi=0$.

The results are shown in Fig. 16.

(a1)

(a2)

(a3)

(b1)

(b2)

(b3)

Figure 16:  Rotated minimum 2d problem with SmReLU activation function. Here
$iteration=300,000$, $\eta=10^{-3}$ and $\varepsilon=\frac{0.1}{32}$. (a)
$\gamma=0.5$ and $\phi=\frac{\pi}{4}$; (b) $\gamma=0.5$ and
$\phi=\frac{\pi}{8}$. left: $u_{x}$; middle: $u_{y}$; right:
$\cos(\phi)u_{x}+\sin(\phi)u_{y}$.

(a)

(b)

(c)

Figure 17:  Different activation functions in the 2d regularized problem.
Here $iteration=300,000$, $\eta=10^{-3}$ and $\varepsilon=\frac{0.1}{32}$.
(a)Tanh activation function; (b) ReLU activation function; (c) SmReLU
activation function.

##  4 Conclusions and Discussions

This work makes use of the Deep Ritz Method for the first time in order to
solve a nonconvex gradient-energy minimization problem from materials science.
Such problems are challenging in two ways. First, in general they do not
possess a global minimum but only minimizing sequences, as well as a complex
large set of local minima. Second, minimizers, including local ones, are
characterized by complicated geometric microstructures, comprising multiple
layers separated by twin interfaces, which are in effect free boundaries in
the form of gradient discontinuities. These exhibit needle tapering and tip-
splitting topological transitions near boundaries of the domain [20, 25, 11,
12, 39, 21].

A DNN represents the minimizer-candidate, which allows representation of the
energy as a function of the weights and biases (the DNN parameters), and
optimization with respect to them. Our results confirm the ability of the Deep
Ritz method to spontaneously capture the above microstructure complexities of
local or global minimizers. A crucial new ingredient of our method is the
activation function proposed here for the first time. To begin with, we
observe that the ReLU activation function captures exact global minimizers in
a simple unregularized 1D problem, due to its piecewise linear character.
Other activation functions, such as Tanh and Sigmoid, struggle and are indeed
unable to approximate these weak solutions with derivative jumps Fig. 3. For
regularized problems, we introduce a new smoothened version of ReLU which we
term SmReLU; see (27). It is quite successful in producing results that agree
with careful FD computations [20] and captures the correct microstructure Fig.
1 (b), (c), while it outperforms other traditional activation functions, which
fail at this task; see Fig. 17.

We find that it is essential for the DNN to have both sufficient width and
depth to capture local minimizers. In general, increasing the number of layers
(depth) increases the number of twin bands appearing in our solution until a
point where further depth increases does not improve the solution. In that
sense, depth plays a role analogous to mesh size in FE. We find that for fixed
depth, it is not necessary to increase width past a certain point, so here
depth plays a more important role. Equipped with our new SmReLU activation
function, our method has no problem capturing both sharp gradient
discontinuities across twin boundaries and smoothened transition layers. It
also approximates curved interfaces successfully; see Fig. 16. This is because
our approach is mesh-free, which is a great advantage over FEM and FD methods,
especially in the presence of interfacial discontinuities carrying surface
energy, where problems can arise due to mesh misorientation compared to
interfaces [35]. Another advantage of this method is its simplicity, ease of
programming, and the ability to naturally capture nonsmooth solutions with
geometrically complex free interfaces in more than one dimension without
special treatment such as discontinuous Galerkin FE for higher gradients in
nonconvex elasticity [18].

## Acknowledgements

The research of X. Chen is supported by the Ministry of Education, Singapore,
under its Research Centre of Excellence award to the Institute for Functional
Intelligent Materials (I-FIM, project No. EDUNC-33-18-279-V12). The research
of Z. Zhang is supported by Hong Kong RGC grant (projects 17300318 and
17307921), National Natural Science Foundation of China (project 12171406),
and Seed Funding for Strategic Interdisciplinary Research Scheme 2021/22
(HKU).

## References

  * [1] R. Abeyaratne and J. Knowles, A continuum model of a thermoelastic solid capable of undergoing phase transitions, Journal of the Mechanics and Physics of Solids, 41 (1993), pp. 541–571.
  * [2] I. Babuška, The finite element method for elliptic equations with discontinuous coefficients, Computing, 5 (1970), pp. 207–213.
  * [3] J. Ball and R. James, Fine phase mixtures as minimizers of energy, in Analysis and continuum mechanics, Springer, 1989, pp. 647–686.
  * [4] L. Bottou, Large-scale machine learning with stochastic gradient descent, in Proceedings of COMPSTAT’2010, Springer, 2010, pp. 177–186.
  * [5] E. Bronstein, E. Faran, and D. Shilo, Analysis of austenite-martensite phase boundary and twinned microstructure in shape memory alloys: The role of twinning disconnections, Acta Materialia, 164 (2019), pp. 520–529.
  * [6] J. Carr, M. Gurtin, and M. Slemrod, Structured phase transitions on a finite interval, Archive for Rational Mechanics and Analysis, 86 (1984), pp. 317–351.
  * [7] X. Chen, J. Duan, and G. Karniadakis, Learning and meta-learning of stochastic advection–diffusion–reaction systems from sparse measurements, European Journal of Applied Mathematics, 32 (2021), pp. 397–420.
  * [8] Z. Chen and J. Zou, Finite element methods and their convergence for elliptic and parabolic interface problems, Numerische Mathematik, 79 (1998), pp. 175–202.
  * [9] C. Chu, I. Graham, and T. Y. Hou, A new multiscale finite element method for high-contrast elliptic interface problems, Math. Comp., 79 (2010), pp. 1915–1955.
  * [10] C. Chu and R. James, Biaxial loading experiments on Cu-Al-Ni single crystals, in Proceedings of the 1993 ASME Winter Annual Meeting, Publ by ASME, 1993, pp. 61–69.
  * [11] S. Conti, M. Lenz, N. Lüthen, M. Rumpf, and B. Zwicknagl, Geometry of martensite needles in shape memory alloys, Comptes Rendus. Mathématique, 358 (2020), pp. 1047–1057.
  * [12] P. Dondl, B. Heeren, and M. Rumpf, Optimization of the branching pattern in coherent phase transitions, Comptes Rendus Mathematique, 354 (2016), pp. 639–644.
  * [13] W. E, J. Han, and A. Jentzen, Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations, Communications in Mathematics and Statistics, 5 (2017), pp. 349–380.
  * [14] W. E and B. Yu, The deep Ritz method: A deep learning-based numerical algorithm for solving variational problems, Communications in Mathematics and Statistics, 6 (2018), pp. 1–12.
  * [15] J. Ericksen, Equilibrium of bars, Journal of elasticity, 5 (1975), pp. 191–201.
  * [16] E. Georgoulis, M. Loulakis, and A. Tsiourvas, Discrete gradient flow approximations of high dimensional evolution partial differential equations via deep neural networks, arXiv preprint arXiv:2206.00290, (2022).
  * [17] Y. Gong, B. Li, and Z. Li, Immersed-interface finite-element methods for elliptic interface problems with nonhomogeneous jump conditions, SIAM Journal on Numerical Analysis, 46 (2008), pp. 472–495.
  * [18] G. Grekas, K. Koumatos, C. Makridakis, and P. Rosakis, Approximations of energy minimization in cell-induced phase transitions of fibrous biomaterials: $\gamma$-convergence analysis, SIAM Journal on Numerical Analysis, 60 (2022), pp. 715–750.
  * [19] J. He, L. Li, J. Xu, and C. Zheng, Relu deep neural networks and linear finite elements, Journal of Computational Mathematics, 38 (2020), pp. 502–527.
  * [20] T. Healey and A. Sipos, Computational stability of phase-tip splitting in the presence of small interfacial energy in a simple two-phase solid, Physica D: Nonlinear phenomena, 261 (2013), pp. 62–69.
  * [21] T. Hou, P. Rosakis, and P. LeFloch, A level-set approach to the computation of twinning and phase-transition dynamics, Journal of Computational Physics, 150 (1999), pp. 302–331.
  * [22] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Physics-informed machine learning, Nature Reviews Physics, 3 (2021), pp. 422–440.
  * [23] Y. Khoo, J. Lu, and L. Ying, Solving parametric PDE problems with artificial neural networks, European Journal of Applied Mathematics, Special Issue 3: Connections between Deep Learning and Partial Differential Equations, 32 (2021), pp. 421–435.
  * [24] D. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv:1412.6980, (2014).
  * [25] R. Kohn and S. Müller, Branching of twins near an austenite twinned martensite interface, Philosophical Magazine A, 66 (1992), pp. 697–715.
  * [26] I. Lagaris, A. Likas, and D. Fotiadis, Artificial neural networks for solving ordinary and partial differential equations, IEEE Trans. Neural Netw., 9 (1998), pp. 987–1000.
  * [27] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, Nature, 521 (2015), p. 436.
  * [28] R. Leveque and Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM Journal on Numerical Analysis, 31 (1994), pp. 1019–1044.
  * [29] B. Li and M. Luskin, Theory and computation for the microstructure near the interface between twinned layers and a pure variant of martensite, Materials Science and Engineering: A, 273 (1999), pp. 237–240.
  * [30] Z. Li, T. Lin, and X. Wu, New cartesian grid methods for interface problems using the finite element formulation, Numerische Mathematik, 96 (2003), pp. 61–98.
  * [31] Y. Liao and P. Ming, Deep nitsche method: Deep ritz method with essential boundary conditions, arXiv preprint arXiv:1912.01309, (2019).
  * [32] Q. Lou, X. Meng, and G. Karniadakis, Physics-informed neural networks for solving forward and inverse flow problems via the Boltzmann-BGK formulation, Journal of Computational Physics, 447 (2021), p. 110676.
  * [33] R. Mattey and S. Ghosh, A novel sequential method to train physics informed neural networks for Allen Cahn and Cahn Hilliard equations, Computer Methods in Applied Mechanics and Engineering, 390 (2022), p. 114474.
  * [34] H. Montanelli and Q. Du, New error bounds for deep ReLU networks using sparse grids, SIAM Journal on Mathematics of Data Science, 1 (2019), pp. 78–92.
  * [35] M. Negri, The anisotropy introduced by the mesh in the finite element approximation of the Mumford-Shah functional, Numerical Functional Analysis and Optimization, 20 (1999), pp. 957–982.
  * [36] C. Peskin, Numerical analysis of blood flow in the heart, Journal of Computational Physics, 25 (1977), pp. 220–252.
  * [37] M. Raissi, P. Perdikaris, and G. Karniadakis, Multistep Neural Networks for Data-driven Discovery of Nonlinear Dynamical Systems, arXiv:1801.01236, (2018).
  * [38] M. Raissi, P. Perdikaris, and G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics, 378 (2019), pp. 686–707.
  * [39] H. Seiner, P. Plucinsky, V. Dabade, B. Benešová, and R. James, Branching of twins in shape memory alloys revisited, Journal of the Mechanics and Physics of Solids, 141 (2020), p. 103961.
  * [40] Y. Shin, J. Darbon, and G. Karniadakis, On the convergence of physics informed neural networks for linear second-order elliptic and parabolic type PDEs, arXiv preprint arXiv:2004.01806, (2020).
  * [41] S. Wang and P. Perdikaris, Deep learning of free boundary and Stefan problems, Journal of Computational Physics, 428 (2021), p. 109914.
  * [42] Z. Wang and Z. Zhang, A mesh-free method for interface problems using the deep learning approach, Journal of Computational Physics, 400 (2020), p. 108963.
  * [43] Z. Zhang, P. Rosakis, T. Y. Hou, and G. Ravichandran, A minimal mechanosensing model predicts keratocyte evolution on flexible substrates, Journal of the Royal Society Interface, 17 (2020), p. 20200175.


# Revisiting a disky origin for the faint branch of the Sagittarius stellar
stream

Pierre-Antoine Oria  Université de Strasbourg, CNRS, Observatoire astronomique
de Strasbourg, UMR 7550, F-67000 Strasbourg, France  Rodrigo Ibata  Université
de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550,
F-67000 Strasbourg, France  Pau Ramos  Université de Strasbourg, CNRS,
Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France
Benoit Famaey Université de Strasbourg, CNRS, Observatoire astronomique de
Strasbourg, UMR 7550, F-67000 Strasbourg, France  Raphaël Errani  Université
de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550,
F-67000 Strasbourg, France

###### Abstract

We investigate ways to produce the bifurcation observed in the stellar stream
of the Sagittarius dwarf galaxy (Sgr). Our method consists in running $N$-body
simulations of Sgr falling into the Milky Way for the last 3 Gyr, with added
test particles on disk orbits that span a wide range of initial positions,
energies, and angular momenta. We find that particles that end up in the faint
branch are predominantly high angular momentum particles that can all
originate from a single plane within the progenitor, nearly perpendicular both
to the orbital plane of the progenitor and to the Milky Way stellar disk.
Their original configuration at the start of the simulation corresponds to
spiral features already present 3 Gyr ago, which could be, e.g., the result of
a disk-like component being tidally perturbed, or the tidal tails of a
satellite being disrupted within Sgr. We then run a simulation including the
self-gravity of this disky component. Despite the remaining ambiguity of its
origin, this disk component of the Sgr dwarf with spiral over-densities
provides a first step towards a working model to reproduce the observed faint
branch of the bifurcated Sgr stream.

galaxies: dwarf — galaxies: kinematics and dynamics — Local Group

††journal: ApJ

##  1 Introduction

Since its discovery (Ibata et al., 1994, 1995), the Sagittarius dwarf galaxy
(Sgr) has been under intense scrutiny as the closest example of an on-going
galactic merging event. The stellar stream generated by its tidal disruption
(Ibata et al., 2001; Majewski et al., 2003) is an extended and complex
kinematic structure in the stellar halo of the Milky Way (MW), and as such
constitutes an invaluable source of information on the gravitational potential
and history of both the MW and the progenitor dwarf galaxy itself.

Over the years, several models have been put forward in order to reproduce the
shape of the stream and its kinematics (e.g. Ibata et al. 1997; Helmi 2004;
Law et al. 2005). Among those, Law & Majewski (2010) reproduced most of the
observational constraints at the time, involving however an unrealistic and
unstable triaxial dark matter halo configuration for the MW. The latest up-to-
date model by-passing this problem is that of Vasiliev et al. (2021, hereafter
V21), in which the Sgr dwarf is infalling in the joint evolving gravitational
potential of the MW and the Large Magellanic Cloud (LMC), yielding a very good
agreement with recent Gaia data.

One of the remaining mysteries about the Sgr stream is the presence of a
bifurcation, in the form of a faint branch running parallel to the main
brighter branch, observed first in the leading arm (Belokurov et al., 2006),
then in the trailing arm (Koposov et al., 2012). More recently, this
bifurcation has been outlined with great precision by Ramos et al. (2021)
using the latest Gaia EDR3 data (Gaia Collaboration et al., 2021).

Fellhauer et al. (2006) proposed an early explanation for the bifurcation, as
the result of the young leading, old leading and trailing wraps overlapping
and being slightly displaced due to the precession of the orbit (Yanny et al.,
2009; Niederste-Ostholt et al., 2010), but this model did not match later
observations of the stream.

Peñarrubia et al. (2010) then proposed a model in which the Sgr dwarf
originally consisted of a rotating stellar disk embedded in a cold dark matter
halo. A disk slightly misaligned with respect to the orbital plane was shown
to produce a bifurcation as observed in the Sgr stream. However, the model
predicted some remnant rotation in the centre of Sgr today, which was not
observed (Peñarrubia et al., 2011).

Although not in the context of the bifurcation, an originally disky Sgr was
also studied by Łokas et al. (2010) in order to explain the elongated shape of
the remnant. This model makes use of the tidal stirring mechanism (Mayer et
al., 2001; Kazantzidis et al., 2011) according to which dwarf spheroidal
galaxies are the outcome of disky satellite galaxies being deformed due to
galactic tides.

In this letter, we re-investigate the production of a bifurcation by
selecting, within simulations of the Sgr stream, particles that end up in the
observed faint branch, and then examine the properties of the initial
conditions.

##  2 Methodology

###  2.1 Reference model

The underlying model that we use for the present work is the $N$-body
simulation proposed by V21 of the Sgr dwarf falling into a joint, evolving MW
and LMC gravitational potential. This model constitutes an ideal basis for our
investigations as it already reproduces many observational constraints (e.g.
positions, proper motions, distances and line-of-sight velocities, with the
inclusion of the LMC being key for the latter two, especially in the leading
arm), leaving us free to focus our efforts on the production of the
bifurcation. In this context, the simulation starts 3 Gyr ago, at which point
Sgr is made of a spherical King distribution stellar component of mass
$2\times 10^{8}$ M⊙, immersed in a spherical dark halo of mass $3.6\times
10^{9}$ M⊙. The stellar and dark matter components are made of $2\times
10^{5}$ particles each. The MW and LMC models are described in detail in V21.

###  2.2 Sagittarius model and simulation

Using the $N$-body code Gyrfalcon (Dehnen, 2000), we first reproduce the
simulation provided by V21. Then, we add test particles to the initial
conditions of the Sgr dwarf (self-gravity will be included in §3.2) to see
which ones are more likely to end up in the faint branch of the Sgr stream by
the end of the simulation.

Given the aforementioned works hinting strongly at the importance of rotation
in Sgr to produce a bifurcation, we choose to populate our simulations with
test particles with wide ranges of angular momenta. The sample of test
particles is produced by generating stellar disks using Agama (Vasiliev, 2019)
and giving them each a different inclination w.r.t. the orbital plane of Sgr.
Each disk has a scale radius $R_{\rm disk}=0.9$ kpc (V21’s King model has a
scale radius of $1$ kpc), a scale height $H_{\rm disk}=0.18$ kpc, central
velocity dispersion $\sigma_{r,0}=4$ kms-1, and is generated through a
QuasiIsothermal distribution function. The full Agama script for generating
the disk is available in the shared data.

We use a right-handed Cartesian coordinate system centered on the MW with the
$xy$ plane being its disk plane, and the $x$-axis pointing along the Sun-
Galactic center direction, with the sun at $(x,y,z)=(-8,0,0)$ kpc. Our disks
are generated in this plane, then we rotate them before launch by probing
inclination angles $i$ every $20^{\circ}$, both around the $x$-axis and around
the $y$-axis using the following matrices respectively:

$\displaystyle R_{x}=\begin{pmatrix}1&0&0\\\ 0&\cos(i)&-\sin(i)\\\
0&\sin(i)&\cos(i)\end{pmatrix},$ (1) $\displaystyle
R_{y}=\begin{pmatrix}\cos(i)&0&-\sin(i)\\\ 0&1&0\\\
\sin(i)&0&\cos(i)\end{pmatrix}.$ (2)

For reference, for such a disk to be in the orbital plane of Sgr (at present),
it would have to be rotated around the $x$-axis with the $R_{x}$ matrix by an
angle of $i\simeq 100^{\circ}$. Before applying a rotation, our disk models
have angular momentum aligned with the positive $z$-axis. Preliminary tests
showed us that particles in the inner regions of the disk would not end up in
the faint branch, but rather end up close to the remnant of the progenitor.
This is understandable as those particles are deeper in the potential well of
the King model and much harder to strip. We thus select from those disks the
$2\times 10^{4}$ outermost test particles out of the $5\times 10^{4}$
particles generated, allowing us to better probe the regions of interest. This
corresponds to a hole in the inner $\simeq 1.5$ kpc of the disks.

###  2.3 Stream selection

Figure 1: The Gaia EDR3 Sgr sample from Ramos et al. (2021) in the background
(grey), with our selection of faint branch test particles from our simulations
from Section 3.1.1 over-plotted (red for leading arm, blue for trailing arm).
The remnant of the progenitor lies in the $275\leq\rm{RA}\leq 300$ region.

In order to pick the particles that best match the observations of the faint
branch of the Sgr stream, we use the polynomial fits proposed by Ramos et al.
(2021, Table 1) for the ($\tilde{\Lambda}_{\odot}$,$\tilde{\beta}_{\odot}$)
coordinate system centred on Sgr, introduced in Majewski et al. (2003) and
representing the latitude and longitude along its stream. We use the slight
sign modification of Belokurov et al. (2014) for this coordinate system, in
which $\tilde{\Lambda}_{\odot}$ increases towards the leading arm. In the
final snapshot of our simulations, we thus select the test particles with
$\tilde{\beta}_{\odot}$ such that
$|\tilde{\beta}_{\odot}/P(\tilde{\Lambda}_{\odot})-1|<0.2$ where

$P(\tilde{\Lambda}_{\odot})=-0.0003819\tilde{\Lambda}_{\odot}^{2}+0.01904\tilde{\Lambda}_{\odot}+6.084$
(3)

applies to the leading arm part of the faint branch, and

$P(\tilde{\Lambda}_{\odot})=-0.001563\tilde{\Lambda}_{\odot}^{2}-0.2805\tilde{\Lambda}_{\odot}-3.040$
(4)

applies to the trailing arm part of the faint branch. We also require that
$|\tilde{\Lambda}_{\odot}|\geq 20$ in order to exclude the progenitor. In
Figure 1, we show the Gaia EDR3 Sgr sample from Ramos et al. (2021), and we
over-plot what our faint branch selection of Section 3.1 based on Equations
(3) and (4) from our simulations with test particles looks like.

##  3 Results

###  3.1 The faint branch as test particles

Figure 2: Fraction of disk particles that end up in the faint branch for a
given disk inclination $i$. Each inclination angle is probed by a holed disk
made of the $2\times 10^{4}$ outermost test particles, with the holed disk
rotated as explained in Section 2 around the $x$-axis (Sun-MW center axis) for
the top panels, and around the $y$-axis for the bottom panels. The arrow
points to our best model.

We trace our selection of disk particles in the faint stream (see Equations
(3) and (4) and Figure 1) back to the initial conditions.

Figure 2 shows the fraction of test particles that end up in the faint branch
based on initial disk inclination and rotation. We only show the exploitable
results: disks rotated around the $x$-axis by an angle $0^{\circ}\leq i\leq
180^{\circ}$ (top-panels) and disks rotated around the $y$-axis by an angle
$180^{\circ}\leq i\leq 360^{\circ}$ (bottom-panels). We find that other
rotations and angles do not lead efficiently to the creation of a faint
branch, with at best $\simeq 2-3$% of particles ending in the desired regions.

####  3.1.1 Best model: rotation around the $y$-axis

The most appealing model consists of a single disk rotated by an angle
$i=280^{\circ}$ around the $y$-axis, which leads to a faint branch with
roughly the same amount of stars in the leading and trailing parts, as can be
seen in Figure 2. We pick this option as the preferred model in this work. In
this case, the disk almost lies in the $yz$ MW plane, making it nearly
perpendicular both to the MW disk ($\sim xy$ MW plane) and to the Sgr orbital
plane ($\sim xz$ MW plane).

Figure 3: Evolution of the faint branch selection, from spirals at the start
of the simulation 3 Gyr ago (top-left panel, in the initial disk plane, close
to the $yz$ MW plane) to present time (bottom-right panel). The evolution is
seen in the $xz$ MW plane, close to the Sgr orbital plane. Particles that will
make up the final leading arm are shown in red, and the final trailing arm in
blue. The black curve represents the progenitor’s orbit, and the orange ball
represents the Sun. Pericentric passages of the Sgr occur around $t\simeq-2.3$
Gyr and $t\simeq-1.1$ Gyr. Similar plots for every snapshot of the simulation
and the corresponding video are available in the shared data.

We thus run another simulation with the single disk rotated by an angle
$i=280^{\circ}$ around the $y$-axis added to the model of V21, still as test
particles, but keeping this time the full disk, made of $10^{5}$ particles.
After this rotation, our disk has angular momentum nearly aligned with the
positive $x$-axis, with a small positive $z$ component. A majority of this
initially full disk ends up in the bright branch or near the progenitor, and
is not part of our selection. However, picking the disk test particles that
end up in the faint branch once more, we are now interested in their
distribution in the initial conditions. We find that our selection picks out
high energy and angular momentum particles of the disk, and traces spiral arm-
like features (Figure 3, top-left panel) in the outer disk, which would be
sufficient to lead to the creation of the faint branch.

In order to highlight the importance of angular momentum, we note in passing
that doing the same exercise (selecting the faint branch and looking back in
initial conditions) with the stellar particles of the King model of V21 does
not lead to any clear signature distribution in position, energy, or angular
momentum.

####  3.1.2 Alternative model: rotation around the $x$-axis

Disks rotated around the $x$-axis with inclination angles $i=60^{\circ}$ and
$i=80^{\circ}$ are also interesting, with $\simeq 16-17$% of the particles
that end up in the faint branch (Figure 2). This is not too surprising: at
such inclinations, the disk plane matches closely the Sgr orbital plane. In
this configuration, stars in the Sgr disk are on prograde orbits with respect
to the orbit of Sgr around the MW. This has been shown in Łokas et al. (2015)
to lead to stars being stripped easily, producing thin streams.

The model with a disk rotated by an angle $i=70^{\circ}$ around the $x$-axis
is not implausible, but produces slightly worse results than our best disk
model rotated by 280 degrees around the $y$-axis: the trailing arm is harder
to populate, and the agreement with Gaia kinematics is not as good. From the
Sgr orbital plane, the plane of such a disk makes an angle $-30^{\circ}$
around the $x$-axis. This value, which emerges naturally from our probing of
the initial conditions when considering rotations around the $x$-axis, is very
close to the value of $-20^{\circ}$ originally proposed by Peñarrubia et al.
(2010). Data and plots for this model are provided in the repository.

###  3.2 Including self-gravity

Figure 4: Final snapshot of the simulation (at present time) of our model with
massive particles replacing some stellar particles of V21’s initial model. Top
panel: The faint branch (red for leading arm, blue for trailing arm) is well
populated over-plotted on V21’s particles (green). Bottom panel: Density plot
of all stream stars (V21’s stellar particles plus our faint branch particles).

We now study whether the results of the previous section, obtained assuming
that disk stars are mass-less tracers of the underlying potential, also hold
when taking into account the self-gravity between disk stars. This will allow
us to put forward a model that creates a faint branch like the one observed in
the Sgr stream, using the initial conditions and gravitational potential of
V21 as a backbone. In the Gaia EDR3 sample of Sgr of Ramos et al. (2021),
stars with probability $\geq 80$% of being part of the faint branch make up
$\simeq 4$% of the total. We thus aim to be close to this ratio, and replace
6600 of the $2\times 10^{5}$ stellar particles in the V21 model by new ones
following our initial disky distribution. In order to include our particles
into the reference model, we give them the same mass as the stellar particles
of V21, and for each particle that we include, we remove one stellar particle
from V21 sitting at the closest radius from that of our particle. Doing so
ensures that we keep the same total mass, and does not alter the non-linear
dynamics too much.

We follow the evolution of our faint branch selection along the simulation in
Figure 3, from initial spirals to eventually forming the faint branch. Our
selection remains largely bound with angular momentum still pointing towards
the positive $x$ direction until the second pericentric passage (around
$t\simeq-1.1$ Gyr), which strips the faint branch particles from the
progenitor.

Figure 4 (top panel) shows the
$\tilde{\Lambda}_{\odot}$-$\tilde{\beta}_{\odot}$ view of the simulation at
present time. The faint branch can be seen to be well populated, although with
self-gravity now playing a role, a few of our particles end up close to the
progenitor.

An issue is that the thick stream of the V21 model extends to the faint branch
region already, resulting in an overly dense faint branch in Figure 4 (bottom
panel). In a complete bifurcation model, the bright branch should ideally be
thinner, which could probably be achieved with e.g. a different initial
density profile, or a non-equilibrium transitional situation. We leave this
exploration to a future study.

We compute the mean line-of-sight velocity in the remnant of the progenitor to
make sure that our faint branch particles did not perturb the spherical model
of V21 by adding significant rotation. We find a gradient of $\sim 10$ km s-1,
similarly to the pressure-supported model of Peñarrubia et al. (2011, Fig. 2).

Interestingly, we note an over-density of our particles in the
$-180^{\circ}\leq\tilde{\Lambda}_{\odot}\leq-130^{\circ}$ region
($70^{\circ}\leq\rm{RA}\leq 120^{\circ}$), for which we do not have Gaia data.
This signature appeared in all our simulations with disks of inclination close
to that of the chosen model. It would be interesting to see if such an over-
density can be observed.

Finally, we compare in Figure 5 our faint branch particles to the faint branch
selection (probability $\geq 80\%$) from the Gaia EDR3 sample of Ramos et al.
(2021). The radial velocities follow the observed trend for the faint branch
in the trailing and leading arms, and are different from those of the bright
branch in agreement with the data. We remind that our initial selection has
been made purely in configuration space, so that this phase-space agreement is
impressive. For proper motions, the difference in trends between the bright
and faint branches is small in the data, as can be seen in Ramos et al. (2021,
Fig. 3). We note in passing that our faint branch simulation has smaller
scatter than the observed data in proper motions, probably due to both
intrinsic dispersion (both in velocity and distance) and observational
uncertainties. However, transforming our model into star particles with
magnitudes and hence Gaia astrometric uncertainties is far beyond the scope of
this contribution.

Figure 5: Comparison of the faint branch (probability $\geq 80\%$) from the
Gaia EDR3 Sgr sample of Ramos et al. (2021, grey) with the faint branch
selection in our $N$-body model including self-gravity (red for leading arm,
blue for trailing arm). From top to bottom: proper motions in
$\tilde{\Lambda}_{\odot}$, in $\tilde{\beta}_{\odot}$, and line-of-sight
velocities (heliocentric reference frame). For comparison, the bottom panel
shows bright branch members (probability $\geq 80\%$, yellow).

##  4 Discussion

Despite the obvious similarities with a full disk model (Peñarrubia et al.,
2010), there are a few differences. For one, this spiral disky distribution
allows to populate only the faint branch and can be combined with a more
massive spherical model to also populate the bright branch. It also alleviates
the issue of requiring substantial rotation in the progenitor’s remnant at
present time (Peñarrubia et al., 2011). While out of the scope of the present
work to produce a full Sgr model, we discuss possibilities that would lead to
the presence of our faint branch selection and how the inner Sgr could be
populated.

A possible scenario for the presence of this disky spiral distribution 3 Gyr
ago could be a disky dwarf perturbed by tidal effects (Gajda et al., 2017),
and/or having been affected by disk-shocking while crossing the MW-disk. In
addition to the spirals, the rest of the disk could be transitioning to a
pressure supported spheroidal galaxy as in the tidal stirring mechanism (Mayer
et al., 2001). In this context, the inner galactic disk often goes through a
bar perturbation (Kazantzidis et al., 2011; Łokas et al., 2014). It is thus
possible that a bar would be present in Sgr 3 Gyr ago. With tidal heating, the
bar transforms into a diffuse spheroid, part of which would then end up in the
bright branch of the stream, and the rest would form the elongated remnant of
the progenitor that is now observed. This model is attractive because both
branches of the Sgr stream would come from the same stellar population,
consistent with the small difference observed in metallicity between the faint
and bright branches (Ramos et al., 2021).

Another possibility would be that Sgr was already having a substantial
spheroidal component and that only a remaining rotating disk was affected.
Indeed, fitting a full exponential density profile from the surface density
profile of our spiral selection and extrapolating it to the inner disk, we
find that a total disk mass of $2\sim 3\times 10^{7}$ M⊙, or $10\sim 15$% of
the mass of the stellar component in V21’s model would be sufficient (the mass
range depending on the proportion of faint branch stars, $\simeq 4$% of stream
stars in the data). Such a minor disk component would produce a very low
rotation signal in the Sgr remnant at present time.

The spiral features could also be the tidal tails or stellar stream of an
accreted globular cluster or dwarf galaxy inside the Sgr system. This is
however less enticing as it would require the stellar populations of Sgr and
the putative satellite to be fairly similar.

It is also interesting to compare the stripping history and geometry of this
faint branch to full models (see e.g. Ramos et al. 2021, their Fig. 7). As
shown in Figure 3 and in the shared material, our particles for the leading
and trailing faint branch are both mostly stripped during the second
pericenter of the simulation ($t\simeq-1.1$ Gyr). In addition, this stripping
produces a single “upper” faint branch that can be paired with another Sgr
component that would produce the parallel bright branch, as opposed to the
undesired “X-shape” (Ramos et al., 2021) that is usually obtained when
considering inner rotation and/or orbital precession.

##  5 Conclusion

We propose a model for the bifurcation of the Sgr stream in which the faint
branch is populated by stars that were distributed in a disky spiral
distribution within the progenitor 3 Gyr ago, in a plane nearly perpendicular
to both the Sgr orbital plane and the MW disk plane. This pattern emerged here
naturally by probing a large range of initial position, energy and angular
momentum distributions for stellar test particles that end up in both the
leading and trailing parts of the observed faint branch. Populating the faint
branch this way opens the possibility of freely pairing this work with other
Sgr components that would produce the parallel bright branch.

In the context of the tidal stirring mechanism studied in detail in
Kazantzidis et al. (2011) for the formation of dwarf spheroidals, Sgr could
previously have been a disky galaxy which 3 Gyr ago already held a bar (Łokas
et al., 2014) and started the transition from a disky galaxy to a more
isotropic and diffuse one. Low rotational velocity would then be observed
today in the remnant, in agreement with del Pino et al. (2021). The spiral
features could be tidally-induced, bar-induced, or the result of disk shocking
when crossing the MW. Although out of scope for the present study, it would be
very interesting to see if this could be turned into a working model for the
entire Sgr stream.

Another interesting albeit less likely possibility would be that this spiral
distribution is the tidal tail or stellar stream caused by the disruption of a
satellite of the Sgr system. Further observations of the stellar populations
and their detailed chemistry in both the bright and faint branch will likely
provide very useful information in deciding this matter.

## Data Availability

We make available our addition to the model of V21 (DOI
10.5281/zenodo.4300977), including the particles leading to the creation of
the faint branch as well as plots, a movie, and all related scripts, at DOI
10.5281/zenodo.6581185.

The authors thank the referee for an extremely helpful report which led to
significant improvements of the paper. The authors acknowledge funding from
the European Research Council (ERC) under the European Unions Horizon 2020
research and innovation programme (grant agreement No. 834148) and from the
Agence Nationale de la Recherche (ANR projects ANR-18-CE31-0006 and
ANR-19-CE31-0017).

## References

  * Belokurov et al. (2006) Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJ, 642, L137, doi: 10.1086/504797
  * Belokurov et al. (2014) Belokurov, V., Koposov, S. E., Evans, N. W., et al. 2014, MNRAS, 437, 116, doi: 10.1093/mnras/stt1862
  * Dehnen (2000) Dehnen, W. 2000, ApJ, 536, L39, doi: 10.1086/312724
  * del Pino et al. (2021) del Pino, A., Fardal, M. A., van der Marel, R. P., et al. 2021, ApJ, 908, 244, doi: 10.3847/1538-4357/abd5bf
  * Fellhauer et al. (2006) Fellhauer, M., Belokurov, V., Evans, N. W., et al. 2006, ApJ, 651, 167, doi: 10.1086/507128
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1, doi: 10.1051/0004-6361/202039657
  * Gajda et al. (2017) Gajda, G., Łokas, E. L., & Athanassoula, E. 2017, ApJ, 842, 56, doi: 10.3847/1538-4357/aa74b4
  * Helmi (2004) Helmi, A. 2004, ApJ, 610, L97, doi: 10.1086/423340
  * Ibata et al. (2001) Ibata, R., Irwin, M., Lewis, G. F., & Stolte, A. 2001, ApJ, 547, L133, doi: 10.1086/318894
  * Ibata et al. (1994) Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Nature, 370, 194, doi: 10.1038/370194a0
  * Ibata et al. (1995) —. 1995, MNRAS, 277, 781, doi: 10.1093/mnras/277.3.781
  * Ibata et al. (1997) Ibata, R. A., Wyse, R. F. G., Gilmore, G., Irwin, M. J., & Suntzeff, N. B. 1997, AJ, 113, 634, doi: 10.1086/118283
  * Kazantzidis et al. (2011) Kazantzidis, S., Łokas, E. L., Callegari, S., Mayer, L., & Moustakas, L. A. 2011, ApJ, 726, 98, doi: 10.1088/0004-637X/726/2/98
  * Koposov et al. (2012) Koposov, S. E., Belokurov, V., Evans, N. W., et al. 2012, ApJ, 750, 80, doi: 10.1088/0004-637X/750/1/80
  * Law et al. (2005) Law, D. R., Johnston, K. V., & Majewski, S. R. 2005, ApJ, 619, 807, doi: 10.1086/426779
  * Law & Majewski (2010) Law, D. R., & Majewski, S. R. 2010, ApJ, 714, 229, doi: 10.1088/0004-637X/714/1/229
  * Łokas et al. (2014) Łokas, E. L., Athanassoula, E., Debattista, V. P., et al. 2014, MNRAS, 445, 1339, doi: 10.1093/mnras/stu1846
  * Łokas et al. (2010) Łokas, E. L., Kazantzidis, S., Majewski, S. R., et al. 2010, ApJ, 725, 1516, doi: 10.1088/0004-637X/725/2/1516
  * Łokas et al. (2015) Łokas, E. L., Semczuk, M., Gajda, G., & D’Onghia, E. 2015, ApJ, 810, 100, doi: 10.1088/0004-637X/810/2/100
  * Majewski et al. (2003) Majewski, S. R., Skrutskie, M. F., Weinberg, M. D., & Ostheimer, J. C. 2003, ApJ, 599, 1082, doi: 10.1086/379504
  * Mayer et al. (2001) Mayer, L., Governato, F., Colpi, M., et al. 2001, ApJ, 547, L123, doi: 10.1086/318898
  * Niederste-Ostholt et al. (2010) Niederste-Ostholt, M., Belokurov, V., Evans, N. W., & Peñarrubia, J. 2010, ApJ, 712, 516, doi: 10.1088/0004-637X/712/1/516
  * Peñarrubia et al. (2010) Peñarrubia, J., Belokurov, V., Evans, N. W., et al. 2010, MNRAS, 408, L26, doi: 10.1111/j.1745-3933.2010.00921.x
  * Peñarrubia et al. (2011) Peñarrubia, J., Zucker, D. B., Irwin, M. J., et al. 2011, ApJ, 727, L2, doi: 10.1088/2041-8205/727/1/L2
  * Ramos et al. (2021) Ramos, P., Antoja, T., Yuan, Z., et al. 2021, arXiv e-prints, arXiv:2112.02105.  https://arxiv.org/abs/2112.02105
  * Vasiliev (2019) Vasiliev, E. 2019, MNRAS, 482, 1525, doi: 10.1093/mnras/sty2672
  * Vasiliev et al. (2021) Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279, doi: 10.1093/mnras/staa3673
  * Yanny et al. (2009) Yanny, B., Newberg, H. J., Johnson, J. A., et al. 2009, ApJ, 700, 1282, doi: 10.1088/0004-637X/700/2/1282