From fd98919bda4ba900559d67ed220f9719d4982969 Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Mon, 29 Aug 2022 14:57:14 +0300
Subject: [PATCH] DOC
---
README.md | 4 +-
doc/discretization.md | 711 ++++++++++++++++++++++++++++++++++++++++++
2 files changed, 713 insertions(+), 2 deletions(-)
create mode 100644 doc/discretization.md
diff --git a/README.md b/README.md
index 146c623..5b76791 100644
--- a/README.md
+++ b/README.md
@@ -6,8 +6,8 @@ Uses a discretization with local gauge invariance.
## Requirements
-- C++ compiler
-- C++ libraries: Eigen, CppAD
+- C++17 compatible compiler
+- C++ libraries: [Eigen](https://eigen.tuxfamily.org/), [CppAD](https://coin-or.github.io/CppAD/)
- Python 3
## Compiling
diff --git a/doc/discretization.md b/doc/discretization.md
new file mode 100644
index 0000000..1a8001d
--- /dev/null
+++ b/doc/discretization.md
@@ -0,0 +1,711 @@
+# Gauge-invariant discretization
+
+We do the discretization of the Usadel equation as follows. We first
+discretize the action $S$, and derive the discrete saddle-point
+equation from the discretized action. This simplifies things somewhat,
+because proper handling of the gauge fields is simpler to do
+in the action formulation.
+
+The action is
+```math
+ S[Q]
+ &=
+ \mathrm{Tr}[D(\hat{\nabla}Q)^2 + F_{ij}Q\hat{\nabla}_iQ\hat{\nabla}_jQ + \Omega Q]
+ \,,
+```
+where $\hat{\nabla}_j Q = \partial_j Q - i[\mathcal{A}_j, Q]$. We'd
+like the discretized functional to be also gauge invariant.
+
+
+
+We subdivide the space to cells, centered on a rectangular lattice
+$\vec{r}_j=(h j_x, h j_y, h j_z)$ with $j_{x,y,z}\in\mathbb{Z}$
+and $h$ are the lattice spacings. See Fig.~\ref{fig:discr}
+
+We choose $Q(j) = Q(\vec{r}_j)$ to be values of $Q$ at the lattice
+sites. We define the Wilson link matrices
+$U_{ij}
+ =U_{ji}^{-1}
+ =U(\vec{r}_{i},\vec{r}_{j})
+ =
+ \mathrm{Pexp}[i\int_{L(\vec{r}_i,\vec{r}_j)}\dd{\vec{r}'}\cdot\vec{\mathcal{A}}(\vec{r}')]$
+%
+where $L(\vec{r}_i,\vec{r}_j)$ is straight line from $\vec{r}_j$ to
+$\vec{r}_i$ and $\mathrm{Pexp}$ is the path-ordered integral.
+The neighbor cells of $j$ are $i\in{}\mathrm{neigh}(j)$
+ie. $i=(j_x\pm1,j_y,j_z)$, $i=(j_x,j_y\pm1,j_z)$,
+$i=(j_x,j_y,j_z\pm1)$. Since the gauge field $\vec{\mathcal{A}}$ is
+fixed, the link matrices can be computed cheaply, and there's no need
+to expand them.
+
+We want a discretized theory that is invariant under the discrete gauge transformation
+%
+\begin{align}
+ Q(j)
+ &\mapsto
+ u(j)
+ Q
+ u(j)^{-1}
+ \,,
+ &
+ U_{ij}
+ \mapsto
+ u(i)U_{ij}u(j)^{-1}
+ \,,
+\end{align}
+%
+and in the continuum limit reduces to the original theory. Such
+discretization has the advantage of being insensitive to the gauge
+choice, and having exact discrete conservation laws associated with
+the symmetry.
+
+We then discretize
+%
+\begin{align}
+ S_2
+ &=
+ \Tr D(\hat{\nabla} Q)^2
+ \mapsto
+ D
+ \frac{2d}{|\mathrm{neigh}|} \frac{h^d}{h^2} \sum_{\mathrm{neigh}(i,j)}\Tr[Q(i) - U_{ij}Q(j)U_{ji}]^2
+ \\
+ &\simeq
+ \frac{h^{d}d}{|\mathrm{neigh}|}
+ \sum_{i}\sum_{j\in{}\mathrm{neigh}(i)}\Tr\Bigl(\frac{Q(i) - Q(j)}{h} - [i\mathcal{A}(\frac{\vec{r}_i+\vec{r}_j}{2}), Q(j)]\Bigr)^2
+ \,,
+\end{align}
+%
+where $d$ is the space dimension. The first line is gauge
+invariant. The second line verifies the reduction to the continuum
+limit.
+
+The link matrices have the usual continuum limit expansion
+%
+\begin{align}
+ U_{ji}
+ &=
+ %\sum_{n=0}^\infty \frac{i^n}{n!} |\vec{r}_i-\vec{r}_j|^n \int_{-1/2}^{1/2} \dd{s_1}\!\!\cdots{}\dd{s_n} P[\mathcal{A}(s_1)\cdots\mathcal{A}(s_n)]
+ %=
+ \sum_{n=0}^\infty i^n |\vec{r}_i-\vec{r}_j|^n \int_{-1/2}^{1/2}\dd{s_1}\int_{-1/2}^{s_1}\dd{s_2}\ldots\int_{-1/2}^{s_{n-1}}\dd{s_n} \mathcal{A}(s_1)\cdots\mathcal{A}(s_n)
+ \label{eq:link-expansion}
+ \\
+ &\approx
+ 1 + i h \mathcal{A}(s=0)
+ +
+ \frac{(ih)^2}{2}\mathcal{A}(s=0)^2
+ +
+ \mathcal{O}(h^3)
+ \,,
+\end{align}
+%
+where $\mathcal{A}(s)=\mathcal{A}[(\frac{1}{2} - s)\vec{r}_i + (\frac{1}{2} + s)\vec{r}_j]$. The second line indicates the midpoint rule,
+found using $\mathcal{A}(s)\simeq\mathcal{A}(\frac{1}{2}) + (s-\frac{1}{2})\mathcal{A}'(\frac{1}{2})+\ldots$
+and observing that and $\frac{d^n}{ds^n}\mathcal{A}(s) = \mathcal{O}(h^n)$.
+
+As well known in lattice QCD, the field strength can be expressed in
+terms of the link matrices $U$ via a plaquette loop. Consider then
+the plaquette (see Fig.~\ref{fig:discr}), with edges $j\leftarrow{}i$ and $k\leftarrow{}l$ in
+direction $\mu$ and $l\leftarrow{}i$, $k\leftarrow{}j$ in direction $\nu$, and expand around
+$\vec{r}_0=\frac{\vec{r}_i+\vec{r}_j+\vec{r}_k+\vec{r}_l}{4}$:
+%
+\begin{align}
+ P_{lkji}
+ &=
+ U_{il}U_{lk}U_{kj}U_{ji}
+ \,,
+ \\
+ U_{ji}
+ &\simeq
+ 1 + i h \mathcal{A}_\mu(\vec{r}_0) + \frac{1}{2}(i h \mathcal{A}_\mu(\vec{r}_0))^2 + i h (-h)\partial_\nu\mathcal{A}_{\mu}(\vec{r}_0)
+ +
+ \ldots
+ \,,
+ \\
+ U_{kj}
+ &\simeq
+ 1 + i h \mathcal{A}_\nu(\vec{r}_0) + \frac{1}{2}(i h \mathcal{A}_\nu(\vec{r}_0))^2 + i h (+h)\partial_\mu\mathcal{A}_{\nu}(\vec{r}_0)
+ +
+ \ldots
+ \,,
+ \\
+ U_{lk}
+ &\simeq
+ 1 - i h \mathcal{A}_\mu(\vec{r}_0) + \frac{1}{2}(-i h \mathcal{A}_\mu(\vec{r}_0))^2 - i h (+h)\partial_\nu\mathcal{A}_{\mu}(\vec{r}_0)
+ +
+ \ldots
+ \,,
+ \\
+ U_{il}
+ &\simeq
+ 1 - i h \mathcal{A}_\nu(\vec{r}_0) + \frac{1}{2}(-i h \mathcal{A}_\nu(\vec{r}_0))^2 - i h (-h)\partial_\mu\mathcal{A}_{\mu}(\vec{r}_0)
+ +
+ \ldots
+ \,,
+ \\
+ \Rightarrow
+ P_{lkji}
+ &\simeq
+ 1
+ - i h^2
+ \bigl(
+ \partial_\mu \mathcal{A}_\nu - \partial_\nu \mathcal{A}_\mu - i[\mathcal{A}_\mu,\mathcal{A}_\nu]
+ \bigr)
+ +
+ \mathcal{O}(h^3)
+ =
+ 1
+ + i h^2 F_{\mu \nu}(\vec{r}_0)
+ +
+ \mathcal{O}(h^3)
+ \\
+ \Rightarrow
+ F_{\mu \nu}(\vec{r}_0)
+ &\simeq
+ -\frac{1}{i h^2}(1 - P_{lkji})
+ \simeq
+ +
+ \frac{1}{i h^2}(1 - P_{jkli})
+\end{align}
+%
+Let us then construct computer algebra for doing the continuum limit expansions to order $h^4$,
+and use it to verify this. First, define the algebra of noncommutative indexed symbols $a$, $Q$,
+with coefficients being polynomials of $h$ with $h^5=0$. The indices on $a$ and $Q$ indicate plain derivatives.
+%
+\begin{sageblock}
+from indexed_letter_monoid import IndexedLetterMonoid, relabel, distributive_op
+from reduce_term import reduce_term, coerce_algebra, divide_monomial, divide_monomial_cyclic
+
+var_props = {'Q': 'involution', 'q': 'involution', 'a': 'symmetric-tail', 'F': '2-idx'}
+GG = IndexedLetterMonoid(letter_init=tuple(var_props.items()))
+QQi = I.parent()
+PP0 = PolynomialRing(QQi, 'h,s')
+PP = PP0.quotient(PP0.ideal(PP0.gens()[0]**5), names='h,s')
+PP.inject_variables()
+FF = GG.algebra(PP)
+FF0 = GG.algebra(PP0)
+def a(idxs): return FF(GG([("a", tuple(idxs))]))
+def Q(idxs=()): return FF(GG([("Q", tuple(idxs))]))
+def q(idxs=()): return FF(GG([("q", tuple(idxs))]))
+def F(idxs=()): return FF(GG([("F", tuple(idxs))]))
+\end{sageblock}
+%
+Define taylor expansions of $a$ and $Q$ matrices around center of the
+plaquette, and compute $U$ from Eq.~\eqref{eq:link-expansion}. We want
+the conditions $U_{ij}U_{ji}=1$ and $Q^2=1$ to hold to order $h^4$.
+%
+\begin{sageblock}
+@mem_cache
+def symtaylor(fun, mu, dx, dy, order=4):
+ res = fun(mu)
+ for n in range(1, order + 1):
+ for ni in range(n + 1):
+ nj = n - ni
+ res = res + PP(QQi(binomial(n, ni)) / factorial(n)
+ * dx**ni * dy**nj) * fun(mu + "i"*ni + "j"*nj)
+ return res
+
+def Qexpansion(dx, dy, order=4):
+ return symtaylor(Q, "", dx, dy, order=order)
+
+def aexpansion(mu, dx, dy, order=4):
+ return symtaylor(a, mu, dx, dy, order=order)
+
+@mem_cache
+def Uexpansion(mu, dx, dy, sgn=1, order=4):
+ # path-ordered integral of midpoint Taylor expansion of a(x,y)
+ if dx == 0:
+ b = aexpansion(mu, sgn * s * h, dy, order=order)
+ elif dy == 0:
+ b = aexpansion(mu, dx, sgn * s * h, order=order)
+ else:
+ raise ValueError()
+
+ res = 1
+ for n in range(1, order + 1):
+ res = res + PSintegral(b, n) * (I*h*sgn)**n
+ return res
+
+def PSintegral(b, n, expr=None):
+ r"""
+ Pexp(\int_{-1/2}^{1/2} ds b(s))
+ """
+ ring = b.parent().base_ring()
+ if expr is None:
+ expr = b
+ ss = SR(s.lift())
+ if n == 1:
+ ig = expr.map_coefficients(lambda z: ring(integrate(SR(z.lift()), ss, -1/2, 1/2)))
+ return ig
+ else:
+ ig = expr.map_coefficients(lambda z: ring(integrate(SR(z.lift()), ss, -1/2, ss)))
+ return PSintegral(b, n-1, b * ig)
+
+Q1val = Qexpansion(-h/2, -h/2)
+Q2val = Qexpansion(h/2, -h/2)
+Q3val = Qexpansion(h/2, h/2)
+Q4val = Qexpansion(-h/2, h/2)
+
+U21 = Uexpansion("i", 0, -h/2, sgn=1)
+U32 = Uexpansion("j", h/2, 0, sgn=1)
+U43 = Uexpansion("i", 0, h/2, sgn=-1)
+U14 = Uexpansion("j", -h/2, 0, sgn=-1)
+
+U12 = Uexpansion("i", 0, -h/2, sgn=-1)
+U23 = Uexpansion("j", h/2, 0, sgn=-1)
+U34 = Uexpansion("i", 0, h/2, sgn=1)
+U41 = Uexpansion("j", -h/2, 0, sgn=1)
+
+Fijval = a("ji") - a("ij") - I*(a("i")*a("j") - a("j")*a("i"))
+Qival = Q("i") - I*(a("i") * Q() - Q() * a("i"))
+Qjval = Q("j") - I*(a("j") * Q() - Q() * a("j"))
+
+def hcoef(expr, order):
+ if order == 0:
+ return expr.map_coefficients(lambda z: z.lift().constant_coefficient())
+ elif order > 0:
+ return expr.map_coefficients(lambda z: z.lift().monomial_coefficient(h.lift()**order))
+ else:
+ raise ValueError(order)
+
+def htrunc(expr, order):
+ return sum(h**n * hcoef(expr, n) for n in range(order + 1))
+
+@distributive_op
+def tr_permute(m, c, algebra):
+ return min((algebra.term(GG(m.value[k:] + m.value[:k]), c) for k in range(len(m))),
+ default=algebra.term(m, c))
+
+def partitions(total, left=(), right=()):
+ total = tuple(total)
+ if not total:
+ yield left, right
+ else:
+ t2 = total[1:]
+ s = total[0]
+ l2 = left + (s,)
+ r2 = right + (s,)
+ yield from partitions(t2, l2, right)
+ yield from partitions(t2, left, r2)
+
+@mem_cache
+def Qideal():
+ # commutation rules from derivatives of Q^2 = 1
+ return [
+ sum(Q(p1) * Q(p2) for p1, p2 in partitions(ptot))
+ for ptot in ["i", "j", "ii", "ij", "jj",
+ "iii", "iij", "ijj", "jjj",
+ "iiii", "iiij", "iijj", "ijjj", "jjjj"]
+ ]
+
+@mem_cache
+def simplify(expr, trace=False):
+ ideal = Qideal()
+ s = tr_permute if trace else None
+ divide = divide_monomial_cyclic if trace else divide_monomial
+ return reduce_term(expr, ideal, algebra=FF0, simplify=s, divide=divide)
+
+@distributive_op
+def gderiv(m, c, algebra, i):
+ assert len(i) == 1
+ new = []
+ for letter, idx in m.value:
+ new.append((letter, tuple(idx) + tuple(i)))
+ expr0 = FF.term(m, c)
+ expr1 = FF.term(GG(new), c)
+ return expr1 - I*(a(i)*expr0 - expr0*a(i))
+
+@mem_cache
+def to_covariant(expr, trace=False, **kw):
+ comm = lambda a, b: a*b - b*a
+ ideal = [Q() - q()]
+ ideal += [gderiv(Q(), i) - q(i) for i in "ij"]
+ ideal += [gderiv(gderiv(Q(), i), j) - q((i,j)) for i in "ij" for j in "ij"]
+ ideal += [gderiv(gderiv(gderiv(Q(), i), j), k) - q((i,j,k)) for i in "ij" for j in "ij" for k in "ij"]
+ ideal += [Fijval - F(("i", "j"))]
+ ideal += [F(("j", "i")) + F(("i", "j"))]
+ ideal += [gderiv(Fijval, k) - F(("i","j",k)) for k in "ij"]
+ ideal += [F(("j","i",k)) + F(("i","j",k)) for k in "ij"]
+ ideal += Qideal()
+ s = tr_permute if trace else None
+ divide = divide_monomial_cyclic if trace else divide_monomial
+ return reduce_term(expr, ideal, max_key=_max_key_cov, algebra=FF0,
+ simplify=s, divide=divide, **kw)
+
+def _max_key_cov(m):
+ cov_badness = sum(letter in ("a", "Q") for letter, idx in m.value)
+ return (cov_badness, m)
+\end{sageblock}
+%
+The plaquette loop:
+%
+\begin{sageblock}
+P4321 = U14 * U43 * U32 * U21
+\end{sageblock}
+%
+Then, to order $\mathcal{O}(h^3)$,
+%
+\begin{sageexample}
+sage: I * htrunc(1 - P4321, 2)
+\end{sageexample}
+%
+When discretizing the $Q\hat{\nabla}_iQ\hat{\nabla}_jQ$ part, we can
+consider a discrete derivative
+%
+\begin{align}
+ D_{ij}
+ \coloneqq
+ Q(i) U_{ij} - U_{ij} Q(j)
+ \,,
+ \qquad
+ \Rightarrow
+ Q(i) D_{ij} = - D_{ij} Q(j)
+ \,,
+\end{align}
+%
+which satisfies a discrete gauge-invariant version of the $Q (\hat{\nabla}Q) =
+-(\hat{\nabla}Q) Q$ relation. Then, in the plaquette of Fig.~\ref{fig:discr} we
+have
+%
+\begin{align}
+ \tr
+ F_{\mu\nu}Q \hat{\nabla}_\mu Q \hat{\nabla}_\nu Q
+ \rvert_{\mu=ji, \nu=li}
+ &\simeq
+ \frac{i}{h^4}
+ \tr
+ (1 - P_{lkji}) U_{ij} D_{ji} Q(i) D_{il} U_{li}
+ =
+ \frac{-i}{h^4}
+ \tr
+ (U_{lk} U_{kj} - U_{li} U_{ij}) D_{ji} Q(i) D_{il}
+\end{align}
+%
+So the structure is $\frac{-i}{h^4}$ $\times$ (gauge factors for
+counterclockwise loop $-$ return back clockwise) $\times$ $DQD$ for
+three sites in counterclockwise order. Since only $Q(i)$, $Q(j)$,
+$Q(l)$ appear here, one can also think of this as a plaquette corner
+with $Q(i)$ being the corner.
+
+Indeed:
+%
+\begin{sageblock}
+Fijval = a("ji") - a("ij") - I*(a("i")*a("j") - a("j")*a("i"))
+Fjival = -Fijval
+Qival = Q("i") - I*(a("i") * Q() - Q() * a("i"))
+Qjval = Q("j") - I*(a("j") * Q() - Q() * a("j"))
+
+D41 = Q4val * U41 - U41 * Q1val
+D34 = Q3val * U34 - U34 * Q4val
+D21 = Q2val * U21 - U21 * Q1val
+D32 = Q3val * U32 - U32 * Q2val
+
+D14 = Q1val * U14 - U14 * Q4val
+D43 = Q4val * U43 - U43 * Q3val
+D12 = Q1val * U12 - U12 * Q2val
+D23 = Q2val * U23 - U23 * Q3val
+
+f_1 = -I*((U43*U32 - U41*U12) * D21*Q1val*D14)
+f_2 = -I*((U14*U43 - U12*U23) * D32*Q2val*D21)
+f_3 = -I*((U21*U14 - U23*U34) * D43*Q3val*D32)
+f_4 = -I*((U32*U21 - U34*U41) * D14*Q4val*D43)
+
+fb_xy = Fijval*Q()*Qival*Qjval * h**4
+fb_yx = Fjival*Q()*Qjval*Qival * h**4
+\end{sageblock}
+\begin{sageexample}
+sage: simplify(f_1 - fb_xy, trace=True)
+sage: simplify(f_2 - fb_yx, trace=True)
+sage: simplify(f_3 - fb_xy, trace=True)
+sage: simplify(f_4 - fb_yx, trace=True)
+\end{sageexample}
+
+To do the sum over the $\mu$, $\nu$ indices, one can average over
+corners of the plaquettes as follows:
+%
+\begin{align}
+ S_H
+ &=
+ \alpha
+ \Tr F_{\mu\nu}Q\hat{\nabla}_\mu Q\hat{\nabla}_\nu Q
+ \\
+ &\mapsto
+ \alpha
+ \frac{h^d d}{4ih^4}
+ \sum_{\mathrm{plaqc}(ijkl; i)}
+ \Tr
+ \Bigl(
+ (U_{lk} U_{kj} - U_{li} U_{ij}) D_{ji} Q(i) D_{il}
+ \Bigr)
+ \\
+ &=
+ \frac{\alpha h^{d-4} d}{4i}
+ \sum_{\mathrm{plaqc}(i; jkl)}
+ \Tr
+ \Bigl(
+ (P_{lkji} - 1)(Q(i) - U_{ij} Q(j) U_{ji} Q(i) U_{il} Q(l) U_{li})
+ \Bigr)
+ \,,
+\end{align}
+%
+where $\mathrm{plaqc}(i;jkl)$ are the plaquette corners. The last line
+follows with direct calculation, and considering the cyclic
+permutations.
+
+Indeed:
+%
+\begin{sageblock}
+@mem_cache
+def get_fc():
+ P4321 = U14 * U43 * U32 * U21
+ P3214 = U43 * U32 * U21 * U14
+ P2143 = U32 * U21 * U14 * U43
+ P1432 = U21 * U14 * U43 * U32
+
+ fc_1 = (P4321-1) * (Q1val - U12*Q2val*U21*Q1val*U14*Q4val*U41)
+ fc_2 = (P1432-1) * (Q2val - U23*Q3val*U32*Q2val*U21*Q1val*U12)
+ fc_3 = (P2143-1) * (Q3val - U34*Q4val*U43*Q3val*U32*Q2val*U23)
+ fc_4 = (P3214-1) * (Q4val - U41*Q1val*U14*Q4val*U43*Q3val*U34)
+
+ return (-I/2) * (fc_1 + fc_2 + fc_3 + fc_4)
+
+fc_tot = get_fc()
+\end{sageblock}
+\begin{sageexample}
+sage: simplify(fc_tot - (fb_xy + fb_yx), trace=True)
+\end{sageexample}
+
+The last term is simple:
+%
+\begin{align}
+ \Tr \Omega Q
+ \mapsto
+ \sum_{i} h^3 \Tr \Omega(j) Q(j)
+ \,.
+\end{align}
+%
+It is gauge invariant if $\Omega$ transforms as
+$\Omega(j)\mapsto{}u(j)\Omega(j)u(j)^{-1}$.
+
+We can also calculate the matrix current $j\mapsto{}i$ between neighboring cells:
+%
+\begin{align}
+ J_{ij}^a
+ &=
+ \frac{\partial}{\partial\xi}
+ S[
+ U_{ij}
+ \mapsto
+ e^{i\xi h T^a} U_{ij}
+ ]\rvert_{\xi=0}
+ =
+ \frac{\partial}{\partial\xi}
+ S[
+ U_{ij}
+ \mapsto
+ U_{ij}
+ +
+ \frac{ih}{2}T^a U_{ij}\xi
+ ]\rvert_{\xi=0}
+\end{align}
+%
+where $T^a$ is an appropriate matrix generator, and the transformation
+is inserted only to one of the links. Indeed, in the continuum limit
+%
+\begin{align}
+ J_{ij}^a
+ &\simeq
+ \frac{\partial}{\partial\xi}
+ S[
+ U_{ij}
+ \mapsto
+ 1
+ +
+ i
+ h
+ \mathcal{A}
+ +
+ ih\xi T^a
+ ]\rvert_{\xi=0}
+ \simeq
+ \frac{\partial}{\partial\xi}
+ S[
+ \mathcal{A}
+ \mapsto
+ \mathcal{A}
+ +
+ \xi T^a
+ ]\rvert_{\xi=0}
+ \,.
+\end{align}
+%
+The local gauge invariance now implies that there is an exact discrete
+continuity equation
+%
+\begin{align}
+ \sum_{j\in{}N_i} J_{ij}^a
+ =
+ \frac{\partial}{\partial\xi}
+ S[Q(i)\mapsto{}e^{i\xi T^a}Q(i)e^{-i\xi T^a}]\rvert_{\xi=0}
+ =
+ \frac{\partial}{\partial\xi}
+ S[\Omega(i)\mapsto{}e^{-i\xi T^a}\Omega(i)e^{i\xi T^a}]\rvert_{\xi=0}
+ =
+ R_{i}^a
+ \,.
+\end{align}
+%
+This is one main advantage of the procedure, in addition to its
+invariance vs. gauge fixing.
+
+However, the definition of the matrix current is \emph{asymmetric}, so
+that $J_{ij}\ne{}-J_{ji}$ as they differ by a quantity of order
+$\mathcal{O}(h^2)$. The asymmetry arises when $T^a$ does not commute
+with $U_{ij}$, so it is an issue only for the SU(2) component. The
+problem is to some degree just in the interpretation of what $J_{ij}$
+means: it is the incoming current measured at site $i$, whereas
+$J_{ji}$ is measured as site $j$. Since the parallel transport of
+spin between these two locations can imply rotation due to the gauge
+field, there's no reason why we should have $J_{ij} = -J_{ji}$, except
+in the continuum limit where the two points become close to each
+other.
+
+\subsection{Saddle point}
+
+We have the variations
+%
+\begin{align}
+ \delta S_2
+ &=
+ \frac{2Dd}{|\mathrm{neigh}|} \frac{h^d}{h^2} \sum_{\mathrm{neigh}(i,j)}
+ \Tr\Bigl(
+ 2 \delta Q(i) [Q(i) - U_{ij}Q(j)U_{ji}]
+ +
+ 2 \delta Q(j) [Q(j) - U_{ji} Q(i) U_{ij}]
+ \Bigr)
+ \\
+ &=
+ \frac{4Dh^{d-2}d}{|\mathrm{neigh}|}
+ \sum_{i}
+ \sum_{j\in\mathrm{neigh}(i,j)}
+ \Tr\Bigl( \delta Q(i) [Q(i) - U_{ij}Q(j)U_{ji}]\Bigr)
+ \,,
+\end{align}
+%
+and
+%
+\begin{align}
+ \delta S_H
+ &\simeq
+ \alpha
+ \frac{h^{d-4}d}{4i}
+ \sum_{\mathrm{plaqc}(i; jkl)}
+ \Tr
+ \Bigl(
+ (P_{lkji} - 1)(\delta Q(i)
+ - U_{ij} \delta Q(j) U_{ji} Q(i) U_{il} Q(l) U_{li}
+ - U_{ij} Q(j) U_{ji} \delta Q(i) U_{il} Q(l) U_{li}
+ \\\notag&\qquad
+ - U_{ij} Q(j) U_{ji} Q(i) U_{il} \delta Q(l) U_{li}
+ )
+ \Bigr)
+ \,,
+ \\
+ &=
+ \alpha
+ \frac{h^{d-4}d}{4i}
+ \sum_{i}
+ \sum_{jkl\in\mathrm{plaqc}(i; jkl)}
+ \Tr
+ \delta Q(i)
+ \Bigl(
+ (P_{lkji} - 1)
+ - U_{il} Q(l) U_{lk} Q(k) U_{kl} (P_{kjil} - 1) U_{li}
+ - U_{il} Q(l) U_{li}(P_{lkji} - 1) U_{ij} Q(j) U_{ji}
+ \\\notag&\qquad
+ - U_{ij} () U_{jk} Q(k) U_{kj} Q(j) U_{ji}
+ \Bigr)
+ \,.
+\end{align}
+%
+Hence, the lattice gauge-invariant Usadel equation is
+%
+\begin{align}
+ [Q(i), Z_2(i) + Z_H(i)] &= 0
+ \,,
+ \\
+ Z_2(i)
+ &=
+ \frac{4Dh^{d-2}d}{|\mathrm{neigh}|} \sum_{j\in\mathrm{neigh}(i,j)}
+ [Q(i) - U_{ij}Q(j)U_{ji}]
+ \,,
+ \\
+ Z_H(i)
+ &=
+ \frac{\alpha h^{d-4}d}{4}
+ \sum_{jkl\in\mathrm{plaqc}(i; jkl)}
+ \Bigl(
+ F_{lkji}
+ \\\notag&\qquad
+ - U_{il} Q(l) U_{lk} Q(k) U_{kl}U_{li} F_{lkji}
+ \\\notag&\qquad
+ - U_{il} Q(l) U_{li} F_{lkji} U_{ij} Q(j) U_{ji}
+ \\\notag&\qquad
+ - F_{lkji} U_{ij} U_{jk} Q(k) U_{kj} Q(j) U_{ji}
+ \Bigr)
+ \,,
+ \qquad
+ \qquad
+ F_{lkji}
+ \coloneqq
+ i(1 - P_{lkji})
+ \,,
+ \\
+ &\hat{=}
+ \frac{\alpha h^{d-4} d}{4}
+ \sum_{jkl\in\mathrm{plaqc}(i; jkl)}
+ \Bigl(
+ F_{lkji} - U_{il} Q(l) U_{lk} Q(k) U_{kl}U_{li} F_{lkji}
+ \\\notag&\qquad
+ + Q(i) F_{lkji} Q(i) - U_{il} Q(l) U_{li} F_{lkji} U_{ij} Q(j) U_{ji}
+ \\\notag&\qquad
+ + F_{lkji} - F_{lkji} U_{ij} U_{jk} Q(k) U_{kj} Q(j) U_{ji}
+ \Bigr)
+\end{align}
+%
+In the discretization of the $F_{\mu\nu}$ term, we added
+$[Q(i),F_{lkji} + Q(i)F_{lkji}Q(i)]=0$ to the equation, so that the
+leading order in $h$ of the continuum expansion of $Z_H$ is canceled.
+
+We can verify that
+%
+\begin{sageblock}
+Fp = I*(1 - P4321)
+
+ZH = (
+ Fp - U14*Q4val*U43*Q3val*U32*U21*Fp
+ + Q1val*Fp*Q1val - U14*Q4val*U41*Fp*U12*Q2val*U21
+ + Fp - Fp*U12*U23*Q3val*U32*Q2val*U21)
+\end{sageblock}
+
+\subsection{Implementation}
+
+The above is sufficient for numerical implementation: it is not
+necessary to find out the saddle point equation symbolically. It can
+be deduced by parametrizing first so that $Q^2=1$ and then using
+automatic differentiation (AD). Only a routine evaluating the value of
+the action is then needed.
+
+When looking for the saddle point in equilibrium, the problem becomes
+even simpler; $Q$ has then additional restrictions, and we are looking
+for the minimum of the free energy $S$.
+
+
+
+\bibliography{socquestion}
+
+\end{document}
+
+%%% Local Variables:
+%%% mode: latex
+%%% TeX-master: t
+%%% End:
--
GitLab