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

There was an error rendering this math block. KaTeX parse error: Expected 'EOF', got '&' at position 10: S[Q] &̲= \mathrm{Tr}…

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) * dxni * 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 - expr0a(i))

@mem_cache def to_covariant(expr, trace=False, **kw): comm = lambda a, b: ab - ba 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*((U43U32 - U41U12) * D21Q1valD14) f_2 = -I*((U14U43 - U12U23) * D32Q2valD21) f_3 = -I*((U21U14 - U23U34) * D43Q3valD32) f_4 = -I*((U32U21 - U34U41) * D14Q4valD43)

fb_xy = FijvalQ()QivalQjval * h**4 fb_yx = FjivalQ()QjvalQival * 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 - U14Q4valU43Q3valU32U21Fp + Q1valFpQ1val - U14Q4valU41FpU12Q2valU21 + Fp - FpU12U23Q3valU32Q2valU21) \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: