# An Accelerated Lyapunov Function
for Polyak’s Heavy-Ball on Convex Quadratics

Antonio Orvieto
ETH Zürich

###### Abstract

In 1964, Polyak showed that the Heavy-ball method, the simplest momentum
technique, accelerates convergence of strongly-convex problems in the vicinity
of the solution. While Nesterov later developed a globally accelerated
version, Polyak’s original algorithm remains simpler and more widely used in
applications such as deep learning. Despite this popularity, the question of
whether Heavy-ball is also globally accelerated or not has not been fully
answered yet, and no convincing counterexample has been provided. This is
largely due to the difficulty in finding an effective Lyapunov function:
indeed, most proofs of Heavy-ball acceleration in the strongly-convex
quadratic setting rely on eigenvalue arguments. Our study adopts a different
approach: studying momentum through the lens of quadratic invariants of simple
harmonic oscillators. By utilizing the modified Hamiltonian of Stormer-Verlet
integrators, we are able to construct a Lyapunov function that demonstrates an
$O(1/k^{2})$ rate for Heavy-ball in the case of convex quadratic problems.
This is a promising first step towards potentially proving the acceleration of
Polyak’s momentum method and we hope it inspires further research in this
field.

##  1 Introduction

The problem of unconstrained continuous convex optimization consists in
finding an element of the set
$\operatorname*{arg\,min}_{q\in\mathbb{R}^{d}}f(q)$, for some lower bounded
convex $f:\mathbb{R}^{d}\to\mathbb{R}$, generally assumed to be regular, e.g.,
twice continuously differentiable: $f\in C^{2}(\mathbb{R}^{d},\mathbb{R})$.

###  1.1 Acceleration in discrete- and continuous-time

In 1979 Nemirovsky and Yudin [17] showed that, if $f$ is convex and
$L$-smooth111A differentiable function $f:\mathbb{R}^{d}\to\mathbb{R}$ is said
to be $L$-smooth if it has $L$-Lipschitz gradients., no gradient-based
optimizer can converge to a solution faster than $O(1/k^{2})$, where $k$ is
the number of gradient evaluations222This lower bound holds just for $k<d$
hence it is only interesting in the high-dimensional setting.. While Gradient
Descent (GD) converges like $O(1/k)$, the optimal rate $O(1/k^{2})$ is
achieved by the celebrated Accelerated Gradient Descent (AGD) method, proposed
by Nesterov in 1982 [19]: starting from $p_{0}=0$ and a random $q_{0}$, the
approximation $q_{k}$ to a problem solution $q^{*}$ is computed iteratively
as333Many similar writings are possible. Here, we consider the particular
version studied in [22] and a physicist notation, where $p_{k}$ is a velocity
variable. This makes the connection to continuous-time cleaner and consistent
with recent work on the geometry of momentum methods [5].

$\begin{cases}q_{k+1}=q_{k}+\beta_{k}hp_{k}-h^{2}\nabla
f(q_{k}+\beta_{k}hp_{k})\\\ p_{k+1}=(q_{k+1}-q_{k})/h\end{cases}.$ (AGD)

where $\beta_{k}=\frac{k-1}{k+2}$ and $h^{2}$ is the step-size (we use the
notation $h^{2}$ instead of the standard $\eta$ for a reason which will become
apparent in the next sections). Interestingly, the different behaviour of GD
and AGD is retained in the continuous-time limit (as the step-size vanishes),
recently studied by Su, Boyd and Candes [22], but already present in the
seminal works of Polyak [20] and Gavurin [6]:

$\dot{q}+\nabla f(q)=0;$ (GD-ODE) $\ddot{q}+\frac{r}{t}\dot{q}+\nabla f(q)=0.$
(AGD-ODE)

Namely, we have that GD-ODE converges like $O(1/t)$ and AGD-ODE (with $r\geq
3$) like $O(1/t^{2})$, where $t>0$ is the time variable. This result gave
researchers a new tool to grasp the baffling essence (see discussion in [2,
22]) of accelerated optimizers, and led to the design of many novel fast
interpretable algorithms [1, 14, 25, 24].

###  1.2 Evaluating gradients at a shifted position

There are two modifications of GD that bring AGD about:

  1. 1.

inclusion of the momentum term (i.e. using $\beta_{k}\neq 0$);

  2. 2.

change in gradient extrapolation point:
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\nabla
f(q_{k})}\to{\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@color@cmyk@fill{0}{1}{0}{0}\nabla
f(q_{k}+\beta_{k}hp_{k})}.$

Questions arise immediately:

Are both these modifications necessary for acceleration?
In particular, is evaluating the gradient at non-iterate points
crucial or even necessary for acceleration?

To put these questions in the right historical context, one has to go back to
Polyak’s 1964 seminal paper [20], where the very first momentum method was
proposed for $C^{2}$ and $\mu$-strongly-convex problems444In the strongly-
convex case, $\beta_{k}$ is not monotonically increasing, but is instead
chosen to be a constant dependent on the strong-convexity modulus $\mu$, that
is $\beta=\left(\frac{\sqrt{L}-\sqrt{\mu}}{\sqrt{L}+\sqrt{\mu}}\right)^{2}$..
Using an elegant functional-analytic argument on multistep methods, Polyak
proved that momentum alone — without shifted gradient evaluation (a.k.a.
Heavy-ball (HB), see equation below) — is able to achieve acceleration555Here
to be intended as a dependency of the rate on the square root of the condition
number $L/\mu$. in a neighborhood of the solution. This local argument becomes
of course global in the quadratic case (for a simplified proof, see
Proposition 1 in [15]).

$\begin{cases}q_{k+1}=q_{k}+\beta_{k}hp_{k}-h^{2}\nabla f(q_{k})\\\
p_{k+1}=(q_{k+1}-q_{k})/h\end{cases}.$ (HB)

Despite the many attempts, nobody in the last 56 years has been able to show
that HB has a global (i.e. for any initialization) accelerated rate — neither
in the strongly-convex case (using a fixed momentum) nor in the non-strongly-
convex case (using an increasing $\frac{k-1}{k+2}$ momentum). Beyond the
technical difficulty, another plausible reason may also be lack of interest,
as the introduction of Nesterov’s globally accelerated method in 1982, that
overshadowed the conceptually simpler method from Polyak.

However, many researchers in the last decade, supported by numerical evidence
and by the success of Heavy-ball in deep learning [13], expressed their belief
that HB is accelerated:

[…] supported by the numerical simulations we envisage that the convergence
factor could be strengthened even further. This is indeed left as a future
work.

– Ghadimi et al. [7], 2015

Despite the long history of this approach, there is still an open question
whether the heavy ball method converges to the optimum globally with
accelerated rate when the objective function is twice continuous
differentiable.

– Gorbunov et al. [8], 2019

Neither the evaluation of the gradient at a shifted position, nor a
specifically engineered damping parameter, as for example proposed in Nesterov
(2004, Sec. 2.2), seem666After talking to the first author, we decided to
replace “are” (as in the original preprint) with “seem”: indeed, the argument
in [16] is asymptotic and therefore somewhat equivalent to the one of Polyak
[20]. necessary.

– Muehlebach and Jordan [16], 2020

Other researcher believe HB is not accelerated:

If we can translate this argument to the discrete case we can understand why
AGD achieves acceleration globally for strongly-convex functions but the
Heavy-ball method does not.

– Shi et al. [21], 2018

While on the theoretical side the opinion is mixed, on the experimental side
no numerical simulation777In [15], the authors show that there exist a
strongly-convex smooth function such that Heavy-ball does not converge.
However, as also pointed out by Ghadimi et al. [7], such function is not
$C^{2}$, and that a big step-size is used — which violates the convergence
conditions of Thm. 4 in [7]. As such, this function does not constitute a
proper counterexample. has been able to show that HB is not accelerated. In
Figure 1, we provide two examples for the non-strongly-convex case (i.e. $\mu$
very small, such that an increasing momentum is preferable, leading $1/k^{2}$
convergence as opposed to $(1-\sqrt{\mu/L})^{k}$). In particular, we show that
HB is comparable to AGD through the lens of the pathological lower-bounding
quadratic example introduced by [17] and used to construct the $O(1/k^{2})$
bound in convex optimization — at least until the effect of non-trivial
strong-convexity becomes dominant (at around $f(q_{k})=10^{-6}$).

(a) HB and AGD on the worst-case (lower bound) quadratic objective from
Nesterov [18].

(b) HB and AGD on ill-conditioned linear regression. The optimal step-size
$1/L$ was used.

Figure 1: For both examples HB with momentum $\frac{k-1}{k+2}$ exhibits an
accelerated $1/k^{2}$ convergence rate, even though AGD with momentum
$\frac{k-1}{k+2}$ is faster in a neighborhood of the optimizer due to strong-
convexity. Instead, GD violates the Nesterov $O(1/k^{2})$ upper bound. We
recall that, while Nesterov’s upper bound holds for all $k>0$, the
$O(1/k^{2})$ lower bound (originally discovered by Nemirovski and Yudin [17])
only holds at $k=d/2$ (for more details, check the discussion in [18]).

###  1.3 Contributions

The purpose of the manuscript at hand is to study the effect of shifts in
gradient extrapolation points on acceleration in convex optimization (i.e. to
study the difference between Heavy-ball and Nesterov’s method). In particular,
the next pages are organized as follows:

  1. 1.

We start from a continuous-time argument: inspired by a recent idea from
Flammarion and Bach [4], in Section 2.1 we show how AGD-ODE with damping $2/t$
can be derived from the equation of a simple harmonic oscillator:
$\ddot{u}=-Au$. By using Lyapunov equations and a simple change of variables,
we retrieve the Lyapunov function proposed by Su, Boyd and Candes [22] to
prove a rate $O(1/t^{2})$ for AGD-ODE. This procedure is principled and leads
to many insights on Lyapunov function design.

  2. 2.

In Section 2.2, we apply the same methodology in discrete time, and show that
HB with momentum $\frac{k-1}{k+1}$ can be derived from the Störmer–Verlet
discretization of the simple harmonic oscillator. Solving again Lyapunov’s
equations, we are able to show an $O(1/k^{2})$ rate for a Heavy-ball argorithm
for convex quadratics. While this rate is already present in [4], our proof
technique is different as it relies on a Lyapunov function as opposed to an
eigenvalue analysis.

  3. 3.

In Section 3, by generalizing the discrete-time Lyapunov function found in
Section 2.2 we derive a modified Heavy-ball method

$q_{k+1}=q_{k}+\frac{k-1}{k+r-1}(q_{k}-q_{k-1})-h^{2}\frac{k+\frac{r-2}{2}}{k+r-1}\nabla
f(q_{k}),$

with a rate of convergence $O(1/k^{2})$ for any $k\geq 2$ and $r\geq 2$. Our
result not only generalizes the theory in [4], but also provides an
interesting connection between the continuous and the discrete — as the used
Lyapunov function converges, in the limit $h\to 0$, to the one used in [22]
for $r\geq 2$.

#### Recent related works.

Very recently, Wang et al. [23] proved that Heavy-ball is accelerated for a
class of functions satisfying the Polyak-Łojasiewicz condition. Instead, here
we provide a Lyapunov function for the non-strongly-convex setting, where the
Polyak-Łojasiewicz constant vanishes. We remark that, for strongly-convex
quadratic potentials, Heavy-ball is already known to achieve acceleration
[15]. However, the eigenvalue argument used in [15] cannot be leveraged in the
non-strongly-convex setting, where the minimum eigenvalue can be arbitrarily
low. As such, our work provides insights on how to construct effective
Lyapunov functions in the non-quadratic case, where Lyanonov arguments are
often the go-to option.

##  2 From quadratic invariants of oscillators to accelerated rates

Our procedure in this section is inspired by a beautiful idea presented by
Flammarion and Bach [4]: it is sometimes possible to translate a time-
dependent convergence rate problem into a time-independent stability problem.
Here we go one step further, and show how, with an additional step
(computation of quadratic invariants), it is possible to derive Lyapunov
functions and rates for the corresponding algorithms. We first illustrate the
idea in continuous-time and then proceed with the discrete-time analysis.

Our starting point is the following ODE:

$\ddot{q}+\frac{2}{t}\dot{q}+\nabla f(q)=0.$ (AGD-ODE2)

From the analysis in [22], we know that on a quadratic
$f(q)=f^{*}+\frac{1}{2}\langle(q-q^{*}),A(q-q^{*})\rangle$, with $A$ positive
semidefinite and $f^{*}\in\mathbb{R}$, the solution converges to
$q^{*}\in\operatorname*{arg\,min}_{x\in\mathbb{R}^{d}}f(q)$ at the rate
$O(1/t^{2})$. To prove this rate, the authors in [22] use the following
Lyapunov function:

$V(q,t)=2t^{2}(f(q)-f^{*})+\|t\dot{q}+(q-q^{*})\|^{2}.$ (1)

We show here a constructive way to derive $V$ (Section 2.1) and then (Section
2.2) we apply the same procedure to get a Lyapunov function for Heavy-ball
(i.e., the discretization). For simplicity, we consider here $f^{*}=0$ and
$q^{*}=0$.

###  2.1 Lyapunov functions from continuous-time invariants

Consider an harmonic oscillator on the potential $f(u)=\frac{1}{2}u^{\top}Au$,
i.e. $\ddot{u}=-Au$. From basic physics, we know that such a system is
marginally stable (bounded dynamics). By choosing $u=tq$ we get
$\dot{u}=q+t\dot{q}$ and $\ddot{u}=\dot{q}+\dot{q}+t\ddot{q}$. This implies

$\dot{q}+\dot{q}+t\ddot{q}=\ddot{u}=-Atq\ \quad\Rightarrow\quad\
\ddot{q}+\frac{2}{t}\dot{q}+Aq=0.$

That is, AGD-ODE can be reconstructed from a simple linearized pendulum. By
introducing the variable $v=\dot{u}$, we can write the pendulum in phase space
as a linear dynamical system

$\begin{pmatrix}\dot{u}\\\ \dot{v}\end{pmatrix}=\begin{pmatrix}0&I\\\
-A&0\end{pmatrix}\begin{pmatrix}u\\\ v\end{pmatrix}.$

Hence, the pendulum has the form $\dot{y}=Fy$, where $y=(u,v)$. We would now
like to get a Lyapunov function for this system. To do this, we recall a
fundamental proposition (check Thm. 4.6. in [12]).

######  Proposition 1 (Continuous-time Lyapunov equations).

The linear system $\dot{y}=Fy$ is Lyapunov stable if and only if for all
positive semidefinite matrices $Q$, there exists a symmetric matrix $P$ such
that

$PF+F^{T}P=-Q.$ (2)

Moreover, $V(y)=y^{T}Py$ is a Lyapunov function and $\dot{V}(y)=-y^{T}Qy$.

Since we know that a pendulum is only marginally stable (i.e., not
asymptotically stable), we can limit ourselves to the choice of a null matrix
$Q$. Hence, we need to solve the Lyapunov equation $PF=-F^{T}P$ for $P$. A
solution to this equation (many exist) is $P=\begin{pmatrix}A&0\\\
0&I\end{pmatrix}$, which implies that

$V(u)=\langle u,Au\rangle+\|v\|^{2}$ (3)

is a quadratic invariant, i.e. $\dot{V}(u)=0$. This is well known, since $V$
is actually twice the total energy (Hamiltonian) of the pendulum. Finally, we
can change variables and get that

$V(q)=t^{2}\langle
q,Aq\rangle+\left\|\frac{{\mathrm{d}}}{{\mathrm{d}}t}(tq)\right\|^{2}=2t^{2}f(q)+\|t\dot{q}+(q-q^{*})\|^{2},$
(4)

is a Lyapunov function for AGD-ODE2, with $\dot{V}(q)=0$. This is precisely
equation 1.

#### From quadratic to convex.

With a small modification (using a factor $r-1$ instead of $1$), it is
possible to get a Lyapunov function that works for AGD-ODE in the more general
convex case.

######  Proposition 2 (Theorem 3 from [22]).

For convex $L$-smooth objectives, AGD-ODE converges at a rate $O(1/t^{2})$.
This follows from the fact that

$V(q,t)=2t^{2}(f(q)-f^{*})+\|t\dot{q}+(r-1)(q-q^{*})\|^{2}.$ (5)

is a Lyapunov function, for $r\geq 3$.

###  2.2 Discrete-time invariants

We apply the construction from the last subsection to the discrete case.
Inspired by Flammarion and Bach [4], we consider at first a slightly modified
HB:

$q_{k+1}=q_{k}+\frac{k-1}{k+1}(q_{k}-q_{k-1})-h^{2}\frac{k}{k+1}\nabla
f(q_{k}).$ (HB$2$)

This algorithm is the discrete-time equivalent of
$\ddot{q}+\frac{2}{t}\dot{q}+\nabla f(q)=0$. As for the continuous-time case,
we start from $f(q)=\frac{1}{2}\langle q,Aq\rangle$. In this case, HB$2$ can
be written as

$(k+1)q_{k+1}=2kq_{k}+(k-1)q_{k-1}-h^{2}A(kq_{k}).$

That is, if we set $u_{k}=kq_{k}$, we get

$u_{k+1}-2u_{k}+u_{k-1}=-h^{2}Au_{k}.$ (6)

With surprise, we recognize that this is the Störmer–Verlet method [11] on
$\ddot{u}=-Au$, with step-size $h$ (that’s why had $h^{2}$ from the very
beginning).

It would be natural, as for the continuous-time case, to consider the total
energy as a quadratic invariant to derive a Lyapunov function. However, it
turns out that, interestingly, the Störmer–Verlet method does not precisely
conserve the total energy: there are small oscillations (see Section 3 in
[9])! Taking into account such small oscillations (Figure 2) is of fundamental
importance — since they lead to a crucial modification of the invariants we
have to use.

Figure 2: The Störmer–Verlet method on a one-dimensional quadratic potential
(i.e., a simplified pendulum) does not conserve the total energy. Details on
this phenomenon can be found in [10, 9].

######  Proposition 3 (Discrete-time Lyapunov equations).

The system $y_{k+1}=Fy_{k}$ is Lyapunov stable if an only if for all positive
semidefinite matrices $Q$, there exists a symmetric matrix $P$ such that

$F^{T}PF-P=-Q.$ (7)

Moreover, $V(y)=y^{T}Py$ is a Lyapunov function and
$V(y_{k+1})-V(y_{k})=-y_{k}^{T}Qy_{k}$ for all $k$.

We apply the theorem above (for $Q=0$) to the linear system

$\begin{pmatrix}u_{k+1}\\\ v_{k+1}\end{pmatrix}=\begin{pmatrix}I-h^{2}A&hI\\\
-hA&I\end{pmatrix}\begin{pmatrix}u_{k}\\\ v_{k}\end{pmatrix}.$

Under the choice $v_{k}=(u_{k}-u_{k-1})/h$, this system is equivalent to
equation 6, i.e., the discretized pendulum we want to find a quadratic
invariant for. Solving the discrete Lyapunov equation gives us

$P=\begin{pmatrix}A&-hA/2\\\ -hA/2&I\end{pmatrix},$

and the associated modified total energy:

$V(u_{k},v_{k})=\underbrace{\langle
u_{k},Au_{k}\rangle+\|v_{k}\|^{2}}_{\text{continuous-time invariant
(energy)}}-\underbrace{h\langle v_{k},Au_{k}\rangle}_{\text{vanishing cross-
term}}.$ (8)

We make the following comments:

  * •

As $h\to 0$, the modified energy approaches the total energy in equation 3.
The purpose of the additional cross-term is to eliminate the small energy
oscillations we see in Figure 2.

  * •

Assuming without loss of generality that $A$ does not have zero eigenvalues,
$P$ is positive semidefinite (i.e., yields a valid Lyapunov function) if and
only if the Schur complement of $I$ (i.e., the block $P_{22}$) in $P$ is
positive semidefinite. That is, we need

$B:=A-\frac{h^{2}}{4}A^{2}=A\left(I-\frac{h^{2}}{4}A\right)\geq 0.$ (9)

Since $A$ and $(I-h^{2}A/4)$ are co-diagonalizable, the product is positive
semidefinite if and only if both $A\geq 0$ and $(I-h^{2}A/4)\geq 0$. This
requires $0\leq A\leq\frac{4}{h^{2}}I$, which in turns implies an upper bound
on the step-size $h^{2}$:

$h^{2}\leq\frac{4}{\lambda_{\text{max}}(A)}=\frac{4}{L}.$

The same condition (note that our step-size is $h^{2}$, not $h$) can be
deduced from the analysis of HB$2$ in [4].

Now it’s time to change variables back: $u_{k}=kq_{k}$. If we set
$p_{k}:=(q_{k}-q_{k-1})/h$, as also done in the introduction, we get

$\displaystyle hv_{k}$ $\displaystyle=h(u_{k}-u_{k-1})/h$
$\displaystyle=kq_{k}-(k-1)q_{k-1}+(k-1)q_{k}-(k-1)q_{k}$
$\displaystyle=(k-1)(q_{k}-q_{k-1})+q_{k}$ $\displaystyle=h(k-1)p_{k}+q_{k}.$

By substituting these formulas in equation 8, we get the following final form
for an effective Lyapunov function for HB$2$ — for the quadratic case:

$\displaystyle V_{k}=\langle u_{k},Au_{k}\rangle+\|v_{k}\|^{2}-\langle
hv_{k},Au_{k}\rangle$ $\displaystyle\Longrightarrow\ $ $\displaystyle
V_{k}=k^{2}\langle
q_{k},Aq_{k}\rangle+\frac{1}{h^{2}}\|h(k-1)p_{k}+q_{k}\|^{2}-k\langle
h(k-1)p_{k}+q_{k},Aq_{k}\rangle.$

To better understand this Lyapunov function, we multiply everything by $h^{2}$
and get

$V_{k}=(kh)^{2}\langle
q_{k},Aq_{k}\rangle+\|h(k-1)p_{k}+q_{k}\|^{2}-h^{2}k\langle
h(k-1)p_{k}+q_{k},Aq_{k}\rangle.$

Recalling that the “time” variable $t$ is defined to be $t_{k}=hk$, this cost
becomes

$V_{k}=t_{k}^{2}\langle
q_{k},Aq_{k}\rangle+\|t_{k-1}p_{k}+q_{k}\|^{2}-ht_{k}\langle
t_{k-1}p_{k}+q_{k},Aq_{k}\rangle.$

This Lyapunov function can be easily generalized by noting that $\langle
q_{k},Aq_{k}\rangle=2(f(q_{k})-f^{*})$ and $Aq_{k}=\nabla f(q_{k})$:

$V_{k}=2t_{k}^{2}(f(q_{k})-f^{*})+\|t_{k-1}p_{k}+q_{k}\|^{2}-t_{k}\langle\nabla
f(q_{k}),t_{k-1}p_{k}+q_{k}\rangle.$ (10)

Finally, note that

  * •

From equation 10 as $h\to 0$ we get equation 4: the continuous-time Lyapunov
function of [22].

  * •

The mixing term is necessary and makes the positive definiteness (see equation
9) of $V_{k}$ non-trivial.

All in all, in this subsection, we proved the following result.

######  Proposition 4.

Let $A\in\mathbb{R}^{d\times d}$ be positive semidefinite,
$f^{*}\in\mathbb{R}$ and $f(q)=f^{*}+\frac{1}{2}(q-q^{*})^{T}A(q-q^{*})$. Let
$(q_{k})_{k\geq 0}$ be the iterates of HB$2$ and $p_{k}:=(q_{k}-q_{k-1})/h$.
If the step-size $h^{2}<\frac{4}{\lambda_{\text{max}}(A)}$, then equation 10
is non-negative and such that $V_{k+1}-V_{k}=0$ along the HB$2$ trajectory,
for all $k$. From this, one can deduce an accelerated rate of $O(1/k^{2})$ in
suboptimality.

Details of the proof are given in the proof of Theorem 1 (more general).

##  3 Accelerated Heavy-ball methods for convex quadratics

In this section, we start to lift the discussion to the convex non-quadratic
setting, by providing a generalization of HB2\. Indeed, we know from the
continuous-time analysis in [22] that $\ddot{q}+\frac{2}{t}\dot{q}+\nabla
f(q)=0$ may not have an accelerated rate for functions which are convex but
not necessarily quadratic. In this case, a rate of $O(1/t^{2})$ only holds
888The case $0<r\leq 3$ was studied by Attouch et al. [3]: a convergence rate
of $O(1/t^{p})$ with $p<2r/3$ is shown in this case. The same result also
holds in discrete time. for

$\ddot{q}+\frac{r}{t}\dot{q}+\nabla f(q)=0,$

with $r\geq 3$. In the same way, we expect that HB$2$ (which is the
discretization for $r=2$) may not have an accelerated rate in the convex non-
quadratic setting and a generalization corresponding to high friction is
therefore necessary.

Our objective in this chapter is to construct such a generalization of HB$2$,
which we name HB$r$.

###  3.1 A generalized Heavy-ball with high friction and guarantees on
quadratics

After a few weeks of intense calculations, we found that this algorithm gives
the desired result (Thm. 1).

$q_{k+1}=\underbrace{q_{k}+\frac{k-1}{k+r-1}(q_{k}-q_{k-1})}_{\text{iterate +
momentum}}-h^{2}\underbrace{\frac{k+\frac{r-2}{2}}{k+r-1}\nabla
f(q_{k})}_{\text{scaled gradient of iterate}}.$ (HB$r$)

First, note that $r=2$ recovers HB$2$ — which we proved to be accelerated in
the last subsection using a novel Lyapunov argument. The second, and perhaps
the most crucial, thing to note is that HB$r$ recalls the high friction
generalization of AGD proposed by [22] (see Theorem 6 in their paper):

$q_{k+1}=\underbrace{q_{k}+\frac{k-1}{k+r-1}(q_{k}-q_{k-1})}_{\text{iterate +
momentum}}-h^{2}\underbrace{\nabla
f\left(q_{k}+\frac{k-1}{k+r-1}(q_{k}-q_{k-1})\right)}_{\text{gradient of
[iterate + momentum]}}.$ (AGD$r$)

Between HB$r$ and AGD$r$ there are a few important differences:

  * •

In AGD$r$ the gradient is evaluated at
$q_{k}+\frac{k-1}{k+r-1}(q_{k}-q_{k-1})$, while in HB$r$ it is evaluated at
$q_{k}$.

  * •

in HB$r$ the effective step-size (i.e. what multiplies the gradient) is
iteration-dependent, and goes from $h^{2}/2$ to $h^{2}$ as $k\to\infty$. We
believe this has not to be regarded as part of the acceleration mechanism: it
is just a small modification needed to make the analysis easier.

  * •

Arguably HB$r$ (neglecting the small correction) is conceptually simpler that
AGD$r$: compared to GD, only a momentum term is added at each iteration — and
this can be thought of as the source of acceleration.

We proceed in proving that HB$r$ is accelerated in the quadratic case.

######  Theorem 1.

Let $A\in\mathbb{R}^{d\times d}$ be positive semidefinite,
$f^{*}\in\mathbb{R}$ and $f(q)=f^{*}+\frac{1}{2}(q-q^{*})^{T}A(q-q^{*})$. Let
$(q_{k})_{k\geq 0}$ be the iterates of HB$r$ ($r\geq 2$) and
$p_{k}:=(q_{k}-q_{k-1})/h$. If $h^{2}\leq\frac{4}{\lambda_{\text{max}}(A)}$,
then
$V_{k}=2(k+r-2)^{2}h^{2}(f(q_{k})-f^{*})+\|h(k-1)p_{k}+(r-1)(q_{k}-q^{*})\|^{2}\\\
-h^{2}(k+r-2)\langle\nabla f(q_{k}),h(k-1)p_{k}+(r-1)(q_{k}-q^{*})\rangle$
(11) is non-negative and such that $V_{k+1}-V_{k}\leq 0$ along the HB$r$
trajectory, for all $k$. Moreover, for any
$h^{2}<\frac{4}{\lambda_{\text{max}}(A)}$, HB$r$ is accelerated. In
particular, if $h^{2}=\frac{2}{\lambda_{\text{max}}(A)}$, we have the rate
$f(q_{k})-f^{*}\leq\frac{\lambda_{\text{max}}(A)V_{0}}{2(k+r-2)^{2}}.$

We note a couple of facts about the Lyapunov function $V_{k}$ in equation 11.

  * •

It reduces to equation 10 in the case $r=2$. For $r>2$, it is a graspable
generalization of equation 10 — which we instead derived in a systematic way
using Lyapunov equations. The term $(r-1)$ is inspired by the continuous-time
limit in equation 5.

  * •

Consider the Lyapunov function above, but without cross term i.e.

$2(k+r-2)^{2}h^{2}(f(q_{k})-f^{*})+\|h(k-1)p_{k}+(r-1)(q_{k}-q^{*})\|^{2}.$

This function works for proving an $O(1/k^{2})$ rate for AGDr (it’s a Lyapunov
function, see Thm. 6 from [22]). Therefore, higher complexity (i.e., an
additional cross term) is needed to study the acceleration of Heavy-ball, when
compared to Nesterov’s method.

  * •

As $h\to 0$, the cross term vanishes $V_{k}$ converges to equation 5 — its
continuous-time equivalent. Indeed, both HBr and AGDr converge to AGD-ODE as
$h\to 0$.

###  3.2 Proof of the theorem

It is useful to simplify equation 11 and to work with variables $q_{k}$ and
$q_{k-1}$ — a natural choice in the discrete setting. We split the Lyapunov
function into two parts: $V_{k}=V^{1}_{k}+V^{2}_{k}$.

$\displaystyle V^{1}_{k}:=2(k+r-2)^{2}h^{2}(f(q_{k})-f^{*})$ (12)
$\displaystyle\quad\quad\quad-h^{2}(k+r-2)\langle\nabla
f(q_{k}),(k-1)(q_{k}-q_{k-1})+(r-1)(q_{k}-q^{*})\rangle.$ $\displaystyle
V^{2}_{k}:=\|(k-1)(q_{k}-q_{k-1})+(r-1)(q_{k}-q^{*})\|^{2}.$ (13)

First, we are going to study $V_{k}^{2}$ in the non-quadratic case, and then
$V_{k}^{1}$ in the quadratic case. Theorem 1 will follow from a combination of
the two corresponding lemmata.

The first lemma shares many similarities with the proof of Theorem 1 in [7].

######  Lemma 1.

For any differentiable function $f:\mathbb{R}^{d}\to\mathbb{R}$ (not
necessarily convex or $L$-smooth) and any sequence of iterates $(q_{k})_{k\geq
0}$ returned by HB$r$, we have:

$\displaystyle V_{k+1}^{2}-V_{k}^{2}\ =\ $
$\displaystyle-h^{2}(r-1)(2k+r-2)\langle\nabla f(q_{k}),q_{k}-q^{*}\rangle$
$\displaystyle-h^{2}(k-1)(2k+r-2)\langle\nabla f(q_{k}),q_{k}-q_{k-1}\rangle$
$\displaystyle+\frac{h^{4}}{4}(2k+r-2)^{2}\|\nabla f(q_{k})\|^{2},$

where $V^{2}_{k}$ is defined in equation 13.

###### Proof.

Let $g_{k}:=(k-1)(q_{k}-q_{k-1})+(r-1)(q_{k}-q^{*})$, then,

$V^{2}_{k+1}-V^{2}_{k}=\|g_{k+1}\|^{2}-\|g_{k}\|^{2}=\langle
g_{k+1}+g_{k},g_{k+1}-g_{k}\rangle.$

We proceed in computing $g_{k+1}-g_{k}$. The algorithm symmetric structure
here is fundamental:

$\displaystyle g_{k+1}-g_{k}=\ $ $\displaystyle
k(q_{k+1}-q_{k})+(r-1)(q_{k+1}-q^{*})-(k-1)(q_{k}-q_{k-1})-(r-1)(q_{k}-q^{*})$
$\displaystyle=\ $
$\displaystyle(k+r-1)q_{k+1}-(k+r-1)q_{k}-(k-1)(q_{k}-q_{k-1})$
$\displaystyle\overset{\text{({HB}$r$)}}{=}$
$\displaystyle-h^{2}\left(k+\frac{r-2}{2}\right)\nabla f(q_{k}).$

Instead, $g_{k+1}+g_{k}$ is slightly more complex.

$\displaystyle g_{k+1}+g_{k}=\ $ $\displaystyle
k(q_{k+1}-q_{k})+(r-1)(q_{k+1}-q^{*})+(k-1)(q_{k}-q_{k-1})+(r-1)(q_{k}-q^{*})$
$\displaystyle=\ $
$\displaystyle(k+r-1)q_{k+1}+(-k+r-1)q_{k}+(k-1)(q_{k}-q_{k-1})-2(r-1)q^{*}$
$\displaystyle\overset{\text{({HB}$r$)}}{=}$
$\displaystyle(k+r-1)q_{k}+(k-1)(q_{k}-q_{k-1})-h^{2}\left(k+\frac{r-2}{2}\right)\nabla
f(q_{k})$ $\displaystyle+(-k+r-1)q_{k}+(k-1)(q_{k}-q_{k-1})-2(r-1)q^{*}$
$\displaystyle=\ $ $\displaystyle
2(r-1)(q_{k}-q^{*})+2(k-1)(q_{k}-q_{k-1})-h^{2}\left(k+\frac{r-2}{2}\right)\nabla
f(q_{k}).$

The proof is concluded by taking the inner product. ∎

We proceed by computing the difference $V^{1}_{k+1}-V^{1}_{k}$. Our
calculations will be very quick, since we can leverage, in the quadratic case,
on a simplified expression for $V^{1}_{k}$.

######  Lemma 2.

Let $V^{1}_{k}$ be defined as equation 12. In the context of Theorem 1, we
have

$V^{1}_{k}=h^{2}(k+r-2)(k-1)\langle q_{k-1}-q^{*},A(q_{k}-q^{*})\rangle$

and

$\displaystyle V_{k+1}^{1}-V_{k}^{1}\ =$ $\displaystyle\ h^{2}(2k+r-2)\langle
q_{k}-q^{*},A(q_{k}-q^{*})\rangle$ $\displaystyle+h^{2}(k-1)(2k+r-2)\langle
q_{k}-q^{*},A(q_{k}-q_{k-1})\rangle$
$\displaystyle-\frac{h^{4}}{2}(2k+r-2)k\|A(q_{k}-q^{*})\|^{2}.$

###### Proof.

From equation 12, we get

$\displaystyle V^{1}_{k}=\ $ $\displaystyle\ (k+r-2)^{2}h^{2}\langle
q_{k}-q^{*},A(q_{k}-q^{*})\rangle$ $\displaystyle-h^{2}(k+r-2)\langle
A(q_{k}-q^{*}),(k-1)(q_{k}-q_{k-1})+(r-1)(q_{k}-q^{*})\rangle$
$\displaystyle=\ $ $\displaystyle\ (k+r-2)^{2}h^{2}\langle
q_{k}-q^{*},A(q_{k}-q^{*})\rangle$ $\displaystyle\ -h^{2}(k+r-2)\langle
A(q_{k}-q^{*}),(k+r-2)(q_{k}-q^{*})-(k-1)(q_{k-1}-q^{*})\rangle$
$\displaystyle=\ $ $\displaystyle\ h^{2}(k+r-2)(k-1)\langle
q_{k-1}-q^{*},A(q_{k}-q^{*})\rangle.$

We proceed computing $V^{1}_{k+1}-V^{1}_{k}$ using this simplified form:

$\displaystyle V^{1}_{k+1}-V^{1}_{k}=\ $ $\displaystyle\ h^{2}(k+r-1)k\langle
q_{k}-q^{*},A(q_{k+1}-q^{*})\rangle$ $\displaystyle-h^{2}(k+r-2)(k-1)\langle
q_{k-1}-q^{*},A(q_{k}-q^{*})\rangle$ $\displaystyle=\ $ $\displaystyle
h^{2}\langle q_{k}-q^{*},A\Delta_{k}\rangle,$

where

$\displaystyle\Delta_{k}:=\ $ $\displaystyle\
(k+r-1)k(q_{k+1}-q^{*})-(k+r-2)(k-1)(q_{k-1}-q^{*}).$

Now, recall the definition of HB$r$:

$(q_{k+1}-q^{*})=(q_{k}-q^{*})+\frac{k-1}{k+r-1}(q_{k}-q_{k-1})-h^{2}\frac{k+\frac{r-2}{2}}{k+r-1}A(q_{k}-q*),$

where we subtracted $q^{*}$ from both sides. By plugging this into
$\Delta_{k}$, we get

$\displaystyle\Delta_{k}=\ $ $\displaystyle\
(k+r-1)k(q_{k}-q^{*})+(k-1)k(q_{k}-q_{k-1})-\frac{h^{2}}{2}(2k+r-2)kA(q_{k}-q*)$
$\displaystyle-(k+r-2)(k-1)(q_{k-1}-q^{*})$ $\displaystyle=\ $ $\displaystyle\
(2k+r-2)k(q_{k}-q^{*})-(2k+r-2)(k-1)(q_{k-1}-q^{*})$
$\displaystyle-\frac{h^{2}}{2}(2k+r-2)kA(q_{k}-q^{*})$ $\displaystyle=\ $
$\displaystyle\
(2k+r-2)(q_{k}-q^{*})+(2k+r-2)(k-1)(q_{k}-q_{k-1})-\frac{h^{2}}{2}(2k+r-2)kA(q_{k}-q^{*}).$

The result follows after taking the inner product $h^{2}\langle
q_{k}-q^{*},A\Delta_{k}\rangle$. ∎

We are finally ready to prove the result.

###### Proof of Theorem 1.

First, we compute
$V_{k+1}^{1}-V_{k}=(V_{k+1}^{1}-V_{k}^{1})+(V_{k+1}^{2}-V_{k}^{2})$ using
Lemma 2 and Lemma 1 (written for quadratic $f$). Next, we show that a certain
condition on the step-size implies positivity of $V_{k}$ and a convergence
rate.

$\displaystyle(V_{k+1}^{1}-V_{k}^{1})+(V_{k+1}^{2}-V_{k}^{2})\ =$
$\displaystyle\ h^{2}(2k+r-2)\langle q_{k}-q^{*},A(q_{k}-q^{*})\rangle$
+h^2(k-1)(2k+r-2)⟨q_k-q^*, A(q_k-q_k-1)⟩
$\displaystyle-\frac{h^{4}}{2}(2k+r-2)k\|A(q_{k}-q^{*})\|^{2}$
$\displaystyle-h^{2}(r-1)(2k+r-2)\langle q_{k}-q^{*},A(q_{k}-q^{*})\rangle$
-h^2(k-1)(2k+r-2)⟨q_k-q^*,A(q_k-q_k-1)⟩
$\displaystyle+\frac{h^{4}}{4}(2k+r-2)^{2}\|A(q_{k}-q^{*})\|^{2}.$

Crucially, note that the terms including $\langle
q_{k}-q^{*},A(q_{k}-q_{k-1})\rangle$ cancel. This is necessary to make our
proof (or, probably, any proof) work, since such inner product between the
gradient and the momentum changes sign (infinitely) many times along the
trajectory, and therefore cannot be easily compared to other quantities. For
the same reason, in the corresponding continuous-time proof from [22], the
terms including $\langle\nabla f(q),p\rangle$ also perfectly cancel out.

All in all, by collecting some terms, we get

$\displaystyle V_{k+1}-V_{k}$ $\displaystyle=-h^{2}(r-2)(2k+r-2)\langle
q_{k}-q^{*},A(q_{k}-q^{*})\rangle$
$\displaystyle\quad\quad+\frac{h^{4}}{4}(2k+r-2)(r-2)\|A(q_{k}-q^{*})\|^{2}$
$\displaystyle=-h^{2}(r-2)(2k+r-2)\langle q_{k}-q^{*},B(q_{k}-q^{*})\rangle,$

where $B:=A-\frac{h^{2}}{4}A^{2}$ is the matrix that we already studied in the
context of Lyapunov equations (see equation 9). Since $r\geq 2$, a sufficient
condition for $V_{k+1}-V_{k}\leq 0$ is $B\geq 0$, which holds under
$h^{2}\leq\frac{4}{\lambda_{\text{max}}(A)}$. As a sanity check, the reader
can appreciate the fact that, if $r=2$, then $V_{k+1}=V_{k}$ — as we already
proved in Proposition 4 (follows from the fact that $V_{k}$ solves the
Lyapunov equations).

Last, we have to translate the fact that $V_{k}$ is non-increasing to a
convergence rate. This is not trivial in our case, since $V_{k}$ also contains
a cross term which is not necessarily positive. Actually, we do not even know
that $V_{k}\geq 0$ yet! Hence, we have to come up with some tricks. We start
from rewriting the (simplified) Lyapunov function:

$V_{k}=h^{2}(k+r-2)(k-1)\langle
q_{k-1}-q^{*},A(q_{k}-q^{*})\rangle+\|(k-1)(q_{k}-q_{k-1})+(r-1)(q_{k}-q^{*})\|^{2}.$

Now, let us add and subtract a term $ch^{2}(k+r-2)^{2}\langle
q_{k}-q^{*},A(q_{k}-q^{*})\rangle$, with $c>0$. We have:

$V_{k}=ch^{2}(k+r-2)^{2}\langle
q_{k}-q^{*},A(q_{k}-q^{*})\rangle+\tilde{V}_{k},$

with

$\displaystyle\tilde{V}_{k}:=$ $\displaystyle-ch^{2}(k+r-2)^{2}\langle
q_{k}-q^{*},A(q_{k}-q^{*})\rangle$
$\displaystyle\quad+h^{2}(k+r-2)(k-1)\langle
q_{k-1}-q^{*},A(q_{k}-q^{*})\rangle$
$\displaystyle\quad+\|(k-1)(q_{k}-q_{k-1})+(r-1)(q_{k}-q^{*})\|^{2}.$

Now, if we show that $\tilde{V}_{k}$ is always positive, then $V_{k+1}\leq
V_{k}$ for all $k$ implies:

$ch^{2}(k+r-2)^{2}\langle q_{k}-q^{*},A(q_{k}-q^{*})\rangle\leq
ch^{2}(k+r-2)^{2}\langle
q_{k}-q^{*},A(q_{k}-q^{*})\rangle+\tilde{V}_{k}=V_{k}\leq V_{0},$

which gives the desired rate:

$f(q_{k})-f^{*}=\frac{1}{2}\langle
q_{k}-q^{*},A(q_{k}-q^{*})\rangle\leq\frac{V_{0}}{2ch^{2}(k+r-2)^{2}}.$

Therefore, we only need to show $\tilde{V}_{k}\geq 0$. To do this, we
introduce two new variables:

$u_{k}:=(k+r-2)(q_{k}-q^{*}),\quad\quad w_{k}:=(k-1)(q_{k-1}-q^{*}),$

and get a simplified form for $\tilde{V}_{k}$

$\displaystyle\tilde{V}_{k}$ $\displaystyle=-ch^{2}\langle
u_{k},Au_{k}\rangle+h^{2}\langle v_{k},Aw_{k}\rangle+\|u_{k}-w_{k}\|^{2}$
$\displaystyle=\langle u_{k},(I-ch^{2}A)u_{k}\rangle+\|w_{k}\|^{2}-2\langle
u_{k},\left(I-\frac{h^{2}}{2}A\right)w_{k}\rangle.$

Hence, we just need to show that

$\tilde{P}=\begin{pmatrix}I-ch^{2}A&-\left(I-\frac{h^{2}}{2}A\right)\\\
-\left(I-\frac{h^{2}}{2}A\right)&I\end{pmatrix}$

is positive definite, for some $c$ and $h^{2}$. Using the Schur
characterization for positive semidefinite matrices, $\tilde{P}\geq 0$ if and
only if

$\displaystyle 0\leq\tilde{B}(c)$
$\displaystyle:=I-ch^{2}A-\left(I-\frac{h^{2}}{2}A\right)^{2}$
$\displaystyle=I-ch^{2}A-I-\frac{h^{4}}{4}A^{2}+h^{2}A$
$\displaystyle=h^{2}A\left(1-c-\frac{h^{2}}{4}A\right).$

It is clear that $\tilde{B}(c)$ is positive semidefinite if and only if
$1-c-\frac{h^{2}}{4}\lambda_{\text{max}}(A)\geq 0$. That is,

$h^{2}\leq\frac{4(1-c)}{\lambda_{\text{max}}(A)}.$

Hence, for any $c\in(0,1)$ we get an acceleration. In particular, in the
theorem, we chose $c=\frac{1}{2}$. ∎

###  3.3 Numerical verification of our Lyapunov function

We verify numerically that the Lyapunov function for HBr proposed in equation
11 works on quadratics. To more clearly show the effect the inner product
correction term, which originated from the quadratic invariant of the
Störmer–Verlet method, we use here a slightly different notation:
$V_{k}=V^{11}_{k}+V^{12}_{k}+V^{2}_{k}$, with
$V^{1}_{k}=V^{11}_{k}+V^{12}_{k}$.

$\displaystyle V^{11}_{k}$ $\displaystyle:=2(k+r-2)^{2}h^{2}(f(q_{k})-f^{*});$
$\displaystyle V^{12}_{k}$ $\displaystyle:=-h^{2}(k+r-2)\langle\nabla
f(q_{k}),(k-1)(q_{k}-q_{k-1})+(r-1)(q_{k}-q^{*})\rangle;$ $\displaystyle
V^{2}_{k}$ $\displaystyle:=\|(k-1)(q_{k}-q_{k-1})+(r-1)(q_{k}-q^{*})\|^{2}.$

We recall that, the term $V^{12}_{k}$ (a.k.a. the cross-term) vanishes as
$h\to 0$, and is indeed not present in the continuous-time limit. We show that
this term, which we derived using Lyapunov equations in Sec. 2, plays a
fundamental role in ensuring $V_{k+1}-V_{k}\leq 0$. In Figure 3 we verify
numerically Thm. 1. In Figure 4 we show the essential role of $V^{12}$. Here
we used $h^{2}=1/L$, but HB$r$ can take larger steps (up to $4/L$), while the
other algorithms become unstable (Figure 5).

Figure 3: Dynamics of the Lyapunov function for HB$r$ on linear regression
(ill-conditioned Hessian, with condition number $\kappa$). Shown is the
behavior for $r=2,3$ with step-size $1/L$. For $r=2$, $V_{k}$ is constant, as
predicted by Prop. 4. For $r=3$, $V_{k}$ is decreasing as predicted by Thm. 1.
Figure 4: Same setting of the second example in Figure 3, but different
candidate Lyapunov function (no cross term). This confirms the cross-term is
necessary. Figure 5: HBr also works for big step-sizes (see conditions in
Theorem 1). Here, used is $h^{2}=3.9/L$.

##  4 Conclusion

In conclusion, the question of whether the Heavy-ball method is globally
accelerated for non-strongly-convex quadratic problems has yet to be fully
answered, and has attracted the attention of recent research [23]. Our study
takes a novel approach by examining momentum through the lens of quadratic
invariants of simple harmonic oscillators, and by utilizing the modified
Hamiltonian of Stormer-Verlet integrators we were able to construct a Lyapunov
function that demonstrates an $O(1/k^{2})$ rate for Heavy-ball in the case of
convex quadratic problems, where eigenvalues can vanish. This is a promising
first step towards potentially proving the acceleration of Polyak’s momentum
method through Lyapunov function arguments.

##  5 Acknowledgements

I would like to extend my deepest gratitude to Prof. Boris Polyak, Prof.
Christian Lubich, and Konstantin Mishchenko for the stimulating discussions.
My appreciation goes to Prof. Aurelien Lucchi and Prof. Thomas Hofmann for
their unwavering support and motivation, which helped me to develop the
project idea in Spring 2020. Lastly, I cannot express enough my gratitude to
Johannes Brahms for his Violinkonzert D-Dur op. 77, which provided the perfect
soundtrack to my late-night calculations, igniting my passion and drive to
push through the toughest moments.

## References

  * [1] Foivos Alimisis, Antonio Orvieto, Gary Bécigneul and Aurelien Lucchi  “A continuous-time perspective for modeling acceleration in Riemannian optimization”  In _International Conference on Artificial Intelligence and Statistics_ , 2020
  * [2] Zeyuan Allen-Zhu and Lorenzo Orecchia  “Linear coupling: An ultimate unification of gradient and mirror descent”  In _arXiv:1407.1537_ , 2014
  * [3] Hedy Attouch, Zaki Chbani and Hassan Riahi  “Rate of convergence of the Nesterov accelerated gradient method in the subcritical case $\alpha\leq$ 3”  In _ESAIM: Control, Optimisation and Calculus of Variations_ EDP Sciences, 2019
  * [4] Nicolas Flammarion and Francis Bach  “From averaging to acceleration, there is only a step-size”  In _Conference on Learning Theory_ , 2015
  * [5] Guilherme França, Michael I Jordan and René Vidal  “On Dissipative Symplectic Integration with Applications to Gradient-Based Optimization”  In _arXiv:2004.06840_ , 2020
  * [6] Mark Konstantinovich Gavurin  “Nonlinear functional equations and continuous analogues of iteration methods”  In _Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika_ Kazan (Volga region) Federal University, 1958
  * [7] Euhanna Ghadimi, Hamid Reza Feyzmahdavian and Mikael Johansson  “Global convergence of the heavy-ball method for convex optimization”  In _Control Conference (ECC), 2015 European_ , 2015  IEEE
  * [8] Eduard Gorbunov et al.  “A stochastic derivative free optimization method with momentum”  In _arXiv:1905.13278_ , 2019
  * [9] Ernst Hairer  “Challenges in geometric numerical integration”  In _Trends in contemporary mathematics_ Springer, 2014
  * [10] Ernst Hairer and Christian Lubich  “Long-term analysis of the Störmer–Verlet method for Hamiltonian systems with a solution-dependent high frequency”  In _Numerische Mathematik_ Springer, 2016
  * [11] Ernst Hairer, Christian Lubich and Gerhard Wanner  “Geometric numerical integration illustrated by the Störmer–Verlet method”  In _Acta numerica_ Cambridge University Press, 2003
  * [12] Hassan K Khalil and Jessy W Grizzle  “Nonlinear systems”  Prentice hall, 2002
  * [13] Diederik P Kingma and Jimmy Ba  “Adam: A method for stochastic optimization”  In _arXiv:1412.6980_ , 2014
  * [14] Walid Krichene, Alexandre Bayen and Peter L Bartlett  “Accelerated mirror descent in continuous and discrete time”  In _Advances in Neural Information Processing Systems 28_ Curran Associates, Inc., 2015
  * [15] Laurent Lessard, Benjamin Recht and Andrew Packard  “Analysis and design of optimization algorithms via integral quadratic constraints”  In _SIAM Journal on Optimization_ SIAM, 2016
  * [16] Michael Muehlebach and Michael I Jordan  “Optimization with momentum: dynamical, control-theoretic, and symplectic perspectives”  In _arXiv:2002.12493_ , 2020
  * [17] Nemirovskii and Yudin  “Problem Complexity and Efficiency of Optimization Methods”  Nauka, Moscow, 1979
  * [18] Yurii Nesterov  “Lectures on convex optimization”  Springer, 2018
  * [19] Yurii E Nesterov  “A method for solving the convex programming problem with convergence rate $O(1/k^{2})$”  In _Dokl. Akad. Nauk SSSR_ , 1983
  * [20] Boris T Polyak  “Some methods of speeding up the convergence of iteration methods”  In _USSR Computational Mathematics and Mathematical Physics_ No longer published by Elsevier, 1964
  * [21] Bin Shi, Simon S Du, Michael I Jordan and Weijie J Su  “Understanding the acceleration phenomenon via high-resolution differential equations”  In _Mathematical Programming_ Springer, 2022
  * [22] Weijie Su, Stephen Boyd and Emmanuel J. Candès  “A Differential Equation for Modeling Nesterov’s Accelerated Gradient Method: Theory and Insights”  In _Journal of Machine Learning Research_ , 2016
  * [23] Jun-Kun Wang, Chi-Heng Lin, Andre Wibisono and Bin Hu  “Provable acceleration of heavy ball beyond quadratics for a class of Polyak-Lojasiewicz functions when the non-Convexity is averaged-out”  In _International Conference on Machine Learning_ , 2022
  * [24] Ashia C Wilson, Lester Mackey and Andre Wibisono  “Accelerating rescaled gradient descent: Fast optimization of smooth functions”  In _Advances in Neural Information Processing Systems_ , 2019
  * [25] Pan Xu, Tianhao Wang and Quanquan Gu  “Accelerated stochastic mirror descent: From continuous-time dynamics to discrete-time algorithms”  In _International Conference on Artificial Intelligence and Statistics_ , 2018


# DFT+DMFT study of the magnetic susceptibility and the correlated electronic
structure in transition-metal intercalated NbS2

Hyowon Park1,2 and Ivar Martin1 1Materials Science Division, Argonne National
Laboratory, Argonne, IL, 60439, USA,
2Department of Physics, University of Illinois at Chicago, Chicago, IL 60607,
USA

###### Abstract

The Co-intercalated NbS2 (Co1/3NbS2) compound exhibits large anomalous Hall
conductance, likely due to the non-coplanar magnetic ordering of Co spins. In
this work, we study the relation between this novel magnetism and the
correlated electronic structure of Co1/3NbS2 by adopting dynamical mean field
theory (DMFT) to treat the correlation effect of Co $d$ orbitals. We find that
the hole doping of Co1/3NbS2 can tune the size of the Nb hole pocket at the
DMFT Fermi surface, producing features consistent with those observed in angle
resolved photoemission spectra [Phys. Rev. B 105, L121102 (2022)]. We also
compute the momentum-resolved spin susceptibility, and correlate it with the
Fermi surface shape. We find that the magnetic ordering wavevector of
Co1/3NbS2 obtained from the peak in spin susceptibility agrees with the one
observed experimentally by neutron scattering; it is compatible with
commensurate non-coplanar $3q$ spin structure. We also discuss how results
change if some other than Co transition metal intercalations are used.

##  I Introduction

Understanding the relation between a novel electronic transport and the
correlated electronic structure of complex materials has been a grand
challenge in condensed matter physics. For instance, cobalt-intercalated NbS2
(Co1/3NbS2) shows a very large anomalous Hall effect [1, 2]; however, the
origin of this phenomenon has remained a subject of debate. While the
‘standard’ anomalous Hall effect originates stems from finite uniform
magnetization in ferromagnets [3], the antiferromagnetic (AFM) ground state of
Co1/3NbS2 implies a different, more exotic origin. One promising scenario is
that topologically non-trivial magnetism such as the non-coplanar spin state
of Co $d$ orbitals can generate a strong fictitious magnetic field, thus
inducing a “topological” Hall effect [4]. Equivalently, in the band language,
the Co spin moments couple to the itinerant Nb bands and produce the
topologically non-trivial band structure resulting in a large Berry curvature
leading to anomalously large Hall currents.

Our previous first-principle calculations based on density functional theory
(DFT) support this scenario: the energy of the non-coplanar $3q$ magnetic
structure is the lowest, compared to other $1q$ or $2q$ states [5]. Moreover,
the Berry phase calculation based on the magnetic band structure supports a
large anomalous Hall conductivity (AHC), comparable to $e^{2}/h$ per
crystalline layer [5]. Unfortunately, the DFT band structure of Co1/3NbS2
fails to capture the angle-resolved photoemission spectra (ARPES), which
requires going beyond the rigid-band shift picture due to the Co intercalation
[6, 7]. One important feature observed in ARPES measurements in Co1/3NbS2 is
the appearance of the broad electron pocket around the high symmetry $K$
point, which is not captured by DFT [6, 7, 8]. Moreover, the effective
electron mass at the electron pocket of Co1/3NbS2 is twice larger than that of
NbS2 [7]. These suggest that for a more accurate picture, one needs to treat
the strong correlation effect of Co $d$ orbitals beyond DFT.

Until recently, experimental support for exotic magnetic states in Co1/3NbS2
had been lacking. In fact, early neutron scattering measurement on Co1/3NbS2
argued that the scattering peak data fits well to the standard commensurate
($1q$) AFM structure, but with multiple magnetic domains [9]. (One should
note, that it is quite difficult to distinguish between multi-domain $1q$ AFM
state and a mono-domain $3q$ state in neutron scattering.) Recently, however,
polarized neutron scattering measurements on the related material Co1/3TaS2,
which has the same structure but with Nb replaced by Ta, convincingly
demonstrated the presence of non-coplanar Co magnetism. Moreover, they
demonstrated connection between the appearance of this magnetic state and
large spontaneous topological Hall effect [10].

What is the physical origin of noncoplanar magnetism? In pure-spin models it
typically requires having multi-spin interactions (either four or six spin).
When interactions are mediated by itinerant electrons, such higher order terms
are generated naturally [11]. From the weak-coupling perspective, having
multiple $q$ orders present simultaneously allows to gap out larger total
sections of Fermi. The key ingredients for noncoplanar order are (1) the large
susceptibility with respect to simple collinear single-$q$ order (e.g. if $q$
connect nearly flat opposite sides of the Fermi surface – nesting effect), and
(2) at least approximate commensuration of the magnetic order and the crystal
lattice. The latter criterion allows for several $q$ orderings to coexist with
each other without suppressing the local amplitude of the magnetic order.

In this work, we focus on the first element, the analysis of magnetic
susceptibility. We go beyond the standard DFT by adopting dynamical mean field
theory (DMFT), which allows to treat the strong correlation effect on the
transition metal ($M$) $d$ orbitals in $M_{1/3}$NbS2. As the first step, we
match the main features of DMFT Fermi surface calculations to be consistent
with the ARPES data. We then calculate the momentum-dependent magnetic
susceptibility $\chi$ from first-principles and investigate the momentum $q$
vector showing the leading instability. We note that DFT alone is not
sufficient to study the correlated electronic structure of Co1/3NbS2 as it
fails to capture the essential features of the APRES measurement. Also,
compared to our previous DFT study on the Co1/3NbS2 that showed $3q$ non-
coplanar structure to be the lowest in energy compared to possible $1q$ or
$2q$ states, now we are allowing for the possibility of instability at a
wavevectors incommensurate with the lattice.

##  II Methods

In this section, we explain computational methods used in the band structure
and the magnetic susceptibility calculations of $M$Nb3S6 ($M$= Co, Fe, and
Ni). We also provide parameters used in the calculations.

###  II.1 DMFT calculation

To study the band structure and the Fermi surface of $M$Nb3S6 ($M$= Co, Fe,
and Ni), we adopt DFT+DMFT treating the strong correlation effect of $M$ ions.
The procedure of the DMFT calculation is as follows. First, we obtain the non-
spin-polarized (nsp) band structure from the experimental $M$Nb3S6 crystal
structures. We adopted the Vienna Ab-initio Simulation Package (VASP) [12, 13]
code to compute the nsp band structure using a $14\times 14\times 4$ $k-$mesh
along with the energy cutoff of 400eV for the plane-wave basis. We used the
Perdew-Burke-Ernzerhof (PBE) functional for the exchange and correlation
energy of DFT. Using the nsp band structure, we construct the following tight-
binding Hamiltonian by adopting the maximally localized Wannier function [14]
as the basis,

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle\sum_{\alpha\beta,ij\sigma}t_{\alpha\beta,ij}\hat{c}_{i\alpha\sigma}^{\dagger}\hat{c}_{j\beta\sigma}+\sum_{\alpha\beta,i}\sum_{\sigma\sigma^{\prime}}U^{\sigma\sigma^{\prime}}_{\alpha\beta}\hat{n}_{i\alpha\sigma}\hat{n}_{i\beta\sigma^{\prime}}$
$\displaystyle+\sum_{ij\sigma}t_{ij}\hat{d}_{i\sigma}^{\dagger}\hat{d}_{j\sigma}+\sum_{\alpha,ij\sigma}(t_{\alpha,ij}\hat{c}_{i\alpha\sigma}^{\dagger}\hat{d}_{j\sigma}+t^{*}_{\alpha,ij}\hat{d}_{j\sigma}^{\dagger}\hat{c}_{i\alpha\sigma}),$

where $t$ is the inter-orbital and inter-site hopping matrix elements
including the whole manifold of the Co $d$ orbitals
($\hat{c}^{\dagger},\hat{c}$) and the Nb $d_{z^{2}}$ orbital
($\hat{d}^{\dagger},\hat{d}$).
$\hat{n}_{i\alpha\sigma}(=\hat{c}^{\dagger}_{i\alpha\sigma}\hat{c}_{i\alpha\sigma})$
is the density operator for the orbital $\alpha$ and the spin $\sigma$ at the
site $i$. $U^{\sigma\sigma^{\prime}}_{\alpha\beta}$ is the local Coulomb
interaction matrix for the on-site Co $d$ orbitals and approximated as the
density-density interaction type.

Using the Hamiltonian in Eq. II.1, we solve the DMFT self-consistent equations
[15] using the continuous-time quantum Monte Carlo (CTQMC) method [16] as the
impurity solver, then obtain the the local self-energy $\Sigma(i\nu_{n})$ for
the Co $d$ orbitals. Here, we parameterize the $U_{\alpha\beta}$ matrix
elements by the Slater integrals using the local Hubbard interaction $U$=5eV
and the Hund’s coupling $J$=0.7eV. The temperature $T$ is set to be 116K. In
DMFT, we use a fine $k-$mesh of 30$\times$30$\times$10\. For all compounds,
the fixed double counting potential scheme is adopted using the following
formula to subtract the double-counting correction from the DFT potential

$V_{DC}=U(n_{d}^{0}-\frac{1}{2})-\frac{J}{2}(n_{d}^{0}-1)$ (2)

where $n_{d}^{0}$ is the nominal occupancy of the transition metal $d$
orbitals, i.e. $n_{d}^{0}=7.0$ for Co2+, $n_{d}^{0}=6.0$ for Fe2+, and
$n_{d}^{0}=8.0$ for Ni2+.

###  II.2 Spin susceptibility $\chi^{sp}(\mathbf{q})$ calculation

A general form of the spin susceptibility $\chi^{sp}$ can be given by the
retarded two-particle Green’s function of two spin operators:

$\chi_{mn}^{sp}(\mathbf{r}^{\prime}-\mathbf{r},t^{\prime}-t)=\Theta(t^{\prime}-t)\langle[\hat{S}^{m}(\mathbf{r}^{\prime},t^{\prime}),\>\hat{S}^{n}(\mathbf{r},t)]_{-}\rangle,$
(3)

where ${\hat{S}^{n}(\mathbf{r},t)}$ is the $n$-th component of spin density
and $\Theta(t^{\prime}-t)$ is the step function that imposes causality. The
spin operator can be expanded using the localized orbital basis set and the
fermionic creation/annihilation operators:

$\hat{S}^{n}(\mathbf{r},t)=\frac{1}{N}\sum_{\sigma\sigma^{\prime}}\hat{s}^{n}_{\sigma\sigma^{\prime}}\sum_{nn^{\prime}}\phi^{*}_{n}(\mathbf{r})\phi_{n^{\prime}}(\mathbf{r})\>\hat{c}^{\dagger}_{n\sigma}(t)\hat{c}_{n^{\prime}\sigma^{\prime}}(t),$
(4)

where $\hat{s}^{i}=g\mu_{B}\hat{\sigma}^{i}/2$ with $\hat{\sigma}^{i}$ being
the $i$-th Pauli matrix and $\phi_{n}(\mathbf{r})$ is the basis function with
the index $n\>(=\\{\mathbf{k},\alpha,\tau_{\alpha}\\})$ for the orbital
character $\alpha$ located at the position $\tau_{\alpha}$ with the momentum
$\mathbf{k}$. One can note that the spin operator can be diagonal for both
spin and orbital basis sets if the spin arrangement is collinear and the spin-
orbit coupling is neglected. Here, we consider the paramagnetic spin symmetry
without the spin-orbit coupling.

For the longitudinal and paramagnetic spin symmetries ($m=n=z$), the spin
susceptibility $\chi_{zz}^{sp}$ can be obtained from

$\displaystyle\chi_{zz}^{sp}(\mathbf{r}^{\prime}-\mathbf{r},t^{\prime}-t)$
$\displaystyle=$
$\displaystyle\frac{\Theta(t^{\prime}-t)}{N^{2}}\cdot\sum_{\sigma\sigma^{\prime}}\sum_{\bar{\sigma}\bar{\sigma}^{\prime}}\hat{s}^{z}_{\sigma\sigma^{\prime}}\hat{s}^{z}_{\bar{\sigma}\bar{\sigma}^{\prime}}$
$\displaystyle\sum_{nn^{\prime}}\sum_{\bar{n}\bar{n}^{\prime}}\phi^{*}_{n}(\mathbf{r^{\prime}})\phi_{n^{\prime}}(\mathbf{r^{\prime}})\phi^{*}_{\bar{n}}(\mathbf{r})\phi_{\bar{n}^{\prime}}(\mathbf{r})$
$\displaystyle\cdot\langle[\hat{c}^{\dagger}_{n\sigma}(t^{\prime})\hat{c}_{n^{\prime}\sigma^{\prime}}(t^{\prime}),\>\hat{c}^{\dagger}_{\bar{n}\bar{\sigma}}(t)\hat{c}_{\bar{n}^{\prime}\bar{\sigma}^{\prime}}(t)]_{-}\rangle.$

Here, the paramagnetic symmetry imposes that the two-particle response
function should be invariant upon the spin flip, i.e.
$\sigma\rightarrow-\sigma$.

###  II.3 Form factor $F$ calculation

Using the paramagnetic symmetry, the momentum and frequency dependent
susceptibility, $\chi_{zz}^{sp}(\mathbf{q},\omega)$ can be simplified from the
Fourier transform of Eq. II.2:

$\chi_{zz}^{sp}(\mathbf{q},\omega)=\frac{(g\mu_{B})^{2}}{2N^{2}}\sum_{nn^{\prime}}\sum_{\bar{n}\bar{n}^{\prime}}F_{nn^{\prime}\bar{n}\bar{n}^{\prime}}(\mathbf{q})\cdot\chi_{nn^{\prime}\bar{n}\bar{n}^{\prime}}\left(\omega\right),$
(6)

where $F(\mathbf{q})$ is the atomic form factor describing the modulation of
the charge density and $\chi(\omega)$ is the orbital-dependent two-particle
response function where the spin dependence is simplified due to the
paramagnetic symmetry. In the susceptibility calculation based on DFT, the
matrix element for $F(\mathbf{q})$ can be typically computed using the Kohn-
Sham wavefunction, i.e. $\psi_{n_{\alpha}}^{\mathbf{k}}(\mathbf{r})$ (the
orbital index $\\{\alpha,\tau_{\alpha}\\}$ changes to the band index
$n_{\alpha}$), and it can be expanded continuum basis set, such as plane waves
[17]. Here, we adopt the maximally localized Wannier function for the form
factor $F(\mathbf{q})$ calculation, which is the same basis function for DMFT
calculations:

$\displaystyle\phi_{n}(\mathbf{r})=\phi^{\mathbf{k}}_{n_{\alpha}}(\mathbf{r})=\frac{1}{\sqrt{N}}\sum_{\mathbf{R}}e^{-i\mathbf{k}\cdot\mathbf{R}}\phi^{\tau_{\alpha}}_{\alpha}(\mathbf{r}-\mathbf{R}),$
(7)

where $\phi^{\mathbf{k}}_{n_{\alpha}}(\mathbf{r})$ is the Wannier function
with the index $n_{\alpha}\>(=\\{\alpha,\tau_{\alpha}\\})$ defined in the
$\mathbf{k}-$space for the primitive unit cell. The index $n_{\alpha}$ runs
over both the orbital character $\alpha$ and the internal atomic position
$\tau_{\alpha}$.

If the complete and orthonormalized basis set is used for the $F(\mathbf{q})$
calculation, one can obtain the product of delta functions imposing the
momentum and orbital conservation:

$\displaystyle F_{nn^{\prime}\bar{n}\bar{n}^{\prime}}(\mathbf{q})$
$\displaystyle=$
$\displaystyle\delta_{\mathbf{k}+\mathbf{q},\mathbf{k^{\prime}}}\cdot\delta_{\bar{\mathbf{k}}^{\prime}+\mathbf{q},\bar{\mathbf{k}}}\cdot\delta_{n_{\alpha},n_{\alpha^{\prime}}}\cdot\delta_{n_{\bar{\alpha}},n_{\bar{\alpha}^{\prime}}}.$
(8)

However, in general, $F(\mathbf{q})$ needs to be modified if one uses the
incomplete basis set. In the case of Co1/3NbS2, the magnetic moments primarily
reside on Co $d$ orbitals, which are the subset of the complete band
structure. As a result, the $F(\mathbf{q})$ expression for only Co $d$
orbitals can be modified as follows:

$\displaystyle F_{nn^{\prime}\bar{n}\bar{n}^{\prime}}(\mathbf{q})$
$\displaystyle=$
$\displaystyle\frac{1}{N^{2}_{\tau}}\biggl{(}\>\sum_{\mathbf{\tilde{k}},\mathbf{\tilde{k}}^{\prime}}\sum_{\vec{\tau}_{\alpha},\vec{\tau}_{\alpha^{\prime}}}e^{i(\tilde{\mathbf{k}}^{\prime}\cdot\vec{\tau}_{\alpha^{\prime}}-\tilde{\mathbf{k}}\cdot\vec{\tau}_{\alpha})}$
(9)
$\displaystyle\cdot\delta_{\alpha\alpha^{\prime}}\cdot\delta_{\\{\mathbf{k}\\},\tilde{\mathbf{k}}}\cdot\delta_{\\{\mathbf{k}^{\prime}\\},\tilde{\mathbf{k}}^{\prime}}\cdot\delta_{\tilde{\mathbf{k}}+\\{\mathbf{q}\\},\tilde{\mathbf{k}}^{\prime}}\biggr{)}$
$\displaystyle\cdot\biggl{(}\>\sum_{\mathbf{\tilde{k}},\mathbf{\tilde{k}}^{\prime}}\sum_{\vec{\tau}_{\bar{\alpha}},\vec{\tau}_{\bar{\alpha}^{\prime}}}e^{i(\tilde{\mathbf{k}}^{\prime}\cdot\vec{\tau}_{\bar{\alpha}^{\prime}}-\tilde{\mathbf{k}}\cdot\vec{\tau}_{\bar{\alpha}})}$
$\displaystyle\cdot\delta_{\bar{\alpha}\bar{\alpha}^{\prime}}\cdot\delta_{\\{\bar{\mathbf{k}}\\},\tilde{\mathbf{k}}}\cdot\delta_{\\{\bar{\mathbf{k}}^{\prime}\\},\tilde{\mathbf{k}}^{\prime}}\cdot\delta_{\tilde{\mathbf{k}}^{\prime}+\\{\mathbf{q}\\},\tilde{\mathbf{k}}}\biggr{)},$

where the momentum $\tilde{\mathbf{k}}$ is defined in an extended Brillouin
zone (BZ) obtained for a single Co ion in triangular lattice,
$\\{\mathbf{k}\\}$ represents the set of the momentum $\mathbf{k}$ vectors
shifted by the reciprocal vectors, $\vec{\tau}_{\alpha}$ is the atomic
position of correlated atoms, and $N_{\tau}$ is the number of correlated atoms
in the primitive unit cell. More detailed derivation of Eq. 9 is given in the
Appendix.

###  II.4 The Bethe-Salpeter equation

The orbital-dependent susceptibility
$\chi_{nn^{\prime}\bar{n}\bar{n}^{\prime}}(\omega)$ in Eq. 6 can be computed
using the following Bethe-Salpeter equation:

$\chi_{nn^{\prime}\bar{n}\bar{n}^{\prime}}(\omega)=\chi_{n_{\alpha}n_{\bar{\alpha}}}^{0}(\mathbf{k},\mathbf{k^{\prime}},\omega)+\chi^{0}*\Gamma^{irr}*\chi,$
(10)

where $\chi^{0}$ is the polarizability obtained from the interacting Green’s
function and $\Gamma^{irr}$ is the irreducible vertex function. One should
note that the orbital indices for $\chi^{0}_{n_{\alpha}n_{\bar{\alpha}}}$ run
over all Co $d$ and Nb $d_{z^{2}}$ orbitals in a unit cell while those on the
$F(\mathbf{q})$ in Eq. 9 account only the correlated Co $d$ orbitals. In
general, $\Gamma^{irr}$ is a complex function depending on momentum,
frequency, spin, site, and orbital degrees of freedom. While the effect of
$\Gamma^{irr}$ is crucial to compare both momentum and frequency dependence of
the susceptibility to the experimental neutron scattering data [18, 19], we
adopt the static interaction type based on the random phase approximation
(RPA), namely assuming that the interaction matrix to be independent on
momentum and frequency. Here, we further approximate that it is independent of
orbitals to account for the average interaction effect and has the spin
rotation symmetry to consider only interactions in the spin channel. As a
result, we consider both the on-site interaction $\bar{U}$ within Co ions and
the inter-site interaction $\bar{V}$ between Co and Nb ions, then study the
effects of $\bar{U}$ and $\bar{V}$ on the susceptibility calculations. For the
on-site $\bar{U}$ value, we used the $\bar{U}$=1.5eV, which is smaller than
the DMFT $U$ value. This is because the static two-particle interaction is
further renormalized within the RPA diagrams while the diagrams within DMFT
take the orbital and dynamical fluctuations into account explicitly. We also
used different inter-site $\bar{V}$ values ($\bar{V}$=0, 0.2, and 0.3 eV) to
explore the effect of $\bar{V}$ on the susceptibility. Our results show that
the inter-site $\bar{V}$ can enhance the momentum dependence of the
susceptibility significantly.

The polarizability $\chi^{0}$ can be given by the product of two Green’s
functions using the Wick’s theorem. Our $\chi^{0}$ is different from the bare
susceptibility since it is obtained from the interacting Green’s function. In
DMFT, the Green’s function is dressed with the local self-energy and the
Matsubara frequency sum over $i\nu$ can be evaluated by performing the contour
integral to obtain the polarizability at $\omega=0^{+}$:

$\displaystyle\chi^{0\>\mathbf{k},\mathbf{k^{\prime}}}_{n_{\alpha}n_{\bar{\alpha}}}$
$\displaystyle=$
$\displaystyle-T\sum_{i\nu}G_{n_{\alpha}n_{\bar{\alpha}}}(\textbf{k},i\nu)\cdot
G_{n_{\bar{\alpha}}n_{\alpha}}(\textbf{k}^{\prime},i\nu+i0^{+})$
$\displaystyle=$ $\displaystyle-\frac{T}{2\pi i}\oint
dz\>\>G_{n_{\alpha}n_{\bar{\alpha}}}(\textbf{k},z)\cdot
G_{n_{\bar{\alpha}}n_{\alpha}}(\textbf{k}^{\prime},z+i0^{+})$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int
d\nu\>[ImG_{n_{\alpha}n_{\bar{\alpha}}}(\textbf{k},\nu)\cdot
G_{n_{\bar{\alpha}}n_{\alpha}}(\textbf{k}^{\prime},\nu+i0^{+})$
$\displaystyle+G_{n_{\alpha}n_{\bar{\alpha}}}(\textbf{k},\nu-i0^{+})\cdot
ImG_{n_{\bar{\alpha}}n_{\alpha}}(\textbf{k}^{\prime},\nu)]f(\nu),$

where $f(\nu)$ is the Fermi function. The real part of the susceptibility at
$\omega=0^{+}$ is given by

$\displaystyle
Re\chi^{0\>\mathbf{k},\mathbf{k^{\prime}}}_{n_{\alpha}n_{\bar{\alpha}}}$
$\displaystyle=$ $\displaystyle\frac{1}{\pi}\int d\nu
f(\nu)[ImG_{n_{\alpha}n_{\bar{\alpha}}}(\textbf{k},\nu)\cdot
ReG_{n_{\bar{\alpha}}n_{\alpha}}(\textbf{k}^{\prime},\nu)$ (12)
$\displaystyle+ReG_{n_{\alpha}n_{\bar{\alpha}}}(\textbf{k},\nu)\cdot
ImG_{n_{\bar{\alpha}}n_{\alpha}}(\textbf{k}^{\prime},\nu)],$

where both $ReG$ and $ImG$ are the real and imaginary parts of interacting
Green’s functions defined on the fine $\mathbf{k}$ mesh and the real frequency
which are obtained from the analytic continuation of the DMFT self-energy
using the maximum entropy method. For the $\chi^{0}$ susceptibility
calculation in Eq. II.4, we performed the summation over the dense
$\mathbf{k}-$grid using 60$\times$60$\times$10 $\mathbf{k}-$points at each
$\mathbf{q}-$ points chosen along the high symmetry points.

Within static theories such as DFT, equivalent to setting the DMFT self-energy
is zero, $\chi^{0}(\mathbf{q},\omega)$ from Eq. 6 is given as the bare
susceptibility and its evaluation using Eq. II.4 and Eq. 8 reduces to the
Lindhard formula for the susceptibility:

$\displaystyle\chi^{0}(\mathbf{q},\omega)$ $\displaystyle=$
$\displaystyle\frac{(g\mu_{B})^{2}}{2N}\sum_{\textbf{k},nm}\sum_{\alpha\beta}\frac{f^{m}_{\mathbf{k+q}}-f^{n}_{\mathbf{k}}}{\omega+\epsilon^{n}_{\mathbf{k}}-\epsilon^{m}_{\mathbf{k+q}}+i0^{+}}$
$\displaystyle\cdot\langle\phi_{\alpha}^{\mathbf{k}}|\psi_{n}^{\mathbf{k}}\rangle\langle\psi_{n}^{\mathbf{k}}|\phi_{\beta}^{\mathbf{k}}\rangle\langle\phi_{\beta}^{\mathbf{k+q}}|\psi_{m}^{\mathbf{k+q}}\rangle\langle\psi_{m}^{\mathbf{k+q}}|\phi_{\alpha}^{\mathbf{k+q}}\rangle,$

where $\epsilon^{n}_{\mathbf{k}}$ is the eigenvalue of the DFT band at
momentum $\mathbf{k}$ and the band index $n$. Therefore, it is expected that
the $\chi^{0}(\mathbf{q})$ will be enhanced near the Fermi surface nesting
$\mathbf{q}$ vector. This susceptibility calculation for multi-orbital systems
based on DFT has been applied for various real materials [20, 17, 21].

##  III Results and discussions

We now present the results for the correlated electronic band structure and
the Fermi surface, then relate them to the momentum dependent magnetic
susceptibility of $M_{1/3}$NbS2 ($M$=Co, Ni, and Fe) computed using DMFT. In
particular, we study the effect of the hole doping on the electronic structure
and magnetism.

###  III.1 Correlated electronic structure of Co1/3NbS2

Figure 1:  (a) The crystal structure of Co1/3NbS2, (b) The DFT band structure
of Co1/3NbS2 projected to different orbital characters (red: Nb $d_{z^{2}}$
band, blue: Co $d$ band)  Figure 2:  (a) The DMFT band structure of Co1/3NbS2,
(b) The DMFT band structure upon the hole doping $\delta$=2.0 (two holes per
CoNb3S6. The green lines represent the DFT band from Fig. 1(b).

The Co1/3NbS2 structure is formed by stacking NbS2 layers with the
intercalation of Co ions at two distinct Nb sites between the layers along the
$c-$axis (see Fig. 1(a)). To study the strong correlation effect of the Co $d$
orbitals on the band structure, we compare the DMFT spectral function
$A(\mathbf{k},\omega)$ (Fig. 2) to the DFT band structure in Fig. 1(b). Here,
we impose the paramagnetic spin symmetry (nonmagnetic state). The DFT band
(Fig. 1(b), the green thin lines in Fig. 2) shows that the hole pocket near
the $\Gamma$ and $A$ points are quite small and the band crossing near $K$ and
$H$ points occur above the Fermi energy. Fig. 1(b) shows the orbital
characters of DFT band structure and the hole pockets are mostly the Nb
$d_{z^{2}}$ character. The Co $d$ bands (blue color) are mostly located at $K$
and $H$ points above the Fermi energy with multiple degeneracy.

The DMFT bands in Fig. 2(a) show the strong modification of Co $d$ bands near
$K$ and $H$ points as they are dressed by the DMFT self-energy. The location
of Co $d$ bands is pushed below the Fermi energy and the spectra become much
broader due to the large self-energy effect. The Nb $d_{z^{2}}$ bands still
show the quasi-particle dispersion along with the feature of the well-defined
Fermi surface. We also study the doping effect by changing the total number of
valence electrons within the DMFT calculation and computing the corresponding
band structure. The hole-doping effect on the DMFT bands in Fig. 2(b) shows
the crossing of Co $d$ bands near the Fermi energy near $K$ and $H$ points due
to the up-shift of bands. Upon the hole doping, the size of the hole pocket
near the $\Gamma$ point increases as the Nb $d_{z^{2}}$ band shifts upward. An
important feature of the hole-doping on the Co $3d$ band structure is the
appearance of the small broad electron pocket at the $K$ point, which is
consistent with the ARPES measurement [8].

In Co1/3NbS2, our DMFT calculation shows that the occupancy of the Co $d$
orbital is close to 7.0 meaning that the Co ion has the valence state of $2+$.
Since the S ion is the $2-$ state, the valence state of Nb is close to
$(10/3)+$. Therefore, the Nb ion has the occupancy of $4d^{1.67}$, which is
larger than the $4d^{1}$ one (Nb4+) of the pure NbS2 layer. In other words,
the hole pocket at the $\Gamma$ point of the pure NbS2 is expected to be
larger than that of Co1/3NbS2 (Fig. 2(a)). Our hole doping effect of
$\delta$=2.0 (two holes per CoNb3S6) means that 2/3 holes are mostly doped to
the Nb ions since the DMFT occupancy of the Co $d$ orbitals still remains
close to 7.0. As a result, Co1/3NbS2 has the strong hybridization between Co
$3d$ and Nb $4d_{z^{2}}$ orbitals resulting the doped holes residing mostly on
the Nb $d_{z^{2}}$ orbital. Therefore, this hole doping effect mainly affects
the size the hole pocket near the $\Gamma$ point.

Figure 3: The DMFT Fermi surface for Co1/3NbS2 as a function of hole-doping
(a) $\delta$=0, (b) 1.0, and (c) 2.0. The hole-doping effect increases the
size of the hole pocket centered at the $\Gamma$ point. The calculated Fermi
surface is consistent with the (d) experimental ARPES measurement [8] when
$\delta$=2.0. The dashed line represents the BZ of Co1/3NbS2 and the solid
line shows the BZ of NbS2. Note that the $M$ point in (d) is defined in the
solid-line BZ, while our $M$ points are defined in the dashed-line BZ.

Fig. 3 shows the Fermi surface of Co1/3NbS2 computed using DFT+DMFT at
different hole doping $\delta$ values ($\delta$= the number of holes per Co
ion). At $\delta$=0, the DMFT hole-pocket centered at the $\Gamma$ point has a
circular shape similarly as measured in ARPES although its size is smaller
than the ARPES one. The outer larger pocket has mostly Co $d$ orbital
character and exhibits much weaker intensity due to the large scattering rate
($Im\Sigma(\omega)$). Upon hole-doping, the smaller hole-pocket gets larger in
size, comparable to the ARPES measurement and the Co $d$ states move closer to
the $K$ point in the BZ. This broad Co $d$ spectra near the $K$ point is also
captured in ARPES. Our DMFT Fermi surface calculation shows that the hole
doping of $\delta$=2.0 makes the size of Nb hole pocket in the Fermi surface
similar to the ARPES measurement.

Since the ARPES measurement is sensitive to the surface state of possibly NbS2
termination layers, the electronic structure will have the hole-doping effect
due to the Co ion deficiency. In the bulk, Co or Nb vacancies can induce the
effect of hole-dopings similarly to the surface state. Our DMFT calculation
shows that this doping effect can tune the size of hole pocket significantly,
possibly affecting the magnetic properties. Previous experimental study also
shows that the AHE of Co1/3NbS2-x can be dramatically changed by the S
deficiency, which can lead to the similar doping effect [8].

###  III.2 Co1/3NbS2 magnetic susceptibility

Figure 4: The polarizability $\chi^{0}$ (Eq. II.4) for Co1/3NbS2 computed at
different hole doping $\delta$, i.e. $\delta$=0 (top panel) and $\delta$=2.0
(bottom panel).

Our previous DFT calculation on Co1/3NbS2 shows that a $3q-$type magnetic
structure is energetically stable and may be responsible for the large
observed anomalous Hall current in this material. This particular $3q-$type
spin structure corresponds to the non-coplanar spins arrangement with the four
spins within the magnetic unit cell of one Co layer pointing towards four
vertices of a tetrahedron in the spin space. Such unusual spin configuration
can be stabilized in a triangular lattice due to the Fermi surface nesting of
itinerant electrons [4]. For every triangular plaquette of the spin lattice,
scalar spin chirality, $\chi_{ijk}={\bf S}_{i}\cdot[{\bf S}_{j}\times{\bf
S}_{k}]$ is constant, corresponding to a uniform Berry flux per plaquette. The
modulation vector $\mathbf{q}$ of this $3q-$type spin structure is half of the
reciprocal lattice vectors of the primitive unit cell, i.e., the high-symmetry
$M$ points in the Brillouine zone. Early neutron scattering experiment [9]
showed indeed the scattering peak at $\mathbf{q}=(1/2,0,0)$ ($M$ point).
However, confirming the existence of the non-coplanar 3$q$ state requires more
complex polarized neutron scattering experiments, as has been successfully
realized in some cases [22, 10].

Neutron scattering measures magnetic susceptibility. We here compute the
momentum-dependent magnetic susceptibility, $\chi(\mathbf{q},\omega=0)$ of
Co1/3NbS2 at different hole-doping $\delta$ values to understand the spin
modulation vector $\mathbf{q}$ of the leading magnetic instability and its
relation to the correlated electronic structure. We first compute the real
part of the bare susceptibility $\chi^{0}$ using DMFT (Eq. II.4) at $\delta=0$
and 2.0, as shown in Fig. 4. In both doping levels, the Co bare spin
susceptibility ($\chi^{0}$) obtained from DMFT shows very weak momentum
dependence due to the the strongly localized nature of Co $d$ orbitals. This
can be understood from the one-particle DMFT spectra of Co $d$ bands showing
no clear evidence of quasi-particle peaks near the Fermi energy, rather very
broad feature of band dispersion without much dependence on momenta. Unlike
the Co $d$ orbitals, the Nb $d_{z^{2}}$ orbitals show a rather strong momentum
dependence of the susceptibility at $\delta=2.0$ due to their itinerant nature
near the Fermi energy. The contribution of the Nb $d_{z^{2}}$ orbitals to the
susceptibility $\chi^{0}$ is larger than that of the Co $d$ orbitals and
depends sensitively on momentum and the the doping levels (see Fig. 4).

We argue that the leading modulation vector $\mathbf{q}$ of $\chi^{0}$ in
Co1/3NbS2 can be mostly determined by the Fermi surface momentum ($2k_{F}$) of
the hole pocket centered at the $\Gamma$ point. As shown in Fig. 3, the size
of the hole pocket can be sensitively dependent on the hole doping levels and
it is consistent with the ARPES measurement at $\delta=2.0$. At $\delta$=0,
the Nb $d_{z^{2}}$ orbital contribution to $\chi^{0}$ has no clear momentum
dependence although the spectra near the $\Gamma$ point are slightly larger
than those at other momenta. This is because the smallest hole-pocket near the
$\Gamma$ point has the largest spectral weight while the other Fermi surfaces
have much smaller spectral weights. As the hole doping increases
($\delta$=2.0), the size of the hole pocket near the $\Gamma$ point also
increases and the contribution of the Nb $d_{z^{2}}$ orbital to the
susceptibility favors the modulation $\mathbf{q}$ vector at the high-symmetry
$M$ point, which is close to the Nb Fermi surface momentum (2$k_{F}$). The
susceptibility near the $K$ point also shows the enhanced peak height at
$|\mathbf{q}|\simeq 2k_{F}$ although the $M$ point shows the maximum peak
height. We note that the peak heights of the susceptibility depend on the
effect of the form factor in Eq. 9 – the peak height at the $M$ point becomes
much closer to that at the $K$ point if the form factor is simplified using
Eq. 8 (see Appendix).

Figure 5: The magnetic susceptibility $\chi$ for Co1/3NbS2 including the RPA-
type interaction at different hole-dopings $\delta=0$ (top panel) and
$\delta=2.0$ (bottom panel).

While the bare magnetic susceptibility $\chi^{0}$ of Co $d$ orbitals does not
have any significant momentum dependence due to the strongly localized nature
of the band structure, the full magnetic susceptibility including the
interaction effect shows the preference for a particular momentum suggesting
the long-range spin ordering of Co $d$ spins coupled via the RKKY interaction
mediated by the itinerant Nb $d_{z^{2}}$ bands. It turns out that the inter-
site interaction $\bar{V}$ plays an important role in mediating the RKKY
interaction. Our RPA susceptibility calculation shows that the local
interaction $\bar{U}$ enhances the absolute value of $\chi(\mathbf{q})$ while
retaining the weak momentum dependence. The increase of $\bar{V}$ results in
the momentum dependence of $\chi(\mathbf{q})$, which is peaked at the $M$
point (the same peak position as $\chi^{0}$) for the RPA susceptibility (see
Fig. 5).

We find that the $\chi$ diverges at the $M$ point near $\bar{U}$=1.7eV and
$\bar{V}$=0.3eV, supporting the occurrence of the magnetic instability. While
this is a direct way to study the instability from the susceptibility, one can
also further analyze the different band contributions to the magnetic
instability by decomposing the product of the $\chi^{0}$ and the
$\Gamma^{irr}$ matrices while solving the Bethe-Salpeter equation in Eq. 10
[21]. Again, the reasonable range of the RPA $\bar{U}$ should be much smaller
than the Hubbard $U$ used in DMFT since it does not account for the orbital
and dynamical screening process. The screened inter-site $\bar{V}$ also should
be much smaller than the on-site $\bar{U}$ value. While determining $\bar{U}$
and $\bar{V}$ quantitatively will be a complicated task, we find that the
qualitative feature of the magnetic susceptiblity (i.e., the momentum
dependence) does not vary depending on $\bar{U}$ and $\bar{V}$ values.

Figure 6: The DMFT Fermi surfaces (left panels) for $M_{1/3}$NbS2 ($M$= (a) Fe
(top panel) and (c) Ni (bottom panel)) at the hole doping of $\delta$= 2.0.
The experimental ARPES [8] measurements are also compared (right panels).

###  III.3 Electronic structure of Fe1/3NbS2 and Ni1/3NbS2

While Co1/3NbS2 shows large anomalous Hall effect likely originating from the
non-coplanar spin structure, such effects have not been seen experimentally
for Fe1/3NbS2 and Ni1/3NbS2. Although our previous DFT calculation showed that
both Fe1/3NbS2 and Ni1/3NbS2 can favor the non-coplanar $3q$ spin structure
energetically, the ground-state magnetic state has been studied only for a
small number of possible $\mathbf{q}$ vectors allowed within a supercell.
Therefore, it is plausible that the leading magnetic $\mathbf{q}-$vector can
vary depending on the intercalated transition metal ions due to the electronic
structure change. Here, we compute the magnetic susceptibility and the Fermi
surface of Fe1/3NbS2 and Ni1/3NbS2 at $\delta=2.0$, similarly as the Co1/3NbS2
case.

We find that the hole doping effects in Fe1/3NbS2 and Ni1/3NbS2 can be quite
different. Fig. 6 shows that the DMFT Fermi surface of Fe1/3NbS2 has a
slightly smaller hole pocket compared to the Co1/3NbS2 one at the same hole
doping ($\delta$=2.0). Our DMFT calculation shows that the hole doping induces
the change of the Fe $d$ occupancy as the Fe valence state becomes close to
Fe2.3+. In Fe1/3NbS2, the hole doping mostly affects the Fe states near the
Fermi energy since the hybridization between the Fe $3d$ and Nb $4d$ orbitals
is rather weak. As a result, the size of the Nb hole pocket in Fe1/3NbS2 is
less sensitive to dopings while it becomes larger upon hole doping for the
Co1/3NbS2 case. Moreover, the Fe $d$ character becomes much weaker and not
visible near the $K$ point, as also seen in the ARPES data [8]. In Ni1/3NbS2,
the DMFT valence of the Ni ion is still close to Ni2+ and the Ni band is
located lower than the Co or Fe bands in the other materials. This means that
the Ni band has the negative charge-transfer effect, similarly to what is
observed in some other rare-earth nickelates [23]. As a result, the hole
doping mostly affects the Nb hole states and the Nb hole pocket in Ni1/3NbS2
becomes the largest among three materials.

Figure 7: The magnetic susceptibility for Fe1/3NbS2, the polarizability
$\chi^{0}$ (top panel) and the RPA $\chi$ (bottom panel)

Figure 8: The magnetic susceptibility for Ni1/3NbS2, the polarizability
$\chi^{0}$ (top panel) and the RPA $\chi$ (bottom panel)

###  III.4 The magnetic susceptibility of Fe1/3NbS2 and Ni1/3NbS2

These notable variations of the DMFT Fermi surface due to the intercalation by
different transition metals can also change the leading magnetic
susceptibility momentum $\mathbf{q}$. In Fe1/3NbS2, the susceptibility
$\mathbf{q}$ peaks near $M$ and $K$ points are almost degenerate as the the
Fermi momentum $2k_{F}$ of the hole pocket in the Fermi surface gets closer to
both $M$ and $K$ points. Similar to the Co1/3NbS2 case, the inter-site
interaction $V$ strongly enhances the momentum dependence of the Fe $d$
susceptibility. In Ni1/3NbS2, the susceptibility peak is slightly higher at
the $K$ point as the leading modulation vector, while the momentum dependence
of the susceptibility is the weakest among three materials. This is because
the $2k_{F}$ of the hole pocket in Ni1/3NbS2 is much larger than the high
symmetry points and, as a result, the susceptibility does not show the
dominant momentum peak.

In both cases of Fe1/3NbS2 and Ni1/3NbS2, the modulation factor (Eq. 14) in
the susceptibility can change the $\chi(\mathbf{q})$ profile. Without the
factor, both susceptibilities enhance the peak near the $K$ point (see
Appendix). It is possible that both Fe1/3NbS2 and Ni1/3NbS2 can have different
magnetic $\mathbf{q}$ instability compared to the Co1/3NbS2 case. Finally, the
spectral weight of the susceptibility is the smallest for Ni1/3NbS2 as the
static magnetic moment of the Ni ion is the smallest compared to the other
materials [24].

##  IV Conclusion

We studied the magnetic susceptibility and the correlated electronic structure
of $M_{1/3}$NbS2 ($M$= Co, Fe, and Ni) using DMFT to treat the strong
correlation effect of transition metal ions. Our DMFT band structure and the
Fermi surface calculations upon hole-doping are consistent with the ARPES
measurements [8] of these compounds. The size of the hole pocket centered at
the $\Gamma$ point is comparable to the ARPES data for all compounds and the
appearance of the electron pocket at the $K$ point in Co1/3NbS2 is correctly
captured. Due to the strong hybridization between Co/Ni $3d$ orbitals and Nb
$4d$ orbitals, the doped holes mostly change the Nb valence states and the
size of the hole pocket centered at the $\Gamma$ point can be tuned upon the
hole doping. In Fe1/3NbS2, the hole-doping effect changes mostly the Fe
valence state.

We also show that the spin susceptibility $\chi(\mathbf{q})$ calculation using
DMFT can help to identify the momentum $\mathbf{q}$ of the leading magnetic
instability in strongly correlated materials. This method allows to avoid the
need for constructing a large supercell to study the magnetic instability at
the arbitrary momentum $\mathbf{q}$ vector. While the spin susceptibility of
Co1/3NbS2 is peaked at $\mathbf{q}=(1/2,0,0)$ (the $M$ point), which is
consistent with the $3q$-type non-coplanar spin structure, the maximum peak
positions change to the $K$ point for Fe1/3NbS2 and Ni1/3NbS2. This suggests
that the magnetic ground state of these two compounds will be distinct from
that of Co1/3NbS2.

## Acknowledgement

We would like to thanks Mike Norman and Chris Lane for useful discussions.
This work was supported by the Materials Sciences and Engineering Division,
Basic Energy Sciences, Office of Science, US Department of Energy. We
gratefully acknowledge the computing resources provided on Bebop, a high-
performance computing cluster operated by the Laboratory Computing Resource
Center at Argonne National Laboratory.

## References

  * Ghimire  _et al._ [2018] N. J. Ghimire, A. S. Botana, J. S. Jiang, J. Zhang, Y.-S. Chen, and J. F. Mitchell, Nature Communications 9, 3280 (2018).
  * 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, Phys. Rev. Res. 2, 023051 (2020).
  * Nagaosa  _et al._ [2010] N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Rev. Mod. Phys. 82, 1539 (2010).
  * Martin and Batista [2008] I. Martin and C. D. Batista, Phys. Rev. Lett. 101, 156402 (2008).
  * Park  _et al._ [2022] H. Park, O. Heinonen, and I. Martin, Phys. Rev. Mater. 6, 024201 (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, 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š, Phys. Rev. B 105, 155114 (2022).
  * 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, Phys. Rev. B 105, L121102 (2022).
  * Parkin  _et al._ [1983] S. S. P. Parkin, E. A. Marseglia, and P. J. Brown, J. Phys. C: Solid State Phys. 16, 2765 (1983).
  * Takagi  _et al._ [2023] H. Takagi, R. Takagi, S. Minami, T. Nomoto, K. Ohishi, M.-T. Suzuki, Y. Yanagi, M. Hirayama, N. D. Khanh, K. Karube, H. Saito, D. Hashizume, R. Kiyanagi, Y. Tokura, R. Arita, T. Nakajima, and S. Seki, Nature Physics 19, 961 (2023).
  * Batista  _et al._ [2016] C. D. Batista, S.-Z. Lin, S. Hayami, and Y. Kamiya, Reports on Progress in Physics 79, 084504 (2016).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Marzari and Vanderbilt [1997] N. Marzari and D. Vanderbilt, Physical review B 56, 12847 (1997).
  * Singh  _et al._ [2021] V. Singh, U. Herath, B. Wah, X. Liao, A. H. Romero, and H. Park, Computer Physics Communications 261, 107778 (2021).
  * Haule [2007] K. Haule, Phys. Rev. B 75, 155113 (2007).
  * Heil  _et al._ [2014] C. Heil, H. Sormann, L. Boeri, M. Aichhorn, and W. von der Linden, Phys. Rev. B 90, 115143 (2014).
  * Park  _et al._ [2011] H. Park, K. Haule, and G. Kotliar, Phys. Rev. Lett. 107, 137007 (2011).
  * Goremychkin  _et al._ [2018] E. A. Goremychkin, H. Park, R. Osborn, S. Rosenkranz, J.-P. Castellan, V. R. Fanelli, A. D. Christianson, M. B. Stone, E. D. Bauer, K. J. McClellan, D. D. Byler, and J. M. Lawrence, Science 359, 186 (2018).
  * Johannes  _et al._ [2006] M. D. Johannes, I. I. Mazin, and C. A. Howells, Phys. Rev. B 73, 205102 (2006).
  * Lane  _et al._ [2023] C. Lane, R. Zhang, B. Barbiellini, R. S. Markiewicz, A. Bansil, J. Sun, and J.-X. Zhu, Communications Physics 6, 90 (2023).
  * Takagi  _et al._ [2018] R. Takagi, J. S. White, S. Hayami, R. Arita, D. Honecker, H. M. Rønnow, Y. Tokura, and S. Seki, Science Advances 4, eaau3402 (2018).
  * Bisogni  _et al._ [2016] V. Bisogni, S. Catalano, R. J. Green, M. Gibert, R. Scherwitzl, Y. Huang, V. N. Strocov, P. Zubko, S. Balandeh, J.-M. Triscone, G. Sawatzky, and T. Schmitt, Nature Communications 7, 13017 (2016).
  * Anzenhofer  _et al._ [1970] K. Anzenhofer, J. van den Berg, P. Cossee, and J. Helle, J. Phys. Chem. Solids 31, 1057 (1970).
  * Ku  _et al._ [2010] W. Ku, T. Berlijn, and C.-C. Lee, Phys. Rev. Lett. 104, 216401 (2010).

## Appendix: The derivation and effect of the Form factor

The Form factor $F(\mathbf{q})$ in Eq. 6 can be given by

$\displaystyle F_{nn^{\prime}\bar{n}\bar{n}^{\prime}}(\mathbf{q})$
$\displaystyle=$ $\displaystyle\int d\mathbf{r}\int
d\mathbf{r^{\prime}}e^{i\mathbf{q}(\mathbf{r^{\prime}}-\mathbf{r})}\phi^{*}_{n}(\mathbf{r^{\prime}})\phi_{n^{\prime}}(\mathbf{r^{\prime}})\phi^{*}_{\bar{n}}(\mathbf{r})\phi_{\bar{n}^{\prime}}(\mathbf{r})$
(14) $\displaystyle=$
$\displaystyle\langle\phi^{\mathbf{k}}_{n_{\alpha}}|e^{i\mathbf{q}\cdot\hat{\mathbf{r}}^{\prime}}|\phi^{\mathbf{k^{\prime}}}_{n_{\alpha^{\prime}}}\rangle\langle\phi^{\mathbf{\bar{k}}}_{n_{\bar{\alpha}}}|e^{-i\mathbf{q}\cdot\hat{\mathbf{r}}}|\phi^{\bar{\mathbf{k}^{\prime}}}_{n_{\bar{\alpha}^{\prime}}}\rangle,$

where the $e^{i\mathbf{q}\cdot\hat{\mathbf{r}}}$ operator can be typically
simplified using the dipole approximation or the constant matrix element
approximation ($e^{i\mathbf{q}\cdot\hat{\mathbf{r}}}\simeq 1$).

For a localized orbital $\phi_{n}(\mathbf{r})$, such as the Wannier function,
the expectation value of the $e^{i\mathbf{q}\cdot\hat{\mathbf{r}}}$ factor can
be performed using the discrete sum of the center of orbitals.

$\displaystyle F_{nn^{\prime}\bar{n}\bar{n}^{\prime}}(\mathbf{q})$
$\displaystyle\simeq$
$\displaystyle\langle\phi^{\mathbf{k}}_{n_{\alpha}}|e^{i\mathbf{q}\cdot\hat{\mathbf{R}}^{\prime}}|\phi^{\mathbf{k^{\prime}}}_{n_{\alpha^{\prime}}}\rangle\langle\phi^{\mathbf{\bar{k}}}_{n_{\bar{\alpha}}}|e^{-i\mathbf{q}\cdot\hat{\mathbf{R}}}|\phi^{\bar{\mathbf{k}^{\prime}}}_{n_{\bar{\alpha}^{\prime}}}\rangle$
(15) $\displaystyle=$
$\displaystyle\langle\phi^{\mathbf{k+q}}_{n_{\alpha}}|\phi^{\mathbf{k^{\prime}}}_{n_{\alpha^{\prime}}}\rangle\langle\psi^{\bar{\mathbf{k}}}_{n_{\bar{\alpha}}}|\psi^{\bar{\mathbf{k}}^{\prime}+\mathbf{q}}_{n_{\bar{\alpha}^{\prime}}}\rangle.$

Since the Wannier function $\phi^{\mathbf{k}}_{n_{\alpha}}$ is complete and
orthonomalized, the Form factor can be generally given by the form of delta
functions obeying the momentum and orbital conservation:

$\displaystyle F_{nn^{\prime}\bar{n}\bar{n}^{\prime}}(\mathbf{q})$
$\displaystyle=$
$\displaystyle\delta_{\mathbf{k}+\mathbf{q},\mathbf{k^{\prime}}}\cdot\delta_{\bar{\mathbf{k}}^{\prime}+\mathbf{q},\bar{\mathbf{k}}}\cdot\delta_{n_{\alpha},n_{\alpha^{\prime}}}\cdot\delta_{n_{\bar{\alpha}},n_{\bar{\alpha}^{\prime}}}.$
(16)

Here, one should note that the orbital index $n_{\alpha}$ used in computing
the magnetic susceptibility runs over the correlated orbitals with finite
magnetic moments, while the generated Wannier functions in the primitive unit
cell contain both correlated and non-correlated orbitals. For example,
Co1/3NbS2 structure contains two Co ions per primitive unit cell, with the
magnetic moments primarily residing on $d$ orbitals of Co ions, while the band
structure calculation can generate both the Co $d$ and the Nb $d_{z^{2}}$
Wannier orbitals as they are strongly hybridized. Once the Form factor
considers only Co ions, one can define a new momentum vector
$\tilde{\mathbf{k}}$ to describe the Wannier function $\tilde{\phi}$ for Co
ions only, which will form effectively a single triangular lattice with the
new periodicity:

$\tilde{\phi}^{\tilde{\mathbf{k}}}_{\alpha}(\mathbf{r})=\frac{1}{\sqrt{NN_{\tau}}}\sum_{\mathbf{R},\vec{\tau}_{\alpha}}e^{-i\tilde{\mathbf{k}}\cdot(\mathbf{R}+\vec{\tau}_{\alpha})}\tilde{\phi}^{\tau_{\alpha}}_{\alpha}(\mathbf{r}-\mathbf{R}-\vec{\tau}_{\alpha}).$
(17)

Here, the $\tilde{\phi}$ is obtained from the Fourier transform using
explicitly the position $\tau_{\alpha}$ of the correlated orbitals.

To evaluate the matrix element in Eq. 15, we expand it using the new Wannier
orbital $\tilde{\phi}$:

$\displaystyle\langle\phi^{\mathbf{k}}_{n_{\alpha}}|e^{i\mathbf{q}\cdot\hat{\mathbf{R}}^{\prime}}|\phi^{\mathbf{k^{\prime}}}_{n_{\alpha^{\prime}}}\rangle$
$\displaystyle=$
$\displaystyle\sum_{\mathbf{\tilde{k}},\tilde{\mathbf{k}}^{\prime}}\sum_{\tilde{\alpha},\tilde{\alpha}^{\prime}}\langle\phi^{\mathbf{k}}_{n_{\alpha}}|\tilde{\phi}^{\tilde{\mathbf{k}}}_{\tilde{\alpha}}\rangle\langle\tilde{\phi}^{\tilde{\mathbf{k}}}_{\tilde{\alpha}}|e^{i\mathbf{q}\cdot\mathbf{\hat{R}^{\prime}}}|\tilde{\phi}^{\tilde{\mathbf{k}}^{\prime}}_{\tilde{\alpha}^{\prime}}\rangle$
(18)
$\displaystyle\cdot\langle\tilde{\phi}^{\tilde{\mathbf{k^{\prime}}}}_{\tilde{\alpha}^{\prime}}|\phi^{\mathbf{k^{\prime}}}_{n_{\alpha^{\prime}}}\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{N_{\tau}}\sum_{\mathbf{\tilde{k}},\mathbf{\tilde{k}}^{\prime}}\sum_{\vec{\tau}_{\alpha},\vec{\tau}_{\alpha^{\prime}}}e^{i(\tilde{\mathbf{k}}^{\prime}\cdot\vec{\tau}_{\alpha^{\prime}}-\tilde{\mathbf{k}}\cdot\vec{\tau}_{\alpha})}\cdot\delta_{\alpha\alpha^{\prime}}$
$\displaystyle\cdot\delta_{\\{\mathbf{k}\\},\tilde{\mathbf{k}}}\cdot\delta_{\\{\mathbf{k}^{\prime}\\},\tilde{\mathbf{k}}^{\prime}}\cdot\delta_{\tilde{\mathbf{k}}+\\{\mathbf{q}\\},\tilde{\mathbf{k}}^{\prime}},$

where the momentum $\tilde{\mathbf{k}}$ is defined in an extended BZ obtained
for a single Co ion in triangular lattice, therefore the overlap between
$\phi^{\mathbf{k}}_{n_{\alpha}}$ and
$\tilde{\phi}^{\tilde{\mathbf{k}}}_{\tilde{\alpha}}$ can produce the
modulation factor obtained from the Co ion positions due to the effective band
unfolding effect [25].

To obtain Eq. 18, we evaluate the overlap between these orbitals as follows:

$\displaystyle\langle\phi_{n_{\alpha}}^{\mathbf{k}}|\tilde{\phi}_{\tilde{\alpha}}^{\tilde{\mathbf{k}}}\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{N\sqrt{N_{\tau}}}\sum_{\mathbf{R},\vec{\tau}_{\tilde{\alpha}}}e^{i(\mathbf{k}\cdot\mathbf{R}-\tilde{\mathbf{k}}\cdot(\mathbf{R}+\vec{\tau}_{\tilde{\alpha}}))}\langle\phi_{\alpha}^{\tau_{\alpha}}|\tilde{\phi}_{\tilde{\alpha}}^{\tau_{\tilde{\alpha}}}\rangle$
(19) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N_{\tau}}}\sum_{\vec{\tau}_{\tilde{\alpha}}}e^{-i\tilde{\mathbf{k}}\vec{\tau}_{\tilde{\alpha}}}\cdot\delta_{\\{\mathbf{k}\\},\tilde{\mathbf{k}}}\cdot\delta_{\alpha,\tilde{\alpha}}$

Similarly,

$\displaystyle\langle\tilde{\phi}_{\tilde{\alpha}}^{\tilde{\mathbf{k}}}|e^{i\mathbf{q}\cdot\mathbf{\hat{R}^{\prime}}}|\tilde{\phi}_{\tilde{\alpha}^{\prime}}^{\tilde{\mathbf{k}}^{\prime}}\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{NN_{\tau}}\sum_{\mathbf{R^{\prime}},\vec{\tau}_{\tilde{\alpha}}}e^{i(\tilde{\mathbf{k}}-\tilde{\mathbf{k}}^{\prime}+\mathbf{q})\cdot(\mathbf{R^{\prime}}+\vec{\tau}_{\tilde{\alpha}})}\langle\tilde{\phi}_{\tilde{\alpha}}^{\tau_{\tilde{\alpha}}}|\tilde{\phi}_{\tilde{\alpha}^{\prime}}^{\tau_{\tilde{\alpha}^{\prime}}}\rangle$
(20) $\displaystyle=$
$\displaystyle\delta_{\tilde{\mathbf{k}}+\\{\mathbf{q}\\},\tilde{\mathbf{k}}^{\prime}}\cdot\delta_{\tilde{\alpha},\tilde{\alpha}^{\prime}}$

While the magnetic susceptibility in the main text is computed using the Form
factor in Eq. 9, the Nb $d_{z^{2}}$ $\chi^{0}$ profile can change notably once
the matrix element effect of the Form factor is simplified using the delta
function form in Eq. 16 for the susceptibility calculation. In Co1/3NbS2, the
notable change occurs at the $M$ point where the peak height is slightly
reduced and becomes almost degenerate with the peak at the $K$ point (see Fig.
9). This shows that the matrix element effect in Eq. 9 can favor the nesting
along the $M$ direction as the Nb $d_{z^{2}}$ orbital modulates along the $M$
direction without the hybridization with Co $d$ orbitals. However, the
localized Co $d$ orbitals have no changes due to the matrix element effect as
their profile is momentum independent and not affected by the hybridization
effect. In both Fe1/3NbS2 and Ni1/3NbS2 cases, the Nb $d_{z^{2}}$
susceptibility peak favors the $K$ point in the BZ when the form factor is
simplified using Eq. 16 (see Fig. 10 and Fig. 11).

Figure 9: The magnetic susceptibility for Co1/3NbS2 at the doping
$\delta$=2.0, the polarizability $\chi^{0}$ (left panel) and the RPA $\chi$
(right panel), computed using the form factor in Eq. 16.

Figure 10: The magnetic susceptibility for Fe1/3NbS2 at the doping
$\delta$=2.0, the polarizability $\chi^{0}$ (left panel) and the RPA $\chi$
(right panel), computed using the form factor in Eq. 16.

Figure 11: The magnetic susceptibility for Ni1/3NbS2 at the doping
$\delta$=2.0, the polarizability $\chi^{0}$ (left panel) and the RPA $\chi$
(right panel), computed using the form factor in Eq. 16.