Machine-learning the spectral function of a hole in a quantum antiferromagnet

The $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model is described by the following Hamiltonian Yin et al. (1998):

in the standard notation of the constrained fermionic operators: $\tilde{c}^{\dagger}_{i\sigma}$ creates an electron with the spin index $\sigma$ (either $\uparrow$ or $\downarrow$) at site $i$—with the constraint of no double occupancy at any site—and $\tilde{c}_{i,\sigma}$ annihilates it. The spin operators $\mathbf{S}_{i}$ expressed in the matrix form are given by $\left(\mathbf{S}_{i}\right)_{\sigma\sigma^{\prime}}=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}{\tilde{c}^{\dagger}_{i\sigma}\hat{\tau}_{\sigma\sigma^{\prime}}\tilde{c}_{i\sigma^{\prime}}}$ where $\hat{\tau}=(\hat{\tau}^{x},\hat{\tau}^{y},\hat{\tau}^{z})$ are the $2\times 2$ Pauli matrices. The angle brackets denote the first ($\langle i,j\rangle_{1}$), second ($\langle i,j\rangle_{2}$), and third ($\langle i,j\rangle_{3}$) neighbor sites, respectively. Thus, the $J$ term describes the Heisenberg interaction between nearest-neighboring quantum spins; the $t$, $t^{\prime}$, and $t^{\prime\prime}$ terms describe the electron hopping to nearest, second nearest, and third nearest sites, respectively.

The angle-resolved spectral function of a doped hole with momentum $\mathbf{k}$ and energy $\omega$ is given by

with the retarded Green’s function of the single hole being

where $E_{0}$ and $|\Psi_{0}\rangle$ are the ground-state energy and wave function of the undoped system, respectively, thus $H|\Psi_{0}\rangle=E_{0}|\Psi_{0}\rangle$. Or equivalently,

where $|\nu\rangle$ is an eigenstate of $H$ with one less electron and $E_{\nu}$ is the corresponding eigen-energy satisfying $H|\nu\rangle=E_{\nu}|\nu\rangle$, as the Dirac delta function $\delta(\omega)$ is related to a Lorentzian by $\delta(\omega)=\lim_{\eta\to 0^{+}}\frac{1}{\pi}\frac{\eta}{\omega^{2}+\eta^{2}}=\lim_{\eta\to 0^{+}}-\frac{1}{\pi}\mathrm{Im}\,\frac{1}{\omega+i\eta}$.

The angle-unresolved spectral function is given by

which is also called the density of states (DOS). To obtain the DOS in the normal procedure of theoretical calculations, one needs to have first calculated out $A(\mathbf{k},\omega)$ using Eq. (4) for a dense mesh of $\mathbf{k}$ points, and then sum the results over $\mathbf{k}$ using Eq. (5). This implies that if DOS can be accurately predicted from known DOS data, a significant speedup in evaluating DOS can be achieved, e.g., a four-orders-of-magnitude speedup compared with the normal procedure using a $100\times 100$ $\mathbf{k}$-mesh. More interestingly, we will explore whether the model Hamiltonian parameters can be accurately predicted by machine learning the DOS $A(\omega)$, which is relevant to x-ray photoemission (XPS), or $A(\mathbf{k},\omega)$ with a fixed $\mathbf{k}$, which is relevant to laser-based ARPES where the $\mathbf{k}$ points are most accessible near the zone center $\mathbf{k}=0$.

To obtain the dataset for use in our ML approach, we use the self-consistent born approximation (SCBA) to calculate Green’s function of a hole in the $t$-$t^{\prime}$-$t^{\prime\prime}$-$J$ model Schmitt-Rink et al. (1988); Marsiglio et al. (1991); Martinez and Horsch (1991); Liu and Manousakis (1991, 1992); Yin and Gong (1997, 1998); Manousakis (2007a, b); Valla et al. (2007) (see Appendix A for details). This approximation produces quantitatively accurate results for the hole Green’s function compared with exact diagonalization on small systems Leung and Gooding (1995) and Monte Carlo simulations Diamantis and Manousakis (2021).

We note that although the Hamiltonian has four parameters $(t,t^{\prime},t^{\prime\prime},J)$, all the data can be scaled with respect to $t$, e.g.

where $a$ is an arbitrary positive real number. Setting $t$ as the energy unit ($t=1$) reduces the ML complexity by one dimension, which is of significant advantage in high-throughput computation and big data management. Thus, the Green’s functions are generated in a grid of $t^{\prime}\in[-0.5,0.5]$, $t^{\prime\prime}\in[-0.5,0.5]$ and $J\in[0.2,1.0]$, with each parameter sampled on a 51-point uniform grid.

For each combination of $t^{\prime}$, $t^{\prime\prime}$, and $J$, the calculation of the Green’s function $G(\mathbf{k},\omega)$ is performed by using a $128\times 128$ mesh for the $\mathbf{k}$ points, $\omega\in[-6,6]$ with the step (i.e., energy resolution) being 0.01, and $\eta=0.01$. Then, $\eta=0.1$ is used to broaden the resulting spiky DOS and a uniform grid of 301 $\omega$ points is used to sample the DOS. Therefore, our dataset for the DOS consists of $51^{3}=132,651$ pairs $\left(\textbf{x}^{(i)},\textbf{y}^{(i)}\right)$, with $\mathbf{x}^{(i)}=(t^{\prime},t^{\prime\prime},J)^{(i)}$ being the 3-dimensional vector representation of the $i$th model parameter set and $\mathbf{y}^{(i)}=(A(\omega_{1}),A(\omega_{2}),\dots,A(\omega_{301}))^{(i)}$ the corresponding 301-dimensional vector representation of the DOS. For the forward problem, $\mathbf{x}^{(i)}$ are the input and $\mathbf{y}^{(i)}$ are the output. For the inverse problem, the definitions of $\textbf{x}^{(i)}$ and $\mathbf{y}^{(i)}$ are switched, i.e., $\mathbf{x}^{(i)}=(A(\omega_{1}),A(\omega_{2}),\dots,A(\omega_{301}))^{(i)}$ and $\mathbf{y}^{(i)}=(t^{\prime},t^{\prime\prime},J)^{(i)}$. We then randomly partitioned the dataset into an 80/10/10 training ($\mathbb{T}$), validation ($\mathbb{V}$), and testing $\mathcal{T}$ split. Here we use the computationally generated testing sets to demonstrate without any ambiguity that the ML methods work well for the present baseline problems.

Training a ML model consists of an optimization procedure in which a loss function encoding the difference between predicted and ground-truth outputs is minimized on a training set. In addition, a set of hyperparameters of the ML model is tuned during cross-validation to achieve high accuracy on the validation set. Hyperparameters are untrained parameters that include, but are not limited to, training time, network architecture and activation functions. Ultimately, final results are presented on the testing set in order to provide an unbiased estimate of model performance. Here we use the total mean squared error (MSE) as the loss function. Given the training set $\mathbb{T}=\left\{\left(\mathbf{x}^{(i)},\mathbf{y}^{(i)}\right)\right\}$ of size $|\mathbb{T}|$, i.e., $i=1,2,3,...,|\mathbb{T}|$, for an $n$-dimensional input vector $\mathbf{x}^{(i)}$, the corresponding ground-truth output is a $m$-dimensional vector $\mathbf{y}^{(i)}$; if the ML model predicts $\hat{\mathbf{y}}^{(i)}$, then the individual MSE for that training example is given by

and the total MSE score of the ML method is given by

We now introduce the two ML methods used in this work: K-nearest neighbors (KNN) and feed-forward neural network (FFNN).

To analyze the quality of our dataset and visualize the potential of applying ML algorithms to the data, we performed the following principal component analysis (PCA) Wold et al. (1987) (see Appendix B for details).

The inverse problem, which involves predicting the model’s parameters from observable quantities, has important experimental implications. Since ARPES experiments produce the spectral function and the DOS, the final goal of inverse modeling would be to predict the Hamiltonian parameters from this available experimental data. In order to make our DOS dataset more experimentally relevant, we shifted every DOS with respect to the top of the quasiparticle valence band [see Fig. 1(b)]. This better mimics experimental data, where absolute energies are not measured, but are instead found relative to the Fermi level [see Fig. 1(a)]. We also limited the dataset to include only examples with $t^{\prime}<0$ and $t^{\prime\prime}>0$, i.e, the hole doped case (the case of $t^{\prime}>0$ and $t^{\prime\prime}>0$ corresponds to electron doping). As different DOS in our dataset requires shifting of different amounts, this demands an expansion of our energy window. As a result, we now use a 354-point linear grid to sample the shifted DOS. The input is $\mathbf{x}^{(i)}=(A(\omega_{1}),A(\omega_{2}),\dots,A(\omega_{354}))^{(i)}$ and the output is $\mathbf{y}^{(i)}=(t^{\prime},t^{\prime\prime},J)^{(i)}$, where $i=1,2,3,\dots,\sim 51^{3}/4$.

In our above analysis, we have produced and analyzed DOS with $t$ being the energy unit, i.e., $t=1$. However, ARPES experiments produce DOS that are measured in terms of absolute energy. We thus proceed to analyze the feasibility of obtaining ground truth $t$ (referred to as $t_{\mathrm{truth}}$) from more experimentally realistic DOS. To this end, we examine the following simulation: We start with an SCBA-generated DOS with $t=1$ from our dataset, shifted with respect to the top of the valence band. We then scale the DOS by using $A(\omega)\to A(\omega\,t_{\mathrm{truth}})/t_{\mathrm{truth}}$, thus producing an “experimental” DOS in units of absolute energy. The task is to find $t_{\mathrm{truth}}$ from this “experimental” DOS while pretending that we do not know this $t_{\mathrm{truth}}$.

We propose the following algorithm for this task: We add to the ML methods presented in Section III.3 an outermost loop over various guesses of $t_{\mathrm{truth}}$. Specifically, for each $t_{\mathrm{guess}}$, we rescale the experimental DOS by using $A(\omega)\to A(\omega/t_{\mathrm{guess}})\,t_{\mathrm{guess}}$, producing the DOS with $t_{\mathrm{guess}}$ being the energy unit; thus, we arrive at the same inverse problem of predicting $t^{\prime}$, $t^{\prime\prime}$, and $J$ with $t=1$ studied in Section III.3. Then, we resample the rescale DOS with the the same 354 point energy grid as for the previous inverse problem, using cubic spline interpolation. After that, we run the trained neural network to produce $t^{\prime}$, $t^{\prime\prime}$, and $J$ from the rescaled DOS and calculate the MSE.

The results of this procedure are shown in Fig. 7. We find that this algorithm is able to predict $t_{truth}$ very accurately, as seen by the steep drop in the mean squared error when $t_{\mathrm{guess}}=t_{\mathrm{truth}}$. Adding such one outermost loop takes advantage of the scalability of the spectral functions [Eq. (6)] and reduces the dimensionality of the model parameter space from 4 to 3, a significant improvement in coping with the curse of dimensionality in big-data ML research.