From bf6a3a887de0da7f37e7d80eddc37812c3df446b Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Mon, 29 Aug 2022 15:30:19 +0300
Subject: [PATCH] DOC

---
 doc/discretization.md            | 669 ++++++-------------------------
 doc/{ => img}/discretization.svg |   0
 doc/pandoc/convert.sh            |   8 +
 doc/pandoc/fix-links.lua         |  44 ++
 doc/pandoc/gitlab-math.lua       |  18 +
 5 files changed, 184 insertions(+), 555 deletions(-)
 rename doc/{ => img}/discretization.svg (100%)
 create mode 100755 doc/pandoc/convert.sh
 create mode 100644 doc/pandoc/fix-links.lua
 create mode 100644 doc/pandoc/gitlab-math.lua

diff --git a/doc/discretization.md b/doc/discretization.md
index 1a8001d..b93102c 100644
--- a/doc/discretization.md
+++ b/doc/discretization.md
@@ -1,150 +1,118 @@
 # 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
+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]
+  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
+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}
+$`\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 the figure below:
+
+![Lattice discretization](img/discretization.svg "FIG 1. Lattice discretization.")
 
-We choose $Q(j) = Q(\vec{r}_j)$ to be values of $Q$ at the lattice
+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_{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
+  \mathrm{Pexp}[i\int_{L(\vec{r}_i,\vec{r}_j)}d\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}
-  \,,
-  &
+```math
+Q(j) \mapsto u(j) Q u(j)^{-1} \,,
+```
+and
+```math
   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.
 
+## Derivative term
+
 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.
+```math
+S_2
+=
+\mathrm{Tr} D(\hat{\nabla} Q)^2
+\mapsto
+D
+\frac{2d}{|\mathrm{neigh}|} \frac{h^d}{h^2} \sum_{\mathrm{neigh}(i,j)}\mathrm{Tr}[Q(i) - U_{ij}Q(j)U_{ji}]^2
+```
+where $`d`$ is the space dimension.
+
+Expanding in small $h$ we have
+```math
+S_2
+\simeq
+\frac{h^{d}d}{|\mathrm{neigh}|}
+\sum_{i}\sum_{j\in{}\mathrm{neigh}(i)}\mathrm{Tr}\Bigl(\frac{Q(i) - Q(j)}{h} - [i\mathcal{A}(\frac{\vec{r}_i+\vec{r}_j}{2}), Q(j)]\Bigr)^2
+\,,
+```
+which verifies the reduction to the continuum limit.
+
+## Field strength
 
 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
+```math
+U_{ji}
+=
+\sum_{n=0}^\infty i^n |\vec{r}_i-\vec{r}_j|^n \int_{-1/2}^{1/2}ds_1\int_{-1/2}^{s_1}ds_2\ldots\int_{-1/2}^{s_{n-1}}ds_n\, \mathcal{A}(s_1)\cdots\mathcal{A}(s_n)
+```
+```math
+U_{ji}
+\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)$.
+```
+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}
+terms of the link matrices $`U`$ via a plaquette loop.  Consider then
+the plaquette (see Fig. 1), 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}`$:
+```math
   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
-  \,,
-  \\
+```
+```math
   \Rightarrow
   P_{lkji}
-  &\simeq
+  \simeq
   1
   - i h^2
   \bigl(
@@ -157,340 +125,88 @@ $\vec{r}_0=\frac{\vec{r}_i+\vec{r}_j+\vec{r}_k+\vec{r}_l}{4}$:
   + 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))
+```
+which then allows expressing $`F_{ij}`$ in terms of the link matrices.
 
-@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)
+## Hall term
 
-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}
+```math
   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
+```
+which satisfies a discrete gauge-invariant version of the $`Q (\hat{\nabla}Q) =
+-(\hat{\nabla}Q) Q`$ relation.  Then, in the plaquette of Fig. 1 we
 have
-%
-\begin{align}
-  \tr
+```math
+  \mathrm{tr}
   F_{\mu\nu}Q \hat{\nabla}_\mu Q \hat{\nabla}_\nu Q
   \rvert_{\mu=ji, \nu=li}
-  &\simeq
+  \simeq
   \frac{i}{h^4}
-  \tr
+  \mathrm{tr}
   (1 - P_{lkji}) U_{ij} D_{ji} Q(i) D_{il} U_{li}
+```
+```math
   =
   \frac{-i}{h^4}
-  \tr
+  \mathrm{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}
+corners of the plaquettes as follows, which gives the final discretization
+of the Hall term:
+```math
   S_H
-  &=
+  =
   \alpha
-  \Tr F_{\mu\nu}Q\hat{\nabla}_\mu Q\hat{\nabla}_\nu Q
-  \\
-  &\mapsto
+  \mathrm{Tr} F_{\mu\nu}Q\hat{\nabla}_\mu Q\hat{\nabla}_\nu Q
+```
+```math
+  S_H
+  \mapsto
   \alpha
   \frac{h^d d}{4ih^4}
   \sum_{\mathrm{plaqc}(ijkl; i)}
-  \Tr
+  \mathrm{Tr}
   \Bigl(
   (U_{lk} U_{kj} - U_{li} U_{ij}) D_{ji} Q(i) D_{il}
   \Bigr)
-  \\
-  &=
+```
+```math
+  S_H
+  =
   \frac{\alpha h^{d-4} d}{4i}
   \sum_{\mathrm{plaqc}(i; jkl)}
-  \Tr
+  \mathrm{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}$.
+## Matrix current
 
 We can also calculate the matrix current $j\mapsto{}i$ between neighboring cells:
-%
-\begin{align}
+```math
   J_{ij}^a
-  &=
+  =
   \frac{\partial}{\partial\xi}
   S[
     U_{ij}
@@ -506,42 +222,14 @@ We can also calculate the matrix current $j\mapsto{}i$ between neighboring cells
     +
     \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}
-%
+is inserted only to one of the links.
+
 The local gauge invariance now implies that there is an exact discrete
 continuity equation
-%
-\begin{align}
+```math
   \sum_{j\in{}N_i} J_{ij}^a
   =
   \frac{\partial}{\partial\xi}
@@ -552,160 +240,31 @@ continuity equation
   =
   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
+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
+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
+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}
+## 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.
+the action is then needed, in terms of the AD dual variables.
 
 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:
diff --git a/doc/discretization.svg b/doc/img/discretization.svg
similarity index 100%
rename from doc/discretization.svg
rename to doc/img/discretization.svg
diff --git a/doc/pandoc/convert.sh b/doc/pandoc/convert.sh
new file mode 100755
index 0000000..1c4608d
--- /dev/null
+++ b/doc/pandoc/convert.sh
@@ -0,0 +1,8 @@
+#!/bin/bash
+PDIR=$(dirname "$0")
+pandoc --from markdown --to html \
+       --embed-resources --standalone \
+       --lua-filter $PDIR/gitlab-math.lua \
+       --lua-filter $PDIR/fix-links.lua \
+       --katex \
+       -o "${1/.md/.html}" "$1"
diff --git a/doc/pandoc/fix-links.lua b/doc/pandoc/fix-links.lua
new file mode 100644
index 0000000..36ac7dc
--- /dev/null
+++ b/doc/pandoc/fix-links.lua
@@ -0,0 +1,44 @@
+-- see if the file exists
+function file_exists(file)
+  local f = io.open(file, "rb")
+  if f then f:close() end
+  return f ~= nil
+end
+
+-- 1. Replaces .md with .html
+-- 2. Replaces absolute paths with relative ones
+-- 3. Replaces folder links with links to their index / readme
+function fix_link (url) 
+  -- Replace md with html
+  fixed_url = url:gsub("%.md", ".html")
+  -- Replace project-root-relative (i.e. /path/to/thing.md) paths with relative ones
+  fixed_url = fixed_url:gsub("^/", (
+    function(s)
+      basedir = io.popen("realpath --relative-to=\"$(dirname \"" .. s:sub(2) .. "\")\" .", 'r'):read('*a')
+      return basedir:gsub("%s", "") .. "/" .. s:sub(2)
+    end
+  ))
+  -- Allowed readme / index names
+  dir_name = {} -- values are functions that return the path to the new link
+  dir_name["index.md"] = (function(s) return s:gsub("/?$", "/index.html") end)
+  dir_name["INDEX.md"] = (function(s) return s:gsub("/?$", "/INDEX.html") end)
+  dir_name["readme.md"] = (function(s) return s:gsub("/?$", "/readme.html") end)
+  dir_name["README.md"] = (function(s) return s:gsub("/?$", "/README.html") end)
+  -- Is the url pointing to a folder?
+  file_or_folder = io.popen("[ -d '" .. fixed_url .. "' ] && echo 'folder'", 'r'):read('*a')
+  -- If it's a folder, find the readme / index and replace the link with that. 
+  if file_or_folder:gsub("%s", "") == "folder" then
+    for index,mk_path in pairs(dir_name) do
+      if file_exists(fixed_url .. "/" .. index) then
+        fixed_url = mk_path(fixed_url)
+        break
+      end
+    end
+  end
+  return fixed_url
+end
+
+
+function Link(link) link.target = fix_link(link.target); return link end
+function Image(img) img.src = fix_link(img.src); return src end
+-- return {{Meta = Meta}, {Link = Link, Image = Image}}
\ No newline at end of file
diff --git a/doc/pandoc/gitlab-math.lua b/doc/pandoc/gitlab-math.lua
new file mode 100644
index 0000000..6161031
--- /dev/null
+++ b/doc/pandoc/gitlab-math.lua
@@ -0,0 +1,18 @@
+function Math(el)
+    if el.mathtype == "InlineMath" then
+        if el.text:sub(1,1) == '`' and el.text:sub(#el.text) == '`' then
+            local text = el.text:sub(2,#el.text-1)
+            return pandoc.Math(el.mathtype, text)
+        else
+            local cont = pandoc.read(el.text)
+            return el
+            --return { pandoc.Str("$") } .. cont.blocks[1].content .. { pandoc.Str("$") }
+        end
+    end
+end
+
+function CodeBlock(el)
+    if el.classes[1] == "math" then
+        return pandoc.Para({ pandoc.Math("DisplayMath", el.text) })
+    end
+end
-- 
GitLab