Current affiliation: ]Physik-Institut, Universität Zürich, Winterthurerstrasse
190, CH-8057 Zürich, Switzerland

# Intercalation-induced states at the Fermi level and the coupling of
intercalated magnetic ions to conducting layers in Ni1/3NbS2

Yuki Utsumi Boucher<EMAIL_ADDRESS>Institute of Physics, Bijenička c. 46,
10000 Zagreb, Croatia  Izabela Biało  [  AGH University of Krakow, Faculty of
Physics and Applied Computer Science, 30-059 Krakow, Poland  Mateusz A. Gala
AGH University of Krakow, Faculty of Physics and Applied Computer Science,
30-059 Krakow, Poland  Wojciech Tabiś  AGH University of Krakow, Faculty of
Physics and Applied Computer Science, 30-059 Krakow, Poland  Marcin Rosmus
Solaris National Synchrotron Radiation Centre, Jagiellonian University,
Czerwone Maki 98, 30-392 Krakow, Poland  Faculty of Physics, Astronomy, and
Applied Computer Science, Jagiellonian University, Łojasiewicza 11, 30-348
Krakow, Poland  Natalia Olszowska  Solaris National Synchrotron Radiation
Centre, Jagiellonian University, Czerwone Maki 98, 30-392 Krakow, Poland
Jacek J. Kolodziej  Solaris National Synchrotron Radiation Centre,
Jagiellonian University, Czerwone Maki 98, 30-392 Krakow, Poland  Faculty of
Physics, Astronomy, and Applied Computer Science, Jagiellonian University,
Łojasiewicza 11, 30-348 Krakow, Poland  Bruno Gudac  Department of Physics,
Faculty of Science, University of Zagreb, Bijenička c. 32, 10000 Zagreb,
Croatia  Mario Novak  Department of Physics, Faculty of Science, University of
Zagreb, Bijenička c. 32, 10000 Zagreb, Croatia  Naveen Kumar Chogondahalli
Muniraju  Institute of Physics, Bijenička c. 46, 10000 Zagreb, Croatia
Institute of Nuclear Physics PAN, Radzikowskiego 152, 31-342 Kraków, Poland
Ivo Batistić  Department of Physics, Faculty of Science, University of Zagreb,
Bijenička c. 32, 10000 Zagreb, Croatia  Neven Barišić  Department of Physics,
Faculty of Science, University of Zagreb, Bijenička c. 32, 10000 Zagreb,
Croatia  Institute of Solid State Physics, TU Wien, 1040 Vienna, Austria
Petar Popčević<EMAIL_ADDRESS>Institute of Physics, Bijenička c. 46, 10000
Zagreb, Croatia  Eduard Tutiš<EMAIL_ADDRESS>Institute of Physics, Bijenička
c. 46, 10000 Zagreb, Croatia

###### Abstract

The diversity of magnetic orders that appear in layered magnetic materials is
of great interest from the fundamental point of view and for applications. In
particular, the magnetic sublayers, introduced by intercalation into van der
Waals gaps of the host transition-metal dichalcogenide (TMD), are known to
produce various magnetic states, with some being tunable by pressure and
doping. The magnetic sublayers and their magnetic ordering strongly modify the
electronic coupling between layers of the host compound. Understanding the
roots of this variability, starting from the underlying electronic structure,
is a significant challenge. Here we employ the angle-resolved photoelectron
spectroscopy at various photon energies, the ab initio electronic structure
calculations, and modeling to address the particular case of Ni-intercalate,
Ni1/3NbS2. We find that the bands around the Fermi level bear the signature of
a strong yet unusual hybridization between NbS2 conduction band states and the
Ni 3$d$ orbitals. The hybridization between metallic NbS2 layers is almost
entirely suppressed in the central part of the Brillouin zone, including the
part of the Fermi surface around the $\mathrm{\Gamma}$ point. Simultaneously,
it gets very pronounced towards the zone edges. It is shown that this behavior
is the consequence of the rather exceptional, symmetry imposed, spatially
strongly varying, zero total hybridization between relevant Ni magnetic
orbitals and the neighboring Nb orbitals that constitute the metallic bands.
We also report the presence of the so-called $\beta$-feature, discovered only
recently in two other magnetic intercalates with very different magnetic
orderings. In Ni1/3NbS2, the feature shows only at particular photon energies,
indicating its bulk origin. Common to prior observations, it appears as a
series of very shallow electron pockets at the Fermi level, positioned along
the edge of the Brillouin zone. Unforeseen by ab initio electronic
calculations, and its origin still unresolved, the feature appears to be a
robust consequence of the intercalation of 2H-NbS2 with magnetic ions.

ARPES; Ni intercalation; band structure; transition-metal dichalcogenides;
resonance

##  I Introduction

For decades, the interaction between localized magnetic moments and itinerant
electronic subsystems has been central to condensed matter research. In
particular, it has been regularly and extensively discussed in the context of
heavy fermions [1] and high-temperature superconducting cuprates, whose
discovery has posed many puzzles to the research community [2, 3, 4, 5, 6].
Even earlier, various materials, where the same interaction may be expected,
have been discovered and cataloged after inquiries appropriate for the times.
The magnetic intercalates of transition metal dichalcogenides (TMDs) are one
such family of compounds [7, 8, 9]. Probably the best known and most numerous
of those are based on 2H-NbS2, with formula M1/3NbS2 (M: 3d transition metal).
These compounds are particularly attractive due to the high quality of the
crystals achieved at the 1/3 concentration of intercalated atoms and the range
of exciting intertwined magnetic and transport properties they exhibit [9, 10,
11, 12]. In recent times, along with the development of new experimental
techniques and theoretical concepts, there has been a revived interest in the
2H-NbS2 intercalates. With the interaction mentioned above taking the
limelight, the question of the electronic properties of these materials seems
far from being settled. Apart from the common questions regarding the charge
transfer between metallic layers of the host compound and the intercalated
species, and the resulting spin and charge states of the latter, the other
questions remain to be addressed: How do the magnetic ions connect the
otherwise weakly connected metallic layers, and how their magnetic states and
fluctuations affect the electronic structure in the vicinity of the Fermi
level. Some of these questions are to be addressed here.

The host compound 2H-NbS2 is a superconductor below 6 K [13]. The
intercalation of transition metal ions between NbS2 layers suppresses the
superconductivity and often results in complex magnetic orderings. In
M1/3NbS${}_{2\ }$the transition metal ions occupy the octahedral voids between
NbS2 layers and form a $\mathrm{\sqrt{3}}\mathrm{\times}$$\mathrm{\sqrt{3}}$
superlattice rotated by 30∘ in respect to the 1$\mathrm{\times}$1 unit cell of
2H-NbS2. Intercalation of the late transition metals Fe, Co, and Ni leads to
predominant antiferromagnetism, whereas the early transition metals V, Mn, and
Cr tend to order ferromagnetically [9, 10, 11, 12]. In some cases, the type of
long-range magnetic order is found tunable by small changes in the
concentration of the intercalate [14]. Recent studies have revealed the chiral
helimagnetic structure in Cr1/3NbS2 and Mn1/3NbS2 [15, 16], whereas weak
ferromagnetism along the c-axis was found to accompany the dominant
antiferromagnetic ordering in Co1/3NbS2 [17]. Furthermore, a large anomalous
Hall effect has been reported in the latter compound, the origin of which is
still debated [17, 18, 19, 20, 21].

The angle-resolved photoelectron spectroscopy (ARPES) is expected to help
significantly in understanding the electronic properties of these crystals.
The current study uses ARPES to address the case of Ni1/3NbS2. Along the way,
we will contrast this compound to Co1/3NbS2, recently studied by several
groups [21, 22, 23]. At first, going from Co-intercalated to Ni-intercalated
material can appear equivalent to doping, imposing a shift of the Fermi level
within the almost-rigid electronic band structure of the host compound. A
deeper comparison finds more significant and intricate differences, and the
difference in band fillings in the two compounds is not one of them. Notably,
the antiferromagnetic orderings in Ni1/3NbS2 and Co1/3NbS2 differ
significantly. The magnetic structure of Co1/3NbS2 was solved previously [12]
and confirmed recently [19]. The reported antiferromagnetic “hexagonal
ordering of the first kind” breaks the 120${}^{\circ}\ $rotational symmetry of
the crystal and doubles the unit cell. The only attempt of refining the
magnetic structure in Ni1/3NbS2 resulted in the claim that no additional
magnetic peaks are observed in the antiferromagnetically ordered phase [24].
Such a finding is in accordance with the theoretically proposed magnetic
structure [25] where Ni magnetic moments are ferromagnetically arranged within
Ni-layers, whereas their direction alternates between subsequent Ni-layers. In
that way, the magnetic unit cell stays the same as the crystallographic one,
as the crystallographic unit cell already comprises two Ni layers. The same
magnetic structure will be used in our calculations of the electronic
structure of Ni1/3NbS2, to be presented further below. Also, in contrast to
Co1/3NbS2, no evidence of the out-of-plane ferromagnetic moment nor anomalous
Hall effect has been reported for Ni1/3NbS2. Other differences and
similarities should be sought in the electronic structure, which we explore
here. We will mainly focus on the states near the Fermi level, dominantly
influencing the low-energy electronic process. Motivated by the recent finding
of a new "$\beta$-feature" in the vicinity of the Fermi level, around the
$\mathrm{K}$-point of the first Brillouin zone in Cr1/3NbS2 and Co1/3NbS2 [26,
27, 21, 22, 23, 28], we also explore its possible existence in Ni1/3NbS2.
Experimentally detected, the feature is unforeseen by any electronic structure
calculation, whereas the arguments in favor of its origin in strong electronic
correlations have been presented recently [23]. Despite the mechanism of its
appearance not being fully understood, the $\beta$-feature has been observed
only in compounds with intercalated magnetic ions, with the orbitals of the
intercalated magnetic ions participating in the phenomenon [26, 27, 23].

In this work, we devote special attention to the influence of intercalated
ions on the Nb 4d conduction bands. Recently a large amplification of the
interlayer hybridization has been found in Co1/3NbS2 and Cr1/3NbS2 [23, 28],
effectively changing the character of these bands near the $\mathrm{\Gamma}$
point from quasi-two to three-dimensional. However, it will turn out that the
coupling between NbS2 layers and the intercalated magnetic ions works very
differently in Co1/3NbS2 and Ni1/3NbS2.

##  II Experiment and Methods

Single crystals of Ni1/3NbS2 were grown by chemical vapor transport method
using iodine as a transport agent [9]. The grown single-crystals were
characterized by electrical resistivity measurement, x-ray diffraction, and
magnetic susceptibility. Electrical resistivity and magnetic susceptibility
measurements revealed the antiferromagnetic ordering temperature, $T_{N}=$95
K, higher than earlier reported 90 K [9, 10]. The residual resistivity ratio
(R(300 K)/R(2 K)) of 19, relatively high for this class of materials [9, 11],
also indicates the high quality of crystals. The measured room temperature
lattice parameters are a = 5.7609 (5) Å, and c = 11.9014 (7) Å and agree well
with previously reported values [8]. Using the Laue diffraction, samples were
oriented along the high symmetry $\mathrm{\Gamma}$M0 direction (which
corresponds to $\mathrm{\Gamma}$M direction of 2H-NbS2 Brillouin zone) and
glued to Ti flat plates. The pristine surfaces were obtained by in situ
cleaving. A sample was placed within the preparation chamber, cleaved under
ultra-high vacuum (better than 2$\mathrm{\times}$10-10 mbar), and immediately
transferred to the measurement chamber with the base pressure of
8$\mathrm{\times}$10-11 mbar. The experiment was performed at the UARPES
(currently URANOS) beamline of the SOLARIS synchrotron, using the photons in
the 34–82 eV energy range. The measurements were conducted at 20 K, well below
the magnetic ordering temperature. The photoelectrons were collected using the
hemispherical analyzer (Scienta-Omicron DA30) and the multichannel plate
coupled to a CCD detector. The overall energy resolution was set to
$\mathrm{\sim}$30 meV. Linearly polarized light was used in the experiment.
The electric field vector of the incident light was parallel to the incident
plane defined by the light propagation vector and was normal to the sample
surface. The details of the experimental setup can be found in Ref. 23.

The band structure for Ni1/3NbS2 was calculated using the Quantum ESPRESSO
(QE) Density Functional Theory package [29, 30]. We have used the kinetic
energy cut-off of 160 Ry for wavefunctions and 800 Ry for charge density and
potentials. To some level, the electronic correlations on Ni ions were
accounted for using the DFT+$U$ approach proposed by Cococcioni and Gironcoli
[31]. The Hubbard interaction, $U$ = 7.7 eV, was determined through first-
principles calculation [32]. The ultrasoft pseudopotentials from
pslibrary-0.3.1 [33] are based on Perdew-Burke-Ernzerhof exchange-correlation
functional revised for densely packed solids [34, 35]. The k-mesh within the
Brillouin zone was $10\times 10\times 5$, without shift, whereas the Fermi-
energy discontinuity was smeared as proposed by Marzari-Vanderbilt [36] using
the broadening of 0.01 Ry.

##  III 40 eV ARPES and DFT spectra

Figure 1: (Color online) The ARPES signal in $k$-space. Here and throughout
the paper, $\mathrm{\Gamma}$, K0, and M0 denote the usual characteristic
points of the first Billiouin zone of the 2H-NbS2, whereas K stands for the
usual characteristic point of the three-fold smaller first Billiouin zone of
Ni1/3NbS2. $\overline{\mathrm{\Gamma}}$, ${\overline{\mathrm{M}}}_{0}$,
${\overline{\mathrm{K}}}_{0}$, and ${\overline{\mathrm{K}}}$ represent the
corresponding points of the surface Brillouin zones. (a) ARPES intensity plot
of Ni1/3NbS2 along $\overline{\mathrm{\Gamma}}{\overline{\mathrm{M}}}_{0}$
direction measured with h$\nu$ = 40 eV is shown on the right. The measured
spectra are compared to the unfolded DFT-calculated spectra for magnetically
ordered Ni1/3NbS2 (shown on the left). The vertical dashed lines represent the
boundaries of the Ni1/3NbS2 surface Brillouin zone. The red (light green) dots
on the left panel mark the electronic states with predominant contributions
from niobium (sulfur) orbitals. (b) Fermi surface map of Ni1/3NbS2, also
obtained by using h$\nu$ = 40 eV. The black hexagon and the smaller blue
hexagon shown in the dashed line represent the first Brillouin zone of 2H-NbS2
and Ni1/3NbS2, respectively. The inset shows another scan that extends to the
boundary of the (big) first Brillouin zone of 2H-NbS2

Figure 1 shows the ARPES intensity plot obtained using the incident photon
energy (h$\nuup$) of 40 eV, along with the results of the performed density
functional theory (DFT) calculations. The right-hand side of Fig. 1 (a)
presents the ARPES intensity plot of Ni1/3NbS2 measured along
$\overline{\mathrm{\Gamma}}{\overline{\mathrm{M}}}_{0}$ direction. The
notation regarding the $k$-space points is explained in the caption of Fig. 1.
In agreement with other ARPES results in the M1/3NbS2 series [37, 26, 21, 22,
23], no evidence of the band splitting can be observed at the boundary of
Ni1/3NbS2 Brillouin zone (vertical dashed line in Fig. 1 (a)). However, the
part of the conduction band signal weakens significantly after crossing the
dashed boundary.

The DFT-calculated spectra are shown on the left-hand side of Fig. 1 (a)
111The conventional illustration of the calculated band structure of Ni1/3NbS2
can be found in Supplemental Material.. These spectra are shown as unfolded
into the Brillouin zone of 2H-NbS2 for a more straightforward comparison with
the experimental data. Two approaches to unfolding were used. In the most
direct approach, we use the unfold-x procedure [39, 40, 41] on the Quantum
Espresso results. The second approach is based on the tight-binding (TB)
parametrization of the DFT-calculated bands. The TB parametrization, also used
elsewhere in the paper, is based on the Wannier90 code [42]. The two
approaches produce the same results regarding the spectra in Fig. 1(a).
However, the TB approach may also be exploited to get insight into the orbital
content of particular electronic states. In that respect, the red and light
green dots in Fig. 1(a) mark the energies of electronic states with the
predominant contribution of niobium and sulfur orbitals, respectively. The two
red dots closest to the Fermi level mark the Nb 4$d$ bands. These signals
resemble the anti-bonding and the bonding Nb 4$d$ bands of 2H-NbS${}_{2\ }$
studied previously, both experimentally and through DFT calculations [26, 43,
44, 23].

The band fillings are notoriously hard to assess precisely from ARPES scans.
Even under the assumption of Luttinger’s theorem being applicable in a
particular system [45], the need arises to perform ARPES at the dense mesh of
$k$-points in three dimensions instead of addressing only the high symmetry
cross-sections. Including the several high-symmetry cross-sections, comparable
to $\mathrm{\Gamma}-\mathrm{M}_{0}$ cross-section in Fig. 1(a) helps only
slightly [23], and the experimental-only means to assess the bands-filling may
require extensive soft X-ray ARPES scans, providing more precise $k_{z}$
mapping of the electronic spectra [46]. In intercalated systems, the
observation of the full Fermi surface in the $k$-space is further obscured by
a large number of Fermi surface segments, a weak signal of some bands within
the first Brillouin zone of the intercalated material, coming from ARPES
matrix elements. In addition, and not unlike our system, the strong electron
correlations tend to add complications that go beyond a simple band structure.
In practice, the best course in estimating the band filling in the
intercalated system partly relies on comparing the ARPES spectra in the
intercalated and the host compound and comparing measured and DFT-calculated
spectra. In the intercalated systems, the comparison between the latter two is
further complicated by the fact that the intercalation is expected to perturb
the spectrum of the host material extending over the Brillouin zone, several
times bigger than the one calculated for the intercalated system. Thus, there
is a need for unfolding the calculated spectra to the original Brillouin zone
of the host material before making the comparison. In addition, in our
particular case, the problem with the precise determination of the Fermi
surface volume arises from the $\beta$-feature appearing in the measured
spectra, to be addressed in the following section, which lacks its DFT-
calculated counterpart.

The DFT results are used to calculate the charge of the Ni ion by integrating
the charge density around the Ni site, with the result of +1.6 $e$. The
numerical evaluation of the conduction bands filling from DFT results produces
the result of 5/6. Compared to the host compound 2H-NbS2, the extra 2/3
electrons per NbS2 formula unit is found in conduction bands. Within the rigid
band picture, this would require the charge transfer of precisely 2 electrons
from Ni ion into NbS2 layers. Thus, the result for the conduction bands
filling seems to oppose the result for the charge of the Ni ion of +1.6 $e$,
mentioned above. The resolution of this small puzzle lies in the fact that the
situation in Ni1/3NbS2 very much departs from the rigid band picture. The
strong hybridization between intercalated ions and the layers of the host
compound invalidates the simple relation between conduction bands filling and
the charge state of the intercalated ion. Further details are provided in
Appendix B.

Going to higher binding energies, we first address the band that separates
from the conduction bands approximately halfway between
${\overline{\mathrm{M}}}_{0}$ and $\overline{\mathrm{\Gamma}}$ points. In the
experimental side of Fig. 1(a), this band reaches the $\mathrm{\Gamma}$ point
at the binding energy ($E_{B}$) of 0.4 eV, i.e., at energy $E=E_{\rm
F}-E_{B}=-0.4$ eV, below the Fermi level ($E_{\rm F}$).

The band that shows similar behavior in the calculated spectrum reaches
$\mathrm{\Gamma}$ point at $E_{B}\approx$ 0.65 eV (marked by the higher light-
green dot in Fig. 1(a)). This band is the sulfur 3$p_{z}$ band, which is
particularly sensitive, especially around the $\mathrm{\Gamma}$ point, to
minor changes in the $c$-axis lattice constant. This high sensitivity has been
documented in the recent numerical studies in 2H-NbS2, and its intercalates
[43, 23]. Thus, the 0.25 eV shift in energy at the $\mathrm{\Gamma}$ point
should not prevent us from making the correspondence between the observed band
and the calculated band, otherwise behaving similarly. A similar shift in size
and direction can be observed for the sulfur 3$p$ signal around 1.6 eV binding
energy, marked by the lower light-green dot in Fig. 1(a). Finally, the lowest
red dot in Fig. 1(a) indicates the Nb-prevailing state around 3 eV at
$\mathrm{\Gamma}$ point and marks the only region in the lower part of the
spectrum where the contribution from sulfur orbitals is not dominating.

Further similarities and differences between ARPES results for Ni1/3NbS2 and
2H-NbS2 [26, 43] can be discussed in some detail, as it was done in the case
of Co1/3NbS2 [23]. Instead, we focus on several essential features the Ni
intercalation brings to the electronic spectrum.

Fig. 1 (b) shows the Fermi surface of Ni1/3NbS${}_{2\ }$measured using $h\nu$
= 40 eV. The large black hexagon and the small blue hexagon in the dashed line
represent the surface Brillouin zones of 2H-NbS2 and Ni1/3NbS2, respectively.
As in the previous ARPES studies of Ni1/3NbS2 [37, 21], the part of the Fermi
surface centered at the $\mathrm{\Gamma}$ point is present in the form of a
hexagonally shaped hole pocket. Fig. 1 (b) does not extend wide enough to make
visible the part of the Fermi surface centered around ${\overline{K}}_{0}$
points at the edge of the Brillouin zone of 2H-NbS2. This part of the Fermi
surface, evolving from the similar sections in 2H-NbS2, is presented in the
inset as a separate scan. On the other hand, keeping in mind our interest for
possible non-NbS2${}_{\ }$features appearing around the $\overline{K}$-point,
at the edge of the Ni1/3NbS2 Brillouin zone, we note that no such signal is
found in Fig. 1 (b). However, as pointed out earlier [23], the search is to be
continued by scans at other photon energies. This is precisely what we do
next, by measuring the ARPES spectra in the 34 to 82 eV photon energy range in
a few eV steps.

Figure 2: (Color online) ARPES intensity plots of Ni1/3NbS2 measured along
$\overline{\mathrm{\Gamma}}{\overline{\mathrm{M}}}_{0}$ direction using
various incident photon energies from 34 to 82 eV. The vertical dashed lines
represent the boundaries of the Ni1/3NbS2 Brillouin zone. The total intensity
of each ARPES spectrum, extending over a measured span of electronic momenta
and the binding energy between 0 and 4 eV, is subsequently normalized to the
same value. The signal very narrow in energy around the Fermi level, showing
at the left and the right edges of the 58 eV panel, is the first trace of the
$\beta$-feature, to be followed in other figures.

##  IV The $\beta$-feature in the nickel intercalate

Fig. 2 shows ARPES intensity plots of Ni1/3NbS2 measured along
$\overline{\mathrm{\Gamma}}{\overline{\mathrm{M}}}_{0}$ line using various
photon energies. What is immediately apparent are the changes between the
panels that occur in the S 3$p$ bands located between $E_{B}=1$ and 1.8 eV.
These changes indicate a strong $k_{z}$ dispersion of these bands. In
contrast, the Fermi level crossings of the Nb 4$d$ bands are almost
insensitive to the variation of the photon energy, pointing at the two-
dimensional character of the Fermi surface, as expected for very anisotropic
host material 2H-NbS2 [43].

The more important feature to be spotted in these spectra is the appearance of
the shallow electron pocket that shows at the Fermi level at the boundary of
Ni1/3NbS2 Brillouin zone (dashed vertical line), especially evident in
$h\nu$=58 eV panel. These shallow electron pockets resemble those observed in
Co1/3NbS2 [23, 21, 22], following the Brillouin zone boundary of the
intercalate, and exceptionally well pronounced around its $\overline{K}$-point
corners. Apart from Co1/3NbS2, a similar feature has been observed also in
Cr1/3NbS${}_{2\ }$[26, 27, 28]. The signal/pocket observed around the
$\overline{K}$-point was given the name “$\beta$-feature” [27, 23], which we
also adopt here.

Figure 3: (Color online) ARPES intensity plots along $k_{z}$ direction at (a)
the Fermi level and (b) 1.1 eV below the Fermi level along the cuts indicated
by the black dashed lines in Fig. 2. The vertical dashed lines represent the
boundary of the Brillouin zone of Ni1/3NbS2. The horizontal green lines and
labels ($\mathrm{\Gamma}_{n+1}$) mark the $k_{z}$ levels corresponding to
$\mathrm{\Gamma}_{n+1}$ $\equiv$ (0, 0, 2$\pi/c\times n$). The horizontal blue
lines and labels $\mathrm{A_{n+1}}$ mark the $k_{z}$ levels corresponding to
$\mathrm{A}_{n+1}\equiv(0,0,\pi/c+2\pi/c\times n$).

It is somewhat easier to identify the $\beta$-feature in Fig. 3, which
presents ARPES intensity plots at constant energy along the kz ($c$-axis)
direction. The same data presented in Fig. 2 are now used to plot the kz band
dispersion shown in Fig. 3, with the $k_{z}$ values obtained from the photon
energy using the usual formula,
$k_{z}=\sqrt{(E_{k}{{\mathrm{cos}}^{2}\theta+V_{0})\ }2m/{{\hslash}^{2}}}$.
Here $E_{k}$, $\theta$, m, and V0 denote the kinetic energy of a
photoelectron, the angle between the surface normal and photoelectron’s wave-
vector, the free electron mass, and the inner potential, respectively [47].
The inner potential was set to V0 = 14 eV, previously used in cases of
Cr1/3NbS2 and Co1/3NbS2 [26, 23]. It is worth to mention that $V_{0}$ might
vary among different intercalates. The variation between ionic charges of
magnetic ions is prone to strongly affect the electrostatic potential energy
drop over the dipolar layer at the crystal surface. The versions of Fig.3 done
for the inner potential $V_{0}$ range from 14 eV to 19 eV are provided in the
Supplemental material. The figures do not differ very much, with the
differences being mostly perceived as shifts in $k_{z}$.

In the Fermi level scan in Fig. 3 (a), the $\beta$-feature appears close to
the Ni1/3NbS2 Brillouin zone boundary, marked by dashed vertical lines at
$k_{\parallel}$ ${}_{\ }$= $\mathrm{\pm}$0.727 Å-1
$(\overline{\mathrm{K}}$-point). This signal can be followed within the
kz${}_{\ }$= 4.0-4.2 Å-1 range around ${\mathrm{\Gamma}}_{9}$. Apart from this
$k_{z}$ dependence, the quasi-2D character of the Fermi surface originating
from “ordinary” Nb 4d bands is evident from Fig. 3 (a). The signals
corresponding to Nb 4d bands cross the Fermi level at $k_{\parallel}$${}_{\
}$= $\mathrm{\pm}$0.4 Å-1${}^{\ }$ throughout the measured $k_{z}$ range.

Fig. 3 (b) shows the constant energy surface cut at $E_{B}=1.1$ eV. An oval-
shaped constant energy surface elongated along the $k_{z}$ direction can be
clearly distinguished with a pronounced $k_{z}$ dispersion. The signal
probably combines contributions of more than one sulfur $3p$ band. A more
direct insight may require an experiment covering a broader $k_{z}$ momentum
range, using the incident photon energy extending into the soft X-ray region.

Motivated by findings related to Fig. 3 (a) and to further examine the
$\beta$-feature, we performed another Fermi-level scan using a photon energy
of 58 eV. The result of the scan is shown in Fig. 4. (a). Along with the
hexagonal hole pocket centered at the $\overline{\mathrm{\Gamma}}$ point, the
$\beta$-feature shows as smaller triangular electron pockets at each
$\overline{\mathrm{K}\ }$corner of the small Ni1/3NbS2 Brillouin zone. The
shallowness of these pockets can be appreciated upon inspecting Figs. 4 (b)
and (c). These panels show the ARPES intensity profiles along
$\overline{\mathrm{K}}$-$\overline{\mathrm{M}}$-$\overline{\mathrm{K}}$ and
${\overline{\mathrm{M}}}_{0}$-$\overline{\mathrm{\Gamma}}$-${\overline{\mathrm{M}}}_{0}$
directions, indicated by light green arrows in Fig. 4 (a).

The feature showing around $\overline{\mathrm{K}\ }$points in Fig. 4 (a) is
somewhat weaker and less structured than the one recently found in Co1/3NbS2
[21, 22, 23]. Nevertheless, it shares several qualities with the features
observed recently in Co1/3NbS2 and Cr1/3NbS2 [26, 27, 21, 22, 23, 28].
Centered around the same points in the $(k_{x},k_{y})$ plane, the
$\beta$-feature has always been found as a relatively weak signal, confined
around the Fermi level and occupying a very narrow energy window relative to
other bands. These similarities appear despite the three materials having
different magnetic ordering types and significant electronic structure
differences, some of which are to be addressed further below. Our finding
indeed supports the possible universality of the $\beta$-feature in magnetic
intercalates of TMDs, as suggested in Ref. 23. As elaborated there, this
universality may stem from strong electronic correlations produced by the
coupling between magnetic ions and the metallic subsystem. The strong electron
correlations are expected to arise primarily because the low energy
excitations in magnetic sublayers tend to dynamically reconfigure the charge
and spin exchange between metallic subsystems [48]. In turn, and on a slower
scale, the itinerant electrons can affect the state of the magnetic subsystem,
thus introducing the "long-term memory" effects for the itinerant subsystem.
It should be noted, however, that the recent ARPES study of the M1/3NbS2 (M=
Fe, Co, Ni) series, using $h\nu$ = 120 eV, reveals the $\beta$-feature around
$\overline{\mathrm{K}}$-point in Co1/3NbS2 but not in Fe1/3NbS2 and Ni1/3NbS2
[21]. Part of the reason for the latter discrepancy from our results may lie
in the values of the incident photon energy used. The earlier study in
Co1/3NbS2 [23] and here reported dependence of the signal intensity on photon
energy speak in favor of such a possibility.

In principle, the dependence of ARPES spectra on photon energy offers
information on band dispersion along the kz direction. In Co1/3NbS2 the kz
band dispersion has been reported for the Co-derived band near the Fermi level
[22, 21]. The kz dependence of the $\beta$-feature may be considered a
possible reason for it being seen only within a particular photon energy
range. However, other factors can also affect the ARPES signal, such as the
photon energy-dependent transition matrix element influenced by the symmetry
of the underlying electronic wavefunctions. The $k_{z}$ range covered in our
experiment by varying the photon energy exceeds the size of a single Brillouin
zone, as evident from Fig. 3, which includes several replicas of the
$\mathrm{\Gamma}$ point along the $k_{z}$-axis direction. More precisely, the
$9^{th}$ replica of the first Brillouin zone in $k_{z}$-direction extends
between points $\mathrm{A}_{8}$ and $\mathrm{A}_{9}$ in Fig. 3, with the point
$\mathrm{\Gamma}_{9}$ in its center. Likewise, the $8^{th}$ and $10^{th}$
replicas are centered around $\mathrm{\Gamma}_{8}$ and $\mathrm{\Gamma}_{10}$.

The $\beta$-feature in Fig. 3(a) appears around ${\mathrm{\Gamma}}_{9}$, but
does not show around ${\mathrm{\Gamma}}_{8}$ and ${\mathrm{\Gamma}}_{10}$.
Accordingly, the mere $k_{z}$ dispersion is probably not the main reason for
the $\beta$-feature showing in a particular photon energy range. Further ARPES
experiments over an extended photon energy range, possibly using soft X-rays,
may help to resolve the question entirely.

Figure 4: (Color online) (a) Fermi surface map of Ni1/3NbS2 measured with
h$\nu$ = 58 eV. The pockets centered around $\overline{\mathrm{K}}$ and
equivalent points in the Brillouin zone appear well pronounced here,
witnessing the $\beta$-feature in Ni-intercalate. The ARPES intensity plots
along (b)
$\overline{\mathrm{K}}$-$\overline{\mathrm{M}}$-$\overline{\mathrm{K}}$ and
(c)
${\overline{M}}_{0}$-$\overline{\mathrm{\Gamma}}$-${\overline{\mathrm{M}}}_{0}$
direction, indicated by light green allows in panel (a).

##  V Qualitative differences in couplings between different magnetic ions and
TMD layers

The second rather significant difference between the electronic structures of
2H-NbS2 and Co1/3NbS2 is related to Nb 4d conduction bands [23]. A much larger
splitting in energy between Nb 4d bonding and anti-bonding bands was found in
the central part of the Brillouin zone in Co1/3NbS2 than in 2H-NbS2. In
2H-NbS2, both bands cross the Fermi level upon approaching the
$\mathrm{\Gamma}$ point. In Co1/3NbS2, only the anti-bonding band shows such
crossing, while the bonding band is found submerged 0.3 eV below the Fermi
level at the $\mathrm{\Gamma}$ point. The latter was primarily attributed to
strongly amplified hybridization between NbS2 layers, with the intercalated Co
ions serving as hybridization bridges [23]. The electron transfer from Co ions
into NbS2 layers was found to play a minor role in positioning the bonding
band in Co1/3NbS2 below the Fermi level. No such splitting caused by
intercalation-amplified hybridization between NbS2 layers is visible in
Ni1/3NbS2. No Nb 4d bonding band is identified below the Fermi level at the
$\mathrm{\Gamma}$ point. The splitting between conduction bands in Ni1/3NbS2
is minimal at their crossing of the Fermi level. On the other hand, the
spotting between Nb-bands gets very pronounced upon approaching the zone
boundary, as evidenced by Fig. 1(a).

The DFT calculations agree with such a finding: The calculated electronic
spectra show strong intercalation-induced splitting between the bonding and
the anti-bonding bands in Co1/3NbS2 [23], and very small splitting in the
vicinity of the $\mathrm{\Gamma}$ point in Ni1/3NbS2 (see. Fig.1(a)). This
large difference between the two compounds is surprising since the Ni-Nb and
Co-Nb hybridizations are expected to be similar in size, based on listed Ni
and Co ionic radii [49] and the similar c-axis lattice constants found in two
compounds. The reasons for the discrepancy between the two systems can be
explored by examining the Wannier90 tight-binding parametrization of the DFT-
calculated spectra in terms of maximally localized orbitals and the
hybridization integrals between them. The Wannier90 parametrization in
Ni1/3NbS2 produces the "maximally localized" Wannier functions that resemble
the Ni 3$d$ and S 3$p$ atomic orbitals. As for the Nb-dominated conduction
bands, one Wannier function per Nb site is produced. These Wannier functions
are dominantly composed of Nb 4$d$ atomic orbitals, with three neighboring Nb
ions equally contributing to each Wannier function. These effective three-
legged, Y-shaped "Nb orbitals" are centered between three contributing Nb ions
and extend toward their nuclei, as shown in Fig. 5(c). For simplicity, in the
remaining, we refer to them as "Nb Wannier orbitals", "Nb($Y$) orbitals", or
just as the "Nb orbitals". Note that the conduction bands in any 2H-NbS2-based
systems are dominantly composed of such Nb($Y$) orbitals.

The participation of Ni ions in building the conduction bands comes through
two minority-spin 3$d$ $e_{g}$ orbitals, $d_{1}$ and $d_{2}$. These two
orbitals are degenerate in energy, with their energy levels positioned less
than 1 eV above the energy levels of the Nb orbitals and the Fermi level. To a
large degree, $d_{1}$ and $d_{2}$ resemble the standard 3$d_{z^{2}_{1}}$ and
3$d_{x^{2}_{1}-y^{2}_{1}}$ orbitals, with the coordinate system
$(x_{1},y_{1},z_{1})$ set by the ideal octahedron with the sulfur atoms at its
vertices. In essence, we find that the effective hybridization between any of
the orbitals $d_{1}$ and $d_{2}$ and the Nb orbitals is comparable in size to
analogous hybridization integrals found in Co1/3NbS2 [23]. However, contrary
to the Co system, the effective hybridization integrals between a Ni orbital
and the Nb orbitals within a given Nb-plane vary in size and sign, as
illustrated in Fig. 5. Panels (a) and (b) in Fig. 5 show the hybridization
integrals ${t_{\mathrm{Ni}(d_{1,2})-\mathrm{Nb(\mathit{Y})}_{j}}}$ between the
Ni $3d$ $e_{g}$ orbitals $d_{1}$ and $d_{2}$, and the Nb($Y$) orbitals at
position $j$ in the Nb layer below the Ni ion. The sizes and signs of the
hybridization integrals are shown through the areas of the circular surfaces
and the colors of the disks centered at the positions of the Nb($Y$) orbitals.
Only the hybridization integrals for the three nearest-neighboring and the
three next-to-nearest-neighboring Nb orbitals are shown. The hybridization
integrals to further Nb orbitals are very small in comparison. It may be noted
that the sign of the hybridization integral depends on the position of the Nb
orbital. Also, the absolute values of hybridization integrals for next-to-
nearest-neighbors tend to be much bigger than for the nearest neighbors. Most
importantly, it turns out that the sum of the hybridization integrals between
any of the two Ni $3d$ $e_{g}$ orbitals and all the Nb Wannier orbitals
Nb($Y$)j of a particular Nb plane,
$\sum_{j}{t_{\mathrm{Ni}(d_{1,2})-\mathrm{Nb(\mathit{Y})}_{j}}}$, is
numerically indistinguishable from zero. This vanishingly small sum is not
accidental. Instead, it is the consequence of the exact symmetry present in
the system. The claim can be exactly proved by constructing two linear
combinations of two wavefunctions, $d_{o1,o2}\propto d_{1}\pm id_{2}$, with
$d_{o1}$ being shown in Fig. 5(c). By construction (see also Appendix C),
these orbitals conform to the crystal symmetry, contrary to the orbitals
$d_{1}$ and $d_{2}$. Upon their rotation by $2\pi/3$ around the $c$-axis of
the crystal, only their complex phase changes by $\pm 2\pi/3$. The
corresponding sets of hybridization integrals possess the same symmetry: The
three hybridization integrals to Nb orbitals connected by the $\pm 2\pi/3$
rotation around the $c$-axis passing through the Ni-site have the same
absolute value. Only the complex phase of each hybridization integral is
shifted by $\pm 2\pi/3$ relative to the other two. Consequently, their partial
sum is precisely zero. The situation is illustrated in Fig. 5(d), where six
hybridization integrals are positioned in the complex plane, and within the
complex phase color wheel, forming two equilateral triangles. The symmetry
argument appears in further detail in Appendix C.

The vanishing sum of the hybridization integrals implies the zero value of the
$(k_{x},k_{y})=0$ Fourier transform of the Ni-Nb hybridization. Consequently,
it also implies zero coupling between any of two Ni 3d $e_{g}$ orbitals and
the $(k_{x},k_{y})=0$ states of the Nb-layers. The Fourier transform of the
hybridization integral and the coupling to Nb states remains finite for non-
zero $(k_{x},k_{y})$. As a result, the splitting between the bonding and the
anti-bonding bands produced by Ni intercalation vanishes at the
$\mathrm{\Gamma}$ point and remains small in the central part of the first
Brillouin zone. It can be concluded that the striking difference in how Ni and
Co ions hybridize/interact with the NbS2 layers is mainly related to their
spin and charge states acquired in similar local environments. The respective
"starting/ideal" spin/charge states are Ni+2 3$d^{8}$ and Co+2 3$d^{7}$.
Different symmetries of hybridization configurations in Ni and Co compounds
are set entirely by these "ideal" Ni/Co states. The substantial hybridization
amplitude between Ni $3d$ $e_{g}$ orbitals and the conduction band Nb orbitals
leads to Ni charge and spin state deviations from the ideal Ni+2 3$d^{8}$
state, but does affect the established symmetry. These further deviations in
charge are discussed in more detail in Appendix B.

Figure 5: (Color online) (a) The crystal structure of 2H-NbS2 as viewed along
the $c$-axis direction. Grey balls represent the Nb atoms, and the sulfur
atoms are shown in green. The three-legged Nb Wannier orbitals, each joining
three neighboring Nb sites, are shown in dark-red (positive sign) and blue
(negative sign). The red circle in the center marks the position of the
magnetic Ni ion above the NbS2 layer in Ni1/3NbS2. (b), (c) The Ni $3d$
$e_{g}$ orbitals with energy above the Fermi level. The orbitals are shown in
the usual way, with their absolute value being proportional to the distance of
the surface point from the center of the Ni atom, and the sign is represented
through color (red for +1, cyan for -1). The cylinders below each orbital
represent the hybridization integrals between the orbital and the six Nb($Y$)
Wannier orbitals, labeled by 1 to 6 in (a), in the Nb plane below the Ni ion.
The cylinders are positioned at the centers of the Nb Wannier orbitals. The
sizes and colors of the cylinders represent the sizes and signs of the
hybridization integrals. The numbers next to some cylinders illustrate the
absolute values of respective hybridization integrals in eV. Note the
variation in absolute values and signs. (d) The linear combination,
$d_{o1}\sim d_{1}+id_{2}$, of the two $e_{g}$ orbitals. The complex phase of
the wavefunction is coded in color. The complex phase color wheel is shown in
(e). Contrary to $d_{1}$ and $d_{2}$, $d_{o1}$ is the eigenfunction of the
$2\pi/3$ crystal symmetry rotation. The absolute values of hybridization
integrals to Nb($Y$) orbitals at positions equivalent by crystal symmetry are
the same in this representation, whereas their phases (the same color coding
applies) differ by $\pm 2\pi/3$. (e) The same six cylinders (disks) are
positioned within the complex phase color wheel, showing the equilateral
triangular arrangement.

##  VI Conclusion

Using various photon energies, we have performed ARPES measurements in
Ni1/3NbS2. The calculated electronic spectra provide good guidelines for the
analysis of our experimental findings. The main elements of the observed band
dispersion of Ni1/3NbS2 can be traced back to those of the host 2H-NbS2. The
electron transfer from Ni atoms to NbS2 layers is mainly responsible for the
observed Fermi level shift within the 2H-NbS2 conduction bands. However, the
2H-NbS2 rigid band picture cannot describe the electronic structure of
Ni1/3NbS2. Interestingly, the deviations from the rigid band picture
qualitatively differ from those recently reported for Co1/3NbS2 [23]. In
contrast to the cobalt case, the intercalation by nickel does not produce a
strong splitting of the Nb 4$d$ conduction bands in the center of the
Brillouin zone. However, we argue that both behaviors are consistent with
strong hybridization between 3$d$ orbitals of the intercalated ion and Nb
metallic states. We find that the hybridization between individual 3$d$
orbitals and neighboring Nb orbitals is comparable in size in both Ni and Co
systems. What differentiates the two systems are the symmetries of the
relevant 3$d$ states, the ones closest to the Fermi level. In the cobalt case
the strong hybridization is reflected in the large conduction band splitting
around the Brillouin zone center. Contrary to that, in the nickel case, the
large conduction bands splitting caused by Ni orbitals occurs far from the
center of the Brillouin zone but turns out exactly zero at the center of the
Brillouin zone, and remains small around the zone center. The latter behavior
seen in our ARPES data is understood through our DFT calculations and the
symmetry analysis.

Regarding the newly discovered $\beta$-feature in Ni1/3NbS2, the pronounced
differences between Ni1/3NbS2, Co1/3NbS2, and Cr1/3NbS2 speak in favor of the
anticipated universality of the $\beta$-feature [23]. Again, as in the case of
Co intercalated compound, the $\beta$-feature is particularly pronounced
around $\overline{\mathrm{K}}$-points of the Brillouin zone of the Ni
intercalated material. Nothing similar is expected from DFT calculations
performed for 2H-NbS2 or Ni1/3NbS2. The $\beta$-feature signal’s strength in
Ni1/3NbS2 is found to depend strongly on the photon energy used in ARPES,
pointing to the structure of its wavefunction along the $c$-axis, or its
energy dispersion along the same direction. The previously proposed [23] role
of strong electron correlations in creating the $\beta$-feature still remains
to be confirmed.

## Acknowledgements

The Croatian Science Foundation has partly supported this work under the
project numbers IP 2020-02-9666, UIP 2019-04-2154, and by the project
Cryogenic Centre at the Institute of Physics – KaCIF co-financed by the
Croatian Government and the European Union through the European Regional
Development Fund-Competitiveness and Cohesion Operational Programme (Grant No.
KK.01.1.1.02.0012). The research has been performed at the National
Synchrotron Radiation Centre SOLARIS at the UARPES beamline (proposal No.
201046). This publication was partially developed under the provision of the
Polish Ministry and Higher Education project Support for research and
development with the use of research infra-structure of the National
Synchrotron Radiation Centre SOLARIS under contract nr 1/SOL/2021/2. The work
at the TU Wien was supported by the European Research Council (ERC
Consolidator Grant No 725521). The work at the University of Zagreb was
supported by the Croatian Science Foundation under the project IP 2018 01 8912
and CeNIKS project co-financed by the Croatian Government and the European
Union through the European Regional Development Fund - Competitiveness and
Cohesion Operational Programme (Grant No. KK.01.1.1.02.0013). The work at AGH
University was supported by the National Science Centre, Poland, Grant OPUS:
2021/41/B/ST3/03454; the Polish National Agency for Academic Exchange under
“Polish Returns 2019” Programme No. PPN/PPO/2019/1/00014; and the “Excellence
Initiative–Research University” program for AGH University of Krakow. I.B.
acknowledges support from the Swiss Confederation through Government
Excellence Scholarships. E.T. acknowledges the support of the Croatian Science
Foundation Project IP-2022-10-9423.

## Appendix A. Unfolding from the Wannier90 representation

The superstructure introduced by intercalation leads to the small first
Brillouin zone that corresponds to the supercell in real space. The DFT
calculations can provide the energy eigenvalues $E_{n,k}$ for the Kohn-Sham
band $n$ and any wave vector $k$ within the small first Brillouin zone. The
ARPES spectra, however, are not confined to the wave vectors of the small
Brillouin zone, or any Brillouin zone for that matter, but address the whole
reciprocal space. The DFT results being "folded" into a small part of the
reciprocal space, makes the comparison with measured spectra very difficult.
The first step in comparing the DFT-calculated results with the spectra
measured in ARPES goes in the direction of projecting the wavefunctions
$\psi_{n,k}$ onto plain waves $e^{ik^{\prime}r}$ at arbitrary wavevector
$k^{\prime}$.

The approach to "unfolding" that we use in the paper uses the Wannier90
parameterization of the DFT results [42], relying on a handful of maximally
localized Wannier wavefunctions ("orbitals") $\phi_{m}$ in the unit cell. The
orbitals are usually centered around particular atoms and often acquire the
form of the textbook atomic orbitals. Within the Wannier90 approach, the Kohn-
Sham wave function $\psi_{n,k}$ follows the tight binding form for the crystal
with several atoms per unit cell.

$\psi_{n,k}(r)=\frac{1}{\sqrt{N}}\sum_{R,m}e^{ikR}\alpha_{n,k}^{(m)}\phi_{m}(r-R-
d_{m}).$ (1)

Here $R$ stands for the vector of the Bravais lattice of the crystal, running
over the positions of its unit cells. $N$ stands for the number of unit cells
in the large-yet-finite crystal. The Wannier90 approach also provides the
energy levels of the orbitals, $E_{m}$ as well as the overlap integrals
$t_{m,m^{\prime}}(R)$ between orbitals $m$ and $m^{\prime}$ within the same
unit cell ($R=0$), or within the unit cells separated by $R$. Within the usual
scheme of the tight-binding approximation, these parameters lead to the energy
bands $E_{n,k}$ that reproduce the DFT- calculated spectra very well, and
provide the orbital amplitudes $\alpha_{n,k}^{(m)}$ that enter the wave
function $\psi_{n,k}$ in Eq.(1). The projection of the wave function
$\psi_{n,k}$ on the plane wave $\left\langle
r\middle|k^{\prime}\right\rangle=\frac{1}{\sqrt{V}}e^{Ik^{\prime}r}$, is easy
to calculate, yielding

$\left\langle
k^{\prime}\middle|\psi_{n,k}\right\rangle=\sum_{m,G}\delta_{k^{\prime},k+G}\alpha_{n,k}^{(m)}g_{m,k^{\prime}}e^{ik^{\prime}d_{m}}$
(2)

where $g_{m,k^{\prime}}$ stands for the projection of the local orbital on the
plane wave,

$g_{m,k^{\prime}}\equiv\frac{1}{\sqrt{v_{1}}}\int
dre^{ik^{\prime}r}\phi_{m}(r).$ (3)

and $G$ is the vector of the crystal’s reciprocal lattice and $v_{1}$ denotes
the volume of its unit cell. The absolute value of the projection determines
the intensity of the unfolded spectrum at point $(k^{\prime},E=E_{n,k})$,

$w(k^{\prime},E)=\sum_{n,k}\left|\left\langle
k^{\prime}\middle|\psi_{n,k}\right\rangle\right|^{2}\delta(E_{n,k}-E),$ (4)

or,

$w(k^{\prime},E)=\sum_{n,k}\left|\sum{m}\alpha_{n,k}^{(m)}g_{m,k^{\prime}}e^{ik^{\prime}d_{m}}\right|^{2}\delta(E_{n,k}-E).$
(5)

Here $k^{\prime}$ stands for whatever vector in the reciprocal space, whereas
$k=k^{\prime}-G$ stands for its counterpart in the first Brillouin zone of the
crystal. In order to make the variation of the weights visible in the image
plot, one needs to replace the function $\delta(E)$ by the test function
$\gamma(E,\eta)$ of a finite width $\eta$. In this paper we choose to use the
Gaussian function for $\gamma(E,\eta)$ ,
$\gamma(E,\eta)=\frac{1}{{\gamma\sqrt{2\pi}}}e^{-E^{2}/2\gamma}$, and
$\eta=0.03$ eV. The advantage of the unfolding being done through the
Wannier90 parametrization of the spectra, allows us to get additional insights
by making the various choices for $g_{m}(m,k^{\prime})$, instead of the one
prescribed by Eq.(3). First, it should be remarked that the spectrum
$w(k^{\prime},E)$ as defined in Eq.(5) is not expected to simulate the ARPES
spectra properly. In the simplest three-step approximation to ARPES [50], the
ARPES matrix elements include the polarization wave vector $\epsilon$ of the
photon beam, leading to

$g_{m,k^{\prime},\epsilon}\equiv\frac{1}{\sqrt{v_{1}}}\int
dre^{ik^{\prime}r}(\epsilon\cdot r)\phi_{m}(r)$ (6)

to appear in Eq.(5) instead of $g_{m,k^{\prime}}$ [50]. Second, the effect of
the crystal surface is to be taken into account, leading to further
modification of the spectrum [46, 50, 51]. In fact, the simplest way to
compare the DFT-calculated and ARPES spectra is to show all calculated bands
in their unfolded variant, without some bands getting suppressed relative to
others due to the particular symmetry or the size of the orbitals involved, or
the photon polarization. This choice is the closest one to the usual plots of
bands dispersion in the extended zone in the $k$-space, however with weights
of particular points being modulated by the variations of atomic positions,
orbital energy levels $E_{n}$ and the orbital overlap integrals
$t_{m,m^{\prime}}(R)$ being modulated by the superstructure. This choice is
made in preparing the DFT-unfolded spectra in Fig. 1 (a). It corresponds to
taking $g_{m,k^{\prime}}=1$ in Eq.(5). Another way to profit from a particular
choice of $g_{m}$’s, is to explore the contributions of particular
orbitals/atoms by choosing $g_{m}=0$ for all the other orbitals. We use this
possibility to focus on the contribution of the Nb orbitals to bands around
the Fermi level. In Fig. 6 we show and compare the DFT+unfolded spectra in
Ni1/3NbS2 and Co1/3NbS2 featuring only the contributions from Nb orbitals.
Notably, contrary to the Co1/3NbS2 case, the Nb contribution close below the
Fermi level is not observed around the $\Gamma$ point in Ni1/3NbS2. Another
way to plot the contribution of various orbitals to the unfolded spectra is
provided in Supplemental Material, Fig. S2 (b).

Figure 6: (Color online) The unfolded spectra of the DFT-calculated results
for Ni (a) and Co (b) intercalated systems, featuring only the contribution of
Nb atoms. The great difference between the two figures is the result of very
different ways in which Ni and Co ions connect the neighboring Nb layers,
discussed in the main text. In particular, the big separation in energy
between bands is visible at the $\Gamma$ point in the Co-based system, whereas
almost no separation is found in the Ni-based system. The qualitative
differences are also noticeable throughout the (large) Brillouin zone of
2H-NbS2.

## Appendix B. Bands filling

Here, we discuss the effects of the intercalation on Ni charge and the
conduction bands filling. The direct approach, based on DFT-calculated
"projected density of states" (PDOS), or the approach starting from the
Wannier90 parametrization of the DFT spectra, may be used, giving the same
results.

The integration of the density of states (DOS) obtained from the Quantum
Espresso DFT calculation provides the value of the conduction bands filling
factor of precisely 5/6. This result is evident in the rigid-band picture,
assuming that each Ni atom donates $q=2$ electrons into NbS2 conduction band,
resulting in the Ni+2 ion. For the $c_{i}=1/3$ concentration of intercalated
atoms per NbS2 formula unit, the conduction bands filling changes from
$f_{0}=1/2$ to $f=f_{0}+q\times c_{i}\times{1/2}=5/6$. The factor $1/2$
appears as the consequence of the spin degeneracy of band states. In our case,
the two electrons to be transferred to NbS2 conduction bands reside in Ni 4$s$
orbitals of the standalone magnetic Ni atoms. The transfer occurs upon
positioning these atoms at the centers of sulfur octahedral voids.
Simultaneously, the lowest Ni orbitals above the Fermi level become the local
minority-spin Ni $3d$ $e_{g}$ orbitals, $d_{1}\approx d_{z_{1}^{2}}$ and
$d_{2}\approx d_{x_{1}^{2}-y_{1}^{2}}$. The labels correspond to the local
coordinate systems set by the sulfur octahedron (see also Fig. 5 and Appendix
B). The final ingredient to account for is the hybridization between the
intercalated atom’s orbitals and the orbitals of the host material. It may not
be evident that hybridization invalidates the above relation between the
charge of the intercalated ion and the conduction band filling. The
hybridization does not affect the band filling, as the number of states in
conduction bands and the number of electrons hosted by those bands remain
unchanged. Conversely, the charge of the intercalated atom does change as its
orbitals participate in the filled conduction bands states. Our DFT
calculation gives the value around 0.37 for the occupancy per Ni ion of its
3$d_{1}$ and 3$d_{2}$ orbitals. This calculation is done by integrating the
partial density of states, derived by projecting the conduction bands
wavefunctions to these atomic states (the Löwdin method [52]). As already
stated in the main text, the standard QE calculation of ionic charges,
performed by integrating the spatial charge density distribution within the
non-overlapping spheres around atomic nuclei, provides a compatible value of
$+1.6(\approx 2-0.37)$ $e$ for the Ni ionic charge 222A more elaborate
discussion of the charge allocation to atoms within Quantum Espresso, using
the Löwdin method, would also address the charge density and wavefunctions
further away from the Ni nucleus. The related calculations tend to attribute
some charge to Ni 4$s$ and Ni 4$p$ orbitals, with energy levels further above
the Fermi level, through their hybridization with the states/bands deep below
the conduction bands, presumably sulfur-dominated. Similarly, the
hybridization of 3$d_{x_{1}^{2}-y_{1}^{2}}$ and 3$d_{z_{1}^{2}}$ orbitals with
the states deep below the conduction band also provides some additional charge
in these orbitals. We consider such details to be out of the scope of the
present paper. The message of this appendix is illustrated in Fig. 7.

Figure 7: (Color online) The schematic view of Ni and NbS2 subsystems. The Ni
3$d$ orbital levels are marked in red. The electrons and the filled portion in
conduction bands are painted green. The NbS2 empty states are shown in light
blue. (a) The standalone magnetic Ni atom and NbS2 crystal. (b) The transfer
of 2 electrons from one Ni atom per 3 NbS2 formula units, filling the
conduction band to 5/6. (c) The hybridization mixes the Ni orbitals and the
NbS2 conduction bands states. In turn, the NbS2 orbitals participate in the
Ni-dominated bands. The number of states in the conduction bands and the
number of electrons in those bands do not change. The band filling remains
5/6. The charge of the Ni ion profits from participating in occupied
conduction band states by receiving extra $\delta$ electrons per Ni atom. The
charge of the Ni ion is $+(2-\delta)$.

## Appendix C. The hybridization of Ni orbitals with the Nb layers

The Wannier90 parametrization of the DFT-calculated wavefunctions and
electronic spectra in Co1/3NbS2 reveals only one relevant Co 3d orbital close
to the Fermi level that effectively (via intermediate sulfur orbitals)
hybridizes with the Nb($Y$) orbitals [23]. This orbital is approximately
described by the Co 3$d_{z^{2}}$ electronic state, oriented along the
crystal’s $c$-axis. The Co-Nb hybridization integral maximizes at
approximately 0.24 eV for Nb($Y$)-orbitals nearest to the Co ion. The total
hybridization,
$\sum_{j}{t_{\mathrm{Co}(d_{z^{2}})-\mathrm{Nb(\mathit{Y})}_{j}}}$, between
the relevant Co orbital and all the Nb orbitals of the closest Nb layer
amounts to approximately three times that value. In Ni1/3NbS2, the similar
Wannier90 parametrization of the bands around the Fermi level produces two
maximally-localized, real-numbered Ni orbitals strongly hybridizing with the
Nb($Y$) orbitals. The orbitals $d_{1}$ and $d_{2}$ are the usual $e_{g}$ pair
of states emerging in the octahedral environment provided by the sulfur ions.
As expected, they appear degenerate in energy for the ideal octahedron or the
octahedron dilated along the $c$-axis, perpendicular to its two faces, as
found in our crystal. The resemblance between orbitals $d_{1}$ and $d_{2}$ and
their counterparts within the ideal octahedron is very pronounced,
$d_{1}\approx d_{z_{1}^{2}}$ and $d_{2}\approx d_{x_{1}^{2}-y_{1}^{2}}$. The
coordinate system $(x_{1},y_{1},z_{1})$ implied here points along the axes of
the imaginary "ideal sulfur octahedron".

The calculated Ni($d_{1}$)-Nb($Y$)j and Ni($d_{2}$)-Nb($Y$)j hybridization
integrals maximize around 0.2 eV in absolute value for specific Nb($Y$)
orbitals. However, for both Ni orbitals, the hybridization integrals appear
with positive and negative signs. Their sum, per single Ni orbital and single
Nb layer, turns out numerically indistinguishable from zero within the
numerical accuracy of the DFT+Wannier90 calculation. The vanishing sum implies
zero coupling between any of two Ni 3d orbitals and the $(k_{x},k_{y})=0$
states of the Nb-layers. Consequently, the splitting between bonding and anti-
bonding bands produced by Ni intercalation vanishes at the $\mathrm{\Gamma}$
point and remains very small in the central part of the Brillouin zone.

It should also be remarked that the sums of absolute values of hybridization
integrals, per single Ni orbital and all Nb orbitals of the nearby layer, are
only some 40% smaller in Ni1/3NbS2 than in Co1/3NbS2. Certain Ni-Nb
hybridization integrals are indeed strongly pronounced, partly justifying the
naive preliminary expectation regarding the comparable hybridization strengths
in Co and Ni intercalates.

The strictly vanishing sum of Ni-Nb hybridization integrals between relevant
Ni orbitals and the metallic layer is not likely to be a mere numerical
coincidence. As we show now, it appears as the consequence of the symmetry
that involves the $3d^{8}$ electronic configuration of the magnetic ion in the
octahedral environment and the $2\pi/3$ rotation symmetry of the surrounding
crystal.

The reasons for the vanishing average of the hybridization to nearby layers
are not evident for Wannier90-supplied maximally localized orbitals. It
becomes more evident upon choosing the equivalent orbitals, $d_{o1}$ and
$d_{o2}$, fully obeying the symmetry of the Ni ion environment, although
slightly more extended in space. The crystal symmetry operator,
$\mathcal{R}\mathrm{\equiv}\mathrm{Rot[}c,n\mathrm{\times}\mathrm{\
2}\pi\mathrm{/}$3]$\mathrm{,\ }$ denoting the $2\pi\mathrm{/3}$ rotation
around the $c$-axis of the crystal is essential in that respect. Indeed, the
wavefunctions ${\mathcal{R}}^{n}d_{1}$ and ${\mathcal{R}}^{n}d_{2}$ ($n$
denoting an integer) represent the choices equivalent to $d_{1}$ and
$d_{2}\mathrm{\ }$for electronic structure description. The same applies to
the two linear combinations,

$\displaystyle d_{o1}\propto\sum_{n=0}^{2}{e^{+in\frac{2\pi}{3}}\
\mathcal{R}^{n}\ d_{1}},$ $\displaystyle
d_{o2}\propto\sum_{n=0}^{2}{e^{-in\frac{2\pi}{3}}\ \mathcal{R}^{n}\ d_{2}},$

which are the eigenfunctions of the $\mathcal{R}$ operator by construction,
and mutually orthogonal. Upon normalization, and for the ideal octahedron,
these two orbitals can be written as

$\displaystyle d_{o1}=\frac{1}{\sqrt{2}}\left(d_{1}+id_{2}\right),$
$\displaystyle d_{o2}=\frac{1}{\sqrt{2}}\left(d_{1}-id_{2}\right).$

The choice of $d_{o1}$ and $d_{o2}$ for the Ni orbitals restores the crystal
symmetry of the electronic structure parametrization, missing in the Wannier90
maximally localized approach.

Finally, we turn to $\mathrm{Ni}\left(d_{o1}\right)-\mathrm{Nb}(Y)_{j}$ and
$\mathrm{Ni}\left(d_{o2}\right)-\mathrm{Nb}(Y)_{j}$ hybridization integrals:
For a given Ni ion, these can be grouped into groups of three, each group
containing Nb($Y$) orbitals connected by $\mathcal{R}$ rotations around the Ni
ion. Through crystal symmetry and our construction of $d_{o1}$ and $d_{o2}$,
the related hybridization integrals are equal in magnitude and bear the phases
proportional to unity, $e^{+i2\pi/3}$ and $e^{-i2\pi/3}$. Thus, the sum of
hybridization integrals exactly vanishes for each group, rendering to zero the
sum of the hybridization integrals between $d_{o1}$ ($d_{o2}$) orbital and any
Nb layer. However, as the hybridization varies from one Nb orbital position to
another, the Fourier transform of the hybridization integral remains finite at
finite $(k_{x},k_{y})$, being precisely zero only at $(k_{x},k_{y})=0$ point.

The illustration that involves the Ni ($d_{o1}$ orbital and six orbital
positions in one Nb layer is provided in Fig. 5. Further illustrations can be
found in the Supplemental Material.

## References

  * Steglich  _et al._ [1979] F. Steglich, J. Aarts, C. D. Bredl, W. Lieke, D. Meschede, W. Franz, and H. Schäfer, Superconductivity in the presence of strong Pauli paramagnetism: CeCu2Si2, Physical Review Letters 43, 1892 (1979).
  * Bednorz and Müller [1986] J. G. Bednorz and K. A. Müller, Possible highTc superconductivity in the Ba-La-Cu-O system, Zeitschrift für Physik B Condensed Matter 64, 189 (1986).
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, From quantum matter to high-temperature superconductivity in copper oxides, Nature 518, 179 (2015).
  * Kamihara  _et al._ [2006] Y. Kamihara, H. Hiramatsu, M. Hirano, R. Kawamura, H. Yanagi, T. Kamiya, and H. Hosono, Iron-based layered superconductor: LaOFeP, Journal of the American Chemical Society 128, 10012 (2006), pMID: 16881620, https://doi.org/10.1021/ja063355c .
  * de la Cruz  _et al._ [2008] C. de la Cruz, Q. Huang, J. W. Lynn, J. Li, W. R. II, J. L. Zarestky, H. A. Mook, G. F. Chen, J. L. Luo, and N. L. Wang, Magnetic order close to superconductivity in the iron-based layered LaO1-xFxFeAs systems, Nature 453, 899 (2008).
  * Barišić and Sunko [2022] N. Barišić and D. K. Sunko, High-Tc cuprates: a story of two electronic subsystems, Journal of Superconductivity and Novel Magnetism 35, 1781 (2022).
  * van den Berg and Cossee [1968] J. van den Berg and P. Cossee, Structural aspects and magnetic behaviour of NbS2 and TaS2 containing extra metal atoms of the first transition series, Inorganica Chimica Acta 2, 143 (1968).
  * Anzenhofer  _et al._ [1970] K. Anzenhofer, J. M. van den Berg, P. Cossee, and J. N. Helle, The crystal structure and magnetic susceptibilities of MnNb3S6, FeNb3S6, CoNb3S6 and NiNb3S6, Journal of Physics and Chemistry of Solids 31, 1057 (1970).
  * Friend  _et al._ [1977] R. H. Friend, A. R. Beal, and A. D. Yoffe, Electrical and magnetic properties of some first row transition metal intercalates of niobium disulphide, The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 35, 1269 (1977).
  * Parkin and Friend [1980a] S. S. P. Parkin and R. H. Friend, 3d transition-metal intercalates of the niobium and tantalum dichalcogenides. I. Magnetic properties, Philosophical Magazine B 41, 65 (1980a).
  * Parkin and Friend [1980b] S. S. P. Parkin and R. H. Friend, 3d transition-metal intercalates of the niobium and tantalum dichalcogenides. II. Transport properties, Philosophical Magazine B 41, 95 (1980b).
  * Parkin  _et al._ [1983] S. S. P. Parkin, E. A. Marseglia, and P. J. Brown, Magnetic structure of Co1/3NbS2 and Co1/3TaS2, Journal of Physics C: Solid State Physics 16, 2765 (1983).
  * Van Maaren and Schaeffer [1966] M. Van Maaren and G. Schaeffer, Superconductivity in group va dichalcogenides, Physics Letters 20, 131 (1966).
  * Wu  _et al._ [2022] S. Wu, Z. Xu, S. C. Haley, S. F. Weber, A. Acharya, E. Maniv, Y. Qiu, A. A. Aczel, N. S. Settineri, J. B. Neaton, J. G. Analytis, and R. J. Birgeneau, Highly tunable magnetic phases in transition-metal dichalcogenide ${\mathrm{Fe}}_{1/3+\delta}{\mathrm{NbS}}_{2}$, Phys. Rev. X 12, 021003 (2022).
  * Togawa  _et al._ [2012] Y. Togawa, T. Koyama, K. Takayanagi, S. Mori, Y. Kousaka, J. Akimitsu, S. Nishihara, K. Inoue, A. S. Ovchinnikov, and J. Kishine, Chiral magnetic soliton lattice on a chiral helimagnet, Physical Review Letters 108, 107202 (2012).
  * Dai  _et al._ [2019] Y. Dai, W. Liu, Y. Wang, J. Fan, L. Pi, L. Zhang, and Y. Zhang, Critical phenomenon and phase diagram of Mn-intercalated layered MnNb3S6, Journal of Physics: Condensed Matter 31, 195803 (2019).
  * Ghimire  _et al._ [2018] N. J. Ghimire, A. S. Botana, J. S. Jiang, J. Zhang, Y.-S. Chen, and J. F. Mitchell, Large anomalous Hall effect in the chiral-lattice antiferromagnet CoNb3S6, Nature Communications 9, 3280 (2018).
  * Šmejkal  _et al._ [2020] L. Šmejkal, R. González-Hernández, T. Jungwirth, and J. Sinova, Crystal time-reversal symmetry breaking and spontaneous hall effect in collinear antiferromagnets, Science Advances 6, eaaz8809 (2020), https://www.science.org/doi/pdf/10.1126/sciadv.aaz8809 .
  * Tenasini  _et al._ [2020] G. Tenasini, E. Martino, N. Ubrig, N. J. Ghimire, H. Berger, O. Zaharko, F. Wu, J. F. Mitchell, I. Martin, L. Forró, and A. F. Morpurgo, Giant anomalous hall effect in quasi-two-dimensional layered antiferromagnet ${\mathrm{Co}}_{1/3}{\mathrm{NbS}}_{2}$, Phys. Rev. Res. 2, 023051 (2020).
  * Mangelsen  _et al._ [2021] S. Mangelsen, P. Zimmer, C. Näther, S. Mankovsky, S. Polesya, H. Ebert, and W. Bensch, Interplay of sample composition and anomalous hall effect in CoxNbS2, Physical Review B 103, 184408 (2021).
  * Tanaka  _et al._ [2022] H. Tanaka, S. Okazaki, K. Kuroda, R. Noguchi, Y. Arai, S. Minami, S. Ideta, K. Tanaka, D. Lu, M. Hashimoto, V. Kandyba, M. Cattelan, A. Barinov, T. Muro, T. Sasagawa, and T. Kondo, Large anomalous hall effect induced by weak ferromagnetism in the noncentrosymmetric antiferromagnet ${\mathrm{CoNb}}_{3}{\mathrm{S}}_{6}$, Phys. Rev. B 105, L121102 (2022).
  * Yang  _et al._ [2022] X. P. Yang, H. LaBollita, Z.-J. Cheng, H. Bhandari, T. A. Cochran, J.-X. Yin, M. S. Hossain, I. Belopolski, Q. Zhang, Y. Jiang, N. Shumiya, D. Multer, M. Liskevich, D. A. Usanov, Y. Dang, V. N. Strocov, A. V. Davydov, N. J. Ghimire, A. S. Botana, and M. Z. Hasan, Visualizing the out-of-plane electronic dispersions in an intercalated transition metal dichalcogenide, Phys. Rev. B 105, L121107 (2022).
  * Popčević  _et al._ [2022] P. Popčević, Y. Utsumi, I. Biało, W. Tabis, M. A. Gala, M. Rosmus, J. J. Kolodziej, N. Tomaszewska, M. Garb, H. Berger, I. Batistić, N. Barišić, L. Forró, and E. Tutiš, Role of intercalated cobalt in the electronic structure of ${\mathrm{Co}}_{1/3}{\mathrm{NbS}}_{2}$, Phys. Rev. B 105, 155114 (2022).
  * Van Laar  _et al._ [1971] B. Van Laar, H. Rietveld, and D. Ijdo, Magnetic and crystallographic structures of MexNbS2 and MexTaS2, Journal of Solid State Chemistry 3, 154 (1971).
  * Polesya  _et al._ [2019] S. Polesya, S. Mankovsky, and H. Ebert, Electronic and magnetic properties of the 2H-NbS2 intercalated by 3d transition metal atoms, Zeitschrift für Naturforschung B 74, 91 (2019).
  * Sirica  _et al._ [2016] N. Sirica, S.-K. Mo, F. Bondino, I. Pis, S. Nappini, P. Vilmercati, J. Yi, Z. Gai, P. C. Snijders, P. K. Das, I. Vobornik, N. Ghimire, M. R. Koehler, L. Li, D. Sapkota, D. S. Parker, D. G. Mandrus, and N. Mannella, Electronic structure of the chiral helimagnet and $3d$-intercalated transition metal dichalcogenide Cr1/3NbS2, Physical Review B 94, 075141 (2016).
  * Sirica  _et al._ [2020] N. Sirica, P. Vilmercati, F. Bondino, I. Pis, S. Nappini, S. K. Mo, A. V. Fedorov, P. K. Das, I. Vobornik, J. Fujii, L. Li, D. Sapkota, D. S. Parker, D. G. Mandrus, and N. Mannella, The nature of ferromagnetism in the chiral helimagnet Cr1/3NbS2, Communications Physics 3, 1 (2020).
  * Qin  _et al._ [2022] N. Qin, C. Chen, S. Du, X. Du, X. Zhang, Z. Yin, J. Zhou, R. Xu, X. Gu, Q. Zhang, W. Zhao, Y. Li, S.-K. Mo, Z. Liu, S. Zhang, Y. Guo, P. Tang, Y. Chen, and L. Yang, Persistent exchange splitting in the chiral helimagnet ${\mathrm{Cr}}_{1/3}{\mathrm{NbS}}_{2}$, Phys. Rev. B 106, 035129 (2022).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with q uantum espresso, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Cococcioni and de Gironcoli [2005] M. Cococcioni and S. de Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the $\mathrm{LDA}+\mathrm{U}$ method, Phys. Rev. B 71, 035105 (2005).
  * Timrov  _et al._ [2018] I. Timrov, N. Marzari, and M. Cococcioni, Hubbard parameters from density-functional perturbation theory, Phys. Rev. B 98, 085127 (2018).
  * Corso [2014] A. D. Corso, Pseudopotentials periodic table: From H to Pu, Computational Materials Science 95, 337 (2014).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * Perdew  _et al._ [2009] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Erratum: Restoring the density-gradient expansion for exchange in solids and surfaces [phys. rev. lett. 100, 136406 (2008)], Phys. Rev. Lett. 102, 039902 (2009).
  * Marzari  _et al._ [1999] N. Marzari, D. Vanderbilt, A. DeVita, and M. C. Payne, Thermal contraction and disordering of the Al(110) surface, Physical Review Letters 82, 3296 (1999).
  * Battaglia  _et al._ [2007] C. Battaglia, H. Cercellier, L. Despont, C. Monney, M. Prester, H. Berger, L. Forró, M. G. Garnier, and P. Aebi, Non-uniform doping across the fermi surface of NbS2 intercalates, The European Physical Journal B 57, 385 (2007).
  * Note [1] The conventional illustration of the calculated band structure of Ni1/3NbS2 can be found in Supplemental Material.
  * Popescu and Zunger [2012] V. Popescu and A. Zunger, Extracting e versus k - effective band structure from supercell calculations on alloys and impurities, Physical Review B - Condensed Matter and Materials Physics 85, 085201 (2012).
  * Pacilè  _et al._ [2021] D. Pacilè, C. Cardoso, G. Avvisati, I. Vobornik, C. Mariani, D. A. Leon, P. Bonfà, D. Varsano, A. Ferretti, and M. G. Betti, Narrowing of d bands of feco layers intercalated under graphene, Applied Physics Letters 118, 121602 (2021).
  * [41] See also P. Bonfà, https://bitbucket.org/bonfus/unfold-x for the unfold-x code.
  * Pizzi  _et al._ [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, Wannier90 as a community code: new features and applications, Journal of Physics: Condensed Matter 32, 165902 (2020).
  * El Youbi  _et al._ [2021] Z. El Youbi, S. W. Jung, C. Richter, K. Hricovini, C. Cacho, and M. D. Watson, Fermiology and electron-phonon coupling in the 2H and 3R polytypes of NbS2 , Physical Review B 103, 155105 (2021).
  * Heil  _et al._ [2017] C. Heil, S. Poncé, H. Lambert, M. Schlipf, E. R. Margine, and F. Giustino, Origin of superconductivity and latent charge density wave in ${\mathrm{NbS}}_{2}$, Phys. Rev. Lett. 119, 087003 (2017).
  * Luttinger [1960] J. M. Luttinger, Fermi surface and some simple equilibrium properties of a system of interacting fermions, Phys. Rev. 119, 1153 (1960).
  * Strocov [2003] V. Strocov, Intrinsic accuracy in 3-dimensional photoemission band mapping, Journal of Electron Spectroscopy and Related Phenomena 130, 65 (2003).
  * Himpsel [1983] F. Himpsel, Angle-resolved measurements of the photoemission of electrons in the study of solids, Advances in Physics 32, 1 (1983).
  * Popčević  _et al._ [2023] P. Popčević, I. Batistić, A. Smontara, K. Velebit, J. Jaćimović, I. Živković, N. Tsyrulin, J. Piatek, H. Berger, A. Sidorenko, H. Rønnow, L. Forró, N. Barišić, and E. Tutiš, Electronic transport and magnetism in the alternating stack of metallic and highly frustrated magnetic layers in ${\mathrm{Co}}_{1/3}\mathrm{Nb}{\mathrm{S}}_{2}$, Phys. Rev. B 107, 235149 (2023).
  * Shannon [1976] R. D. Shannon, Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides, Acta Crystallographica Section A 32, 751 (1976).
  * Moser [2017] S. Moser, An experimentalist’s guide to the matrix element in angle resolved photoemission, Journal of Electron Spectroscopy and Related Phenomena 214, 29 (2017).
  * Strocov  _et al._ [2012] V. N. Strocov, M. Shi, M. Kobayashi, C. Monney, X. Wang, J. Krempasky, T. Schmitt, L. Patthey, H. Berger, and P. Blaha, Three-dimensional electron realm in ${\mathrm{VSe}}_{2}$ by soft-x-ray photoelectron spectroscopy: Origin of charge-density waves, Phys. Rev. Lett. 109, 086401 (2012).
  * Löwdin [1950] P. Löwdin, On the non-orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals, The Journal of Chemical Physics 18, 365 (1950), https://doi.org/10.1063/1.1747632 .
  * Note [2] A more elaborate discussion of the charge allocation to atoms within Quantum Espresso, using the Löwdin method, would also address the charge density and wavefunctions further away from the Ni nucleus. The related calculations tend to attribute some charge to Ni 4$s$ and Ni 4$p$ orbitals, with energy levels further above the Fermi level, through their hybridization with the states/bands deep below the conduction bands, presumably sulfur-dominated. Similarly, the hybridization of 3$d_{x_{1}^{2}-y_{1}^{2}}$ and 3$d_{z_{1}^{2}}$ orbitals with the states deep below the conduction band also provides some additional charge in these orbitals. We consider such details to be out of the scope of the present paper.


# Bringing PDEs to JAX with forward and reverse modes automatic
differentiation

Ivan Yashchuk
Department of Computer Science, Aalto University
VTT Technical Research Centre of Finland
ivan.yashchuk@{aalto,vtt}.fi

###### Abstract

Partial differential equations (PDEs) are used to describe a variety of
physical phenomena. Often these equations do not have analytical solutions and
numerical approximations are used instead. One of the common methods to solve
PDEs is the finite element method. Computing derivative information of the
solution with respect to the input parameters is important in many tasks in
scientific computing. We extend JAX automatic differentiation library with an
interface to Firedrake finite element library. High-level symbolic
representation of PDEs allows bypassing differentiating through low-level
possibly many iterations of the underlying nonlinear solvers. Differentiating
through Firedrake solvers is done using tangent-linear and adjoint equations.
This enables the efficient composition of finite element solvers with
arbitrary differentiable programs. The code is available at
github.com/IvanYashchuk/jax-firedrake.

##  1 Introduction

There is a growing interest in differentiable physical simulators and
including them as components in machine learning systems. Recent work in this
field includes either implementing physical simulators from scratch in popular
deep learning frameworks (Holl et al., 2020) or developing new software
specifically for physical simulations (Hu et al., 2019). In the case of
simulators consisting of the partial differential equation (PDE) solving
differentiating through all the operations of the PDE solver is not scalable
and heavy memory-bound. For PDEs, explicit expressions are available for
deriving efficient functions needed for automatic differentiation (AD). In
this work, we integrate Firedrake PDE library with JAX AD library. As a
result, it is possible to compose Firedrake models with arbitrary programs
that JAX can differentiate.

##  2 Background

The variational formulation also known as weak formulation allows finding the
solution to problems modeled with PDEs using an integral form. Having the
integral form gives simpler equations to solve using linear algebra methods
over a vector space of infinite dimension or functional space. In finite
element methods (FEM), the variational formulation of the problem is
transformed into a system of nonlinear equations for unknown finite element
function coefficients that can be solved numerically.

Let a system of PDEs that describe the physics of the problem of interest be
written as

$F(u,m)=0,$ (1)

where $u$ is the solution and $m$ represents parameters that affect the
solution. The solution $u$ can be written as an implicit function $u(m)$ of
the parameters $m$, then we can formulate its derivative $\frac{du}{dm}$. In
this work we consider the non-linear systems of PDEs that are discretized
using the finite element method.

Let $J(u,m)$ be a functional of interest, it represents any quantity that
depends on $u$ and $m$. Then the problem is computing the derivatives
$\frac{dJ(u(m),m)}{dm}$. It can be done using finite difference methods,
tangent-linear or adjoint approaches.

The tangent-linear system is the same idea as the forward mode of automatic
differentiation, while the adjoint approach corresponds to the reverse mode
automatic differentiation. In (Farrell et al., 2013), authors describe the
techniques that can be used to derive possibly time-dependant tangent-linear
and adjoint models. Having the solution to the tangent-linear system we can
evaluate the gradient of any functional. However, in many applications the
functional is fixed and the goal is to calculate the derivative with respect
to any parameter the chosen functional depends on. In this situation, the
better approach is to use the adjoint equations for derivative computation.
Higher-order derivatives can be derived by composing tangent-linear and
adjoint models (Maddison et al., 2019).

### Firedrake library

Firedrake is an automated system for the solution of partial differential
equations using the finite element method. It uses the UFL language to express
variational problems (Rathgeber et al., 2016; Alnaes et al., 2014). Then this
high-level problem representation is compiled into low-level C code for
assembling vectors and matrices for the nonlinear solver. Nonlinear and linear
solvers are based on PETSc library (Balay et al., 2019). UFL representation of
the residual equation makes it straightforward to generate both tangent-linear
and adjoint equations using built-in automatic differentiation. However, pure
UFL implementation of derivative calculation is restricted only to one
variational problem. Time-dependent problems are solved as a sequence of
variational problems and to differentiate through this sequence program
execution should be traced. Two libraries interface directly with Firedrake
for automated derivative computation: dolfin-adjoint (Mitusch et al., 2019)
and tlm_adjoint (Maddison et al., 2019). So it is not necessary to use dolfin-
adjoint and relying only on UFL for adjoint and tangent-linear derivation is
possible if the Firedrake model is a single variational problem or time-
stepping is implemented in JAX. However, we chose to depend on dolfin-adjoint
library here as it also makes possible to calculate shape derivatives for
domain optimization and boundary conditions derivatives.

##  3 Automatically differentiating PDE solvers in JAX

JAX is a numerical computing library that includes forward and reverse mode
automatic differentiation. JAX uses terminology and concepts from differential
geometry, namely pushforward map for the forward mode AD and pullback map for
the reverse mode AD. These maps can also be described in terms of the
Jacobian: The pushforward is Jacobian-vector product (JVP), and pullback is
Jacobian-transpose-vector product, or vector-Jacobian product (VJP).

Given a function $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$, the Jacobian matrix of
$f$ evaluated at an input point $x\in\mathbb{R}^{n}$, denoted $\partial f(x)$,
is often thought of as a matrix of partial derivatives of size $m\times n$.
Alternatively, $\partial f(x)$ represents a linear map, which maps the tangent
space of the domain of $f$ at the point $x$ to the tangent space of the
codomain of $f$ at the point $f(x)$:

$\partial f(x):\mathbb{R}^{n}\to\mathbb{R}^{m}.$

This map is called the pushforward map of $f$ at $x$. Given input point
$x\in\mathbb{R}^{n}$ and a tangent vector $v\in\mathbb{R}^{n}$, we get back an
output tangent vector in $\mathbb{R}^{m}$. This mapping, from $(x,v)$ pairs to
output tangent vectors, is called the Jacobian-vector product, and written as

$(x,v)\mapsto\partial f(x)v.$

The pullback map of $f$ at $x$ is

$\partial f(x)^{*}:\mathbb{R}^{m}\to\mathbb{R}^{n}.$

The Jacobian-transpose-vector product is then

$(x,v)\mapsto\partial f(x)^{*}v.$

Referring back to the tangent-linear model in the context of PDEs, there
Jacobian-vector product corresponds to

$(m,v)\mapsto\frac{du}{dm}v:=\dot{u}_{v}\quad\in\mathbb{R}^{m},$ (2)

where $\dot{u}_{v}$ is the solution to the following tangent-linear equation
with right hand side defined as the derivative of $F(u,m)$ with respect to $m$
in the direction $v\in\mathbb{R}^{n}$:

$\frac{\partial F(u,m)}{\partial u}\dot{u}_{v}=-\frac{\partial
F(u,m)}{\partial m}\cdot v.$ (3)

Note that compared to the standard tangent-linear equations the unknown in the
equation 3 is not the full solution Jacobian matrix but its multiplication
with a vector, making the problem easier to solve.

Jacobian-transpose-vector product function is defined as

$(m,v)\mapsto\frac{du}{dm}^{*}v:=-\frac{\partial F(u,m)}{\partial
m}^{*}\cdot\lambda_{v}\quad\in\mathbb{R}^{n},$

where $\lambda_{v}$ is the solution to the following adjoint equation with
right hand side given by the cotangent vector $v\in\mathbb{R}^{m}$:

$\frac{\partial F(u,m)}{\partial u}^{*}\lambda_{v}=v.$ (4)

Compared to the standard adjoint equations here we do not need any information
about the functional and it is AD’s system responsibility to compose VJPs for
reverse pass and supply cotangent vectors.

Internally, JAX functions are called Primitives and for each Primitive
Jacobian-vector products (JVP), vector-Jacobian products (VJP), batching and
just-in-time compilation rules are implemented. To make JAX work with external
functions we need to implement new Primitives together with functions that
define vector-Jacobian and Jacobian-vector products.

We have implemented a new Primitive function for the Firedrake solver that
transforms inputs and outputs to appropriate data types and registers
associated JVP and VJP functions that solve tangent-linear and adjoint
equations respectively.

##  4 Numerical examples

For demonstrating the implementation we consider the Poisson equation as the
forward model problem. Let $\Omega\subset\mathbb{R}^{n},n\in\\{1,2,3\\}$ be an
open, bounded domain and consider the following problem:

$\begin{split}\quad-\nabla\cdot(\kappa\nabla u)&=f\text{ in }\Omega,\\\
u&=0\text{ on }\partial\Omega,\end{split}$ (5)

where $u:\Omega\to\mathbb{R}$ is the unknown temperature,
$\kappa\in\mathbb{R}$ is the thermal conductivity, $f:\Omega\to\mathbb{R}$ is
the source term ($f(x)>0$ corresponds to heating and $f(x)<0$ corresponds to
cooling).

The variational form of the equation 5: Find $u\in H_{0}^{1}(\Omega)$ such
that

$(\kappa\nabla u,\nabla v)_{L^{2}(\Omega)}-(f,v)_{L^{2}(\Omega)}=0,\text{ for
all }v\in H_{0}^{1}(\Omega),$ (6)

where $H_{0}^{1}(\Omega)$ is the space of functions vanishing on
$\partial\Omega$ with square integrable derivatives.
$(\cdot\,,\cdot)_{L^{2}(\Omega)}$ denotes the $L^{2}$-inner product.

### Optimal control of the Poisson equation

We solve the standard problem in PDE-constrained optimization: the optimal
control of the Poisson equation. Physically, this problem can be interpreted
as finding the best heating or cooling of a surface to achieve a desired
temperature profile. The problem is to minimize the following functional

$\min_{f}J(f):=\frac{1}{2}\int_{\Omega}(u-u_{d})^{2}\,dx+\frac{\gamma}{2}\int_{\Omega}f^{2}\,ds,$
(7)

where $u$ is the solution to the Poisson equation 6, $u_{d}$ is the desired
temperature profile, $f$ is the unknown control function, $\gamma$ is the
regularization parameter. Additionally, $f$ is constrained to bounds $a,b$
such that $a\leq f\leq b$.

The unknown function $f$ is discretized in a finite element space such that
values of $f$ at each cell of the mesh are the optimization parameters. As the
mesh is refined the number of parameters in the optimization problem
increases.

For our example we take the desired temperature to be
$u_{d}=\frac{1}{2\pi^{2}}\frac{1}{1+4\gamma\pi^{4}}\sin(\pi x)\sin(\pi y)$,
$\gamma=10^{-6}$, and $f$ is bounded between $0$ and $0.8$. L-BFGS-B optimizer
from SciPy library (Virtanen et al., 2020) is used with the gradient values
calculated by JAX. Convergence is achieved after 38 iterations with gradient
norm tolerance $10^{-10}$.

### Coefficient field inversion with neural network representation

Representing the model parameters at each point in space quickly leads to a
large number of model parameters. The neural network can be used as an
approximation to the spatially varying coefficients characterized by the
weights of the neural network. The optimization problem is then posed in the
space of network-parameters rather than at each cell of the computational
grid. In the task of topology optimization, it was shown that neural network
representation of the solution to the optimization problem helps to find a
better design in many cases (Hoyer et al., 2019). In (Berg & Nyström, 2017),
the authors studied neural network representation stability for inverse
problems.

In this example, we demonstrate the use of neural nets to parameterize inputs
of Firedrake solver and demonstrate the differentiability of the pipeline. The
problem is to minimize the following functional

$\min_{\kappa}J(\kappa):=\frac{1}{2}\int_{\Omega}(u-u_{m})^{2}\,dx+\frac{\gamma}{2}\int_{\Omega}|\nabla\kappa|^{2}\,dx,$
(8)

where $u$ is the solution to the Poisson equation 6, $u_{m}$ is the noisy
temperature measurement, $\kappa$ is the unknown material coefficient field,
$\gamma$ is the regularization parameter.

Here, for the simplicity, we choose feed-forward neural network with single
hidden layer with 10 neurons and tanh activation function. We set up synthetic
measurements temperature data with spatially varying coefficient
$\kappa=1+x+y$ and adding small noise to the true temperature state.

| # parameters | # iterations | $\|\kappa_{opt}-\kappa_{true}\|$ | $\|u_{opt}-u_{true}\|$
---|---|---|---|---
FEM | 981 | 31 | 1.356112e-1 | 2.711260e-4
NN | 47 | 89 | 5.362389e-2 | 1.999892e-4
Table 1: Performance comparison between FEM and NN representation of the
coefficient.

L-BFGS-B optimizer is used with gradient norm tolerance $10^{-6}$. Number of
iterations and $L^{2}$-norm of the difference between the optimal solution and
the true one are summarized in Table 1. In this problem, finite element
representation overfits to the noise present in the measurements while neural
network representation can recover material coefficient function closer to the
true function and visually seems to be unaffected by the noise (Figure 1).

Noisy temperature field

Inverted $\kappa$, FEM

Inverted $\kappa$, NN

Figure 1: Conductivity field inversion from noisy measurements using FEM and
NN representations of optimization parameters.

##  5 Conclusions

We describe an extension to JAX that allows a seamless inclusion of PDE
solvers written in Firedrake into arbitrary JAX differentiable programs,
including, but not limited to, deep learning models, Bayesian probabilistic
models. Ongoing work targets further integration of Firedrake code with JAX
allowing just-in-time compilation and GPU computing and better
interoperability of JAX parallelism with Firedrake’s MPI parallelism.

## References

  * Alnaes et al. (2014) Martin S. Alnaes, Anders Logg, Kristian B. Oelgaard, Marie E. Rognes, and Garth N. Wells.  Unified form language: A domain-specific language for weak formulations of partial differential equations.  _ACM Trans. Math. Softw._ , 40(2), March 2014.  ISSN 0098-3500.  doi: 10.1145/2566630.  URL https://doi.org/10.1145/2566630.
  * Balay et al. (2019) Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Lisandro Dalcin, Alp Dener, Victor Eijkhout, William D. Gropp, Dmitry Karpeyev, Dinesh Kaushik, Matthew G. Knepley, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Todd Munson, Karl Rupp, Patrick Sanan, Barry F. Smith, Stefano Zampini, Hong Zhang, and Hong Zhang.  PETSc Web page.  https://www.mcs.anl.gov/petsc, 2019.  URL https://www.mcs.anl.gov/petsc.
  * Berg & Nyström (2017) Jens Berg and Kaj Nyström.  Neural network augmented inverse problems for pdes, 2017.
  * Farrell et al. (2013) Patrick E. Farrell, David A. Ham, Simon W. Funke, and Marie E. Rognes.  Automated derivation of the adjoint of high-level transient finite element programs.  _SIAM Journal on Scientific Computing_ , 35(4):C369–C393, January 2013.  doi: 10.1137/120873558.  URL https://doi.org/10.1137/120873558.
  * Holl et al. (2020) Philipp Holl, Vladlen Koltun, and Nils Thuerey.  Learning to control pdes with differentiable physics, 2020.
  * Hoyer et al. (2019) Stephan Hoyer, Jascha Sohl-Dickstein, and Sam Greydanus.  Neural reparameterization improves structural optimization, 2019.
  * Hu et al. (2019) Yuanming Hu, Luke Anderson, Tzu-Mao Li, Qi Sun, Nathan Carr, Jonathan Ragan-Kelley, and Frédo Durand.  Difftaichi: Differentiable programming for physical simulation, 2019.
  * Maddison et al. (2019) James R. Maddison, Daniel N. Goldberg, and Benjamin D. Goddard.  Automated calculation of higher order partial differential equation constrained derivative information.  _SIAM Journal on Scientific Computing_ , 41(5):C417–C445, January 2019.  doi: 10.1137/18m1209465.  URL https://doi.org/10.1137/18m1209465.
  * Mitusch et al. (2019) Sebastian Mitusch, Simon Funke, and Jørgen Dokken.  dolfin-adjoint 2018.1: automated adjoints for FEniCS and firedrake.  _Journal of Open Source Software_ , 4(38):1292, June 2019.  doi: 10.21105/joss.01292.  URL https://doi.org/10.21105/joss.01292.
  * Rathgeber et al. (2016) Florian Rathgeber, David A. Ham, Lawrence Mitchell, Michael Lange, Fabio Luporini, Andrew T. T. Mcrae, Gheorghe-Teodor Bercea, Graham R. Markall, and Paul H. J. Kelly.  Firedrake.  _ACM Transactions on Mathematical Software_ , 43(3):1–27, December 2016.  doi: 10.1145/2998441.  URL https://doi.org/10.1145/2998441.
  * Virtanen et al. (2020) Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, CJ Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake Vand erPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1. 0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 2020.  doi: https://doi.org/10.1038/s41592-019-0686-2.

##  Appendix A Derivations of tangent-linear and adjoint equations

Begin with the non-linear system of equations given by equation 1. Then,
taking the total derivative of the above equation with respect to the
parameter $m$ gives

$\frac{d}{dm}F(u,m)=\frac{\partial F(u,m)}{\partial
u}\frac{du}{dm}+\frac{\partial F(u,m)}{\partial m}\frac{dm}{dm}=0$ (9)

and re-arranging gives the tangent-linear equation associated with the PDE
equation 1

$\frac{\partial F(u,m)}{\partial u}\frac{du}{dm}=-\frac{\partial
F(u,m)}{\partial m}.$ (10)

The tangent equation is always linear.

Consider a functional $J(u,m)$. Let $\widehat{J}(m):=J(u(m),m)$ be a pure
function of $m$. Applying the chain rules gives the expression for the
gradient

$\frac{d}{dm}\widehat{J}(m)=\frac{\partial J(u,m)}{\partial
u}\frac{du}{dm}+\frac{\partial J(u,m)}{\partial m}.$ (11)

Having the solution to the tangent-linear equation 10 we can evaluate the
gradient of any functional $\widehat{J}$. However, in many applications the
functional is fixed and the goal is to calculate the derivative with respect
to any parameter the chosen functional depends on. In this case, the
alternative is the adjoint approach.

Suppose the tangent-linear system is invertible. Then rewrite the solution to
the equation 10 as

$\frac{du}{dm}=-\left(\frac{\partial F(u,m)}{\partial
u}\right)^{-1}\frac{\partial F(u,m)}{\partial m}$ (12)

and substitute $\frac{du}{dm}$ into the expression for the gradient of
$\widehat{J}$:

$\frac{d\widehat{J}(m)}{dm}=-\frac{\partial J(u,m)}{\partial
u}\left(\frac{\partial F(u,m)}{\partial u}\right)^{-1}\frac{\partial
F(u,m)}{\partial m}+\frac{\partial J(u,m)}{\partial m}.$ (13)

Take the adjoint of the above equation:

$\frac{d\widehat{J}(m)}{dm}^{*}=-\frac{\partial F(u,m)}{\partial
m}^{*}\left(\frac{\partial F(u,m)}{\partial u}\right)^{-*}\frac{\partial
J(u,m)}{\partial u}^{*}+\frac{\partial J(u,m)}{\partial m}^{*}.$ (14)

Now define the new variable as

$\lambda=\left(\frac{\partial F(u,m)}{\partial u}\right)^{-*}\frac{\partial
J(u,m)}{\partial u}^{*}.$ (15)

This new variable is called the adjoint variable and it is the solution of the
adjoint equation:

$\frac{\partial F(u,m)}{\partial u}^{*}\lambda=\frac{\partial J(u,m)}{\partial
u}^{*}.$ (16)

The adjoint equation is always linear.