# Efficient core-excited state orbital perspective on calculating X-ray
absorption transitions in determinant framework

Subhayan Roychoudhury<EMAIL_ADDRESS>The Molecular Foundry, Lawrence Berkeley
National Laboratory, Berkeley CA 94720, USA  David Prendergast
<EMAIL_ADDRESS>The Molecular Foundry, Lawrence Berkeley National
Laboratory, Berkeley CA 94720, USA

###### Abstract

X-ray absorption spectroscopy (XAS) is an explicit probe of the unoccupied
electronic structure of materials and an invaluable tool for fingerprinting
various electronic properties and phenomena. Computational methods capable of
simulating and analysing such spectra are therefore in high demand for
complementing the experimental results and for extracting valuable insights
therefrom. In particular, a recently proposed first-principles approach titled
Many-Body XAS (MBXAS), which approximates the final (initial) state as a
Slater determinant constructed from Kohn-Sham (KS) orbitals optimized in
absence (presence) of the relevant core-hole has shown promising prospects in
evaluating the transition amplitudes. In this article, we show that the MBXAS
approach can be rederived using a transition operator expressed entirely in
the basis of core-excited state KS orbitals and that this reformulation offers
substantial practical and conceptual advantages. In addition to circumventing
previous issues of convergence with respect to the number of unoccupied
ground-state orbitals, the aforementioned representation reduces the
computational expense by rendering the calculation of such orbitals
unnecessary altogether. The reformulated approach also provides a direct
pathway for comparing the many-body approximation with the so-called single-
particle treatment and indicates the relative importance in observed XAS
intensity of the relaxation of the valence occupied subspace induced by the
core excitation. Finally, using the core-excited state basis, we define
auxiliary orbitals for x-ray absorption and demonstrate their utility in
explaining the spectral intensity by contrasting them with single-particle
approximations to the excited state.

##  I Introduction

Several major aspects of modern technology are heavily reliant on our ability
to investigate and engineer the electronic structure of functional materials.
Therefore, a considerable portion of the advancements in electronics and
computing Coropceanu  _et al._ (2013); Messmer (1985); Hasnip  _et al._
(2014), energy storage Eng  _et al._ (2022); He  _et al._ (2019), renewable
energy harvesting Eich  _et al._ (2014); Shi and Yun (2018), drug discovery
Schaduangrat  _et al._ (2020); Arodola and Soliman (2017); Ye  _et al._
(2022), discovery of magnetic materials Sanvito  _et al._ (2017), etc. can be
attributed to the advancements in the theoretical and experimental research on
electronic structure theory. On the experimental front, electronic structure
of materials is usually explored by creating electronic excitations. Among
such methods, X-ray absorption spectroscopy (XAS) de Groot (2001); Smith and
Saykally (2017); Prendergast and Galli (2006); Li  _et al._ (2016); Freiwald
_et al._ (2004); Liu and Weng (2016); Pascal  _et al._ (2015); Wan  _et al._
(2016), which offers an element-specific probe, typically with good
resolution, has emerged as one of the most popular and effective techniques
over the last few decades. Besides providing a relatively straightforward
avenue for investigating, explicitly, the electronic structure (mostly of the
unoccupied subspace) of materials, the XAS process is widely employed for
fingerprinting various chemical attributes and phenomena of interest Qiao  _et
al._ (2014); Roychoudhury  _et al._ (2021a); Horikawa  _et al._ (2012); Qiao
_et al._ (2015); Branci  _et al._ (2000). Within a quasielectron framework,
the XAS process amounts to excitation of the core electrons of the desired
atomic species to the unoccupied levels and measures, albeit typically through
indirect means, the cross-section of such transitions. As with most
experimental techniques, the interpretation of XAS data and the extraction of
information related to electronic structure and properties therefrom is often
facilitated by complementary theoretical and computational research
Roychoudhury  _et al._ (2021b). Development of efficient and accurate
theoretical models for simulating XAS is, therefore, an active and dynamic
field of research.

Among the various first-principles methods employed in computational
electronic structure theory, Kohn-Sham (KS) density-functional theory (DFT)
has witnessed unprecedented popularity in a myriad of applications owing to
its favorable system-size scaling, modest computational requirements and
typically reasonable accuracy. Consequently, in the field of X-ray
spectroscopy, substantial effort has been invested in developing and applying
theoretical models built upon KS DFT. For instance, the response based
approaches, such as linear-response (LR) time-dependent DFT (TDDFT) Besley and
Asmuruf (2010); Fronzoni  _et al._ (2005); Hait  _et al._ (2022); Lopata  _et
al._ (2012) or the Bethe-Salpeter equation (BSE) Shirley (1998); Laskowski and
Blaha (2010); Olovsson  _et al._ (2009), aim to evaluate the absorption cross-
section from the response of the electron density to a change in the
potential. Such response, in turn, is obtained most commonly by approximating
the excited states as a linear combination of singles excitations of the KS
ground state (GS). However, the response-based technique can prove to be
expensive, especially when it involves explicit solution of the relevant
equation (LR-TDDFT/BSE) within an all-electron calculation Laskowski and Blaha
(2010); Olovsson  _et al._ (2009). Additionally, the response-based approach,
within the commonly used adiabatic approximation can incur substantial errors
in core-level spectroscopy, due to the inadequate description of the electron-
hole (e-h) interaction Roychoudhury  _et al._ (2022); Hait  _et al._ (2022). A
viable inexpensive alternative is the constrained-occupation approach, whereby
the excited state is modeled by maintaining a hole in the core system with the
help of additional constraints. Interestingly, such a treatment can be used in
a pseudopotential-based calculation by representing the core of the excited
atom by a modified pseudopotential that encodes the effects of the core-hole
Taillefumier  _et al._ (2002); Prendergast and Galli (2006); Liang  _et al._
(2017). The approximate KS orbitals corresponding to the core-excited state
can then be obtained from a self-consistent field (SCF) calculation employing
the aforementioned core-hole pseudopotential either within the Full Core-Hole
(FCH) treatment, which minimizes the energy of the core-ionized system with a
net positive charge or within the eXcited-electron and Core-Hole (XCH)
treatment, which minimizes the energy of the neutral core-excited state.

Within the so-called single-particle approximation, one can estimate the
absorption amplitude as that of a transition of a non-interacting electron
from the core-orbital to the relevant unoccupied KS orbital. Such a simplified
calculation, however, neglects the explicit effects of the relaxation of the
occupied orbitals in response to the core-excitation, resulting often in
severe inaccuracies in the absorption intensities. To mitigate these errors,
the recently proposed Many-Body X-ray Absorption Spectroscopy (MBXAS) method
relies on two separate SCF calculations:

  1. 1.

a GS calculation generating the orbitals for the initial state, and

  2. 2.

a FCH/XCH calculation generating the orbitals for the core-excited state.

Each state is then expressed as a Slater determinant (SD) composed of the
corresponding orbitals and the transition amplitude is obtained, typically
within the dipole approximation, as a many-body matrix element of the
transition operator between the initial and the final state SDs. This
determinant-based technique, which only requires simple pseudopotential-based
KS DFT calculations and consequently retains their low computational cost and
scaling, has been shown to significantly improve the peak-intensities of the
single-particle approximation, especially for the O K-edge of transition metal
oxides, and to obtain very good agreement with the experimental spectra.

Building upon this development, in this article we reformulate the MBXAS
formalism by recasting the relevant matrix elements with respect to the Full
core-hole Basis (FB) which improves the computational efficiency, eliminates
issues of convergence and enables direct comparison with the single-particle
approximation. Additionally, we develop the concept of auxiliary orbitals for
XAS which helps with the association of physically-motivated single-particle
orbitals to the excited many-electron system.

The rest of the article is organized as follows. We introduce the reformulated
treatment of MBXAS in the next section, which is followed by a prescription,
inspired from Ref. Liang and Prendergast (2018), for achieving computational
simplification. The following subsection, where we present and discuss the
results of our calculation on some representative systems, starts with a
discussion, along with some illustrative examples, on the computational
advantages of recasting the expressions to the core-hole basis set. This is
then followed by a separation of the terms contributing to the MBXAS amplitude
into categories of different physical origins which facilitate a direct
comparison with the single-particle amplitude, as discussed in the next
subsection. Finally, we move on to the discussion of auxiliary orbitals
corresponding to X-ray excitation.

##  II Method

###  II.1 General Background

Within a quasiparticle-inspired framework, X-ray absorption results in the
neutral excitation of an electronic system typically from its ground state to
an excited final state containing a core-hole and a conduction electron. In
order to simplify the terminology, without loss of generality, in the
following, the term “core” is reserved only to denote the electron/hole
corresponding to the specific orbital emptied as a consequence of X-ray
absorption. Denoting by $a^{\dagger}_{i}$ ($a^{\dagger}_{0}$) the creation
operator corresponding to the $i$-th valence level $\ket{\phi_{i}}$ (core
level $\ket{\phi_{0}}$), the ground state (GS) of the system can be expressed
as the Slater determinant (SD)

$\displaystyle\ket{\Psi_{\rm{GS}}}=\left(\prod_{i=1}^{N}a^{\dagger}_{i}\right)a_{0}^{\dagger}\ket{\emptyset},$
(1)

where $\ket{\emptyset}$ denotes the null state. In contrast, a final state can
be constructed by populating the KS orbitals obtained from a SCF calculation
performed in absence of the relevant core electron, which, in accordance with
the standard convention, will be called a full core-hole (FCH) SCF hereafter.
Using a $\sim$ to indicate quantities corresponding to such a core-excited
SCF, the core-ionized system in its lowest energy state will be given by

$\displaystyle\ket{\Psi_{\rm{FCH}}}=\left(\prod_{\mu=1}^{N}\tilde{a}^{\dagger}_{\mu}\right)\ket{\emptyset},$
(2)

with the explicit absence of the core orbital, indexed as 0. Then, within the
singles approximation, any neutral final state $\ket{\Psi_{f}}$ in which, in
addition to the lowest $N$ valence orbitals, the conduction orbital
$\ket{\tilde{\phi}_{f}}=\tilde{a}_{f}^{\dagger}\ket{\emptyset}$ is now
occupied, can be constructed from $\ket{\Psi_{\rm{FCH}}}$ as

$\displaystyle\ket{\Psi_{f}}=\ket{\Psi_{\rm{FCH}}^{+f}}=\tilde{a}_{f}^{\dagger}\ket{\Psi_{\rm{FCH}}}.$
(3)

###  II.2 MBXAS in core-excited state basis

The process of electronic transition is governed by the many-body transition
operator $\hat{O}$ which, in the earlier works, has been represented in terms
of the GS orbitals as

$\displaystyle\hat{O}=\sum_{i,j}\braket{\phi_{i}}{\hat{o}}{\phi_{j}}a^{\dagger}_{i}a_{j}.$
(4)

with $\hat{o}$ standing for the single-particle transition operator and the
sum extending over all pairs of orbitals of the GS system.

However, as with any operator, $\hat{O}$ can, in fact, be expanded in any
complete single-particle Hilbert space representation. This prompts us to
express the transition operator with respect to the FCH state orbitals as

$\displaystyle\hat{O}$
$\displaystyle{}=\sum_{l,h}\braket{\tilde{\phi}_{l}}{\hat{o}}{\tilde{\phi}_{h}}\tilde{a}_{l}^{\dagger}\tilde{a}_{h}$
$\displaystyle{}=\sum_{l,h}\tilde{o}_{lh}\tilde{a}_{l}^{\dagger}\tilde{a}_{h},$
(5)

which, as discussed in the following sections, offers substantial advantage in
terms of computation and interpretation of the X-ray absorption spectra. By
combining Eqs. 1, 2, 3, and II.2, we obtain the expression for the transition
amplitude

$\displaystyle A_{\rm{GS}\rightarrow f}$
$\displaystyle{}=\braket{\Psi_{f}}{\hat{O}}{\Psi_{\rm{GS}}}$
$\displaystyle{}=\sum_{l,h}\tilde{o}_{lh}\braket{0}{\left(\prod_{\mu=1}^{N}\tilde{a}_{\mu}\right)\tilde{a}_{f}\tilde{a}_{l}^{\dagger}\tilde{a}_{h}\left(\prod_{i=1}^{N}a_{i}^{\dagger}\right)a_{0}^{\dagger}}{0}$
$\displaystyle{}=\sum_{l=1..N,f}\tilde{o}_{l}\braket{0}{\left(\prod_{\mu=1}^{N}\tilde{a}_{\mu}\right)\tilde{a}_{f}\tilde{a}_{l}^{\dagger}\tilde{a}_{0}\left(\prod_{i=1}^{N}a_{i}^{\dagger}\right)a_{0}^{\dagger}}{0}$
$\displaystyle{}=\sum_{l=1..N,f}\tilde{o}_{l}\braket{\Psi_{f}}{\tilde{a}_{l}^{\dagger}\tilde{a}_{0}}{\Psi_{\rm{GS}}},$
(6)

where, for the sake of simplicity, we have used the notation
$\tilde{o}_{l}=\tilde{o}_{l0}$. The third line follows from the observation
that, in the sum over $h$, all terms except for that corresponding to $h=0$
(i.e, where $h$ is essentially the core-level) vanish. The approximation
$\ket{\phi_{0}}\approx\ket{\tilde{\phi}_{0}}$ is justified due to the
extremely localized nature of the core orbital Roychoudhury  _et al._ (2022).

Now, for any $f>N$ and $l$ in the range $1..N$, we introduce the many-electron
state

$\displaystyle\ket{\Psi_{\rm{FCH}}^{+f-l+0}}=(-1)^{\gamma_{l}}\left(\tilde{a}_{0}^{\dagger}\tilde{a}_{l}\ket{\Psi_{f}}\right),$
(7)

which has an index-ordered set of occupied orbitals, where the factor
$(-1)^{\gamma_{l}}$, with $\gamma_{l}=(2N-l+1)$, originates from the
anticommutation relations of the creation and annihilation operators. More
specifically, in any expression involving the creation and annihilation
operators, an exchange of two consecutive operators generates a factor of
$-1$. Thus, given that, in order to obtain $\ket{\Psi_{\rm{FCH}}^{+f-l+0}}$
from $\ket{\Psi_{f}}$, the operators $\tilde{a}_{l}$ and
$\tilde{a}_{0}^{\dagger}$ must be moved $N+1-l$ and $N$ places, respectively,
a resultant factor of $(-1)^{2N+l-1}$ is introduced. Eq. II.2 can then be
rewritten as

$\displaystyle A_{\rm{GS}\rightarrow f}=$
$\displaystyle\sum_{l=1..N}(-1)^{\gamma_{l}}\tilde{o}_{l}\braket{\Psi_{\rm{FCH}}^{+f-l+0}}{\Psi_{\rm{GS}}}$
(8)
$\displaystyle+(-1)^{\gamma_{f}}\tilde{o}_{f}\braket{\Psi_{\rm{FCH}}^{+0}}{\Psi_{\rm{GS}}},$

where $\gamma_{f}=N$, and we have used the notation
$\ket{\Psi_{\rm{FCH}}^{+f-f+0}}=\ket{\Psi_{\rm{FCH}}^{+0}}=(-1)^{N}\tilde{a}_{0}^{\dagger}\tilde{a}_{f}\ket{\Psi_{f}}$.
Now, noting that an overlap of two SDs can be expressed as a determinant
composed of the constituent orbital overlaps, the term
$\braket{\Psi_{\rm{FCH}}^{+f-l+0}}{\Psi_{\rm{GS}}}$ can be expressed as

$\displaystyle\braket{\Psi_{\rm{FCH}}^{+f-l+0}}{\Psi_{\rm{GS}}}=\det\begin{bmatrix}\xi^{*}_{1,1}&\ldots&\xi^{*}_{1,N}\\\
\vdots&\ddots&\vdots\\\ \xi^{*}_{l-1,1}&\ldots&\xi^{*}_{l-1,N}\\\
\xi^{*}_{l+1,1}&\ldots&\xi^{*}_{l+1,N}\\\ \vdots&\ddots&\vdots\\\
\xi^{*}_{N,1}&\ldots&\xi^{*}_{N,N}\\\
\xi^{*}_{f,1}&\ldots&\xi^{*}_{f,N}\end{bmatrix},$ (9)

where $\xi_{i,j}=\braket{\phi_{j}}{\tilde{\phi}_{i}}$ denotes the overlap
between the GS and the FCH KS orbitals. Thus, with the help of Eq. 8 and 9,
the MBXAS method approximates the transition amplitude for any final state
$\ket{\Psi_{\rm{FCH}}^{+f}}$ from the results of two converged SCF
calculations of KS DFT : one for the ground-state and one for the core-ionized
state.

###  II.3 Computational simplification

Introducing the notation

$\displaystyle\alpha_{q}=\left[\xi_{q,1}\textrm{
}\xi_{q,2}\ldots\xi_{q,N}\right],$ (10)

a vector of overlaps, we have

$\displaystyle\braket{\Psi_{\rm{GS}}}{\Psi_{\rm{FCH}}^{+f-l+0}}=\alpha_{1}\wedge\ldots\wedge\alpha_{l-1}\wedge\alpha_{l+1}\wedge\ldots\wedge\alpha_{N}\wedge\alpha_{f}$
(11)

Givn that this $N$-dimensional determinant defines a finite $N$-dimensional
vector space, then, for any $f>N$, $\alpha_{f}$ must be linearly dependent on
the original $N$ vectors, $\alpha_{1}\dots\alpha_{N}$, which can be expressed
in the form of a linear equation

$\displaystyle\alpha_{f}=\sum_{k=1}^{N}\kappa_{f,k}\alpha_{k},$ (12)

with $\\{\kappa_{f,k}\\}$ denoting a set of unknown coefficients which can be
computed by expressing the linear relations of Eq. 12, for $f$ in the range
$N+1\ldots M$, in the form of a matrix equation as

$\displaystyle\begin{bmatrix}\alpha_{N+1}\\\ \vdots\\\
\alpha_{M}\end{bmatrix}=\begin{bmatrix}\kappa_{N+1,1}&\ldots&\kappa_{N+1,N}\\\
\vdots&\ddots&\vdots\\\
\kappa_{M,1}&\ldots&\kappa_{M,N}\end{bmatrix}\begin{bmatrix}\alpha_{1}\\\
\vdots\\\ \alpha_{N}\end{bmatrix}$ (13)

Denoting the first, second and third matrices in Eq. 13 as
$\rm{A^{\prime}_{mat}}$, $\rm{K_{mat}}$, and $\rm{A_{mat}}$, respectively,
$\rm{K_{mat}}$ can be calculated from the equation

$\displaystyle\rm{K_{mat}}=\rm{A^{\prime}_{mat}}.\left(\rm{A_{mat}}\right)^{-1}$
(14)

Now, using the identity $a_{i}\wedge a_{j}=(-a_{j}\wedge a_{i}).\delta_{i,j}$,
the complex-conjugate of the determinantal overlap of Eq. 9 can be rewritten
as

$\displaystyle{}\braket{\Psi_{\rm{GS}}}{\Psi_{\rm{FCH}}^{+f-l+0}}$
$\displaystyle{}=\alpha_{1}\wedge\ldots\wedge\alpha_{l-1}\wedge\alpha_{l+1}\wedge\ldots\wedge\alpha_{N}\wedge\left(\sum_{k=1}^{N}\kappa_{f,k}\alpha_{k}\right)$
$\displaystyle{}=(-1)^{N-l}\kappa_{f,l}\left(\alpha_{1}\wedge\ldots\wedge\alpha_{N}\right)$
$\displaystyle{}=(-1)^{N-l}\kappa_{f,l}\rm{A_{det}},$ (15)

where

$\displaystyle\rm{A_{det}}=\left(\alpha_{1}\wedge\ldots\wedge\alpha_{N}\right)=\braket{\Psi_{\rm{GS}}}{\Psi_{\rm{FCH}}^{+0}}.$
(16)

Then, from Eq. 8, the transition amplitude can be calculated as

$\displaystyle\braket{\Psi_{f}}{\hat{O}}{\Psi_{\rm{GS}}}$
$\displaystyle{}=(-1)^{N}\left[\tilde{o}_{f}-\sum_{l=1..N}\kappa^{*}_{f,l}\tilde{o}_{l}\right]\rm{A^{*}_{det}}.$
(17)

##  III Discussions

###  III.1 Computational Advantage

Figure 1: (a) Methanol C K-edge, (b) Acetone O K-edge and (c) CuO O K-edge
absorption spectra. The number in the bracket denotes the total number of FCH
KS orbitals (excluding the core orbitals) used in the calculation. The
“scaled” curve in each panel is obtained by multiplying the dark-green one by
a constant factor in order to match the first-peak intensity with that of the
broken red one.

It is instructive to compare the MBXAS formalism in the FB Approach as
presented above with the original Ground state Basis (GB) Approach proposed in
Ref. Liang and Prendergast (2018). Both of the approaches aim to compute the
transition amplitude from the identical expression $A_{\rm{GS}\rightarrow
f}=\braket{\Psi_{f}}{\hat{O}}{\Psi_{\rm{GS}}}$ with the state
$\ket{\Psi_{\rm{GS}}}$ and $\ket{\Psi_{f}}$ being given by Eq. 1 and Eq. 3,
respectively. Therefore, at convergence, the spectra calculated using the FB
and the GB approaches should, in principle, be the same. In practice, however,
there can be substantial differences.

In the GB treatment, the transition amplitude within the singles approximation
is simplified to

$\displaystyle A_{\rm{GS}\rightarrow f}$
$\displaystyle{}=\sum_{c>N}\braket{\Psi_{\rm{FCH}}^{+f}}{\Psi_{\rm{GS}}^{+c-0}}\braket{\Psi_{\rm{GS}}^{+c-0}}{\hat{O}}{\Psi_{\rm{GS}}}$
$\displaystyle{}=\sum_{c>N}o_{c}\braket{\Psi_{\rm{FCH}}^{+f}}{\Psi_{\rm{GS}}^{+c-0}},$
(18)

where $o_{c}=\braket{\phi_{c}}{\hat{o}}{\phi_{0}}$ and

$\displaystyle\ket{\Psi_{\rm{GS}}^{+c-0}}=(-1)^{N}a_{c}^{\dagger}a_{0}\ket{\Psi_{\rm{GS}}},$
(19)

The relevant many-body overlap is then given by

$\displaystyle\braket{\Psi_{\rm{FCH}}^{+f}}{\Psi_{\rm{GS}}^{+c-0}}=(-1)^{N}\begin{bmatrix}\xi^{*}_{1,1}&\ldots&\xi^{*}_{1,N}&\xi^{*}_{1,c}\\\
\vdots&\ddots&\vdots&\vdots\\\
\xi^{*}_{N,1}&\ldots&\xi^{*}_{N,N}&\xi^{*}_{N,c}\\\
\xi^{*}_{f,1}&\ldots&\xi^{*}_{f,N}&\xi^{*}_{f,c}\end{bmatrix},$ (20)

Note that, in principle, the sum in Eq. III.1 extends over all unoccupied GS
orbitals - bound and unbound both. Thus, if one fails to include any GS
unoccupied orbital which has a non-zero overlap with any of the valence
occupied FCH-state orbitals (i.e.
$\ket{\tilde{\phi}_{1}},\ldots,\ket{\tilde{\phi}_{N}}$), then, in principle,
$A_{\rm{GS}\rightarrow f}$ will be under-converged for any $f>N$.
Additionally, if, for a given $f>N$, $\ket{\tilde{\phi}_{f}}$ overlaps with
any missing unoccupied GS orbitals, then the corresponding
$A_{\rm{GS}\rightarrow f}$ will be under-converged. Therefore, convergence in
Eq. III.1 can be particularly difficult to achieve if one or more of these
FCH-state orbitals has appreciable overlap with the unbound unoccupied
orbitals of the GS Forestell and Marsiglio (2015).

As an aside, we may note that similar convergence issues should plague any
computational method for simulating XAS that is rooted in a GS orbital
representation of the XAS transition amplitude. Prominent examples include
both the aforementioned LR-TDDFT Besley and Asmuruf (2010); Fronzoni  _et al._
(2005); Hait  _et al._ (2022); Lopata  _et al._ (2012) and BSE Shirley (1998);
Laskowski and Blaha (2010); Olovsson  _et al._ (2009) approaches, which are
both expressed within a GS orbital basis that includes unoccupied orbitals.

As illustrative examples, in Fig. 1, we present comparative studies of the GB
and the FB treatments for the C K-edge spectrum of methanol ($\rm{CH_{3}OH}$),
O K-edge spectrum of acetone ($\rm{C_{3}H_{6}O}$) and the O K-edge spectrum of
cupric oxide ($\rm{CuO}$). The GB spectra are seen to be heavily dependent on
the number of unoccupied GS orbitals employed in the calculation. The fact
that the large GB spectra do not coincide with the scaled small GB
counterparts indicates that the under-converged GB spectra differ inherently
in their lineshapes.

Since the sum in Eq. II.2 is limited to a finite number $(N+1)$ of terms
anyway, this complication is circumvented entirely in the FB treatment as
evident from Fig. 1, which shows, as expected, that the FB spectra do not
depend on the number of orbitals included in the calculation, as long as all
FCH-state orbitals within the desired energy range are included. In addition
to not suffering from convergence issues, for any given energy range of the
spectrum, the FB treatment is computationally favourable over the GB treatment
since, unlike the latter, the former does not require the computation of any
unoccupied GS orbital.

For a simulation involving $N_{\rm{tot}}$ number of GS orbitals (including the
valence occupied orbitals as well as the unoccupied ones), within the GB
approach, a practical MBXAS calculation requires the execution of the
following steps separately for each of the GS and the FCH systems:

  1. 1.

a SCF calculation to obtain the electron density,

  2. 2.

a non self-consistent field (NSCF) calculation to obtain the unoccupied KS
orbitals, which

  3. 3.

for computational simplification, are often transformed into an optimal basis-
representation to facilitate faster Brillouin zone integration.

Additionally, for the GS system, we

  1. 4.

calculate the single-particle transition matrix elements $o_{c}$ for all
unoccupied GS orbitals $\ket{\phi_{c}}$.

This is followed by

  1. 5.

calculation of the $(N_{\rm{tot}}\times N_{\rm{tot}})$ overlap matrix between
the ground and the FCH state orbitals,

and finally

  1. 6.

the computation of the MBXAS amplitudes from Eq. III.1.

On the other hand, in the FB treatment, for the GS system, step (2) is not
required at all while step (3) can become substantially cheaper due to the
absence of the unoccupied orbitals. The same is true for step (5) which now
calculates the $(N\times N_{\rm{tot}})$ overlap matrix $\xi_{i,j}$ $\forall
i\leq N_{\rm{tot}}\textrm{ and }j\leq N$ and for step (6) which, now follows
Eq. 8 and thereby involves a sum over $(N+1)$ terms. Note that, in the FB
approach, step (4) is replaced by a calculation of the single-particle
transition matrix elements $\tilde{o}_{l}$ over all of the FCH-state orbitals,
occupied and unoccupied alike, and therefore, should be a little more
expensive than the corresponding GB calculation. However, typically, for a
calculation involving a large number of unoccupied orbitals, steps (2) and (3)
are, by a large margin, the most expensive ones computationally. Therefore,
assuming that, for a given system, they have comparable costs for the GS and
the FCH state, the FB approach is expected to reduce the total computational
cost approximately by half. This, in fact, can be seen from Table 1, which
shows the estimates of the computational time required for methanol C K-edge,
acetone O K-edge and CuO O K-edge calculations with a relatively large number
of unoccupied orbitals.

System |  Core hours (GB) |  Core hours (FB)
---|---|---
Methanol C K (2750) |  95 |  46
Acetone O K (3200) |  166 |  82
CuO O K (3260) |  1056 |  597
Table 1: Table showing for methanol C K-edge, acetone O K-edge and CuO O
K-edge, the estimates of core hours used in the XAS calculations using the GB
and the FB approaches. The number in the bracket denotes the total number of
FCH KS orbitals (excluding the core orbitals) used in the calculation.

###  III.2 Separation of Spectral Contributions

Note from Eq. 8, that the MBXAS amplitude $A_{\rm{GS}\rightarrow f}$ is
essentially a linear-combination, for $l=1,\ldots,N$ and $l=f$, of single-
particle transition amplitudes $\tilde{o}_{l}$ from the core orbital to the
orbital $\ket{\tilde{\phi}_{l}}$, with the factor
$\braket{\Psi_{\rm{FCH}}^{+f-l+0}}{\Psi_{\rm{GS}}}$ serving as the
coefficient. We can naturally separate these terms according to their physical
origins as follows:

  1. 1.

The $l=f$ term essentially corresponds to the expected quasiparticle
transition, from the core orbital to orbital $\ket{\tilde{\phi}_{f}}$, which
is formally unoccupied in $\ket{\Psi_{\rm{FCH}}}$, but indirectly includes
some perturbative effects of the presence of the core-hole (since it is
calculated based on the core-excited SCF).

  2. 2.

The $l=1,\ldots,N$ terms accounts for the collective response of the rest of
the electron density to the perturbation caused by the core hole by
considering single-particle contributions from each of the occupied valence
orbitals of $\ket{\Psi_{\rm{FCH}}}$.

This separation is just a natural consequence of the many-body transition
operator coupling the $N$-electron ground and core-excited states. However,
another way to understand it is to consider the equivalence of absorption
amplitude to the (stimulated) emission amplitude. We can rewrite Eq. 8 as

$\displaystyle{}A_{\rm{GS}\rightarrow f}$
$\displaystyle{}=(-1)^{\gamma_{l}}\left[\sum_{l=1,\ldots,N,f}\braket{\Psi_{\rm{FCH}}^{+f}}{\hat{O}}{\Psi_{\rm{FCH}}^{+f-l+0}}\braket{\Psi_{\rm{FCH}}^{+f-l+0}}{\Psi_{\rm{GS}}}\right].$
(21)

Now we can interpret this as a de-excitation from $\ket{\Psi_{\rm{FCH}}^{+f}}$
to $\ket{\Psi_{\rm{GS}}}$, where the $l$-th representative term can be
recognized as the amplitude
$\braket{\Psi_{\rm{FCH}}^{+f}}{\hat{O}}{\Psi_{\rm{FCH}}^{+f-l+0}}$ of a non-
selfconsistent de-excitation from the state $\ket{\Psi_{\rm{FCH}}^{+f}}$ to
the core-filled state $\ket{\Psi_{\rm{FCH}}^{+f-l+0}}$, multiplied, for the
sake of restoring self-consistency, by the projection of this non-self
consistent state onto the actual final state $\ket{\Psi_{\rm{GS}}}$. Thus,
each of the $(N+1)$ terms in the sum corresponds to de-excitation of an
electron from one of the occupied levels of $\ket{\Psi_{\rm{FCH}}^{+f}}$ to
the core level.

###  III.3 Comparison with the single-particle formalism

Armed with the numerically converged FB treatment, we can now examine the
relative importance of its various contributions. According to the expansion
shown in Eq. 8, $\tilde{o}_{l}$ is one among the $N+1$ terms
($l=1,\ldots,N,f$) contributing to the transition amplitude. This affords a
straightforward avenue for comparison of the MBXAS formalism with the more
common single-particle formalism of XAS, in which the transition amplitude
$A_{\rm{GS}\rightarrow f}$ is approximated simply by $\tilde{o}_{f}$.

Mathematically, it can be shown that this truly single-particle approximation
would be true if the valence occupied subspace of the ground and the FCH state
were identical. Then
$\braket{\Psi_{\rm{FCH}}^{+f-l+0}}{\Psi_{\rm{GS}}}=\delta_{f,l}$ and
consequently, $A_{\rm{GS}\rightarrow f}=\tilde{o}_{f}$ (possibly with an
inconsequential change in sign). We could now ask whether the degree to which
the valence occupied subspaces differ is a useful metric to indicate a larger
difference between $A_{\rm{GS}\rightarrow f}$ and $\tilde{o}_{f}$. We explore
this possibility through a detailed analysis of the first peak of the C and O
K-edge absorption spectra of an isolated methanol molecule.

Figure 2: Panels (a) and (b) show, for the methanol molecule, the C K-edge XAS
spectrum for incident polarization along the x-axis and the O K-edge XAS
spectrum for incident polarization along the y-axis, respectively (see
molecular polarizaton in Fig. 3). Panels (c) and (d) analyse, in detail, the
orbital contributions responsible for MBXAS transition amplitude of the most
intense peaks seen in panels (a) and (b), respectively. In panel (c) (panel
(d)) the green, black and blue lines correspond to
$\left|A_{\rm{GS}\rightarrow 11}\right|$, $\left|\tilde{o}_{11}\right|$ and
$\left|C_{11}^{(11)}\right|$ ($\left|A_{\rm{GS}\rightarrow 9}\right|$,
$\left|\tilde{o}_{9}\right|$ and $\left|C_{9}^{(9)}\right|$), respectively.
Figure 3: Isovalue plots of occupied KS orbitals of the FCH-states of the
methanol molecule resulting from core-ionization of specific atoms. The
strongly contributing orbitals to the XAS intensity are indexed according to
the energy ordering of their KS eigenvalues as (a) 3, (b) 5, (c) 11 (LUMO+2)
for the C K edge and (d) 2, (e) 5, (f) 9 (LUMO) for the O K edge. The core-
ionized atom is marked by a “+” sign.

As stated concisely in Eq. 8, the MBXAS amplitude is essentially a sum of the
terms

$\displaystyle
C_{l}^{(f)}=(-1)^{\gamma_{l}}\tilde{o}_{l}\braket{\Psi_{\rm{FCH}}^{+f-l+0}}{\Psi_{\rm{GS}}}\textrm{
for }l=1,\ldots,N,f,$ (22)

for a given core-excited state $\ket{\Psi_{\rm{FCH}}^{+f}}$ labeled by core-
excited (formerly unoccupied) FCH orbital index $f$. We will now assess which
$C_{l}^{(f)}$ are strong contributors to specific X-ray absorption transitions
in methanol, by reference to Figs. 2 and 3.

At the C K-edge in Fig. 2(a), the most intense XAS peak for photons polarized
along the x-axis (see Fig. 3 for the molecular orientation) occurs just below
290 eV and corresponds to transitions to $\ket{\Psi_{\rm{FCH}}^{+11}}$ for the
core-excited C. As expected, the dominant contribution to this transition is
$C_{11}^{(11)}$, defined by the $f=11$ orbital (FCH LUMO+2) that labels this
excited state. However, in addition, we see in Fig. 2(c) that there are also
significant transition amplitude contributions from $C_{3}^{(11)}$ and
$C_{5}^{(11)}$ which involve occupied orbitals 3 and 5 from this FCH SCF. As
expected, each of these orbitals (3, 5, 11) display noticeable
$p_{x}$-character on the core-excited C atom, as shown in Fig. 3(a-c).

On the other hand, at the O K-edge in Fig. 2(b), the most intense XAS peak for
photons polarized along the y-axis occurs above 536 eV and corresponds to
transitions to $\ket{\Psi_{\rm{FCH}}^{+9}}$ for the core-excited O. As before,
the dominant contribution to this transition is $C_{9}^{(9)}$, defined by the
$f=9$ orbital (FCH LUMO) that labels this excited state. However, in addition,
we see in Fig. 2(d) that there are also significant transition amplitude
contributions from $C_{2}^{(9)}$ and $C_{5}^{(9)}$ which involve occupied
orbitals 2 and 5 from this FCH SCF. As expected, each of these orbitals (2, 5,
9) display noticeable $p_{y}$-character on the core-excited O atom, as shown
in Fig. 3(d-f).

The computed absolute values of the valence-overlap between
$\ket{\Psi_{\rm{FCH}}}$ and $\ket{\Psi_{\rm{GS}}}$ are 0.94 and 0.83 for the C
and the O K-edge of methanol, respectively. Since the excitation is expected
to induce only negligible change in the tightly-bound core orbital, the
aforementioned valence-overlap can be used as a good approximation for
$\braket{\Psi_{\rm{FCH}}^{+0}}{\Psi_{\rm{GS}}}$. As expected, owing to a
higher overlap, the difference between $|C_{11}^{(11)}|$ and
$|\tilde{o}_{11}|$ is smaller at the C K edge than that between
$|C_{9}^{(9)}|$ and $|\tilde{o}_{9}|$ at the O K edge. However, the relative
difference between the absolute values of the MBXAS and the single-particle
amplitude, i.e., $\frac{|A_{\rm{GS}\rightarrow
f}|-|\tilde{o}_{f}|}{|\tilde{o}_{f}|}$ (this is the relative difference
between the green and black horizontal lines in Fig. 2 (c) and (d)) is much
higher for C (0.15) than it is for O (0.07), demonstrating that a lower
overlap between the valence and the FCH occupied subspaces does not
necessarily imply a larger difference between $|A_{\rm{GS}\rightarrow f}|$ and
$|\tilde{o}_{f}|$. In general, the formally more accurate option is clearly to
include all additional valence occupied contributions, although it is clear
that the XAS is somewhat dominated by its single-particle contributions in
this case. We will explore a counter-example below.

###  III.4 Auxiliary Orbital

Figure 4: Panel (a) shows the O K-edge XAS of CrO2 simulated using the single-
particle (S-P) and MBXAS methods. Panel (b) shows separately the up and down
spin contributions to spectra of both methods. The small orange (blue) circle
corresponds to the first peak in the up-spin spectrum and is analysed using
the isovalue plot of panel (c) (panel (d)) which shows, at the $\Gamma$ point,
the orbital $\tilde{\phi}_{f}(\mathbf{r})$
($\tilde{\phi}^{\rm{aux}}_{f}(\mathbf{r})$) that contributes most strongly to
this peak. In this case, $f$ corresponds to the lowest unoccupied level of the
FCH state. The core-ionized oxygen atom is marked by a “+” symbol in both
lower panels.

It is important to note that the excitation described by MBXAS is an
inherently many-body process with no physical single-particle analogue.
However, the transition amplitude given by Eq. 8 inspires the formulation of a
hypothetical single-particle excitation of a core electron to an auxiliary
orbital given by

$\displaystyle\ket{\tilde{\phi}_{\rm{aux}}^{f}}=\sum_{l=1,\ldots,N,f}(-1)^{\gamma_{l}}\braket{\Psi_{\rm{GS}}}{\Psi_{\rm{FCH}}^{+f-l+0}}\ket{\tilde{\phi}_{l}}$
(23)

such that the amplitude of this excitation equals that of the actual many-body
one, i.e.,

$\displaystyle\tilde{o}^{\rm{aux}}_{f}=\braket{\tilde{\phi}_{\rm{aux}}^{f}}{\hat{o}}{\phi_{0}}=A_{\rm{GS}\rightarrow
f},$ (24)

as can be seen by comparison with Eq. II.2. The utility of the auxiliary
orbital can be more easily appreciated by looking at the O K-edge spectrum of
CrO2 which was discussed in detail in Ref Liang and Prendergast (2018).

As shown in Fig. 4(b), for the up-spin channel, the intensity of the O K-edge
first peak at 529 eV is severely underestimated by the single-particle
approximation. This underestimation can be understood by examining, at the
$\Gamma$ point, the corresponding excited single-particle orbital in the
isovalue plots. The filled FCH orbital $f$, in this case, sits at the
conduction band minimum (CBM), or first unoccupied orbital (analogous to LUMO
at the $\Gamma$-point) of the FCH state. We see from Fig. 4(c) that
$\tilde{\phi}_{\rm{LUMO}}(\mathbf{r})$ has apparently no visible presence on
the core-ionized O atom, which explains the almost vanishing intensity of this
first peak in the single-particle spectrum. On the other hand, the associated
auxiliary orbital $\tilde{\phi}^{\rm{aux}}_{\rm{LUMO}}(\mathbf{r})$ shown in
Fig. 4(d) exhibits substantial O-p character on the core-ionized atom, due to
inclusion of contributions from occupied valence orbitals in Eq. 23, and hence
gives rise to a drastically higher (and more accurate) MBXAS intensity.

Fig. 5(a) shows, as a function of the FCH orbital $\tilde{\phi}_{l}$, the
absolute value of $C_{l}^{(\mathrm{LUMO})}$ (see Eq. 22), which is essentially
the contribution to the first-peak MBXAS amplitude (Eq. 8) of the single-
particle transition associated with the orbital $\tilde{\phi}_{l}$. On the
other hand, Fig. 5(b) shows the absolute value of
$\mathbb{C}_{l}^{(\rm{LUMO})}=(-1)^{\gamma_{l}}\braket{\Psi_{\rm{GS}}}{\Psi_{\rm{FCH}}^{+{\rm{LUMO}}-l+0}}$,
which, as seen from Eq. 23, serves as the contribution of $\tilde{\phi}_{l}$
to the auxiliary orbital $\tilde{\phi}^{\rm{aux}}_{\rm{LUMO}}$. As expected,
due to the factor of $\tilde{o}_{l}$, these two contributions do not follow
the same trend. Panel (a) of Fig. 5 clearly shows that the largest
contribution to the MBXAS amplitude arises from the 148-th FCH orbital
$\tilde{\phi}_{148}$, situated deep in the valence subspace. Predictably,
$\tilde{\phi}_{148}$ exhibits strong $p$ character on the excited oxygen atom,
as seen from panel (e) of Fig. 6, which shows the isovalue plots of FCH
orbitals with appreciable contribution in
$\tilde{\phi}^{\rm{aux}}_{\rm{LUMO}}$. Conversely, as evident from Fig. 5(b),
the largest contribution to the auxiliary orbital comes from the FCH LUMO,
which, as already seen from Fig. 4(c), has negligible $p$-character on the
core-excited oxygen atom. For the auxiliary orbital
$\tilde{\phi}_{\rm{LUMO}}^{\rm{aux}}$, the noticeable $p$ character on the
excited atom can then be attributed to the cumulative contributions from the
O-p character of various valence orbitals having appreciable
$\mathbb{C}_{l}^{(\rm{LUMO})}$.

Figure 5: For O K-edge electronic transition at the $\Gamma$ point of CrO2,
panels (a) and (b) show, as a function of the FCH KS orbital indexed by $l$,
the absolute values $\left|C_{l}^{(337)}\right|$ (see Eq. 22)and
$\left|\mathbb{C}_{l}^{(337)}\right|$, respectively, where the 337-th FCH
orbital is the lowest unoccupied one.

Note, from Eq. III.1 that, auxiliary orbitals can also be defined in the GB
approach as

$\displaystyle\ket{\tilde{\phi}_{\rm{aux}}^{f}}=\sum_{c>N}\braket{\Psi_{\rm{FCH}}^{+f}}{\Psi_{\rm{GS}}^{+c-0}}\ket{\phi_{c}}.$
(25)

However, due to the convergence issues discussed above, this might require the
evaluation and manipulation of a large number of GS unoccupied orbitals,
rendering the calculation formidably expensive.

Figure 6: Isovalue plots of $\Gamma$-point FCH KS orbitals contributing most
strongly to the auxiliary orbital $\tilde{\phi}^{\rm{LUMO}}_{f}(\mathbf{r})$
shown in Fig. 4(d). Panels (a) through (f) indicate decreasing contributions
$\left|\mathbb{C}_{l}^{\rm{LUMO}}\right|$ starting from the largest: (a)
orbital 337 (LUMO), (b) 335 (HOMO-1), (c) 336 (HOMO), (d) 282, (e) 148 and (f)
260. The corresponding contributions can be found in Fig. 5(b).

###  III.5 The $f^{(n)}$ contribution

At this point, it is worth mentioning briefly the contribution of the so-
called $f^{(n)}$ term Liang and Prendergast (2018) within the MBXAS formalism,
whereby the final state

$\displaystyle\ket{\Psi_{f}^{(n)}}$
$\displaystyle{}=\ket{\Psi_{\rm{FCH}}^{+f_{n-1}+\ldots+f_{1}+f-h_{n-1}-\ldots-
h_{1}}}$
$\displaystyle{}=(-1)^{\Lambda}\tilde{a}^{\dagger}_{f_{n-1}}\ldots\tilde{a}^{\dagger}_{f_{1}}\tilde{a}^{\dagger}_{f}\tilde{a}_{h_{n-1}}\ldots\tilde{a}_{h_{1}}\ket{\Psi_{\mathrm{FCH}}},$
(26)

with $\Lambda=(nN-h_{1}-\ldots-h_{n-1})$, is obtained by generating, in
addition to the core-excitation, $n-1$ number of electron-hole pairs outside
the core-region, with the $m$-th conduction electron (valence hole) being
indexed by $f_{m}$ ($h_{m}$). Without loss of generality, the orbital-indices
are made to follow the order $f_{n-1}>\ldots>f_{1}>f>h_{n-1}>\ldots>h_{1}$.
The core-excited state shown in Eq. 3, which contains no valence-hole, can now
be recognized as $\ket{\Psi_{f}^{(1)}}$.

Within the FB approach, the transition amplitude associated with
$\ket{\Psi_{f}^{(n)}}$ can then be expressed as

$\displaystyle{}A_{\rm{GS}\rightarrow f}^{(n)}$
$\displaystyle{}=\braket{\Psi_{f}^{(n)}}{\sum_{l}\tilde{o}_{l}\tilde{a}_{l}^{\dagger}\tilde{a}_{0}}{\Psi_{\mathrm{GS}}}$
$\displaystyle{}=\sum_{l}(-1)^{\gamma_{l}^{(n)}}\tilde{o}_{l}\braket{\Psi_{\mathrm{FCH}}^{+f_{n-1}+\ldots+f_{1}+f-h_{n-1}-\ldots-
h_{1}-l+0}}{\Psi_{\mathrm{GS}}},$ (27)

where

$\displaystyle\ket{\Psi_{\mathrm{FCH}}^{+f_{n-1}+\ldots+f_{1}+f-h_{n-1}-\ldots-
h_{1}-l+0}}=(-1)^{\gamma_{l}^{(n)}}\tilde{a}_{l}\tilde{a}^{\dagger}_{0}\ket{\Psi_{f}^{(n)}},$
(28)

with the term $\gamma_{l}^{(n)}$ depending on the relative position of $l$
with respect to $h_{1}\ldots h_{n-1}$. Eq. III.5 and 28 reveal that, within
the FB approach, the GS unoccupied orbitals are not needed to compute the
transition amplitude $A^{(n)}_{\rm{GS}\rightarrow f}$ for any value of $n$.
Similar to the special case of $f^{(1)}$ (see Eq. 23), a general transition
arising from the $f^{(n)}$ term can be associated with the auxiliary orbital

$\displaystyle{}\ket{\tilde{\phi}^{f^{(n)}}_{\rm{aux}}}$ (29)
$\displaystyle{}=\sum_{l}(-1)^{\gamma_{l}^{(n)}}\braket{\Psi_{\rm{GS}}}{\Psi_{\rm{FCH}}^{+f_{n-1}+\ldots+f_{1}+f-h_{n-1}-\ldots-
h_{1}-l+0}}\ket{\tilde{\phi}_{l}}$ (30)

As originally defined, we reiterate that within the MBXAS formalism, the
contribution of $f^{(n)}$ for $n>1$ ignores interactions between the
subsequent electron-hole pairs and, therefore, is really only accurate for
studies of metallic systems. It is especially inaccurate for
insulating/semiconducting systems which should exhibit significant exciton
binding for additional valence electron-hole pairs.

##  IV Conclusion

In summary, in this article, we have presented an efficient representation of
the determinant-based many-body formalism for simulating x-ray absorption
spectroscopy, called MBXAS. We have shown that by recasting the relevant
equations in terms of the FB orbitals, one can circumvent the need for
calculating the unoccupied orbitals of the ground-state, leading to
appreciable reduction in the computational expenses. Additionally, the
reformulated expression for the transition amplitude does not suffer from
convergence issues since it is now limited to a sum of a finite number of
terms. Within the FB approach, the amplitude is expanded as a linear
combination of independent-electron transitions from various core-excited
state orbitals. Therefore, such a treatment facilitates a direct comparison of
the MBXAS formalism with the so-called single-particle approximation. As we
show analytically and with a representative example, the aforementioned
comparison motivates the association of a given MBXAS transition to a purely
single-particle electronic excitation from the core level to a hypothetical
auxiliary orbital. A visual representation of this auxiliary orbital offers
intuitive justification for the relevant absorption intensity.

##  V Acknowledgement

This work was carried out at the Molecular Foundry (TMF), Lawrence Berkeley
National Laboratory (LBNL). Financial support for this work was provided by
the Office of Science, Office of Basic Energy Sciences, of the U.S. Department
of Energy under Contract No. DE-AC02-05CH11231. Computational works were
performed using the supercomputing resources of the National Energy Research
Scientific Computing Center (NERSC) and the TMF clusters managed by the High
Performance Computing Services Group, at LBNL.

## References

  * Coropceanu  _et al._ (2013) V. Coropceanu, H. Li, P. Winget, L. Zhu, and J.-L. Brédas, Annual Review of Materials Research 43, 63 (2013), https://doi.org/10.1146/annurev-matsci-071312-121630 .
  * Messmer (1985) R. P. Messmer, Journal of Non-Crystalline Solids 75, 285 (1985), proceedings of the international conference on the theory of the structures of non-crystalline solids.
  * Hasnip  _et al._ (2014) P. J. Hasnip, K. Refson, M. I. J. Probert, J. R. Yates, S. J. Clark, and C. J. Pickard, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20130270 (2014).
  * Eng  _et al._ (2022) A. Y. S. Eng, C. B. Soni, Y. Lum, E. Khoo, Z. Yao, S. K. Vineeth, V. Kumar, J. Lu, C. S. Johnson, C. Wolverton, and Z. W. Seh, Science Advances 8, eabm2422 (2022), https://www.science.org/doi/pdf/10.1126/sciadv.abm2422 .
  * He  _et al._ (2019) Q. He, B. Yu, Z. Li, and Y. Zhao, ENERGY & ENVIRONMENTAL MATERIALS 2, 264 (2019), https://onlinelibrary.wiley.com/doi/pdf/10.1002/eem2.12056 .
  * Eich  _et al._ (2014) F. G. Eich, M. Di Ventra, and G. Vignale, Phys. Rev. Lett. 112, 196401 (2014).
  * Shi and Yun (2018) J. Shi and S. Yun, First-principles dft calculations for perovskite solar cells, in  _Counter Electrodes for Dye‐sensitized and Perovskite Solar Cells_ (John Wiley & Sons, Ltd, 2018) Chap. 19, pp. 487–509.
  * Schaduangrat  _et al._ (2020) N. Schaduangrat, S. Lampa, S. Simeon, M. P. Gleeson, O. Spjuth, and C. Nantasenamat, Journal of Cheminformatics 12, 9 (2020).
  * Arodola and Soliman (2017) O. A. Arodola and M. E. Soliman, Drug design, development and therapy 11, 2551 (2017), 28919707[pmid].
  * Ye  _et al._ (2022) N. Ye, Z. Yang, and Y. Liu, Drug discovery today 27, 1411 (2022), 34954327[pmid].
  * Sanvito  _et al._ (2017) S. Sanvito, C. Oses, J. Xue, A. Tiwari, M. Zic, T. Archer, P. Tozman, M. Venkatesan, M. Coey, and S. Curtarolo, Science Advances 3, e1602241 (2017), https://www.science.org/doi/pdf/10.1126/sciadv.1602241 .
  * de Groot (2001) F. de Groot, Chemical Reviews 101, 1779 (2001), pMID: 11709999, https://doi.org/10.1021/cr9900681 .
  * Smith and Saykally (2017) J. W. Smith and R. J. Saykally, Chemical Reviews 117, 13909 (2017), pMID: 29125751, https://doi.org/10.1021/acs.chemrev.7b00213 .
  * Prendergast and Galli (2006) D. Prendergast and G. Galli, Phys. Rev. Lett. 96, 215502 (2006).
  * Li  _et al._ (2016) Q. Li, R. Qiao, L. A. Wray, J. Chen, Z. Zhuo, Y. Chen, S. Yan, F. Pan, Z. Hussain, and W. Yang, Journal of Physics D: Applied Physics 49, 413003 (2016).
  * Freiwald  _et al._ (2004) M. Freiwald, S. Cramm, W. Eberhardt, and S. Eisebitt, Journal of Electron Spectroscopy and Related Phenomena 137-140, 413 (2004), iCESS-9 Proceedings of the 9th International Conference on Electronic Spectroscopy and Structure.
  * Liu and Weng (2016) X. Liu and T.-C. Weng, MRS Bulletin 41, 466–472 (2016).
  * Pascal  _et al._ (2015) T. A. Pascal, C. D. Pemmaraju, and D. Prendergast, Phys. Chem. Chem. Phys. 17, 7743 (2015).
  * Wan  _et al._ (2016) L. F. Wan, J. Wright, B. R. Perdue, T. T. Fister, S. Kim, C. A. Apblett, and D. Prendergast, Phys. Chem. Chem. Phys. 18, 17326 (2016).
  * Qiao  _et al._ (2014) R. Qiao, I. T. Lucas, A. Karim, J. Syzdek, X. Liu, W. Chen, K. Persson, R. Kostecki, and W. Yang, Advanced Materials Interfaces 1, 1300115 (2014), https://onlinelibrary.wiley.com/doi/pdf/10.1002/admi.201300115 .
  * Roychoudhury  _et al._ (2021a) S. Roychoudhury, R. Qiao, Z. Zhuo, Q. Li, Y. Lyu, J.-H. Kim, J. Liu, E. Lee, B. J. Polzin, J. Guo, S. Yan, Y. Hu, H. Li, D. Prendergast, and W. Yang, ENERGY & ENVIRONMENTAL MATERIALS 4, 246 (2021a), https://onlinelibrary.wiley.com/doi/pdf/10.1002/eem2.12119 .
  * Horikawa  _et al._ (2012) Y. Horikawa, H. Arai, T. Tokushima, and S. Shin, Chemical Physics Letters 522, 33 (2012).
  * Qiao  _et al._ (2015) R. Qiao, L. A. Wray, J.-H. Kim, N. P. W. Pieczonka, S. J. Harris, and W. Yang, The Journal of Physical Chemistry C 119, 27228 (2015), https://doi.org/10.1021/acs.jpcc.5b07479 .
  * Branci  _et al._ (2000) C. Branci, M. Womes, P. Lippens, J. Olivier-Fourcade, and J. Jumas, Journal of Solid State Chemistry 150, 363 (2000).
  * Roychoudhury  _et al._ (2021b) S. Roychoudhury, Z. Zhuo, R. Qiao, L. Wan, Y. Liang, F. Pan, Y.-d. Chuang, D. Prendergast, and W. Yang, ACS Applied Materials & Interfaces 13, 45488 (2021b).
  * Besley and Asmuruf (2010) N. A. Besley and F. A. Asmuruf, Phys. Chem. Chem. Phys. 12, 12024 (2010).
  * Fronzoni  _et al._ (2005) G. Fronzoni, R. De Francesco, and M. Stener, The Journal of Physical Chemistry B 109, 10332 (2005), pMID: 16852252, https://doi.org/10.1021/jp050755y .
  * Hait  _et al._ (2022) D. Hait, K. J. Oosterbaan, K. Carter-Fenk, and M. Head-Gordon, The Journal of Chemical Physics 156, 201104 (2022), https://doi.org/10.1063/5.0092987 .
  * Lopata  _et al._ (2012) K. Lopata, B. E. Van Kuiken, M. Khalil, and N. Govind, Journal of Chemical Theory and Computation 8, 3284 (2012), pMID: 26605735, https://doi.org/10.1021/ct3005613 .
  * Shirley (1998) E. L. Shirley, Phys. Rev. Lett. 80, 794 (1998).
  * Laskowski and Blaha (2010) R. Laskowski and P. Blaha, Phys. Rev. B 82, 205104 (2010).
  * Olovsson  _et al._ (2009) W. Olovsson, I. Tanaka, T. Mizoguchi, P. Puschnig, and C. Ambrosch-Draxl, Phys. Rev. B 79, 041102 (2009).
  * Roychoudhury  _et al._ (2022) S. Roychoudhury, L. A. Cunha, M. Head-Gordon, and D. Prendergast, Phys. Rev. B 106, 075133 (2022).
  * Taillefumier  _et al._ (2002) M. Taillefumier, D. Cabaret, A.-M. Flank, and F. Mauri, Phys. Rev. B 66, 195107 (2002).
  * Liang  _et al._ (2017) Y. Liang, J. Vinson, S. Pemmaraju, W. S. Drisdell, E. L. Shirley, and D. Prendergast, Phys. Rev. Lett. 118, 096402 (2017).
  * Liang and Prendergast (2018) Y. Liang and D. Prendergast, Phys. Rev. B 97, 205127 (2018).
  * Forestell and Marsiglio (2015) L. Forestell and F. Marsiglio, Canadian Journal of Physics 93, 1009 (2015).


# Spanning Trees of Recursive Scale-Free Graphs

C. Tyler Diggans  Clarkson Center for Complex Systems Science, Clarkson
University, Potsdam, NY 13699  Department of Physics, Clarkson University,
Potsdam, NY 13699  Air Force Research Laboratory: Information Directorate,
Rome, NY 13441  Erik M. Bollt  Clarkson Center for Complex Systems Science,
Clarkson University, Potsdam, NY 13699  Department of Electrical and Computer
Engineering, Clarkson University, Potsdam, NY 13699  Daniel ben-Avraham
Clarkson Center for Complex Systems Science, Clarkson University, Potsdam, NY
13699  Department of Physics, Clarkson University, Potsdam, NY 13699

###### Abstract

We present a link-by-link rule-based method for constructing all members of
the ensemble of spanning trees for any recursively generated, finitely
articulated graph, such as the DGM net. The recursions allow for many large-
scale properties of the ensemble of spanning trees to be analytically solved
exactly. We show how a judicious application of the prescribed growth rules
selects for certain subsets of the spanning trees with particular desired
properties (small-world, extended diameter, degree distribution, etc.), and
thus approximates and/or provides solutions to several optimization problems
on undirected and unweighted networks. The analysis of spanning trees enhances
the usefulness of recursive graphs as sophisticated models for everyday life
complex networks.

##  I Introduction

Recursive, finitely articulated models of complex networks mimic many key
properties of real-world networks. At the same time, these models are amenable
to exact analysis, providing important insights into the nature of everyday
life complex networks. In this article, we revisit the subject of spanning
trees for recursive scale-free nets, and provide a method for generating all
spanning trees, enabling a more complete analysis of this important set of
subgraphs.

The link-based method presented here is applicable to any recursive generative
network model that produces finitely articulated graphs, but we focus on the
well-known DGM net DGM  as a representative example of this class. This
particular family of networks was introduced in DGM  as a Pseudo-fractal
Scale-free Web (and was referred to as the PSW in subsequent literature), but
is equivalent to the $(u,v)$-flower graph Rozenfeld06  when $u=1$ and $v=2$.
We use our notation, $F_{\nu}(1,2)$, to denote the $\nu$-th generation of this
network Rozenfeld06 , for which there are two equivalent ways to construct it
recursively: the link-by-link approach and the hub-pasting approach. Both
recursive growth processes create the same set of small-world complex networks
with hierarchical structure, a scale-free degree distribution with power law
exponent $\gamma\approx 2.57$, and a high clustering coefficient of $C=4/5$,
matching several characteristic benchmarks of real-life complex nets. The
recursive nature and the finite articulation have enabled exact analytical
studies of the statistics of cycles Rozenfeld04 , diffusion Bollt05 ,
percolation Rozenfeld07 , spectral properties Xie16 , and minimum dominating
sets Shan17  for this network, shedding light on the analogous properties in
real-life complex nets. Previous work has exploited the hub-pasting approach
of constructing this and similar models to enumerate the spanning trees (and
calculate the subsequent tree entropy Lyons05  of the ensembles) Zhang10_PSW ;
Zhang11_trees ; Zhang14_trees ; Sun16_trees , but additional results related
to spanning trees of these networks are now obtained by considering a link-by-
link approach.

After providing a brief description of both methods of constructing the
original $F_{\nu}(1,2)$ in section II, a few new results on the spanning trees
of these graphs will be provided from the hub-pasting viewpoint in section
III, before turning to the link-based approach. In section IV, we provide an
explicit process for efficiently generating the whole ensemble of spanning
trees. We use the same recursive link-by-link process that creates the
original network, except that the links are assigned an $R$ or $G$, which
stand for “real” or “ghost” respectively. The term ghost indicates that a link
is part of the original graph, but not part of the spanning tree, and we
represent these links by dashed lines in all the figures.

Using the link-based approach, we first verify previous results for the
enumeration of the labeled spanning trees of $F_{\nu}(1,2)$ in a new
simplified way and calculate the tree entropy of the general $(u,v)$-flower
for comparison with previous results. We then show that this approach finds
added utility in solving optimization problems on undirected and unweighted
networks of this type. While the well known Minimum Spanning Tree (MST)
problem on weighted networks has garnered a lot of attention Graham85 ,
because of its relevance to the famed traveling salesman problem Held70 ;
Magnanti84  and its useful applications in data science Shapiro95 ; Xu02 ;
Varma04 ; He08 ; Galluccio12 ; Maeng12 ; Rainbolt17 ; the MST is not clearly
defined for unweighted networks. Common generalizations do exist, however, for
the undirected and unweighted case including the dense, sparse, and minimum
routing cost spanning tree problems, as well as many other lesser studied
variations Wu04 ; Dobrynin01 ; Kim04 ; Salamon08 ; Cerulli09 ; Silva14 ;
Szekely16 ; Li17 ; Silvestri17 ; Merabet18 ; Ozen20 ; Bozovic21 . The
solutions to these problems are relevant to an ever widening diversity of
applications from optical network design Gargano02 ; Marin15  to networked
oscillator synchronization Diggans21 . We show that while these combinatorial
problems are often NP-hard in general Gargano02 ; Salamon08 ; Silvestri17 ,
for recursive, finitely articulated networks, such optimal spanning trees (and
indeed the whole ensemble of spanning trees) can now be created alongside the
original network with minimal overhead. We conclude this article by describing
a process for selecting such optimal spanning trees through imposing
restrictions on the application of the link-based rules, and provide an
interesting open problem that we hope this work sheds light on.

##  II Constructing $F_{\nu}(u,v)$

Starting from an initial seed graph of the complete graph on two vertices,
i.e. $K_{2}$, in each generation of the growth process, we replace each link
from the previous generation with a pair of paths of length $u$ and $v$, in
parallel, between the nodes of the original link. For the case of $u=1$ and
$v=2$, this is equivalent to just adding a new path of length two in parallel
to the existing link as shown in Fig. 1 (a), which illustrates the first three
generations of $F_{\nu}(1,2)$.

For this case, it follows that for each link from the previous generation, one
new node is added to the network, while $w=u+v=3$ links replace each old link
in the new generation. Thus, the size (total number of links), $M_{\nu}$, and
order (total number of nodes), $N_{\nu}$, of the $\nu$-th generation are given
by $M_{\nu}=w^{\nu}=3^{\nu}$ and

$N_{\nu}=\frac{w-2}{w-1}w^{\nu}+\frac{w}{w-1}=\frac{1}{2}3^{\nu}+\frac{3}{2},$
(1)

where detailed derivations through recursion relations are provided elsewhere
DGM ; Rozenfeld06 ; Diggans20 .

Alternately, one can construct $F_{\nu}(1,2)$ using a hub-pasting approach by
taking $w=u+v=3$ copies of $F_{\nu-1}(1,2)$ and identifying one hub vertex
(those of highest degree) on each copy with another hub vertex on each other
copy resulting in $F_{\nu}(1,2)$. This process is shown in Fig. 1 (b) for the
creation of $F_{3}(1,2)$ by hub-pasting three copies of $F_{2}(1,2)$, where
the identified hub nodes are indicated by dashed circles.

|
---|---
(a) | (b)
Figure 1: Two equivalent methods of constructing $F_{\nu}(1,2)$ from
$F_{\nu-1}(1,2)$: (a) the first three generations of $F_{\nu}(1,2)$ are
obtained by adding a path of length $2$ in parallel to each link from the
previous generation and (b) the construction of $F_{3}(1,2)$ is obtained by
identifying the hub vertices of $u+v=3$ copies of $F_{2}(1,2)$; the three
nodes resulting from the hub identifications (the new hubs) are indicated by
dashed circles for clarity.

Prior to exploring the main contribution of this work within the context of
the link-by-link approach, we first mention two points that have not
previously been discussed, but are best described within the hub-pasting
paradigm.

##  III Hub-pasting Construction and the Typical Diameter

For generating spanning trees of $F_{\nu}(1,2)$, one can begin by selecting
(with replacement) three members from the ensemble of trees from generation
$\nu-1$ and identifying them at the vertices that would be considered hubs of
$F_{\nu-1}(1,2)$ as shown in Fig 2 (a). We will use the italicized version of
the word hub to indicate these hub nodes of the original network when in
reference to a spanning tree, whether those nodes in fact have the largest
degree in the tree or not. This pasting process is done explicitly by
temporarily adding a link between the hub vertices from adjacent spanning
trees (chosen from the previous generation), and then contracting the graph
along these temporary links to identify the three pairs of hub nodes into the
hubs for the new generation. Due to the connectedness of the $3$ individual
spanning trees, the identification process will form a single cycle through
the new hub nodes (shown in bold). To obtain a spanning tree of the new
generation, remove one of the links from this cycle. The ambiguity of this
link removal step complicates the enumeration of spanning trees, which led to
a more arduous accounting process used in previous work Zhang10_PSW . However,
this simplified description provides insight into the typical diameter of the
spanning trees of $F_{\nu}(1,2)$.

In Rozenfeld04 , the typical cycle length of $F_{\nu}(1,2)$ was shown to scale
as $2^{\nu}$ in the thermodynamic limit of large graphs. This implies that the
typical spanning tree of $F_{\nu}(1,2)$ contains a path having length on the
order of $2^{\nu}$. Given that the maximum diameter of any tree is also of
order $2^{\nu}$ (see section V), then the typical spanning tree must have a
diameter of order $2^{\nu}$. It is interesting that while the $(1,2)$-flower
is a small-world network, the diameter of the typical spanning tree of this
network scales like $\sim 2^{\nu}\sim N^{\log_{3}2}$ (where the order of
$F_{\nu}(1,2)$ from Eq. (1) is $N=\frac{1}{2}3^{\nu}+\frac{3}{2}\sim 3^{\nu}$)
and is therefore not small-world. However, section V shows how to select a
small-world spanning tree from the (typically non small-world) ensemble.

Additionally, we note that the “copy machine” algorithm outlined in Appendix C
of Bollt05 , which is based on the hub-pasting approach, can be adapted to
efficiently generate adjacency matrices of spanning trees of $F_{\nu}(1,2)$ in
a manner similar to what was outlined for ensembles in Diggans20 .

##  IV Link-by-link Construction and Tree Entropy

The link-by-link approach to spanning tree construction will mirror the
general link-by-link construction process of the recursive net, while keeping
track of two different kinds of links: “real” links ($R$) that belong to the
spanning tree, and “ghost” links ($G$) that belong to the original network,
but not the spanning tree.

For $F_{\nu}(1,2)$, in each generation, an $R$-link is either replaced by a
$G$-link and two $R$-links (rule # 1 in Fig. 2 (b)), or the $R$-link is
replaced by an $R$-link and a path of one $R$-link and one $G$-link in either
order (rule # 2). Additionally, in the process going forward, a $G$-link is
always replaced by a $G$-link and a path of one $R$-link and one $G$-link
(rule # 3). Other recursive scale-free network models such as the generic
$(u,v)$-flower graph Rozenfeld06  or a similar construction described in
Zhang07  may require a larger set of rules with more cases to consider, making
the subsequent analysis more tedious, but the process is, in essence, the
same.

---
(a)
(b)
Figure 2: As with the regular flower graph, there are two methods of
constructing spanning trees of $(1,2)$-flower graphs: (a) the hub-pasting
approach (here requiring an additional step) entails choosing $1+2=3$ elements
from the previous generation of spanning trees and identifying pairs of hub
nodes as shown; this, however, results in a single cycle (bolded), which must
be broken by removing one of its links. (b) The link-by-link approach involves
applying one of three rules to each link in each generation recursively. Each
sequence of valid rule operations results in a unique spanning tree, allowing
for the enumeration and generation of all spanning trees in the ensemble.

Denoting by $R_{\nu}$ ($G_{\nu}$) the number of real (ghost) links in
generation $\nu$, the construction rules in Fig. 2 imply

$\displaystyle R_{\nu+1}$ $\displaystyle=2R_{\nu}+G_{\nu},\qquad R_{0}=1$ (2)
$\displaystyle G_{\nu+1}$ $\displaystyle=R_{\nu}+2G_{\nu},\qquad G_{0}=0,$

leading to $R_{\nu}=\frac{1}{2}\left(3^{n}+1\right)$ and
$G_{\nu}=\frac{1}{2}\left(3^{n}-1\right)$. This is consistent with the fact
that $R_{\nu}=N_{\nu}-1$ and $G_{\nu}+R_{\nu}=M_{\nu}$, where $M_{\nu}$ and
$N_{\nu}$ are as defined in Eq. (1).

Each sequence of link updatings results in a single unique spanning tree;
moreover, each updating rule is uniquely reversible. Finally, the reverse
rules applied (in an appropriate manner) to an arbitrary spanning tree of
$F_{\nu}(1,2)$ generate a spanning tree of $F_{\nu-1}(1,2)$ (i.e. no loops are
created, and neither is the tree disconnected by the reverse rule
application). It follows that there is a one-to-one correspondence between
spanning trees and valid generating sequences of rules, and therefore, these
rules generate the complete set of labeled spanning trees.

This approach enables an independent verification of the number of labeled
spanning trees of generation $\nu$ for the DGM net Zhang10_PSW . Denote the
number of spanning trees by $T_{\nu}$, and consider the recursion

$T_{\nu+1}=\left(3^{R_{\nu}}\cdot 2^{G_{\nu}}\right)T_{\nu},$ (3)

where the term in parenthesis accounts for the 3 (2) ways of updating $R$
($G$) links. Together with the expressions for $R_{\nu}$ and $G_{\nu}$ and the
initial condition of $T_{0}=1$, we get

$T_{\nu}=\left(\frac{3}{2}\right)^{\frac{\nu}{2}}6^{\left(3^{\nu}-1\right)/4}.$

In terms of the order of the graph, $N$,

$T_{\nu}=\left(\sqrt{\frac{3}{2}}\right)^{\log_{3}{(2N-3)}}6^{\left(\frac{N}{2}-1\right)}\sim\frac{2^{C}}{6}~{}N^{C}~{}\sqrt{6}^{N},$

where $C=(1-\log_{3}(2))/2$ is a constant. Although a different numbering
convention for the generation was used here, this confirms the result found in
Zhang10_PSW .

For comparison, the total number of all labeled trees of order $N$ is
$N^{N-2}$ Cayley1889 . Amusingly, this functional dependence of $T_{\nu}$ with
$N$ is closer to that of unlabeled trees of order $N$ Otter48 .

For further comparison, this approach can be adapted to the general
$(u,v)$-flower by letting $w=u+v$. The recursion for the enumeration of all
labelled spanning trees [analogous to Eq. (3)] is then:
$T_{\nu+1}=w^{R_{\nu}}\cdot(uv)^{G_{\nu}}T_{\nu}$ with $T_{0}=1$, since there
are $w$ ways that any $R$-link can be propagated (the spanning trees of a $w$
cycle) and $uv$ ways that any $G$-link can be propagated (one link in each
path having lengths $u$ and $v$ must be a ghost link). A recursion similar to
Eq. (2) then gives $R_{\nu}=\frac{w-2}{w-1}w^{\nu}+\frac{1}{w-1}$ and
$G_{\nu}=\frac{1}{w-1}w^{\nu}-\frac{1}{w-1}$, which leads to the closed form
solution

$T_{\nu}=w^{\frac{w-2}{w-1}\left(\frac{w^{\nu}-1}{w-1}+\frac{\nu}{w-2}\right)}\cdot(uv)^{\frac{1}{w-1}\left(\frac{w^{\nu}-1}{w-1}-\nu\right)}.$

This gives the tree entropy Lyons05  of the $(u,v)$-flower as

$\displaystyle\lim_{N\rightarrow\infty}\frac{\ln{T_{\nu}}}{N}$
$\displaystyle=\frac{ln{(w)}}{w-1}+\frac{\ln{(uv)}}{(w-1)(w-2)},$

So, while the $(1,2)$-flower indeed has a small tree entropy when compared
with other networks having an average degree of $\langle k\rangle=4$ as
described in Zhang10_PSW , it is also true that larger values of $w$ will
result in smaller tree entropies tending toward zero in the limit of large
$w$; although the average degrees of flower graphs will tend toward $\langle
k\rangle=2$ as $w\rightarrow\infty$. Finally, we note that the tree entropy
for the $w$-cycle is simply $\frac{\ln{(w)}}{w}$, meaning the tree entropy of
the flower graphs are increased by an additional term that is quadratic in
$\frac{1}{w}$ due to the overlapping $w$-cycles that define them.

##  V Spanning Tree Subsets having optimal characteristics

In many applications of spanning trees, one is not concerned with simply the
number of trees, but one seeks an optimal spanning tree for some purpose,
e.g., minimal branching vertex spanning tree (MBVST), dense spanning tree
(DST), etc. Within the process outlined, due to the finite articulation of the
models being addressed, large-scale properties can be defined in terms of
recursion relations allowing such optimal trees to be selected from various
subsets of the ensemble of spanning trees that share certain properties. This
selection is done by restricting how the rules are applied, which alters the
recursions that define these properties.

To illustrate, consider the maximum diameter of any spanning tree from the
ensemble in generation $\nu$, which we denote by $s_{\nu}^{\textrm{max}}$. If
rule # 1 is applied to all $R$ links, then for each $R$ link in the path that
defined the diameter of the previous generation, we will have two $R$ links in
the path that contribute to the diameter of the new generation. At least for
generation $\nu>2$, there will also be two additional $G$ links that will
enable the diameter to be extended further, meaning we have the recursion
$s_{\nu}^{\textrm{max}}=2s_{\nu-1}^{\textrm{max}}+2$, with initial condition
$s_{2}^{\textrm{max}}=5$, leading to the formula
$s_{\nu}^{\textrm{max}}=7\cdot 2^{\nu-2}-2$. So, applying rule # 1 to all
$R$-links (excluding rule # 2 entirely) leads to the diameter scaling $\sim
2^{\nu}$, and in fact, the application of any finite proportion of rule # 1 to
$R$-links will result in a non small-world spanning tree of these small-world
networks.

These maximal diameter spanning trees are solutions to both the minimal
branching vertex and the sparse spanning tree problems (MBVST and SST
respectively) Gargano02 ; Cerulli09 ; Ozen20 ; Bozovic21 . So by applying rule
# 1 to all $R$-links and reducing the ambiguity of rule #3 by choosing the $R$
link to extend the diameter whenever possible (e.g. choosing the $R$ link to
be adjacent to nodes with smaller degree), a member of the subset of spanning
trees with maximal diameter is obtained, see Fig. 3 (a) for one of the subset
of eight such spanning trees of $F_{4}(1,2)$ having $s_{4}^{\textrm{max}}=26$
that are in fact solutions to the MVBST and SST problems.

---
(a)
(b)
(c)
Figure 3: Representative examples of subsets of spanning trees that can be
selected from within the set of spanning trees of $F_{\nu}(1,2)$: (a) using
rule #’s 1 and 3 only, where the ambiguity of rule # 3 is reduced by
maximizing the diameter when possible (breaking ties randomly) resulting in a
solution to the MBVST problem; (b) using rules #’s 2 and 3 only with $p=1$
resulting in a solution to the DST problem; and (c) using rule #’s 2 and 3
only with $p=0$ leads to a balanced small world tree with a more homogeneous
degree distribution that is a candidate for a solution to the $k$-MBVST
problem.

Contrast this spanning tree with a spanning tree that is created by never
applying rule # 1, which is a small-world spanning tree since the diameter
then scales as the logarithm of the order. In the class of trees where only
rule #’s 2 and 3 are applied, if one chooses $R$-links to be adjacent to the
largest degree node (breaking ties randomly), then one obtains a member of the
subset of spanning trees with the smallest diameter, i.e. a solution to the
dense spanning tree problem, which has minimal Wiener index Wiener47 ;
Szekely16 ; Ozen20  (see Fig. 3 (b)). In particular, we obtain a spanning tree
having the minimum diameter, defined by the recursion
$s_{\nu+1}^{\textrm{min}}=s_{\nu}^{\textrm{min}}+1$, since the ambiguity of
rule # 2 prevents any increase in the span, but there is one ghost link in
each generation that must add to the span. The initial condition of
$s_{1}^{\textrm{min}}=2$ leads to $s_{\nu}^{\textrm{min}}=\nu+1$, meaning in
terms of the order of the graph, we have $s_{\nu}^{\textrm{min}}\sim
log_{3}(N)$ in the limit of large graphs.

Letting a parameter $p$ control the proportion of instances (when applying
rule #’s 2 and 3) where the $R$-link is attached to the largest degree node,
the tree in Fig. 3 (b) is a representative of the set created using $p=1$; and
the tree in Fig. 3 (c) is representative of using $p=0$. The parameter $p$,
then alters the degree distribution while retaining the small-world property,
and values of $p\in(0,1)$ then define probability distributions over this
subset of trees having small Wiener indices, providing a useful tool for
approaching the Minimal $k$-Branching Vertex problem (k-MBVST) Merabet18 ;
Bozovic21 .

##  VI Conclusion

In summary, a general approach for obtaining spanning trees of recursive,
finitely articulated graphs has been presented through the example of the DGM
net, and basic properties of the ensemble of all spanning trees, such as the
typical/maximum diameter and the number of labeled trees are obtained exactly.
Additionally, one can now select a spanning tree from subsets of the ensemble
of spanning trees that share desirable properties such as small-world
features, certain degree distributions, and maximal/minimal diameter,
effectively solving many optimal spanning tree problems and approximating
solutions to others on such finitely articulated models.

Perhaps the most interesting related open question is to find the number of
unlabelled spanning trees of recursive finitely articulated graphs. While this
is, in general, an intractable problem, we hope that the recursive properties
of these graphs might provide a convenient foothold to approach such a
problem.

##  VII Acknowledgments

CTD has recieved funding from the Air Force Office of Scientific Research
(20RICOR010)

EB has received funding from the Army Research Office (N68164-EG) and also
DARPA

## References

  * [1] Sergey N. Dorogovtsev, Alexander V. Goltsev, and José F. F. Mendes.  “pseudofractal scale-free web”.  Phys. Rev. E, 65(6):066122, jun 2002.
  * [2] Hernán D. Rozenfeld, Shlomo Havlin, and Daniel ben-Avraham.  Fractal and transfractal recursive scale-free nets.  New J. Phys., 9, 12 2006.
  * [3] Hernán Rozenfeld, Joseph E. Kirk, Erik M. Bollt, and Daniel ben-Avraham.  Statistics of cycles: How loopy is your network?  J. Phys. A, 38, 03 2004.
  * [4] Erik M. Bollt and Daniel ben-Avraham.  “what is special about diffusion on scale-free nets”.  New J. Phys., 2005.
  * [5] Hernán D. Rozenfeld and Daniel ben Avraham.  Percolation in hierarchical scale-free nets.  Phys. Rev. E, 75:061102, 2007.
  * [6] Pinchen Xie, Zhong-Zhi Zhang, and Francesc Comellas.  On the spectrum of the normalized laplacian of iterated triangulations of graphs.  Appl. Math. Computation, 273:1123–1129, 2016.
  * [7] Liren Shan, Huan Li, and Zhong-Zhi Zhang.  Domination number and minimum dominating sets in pseudofractal scale-free web and sierpinski graph.  Theor. Comput. Sci., 677:12–30, 2017.
  * [8] Russell Lyons.  Asymptotic enumeration of spanning trees.  Combinatorics, Probability and Computing, 14(4):491–522, 2005\.
  * [9] Zhongzhi Zhang, Hongxiao Liu, Bin Wu, and Shuigeng Zhou.  Enumeration of spanning trees in a pseudofractal scale-free web.  EPL (Europhysics Letters), 90(6):68002, jun 2010.
  * [10] Zhongzhi Zhang, Hongxiao Liu, Bin Wu, and Tao Zou.  Spanning trees in a fractal scale-free lattice.  Physical Review E, 83(1):016116, 2011.
  * [11] Zhongzhi Zhang, Bin Wu, and Francesc Comellas.  The number of spanning trees in apollonian networks.  Discrete applied mathematics, 169:206–213, 2014.
  * [12] Weigang Sun, Shuai Wang, and Jingyuan Zhang.  Counting spanning trees in prism and anti-prism graphs.  J. Appl. Anal. Comput, 6(1):65–75, 2016.
  * [13] R.L. Graham and Pavol Hell.  On the history of the minimum spanning tree problem.  Annals of the History of Computing, 7(1):43–57, 1985.
  * [14] Michael Held and Richard M. Karp.  The traveling-salesman problem and minimum spanning trees.  Operations Research, 18(6):1138–1162, 1970.
  * [15] T. L. Magnanti and R. T. Wong.  Network design and transportation planning: Models and algorithms.  Transportation Science, 18(1):1–55, 1984.
  * [16] Larry S. Shapiro.  Clustering using maximum affinity spanning trees, page 61–99.  Distinguished Dissertations in Computer Science. Cambridge University Press, 1995.
  * [17] Xu Y., Olman V., and Xu D.  Clustering gene expression data using a graph-theoretic approach: an application of minimum spanning trees.  Bioinformatics, 4:536–45, 2002.
  * [18] Sudhir Varma and Richard Simon.  Iterative class discovery and feature selection using minimal spanning trees.  BMC Bioinformatics, 5(126), 2004.
  * [19] Qiang He and Chee-Hung Henry Chu.  A new image segmentation technique using maximum spanning tree.  In Valentin E. Brimkov, Reneta P. Barneva, and Herbert A. Hauptman, editors, Combinatorial Image Analysis, pages 197–204, Berlin, Heidelberg, 2008. Springer Berlin Heidelberg.
  * [20] Laurent Galluccio, Olivier Michel, Pierre Comon, and Alfred O. Hero.  Graph based k-means clustering.  Signal Processing, 92(9):1970–1984, 2012.
  * [21] Seong Eun Maeng, Hyung Wooc Choi, and Jae Woo Lee.  Complex networks and minimal spanning trees in international trade network.  Intl. J. of Modern Physics: Conference Series, 16:51–60, 2012.
  * [22] J. Lovelace Rainbolt and M. Schmitt.  The use of minimal spanning trees in particle physics.  Journal of Instrumentation, 12(02):P02009–P02009, feb 2017.
  * [23] B.Y. Wu and K.-M Chao.  Minimum routing cost spanning trees, pages 41–84.  Discrete Mathematics and its Applications. CRC Press, 2004.
  * [24] Andrey A. Dobrynin, Roger Entringer, and Ivan Gutman.  Wiener index of trees: Theory and applications.  Acta Applicandae Mathematicae, 66(3):211, 05 2001.
  * [25] Dong-Hee Kim, Jae Dong Noh, and Hawoong Jeong.  Scale-free trees: The skeletons of complex networks.  Phys. Rev. E, 70:046126, Oct 2004.
  * [26] Gábor Salamon and Gábor Wiener.  On finding spanning trees with few leaves.  Information Processing Letters, 105(5):164–169, 2008.
  * [27] Raffaele Cerulli, Monica Gentili, and Antonio Iossa.  Bounded-degree spanning tree problems: Models and new algorithms.  Computational Optimization and Applications, 42:353–370, 04 2009\.
  * [28] Ricardo Silva, Diego Silva, Mauricio Resende, Geraldo Mateus, José Gonçalves, and Paola Festa.  An edge-swap heuristic for generating spanning trees with minimum number of branch vertices.  Optimization Letters, 8, 04 2014.
  * [29] László Székely, Stephan Wagner, and Hua Wang.  Problems related to graph indices in trees, pages 3–30.  Springer International Publishing, 2016.
  * [30] Tao Li, Yingqi Gao, Qiankun Dong, and Hua Wang.  Degree sums and dense spanning trees.  PLOS ONE, 12(9):1–12, 09 2017.
  * [31] Selene Silvestri, Gilbert Laporte, and Raffaele Cerulli.  A branch-and-cut algorithm for the minimum branch vertices spanning tree problem.  Computers & Operations Research, 81:322–332, 2017.
  * [32] Massinissa Merabet, Jitamitra Desai, and Miklos Molnar.  A generalization of the minimum branch vertices spanning tree problem.  In Jon Lee, Giovanni Rinaldi, and A. Ridha Mahjoub, editors, Combinatorial Optimization, pages 338–351, Cham, 2018. Springer International Publishing.
  * [33] Mustafa Ozen, G. Lesaja, and Hua Wang.  Globally optimal dense and sparse spanning trees, and their applications.  Statistics, Optimization & Information Computing, 8:328–345, 05 2020.
  * [34] Vladimir Božović, Z̆ana Kovijanić Vukićević, Goran Popivoda, Rong-Ying Pan, and Xiao-Dong Zhang.  Extreme wiener indices of trees with given number of vertices of maximum degree.  Discrete Applied Mathematics, 304:23–31, 2021.
  * [35] Luisa Gargano, Pavol Hell, Ladislav Stacho, and Ugo Vaccaro.  Spanning trees with bounded number of branch vertices.  In Peter Widmayer, Stephan Eidenbenz, Francisco Triguero, Rafael Morales, Ricardo Conejo, and Matthew Hennessy, editors, Automata, Languages and Programming, pages 355–365, Berlin, Heidelberg, 2002. Springer Berlin Heidelberg.
  * [36] Alfredo Marín.  Exact and heuristic solutions for the minimum number of branch vertices spanning tree problem.  European Journal of Operational Research, 245(3):680–689, 2015\.
  * [37] C. Tyler Diggans, Jeremie Fish, Abd AlRahman R. AlMomani, and Erik M. Bollt.  The essential synchronization backbone problem.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(11):113142, 2021.
  * [38] C. Tyler Diggans, Erik M. Bollt, and Daniel ben Avraham.  Stochastic and mixed flower graphs.  Phys. Rev. E, 101:052315, May 2020.
  * [39] Zhong-Zhi Zhang, Li-Li Rong, and Shuigeng Zhou.  A general geometric growth model for pseudofractal scale-free web.  Physica A, 377(1):329 – 339, 2007.
  * [40] Arthur Cayley.  A theorem on trees, volume 13 of Cambridge Library Collection - Mathematics, page 26–28.  Cambridge University Press, 2009.
  * [41] Richard Otter.  The number of trees.  Annals of Mathematics, 49(3):583–599, 1948.
  * [42] Harry Wiener.  Structural determination of paraffin boiling points.  Journal of the American chemical society, 69(1):17–20, 1947.