Physics-Informed Neural Networks for Discovering Localised Eigenstates
  in Disordered Media
The Schr\"{o}dinger equation with random potentials is a fundamental model
for understanding the behaviour of particles in disordered systems. Disordered
media are characterised by complex potentials that lead to the localisation of
wavefunctions, also called Anderson localisation. These wavefunctions may have
similar scales of eigenenergies which poses difficulty in their discovery. It
has been a longstanding challenge due to the high computational cost and
complexity of solving the Schr\"{o}dinger equation. Recently, machine-learning
tools have been adopted to tackle these challenges. In this paper, based upon
recent advances in machine learning, we present a novel approach for
discovering localised eigenstates in disordered media using physics-informed
neural networks (PINNs). We focus on the spectral approximation of Hamiltonians
in one dimension with potentials that are randomly generated according to the
Bernoulli, normal, and uniform distributions. We introduce a novel feature to
the loss function that exploits known physical phenomena occurring in these
regions to scan across the domain and successfully discover these eigenstates,
regardless of the similarity of their eigenenergies. We present various
examples to demonstrate the performance of the proposed approach and compare it
with isogeometric analysis.
Liam Harcombe, Quanling Deng
# Physics-Informed Neural Networks for Discovering Localised Eigenstates in Disordered Media

###### Abstract

The Schrodinger equation with random potentials is a fundamental model for understanding the behavior of particles in disordered systems. Disordered media are characterised by complex potentials that lead to the localisation of wavefunctions, also called Anderson localisation. These wavefunctions may have similar scales of eigenenergies which poses difficulty in their discovery. It has been a longstanding challenge due to the high computational cost and complexity of solving the Schrodinger equation. Recently, machine-learning tools have been adopted to tackle these challenges. In this paper, based upon recent advances in machine learning, we present a novel approach for discovering localised eigenstates in disordered media using physics-informed neural networks (PINNs). We focus on the spectral approximation of Hamiltonians in one dimension with potentials that are randomly generated according to the Bernoulli, normal, and uniform distributions. We introduce a novel feature to the loss function that exploits known physical phenomena occurring in these regions to scan across the domain and successfully discover these eigenstates, regardless of the similarity of their eigenenergies. We present various examples to demonstrate the performance of the proposed approach and compare it with isogeometric analysis.

keywords: Schrodinger equation, Hamiltonian, Anderson localisation, eigenvalues and eigenstates, neural networks, isogeometric analysis +
Footnote †: journal: arXiv

## 1 Introduction and Background

Partial differential equations (PDEs) are powerful tools for studying dynamic systems and have applications across various fields. Despite their usefulness, these equations are notoriously difficult to solve explicitly, leading to an increased interest in numerical approximations. In addition to classical numerical
methods, recent research has shown the potential of using neural networks to estimate solutions to differential equations [26; 27; 31]. There are several advantages of using neural networks to solve PDEs over traditional numerical methods. Firstly, numerical errors are usually not accumulated [33]. Secondly, they are more robust against the "curse of dimensionality," which refers to the difficulty of solving PDEs in high-dimensional spaces using traditional numerical methods [19]. Furthermore, they can be trained to solve a wide range of PDEs, including those with complex boundary conditions or nonlinear operators [16]. Once trained, neural networks can quickly produce solutions for new PDEs that are similar to those encountered during training, potentially without requiring further adjustments or recalibration of parameters [13]. Overall, neural networks can provide a powerful and flexible tool for solving PDEs, especially in cases where traditional numerical methods may be insufficient or computationally expensive.

Data-driven supervised networks [15] and data-free unsupervised networks [23] have been shown to produce efficient approximations of differential eigenvalue problems. Supervised networks aim to discover patterns in labelled datasets to infer a general model that is able to predict examples outside of the given dataset. The training process involves calculating approximations for given inputs and computing how much these estimations deviate from the true outputs using a loss function. Then, the networks work to minimise this loss function by adjusting its estimation process. Artificial Neural Networks (ANNs) are comprised of layers of interconnected nodes with weights and biases that calculate the estimations. Convolutional Neural Networks (CNNs) are a newer form of neural networks that are especially known for their ability to identify patterns faster and more reliably than ANNs in solving eigenvalue problems using labelled datasets [15]. They introduce a convolutional layer built up of windows that slide around the input data, searching for patterns in whole regions at a time. They are commonly used to pass information on to ANNs and are optimised similarly. In situations where training data is unavailable or difficult to obtain, unsupervised networks are employed as a solution, operating without a labelled dataset. In these situations, more attention is drawn to the design of the network's loss function, as there is no labelled data to compare predictions with [7].

In the case of finding the eigenstates of Hamiltonians, we leverage known physical phenomena to design our loss function. This kind of model is an example of a Physics Informed Neural Network (PINN, c.f., [4; 35]), one that embeds the knowledge of physical laws governing the system into the learning process to drive the network to an admissible approximation.

Differential eigenvalue problems occur in a wide range of problems in physics and applied mathematics, such as quantum energy problems. In early works, Lagaris et al. [27] proposed an ANN to discover the solu
tions to differential eigenvalue problems, tested on various problems in quantum mechanics. More recently, Finol et al. [15] presented a supervised ANN to solve eigenvalue problems in mechanics, showing that CNNs were outperforming traditional ANNs in contexts with labelled datasets. Sirignano and Spiliopoulos [36] developed a deep learning network to accurately solve partial differential equations in dimensions as high as 200. They also proposed a mesh-free algorithm, which was desirable as meshes become infeasible in such high dimensions. Chen et al. [6] adopted the emerging PINNs to solve representative inverse scattering problems in photonic metamaterials and nano-optic technologies. Their method was also mesh-free and was tested against numerical simulations based on the finite element method. These ANN models are excellent at learning time series data. The work [32] demonstrated its application in wave packet displacements in localised and delocalised quantum states. Another work is [25], where they adopted CNNs to solve a classification problem. The authors constructed a supervised model learning from experimental data to distinguish between the dynamics of an Anderson insulator and a many-body localised phase. Yang et al. [38] proposed a combination of PINNs with Bayesian Neural Networks, which solved both forward and inverse nonlinear problems, obtaining more accurate predictions than PINNs in systems with large noises.

Jin et al. [23] have shown that PINNs can solve single and multiple square well problems, calculating the eigenfunction and eigenvalue simultaneously. Their network employs a scanning mechanism that pushes the network to search for higher eigenvalues as the training process evolves, while the network updates the eigenfunction prediction accordingly. They store found eigenfunctions and exploit the orthogonality of wave functions to search for higher eigenstates. Grubisic et al. [18] applied PINNs to identify localised eigenstates of operators with random potential. They studied the effective potential of the operator, whose local minima correspond to the different eigenstates. They first built a PINN to solve these eigenvalue problems in one dimension, then generalised to a deeper model that solves higher dimensional problems. Effective potentials provide a neat way of locating the localised eigenstates, however, are limited in the number of locations they can identify, and have difficulty differentiating between eigenstates of identical eigenvalue. These limitations occur in Bernoulli distributed potentials, which we aim to solve with our model.

We build on the design of [23] to approximate these eigenstates for Hamiltonians whose potential is distributed randomly, leading to localisation of the solution to specific regions of the domain. In particular, the core of our work is in finding these eigenstates where the eigenenergies are nearly identical, a task that current models are unable to accomplish. This is most prevalent with potentials distributed according to the Bernoulli distribution but can occur in other distributions.

The rest of the paper is organised as follows. In Section 2, we state the Schrodinger differential eigen
value problem and present the isogeometric analysis, followed by a discussion on the challenges of its spectral approximation. In section 3, we propose our novel loss function design for the PINN for solving the Schrodinger equation. Section 4 collects and discusses various numerical tests to demonstrate the performance of the proposed method. Concluding remarks are presented in section 5.

## 2 Problem Statement and Isogeometric Analysis

In this section, we first introduce the modeling problem, followed by a presentation of a classic numerical method, namely isogeometric analysis, used to solve the problem. We then provide a general discussion on the challenges associated with the numerical computation of the model.

### Problem Statement

We study the time-independent Schrodinger equation: Find the eigenpair \((E,u)\) such that

\[\begin{split}\left[-\frac{\hbar^{2}}{2m}\Delta+V\right]u& =Eu\quad\text{in}\quad\Omega,\\ u&=0\quad\text{on}\quad\partial\Omega,\end{split} \tag{2.1}\]

where \(\Delta=\nabla^{2}\) is the Laplacian, \(V=V(x)\in L^{2}(\Omega)\) specifies the potential and is a non-negative function, and \(\Omega\subset\mathbb{R}^{d},d=1,2,3,\) is a bounded open domain with Lipschitz boundary \(\partial\Omega\). Throughout the paper, we will focus on the case with \(d=1\), followed by an example of an extension to \(d=2\) in Section 4.3. The differential operator is referred to as the Hamiltonian, i.e., \(\mathcal{H}=-\frac{\hbar^{2}}{2m}\Delta+V\). Here, \(\hbar\) is the reduced Planck constant and \(m\) is the mass of the particle. \(\frac{\hbar^{2}}{2m}\) is a constant that we can divide throughout equation 2.1 and absorb into \(V\) and \(E\). Thus, finding the eigenpairs in equation 2.1 is equivalent to finding those for the following equation, up to a scalar:

\[\begin{split}-\Delta u+Vu&=Eu\quad\text{in}\quad \Omega,\\ u&=0\quad\text{on}\quad\partial\Omega.\end{split} \tag{2.2}\]

This problem is a Sturm-Liouville eigenvalue problem (\(-\Delta+V\) is a diffusion-reaction operator, see, for example, [14; 37]) which has a countable infinite set of positive eigenvalues \(E_{j}\in\mathbb{R}^{+}\)

\[0<E_{1}<E_{2}\leq\cdots\leq E_{j}\leq\cdots \tag{2.3}\]

with an associated set of orthonormal eigenfunctions \(u_{j}\)

\[(u_{j},u_{k})=\int_{\Omega}u_{j}(x)u_{k}(x)\;\text{d}\mathbf{x}=\delta_{jk}, \tag{2.4}\]
where \(\delta_{jk}\) is the Kronecker delta which is equal to 1 when \(j=k\) and 0 otherwise. The set of all the eigenvalues is the spectrum of the Hamiltonian. We normalise the eigenfunctions in the \(L^{2}\) space; hence, the eigenfunctions are orthonormal under the scalar inner product. Let us define two bilinear forms

\[a(w,v)=\int_{\Omega}\nabla w\cdot\nabla v+Vwv\;\text{d}\mathbf{x}\quad\text{and} \quad b(w,v)=(w,v)=\int_{\Omega}wv\;\text{d}\mathbf{x},\quad\forall w,v\in H^{1}_{ 0}(\Omega), \tag{2.5}\]

where \(H^{1}_{0}(\Omega)\) is a Sobolev space with functions vanishing at the boundary \(\partial\Omega.\) Using this notation, the eigenfunctions are also orthogonal with each other with respect to the energy's inner product, i.e.,

\[a(u_{j},u_{k})=E_{j}(u_{j},u_{k})=E_{j}\delta_{jk}. \tag{2.6}\]

We remark that these orthogonalities are critical in the development of the proposed neural networks in Section 3.

### Isogeometric Analysis

At the continuous level, the weak formulation for the eigenvalue problem (2.2) is: Find all eigenvalues \(E\in\mathbb{R}^{+}\) and eigenfunctions \(u\in H^{1}_{0}(\Omega)\) such that,

\[a(w,u)=Eb(w,u),\quad\forall\;w\in H^{1}_{0}(\Omega), \tag{2.7}\]

while at the discrete level, the isogeometric analysis (IGA, c.f., [8; 10; 12; 22]) for the eigenvalue problem (2.2) is: Find all eigenvalues \(E_{h}\in\mathbb{R}^{+}\) and eigenfunctions \(u_{h}\in W_{h}\) such that,

\[a(w_{h},u_{h})=E_{h}b(w_{h},u_{h}),\quad\forall\;w_{h}\in W_{h}(\Omega), \tag{2.8}\]

where \(W_{h}\subset H^{1}_{0}(\Omega)\) is the trial and test space spanned by the B-spline or non-uniform rational basis spline (NURBS) basis functions [9].

In this paper, for the purpose of comparison with the state of the art, we adopt the recently-developed soft isogeometric analysis (softlGA, c.f., [11; 29]). SoftlGA has a similar variational formulation as (2.8) which leads to the matrix eigenvalue problem

\[\mathbf{K}\mathbf{U}=E_{h}\mathbf{M}\mathbf{U}, \tag{2.9}\]

where \(\mathbf{K}_{jk}=a(\phi_{j},\phi_{k}),\mathbf{M}_{jk}=b(\phi_{j},\phi_{k}),\) and \(\mathbf{U}\) contains the coefficients of the eigenfunction \(u_{h}\), for its representation in the linear combination of the B-spline basis functions. For simplicity, the matrix \(\mathbf{K}\) (although it contains a scaled mass) is referred to as the stiffness matrix while the matrix \(\mathbf{M}\) is referred as the mass matrix, and \((E_{h},\mathbf{u}_{h})\) is the unknown eigenpair. We refer to [11; 29] for details.
### A Few Challenges on Numerical Computation

There are a few challenges in the numerical computation, such as using IGA or softIGA, of the Schrodinger equation with random potentials. Firstly, there is a non-uniqueness of solutions. The Schrodinger equation with random potentials can have multiple solutions that correspond to the same (or very similar) energy level. This non-uniqueness can make it difficult to accurately identify the correct eigenstates, particularly when dealing with complex potential landscapes. Secondly, it can be a multi-scale problem with large sampling/discretisation errors. Random potentials can vary over multiple length scales, which leads to inaccurate numerical simulations and can make it challenging to choose an appropriate discretisation size and mesh. Thirdly, the Schrodinger equation for many-body problems in multiple dimensions suffers from the curse of dimensionality. The number of parameters needed to describe the random potential can be very large. As the number of dimensions increases, the computational cost of solving the Schrodinger equation can increase exponentially.

Last but not least, the Schrodinger equation with random potentials is a complex mathematical problem, and solving it numerically can be computationally expensive. This is particularly true for large system sizes and high disorder strengths, which require a large number of numerical simulations. Moreover, solving the Schrodinger equation with a different random potential requires constructing a new matrix eigenvalue problem (2.9), which can be a time-consuming process. This can become particularly challenging when solving the equation for a large number of random potentials. In such cases, one potential solution is to train a neural network that takes the random potential as input and produces the corresponding eigenstates as outputs. For instance, to solve the Schrodinger equation with \(N\) different random potentials, one could use a classic method like softIGA and apply it \(N\) times, or alternatively, train a neural network using data from a subset of cases and then use the well-trained neural network to solve the remaining cases. This approach can save a significant amount of computational time, particularly when \(N\) is large. With these challenges in mind, we introduce the following neural-network-based method as an alternative method to solve the time-independent Schrodinger equations.

## 3 The Physics-Informed Neural Network

Deep neural networks are a type of machine learning algorithm inspired by the structure and function of the human brain. They are composed of multiple layers of artificial neurons, each performing a nonlinear transformation of the input data. The output of each layer is fed as input to the next layer, allowing the network to learn increasingly complex features and relationships in the data.
Some common types of deep neural networks include convolutional neural networks (CNNs) for image processing and computer vision tasks [28], long short-term memory models (LSTMs) for sequential data processing [21], deep belief networks (DBNs) for unsupervised learning [20], and generative adversarial networks (GANs) for image and video synthesis [17]. Deep neural networks have demonstrated impressive performance in many applications. In this paper, we adopt the physics-informed neural networks (PINN) with a novel loss function design to solve the Schrodinger equations.

### Pinn

The Physics-Informed Neural Networks (PINNs) are a type of feed-forward neural networks (FNNs) that incorporate physics laws, typically in the form of partial differential equations (PDEs), into their loss functions [30; 35]. This combination of neural networks and PDEs makes PINNs a powerful tool for solving complex PDEs in various scientific and engineering applications. The overall process of using PINNs to solve a PDE can be summarized as follows. First, a neural network architecture is constructed. This involves defining the structure of the network, including the number and activation functions of its layers, as well as the number of neurons in each layer. Next, the loss function is defined. The loss function quantifies the discrepancy between the neural network's predictions and the actual solution to the PDE. In the case of PINNs, the loss function incorporates terms that enforce the PDE itself, as well as any associated initial or boundary conditions. Once the neural network and loss function are established, the network is trained using an optimization algorithm, such as the Adam optimiser [24]. During the training process, the neural network adjusts its weights and biases to minimise the loss function, gradually improving its ability to approximate the PDE solution. After training, the accuracy of the neural network's predictions can be evaluated by comparing them against calculations made by current numerical methods, such as isogeometric analysis [8]. If the predictions are not satisfactory, further refinement can be performed. This may involve modifying the network architecture, adjusting the loss function, or selecting a different optimization algorithm. The refined network is then retrained to improve its performance. Our PINN implementation is summarised in Algorithm 1.

Typically, the loss function in PINN is defined as

\[L=L_{\text{de}}+L_{\text{reg}}, \tag{3.10}\]

where \(L_{\text{de}}\) specifies the loss associated with the PDE and \(L_{\text{reg}}\) specifies the loss associated with the regularisation such as the boundary conditions. The loss function term \(L_{\text{de}}\) involves derivatives which are usually evaluated by the automatic differentiation by Tensorflow [1] or by Pytorch [34]. In the following subsection,
we present a goal-oriented loss function for the best performance in solving the Schrodinger equations with random potentials. A scale factor is assigned to each regularisation term in \(L_{\text{reg}}\). These scale factors are considered hyperparameters and are tuned to emphasise certain constraints that the network should adhere to in its prediction.

### Loss Function Design

To demonstrate the main idea, we simplify our focus to the one-dimensional case. We assume that the potential \(V(x)\) is randomly distributed over \(m\) regions over the domain \([0,1]\), and each region has a random value. Each set of random values then characterises the individual differential eigenvalue problem to be solved for Anderson localised states. We partition the interval \([0,1]\) into \(n+1\) equally spaced points \(x_{j}=jh,j=0,1,\cdots,n\) where \(h=1/n\) is the grid/mesh size. We denote by \(\mathbf{u}_{h}=(u_{h}(x_{0}),u_{h}(x_{1}),\cdots,u_{h}(x_{n}))^{T}\) as a vector of the approximate values of the true solution \(u(x)\) at each of these coordinates \(x_{j}\). Similarly, we denote by \(\mathbf{V}=(V(x_{0}),V(x_{1}),\cdots,V(x_{n}))^{T}\) as a vector of the approximate values of the potential \(V(x)\). In our PINN model, we choose \(n\) such that each region has exactly 5 nodes each for accuracy as well as for training efficiency. To approximate the second derivative vector \(u^{\prime\prime}(x)\), we use the center finite difference method with an accuracy order of six, i.e., the error is of \(\mathcal{O}(h^{6})\).

We construct a PINN that splits into two separate networks: one that calculates the eigenvalue \(E\), and the other one produces the eigenfunction \(u(x)\) at certain nodes. These two networks are trained and optimised simultaneously, with their outputs used in the same loss function described below. We consider the combination of these two networks as one network and apply it to discover the eigenstates. Figure 3.1 draws this model structure, highlighting the split into separate networks. The network section finding the eigenvalue has two hidden layers with \(n\) nodes, while the network for finding the eigenfunction has three hidden layers with \(n\) nodes. This allows the network to scale with the size of the mesh. Each layer in both networks uses the ReLU activation function \(f(x)=\text{max}(x,0)\), and the loss minimisation is done using the Adam optimiser with a learning rate of \(5\times 10^{-4}\).

The driving force of this network towards a correct solution is the design of its loss function, which is built of multiple terms relating to different aspects of a correct solution. The neural network's optimisation function works to minimise the loss function, so we choose terms that are positive errors for incorrect predictions of \(E\) and \(u(x)\), and are zero for correct predictions. Thus, the minimisation of these terms leads the network to produce correct predictions for \(E\) and \(u(x)\). During training, we choose to input an initial vector of ones (of length \(n\)), which the network uses as a basis to build its predictions of \(E\) and \(u(x)\) from.
The central term of our loss function is the term:

\[L_{\text{de}}=\sqrt{\frac{1}{n+1}\sum_{i=0}^{n}\left(E_{h}u_{h}(x_{i})+u_{h}^{ \prime\prime}(x_{j})-V(x_{j})u_{h}(x_{j})\right)^{2}} \tag{3.11}\]

which is always non-negative and is zero for solutions \(E\) and \(u(x)\) that satisfy equation 2.2. To encourage the network to satisfy the boundary conditions \(u(0)=u(1)=0\), we create the term:

\[L_{\text{bound}}=u_{h}^{2}(x_{0})+u_{h}^{2}(x_{n}) \tag{3.12}\]

which the minimisation of the loss function leads to \(u_{h}(x_{0})=u_{h}(x_{n})=0\). The major issue with using just these terms is that the trivial solution \(u(x)=0\) satisfies these terms (reducing them to zero). Thus, we need the following term to penalise the approximate solution:

\[L_{\text{norm}}=\left(\int_{0}^{1}u_{h}^{2}(x)\text{d}\mathbf{x}-1\right)^{2}. \tag{3.13}\]

The minimisation of the loss function leads to approximate solutions of \(L^{2}\) norm 1. Integration is calculated using the mid-point rule and this calculation is majorly local. In the PINN model developed in [23], this

Figure 3.1: Neural network structure for the prediction of eigenstates of Hamiltonians.
condition was imposed using the terms like \(1/u^{2}(x)\) and \(1/E^{2}\). We point out that since the eigenstates are localised, an eigenstate \(u(x)=0\) for most of the domain, leads to very large and even infinite values due to the design of the term \(1/u^{2}(x)\). Consequently, this leads to difficulties in minimising the loss function and the overall neural network training.

These terms are sufficient to push the network to generate an admissible eigenstate. However, the model often produces the same eigenstates, preventing the discovery of higher states. To overcome this issue, once an eigenstate \((E_{h},u_{h}(x))\) is discovered (using a patience condition described below), we add the eigenstate \(u_{h}(x)\) to a list of eigenstates \(S\). Then, we stop the training, recompile the network and begin training again, where we take advantage of the known physical fact that the eigenfunctions of Hamiltonians are orthogonal and incorporate the following term into the loss function:

\[L_{\text{orth}}=\sum_{\hat{u}_{h}(x)\in S}\left(\int_{0}^{1}\hat{u}_{h}(x)u_{h }(x)dx\right)^{2}, \tag{3.14}\]

where \(u_{h}(x)\) is the network's current state to be produced. We iteratively repeat this process of adding eigenstates to \(S\) and recompiling the network. This iteration is automated and can be done an arbitrary number of times. Sometimes the network will skip an eigenstate and converge to a higher one, however, with enough iterations the network tends to find these states eventually. To decide whether a prediction is admissible, we determine whether \(L_{\text{de}}\) is below a certain threshold, as well as a patience condition as in [23] that checks if the absolute value of the change in \(L_{\text{de}}\) and that of the eigenvalue is below another threshold. This allows the network to converge to solutions with \(L_{\text{de}}\) much lower than the chosen threshold, given that it is converging fast enough.

The loss function we have built so far performs adequately on problems with enough spread in the eigenvalues between different states. However, in situations where the eigenstates have nearly identical eigenvalues, the network is unable to distinguish between the states, converging to an undesirable linear combination of the eigenfunctions. This occurs because eigenspaces form vector spaces, so a linear combination of eigenfunctions of the same eigenvalue \(E\) will be another eigenfunction of eigenvalue \(E\), giving another eigenstate that would satisfy \(L_{\text{de}}\) and all the regularisation terms defined. When studying Hamiltonians with potentials generated randomly, there is a disorder for Anderson localisation to occur. We leverage this fact in our loss function with the following procedure.

We let the model run with the above loss function for a certain number of epochs \(q\), which is a hyperparameter. The goal is to set this parameter so that at epoch \(q\), the network has converged to a linear
combination of the eigenfunctions, for example in Figure 3.2.

Since we know the eigenstates of the Hamiltonian are localised, each spike in 3.2 represents a separate localised state. We thus scan through the eigenfunction and generate a list \(K\) of the regions of the domain where the norm of the eigenfunction is greater than some set hyperparameter, splitting apart the regions where the function returns to zero to separate the localised spikes. Then, in our iterative process of discovering eigenvalues, at each step, we choose one of the intervals in \(K\) (say, \([x_{a},x_{b}]\)) and add the following term to the loss function:

\[L_{\text{loc}}=\left(\int_{0}^{1}u_{h}^{2}(x)\text{d}\mathbf{x}-\int_{x_{a}}^{x_{b }}u_{h}^{2}(x)\text{d}\mathbf{x}\right) \tag{3.15}\]

which encourages the network to set the eigenfunction to zero outside of the localisation interval \([x_{a},x_{b}]\). Equivalently, we could have added terms like \(L_{\text{end}}\) for each nodal point outside of \([x_{a},x_{b}]\). We remark that this term is non-negative. With all these loss terms in mind, we have

\[L_{\text{reg}}=\alpha_{\text{bound}}L_{\text{bound}}+\alpha_{\text{norm}}L_{ \text{norm}}+\alpha_{\text{orth}}L_{\text{orth}}+\alpha_{\text{loc}}L_{\text {loc}}, \tag{3.16}\]

which is added to \(L_{\text{de}}\) to give the overall loss function \(L\) in (3.10). The \(\alpha\) terms in \(L_{\text{reg}}\) scale each term in the loss function, providing emphasis to certain physical constraints. Through hyperparameter tuning, we found that the following values provided adequate predictions for the number of nodes we used (order of magnitudes \(10^{2}\) to \(10^{3}\))

\[\alpha_{\text{bound}}=n^{2},\quad\alpha_{\text{norm}}=n^{3},\quad\alpha_{ \text{orth}}=n^{3},\quad\alpha_{\text{loc}}=n^{4}. \tag{3.17}\]

Figure 3.2: Left: An example of a linear combination of eigenstates where their eigenenergies are approximately equal. Right: The corresponding potential \(V(x)\) with \(m=80\) uniform elements, where \(V(x)\) was randomly generated in each element according to the normal distribution with mean 1 and standard deviation 0.3, scaled by \(10^{6}\).
The idea of having a large weight on the regularisation terms of the loss function is to ensure the network prioritises minimising these terms before working to a solution satisfying the differential equation (i.e. before minimising \(L_{\text{de}}\)). This also penalises the network heavily for straying its solution from the regularisations in \(L_{\text{reg}}\), working to project the parameter search space to a subspace adhering to the regularisations. The method described above is presented in Algorithm 1. The purpose of this algorithm is to discover the eigenstates for a given potential. Throughout the paper, we use this loss function in our proposed PINN with the structure shown in Figure 3.1.

## 4 Numerical Experiments

We test our PINN model on various Hamiltonians with randomly distributed potential and demonstrate the network's robustness against potential with different distributions. In particular, we consider potentials with distributions that cause the Hamiltonian to admit eigenstates with almost identical eigenenergies. This is most prevalent in potentials distributed according to the Binomial distribution. We test the network for different numbers of randomly generated regions (\(m\)), and choose the number of nodes in the mesh (\(n\)) such that each region contains 5 nodes (that is, we choose \(n=5m\)).

The localised states usually occur in the region where the potential takes its smallest value. When there are multiple disjoint regions at (or very close to) this minimum value, each of these regions can admit a localised state, where they all can have similar eigenenergy. Figure 4.3 shows the first four PINN approximated eigenstates where the potential is shown in the right plot of Figure 3.2. The network can distinguish between the first three eigenstates, regardless of how similar their eigenvalues are.

We observe that if \(u(x)\) satisfies equation 2.2, then \(-u(x)\) also does. \(-u(x)\) also has the same \(L^{2}\) norm as \(u(x)\), so it satisfies all other terms in our loss function. This is relevant as the network will choose to converge to either \(u(x)\) or \(-u(x)\). For consistency, we plot all eigenfunctions in their positive form.

Figure 4.4 shows a potential with a Bernoulli distribution, while Figure 4.5 shows the first four eigenstates, the first three of which have almost identical eigenvalues. We expect these eigenstates to localise in the "troughs" of the potential, which is what we observe in the eigenfunction plots. The states that localise in troughs of similar length tend to have almost the same eigenvalue, which our model captures.

Figure 4.6 shows the model's prediction of the eigenvalue as the training process evolves. After 2000 epochs, the model adds the \(L_{\text{loc}}\) term. The model produces converging eigenvalues until around epoch 9000. At this point, the model's prediction has met the convergence criteria and saves the eigenstate to use in \(L_{\text{orth}}\) and it resets its weights to search for the next eigenstate. In our implementation, we include a function that
```
1: Let \(k\) be the epoch where \(L_{\text{loc}}\) is added to the loss function
2: Let \(j=1\)
3:while\(j\leq\) num_eigenstates do
4: Initiate the model
5:while training do
6: Calculate prediction of \(E\) and \(u(x)\)
7: Compute \(L_{\text{de}}\), \(L_{\text{bound}}\), \(L_{\text{norm}}\)
8: Compute \(L_{\text{orth}}\) using all stored eigenfunctions
9:if current epoch \(=k\)then
10: Scan predicted eigenfunction for localisation sections
11: Choose the localisation section \(S\) with greatest norm
12: Reset model weights
13:elseif current epoch \(>k\)then
14: Compute \(L_{\text{loc}}\) across \(S\)
15:endif
16: Backpropagate and step
17:if patience condition and \(L_{\text{de}}<\) threshold then
18: Store prediction of \(E\) and \(u(x)\) for future \(L_{\text{orth}}\) calculation
19:\(j=j+1\)
20:break while
21:endif
22:endwhile
23:endwhile
```

**Algorithm 1** Eigenstate Discovery Using PINN
generates a video to display the network's eigenfunction prediction as the training process evolves, which is available at our GitHub1 page.

Footnote 1: [https://github.com/liamharcombe4/eigen-network](https://github.com/liamharcombe4/eigen-network)

Table 1 displays the PINN eigenvalue errors where we used the reference eigenvalues generated by the softIGA method with a very fine mesh (2000 nodes). We also test the model's robustness against the number of nodes. We note that the model achieves slightly better results for the potential distribution according to the normal distribution compared to the uniform distribution. However, the results for the Bernoulli distributed potentials are much better, all with relative errors below 1%. This demonstrates the success of our approach to the problem of eigenstates with similar eigenvalue. Table 2 displays the eigenvalue prediction results for a single potential generated by the Bernoulli distribution with a different number of mesh elements. As the number of elements increases, there is a clear improvement in the accuracy of the PINN-predicted

Figure 4.3: The first four PINN approximated eigenstates with the potential defined in the right plot of Figure 3.2. Here, ev refers to the eigenvalue associated with the plotted eigenfunction.
eigenvalues. The convergence rate is of order approximately 1 but the theoretical analysis of this convergence is subject to future study.

### Loss Terms Analysis

Another potential based on the Bernoulli distribution is plotted in Figure 4.7, with the network's eigenstate predictions in Figure 4.8. During the training process of the network for these approximations, the value of each term of the loss function is saved every 100 epochs. Figure 4.9 displays the evolution of these terms throughout the training process.

\begin{table}
\begin{tabular}{||c c c|c||} \hline Potential Distribution & Nodes (\(n+1\)) & Regions (\(m\)) & Eigenvalue Error (\(\%\)) \\ \hline \hline Normal & 400 & 80 & 0.863 \\ Normal & 200 & 40 & 1.16 \\ Normal & 100 & 20 & 0.123 \\ \hline Bernoulli & 200 & 40 & 0.351 \\ Uniform & 200 & 40 & 1.92 \\ \hline \end{tabular}
\end{table}
Table 1: The relative eigenvalue errors of the PINN-approximated first eigenenergy of Hamiltonians with potentials generated according to various distributions.

Figure 4.4: The potential \(V(x)\) with \(m=40\) uniform elements, where \(V(x)\) was randomly generated in each element according to the Bernoulli distribution with probability 0.5, scaled by \(10^{6}\) and shifted up by \(10^{5}\).
The spikes in Figure 4.9 represent places in the training process where either the \(L_{\text{loc}}\) term is added to the loss function or the network's parameters are reset after successfully finding an eigenstate and preparing to search for another. We highlight that the regularisation terms are minimised almost instantly, shown by how sharp the spikes are in their plots over training time. This demonstrates that our scaling of the regularisation

\begin{table}
\begin{tabular}{||c|c||} \hline Nodes (\(n+1\)) & Eigenvalue Error (\(\%\)) \\ \hline \hline
100 & 0.783 \\
200 & 0.343 \\
300 & 0.0512 \\ \hline \end{tabular}
\end{table}
Table 2: The relative eigenvalue errors of the PINN-approximated first eigenenergy of Hamiltonian with a fixed Bernoulli distributed potential over \(m=20\) regions.

Figure 4.5: The first four PINN approximated eigenstates with the potential defined in Figure 4.4. Here, ev refers to the eigenvalue associated with the plotted eigenfunction.
terms is successfully restricting the network to a search space of solutions that satisfy these regularisation terms, before working to minimise \(L_{\text{de}}\).

### Comparison with SoftIGA

Figure 4.11 plots the PINN's eigenstate prediction compared to the SoftIGA calculation for the case with a potential plotted in 4.10. We can see the PINN's eigenfunction matching well the ones calculated by SoftIGA, and the eigenvalues each with percentage errors less than 0.5% from the SoftIGA eigenvalue (which is the same for all four of these eigenstates). However, Figure 4.12 shows a significant difference between the PINN and SoftIGA calculation of the fifth and sixth eigenstates for this potential. The PINN's eigenvalue predictions are within 2% of the SoftIGA calculations, however, the SoftIGA eigenfunctions are not localised to single regions. Instead, the eigenfunctions appear to be linear combinations of localised states, while the PINN's calculations are individual localised states. We believe that the individual localised states that make up the linear combination that SoftIGA finds, and the localised states that the PINN finds are all within the same eigenspace (have the same eigenvalue), and thus their linear combinations are also eigenstates with this eigenvalue. We believe this eigenspace is the span of all the localised states in the troughs of the potential that are one segment wide (for example, the intervals \([0,0.02]\), \([0.08,0.1]\), \([0.26,0.28]\), \([0.38,0.4]\), etc., in Figure 4.10). We also believe that the eigenstates in Figure 4.11 are all part of the eigenspace spanned

Figure 4.6: The eigenvalue prediction with respect to the epoch (training process) using the potential in Figure 4.4. The dotted lines represent the true eigenvalues.
by the localised states in the troughs of the potential that are two segments wide (\([0.16,0.2]\), \([0.48,0.52]\), \([0.78,0.82]\), \([0.94,0.98]\)). Since our PINN works by scanning through each segment to specifically detect a localised state, the network does not pick up linear combinations of separate localised states, such as the SoftIGA calculations in Figure 4.12. Thus, our network is capable of separating these eigenspaces into their localised basis elements, while SoftIGA is not, displaying an area where PINNs can outperform traditional computational methods.

### An Example in Two Dimensions

In 2D, we consider a special case where the potential can be decomposed as \(V(x,y)=V_{1}(x)+V_{2}(y)\) in a rectangular domain \(\Omega=[0,1]^{2}\). The Schrodinger operator can be rewritten as

\[-\Delta+V(x,y)=(-\partial_{xx}+V_{1}(x))+(-\partial_{yy}+V_{2}(y)),\]

which allows the decomposition in the weak solution, leading to an equivalence of solving two 1D eigenvalue problems. In particular, in the finite and isogeometric element setting with tensor-product meshes, the 2D matrix eigenvalue problem can be also decomposed into tensor-product structures, leading to 1D solvers. In this setting, a base function in 2D can be rewritten as a product, i.e., \(\phi(x,y)=\phi_{1}(x)\phi_{2}(y)\). Following the derivations in [3; 5], for the eigenvalue problem (2.8) in 2D, let \(W_{h}=\text{span}\{\phi_{1,j}(x)\phi_{2,k}(y)\}\). Then, the bilinear forms can be decomposed as

Figure 4.7: The potential \(V(x)\) with \(m=40\) uniform elements, where \(V(x)\) was randomly generated in each element according to the Bernoulli distribution with probability 0.5, scaled by \(10^{6}\) and shifted up by \(10^{5}\).
\[a(\phi_{1,k}(x)\phi_{2,k}(y),\phi_{1,j}(x)\phi_{2,j}(y)) =a_{x}(\phi_{1,k}(x),\phi_{1,j}(x))\cdot b_{y}(\phi_{2,k}(y),\phi_{ 2,j}(y)) \tag{4.18}\] \[\quad+b_{x}(\phi_{1,k}(x),\phi_{1,j}(x))\cdot a_{y}(\phi_{2,k}(y),\phi_{2,j}(y)),\] \[b(\phi_{1,k}(x)\phi_{2,k}(y),\phi_{1,j}(x)\phi_{2,j}(y)) =b_{x}(\phi_{1,k}(x),\phi_{1,j}(x))\cdot b_{y}(\phi_{2,k}(y),\phi_ {2,j}(y)),\]

where the 1D bilinear forms are defined as

\[a_{x}(\phi_{1,k}(x),\phi_{1,j}(x)) =\int_{0}^{1}\left(\frac{\mathsf{d}}{\mathsf{d}x}\phi_{1,k}(x) \frac{\mathsf{d}}{\mathsf{d}x}\phi_{1,j}(x)+V_{1}(x)\phi_{1,k}(x)\phi_{1,j}(x) \right)\mathsf{d}x, \tag{4.19}\] \[a_{y}(\phi_{2,k}(y),\phi_{2,j}(y)) =\int_{0}^{1}\left(\frac{\mathsf{d}}{\mathsf{d}y}\phi_{2,k}(y) \frac{\mathsf{d}}{\mathsf{d}y}\phi_{2,j}(y)+V_{2}(x)\phi_{2,k}(y)\phi_{2,j}(y) \right)\mathsf{d}y,\] \[b_{x}(\phi_{1,k}(x),\phi_{1,j}(x)) =\int_{0}^{1}\phi_{1,k}(x)\phi_{1,j}(x)\mathsf{d}x,\] \[b_{y}(\phi_{2,k}(y),\phi_{2,j}(y)) =\int_{0}^{1}\phi_{2,k}(y)\phi_{2,j}(y)\mathsf{d}y.\]

This is also referred to as variational separability. With this in mind, we can rewrite the matrix eigenvalue

Figure 4.8: The first four PINN approximated eigenstates with the potential defined in Figure 4.7.
problem (2.9) in 2D as a problem with 1D tensor-products:

\[(K_{x}\otimes M_{y}+M_{x}\otimes K_{y})U=E_{h}(M_{x}\otimes M_{y})U, \tag{4.20}\]

which allows

\[\big{(}(K_{x}-E_{x}M_{x})\otimes M_{y}+M_{x}\otimes(K_{y}-E_{y}M_{y})\big{)}U=0, \tag{4.21}\]

where \(E_{h}=E_{h,x}+E_{h,y}\). As the dimensions are separated, this leads to two 1D eigenvalue problems (see [2; 5] for more details).

\[\begin{array}{l}K_{x}U_{x}=E_{x}M_{x}U_{x},\\ K_{y}U_{y}=E_{y}M_{y}U_{y}.\end{array} \tag{4.22}\]

With this insight of dimension decomposition in mind, we apply the PINN in \(x\) and \(y\) dimensions to find eigenstates for the 2D problem. The eigenvalues in 2D are the sum of the 1D eigenvalues while the eigen

Figure 4.9: The values of the individual terms in the loss function as the training process evolves.
states in 2D are the tensor-products of the 1D eigenstates. This significantly reduces the computational costs. As an example, Figure 4.13 shows the surface plot of one of these two-dimensional potentials, with the PINN eigenstates displayed in Figure 4.14. When the two 1D potentials admit localised eigenstates, the eigenstates of the constructed 2D potential are also localised, as \(u(x,y)=u_{1}(x)u_{2}(y)\) is only nonzero on the small rectangular region where both \(u_{1}(x)\) and \(u_{2}(y)\) are nonzero. This is demonstrated by Figure 4.14, showing the eigenstates localised to the rectangular regions.

## 5 Concluding Remarks

Recently, there has been a growing interest in using machine learning and neural networks to solve differential equations, particularly in the study of Hamiltonians. In this paper, we extend current models by constructing a neural network capable of distinguishing eigenstates from Hamiltonians with randomly distributed potentials that lead to localisation in the states. Our approach involves incorporating a normalisation loss term to avoid trivial solutions and an orthogonality loss term to search for higher eigenstates. Our method eliminates the need for driving terms, which can hinder a network's ability to converge to an admissible solution. The novelty of our approach lies in a new localisation loss term, which encourages the network to converge to a solution that is localised to a specific region. This idea is based on the physical fact that localised eigenstates occur in randomly distributed potentials. We utilise the network's initial at

Figure 4.10: The potential \(V(x)\) with \(m=40\) uniform elements using \(n=200\) nodes, where \(V(x)\) was randomly generated in each element according to the Bernoulli distribution with probability 0.5, scaled by \(10^{6}\) and shifted up by \(10^{5}\).
tempt at convergence to identify regions of localised states before adding in this term to achieve the desired approximation. We demonstrate the effectiveness of our approach in discovering eigenstates for potentials generated randomly according to different distributions, with the highlight being the Bernoulli distribution, which leads to eigenstates with identical eigenvalue. Our model successfully predicts these eigenstates with good accuracy, overcoming a challenge faced by current models.

There are several potential avenues for future work. Some potential future directions in this research area include: (1) exploring the use of more complex loss functions that can capture additional physical constraints, such as symmetries or conservation laws; (2) investigating the use of different network architectures, such as convolutional neural networks or attention-based models, to improve accuracy and scalability; (3) generalising the approach to higher-dimensional systems or non-linear differential equations; and (4) investigating the performance of the approach with larger and more complex systems, including systems with

Figure 4.11: The first four PINN predicted eigenstates compared to the SoftIGA calculated states with the potential defined in Figure 4.10. Here, ev refers to the eigenvalue calculated by the PINN for each state, with the SoftIGA calculated eigenvalue of 103679.96 for each of these four states.
a larger number of particles or more complex interactions. Overall, the use of ANNs to solve differential equations and approximate eigenstates of Hamiltonians has shown to be a promising approach, and there are many potential directions for future research in this area.
