How does Gauge Cooling Stabilize Complex Langevin?

TL;DR

This paper investigates how gauge cooling stabilizes the complex Langevin method by analyzing exact solutions and stochastic equations, demonstrating its effectiveness in localizing distributions in an $SU(2)$ model.

Contribution

It provides a detailed analysis of the gauge cooling mechanism, deriving stochastic equations and showing how it prevents excursions in the complex Langevin process.

Findings

Gauge cooling minimizes Frobenius norm of link variables.

It helps form localized distributions in the Polyakov loop model.

Gauge cooling prevents excursions away from the real axis.

Abstract

We study the mechanism of the gauge cooling technique to stabilize the complex Langevin method in the one-dimensional periodic setting. In this case, we find the exact solutions for the gauge transform which minimizes the Frobenius norm of link variables. Thereby, we derive the underlying stochastic differential equations by continuing the numerical method with gauge cooling, and thus provide a number of insights on the effects of gauge cooling. A specific case study is carried out for the Polyakov loop model in $SU(2)$ theory, in which we show that the gauge cooling may help form a localized distribution to guarantee there is no excursion too far away from the real axis.

Tables

Table 1. Table 1. Numerical results for the observable ⟨ O k ⟩ delimited-⟨⟩ subscript 𝑂 𝑘 \langle O_{k}\rangle

			$N=16$				$N=32$
$k$	Exact		$\alpha =0.4$	$\alpha =1.0$	Optimal		$\alpha =0.4$	$\alpha =1.0$	Optimal
$1$	$2.0957$		$2.1001$	$2.1042$	$2.1017$		$2.1394$	$2.0900$	$2.0966$
$-1$	$2.1026$		$2.1070$	$2.1109$	$2.1083$		$2.1422$	$2.0958$	$2.1033$
$2$	$0.3761$		$0.3850$	$0.3894$	$0.3775$		$0.6548$	$0.3583$	$0.3750$
$-2$	$0.4092$		$0.4182$	$0.4224$	$0.4100$		$0.5463$	$0.3874$	$0.4078$
$3$	$-0.5269$		$-0.5270$	$-0.5272$	$-0.5482$		$-9.0467$	$-0.5514$	$-0.5247$
$-3$	$-0.4800$		$-0.4803$	$-0.4802$	$-0.5014$		$-2.2547$	$-0.4938$	$-0.4778$

Table 2. Table 2. Numerical results for the observable ⟨ O k ⟩ delimited-⟨⟩ subscript 𝑂 𝑘 \langle O_{k}\rangle

		$A=1,B=0.2$			$A=1,B=2$
$k$		Exact	Numerical		Exact	Numerical
$1$		$0.8759+0.1300i$	$0.8569+0.1332i$		$1.7449+0.5495i$	$1.5447+0.5269i$
$2$		$-0.6017+0.1304i$	$-0.6304+0.1298i$		$0.2936+1.7642i$	$-0.1957+0.9220i$
$3$		$-0.7562-0.0652i$	$-0.7543-0.0738i$		$-2.1127+1.2080i$	$20.6438-34.4723i$
		$A=5,B=1$			$A=5,B=10$
$k$		Exact	Numerical		Exact	Numerical
$1$		$1.7189+0.0544i$	$1.7187+0.0545i$		$1.9390+0.1188i$	$1.9395+0.1195i$
$2$		$1.0021+0.1669i$	$1.0013+0.1674i$		$1.7388+0.4511i$	$1.7406+0.4543i$
$3$		$0.1531+0.2341i$	$0.1513+0.2345i$		$1.3565+0.9231i$	$1.3601+0.9318i$

Equations

$$\mathbf{x}=(i,j,k,t)a,i,j,k,t\in \mathbb{Z},$$

$$Z=\int [\mathrm{d}U]\det M(\{U\}){\mathrm{e}}^{-S(\{U\})},$$

$$⟨O⟩=\frac{1}{Z}\int [\mathrm{d}U]O(\{U\})\det M(\{U\}){\mathrm{e}}^{-S(\{U\})}.$$

$$Z_0=\int [\mathrm{d}U]\det M_0(\{U\}){\mathrm{e}}^{-S(\{U\})}$$

$$Z=Z_0{⟨\frac{\det M(\{U\})}{\det M_0(\{U\})}⟩}_0,$$

$$Z=\int [\mathrm{d}U]\exp \left(\mathrm{tr}({\beta }_1{U}_1\cdots U_N+{\beta }_2{U}_N^{-1}\cdots U_1^{-1})\right),$$

$$⟨O⟩=\frac{1}{Z}\int [\mathrm{d}U]O(\{U\}){\mathrm{e}}^{-S(\{U\})},$$

$$S(\{U\})=-\mathrm{tr}({\beta }_1{U}_1\cdots U_N+{\beta }_2{U}_N^{-1}\cdots U_1^{-1}).$$

$$S(\{U\})=-2{\beta }_1\mathrm{Re}\mathrm{tr}(U_1\cdots U_N)\in \mathbb{R},$$

$$P(\{U\})=\frac{1}{Z}\exp (-S(\{U\}))$$

$$⟨O⟩\approx \frac{1}{K}\sum _{k=1}^{K}O(\{U^{(k)}\}).$$

$$\mathrm{d}{U}_k=-\sum _{a=1}^{n^2-1}\mathrm{i}{\lambda }_a[U_k{D}_{ak}S(\{U\})\mathrm{d}t+U_k\circ \mathrm{d}{w}_{ak}],k=1,\cdots ,N.$$

$$\mathrm{tr}({\lambda }_a{\lambda }_b)=2{\delta }_{ab},a,b=1,\cdots ,n^2-1.$$

$$D_{ak}f(\{U\})=\underset{ϵ\to 0}{\lim }\frac{f(\{U^{ϵ}\})-f(\{U\})}{ϵ},$$

$$\frac{\partial p}{\partial t}-D_{ak}[(D_{ak}S)p]=D_{ak}^{2}p,$$

$$U_k\leftarrow \exp \left(-\sum _{a=1}^{n^2-1}\mathrm{i}{\lambda }_a\left(D_{ak}S(\{U\})\Delta t+{\eta }_{ak}\sqrt{2\Delta t}\right)\right)U_k,k=1,\cdots ,N,$$

$$[SL(n,\mathbb{C})]^{N}=\left\{\{U\}\in(\mathbb{C}^{n\times n})^{N}\,\Bigg{|}\,U_{k}=\exp\left(\sum_{a=1}^{n^{2}-1}\mathrm{i}U_{ak}\lambda_{a}\right),~{}~{}U_{ak}\in\mathbb{C},~{}~{}k=1,\cdots,N\right\}.$$

$${\tilde{U}}_k:=V_k^{-1}{U}_k{V}_{k+1},k=1,\cdots ,N,$$

$$O(\{U\})=O(\{\tilde{U}\}),S(\{U\})=S(\{\tilde{U}\}).$$

$$\Vert \{U\}\Vert ={\left(\sum _{k=1}^N\Vert U_k{\Vert }_F^2\right)}^{1/2},\forall \{U\}\in [{\mathbb{C}}^{n\times n}{]}^N,$$

$$\Vert \{U\}\Vert ⩾\sqrt{Nn},\forall \{U\}\in [SL(n,\mathbb{C}){]}^N,$$

$$\frac{1}{n}\sum _{j=1}^n{\nu }_k^{(j)}⩾{\left(\prod _{j=1}^n{\nu }_k^{(j)}\right)}^{1/n},$$

$$\underset{\{V\}\in [SL(n,\mathbb{C}){]}^N}{\mathrm{argmin}}\Vert \{\tilde{U}\}{\Vert }^2.$$

$$V_k=\exp \left(\mathrm{i}\sum _{a=1}^{n^2-1}{V}_{ak}{\lambda }_a\right),$$

$$\left\{\begin{aligned} &\frac{\partial}{\partial X_{ak}}\|\{\widetilde{U}\}\|^{2}\bigg{|}_{\{\widetilde{U}\}=\{U\}}=0,\\ &\frac{\partial}{\partial Y_{ak}}\|\{\widetilde{U}\}\|^{2}\bigg{|}_{\{\widetilde{U}\}=\{U\}}=2\operatorname{tr}[\lambda_{a}(U_{k}U_{k}^{\dagger}-U_{k-1}^{\dagger}U_{k-1})],\end{aligned}\right.\qquad\forall a=1,\cdots,n^{2}-1,\quad\forall k=1,\cdots,N,$$

$$\left\{\begin{aligned} &\frac{\partial}{\partial X_{a}}\|\widetilde{U}\|^{2}\bigg{|}_{\widetilde{U}=U}=0,\\ &\frac{\partial}{\partial Y_{a}}\|\widetilde{U}\|^{2}\bigg{|}_{\widetilde{U}=U}=2\operatorname{tr}[\lambda_{a}(UU^{\dagger}-U^{\dagger}U)],\end{aligned}\right.\qquad\forall a=1,\cdots,n^{2}-1,$$

$$\mathrm{tr}[{\lambda }_a(U{U}^{†}-U^{†}U)]=0,\forall a=1,\cdots ,n^2-1.$$

$$U{U}^{†}-U^{†}U=\sum _{a=1}^{n^2-1}{\alpha }_a{\lambda }_a+\beta I.$$

$$\mathrm{tr}[{\lambda }_b(U{U}^{†}-U^{†}U)]=\sum _{a=1}^{n^2-1}{\alpha }_a\mathrm{tr}({\lambda }_b{\lambda }_a)+\beta \mathrm{tr}({\lambda }_b)={\alpha }_b,b=1,\cdots ,n^2-1,$$

$$\det (U{U}^{†})=\det (U^{†}U+\beta I)⩾\det (U^{†}U)+\beta \det (I)=\det (U{U}^{†})+\beta ,$$

Full text

Zhenning Cai’s work was supported by National University of Singapore Startup Fund under grant No. R-146-000-241-133. Yana Di’s work was supported by National Natural Science Foundation of China Nos. 91630208, 91641107, and 11771437. The authors would like to thank Dr. Lei Zhang in Department of Mathematics, National University of Singapore for useful discussions.

1. Introduction

In quantum chromodynamimcs (QCD), the renormalization of the coupling constant depends on the energy scale. As the energy scale increases, the coupling constant decays to zero. Therefore the perturbative theory works well for high-energy scattering. However, when studying QCD at small momenta or energies (less than 1GeV), the coupling constant is comparable to 1 and the perturbative theory is no longer accurate [13]. In this case, one of the important methods is the path integral formulation, in which people usually employ the lattice gauge theory to perform calculations. In lattice QCD, the degrees of freedom for both gluons and quarks are discretized on a four-dimensional lattice, whose grid points are

[TABLE]

where $a$ is the size of the lattice. Gluons on the lattice are represented by link variables between lattice points, which are matrices $U_{\mu}(\boldsymbol{x})\in SU(3)$, denoting the link between lattice points $\boldsymbol{x}$ and $\boldsymbol{x}+\boldsymbol{e}_{\mu}$. Since the degrees of freedom for quarks can usually be integrated out explicitly, the final form of the path integral is given by the following partition function:

[TABLE]

where $\int[\mathrm{d}U]$ stands for the integral with respect to all link variables $U_{\mu}(\boldsymbol{x})$ defined on the Haar measure of $SU(3)$, and $\{U\}$ represents the collection of all link variables. In the integrand, the matrix $M(\{U\})$ is the fermion Green’s function, and $S(\{U\})$ is the Euclidean action for the gluons. Thus, given an observable $O(\{U\})$, its expected value can be calculated by

[TABLE]

Monte Carlo methods such as the Metropolis algorithm and the Langevin algorithm can be applied to evaluate this integral.

When we consider the system with quark chemical potential, the term $\det M(\{U\})$ may be non-positive, and thus (1) is not a valid partition function [8]. In such a circumstance, the “reweighting” technique is required to carry out the Monte Carlo simulation. Such a method introduces another partition function

[TABLE]

and rewrite $Z$ as

[TABLE]

where $\langle\cdot\rangle_{0}$ denotes the expectation of $\cdot$ based on the partition function $Z_{0}$. However, due to the rapid change of sign in $\det M(\{U\})/\det M_{0}(\{U\})$, significant numerical sign problem may appear, causing large deviation in the numerical integration [7].

To relax the sign problem, numerical methods such as Lefschetz thimble method [6] and complex Langevin method (CLM) [20] have been introduced. This paper focuses on the CLM, which can be considered as a straightforward complexification of the real Langevin method. While the complex Langevin method effectively relaxes the sign problem for some problems, the behavior of this method seems quite unpredictable. Although in a lot of cases, the method produces correct integral values, sometimes it provides incorrect integral values, or even generates diveregent dynamics [4]. A number of efforts have been made to figure out the reason of failure and find a theory for its correct convergence [9, 10, 3, 17, 21]. Although a complete theory has not been found, the problem has been much better understood.

On the other hand, people are trying to stabilize CLM by various numerical techniques such as dynamical stabilization [5] and gauge cooling [23]. The dynamical stabilization adds a regularization term to the dynamics, which avoids divergence by introducing small biases; the gauge cooling technique makes use of the gauge symmetry of the system to stabilize CLM, and therefore does not introduce any biases. The gauge cooling technique has been successfully applied to a number of problems [24, 18, 1], which shows a considerable enhancement of the stability of CLM. Compared with a large number of studies for the original CLM, the analysis for the gauge cooling technique is still rarely seen. The justification of such a technique has been studied in [17, 19], while they do not explain how the complex Langevin process is stabilized. In this paper, we would like to carry out some numerical analysis for CLM with gauge cooling based on the one-dimensional $SU(n)$ theory. By such analysis, we would like to make the stabilizing effect more explicit. A specific analysis will be carried out for the Polyakov loop model. This model describes an infinitely heavy quark propagating only in time direction. Thus all the link variables can be labelled by a single subscript: $U_{k}$, $k=1,\cdots,N$, and we have periodic boundary conditions. Here we follow the reference [23] and choose the partition function to be

[TABLE]

where $\beta_{1}$ and $\beta_{2}$ are constants which are allowed to be complex. We are going to reveal some insights for the gauge cooling technique for this 1D model. Especially, when $n=2$, we will show that the complex Langevin process is localized when $\beta_{1}+\beta_{2}$ is not too large.

The rest of this paper is organized as follows. In Section 2, we briefly review the complex Langevin method and the gauge cooling technique. In Section 3, the optimal gauge transform is solved for one-dimensional periodic models. Section 4 is devoted to the study of gauge cooling after continuation of the numerical method, and some numerical experiments are carried out in Section 5. Finally, the paper ends with some concluding remarks in Section 6.

2. Complex Langevin method and gauge cooling

In this section, we provide a brief review of the CLM and the gauge cooling technique. For the sake of simplicity, we only introduce these methods for the one-dimensional lattice gauge theory,

2.1. Complex Langevin method

For the Polyakov chain model with partition function (2), the expected value of the observable is

[TABLE]

where $\{U\}=(U_{1},U_{2},\cdots,U_{N})\in[SU(n)]^{N}$ with $N$ being the number of lattice points. Here the Euclidean action $S(\{U\})$ is defined by

[TABLE]

Since $SU(n)$ has $n^{2}-1$ infinitesimal generators, the above integral is an $(n^{2}-1)N$-dimensional integral, which needs to be done by Monte Carlo method. When $\beta_{1}$ and $\beta_{2}$ are real and equal, we have

[TABLE]

which means that

[TABLE]

is a probability density function. In this case, to evaluate (3), we just need to draw $K$ samples $\{U^{(k)}\}$, $k=1,\cdots,K$ from the above probability density function, and then compute the average of the observable by

[TABLE]

Thus it remains only to provide a practical algorithm to draw the samples $\{U^{(k)}\}$.

A classical method to draw samples from a partition function is the Langevin method (see e.g. [13, Section 7.1.10]), which turns the sampling problem to solving a stochastic ordinary differential system, called Langevin dynamics, whose invariant measure is the desired probability density function. In our case, the Langevin dynamics has been given in [11], which reads

[TABLE]

Here $\{w_{ak}\}$ denotes the $N(n^{2}-1)$-dimensional Brownian motion, and the symbol “${}\!\circ\mathrm{d}w_{ak}$” stands for the Stratonovich interpretation of the stochastic differential equation (SDE). The matrices $\lambda_{a}$ denotes the infinitesimal generators of $SU(n)$ normalized by

[TABLE]

When $n=2$, they are usually chosen as Pauli matrices; when $n=3$, Gell-Mann matrices are a common choice. In the drift term, $D_{ak}$ is the left Lie derivative operator defined by

[TABLE]

where $\{U^{\epsilon}\}=(U_{1}^{\epsilon},\cdots,U_{N}^{\epsilon})$ and $U_{k}^{\epsilon}=\exp(\mathrm{i}\epsilon\delta_{kl}\lambda^{a}U_{k})$. Let $p(t,\{U\})$ be the probability density function at time $t$. Then $p(t,\{U\})$ satisfies the following Fokker-Planck equation:

[TABLE]

from which one can see that $P(\{U\})$ defined in (5) is a stationary solution of the above equation. Thus, when the stochastic process (7) is ergodic, i.e., there exists a unique invariant measure, then the sampling of $\{U\}$ can be implemented as follows:

Choose an arbitrary initial value $\{U\}=(U_{1},\cdots,U_{N})\in[SU(n)]^{N}$ and a proper time step $\Delta t$. Let $t=0$. Choose $T>0$ which indicates the time at which the probability density function $p(t,\cdot)$ is sufficiently close to the invariant measure $P(\cdot)$. 2. 2.

Evolve the solution by one time step:

[TABLE]

where $\eta_{ak}$ obeys the standard normal distribution. Set $t\leftarrow t+\Delta t$. 3. 3.

If $t<T$, return to step 2; otherwise, choose $\Delta T>0$ as the time interval between two adjacent samplings, and let $T\leftarrow T+\Delta T$. 4. 4.

Evolve the solution by one time step according to (10). 5. 5.

If $t>T$, take the current $\{U\}$ as a new sample and set $T\leftarrow T+\Delta T$. Stop if the number of samples is sufficient. 6. 6.

Return to step 4.

Here the equation (10) is the forward Euler discretization of the Langevin equation (7). Usually, when implementing the above algorithm, $O(\{U\})$ is immediately evaluated and added to the sum in (6), so that it is unnecessary to record all the samples.

In this paper, we are more interested in the case in which $\beta_{1}$ and $\beta_{2}$ are not equal, or even not real. In such a circumstance, the above derivation is no longer valid since $P(\{U\})$ is not a probability function. It has been introduced in the previous section that numerical sign problem may appear if the reweighting technique is used as a remedy. Interestingly, when the action goes complex, the algorithm introduced in the previous paragraph can still be applied if both the action $S(\cdot)$ and the observable $O(\cdot)$ can be holomorphically extended to the following domain of definition:

[TABLE]

Note that on the right-hand side of the above equation, if $U_{ak}\in\mathbb{C}$ is replaced by $U_{ak}\in\mathbb{R}$, then the set is exactly $[SU(n)]^{N}$. Therefore $[SL(n,\mathbb{C})]^{N}$ can be considered as a complexification of $[SU(n)]^{N}$. In this case, the algorithm with exactly the same procedure is called “complex Langevin method”. It can be formally shown that the complex Langevin method also produces an approximation of $\langle O\rangle$ [3, 17, 19, 21]. However, it has been demonstrated in [23, 22] that the complex Langevin method does not provide correct results when $N$ is large. This is also confirmed in our numerical experiments presented in Section 5.

2.2. Gauge cooling technique

It is generally accepted that when the stochastic variable $\{U\}$ drifts too far away from the space $[SU(n)]^{N}$, the computation tends to provide incorrect results or even fails to converge [17, 22]. In [23], a technique called “gauge cooling” was proposed to relax such a problem without introducing any biases to the system. The key idea of gauge cooling is to apply a proper gauge transform after every time step to keep $\{U\}$ close to $[SU(n)]^{N}$. Such a gauge transform is defined in $[SL(n,\mathbb{C})]^{N}$, which is called the “complexified gauge transform”. Due to the gauge symmetry of the lattice gauge theory, such a gauge transform changes neither the action nor the observable. Therefore formally, the process still gives correct expected value. This technique is applicable for the full four-dimensional case. In this paper, we only provide the details for the one-dimensional case with periodic boundary condition.

For the Polyakov chain model, for any field $\{U\}\in[SL(n,\mathbb{C})]^{N}$, the complexified gauge transform reads

[TABLE]

where $\{\widetilde{U}\}=(\widetilde{U}_{1},\cdots,\widetilde{U}_{N})$ is the field after transformation, and $\{V\}=(V_{1},\cdots,V_{N})\in[SL(n,\mathbb{C})]^{N}$ defines the gauge transform. Due to the periodic boundary condition, in (11), when $k=N$, the matrix $V_{N+1}$ is identical to $V_{1}$. The gauge symmetry says that for any $\{U\},\{V\}\in[SL(n,\mathbb{C})]^{N}$, we always have

[TABLE]

For any $\{U\}$ obtained from equation (10), our aim is to find $\{V\}\in[SL(n,\mathbb{C})]^{N}$ such that $\{\widetilde{U}\}$ is close to $[SU(n)]^{N}$. Different measurement of the closeness has been used in the literature [23, 17], and here we adopt the one proposed in [23], which is based on the following theorem:

Theorem 1**.**

Define

[TABLE]

where $\|\cdot\|_{F}$ is the Frobenius norm of the matrices defined by $\|U_{k}\|_{F}^{2}=\operatorname{tr}(U_{k}U_{k}^{\dagger})$. Then we have

[TABLE]

and the equality holds if and only if $\{U\}\in[SU(n)]^{N}$.

Proof.

For any $k=1,\cdots,N$, assuming $\nu_{k}^{(1)},\nu_{k}^{(2)},\cdots,\nu_{k}^{(n)}$ are all the eigenvalues of $U_{k}U_{k}^{\dagger}$ where $U_{k}\in SL(n,\mathbb{C})$, we have $\|U_{k}\|_{F}^{2}=\sum_{j=1}^{n}\nu_{k}^{(j)}$ and $\prod_{j=1}^{n}\nu_{k}^{(j)}=\det(U_{k}U_{k}^{\dagger})=\det(U_{k})\det(U_{k}^{\dagger})=1$. By the inequality of arithmetic and geometric means, we know

[TABLE]

and the equality holds if and only if $\nu_{k}^{(1)}=\nu_{k}^{(2)}\cdots=\nu_{k}^{(n)}$. Thus, $\|U_{k}\|_{F}^{2}\geqslant n$ and the equality holds if and only if all the eigenvalues of $U_{k}U_{k}^{\dagger}$ are $1$. Since any Hermitian matrix whose eigenvalues are all equal to $1$ must be the identity matrix, we conlude that $\|U_{k}\|_{F}^{2}=n$ if and only if $U_{k}\in SU(n)$. The conclusion of the theorem is then naturally obtained by summing up the square Frobenius norms of $U_{k}$. ∎

By this theorem, we can use $\sqrt{\|\{U\}\|^{2}-Nn}$ to characterize the distance between $\{U\}$ and $[SU(n)]^{N}$. Thus finding $\{V\}$ requires to solve the following optimization problem:

[TABLE]

Here $\{\widetilde{U}\}$ is defined by (11). This optimization problem is to be solved right after the evolution of the field. More precisely, in the algorithm described in the previous subsection, the following gauge cooling step is to be inserted after steps 2 and 4:

1. (GC)

Solve the optimization problem (12) and set $U_{k}\leftarrow V_{k}^{-1}U_{k}V_{k+1}$ for all $k=1,\cdots,N$.

Such a gauge cooling process has been formally justified in [19]. To solve the (12) numerically, the gradient descent method is proposed in [23]. By assuming

[TABLE]

where $V_{ak}=X_{ak}+\mathrm{i}Y_{ak}$ and $X_{ak},Y_{ak}\in\mathbb{R}$ for any $a$ and $k$, straightforward calculation yields

[TABLE]

using which the gradient descent method can be developed. In our work, we are going to show that in the one-dimensional periodic case, such an optimization problem can be solved exactly.

3. Gauge cooling in the one-dimensional case

This section is devoted to the exact solution of the optimization problem (12). We are going to start from the one-link case with $N=1$, and then generalize the result to the general one-dimensional case.

3.1. Introductory study: one-link case

If $N=1$, there is only one matrix in $\{U\}$. Hence the index $k$ for the matrix, which is always $1$, will be omitted. In this case, the gauge transform (11) turns out to be a similarity transform, and the gradient equation (13) can be simplified as

[TABLE]

When the optimal gauge transform is chosen, the above gradient must equal zero for every $a=1,\cdots,n^{2}-1$:

[TABLE]

Since $UU^{\dagger}-U^{\dagger}U$ is a Hermitian matrix, there exist real numbers $\alpha_{1},\cdots,\alpha_{n^{2}-1},\beta$ such that

[TABLE]

To determine the coefficients $\alpha_{a}$, we apply the property (8) to get

[TABLE]

from which one can see that $\alpha_{a}=0$ for all $a=1,\cdots,n^{2}-1$ when (14) is fulfilled. Meanwhile, we also have $\beta=0$ since if $\beta$ were positive, we would get

[TABLE]

indicating that $\beta\leqslant 0$, which is a contradiction; a similar contradiction can be obtained if we assume $\beta<0$. Conclusively, after the optimal gauge transform, the link variable $U$ satisfies $UU^{\dagger}=U^{\dagger}U$, which means that $U$ is a normal matrix.

A matrix is a normal matrix if and only if it can be diagonalized by a unitary matrix. Therefore, for any $U\in SL(n,\mathbb{C})$, when $V^{-1}UV$ gives the optimal gauge transform, there exists a diagonal matrix $\Lambda$ and a unitary matrix $P$ such that

[TABLE]

Note that the diagonal entries of $\Lambda$ are the eigenvalues of $U$. When $U$ is diagonalizable, i.e., there exists $Q\in SL(n,\mathbb{C})$ such that $U=Q\Lambda Q^{-1}$, one can find that $V=QP$. Conversely, for a diagonalizable $U$, we can choose an arbitrary unitary matrix $S$, and the gauge transform given by $V=QP$ always transforms $U$ to a normal matrix. A simple choice is $P=I$, for which $\widetilde{U}=V^{-1}UV=\Lambda$.

The above analysis shows that in the one-link case, for every time step, we can choose the gauge transform such that the link $U$ is always a diagonal matrix. Thereby, it is possible to write down a new stochastic process for this diagonal link, and get a clearer picture how gauge cooling helps pull back the complexified links. Before that, we are going to generalize the above results to the 1D case with periodic boundary conditions, and make the analysis more rigorous. Interestingly, in the general 1D case, the link variables $U_{k}$ can still be maintained as diagonal by choosing optimal gauge transforms. This will be detailed in the following subsection.

3.2. General 1D case

Now we are going to restore the index $k$ and study the optimization problem (12) for any $N$. To begin with, we show that equating (13) to zero is equivalent to solving the optimization problem:

Lemma 1**.**

Suppose $\{U\}$ satisfies

[TABLE]

Then for any gauge transform defined by (11), we have $\|\{\widetilde{U}\}\|\geqslant\|\{U\}\|$.

Proof.

It suffices to show that the Hessian matrix is positive semidefinite everywhere, so that the first-order optimality condition (15) gives the global minimum. Again we let $V_{ak}=X_{ak}+\mathrm{i}Y_{ak}$. Because $\partial_{X_{ak}}\|\widetilde{U}\|^{2}$ equal zero at $\widetilde{U}=U$ for any $U$, we only need to calculate the second-order deriatives of $\|\{\widetilde{U}\}\|^{2}$ with respect to the imaginary parts $Y_{ak}$. For any $a,b=1,\cdots,n^{2}-1$ and $k,l=1,\cdots,N$, by straightforward calcuation, we obtain

[TABLE]

To show the convexity of $\|\{\widetilde{U}\}\|^{2}$ as a function of $Y_{ak}$, we just need to show that

[TABLE]

for any $v_{ak}\in\mathbb{R}$ (note that $H_{akbl}$ is real). To prove this, we calculate $R$ as follows:

[TABLE]

For simplicity, we define $M_{k}=\sum_{a=1}^{n^{2}-1}v_{ak}\lambda_{a}$, which is Hermitian, and thus the above equation becomes

[TABLE]

Because of the periodic boundary conditions for $\{U\}$, it is obvious that

[TABLE]

Therefore, we have

[TABLE]

which indicates that $H_{akbl}(\{U\})$ is positive semidefinite. ∎

By the above lemma, we only need to focus on the algebraic equations (15) with $\{U\}$ replaced by $\{\widetilde{U}\}$. In order to apply the idea used in the one-link case, we also write the gauge transform as a simiarity transform by introducing the block matrix $\mathcal{U}$, defined by

[TABLE]

which is an $nN\times nN$ matrix with $N\times N$ blocks. Thereby, the gauge transform

[TABLE]

can be written as

[TABLE]

where $\mathcal{V}=\operatorname{diag}(V_{1},V_{2},\cdots,V_{N})$. Below we claim that similar to the one-link case, the optimal gauge is achieved when $\mathcal{U}$ is a normal matrix.

Theorem 2**.**

The block matrix $\mathcal{U}$ defined by (16) is normal if and only if $\{U\}$ satisfies (15).

Proof.

By the same argument as in the one-link case, we know that (15) is equivalent to

[TABLE]

Since

[TABLE]

the equivalence stated in the theorem is manifest. ∎

To proceed, we again follow the one-link case and try to find the optimal gauge transform by diagonalizing the matrix $\mathcal{U}$. The result is given in the following theorem:

Theorem 3**.**

Suppose $U_{N}U_{1}\cdots U_{N-1}$ can be diagonalized by

[TABLE]

where $\Lambda=\operatorname{diag}(\mu_{1},\cdots,\mu_{n})$ and $\det Q=1$. Then $\mathcal{U}$ can be diagonalized by

[TABLE]

where

[TABLE]

and $\mu_{j}^{(1)},\cdots,\mu_{j}^{(N)}$ are all the distinct $N$th roots of $\mu_{j}$.

Proof.

We first prove that the diagonal entries of $\Lambda_{k}$, $k=1,\cdots,N$ provide all the eigenvalues of $\mathcal{U}$. In other words, we are going to show that $\nu$ is an eigenvalue of $\mathcal{U}$ if and only if $\nu^{N}=\mu_{j}$ for some $j\in\{1,\cdots,n\}$. If $\nu$ is an eigenvalue of $\mathcal{U}$, since $\det\mathcal{U}=(-1)^{N+1}\det(U_{1}U_{2}\cdots U_{N})=(-1)^{N+1}$, we have $\nu\neq 0$. Let $x=(x_{1}^{T},x_{2}^{T},\cdots,x_{N}^{T})^{T}$, where $x_{k}\in\mathbb{C}^{n}$ for all $k=1,\cdots,N$, be the associated eigenvector. By $\mathcal{U}x=\nu x$, we get

[TABLE]

Concatenating these equations yields

[TABLE]

which shows that $\nu^{N}$ is the eigenvalue of $U_{N}U_{1}\cdots U_{N-1}$. On the contrary, for any $j=1,\cdots,n$, suppose $x_{N}\in\mathbb{C}^{n}\setminus\{0\}$ is the eigenvector of $U_{N}U_{1}\cdots U_{N-1}$ associated with the eigenvalue $\mu_{j}$. When $\nu^{N}=\mu_{j}$, we can define $x_{1},\cdots,x_{N-1}$ by equations (21). Let $x=(x_{1}^{T},x_{2}^{T},\cdots,x_{N}^{T})^{T}$, we get $\mathcal{U}x=\nu x$. Therefore, $\nu$ is an eigenvalue of $\mathcal{U}$.

Because $\Lambda_{k}^{N}=\Lambda$ for all $k=1,\cdots,N$ based on the equation (19), it follows from (17) that $U_{N}U_{1}\cdots U_{N-1}Q\Lambda_{k}^{-(N-1)}=Q\Lambda_{k}$. Using these equality, one obtains by straightforward calculation that

[TABLE]

Laying the above equations for all $k=1,\cdots,N$ in a row and dividing the result by $\sqrt{N}$, we get

[TABLE]

which concludes the proof. ∎

By now, it is clear that

[TABLE]

By Theorem 2, we need to find $\mathcal{V}$ such that the above matrix is normal. This requires us to find $\mathcal{V}$ such that $\mathcal{V}^{-1}\mathcal{T}\mathcal{S}$ is a unitary matrix. Fortunately, one can easily achieve this by the following property of $\mathcal{S}$:

Lemma 2**.**

The matrix $\mathcal{S}$ defined by (20) satisfies

[TABLE]

Proof.

It is easy to get

[TABLE]

Now we calculate the matrix blocks. By $\Lambda_{k}^{N}=\Lambda$, we know that $(\Lambda_{k}\Lambda_{k}^{\dagger})=(\Lambda\Lambda^{\dagger})^{\frac{1}{N}}$. For any $p,q=0,-1,\cdots,-(N-1)$ and $p\neq q$,

[TABLE]

and

[TABLE]

On the right-hand side of (23), the sum in the parentheses is

[TABLE]

To calculate it, for all $j=1,\cdots,n$, we write $\mu_{j}$ as $\mu_{j}=|\mu_{j}|\exp(\mathrm{i}\theta_{j})$, and assume

[TABLE]

Thus

[TABLE]

The above equation indicates that $\sum_{k=1}^{N}\Lambda_{k}^{p}(\Lambda_{k}^{q})^{\dagger}=0$. Therefore

[TABLE]

Based on Lemma 2, we provide in the following theorem a series of convenient choices of optimal gauge transforms, which are obviously a generalization of the one-link case:

Theorem 4**.**

Suppose $U_{N}U_{1}\cdots U_{N-1}$ is diagonalizable. For any $P_{k}\in SL(n,\mathbb{C})$, $k=1,\cdots,N$, define

[TABLE]

We have that $\mathcal{V}^{-1}\mathcal{T}\mathcal{S}$ is unitary. In this case,

[TABLE]

Proof.

With $\mathcal{V}=\operatorname{diag}(V_{1},\cdots,V_{N})$, it is easy to calculate

[TABLE]

by Lemma 2 and the equation (18).

Due to $\mathcal{U}\mathcal{U}^{\dagger}=\operatorname{diag}(U_{1}U_{1}^{\dagger},\cdots,U_{N}U_{N}^{\dagger})$, we have

[TABLE]

For $\widetilde{\mathcal{U}}=\mathcal{V}^{-1}\mathcal{T}\mathcal{S}\operatorname{diag}(\Lambda_{1},\cdots,\Lambda_{N})(\mathcal{V}^{-1}\mathcal{T}\mathcal{S})^{\dagger}$ with the equation (22), there is

[TABLE]

Yet,

[TABLE]

In practice, we can simply choose $P_{k}=I$, and thus

[TABLE]

from which we can see that by choosing appropriate gauge transform, all the links are always diagonal. This shows us how the gauge cooling technique helps prevent the field links from excursing too far away from $SU(n)$. Gauge cooling removes redundant degrees of freedom, so that the probability of imaginary excursion is greatly reduced. Therefore gauge cooling can be considered as a method of dimension reduction to enhance the numerical stability. Here we would like to comment that in the multidimensional case, the exact solution no longer has the simple diagonal form as (25). But one can still expect the similar effect. In the next section, we are going to present the effect of gauge cooling on the drift term. Before that, we will first pick up a special case left out in the above discussion.

By now, we have only considered the case when the matrix $U_{N}U_{1}\cdots U_{N-1}$ is diagonalizable. When it is not diagonalizable, we consider its Jordan normal form given by

[TABLE]

where $\mu_{m},m=1,2,\cdots,n$ are eigenvalues of $U_{N}U_{1}\cdots U_{N-1}$, and $\delta_{m},m=1,2,\cdots,n-1$ are [math] or $1$. At least one of $\delta_{m}$, $m=1,2,\cdots,n-1$ is nonzero. In this case, we choose a positive number $\epsilon\ll 1$ and define $Q=\operatorname{diag}\{1,\epsilon,\cdots,\epsilon^{n-1}\}Q^{\prime}$. Thus

[TABLE]

where the middle matrix on the right-hand side is very close to a diagonal matrix. By scaling $Q^{\prime}$ properly such that $\det Q=1$ and choosing $V_{k}$ according to (24), the matrix $\widetilde{\mathcal{U}}$ can be very close to a normal matrix. Since $\epsilon$ can be arbitrarily small, the norm $\|\{\widetilde{U}\}\|$ can be arbitrarily close to $\sqrt{N\left(\sum_{j=1}^{n}|\mu_{j}|^{\frac{2}{N}}\right)}$, and the new links $\{\widetilde{U}\}$ can be arbitrarily close to (25). In fact, this is the situation when the optimal gauge transform does not exist, and $\sqrt{N\left(\sum_{j=1}^{n}|\mu_{j}|^{\frac{2}{N}}\right)}$ is the infimum of all possible values of $\|\{\widetilde{U}\}\|$. However, in $SL(n,\mathbb{C})$, nondiagonalizable matrices only form a submanifold with zero measure. Therefore, in our implementation, we just set the new matrices to be (25) without checking the diagonalizablity of $U_{N}U_{1}\cdots U_{N-1}$.

4. Study of gauge cooling by stochastic differential equations

To further explore the effect of gauge cooling, we will fix the choice of gauge transforms as (25) for all time steps. Thereby, all the links are always diagonal matrices, and in fact, these links have only $n-1$ degrees of freedom, which are $\mu_{1},\cdots,\mu_{n-1}$ (note that $\mu_{1}\cdots\mu_{n}=1$). Thus, it is helpful to study the dynamics of these quantities directly, by taking an infinitesimal time step. Below we will first give the general formulation for the $SU(n)$ case, and then focus on the special action (4) for the $SU(2)$ case.

4.1. General case in $SU(n)$

Our analysis in Section 3.2 shows that optimal gauge cooling can maintain the following structure of $\{U\}$:

[TABLE]

where $\Lambda=\operatorname{diag}(\mu_{1},\cdots,\mu_{N})$. Before deriving the dynamics of $\mu_{j}$, we first review the steps to evolve the links by one time step $\Delta t$:

Let $U_{k}^{*}=\exp\left(-\sum_{a=1}^{n^{2}-1}\mathrm{i}\lambda_{a}(K_{ak}\Delta t+\eta_{ak}\sqrt{2\Delta t})\right)U_{k}$, where $K_{ak}=D_{ak}S(\{U\})$ and $\eta_{ak}$ obeys the standard normal distribution. 2. 2.

Compute the eigenvalues of $U_{N}^{*}U_{1}^{*}\cdots U_{N-1}^{*}$, and use $\mu_{1}^{*},\cdots,\mu_{N}^{*}$ to denote the results. 3. 3.

Update the matrices by setting $U_{1},\cdots,U_{N-1}$ to be $\operatorname{diag}(|\mu_{1}^{*}|^{\frac{1}{N}},\cdots,|\mu_{n}^{*}|^{\frac{1}{N}})$, and setting $U_{N}$ to be $\operatorname{diag}(|\mu_{1}^{*}|^{-\frac{N-1}{N}}\mu_{1}^{*},\cdots,|\mu_{n}^{*}|^{-\frac{N-1}{N}}\mu_{n}^{*})$.

Here the first step is identical to the evolution without gauge cooling, and the second and third steps are an application of equation (25), which performs the optimal gauge cooling. Next, we are going to replace $\Delta t$ in the above equation by an infinitesimal time step $\mathrm{d}t$, and $\eta_{ak}\sqrt{2\Delta t}$ will be replaced by $\mathrm{d}w_{ak}$ correspondingly. Such a limit can be more easily taken if we consider the above process as the Euler-Maruyama method for Itô SDE (see e.g. [14, Section 3.6] for the description of numerical methods for Itô and Stratonovich SDEs). Therefore in this section, we will switch to Itô’s representation, where the circle sign “$\circ$” in front of $\mathrm{d}w_{ak}$ is removed. For easier understanding, below we are going to perform a formal derivation of the SDEs for $\mu_{j}$, $j=1,\cdots,n$, where notations such as $\mathrm{d}t$ and $\mathrm{d}w_{ak}$ will be used directly to replace $\Delta t$ and $\eta_{ak}\sqrt{2\Delta t}$, and all the terms with magnitude smaller than $\mathrm{d}t$ will be discarded without explicit approximation. We claim that such a process can be made rigorous following the classical error analysis for the Euler-Maruyama method without much difficulty [16].

We first write down the infinitesimal version of Step 1 using Itô calculus:

[TABLE]

where we have used $(\mathrm{d}w_{ak})^{2}=2\,\mathrm{d}t$ and

[TABLE]

which can be found in [13]. Step 2 requires us to compute the product of these matrices:

[TABLE]

where

[TABLE]

Next, we need to find all its eigenvalues. Suppose $r_{j}^{*}$ is the unit eigenvector associated with the eigenvalue $\mu_{j}^{*}$:

[TABLE]

By (26), as $U_{N}^{*}U_{1}^{*}\cdots U_{N-1}^{*}$ is a perturbation of the diagonal matrix $\Lambda$, the eigenvalue $\mu_{j}^{*}$ should be the perturbation of the diagonal entry $\mu_{j}$, and $r_{j}^{*}$ should be the perturbation of the canonical basis $e_{j}$. Assume

[TABLE]

Substituting (28) and (26) into (27) yields

[TABLE]

Now we equate the coefficients of the same infinitesimal terms to get

[TABLE]

and

[TABLE]

For the equation (29), we multiply both sides by $e_{j}^{T}$ on the left, and obtain

[TABLE]

Note that $e_{j}^{T}\Lambda=e_{j}^{T}\mu_{j}$. We get $c_{j}^{ak}$ by cancelling out the terms $\mu_{j}e_{j}^{T}h_{j}^{ak}$ on both sides:

[TABLE]

Because $c_{j}^{ak}$ is independent of $k$, we rewrite it as $c_{j}^{a}$ hereafter. If we multiply both sides of (29) by $e_{l}^{T}$ on the left for $l\neq j$, the result is

[TABLE]

When $\mu_{j}\neq\mu_{l}$, we have

[TABLE]

where $(h^{ak}_{j})_{l}$ denotes the $l$th component of $h^{ak}_{j}$. In order to determine the $j$th component of $h^{ak}_{j}$, we impose the following two conditions:

1. (1)

$(r_{j}^{*})^{\dagger}r_{j}^{*}=1$; 2. (2)

The $j$th component of $r_{j}^{*}$ is real.

Since $r_{j}^{*}$ is a perturbation of $e_{j}$, the $j$th component of $r_{j}^{*}$ must be nonzero. Thus the above two constraints fix the eigenvector. By the first constraint, we have

[TABLE]

Therefore,

[TABLE]

The second constraint requires that the imaginary part of $(h_{j}^{ak})_{j}$ be zero as well. Thus $(h_{j}^{ak})_{j}=0$. Now, we can write $h_{j}^{ak}$ explicitly as

[TABLE]

where $k=1,\cdots,N-1$, and

[TABLE]

By these results, we can find $\sigma_{j}$ by multiplying both sides of (30) by $e_{j}^{T}$ on the left:

[TABLE]

where

[TABLE]

To summarize, we insert (31) and (32) into (28), which gives the value of $\mu_{j}^{*}$:

[TABLE]

By now, Step 2 is accomplished.

Step 3 simply means that $\mu_{j}^{*}$ is just $\mu_{j}+\mathrm{d}\mu_{j}$. Hence,

[TABLE]

To make further simplification, we define

[TABLE]

Thus (33) can be rewritten as

[TABLE]

This equation shows that by gauge cooling, the one-dimensional $N$-link model is essentially equivalent to the one-link model, as is known for the exact solution [23]. The optimal gauge cooling automatically utilizes this property, and therefore greatly reduces the difficulty of the simulation. This analysis further confirms the nature of dimension reduction for the gauge cooling technique from another point of view, which also applies in the multi-dimensional cases.

Some other interesting properties can be observed from (34). For example, we can see that only $n-1$ drift terms are active for the Langevin process. This is clear for $n=2$ since the generators of $SU(2)$ are just Pauli matrices (below we add the superscript $(n)$ to matrices $\lambda_{a}$ for clarification):

[TABLE]

and we have $e_{j}^{T}\lambda_{1}^{(2)}e_{j}\equiv 0$ and $e_{j}^{T}\lambda_{2}^{(2)}e_{j}\equiv 0$, which means that both $K_{1}$ and $K_{2}$ are multiplied by zero in (34). When $n=3$, the generators of $SU(3)$ are Gell-Mann matrices:

[TABLE]

from which one sees that $e_{j}^{T}\lambda_{a}^{(3)}e_{j}\equiv 0$ if $a\neq 3$ and $a\neq 8$. Therefore only $K_{3}$ and $K_{8}$ are active. Likewise, for $n\geqslant 3$, the first $n(n-2)$ generators of $SU(n)$ can be chosen as $\lambda_{a}^{(n)}=\operatorname{diag}(\lambda_{a}^{(n-1)},0)$. For the rest generators, the only one with nonzero diagonal entries is [12]

[TABLE]

Therefore, by mathematical induction, the number of active drift terms in (34) equals $n-1$.

In (34), except the term with $K_{a}$, the other two drift terms are independent of the action $S(\{U\})$. These terms also help us stabilize the complex Langevin process, as will be demonstrated by a specific example in the next subsection.

4.2. Polyakov loop model for the $SU(2)$ theory

This subsection will be devoted to a special case study. For easier demonstration, we choose $n=2$ since in this case, there is only one degree of freedom, namely $\mu_{1}$ (or $\mu_{2}$), as stated in the begining of Section 4. By taking this choice, it is much easier to understand the complex Langevin dynamics, since one can explicitly plot the flow field and the distribution of samples on the complex plane, giving a clear illustration of possible issues. Such a method has been widely used in existing investigations of the complex Langevin method [2, 17, 21].

Based on the discussion of the previous subsection, we just need to consider the one-link model with $N=1$. Since $\mu_{1}\mu_{2}=1$, here we only study the dynamics of $\mu_{1}$. For simplicity, we omit the subscript and write $\mu_{1}$ as $\mu$. Thus, by equation (33), we can write down the SDE for $\mu$ as follows:

[TABLE]

For such a representation, $U\in SU(2)$ corresponds to $\mu\in U(1)$. Therefore it is natural to make the change of variable $\mu=\exp(-\mathrm{i}s)$ and consider the SDE for $s$, which is

[TABLE]

The above equation can be verified by Itô’s calculus:

[TABLE]

which agrees with (36). Note that $\mu$ is $2\pi$-periodic with respect to the real part of $s$. Therefore in (37), we need to consider $s(t)$ as a stochastic process defined on the cylinder $(\operatorname{Re}s(t),\operatorname{Im}s(t))\in\mathbb{S}^{1}\times\mathbb{R}$. When $K_{3}=0$, the flow field given by the drift term is plotted in Figure 1, from which we can see that when $s$ is not real, the flow pushes $s$ toward the real axis, meaning that the cotangent term in (37) also helps stabilize the complex Langevin method. This also implies that when $K_{3}$ does not provide a strong repulsion from the real axis, the variable $s$ can stay in a region near the real axis, leading to a localized distribution. Below we are going to show this behavior rigorously for the action (4).

From the action $S(U)$ given by (4) and the definition of the Lie derivative (9), it is easy to get

[TABLE]

For simplicity, we let $\mathrm{d}w=\mathrm{d}w_{3}$ and $\beta=\beta_{1}+\beta_{2}$. Noting that $\mu^{-1}-\mu=2\mathrm{i}\sin s$, we can rewrite the equation (36) as

[TABLE]

To write down the above equation more explicitly, we let $s=x+\mathrm{i}y$ and $\beta=A+\mathrm{i}B$, where $x,y,A,B\in\mathbb{R}$. Then $K_{3}$ can be expanded as

[TABLE]

Below we use $K_{R}$ and $K_{I}$ to denote the real and imaginary parts of $K_{3}$. Then the SDE (38) is equivalent to the following 2D SDEs:

[TABLE]

where we remind the readers that $x$ is a $2\pi$-periodic variable and $w$ is the real Brownian motion.

Below we would like to claim that when $S(U)$ is small, in the complex Langevin process (40), the imaginary part of $s$ cannot excurse to infinity due to the cotangent term in (37) if the initial $y$ is zero. This is detailed in the following theorem:

Theorem 5**.**

Suppose $|A|<\frac{3\sqrt{3}}{2}$. There exists $\epsilon_{A}>0$ such that when $B\in(-\epsilon_{A},\epsilon_{A})$, we can find constants $C_{1}$ and $C_{2}$ satisfying all the following conditions:

• •

$C_{2}>C_{1}\geqslant 0$;

• •

$K_{I}<0$* for all $x\in\mathbb{R}$ and $y\in(C_{1},C_{2})$;*

• •

$K_{I}>0$* for all $x\in\mathbb{R}$ and $y\in(-C_{2},-C_{1})$.*

Here $K_{I}$ is given in the second equation of (40).

Proof.

By the definition of $K_{I}$, we see that it depends on $x,y,A,B$, and it satisfies the following identities:

[TABLE]

Such a symmetry allows us to focus only on the case with $y<0$, $A\geqslant 0$ and $B>0$. Let $\xi=\cos x$ and $\eta=-\tanh y$ which satisfy $-1\leqslant\xi\leqslant 1$ and $0<\eta<1$. Using these new variables, we can rewrite the inequality $K_{I}>0$ in the following equivalent form:

[TABLE]

Below we consider the cases $A=0$ and $A>0$ separately.

If $A=0$, in order that (41) holds for all $x\in\mathbb{R}$, we need to ensure that

[TABLE]

Clearly, the above inequality holds for $B=0$ and any $\eta=1/2$. Thus, by the continuity of both sides of (42), we can find a neighborhood of $(0,1/2)$, say $(-\epsilon_{0},\epsilon_{0})\times(1/2-\delta,1/2+\delta)$, such that (42) still holds, which indicates $K_{I}>0$. Thus, by the relation between $y$ and $\eta$, the constants $C_{1}$ and $C_{2}$ can be found by $C_{1}=\operatorname{arctanh}(1/2-\delta)$ and $C_{2}=\operatorname{arctanh}(1/2+\delta)$.

If $A>0$, the equation (41) can be written as

[TABLE]

Similar to the previous case, we just need to find some $\eta>0$ such that the above inequality holds for $B=0$ and any $\xi\in[-1,1]$. This condition is equivalent to

[TABLE]

When $0<\eta<1$, it is not difficult to work out that $\max\limits_{\xi\in[-1,1]}[\xi^{3}(1-\eta^{2})-\xi]=\frac{2}{3}\sqrt{\frac{1}{3(1-\eta^{2})}}$. Thus a simple rearrangement of (44) yields

[TABLE]

By the assumption $0<A<3\sqrt{3}/2$, we know that the right-hand side of the above equation is positive, and therefore a positive $\eta$ satisfying (44) exists. The existence of $\epsilon_{A}$, $C_{1}$ and $C_{2}$ follows the same argument as in the case $A=0$. ∎

Theorem 5 shows that when $|A|<3\sqrt{3}/2$ and $B\in(-\epsilon_{A},\epsilon_{A})$, the solution of (40) satisfies $y(\eta)\in(-C_{2},C_{2})$ if the initial condition satisfies $y(0)\in(-C_{2},C_{2})$. Recall that the formal justification of the complex Langevin method requires that $y(0)=0\in(-C_{2},C_{2})$. Therefore Theorem 5 does provide us a safe region of parameters, for which the imaginary part cannot excurse out of the striplike area $|y|<C_{2}$. A more precise description of this safe region can be obtained by (43), which inspires us to define

[TABLE]

so that for any given $A\in(-3\sqrt{3}/2,3\sqrt{3}/2)$, the maximum value of $\epsilon_{A}$ satisfies $f(A,\epsilon_{A})=0$. By this argument, the range of $A$ and $B$ which localizes the distribution function can be numerically computed, and we plot such a region in Figure 2.

Such a localizing phenomenon has been studied in [2], where the striplike region is the result of the competition of two parts of the action: one part of the action attracts the Langevin process toward the real axis, while the other part repulses the process away from the real axis. In our work, we have shown that gauge cooling provides an additional attractive drift, which helps the formation of the confining strip even if the action does not contain an attractive part.

One side effect of this advantage of gauge cooling is the existence of singularity at points $s=k\pi$ for $k\in\mathbb{Z}$, corresponding to the case $\mu_{1}=\mu_{2}$. For instance, the SDE around $s=0$ behaves asymptotically like

[TABLE]

Such a velocity field does not even fit in the $L^{2}$ space, causing that most techniques to prove the ergodicity fail to work. Theoretically, it is yet unclear whether such a singularity will lead to the failure of the complex Langevin method, while our numerical experiments seem to suggest its mildness.

5. Numerical examples

In this section, we give two examples to validate our numerical analysis. In the first example, we apply different gauge cooling methods to the $SU(3)$ theory to compare their results; in the second example, we simulate (40) directly for different parameters to examine the convergence behaviors.

5.1. Polyakov loop model for $SU(3)$ theory

In this example, we apply the complex Langevin method with gauge cooling to the Polyakov loop model (4) for $SU(3)$ theory. The aim is to show the effectiveness of optimal gauge cooling in reducing the norm $\|\{U\}\|$. Following [23], we choose $\beta_{1}=\beta+\kappa e^{\mu}$ and $\beta_{2}=\bar{\beta}+\kappa e^{-\mu}$ in (4), and set $\beta=2$, $\kappa=0.1$ and $\mu=1$. In the following numerical experiments, the following three methods will be tested:

• •

Complex Langevin with no gauge cooling.

• •

Complex Langevin with optimal gauge cooling given by (25).

• •

Complex Langevin with gauge cooling implemented by gradient descent method [23].

In the last method, we follow [23] to assume that $V_{ak}$ are purely imaginary ($V_{ak}=\mathrm{i}Y_{ak}$), and the gradient has been given in the second equation of (13). Thus every gradient descending step is given by

[TABLE]

Here $\lambda_{a}$ is the Gell-Mann matrices given in (35), and $\tilde{\alpha}$ indicates the step length. Since gauge cooling is to be done after every time step, and the number of time steps may be large, here we apply neither sophisticated line search technique nor accurate stopping criterion. The step length is simply chosen as $\tilde{\alpha}=\alpha\Delta t$ for some constant $\alpha$. After every time step of complex Langevin, the gradient-descent based gauge cooling (45) is applied five times. Two choices of $\alpha$ ($0.4$ and $1.0$) are tested.

In our numerical tests, we set the Langevin time step to be $\Delta t=2\times 10^{-5}$, and we apply the above methods until $t=10$. The evolution of the norm $\|\{U\}\|$ for both cases $N=16$ and $N=32$ is provided in Figure 3, where the horizontal axis is the Langevin time, and the vertical axis represents the value of $\Delta F:=\|\{U\}\|^{2}-Nn$, which is zero if every $U_{k}$ is unitary. Figure 3 only provides one realization of the complex Langevin process in each case. In our experiments, despite the existence of the randomness, no essential difference can be observed in different runs. Figure 3 shows that if no gauge cooling is used, the numerical method fails due to a rapidly increasing $\Delta F$, and all other three lines show converging results. For $N=16$, the gradient descent methods show a quite good performance. The value of $\Delta F$ are well controlled below $10^{-2}$ for all time steps, and two different step lengths $\alpha=0.4$ and $\alpha=1.0$ show very similar results. However, when $N$ increases to $32$, the step length given by $\alpha=1.0$ shows much better performance. When $t>7$, even for $\alpha=1.0$, the value of $\Delta F$ can reach a number larger than $1$. Contrarily, if the optimal gauge cooling is used, the value of $\Delta F$ is just oscillating around $10^{-6}$, which confirms our statement that the process is essentially independent of $N$ if the optimal gauge cooling is used. These facts indicate that when the number of links is large, one may need to design a more efficient optimization method to suppress the norm $\|\{U\}\|$, so that the power of gauge cooling can be freed maximally.

To verify the correctness of the numerical method, we have calculated the observables $O_{k}(\{U\})=\operatorname{tr}\left[(U_{1}\cdots U_{N})^{k}\right]$ for $k=\pm 1,\pm 2,\pm 3$. Because the Polyakov loop model reduces to a one-link integral analytically [23], we can calculate the exact values of these observables by considering the case $N=1$. For any $U\in SU(3)$, its eigenvalues can be written as $e^{\mathrm{i}\phi_{1}},e^{\mathrm{i}\phi_{2}}$ and $e^{-\mathrm{i}(\phi_{1}+\phi_{2})}$ where $\phi_{1},\phi_{2}\in(-\pi,\pi]$. By Weyl’s integral theorem [15], when $N=1$, we can calculate the integrals $\langle O_{k}\rangle$ as

[TABLE]

where $\eta$ is a normalization constant given by

[TABLE]

and

[TABLE]

Thus, the value of $\langle O_{k}\rangle$ can be obtained by two-dimensional numerical integration, which is used to verify the results from stochastic simulations. In Table 1, some results have been listed. Because the value of $\langle O_{k}\rangle$ is real, we just list the real part of the numerical result. Note that the results for CLM without gauge cooling is not listed, since the computation always fails due to the rapidly growing norm of $\|\{U\}\|$. For $N=16$, all the three methods give reasonable approximations of the exact integral value, while for $N=32$, we observe large deviations from the exact values when the gauge cooling is insufficient ($\alpha=0.4$). In fact, we even experienced divergence in our test runs when $\alpha=0.4$, which again emphasizes the importance to use an effective gauge cooling method.

5.2. Optimal gauge cooling in the $SU(2)$ theory

In this example, we focus on the behavior of equation (40) for different values of $A$ and $B$. For all simulations, we start with the initial value $x=0.5$ and $y=0$, and evolve the process with time step $\Delta t=10^{-5}$. The first $3\times 10^{5}$ steps are considered as “unsteady calculations”, and afterwards, we proceed the evolution by $10^{8}$ steps, in which the samples are drawn at every $10^{4}$ steps.

The observables we compute are still $O_{k}(U)=\operatorname{tr}(U^{k})$, where $U=\operatorname{diag}(\mathrm{e}^{-\mathrm{i}s},\mathrm{e}^{\mathrm{i}s})$. The results for $k=1,2,3$ are provided in Table 2, and we also plot the distributions of the samples $(x,y)$ in Figure 4. The following four choices of $A$ and $B$ are considered in our results:

• •

$A=1$ and $B=0.2$: In this case, the conclusion of Theorem 5 holds, and therefore all the samples are confined in the strip-like region. The simulation is stable and provides reliable values of the observables.

• •

$A=1$ and $B=2$: The conclusion of Theorem 5 does not hold in this case, so the possible values of $y$ can be arbitrary large. The numerical values of the observables also show an interesting behavior. When $k$ increases, the numerical result shows larger and larger deviation from the exact value. In fact, different simulations generate very different results for $k=2$ and $k=3$, while the results for $k=1$ seem stable in different runs. This behavior may be due to insufficient decay in the invariant measure, which is studied in detail in [21].

• •

$A=5$ and $B=1$: Since $A>3\sqrt{3}/2$, Theorem 5 does not hold. However, the value of $B$ is relatively small compared with $A$, which provides only a small channel to allow $y$ to drift far from the real axis. Therefore, in our numerical results, the distribution of the samples still looks as if they are confined. The simulation is again stable and the numerical results look sound.

• •

$A=5$ and $B=10$: This looks like a severe case since we have both large $A$ and large $B$. However, the samples are distributed well around the origin, showing a clear existence of the invariance measure. It shows that the complex Langevin method may work well beyond our theory, which also requires future works.

Another observation worth mentioning is that the singularities introduced by gauge cooling seem not to have a bad effect on the numerical results. The failure in the case $A=1$ and $B=2$ looks more like the result of slow decaying distribution function, since $\mathrm{e}^{\mathrm{i}ks}$ increases faster in the imaginary direction when $k$ is larger, which may yield worse results for the observable $\langle O_{k}\rangle$. However, in the cases $A=1,B=0.2$ and $A=5,B=1$ as shown in Figure 4, we do observe nonsmoothness of distribution functions at the origin, which is clearly the effect of the singular drift.

6. Conclusion

Due to the accessibility of the exact solution for optimal gauge cooling problem in the one-dimensional case, we have presented a clear picture for the stabilizing effect of gauge cooling for CLM. As a summary, we would like to point out the following four effects of gauge cooling:

• •

A large number of redundant degrees of freedom are removed (see Theorem 4);

• •

Some components of the drift velocity no longer take effect (see Section 4.1);

• •

Additional drift velocity toward the non-complexified domain is introduced, causing possible confinement of the samples (see Section 4.2);

• •

Singularities are introduced to the drift velocity (see Section 4.2).

The first three effects clearly help stabilize the CLM, while the outcome of the last one remains unclear. Our analysis clearly emphasizes the importance of a good gauge cooling scheme in the multidimensional case, which is part of our ongoing work. Besides, our numerical results in Section 5.2 suggest that the applicability of CLM with gauge cooling should be far beyond the theory given by Theorem 5, as will also be further studied in our future work.