DOC

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

doc/discretization.md

+# 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.
+
+![Lattice discretization](discretization.svg "Lattice discretization.")
+
+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: