text:
# Lyapunov Exponent Pairing for a Thermostatted Hard-Sphere Gas under Shear in the Thermodynamic Limit

## Abstract

We demonstrate why for a sheared gas of hard spheres, described by the SLLOD equations with an iso-kinetic Gaussian thermostat in between collisions, deviations of the conjugate pairing rule for the Lyapunov spectrum are to be expected, employing a previous result that for a large number of particles $`N`$, the iso-kinetic Gaussian thermostat is equivalent to a constant friction thermostat, up to $`1/\sqrt{N}`$ fluctuations. We also show that these deviations are at most of the order of the fourth power in the shear rate.

Instituut Lorentz, Universiteit Leiden, Postbus 9506, 2300 RA Leiden, The Netherlands
Chemical Physics Theory Group, Chemistry Department, University of Toronto,
80 St. George St. Toronto, Ontario M5S 3H6, Canada

The SLLOD equations of motion, combined with Lees-Edwards boundary condition, were originally proposed in Refs. , and since then they have been convenient tools to calculate the shear viscosity of gases in the bulk by means of non-equilibrium molecular dynamics simulations for many years. These studies consider systems with a large number of mutually interacting particles that are driven by an external shear rate $`\gamma `$ . In these studies, the iso-kinetic Gaussian thermostat is an artificial way to continuously remove the energy generated inside the system due to the work done on it by the external shear field, such that a non-equilibrium steady state, homogeneous in space, can be reached. The Lyapunov spectrum of such systems is of interest since it has been shown that the shear viscosity can be related to the spectrum , which can be numerically obtained as a function of the shear rate. The analysis of the simulation data  indicated that the sum of the largest and the smallest, the sum of the second largest and the second smallest and so on, were the same. The phenomenon of such pairing of the Lyapunov exponents is known as the Conjugate Pairing Rule, or the CPR. Based on this observation, an attempt to prove an exact CPR was made for arbitrary inter-particle potentials and arbitrary $`\gamma `$ , and later studies and better simulation techniques  indicated that for systems obeying the SLLOD equations of motion, the CPR is not satisfied exactly under these general conditions. However, any conclusive theoretical proof regarding the status of an approximate CPR for systems under SLLOD equations of motion is absent in the literature till now, leaving the problem open for a long time.

The SLLOD equations of motion describe the dynamics of a collection of $`N`$ particles constituting a fluid with a macroscopic velocity field $`𝐮(𝐫)=\gamma y\widehat{𝐱}`$. For particles of unit mass, the equations of motion of the $`i`$-th particle, in terms of its position $`𝐫_i`$ and peculiar momentum $`𝐩_i`$, is given by

$`\dot{𝐫}_i=𝐩_i+\gamma y_i\widehat{𝐱},\dot{𝐩}_i=𝐅_i\gamma p_{iy}\widehat{𝐱}\alpha 𝐩_i,`$ (1)
where $`𝐅_i`$ is the force on the $`i`$-th particle due to the other particles in the system. The value of $`\alpha `$, the coefficient of friction representing the iso-kinetic Gaussian thermostat, is chosen such that the total peculiar kinetic energy of the system, $`_ip_i^2/2`$, is a constant of motion in between collisions. In terms of the positions $`𝐫_i`$ and the laboratory momenta $`𝐯_i`$ of the particles, Eq. (1) reads

$`\dot{𝐫}_i=𝐯_i,\dot{𝐯}_i=𝐅_i+\alpha \gamma y_i\widehat{𝐱}\alpha 𝐯_i.`$ (2)
In the present context, the gas particles are hard spheres, which for simplicity are assumed to have unit radius. The dynamics of the gas particles consists of an alternating sequence of flight segments and instantaneous binary collisions. During a flight, the dynamics of the gas particles is described by Eqs. (2) with $`𝐅_i=0`$. At an instantaneous collision between the $`i`$-th and the $`j`$-th sphere, the post-collisional positions and laboratory momenta ($`+`$ subscripts) are related to their pre-collisional values ($``$ subscripts) by

$`𝐫_{i+}`$ $`=`$ $`𝐫_i,𝐫_{j+}=𝐫_j,`$ (3)
$`𝐯_{i+}`$ $`=`$ $`𝐯_i\{(𝐯_i𝐯_j)\widehat{𝐧}_{ij}\}\widehat{𝐧}_{ij}\text{and}`$ (4)
$`𝐯_{j+}`$ $`=`$ $`𝐯_j+\{(𝐯_i𝐯_j)\widehat{𝐧}_{ij}\}\widehat{𝐧}_{ij},`$ (5)
while the positions and the velocities of the rest of the spheres remain unchanged. Here, $`\widehat{𝐧}_{ij}`$ is the unit vector along the line joining the center of the $`i`$-th sphere to the $`j`$-th sphere at the instant of collision. Note that because we applied the iso-kinetic Gaussian thermostat only between collisions, the peculiar kinetic energy changes in individual collisions. These changes are random, both in magnitude and sign, due to the randomness of the collision parameters, and hence it is quite likely that the system would reach a steady state, where the average change of peculiar kinetic energy would be zero.

In terms of the $`3N`$-dimensional vectors $`𝐑=(𝐫_1,𝐫_2,\mathrm{},𝐫_N)`$, $`𝐕=(𝐯_1,𝐯_2,\mathrm{},𝐯_N)`$ and $`\widehat{𝐍}_{ij}`$, whose $`l`$-th entry is given by $`\widehat{𝐍}_{ij}^l=(\delta _{l,i}\delta _{l,j})\widehat{𝐧}_{ij}/\sqrt{2}`$ $`(l=1,2,..,N)`$, Eqs. (2-5) can be compacted to

$`\dot{𝐑}=𝐕,\dot{𝐕}=\alpha \gamma \mathrm{𝐂𝐑}\alpha 𝐕`$ (6)
during a flight segment and

$`𝐑_+=𝐑_{},𝐕_+=𝐕_{}\mathrm{\hspace{0.17em}2}(𝐕_{}\widehat{𝐍}_{ij})\widehat{𝐍}_{ij}`$
at a collision between the $`i`$-th and the $`j`$-th sphere. Here, $`𝐂`$ is a $`3N\times 3N`$ matrix with $`N\times N`$ entries, each of which is a $`3\times 3`$ matrix. In terms of the entry index $`(l,m)`$, in the $`xyz`$-basis, $`𝐂_{lm}=𝐜\delta _{lm}`$ ($`l,m=1,2,..,N`$) and

$`𝐜`$ $`=`$ $`\widehat{𝐱}\widehat{𝐲}=\left[\begin{array}{ccc}\mathrm{0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}& & \\ \mathrm{0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}& & \\ \mathrm{0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}& & \end{array}\right].`$
Having described the dynamics of the infinitesimal deviation $`\delta 𝐗=(\delta 𝐑,\delta 𝐕)`$ between two typical trajectories in the $`6N`$-dimensional phase space for a time $`t`$ as

$`\delta 𝐗(t)`$ $`=`$ $`𝐋(t)\delta 𝐗(0),`$ (8)
the Lyapunov exponents for this system are the logarithms of the eigenvalues of the matrix $`𝚲`$, defined by

$`𝚲=\underset{t\mathrm{}}{lim}\left[\stackrel{~}{}L(t)\right]^{1/(2t)},`$
where $`\stackrel{~}{}L(t)=[𝐋(t)]^\text{T}𝐋(t)`$.

It can be shown that the sufficient condition for the CPR to hold exactly for a dynamical system obeying Eq. (8) is the existence of a constant non-singular matrix $`𝐊`$ satisfying $`𝐊^2𝐈`$, such that

$`[𝐋(t)]^\text{T}𝐊𝐋(t)=\mu 𝐊.`$ (9)
Here, $`\mu `$ is a scalar function of $`t`$. If $`𝐋(t)`$ satisfies Eq. (9), then we call $`𝐋(t)`$ to be “generalized $`\mu `$-symplectic”. It is easy to show from Eq. (9) that if $`\stackrel{~}{L}`$ is an eigenvalue of $`\stackrel{~}{}L(t)`$, then so is $`\mu ^2/\stackrel{~}{L}`$; from which the (exact) CPR follows. For the situations where the CPR has been proved to be exact , only the $`\mu `$-symplecticity case of Eq. (9) (i.e., $`𝐊=𝐉`$, where $`𝐉`$ is the usual symplectic matrix) has been exploited. In this context, we note that despite the similarity between the present problem and the one discussed in Ref. , the elaborate formalism developed therein is not applicable here.

A significant simplification can be achieved by noticing that the coefficient of friction $`\alpha `$, in the non-equilibrium steady state, fluctuates with $`1/\sqrt{N}`$ fluctuations around a fixed value $`\alpha _0`$ in the thermodynamic limit . Thus, to calculate the Lyapunov exponents for large $`N`$, to which we confine ourselves henceforth, $`\alpha `$ can be replaced by $`\alpha _0`$ in Eq. (6), except for a beginning transient time. On average, for small $`\gamma `$, $`\alpha \gamma ^2`$ and so is $`\alpha _0`$. Higher order corrections play a role for larger shear rates.

In the following analysis, we first explore the status of the CPR when the coefficient of friction is a constant, $`\alpha _0`$, and then return to the case where the coefficient of friction represents an iso-kinetic Gaussian thermostat. The detailed derivation of the following results is given elsewhere. At present, we focus only on the main points.

Once $`\alpha _0`$ replaces $`\alpha `$ in Eq. (6), we find that in the time evolution of $`\delta 𝐗`$ over a collision-less flight segment between $`t`$ and $`t+\tau `$ is given by

$`\delta 𝐗(t+\tau )=𝐇(\tau )\delta 𝐗(t).`$ (10)
$`𝐇(\tau )`$ can be decomposed into $`3N\times 3N`$ sub-matrices as

$`𝐇(\tau )=\left[\begin{array}{cc}𝐡^{[1]}(\tau )& 𝐡^{[2]}(\tau )\\ 𝐡^{[3]}(\tau )& 𝐡^{[4]}(\tau )\end{array}\right].`$ (13)
Having further decomposed each of the $`𝐡^{[k]}(\tau )`$ matrices ($`k=1,\mathrm{},4`$) into $`N\times N`$ entries of $`3\times 3`$ matrices as $`𝐡_{lm}^{[k]}(\tau )`$ ($`l`$ and $`m`$ are counted along the row and the column respectively), we have (with $`𝐈`$ as the identity matrix)

$`𝐡_{lm}^{[1]}(\tau )`$ $`=`$ $`\left\{𝐈+\left[\gamma \tau {\displaystyle \frac{\gamma }{\alpha _0}}(1e^{\alpha _0\tau })\right]𝐜\right\}\delta _{lm},`$ (14)
$`𝐡_{lm}^{[2]}(\tau )`$ $`=`$ $`\{{\displaystyle \frac{1e^{\alpha _0\tau }}{\alpha _0}}𝐈`$ (16)
$`+{\displaystyle \frac{\gamma }{\alpha _0^2}}[\alpha _0\tau (1+e^{\alpha _0\tau })2+2e^{\alpha _0\tau }]𝐜\}\delta _{lm},`$
$`𝐡_{lm}^{[3]}(\tau )`$ $`=`$ $`\left\{\gamma [1e^{\alpha _0\tau }]𝐜\right\}\delta _{lm}\text{and}`$ (17)
$`𝐡_{lm}^{[4]}(\tau )`$ $`=`$ $`\left\{e^{\alpha _0\tau }𝐈\gamma \left[\tau +{\displaystyle \frac{1}{\alpha _0}}(1e^{\alpha _0\tau })\right]𝐜\right\}\delta _{lm},`$ (18)
If we now form a $`6N\times 6N`$ matrix $`𝐆`$, \[in the notation of Eq. (13)\] which looks like $`𝐆_{lm}^{[1]}=𝐆_{lm}^{[4]}=Ø`$ and $`𝐆_{lm}^{[2]}=𝐆_{lm}^{[3]}=𝐠\delta _{lm}`$, where

$`𝐠=\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right].`$ (22)
Then $`𝐇(\tau )`$ can be easily shown to satisfy

$`[𝐇(\tau )]^\text{T}𝐆𝐇(\tau )=e^{\alpha _0\tau }𝐆.`$ (23)
Thus, $`𝐇(\tau )`$ is generalized $`\mu `$-symplectic with $`𝐆`$, but it is not $`\mu `$-symplectic, i.e., $`[𝐇(\tau )]^\text{T}𝐉𝐇(\tau )e^{\alpha _0\tau }𝐉`$. The fact that the same analysis \[Eqs. (10-22)\] can be carried out for any constant coefficient of friction (not necessarily $`\alpha _0`$), implies that the CPR is exact for a collision-less gas of point particles obeying Eq. (6) with a constant coefficient of friction. This has been found previously in simulation data .

For the transformation of $`\delta 𝐗`$ over a binary collision between the $`i`$-th and the $`j`$-th sphere, we follow the explicit derivation in Ref. , which in turn is based on the formalism developed simultaneously by Gaspard and Dorfman, and by Dellago and co-workers . The post-collisional infinitesimal deviation vector $`\delta 𝐗_+`$ can be related to its pre-collisional value $`\delta 𝐗_{}`$ by

$`\delta 𝐗_+=𝐌_{ij}\delta 𝐗_{},`$
where the $`6N\times 6N`$ matrix $`𝐌_{ij}`$ can be decomposed into four $`3N\times 3N`$ blocks, having the following structure

$`𝐌_{ij}`$ $`=`$ $`(𝐈\mathrm{\hspace{0.17em}2}\widehat{𝐍}_{ij}\widehat{𝐍}_{ij})\left[\begin{array}{ccc}𝐈& & 0\\ 𝐑& & 𝐈\end{array}\right].`$
Here, $`𝐑`$ is a symmetric matrix. The above expression for $`𝐌`$ implies that the collisions are symplectic, i.e.,

$`𝐌_{ij}^\text{T}𝐉𝐌_{ij}=𝐉,`$
but not generalized symplectic with $`𝐆`$ ($`𝐌_{ij}^\text{T}\mathrm{𝐆𝐌}_{ij}𝐆`$).

We can now express the matrix $`𝐋(t)`$ in terms of the $`𝐇`$ and $`𝐌`$ matrices in the following way: if the dynamics involves free flight segments separated by $`s`$ instantaneous binary collisions at $`t_1`$, $`t_2,\mathrm{},t_s`$ such that $`0<t_1<t_2<\mathrm{}<t_s<t`$, then

$`𝐋(t)=𝐇(\mathrm{\Delta }t_s)𝐌_{i_sj_s}𝐇(\mathrm{\Delta }t_{s1})\mathrm{}𝐌_{i_1j_1}𝐇(\mathrm{\Delta }t_0).`$ (25)
Here, $`\mathrm{\Delta }t_s=tt_s`$ and $`\mathrm{\Delta }t_i=t_{i+1}t_i`$ for $`i=1,\mathrm{},(s1)`$.

The consequences of Eqs. (10-25) can be summarized by the following: for a collection of hard spheres obeying the SLLOD equations of motion with constant coefficient of friction $`\alpha _0`$, (a) the $`𝐇`$ matrices are generalized $`\mu `$-symplectic with $`𝐆`$, but not with $`𝐉`$, and (b) the $`𝐌`$ matrices are symplectic but not generalized $`\mu `$-symplectic with $`𝐆`$. Hence, once the $`𝐇`$ and the $`𝐌`$ matrices are combined together, as in Eq. (25), $`𝐋(t)`$ is seen to be generalized $`\mu `$-symplectic with neither $`𝐆`$ nor $`𝐉`$. This is consistent with the claim that $`𝐋(t)`$ is not generalized $`\mu `$-symplectic (and consequently, the CPR does not hold exactly) for a collection of hard spheres obeying the SLLOD equations of motion with constant coefficient of friction $`\alpha _0`$.

The degree of deviation from an exact CPR must follow from the properties of $`𝐋(t)`$, and to estimate this deviation, we can use either $`𝐊=𝐆`$, or $`𝐊=𝐉`$ in Eq. (9). While the former choice implies that one has to try to estimate the deviation from an exact CPR from the distribution of the unit vectors $`\widehat{𝐍}_{ij}`$’s and the collision angles for different sets of binary collisions in the expression of $`𝐌_{ij}`$’s, the latter choice means that one can make the estimate by using the typical magnitude of a free flight time, i.e., the mean free time $`\tau _0`$. We choose the latter approach, because an estimate of the deviation from the exact CPR can be made at small $`\gamma `$, as a power series expansion in $`\gamma `$. It is important to realize at this point that as the density sets a time scale in the form of the mean flight time $`\tau _0`$ between collisions, the actual dimensionless small parameter corresponding to the shear rate is $`\stackrel{~}{\gamma }=\gamma \tau _0`$.

We begin by constructing another matrix $`𝐇_0(\mathrm{\Delta }t)`$ by setting $`\gamma =0`$ but $`\alpha _00`$ in the explicit form of $`𝐇(\mathrm{\Delta }t)`$ in Eqs. (10-18), i.e.,

$`𝐇_0(\mathrm{\Delta }t)=𝐇(\mathrm{\Delta }t)|_{\alpha _00,\gamma =\mathrm{\hspace{0.17em}0}}.`$
It is easy to show that $`𝐇_0(\mathrm{\Delta }t)`$ satisfies the equation

$`[𝐇_0(\mathrm{\Delta }t)]^\text{T}𝐉𝐇_0(\mathrm{\Delta }t)=e^{\alpha _0\mathrm{\Delta }t}𝐉.`$
Following Eq. (25), we then form the matrix $`𝐋_0(t)`$ as

$`𝐋_0(t)=𝐇_0(\mathrm{\Delta }t_s)𝐌_{i_sj_s}𝐇_0(\mathrm{\Delta }t_{s1})\mathrm{}𝐌_{i_1j_1}𝐇_0(\mathrm{\Delta }t_0),`$ (26)
such that all the $`𝐌_{ij}`$ matrices in Eqs. (25) and (26) are the same. Since both the $`𝐌_{ij}`$ and the $`𝐇_0(\mathrm{\Delta }t)`$ matrices are now $`\mu `$-symplectic with $`𝐉`$, so is $`𝐋_0(t)`$. As a consequence, the logarithms of the eigenvalues of $`\stackrel{~}{}L_0(t)=[𝐋_0(t)]^\text{T}𝐋_0(t)`$ pair exactly. This implies that if we arrange the corresponding Lyapunov spectrum

$`𝚲_0`$ $`=`$ $`\underset{t\mathrm{}}{lim}\left[\stackrel{~}{}L_0(t)\right]^{1/(2t)},`$
in the decreasing order of magnitude as $`\lambda _1^{(0)}\lambda _2^{(0)}\mathrm{}\lambda _{6N}^{(0)}`$, then $`\lambda _i^{(0)}+\lambda _{6Ni+1}^{(0)}=\alpha _0`$.

It is a simple exercise to show that $`𝐇(\mathrm{\Delta }t)𝐇_0(\mathrm{\Delta }t)=O(\stackrel{~}{\gamma }^3)`$, from which we conclude that for $`\mathrm{\Delta }t=\tau =O(\tau _0)`$

$`𝐋(\tau )=𝐋_0(\tau )[𝐈+\stackrel{~}{\gamma }^3𝐁],`$ (27)
where the matrix $`𝐁`$ is of order 1 in $`\stackrel{~}{\gamma }`$ and order 1 in $`N`$. Note that $`𝐁`$ contains higher powers of $`\stackrel{~}{\gamma }`$ as well. Because it involves the matrix $`𝐜`$ and contributions from collisions between spheres, $`𝐁`$ is not proportional to $`𝐈`$, and hence, we cannot regard it simply as a scalar factor (in which case the exact conjugate pairing would be easy to obtain again). Equation (27) implies that for $`\mathrm{\Delta }\stackrel{~}{}L(\tau )\stackrel{~}{}L(\tau )\stackrel{~}{}L_0(\tau )`$

$`\mathrm{\Delta }\stackrel{~}{}L(\tau )`$ $`=`$ $`\stackrel{~}{\gamma }^3[𝐁^\text{T}\stackrel{~}{}L_0(\tau )+\stackrel{~}{}L_0(\tau )𝐁]+\stackrel{~}{\gamma }^6𝐁^\text{T}\stackrel{~}{}L_0(\tau )𝐁.`$ (28)
From Eqs. (27) and (28), we can now see that the differences between $`𝐋(\tau )`$ and $`𝐋_0(\tau )`$, and between $`\stackrel{~}{}L(\tau )`$ and $`\stackrel{~}{}L_0(\tau )`$ are small, by a relative order $`\stackrel{~}{\gamma }^3`$. Therefore the logarithm of the eigenvalues of $`𝐋(\tau )`$ and $`𝐋_0(\tau )`$ also differ by a term of order $`\stackrel{~}{\gamma }^3`$ in an absolute sense. If we now divide the logarithms of these eigenvalues by the time $`\tau `$, we see that the finite time (for time $`\tau `$) Lyapunov exponents, calculated from $`\stackrel{~}{}L_0(\tau )`$ and from $`\stackrel{~}{}L(\tau )`$ (which we denote as $`\lambda _i^{(0)}(\tau )`$ and $`\lambda _i(\tau )`$ respectively, for $`i=1,2\mathrm{}6N`$), differ by a term $`O(\stackrel{~}{\gamma }^3/\tau )=O(\gamma \stackrel{~}{\gamma }^2)`$.

We make one further observation at this stage. The Lyapunov exponents (even the finite time ones) are invariant under $`\gamma \gamma `$, so in a power series expansion in $`\stackrel{~}{\gamma }`$ , the odd powers vanish. Hence, we conclude that the logarithm of the eigenvalues of $`𝐋(\tau )`$ and $`𝐋_0(\tau )`$ must differ by a term of order $`\stackrel{~}{\gamma }^4`$, i.e., the conjugate pairing of $`\lambda _i(\tau )`$’s must be valid up to corrections of the form $`\gamma \stackrel{~}{\gamma }^3`$.

To explicitly extend this formalism to large $`t`$ and thereby obtain a relation between $`\lambda _i`$s and $`\lambda _i^{(0)}`$s, we need to sequentially concatenate a lot of $`𝐋(\tau )`$’s. In general, these matrices neither commute with each other, nor with the $`𝐁`$’s, which prevents us from explicitly demonstrating how the deviation $`[𝐋(t)𝐋_0(t)]`$ is built up. However, we can argue in the following manner: $`\stackrel{~}{}L(t)`$ and $`\stackrel{~}{}L_0(t)`$ are positive definite and symmetric. This allows us to express them in the form $`\stackrel{~}{}L_0(t)=\mathrm{exp}(𝐀_0)`$ and $`\stackrel{~}{}L(t)=\mathrm{exp}(𝐀)`$, where for large $`t`$, both the eigenvalues of $`𝐀_0`$ and $`𝐀`$ must behave $`t`$. From this perspective, the difference between the Lyapunov exponents for $`\stackrel{~}{}L(t)`$ and $`\stackrel{~}{}L_0(t)`$ is related to $`(𝐀𝐀_0)`$. Since the difference between $`\stackrel{~}{}L(t)`$ and $`\stackrel{~}{}L_0(t)`$ has an explicit prefactor of $`\stackrel{~}{\gamma }^3`$, so does $`𝐀𝐀_0`$. Using the symmetry argument that the Lyapunov exponents have to be even functions of $`\gamma `$, we obtain

$`\lambda _i+\lambda _{6Ni+1}=\alpha _0+O(\gamma \stackrel{~}{\gamma }^3).i=1,\mathrm{}6N`$ (29)
For the largest and the most negative Lyapunov exponents, it has been possible to show that they pair to $`\alpha _0`$ plus corrections of $`O(\gamma \stackrel{~}{\gamma }^3)`$ by means of a kinetic theory approach , based on the independence of subsequent collisions of a sphere. Likewise, one expects that in subsequent time-intervals of $`O(\tau _0)`$, the $`𝐋(\tau )`$ matrices are not qualitatively much different from each other. Therefore, we expect that the coefficient of the $`O(\gamma \stackrel{~}{\gamma }^3)`$ term in Eq. (29), to be of the same order as that for a flight time $`\tau =O(\tau _0)`$ \[i.e. of the order of $`𝐁=O(1)`$\], and therefore Eq. (29) to hold.

In summary, for the SLLOD equations with a constant $`\alpha _0`$ thermostat, the finite time Lyapunov exponents obey the CPR up to $`O(\gamma \stackrel{~}{\gamma }^3)`$ when that time is of the order of the mean flight time, and this is expected to hold for the infinite time Lyapunov exponents too. Moreover, the iso-kinetic Gaussian thermostat is equivalent to the constant multiplier thermostat in the thermodynamic limit, and hence one expects that with an iso-kinetic Gaussian thermostat between collisions, the Lyapunov exponent spectrum also exhibits $`O(\gamma \stackrel{~}{\gamma }^3)`$ deviations from the CPR, in the thermodynamic limit. Finally, given that the source of the CPR violation is basically the $`\alpha _0\gamma \mathrm{𝐂𝐑}`$ term in Eq. (6), one can argue that when the gas particles interact with each other by means of a short-ranged, repulsive potential with a constant multiplier thermostat, the violation of the CPR would also be at least of $`O(\stackrel{~}{\gamma }^4)`$ (for gas particles interacting with each other by means of a short-ranged, repulsive potential with an isokinetic Gaussian thermostat, the same results are expected) .

It is a pleasure to thank Prof. J. R. Dorfman and Prof. H. van Beijeren for many useful and motivating discussions regarding this subject. D. P. was supported by the research grants from “Fundamenteel Onderzoek der Materie (FOM)” and that of Prof. J. R. Dorfman, NSF-PHY-9600428. R. v. Z. was supported by the research grant of Prof. H. van Beijeren, by FOM, SMC and NWO Priority Program Non-Linear Systems, and by a grant from the Natural Sciences and Engineering Research Council of Canada.


format: text

text:
# Density Matrix Renormalisation Group Approach to the Massive Schwinger Model

## I INTRODUCTION

The Schwinger model, or quantum electrodynamics in one space and one time dimension, exhibits many analogies with QCD, including confinement, chiral symmetry breaking, charge shielding, and a topological $`\theta `$-vacuum. It is a common test bed for the trial of new techniques for the study of QCD: for instance, several authors have recently discussed new methods of treating lattice fermions using the Schwinger model as an example.

Our purpose in this paper is twofold. First, we aim to explore the physics of this model when an external ‘background’ electric field is applied, as discussed long ago in a beautiful paper by Coleman. Secondly, we wish to demonstrate the application of density matrix renormalisation group methods  to a model of this sort, with long-range, non-local Coulomb interactions.

Coleman showed that the physics of the Schwinger model is periodic in $`\theta =2\pi F/g`$, where $`F`$ is the applied ‘background’ electric field, and $`g`$ is the elementary charge. In the special case $`\theta =\pi `$, some amusing phenomena occur. Whereas the “quarks” in this model are generally confined by a classical linear potential, at $`\theta =\pi `$ it is possible for single quarks to appear deconfined, i.e. move freely, provided they do not cross another quark (“half-asymptotic particles”). Coleman also demonstrated that for $`\theta =\pi `$ a phase transition must occur at some finite value of $`m/g`$, where $`m`$ is the quark mass. These arguments are reviewed more fully in Section II.

There have been very few attempts to verify these predictions numerically, that we are aware of. Hamer, Kogut, Crewther and Mazzolini used finite-lattice techniques to address the problem. They calculated the ground-state energy and ‘string tension’ as functions of $`\theta `$ and $`m/g`$. They located the phase transition at $`\theta =\pi `$ to lie at $`(m/g)_c=0.325(20)`$, with a correlation length index $`\nu =0.9(1)`$. They also attempted to estimate the chiral order parameter; but it was later pointed out  that the chiral order parameter actually suffers from a logarithmic divergence at finite $`m/g`$ in this model. Schiller and Ranft  used Monte Carlo techniques to locate the phase transition at $`(m/g)_c=0.31(1)`$.

Many different numerical methods have been applied to the Schwinger model in zero background field, including strong-coupling series , finite-lattice calculations, Monte Carlo calculations, discrete light-cone quantization and the related light-front Tamm-Dancoff and “fast-moving frame” techniques, a recently proposed “contractor renormalisation group” method, and finally a coupled-cluster expansion. For a recent review, see Sriganesh et al. . Analytic calculations of the mass spectrum have been carried out using mass perturbation theory for small fermion masses, and weak-coupling expansions for large fermion masses.

In this paper, we apply a new technique, namely the density matrix renormalisation group (DMRG), which has been used with great success for lattice spin models and lattice electron models such as the Hubbard model. It was also recently applied to a simple one-particle potential model with asymptotic freedom by Martin-Delgado and Sierra. It might be questioned whether the DMRG can successfully handle a model involving long-range Coulomb interactions such as the Schwinger model. We perform some calculations in zero background field, to show that in fact the approach works extremely well. It gives accurate results for lattices of up to 256 sites, and provides estimates of the continuum limit which are around 50 times more accurate than previous calculations.

The layout of the paper is as follows. In Section II we review Coleman’s discussion  of the Schwinger model in a background field, and the lattice formulation thereof. In Section III the DMRG technique is outlined. In Section IV our main results at background field $`\theta =\pi `$ are presented. Mass gaps and order parameters are calculated, and the critical parameters at the phase transition are accurately estimated, showing that the transition lies in the universality class of the transverse Ising model in (1+1) dimensions. The ‘half-asymptotic’ particles behave in exactly the same way as ‘kinks’ or ‘spinons’ in the transverse Ising model, giving rise to some surprising effects. Finally, in Section V some results at $`\theta =0`$ are presented, and compared with previous results, both analytical and numerical. Our conclusions are summarized in Section VI.

## II THE SCHWINGER MODEL IN A BACKGROUND FIELD

### A Fermion formulation

The Schwinger model Lagrangian density, in standard notation, is

$$=\frac{1}{4}F_{\mu \nu }F^{\mu \nu }+\overline{\psi }(i\overline{)}g\overline{)}Am)\psi $$
(1)
where $`\psi `$ is a 2-component spinor field, since there is no spin in one space dimension. The coupling $`g`$ has dimensions of mass, so the theory is super-renormalisable. Using $`g`$ as the scale of energy, the physical properties of the model are then functions of the dimensionless ratio $`m/g`$. The field strength term is

$$F_{\mu \nu }=_\mu A_\nu _\nu A_\mu .$$
(2)
The equations of motion are the Dirac equation

$$(i\overline{)}g\overline{)}Am)\psi =0$$
(3)
and Maxwell’s equations

$$_\mu F^{\mu \nu }=g\overline{\psi }\gamma ^\nu \psi .$$
(4)
Choosing a time-like axial gauge

$$A_0=0$$
(5)
the field-strength tensor reduces to

$$F^{10}=\dot{A}^1=E$$
(6)
where $`E`$ is the 1-component electric field. Gauss’ law becomes

$$_1E=_1\dot{A}^1=gj^0=g\overline{\psi }\gamma ^0\psi .$$
(7)
The Hamiltonian becomes

$$H=𝑑x\left(i\overline{\psi }\gamma ^1(_1+igA_1)\psi +m\overline{\psi }\psi +\frac{1}{2}E^2\right).$$
(8)
Gauss’ law can be integrated to give

$$E=g𝑑xj^0(x)+F,$$
(9)
showing that $`E`$ is not an independent field, but can be determined in terms of the charge density $`j^0(x)`$, up to the constant of integration $`F`$, which corresponds to a “background field”, as discussed by Coleman . We can think of the background field as created by condenser plates at either end of our one-dimensional universe.

If $`|F|>g/2`$, charged pairs will be produced, and separate to infinity, until the field is reduced within the range $`|F|g/2`$, thus lowering the electrostatic energy per unit length (Fig. 1). Thus physics is periodic in F with period $`g`$, and it is convenient to define an angle $`\theta `$ by

$$\theta =2\pi \frac{F}{g}.$$
(10)
Then we can always choose $`\theta `$ to lie in the interval $`[0,2\pi ]`$.

In the weak-coupling limit $`m/g\mathrm{}`$, the vacuum contains no fermionic excitations, and the vacuum energy density $`ϵ_0`$ is given purely by the electrostatic energy term (we ignore, or renormalise to zero, the energy of the Dirac sea). Hence

$$ϵ_0=\{\begin{array}{cc}\frac{1}{2}F^2=g^2\theta /(8\pi ^2)\hfill & (\theta \pi )\hfill \\ \frac{1}{2}(gF)^2=g^2(2\pi \theta )^2/(8\pi ^2)\hfill & (\pi \theta 2\pi ).\hfill \end{array}$$
(11)
Thus there is a discontinuity in the slope of the energy density, corresponding to a first-order phase transition, at $`\theta =\pi `$. In the strong-coupling limit $`m/g=0`$, on the other hand, chiral invariance demands that the vacuum energy density remains constant as a function of $`\theta `$ (see Sec. II B). Thus we expect a first-order transition at $`\theta =\pi `$ for large $`m/g`$, which terminates at a second order critical point at some finite $`(m/g)_c`$, as illustrated in Fig. 2. This behavior was demonstrated numerically some time ago by Hamer et al. , who located the critical point at $`(m/g)_c=0.325(20)`$, with an associated critical index $`\nu =0.9(1)`$.

Normally, charge is confined in the model: there is a ‘string’ of constant electric field (or flux) connecting any pair of opposite charges . But Coleman  points out that in the very special case $`\theta =\pi `$, or $`F=g/2`$, the peculiar phenomenon of “half-asymptotic” particles arises. In the weak-coupling limit, one can envisage the state shown in Fig. 3. The electric field energy density is the same in between each pair of particles, and they can therefore move freely, as long as they maintain the same ordering, i.e. no pair of fermions interchanges positions.

### B Boson formulation

The one-dimensional fermionic theory can be mapped into an equivalent Bose form . Some of the relevant mappings are

$`:\overline{\psi }\psi :`$ $``$ $`cMN_M\mathrm{cos}(2\sqrt{\pi }\varphi )`$ (12)
$`:i\overline{\psi }\gamma _5\psi :`$ $``$ $`cMN_M\mathrm{sin}(2\sqrt{\pi }\varphi )`$ (13)
$`j^\mu =:\overline{\psi }\gamma ^\mu \psi :`$ $``$ $`{\displaystyle \frac{1}{\sqrt{\pi }}}ϵ_{\mu \nu }^\mu \varphi `$ (14)
The Hamiltonian density in the charge zero sector can then be written

$$=N_M\left[\frac{1}{2}\mathrm{\Pi }^2+\frac{1}{2}(_1\varphi )^2+\frac{1}{2}M^2\varphi ^2cmM\mathrm{cos}(2\sqrt{\pi }\varphi \theta )\right]$$
(15)
Here $`N_M`$ denotes normal ordering with respect to mass $`M`$, $`\varphi `$ is the Bose field and $`\mathrm{\Pi }`$ is its conjugate momentum, and

$$M^2=\frac{g^2}{\pi },$$
(16)
while

$$c=\frac{e^\gamma }{2\pi }$$
(17)
and $`\gamma =0.5774`$ is Euler’s constant. Note that

$$_1\varphi =\sqrt{\pi }j^0$$
(18)
is Gauss’ law, which identifies $`\varphi `$ is proportional to the electric field. Finally, $`\theta `$ is the background field variable, as before. This form of the theory is most convenient for discussing the strong-coupling limit, $`m/g0`$.

A nice discussion of the connection between $`\theta `$ and chiral invariance is given by Creutz . Consider the fermion mass term

$$_m=m\overline{\psi }\psi .$$
(19)
If we now consider a chiral rotation

$$\psi e^{i\theta \gamma _5/2}\psi $$
(20)
the mass term becomes

$$_m=m^{}\overline{\psi }\psi +m_5i\overline{\psi }\gamma _5\psi $$
(21)
where

$`m^{}`$ $`=`$ $`m\mathrm{cos}\theta `$ (22)
$`m_5`$ $`=`$ $`m\mathrm{sin}\theta ,`$ (23)
while the remaining terms of the Lagrangian density are naively invariant. In the Bose form of the theory, the transformed mass term maps into

$$_mcmMN_M\mathrm{cos}(2\sqrt{\pi }\varphi \theta )$$
(24)
which is precisely the form (15), with the chiral rotation parameter $`\theta `$ playing the role of the background field.

It can be seen immediately from (15) that at $`m/g=0`$, the Hamiltonian is independent of the background field variable $`\theta `$ (i.e. is chiral invariant), and reduces to a theory of free, massive bosons, with mass $`M=g/\sqrt{\pi }`$ . There is no sign of the “half-asymptotic” particles found in the weak-coupling limit at $`\theta =\pi `$.

Coleman  gives a very neat semi-classical argument to show how the half-asymptotic particles arise in the Bose formulation. At $`\theta =\pi `$, the Hamiltonian corresponds to an effective potential

$$U(\varphi )=\frac{1}{2}M^2\varphi ^2+cmM\mathrm{cos}(2\sqrt{\pi }\varphi ).$$
(25)
For $`m/g`$ small, there is a unique vacuum at $`\varphi =0`$. For $`m/g`$ large, however, there are two vacua, located at $`\varphi =\pm \frac{1}{2}\sqrt{\pi }`$, and the symmetry $`\varphi \varphi `$ suffers spontaneous breakdown. The two vacua correspond to background field values $`\theta =0`$ or $`2\pi `$, or $`F=0`$ or $`g`$. Creutz notes that in the broken symmetry phase the expectation value of $`\varphi `$ (and hence $`\mathrm{sin}(2\sqrt{\pi }\varphi )`$) should be non-zero, and therefore $`\overline{\psi }\gamma ^5\psi `$ should make a suitable order parameter, by Eq. (14).

Spontaneous breakdown of a $`Z_2`$ symmetry implies that the critical point should belong to the universality class of the (1+1)D or 2D Ising model, with critical indices $`\nu =1`$ and $`\beta =1/8`$, which is consistent with the value for $`\nu `$ found by Hamer et al. .

For a single scalar field in two-dimensional space-time undergoing spontaneous symmetry breaking, there will exist time-independent finite-energy soliton solutions of the classical field equations passing monotonically from one ground state to the other. We may designate the increasing solution a “kink” and the decreasing solution an “antikink” (see Fig. 4). These are the half-asymptotic particles, in Bose language. Kinks and antikinks must alternate with each other when well separated, just as the half-asymptotic fermions do.

At strong couplings, the mass of the single boson excitation at $`\theta =\pi `$ can be estimated semiclassically by fitting a harmonic oscillator to the effective potential at $`\varphi =0`$. Hence one finds

$`{\displaystyle \frac{\mathrm{\Delta }_2}{g}}`$ $`=`$ $`{\displaystyle \frac{M}{g}}\left[1\sqrt{\pi }e^\gamma \left({\displaystyle \frac{m}{g}}\right)+O\left({\displaystyle \frac{m}{g}}\right)^2\right]`$ (26)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{\pi }}}e^\gamma \left({\displaystyle \frac{m}{g}}\right)+O\left({\displaystyle \frac{m}{g}}\right)^2`$ (27)
$`=`$ $`0.5641.78\left({\displaystyle \frac{m}{g}}\right)+O\left({\displaystyle \frac{m}{g}}\right)^2,`$ (28)
where we have denoted this mass by $`\mathrm{\Delta }_2`$, standing for the mass gap in the ‘2-particle’, or electron-positron, sector. A crude linear extrapolation in $`m/g`$ would give the boson mass vanishing at $`m/g=0.317`$, quite close to the critical point found by Hamer et al. .

### C Lattice formulation

We employ the Kogut-Susskind  Hamiltonian spatial lattice formulation of the Schwinger model, with the fermions sited on a “staggered” spatial lattice. Let the lattice spacing be $`a`$, and label the sites of the one-dimensional chain with an integer $`n`$. Define a single-component fermion field $`\varphi (n)`$ at each site $`n`$, and a link variable

$$U(n,n+1)=e^{i\theta (n)}=e^{iagA^1(n)}$$
(29)
on each link***We trust the index $`(n)`$ will differentiate the lattice fields $`\varphi (n)`$ and $`\theta (n)`$ from the quite different fields $`\varphi `$ and $`\theta `$ of Sec. II B. Then the lattice Hamiltonian equivalent to (8) is

$$H=\frac{i}{2a}\underset{n=1}{\overset{N}{}}\left[\varphi ^{}(n)e^{i\theta (n)}\varphi (n+1)\text{h.c.}\right]+m\underset{n=1}{\overset{N}{}}(1)^n\varphi ^{}(n)\varphi (n)+\frac{g^2a}{2}\underset{n=1}{\overset{N}{}}L^2(n)$$
(30)
where the number of lattice sites $`N`$ is even, and the correspondence between lattice and continuum fields is

$$\varphi (n)/\sqrt{a}\{\begin{array}{cc}\psi _{\text{upper}}(x)\hfill & n\text{ even}\hfill \\ \psi _{\text{lower}}(x)\hfill & n\text{ odd}\hfill \end{array}$$
(31)
$`{\displaystyle \frac{1}{ag}}\theta (n)`$ $``$ $`A^1(x)`$ (32)
$`gL(n)`$ $``$ $`E(x).`$ (33)
The $`\gamma `$ matrices are represented by

$$\gamma ^0=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),\gamma ^1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right).$$
(34)
We use a “compact” formulation where the gauge field becomes an angular variable $`0\theta (n)2\pi `$ on the lattice, and $`L(n)`$ is the conjugate spin variable

$$[\theta (n),L(m)]=i\delta _{nm}$$
(35)
so that $`L(n)`$ has integer eigenvalues $`L(n)=0,\pm 1,\pm 2,\mathrm{}`$. In the naive continuum limit $`a0`$, the lattice Hamiltonian (30) reduces to the continuum expression (8).

The Hamiltonian is transcribed to a dimensionless operator

$$W=\frac{2}{ag^2}H=W_0+xV$$
(36)
where

$`W_0`$ $`=`$ $`{\displaystyle \underset{n}{}}L^2(n)+\mu {\displaystyle \underset{n}{}}(1)^n\varphi ^{}(n)\varphi (n),`$ (37)
$`V`$ $`=`$ $`i{\displaystyle \underset{n}{}}\left[\varphi ^{}(n)e^{i\theta (n)}\varphi (n+1)\text{h.c.}\right]`$ (38)
$$\mu =\frac{2m}{g^2a},x=\frac{1}{g^2a^2}.$$
(39)
In the lattice strong-coupling limit $`x0`$, the unperturbed ground state $`|0`$ has

$$L(n)=0,\varphi ^{}(n)\varphi (n)=\frac{1}{2}[1(1)^n],\text{all }n$$
(40)
whose energy we normalise to zero. The lattice version of Gauss’ law is then taken as

$$L(n)L(n1)=\varphi ^{}(n)\varphi (n)\frac{1}{2}[1(1)^n]$$
(41)
which means excitations on odd and even sites create $`1`$ units of flux, corresponding to “electron” and “positron” excitations respectively. Eq. (41) determines the electric field $`L(n)`$ entirely, up to an arbitrary additive constant $`\alpha `$, which then represents the background field. Allowing $`\alpha `$ to be non-zero, the electrostatic energy term is modified to

$$\underset{n}{}L^2(n)\underset{n}{}(L(n)+\alpha )^2.$$
(42)
The physics of the background field then matches precisely with the continuum discussion. If $`|\alpha |>1/2`$, charged pairs will be produced and separate to infinity to lower the electrostatic energy, and bring $`|\alpha |1/2`$. Physics is then periodic in $`\alpha `$ with period 1, and the background field variable is

$$\theta =2\pi \alpha .$$
(43)
In the weak-coupling limit $`m/g\mathrm{}`$ the vacuum contains no fermion excitations, and the ground state energy at $`x=0`$ is

$$\omega _0=N\alpha ^2(|\alpha |<1/2)$$
(44)
corresponding to the ‘zero-loop’ state with $`\{L(n)=0,`$ all $`n\}`$; or

$$\omega _0=N(1\alpha )^2(1/2<\alpha <1)$$
(45)
corresponding to the ‘one-loop’ state with $`\{L(n)=1,\text{ all }n\}`$.

At $`\alpha =1/2`$, these two states are degenerate; and half-asymptotic states consisting of alternating electrons and positrons exist, just as in the continuum model. Let us denote the two states (44) and (45) by $`|\frac{1}{2}`$ and $`|\frac{1}{2}`$, according to the electric fields $`\{L(n)+\alpha =\pm \frac{1}{2},\text{ all }n\}`$ on the links. On a finite lattice at weak coupling, the eigenstates will be the symmetric and antisymmetric combinations

$$|0=\frac{1}{\sqrt{2}}[|\frac{1}{2}+|\frac{1}{2}]$$
(46)
and

$$|0^{}=\frac{1}{\sqrt{2}}[|\frac{1}{2}|\frac{1}{2}].$$
(47)
These states become degenerate in the bulk limit $`N\mathrm{}`$, when spontaneous symmetry-breaking occurs. We shall denote the energy gap between them (or “loop gap”) as $`\mathrm{\Delta }_0`$, standing for the 0-particle sector.

The lowest-energy single-particle excitation (Fig. 14) is

$$|1=\sqrt{\frac{2}{N}}\underset{\begin{array}{cc}n=1& \\ \{n\text{ odd}\}& \end{array}}{\overset{N}{}}\varphi (n)\underset{j=n}{\overset{N}{}}e^{i\theta (j)}|\frac{1}{2}$$
(48)
while the lowest “2-particle” state is

$$|2=\frac{1}{\sqrt{N}}\underset{\begin{array}{cc}n=1& \\ \{n\text{ odd}\}& \end{array}}{\overset{N}{}}\left[\varphi ^{}(n+1)e^{i\theta (n)}\varphi (n)|\frac{1}{2}+\varphi ^{}(n+1)e^{i\theta (n)}\varphi (n+2)|\frac{1}{2}\right].$$
(49)
We shall compute the energy gaps between these states and the ground state, $`\mathrm{\Delta }_1/g`$ and $`\mathrm{\Delta }_2/g`$ respectively.

We will also study two order parameters which can be used to characterize the phase transition at $`\theta =\pi `$. The first one is the average electric field

$$\mathrm{\Gamma }^\alpha =\frac{1}{N}\underset{n}{}(L(n)+\alpha )_0,$$
(50)
which in the weak-coupling limit $`m/g\mathrm{}`$ takes values $`\pm 1/2`$ for the zero-loop and one-loop states respectively. By the Feynman-Hellman theorem, it is proportional to the slope $`E_0/\alpha `$ of the ground-state energy: a glance at Fig. 2 shows that at large $`m/g`$ this slope undergoes a discontinuity at $`\theta =\pi `$, while below $`(m/g)_c`$ it vanishes.

The second order parameter is the axial fermion density suggested by Creutz

$`\mathrm{\Gamma }^5`$ $`=`$ $`i\overline{\psi }\gamma ^5\psi /g_0`$ (51)
$`=`$ $`{\displaystyle \frac{i\sqrt{x}}{N}}{\displaystyle \underset{n}{}}(1)^n\left[\varphi ^{}(n)\varphi (n+1)\text{h.c.}\right]_0`$ (52)
where

$$\gamma ^5=i\gamma ^0\gamma ^1.$$
(53)
Now on a finite lattice there is no spontaneous symmetry breakdown, and the expectation values of the order parameters will remain identically zero. The remedy for this is well-known, however. If the two lowest-energy states which become degenerate in the bulk limit $`N\mathrm{}`$ are denoted $`|0`$ and $`|0^{}`$ respectively, then a simple $`2\times 2`$ matrix calculation shows that the order parameter corresponding to an operator $`Q`$ can be estimated as the overlap matrix element $`0|Q|0^{}`$ on the finite lattice.

In zero background field $`\alpha =0`$, the 2-particle ‘positronium’ excited states of lowest energy in the lattice strong coupling limit are the ‘vector’ state

$$|v=\frac{1}{\sqrt{N}}\underset{n=1}{\overset{N}{}}\left[\varphi ^{}(n)e^{i\theta (n)}\varphi (n+1)+\text{h.c.}\right]|0$$
(54)
and the ‘scalar’ state

$$|s=\frac{1}{\sqrt{N}}\underset{n=1}{\overset{N}{}}\left[\varphi ^{}(n)e^{i\theta (n)}\varphi (n+1)\text{h.c.}\right]|0$$
(55)
where $`|0`$ is the strong-coupling ground state (40).

### D Lattice spin formulation

An equivalent lattice spin formulation can be obtained by a Jordan-Wigner transformation

$`\varphi (n)`$ $`=`$ $`{\displaystyle \underset{l<n}{}}[i\sigma _3(l)]\sigma ^{}(n)`$ (56)
$`\varphi ^{}(n)`$ $`=`$ $`{\displaystyle \underset{l<n}{}}[i\sigma _3(l)]\sigma ^+(n)`$ (57)
giving

$`W_0`$ $`=`$ $`{\displaystyle \underset{n}{}}(L(n)+\alpha )^2+{\displaystyle \frac{\mu }{2}}{\displaystyle \underset{n}{}}(1)^n\sigma _3(n)+{\displaystyle \frac{N\mu }{2}}`$ (58)
$`V`$ $`=`$ $`{\displaystyle \underset{n}{}}\left[\sigma ^+(n)e^{i\theta (n)}\sigma ^{}(n+1)+\text{h.c.}\right]`$ (59)
and

$$L(n)L(n1)=\frac{1}{2}[\sigma _3(n)+(1)^n].$$
(60)
The $`\mathrm{\Gamma }^5`$ order parameter can be written in spin variables

$$\mathrm{\Gamma }^5=\frac{\sqrt{x}}{N}\underset{n}{}(1)^n\left[\sigma ^+(n)e^{i\theta (n)}\sigma ^{}(n+1)+\text{h.c.}\right]_0.$$
(61)
This is the form which we used in the numerical calculations.

## III METHODS

### A The density matrix renormalisation group method

Our results are based on the density matrix renormalisation group (DMRG) method . DMRG has been used primarily to study low dimensional quantum lattice systems in condensed matter physics, and is able to obtain with great accuracy quantities such as the ground and excited state energies, and correlation functions. One of the key features of DMRG is its ability to calculate these quantities for very large system sizes, particularly for systems with a low number of degrees of freedom per site.

The method employed here is the “infinite system” DMRG method, as prescribed by White , used both with open and periodic boundary conditions (OBC and PBC). Due to the presence of the electric field on the links, and the differing nature of the odd and even-numbered sites, some modifications were made to the method, although nothing that changes the spirit of the DMRG. The quantities calculated are the ground state energies, mass gaps and order parameters. We use the form of the Hamiltonian given in (58) and (59).

Let us first look at how the presence of the electric fields affects the implementation of the code. For a particular spin configuration with OBC, if we specify the incoming electric field for the first site of the chain, then according to (60) the electric field for all the links can be deduced. This incoming field can be either $`L_{\text{in}}=0`$ for the zero background field case, or $`L_{\text{in}}=\pm 1/2`$ for the background field case. If PBC are imposed, then as there is no particular link to fix the electric field, we can have loops of electric flux extending throughout the ring. In the presence of a background field we simply add (or subtract) 1/2 from the values of the electric fields due to the spin configuration. A cutoff, or maximum loop value was chosen such that full convergence was reached to machine precision. A loop cutoff of \[-5,5\] was more than sufficient in most cases.

A typical DMRG iteration is shown in Figs. 5 (OBC) and 6 (PBC). The chain (or ring) is split into two blocks and two sites, where blocks contain in general more than one site. In one DMRG iteration, we augment two sites to each block, so that the whole system grows by 4 sites each time. This is more convenient than the standard approach of adding one site at a time, due to the differing nature of the odd and even-numbered sites. For the PBC case this allows one to have the same type of site on each end of the block (i.e. both even or both odd), which considerably simplifies the bookkeeping of the electric fields on the links. In this case both blocks are identical, so that only one augmentation and density matrix is necessary. For the OBC case, we have a different situation. Again augmenting two sites at a time allows us to have the same type of site on each end of the blocks, but in this case the left hand block will always have an odd site on each end, while the right hand block will have an even site on each end. This requires two separate density matrices for each block, as they are not identical.

We start with a 4 site superblock, and continue adding sites until we reach 256 sites. For PBC we keep a maximum of 125 states per spin sector per loop value in a block, which corresponds to a truncation parameter of approximately $`l1100`$, where $`l`$ is the total number of basis states retained in a block. The truncation error for 256 sites can vary from 1 part in $`10^7`$ to $`10^{13}`$ depending on the number of excitations and the parameters. For OBC, again we keep a maximum of 125 states per spin sector, but as there are no electric field loops, here the truncation parameter is $`l300`$. The truncation error for 256 sites is typically 1 part in $`10^{16}`$.

### B Accuracy of the DMRG calculations

The accuracy of the DMRG calculation is strongly dependent on the parameter values $`x`$ and $`m/g`$; in particular we expect the accuracy to worsen near critical regions. We find that OBC provide more accurate results than PBC, which is consistent with previous DMRG studies. We provide details of the accuracy for the worst case scenario, i.e. the parameters closest to the critical regions. This is not representative of the typical accuracies achievable with the DMRG, but gives us an upper bound on the errors.

We examine the point $`(m/g=0.3;x=100;\theta =\pi ;N=256)`$ with PBC, which represents the smallest lattice spacing used for this study, and lies in the critical region. Table I shows the convergence of the DMRG for these parameters as a function of the truncation parameter $`l`$, for the energies and order parameters. The ground state energy can be resolved to 1 part in $`10^6`$, while the two-particle gap is resolved to 1 part in $`10^3`$. The order parameters are subject to round-off errors because they involve the overlap of the two wavefunctions, and hence reduced accuracy is expected: here we have 3 figures. Fig. 7 shows the behaviour of the ground state energy density as a function of the density matrix truncation eigenvalue, which corresponds to the sum of the eigenvalues thrown away from the density matrix. We see a linear relationship between these quantities, which is confirmation that the DMRG is working correctly. Where possible, we obtain the final estimate by taking the results with the two largest truncation parameters, and performing a linear extrapolation to the $`y`$-axis. The error estimate is obtained by taking the difference in the extrapolant and the result with the largest $`l`$.

Using OBC improves the accuracy of the DMRG calculation considerably. Table II shows the convergence of the same quantities with the same parameters as Table I , but with OBC. The ground state energy density is accurate to nearly machine precision, and the 2-particle mass gap to 1 part in $`10^7`$. Again due to round-off errors, we are limited in accuracy for the order parameters, so a similar accuracy to the PBC code is achieved here. A point against using OBC is that the finite-size corrections are much larger for this case. We use OBC to calculate the 1-particle gap and the order parameter estimates in the continuum limit, and PBC otherwise.

## IV RESULTS AT BACKGROUND FIELD $`\theta =\pi `$

### A Analysis of critical behaviour

We use finite-size scaling theory to estimate the position of the critical point, by calculating pseudo-critical points  at each lattice size $`N`$, and each lattice spacing $`x`$. We demand that the scaled energy-gap ratio

$$R_N(m/g)=\frac{N\mathrm{\Delta }_N(m/g)}{(N1)\mathrm{\Delta }_{N1}(m/g)}$$
(62)
is equal to unity at the point $`m/g=(m/g)_N^{}`$, which is the pseudo-critical point. Here $`\mathrm{\Delta }_N(m/g)`$ refers to the energy gap at some finite lattice size $`N`$. In practice we calculate $`R_N(m/g)`$ for a cluster of five points straddling the pseudo-critical point, then use a polynomial interpolation to find $`(m/g)_N^{}`$. The points were chosen at a spacing of $`\mathrm{\Delta }(m/g)=0.02`$ apart, the smallest reasonable spacing which ensures we cover the pseudo-critical point based on the work of Hamer et al. . For this exercise we use the ‘loop gap’ $`\mathrm{\Delta }_0/g`$, which collapses to zero at the pseudo-critical point.

Fig. 8 shows three sample data sets for the pseudo-critical points calculated between $`x=4`$ and $`x=100`$. On a $`1/N^3`$ plot these can be simply extrapolated with a quadratic, with errors in the vicinity of 1 part in $`10^4`$, for all $`x`$. Some numerical instability creeps in for the larger lattice sizes, however, where the change in the gap energy with lattice size $`N`$ becomes so slow that round off errors become appreciable. The continuum limit is then estimated from these bulk critical points, as shown in Fig. 9. A quadratic fit in $`1/\sqrt{x}`$ extracts the continuum limit $`a0`$ or $`x\mathrm{}`$, which we estimate to be

$$\left(\frac{m}{g}\right)_c=0.3335(2).$$
(63)
This is consistent with the previous estimate by Hamer et al. of $`(m/g)_c=0.325(20)`$, or Schiller and Ranft, $`(m/g)_c=0.31(1)`$, but with two orders of magnitude improvement in accuracy.

Finite-size scaling theory  also allows us to estimate the critical indices for the model. The theory tells us that the Callan-Symanzik “beta function”

$$\frac{\beta _N(m/g)}{m/g}=\frac{\mathrm{\Delta }_N(m/g)}{(\mathrm{\Delta }_N(m/g)2(m/g)\mathrm{\Delta }_N^{}(m/g))}$$
(64)
scales at the critical point like $`\beta _N((m/g)_c)N^{1/\nu }`$ as $`N\mathrm{}`$ . Hence we can extract the critical exponent $`\nu `$ from the ratio

$$N\left(1\frac{\beta _N((m/g)_N^{})}{\beta _{N1}((m/g)_{N1}^{})}\right),$$
(65)
since this should approach $`1/\nu `$ as $`N\mathrm{}`$. Alternatively we can use the limiting behavior

$$\frac{\mathrm{ln}[\beta _N((m/g)_N^{})/\beta _{N1}((m/g)_N^{})]}{\mathrm{ln}[N/(N1)]}\frac{1}{\nu }$$
(66)
as $`N\mathrm{}`$. We call the first of these ratios (65) the ‘linear’ estimate, while the second is the ‘logarithmic’ estimate. We calculate these ratios for all lattice sizes $`N`$ and lattice spacings $`x`$. Analysis of the data shows that the ‘logarithmic’ version is more stable numerically, and converges more quickly towards the bulk limit. We show an example of this convergence in Fig. 10 for a particular lattice spacing $`1/\sqrt{x}=ga=0.45`$. We see monotonic convergence towards $`1/\nu 1`$, until the same numerical errors begin to creep in as seen in Fig. 8.

A similar method was used to estimate the critical index $`\beta `$. The order parameters $`\mathrm{\Gamma }_N^\alpha `$ and $`\mathrm{\Gamma }_N^5`$ are expected to scale like $`N^{\beta /\nu }`$ . ‘Linear’ and ‘logarithmic’ ratios are again constructed from these quantities, and analysed in the same way as for the “beta-functions”. Again the logarithmic ratios

$$\frac{\mathrm{ln}[\mathrm{\Gamma }_N((m/g)_N^{})/\mathrm{\Gamma }_{N1}((m/g)_N^{})]}{\mathrm{ln}[N/(N1)]}\frac{\beta }{\nu }$$
(67)
seem to do better than the ‘linear’ ratios, in terms of numerical stability and convergence. Fig. 10 shows the estimate for $`\beta /\nu `$, which uses the electric field order parameter. All final estimates for $`\beta /\nu `$ shown in Table III are calculated using $`\mathrm{\Gamma }_N^\alpha =(L+\alpha )_0`$.

Table III displays our results for the critical exponents, for each different lattice spacing. We see essentially no variation in the exponents with lattice spacing, to within the accuracy of our calculations. Our best estimates for the critical exponents are thus

$`\nu `$ $`=`$ $`0.99(1)`$ (68)
$`\beta /\nu `$ $`=`$ $`0.125(5).`$ (69)
These results provide reasonably conclusive evidence that the Schwinger model transition at $`\theta =\pi `$ lies in the same universality class as the one dimensional transverse Ising model, or equivalently the 2D Ising model, with ($`\nu =1,\beta =1/8`$).

### B Behaviour in the continuum limit

We now turn to estimating continuum limit values for the energy gaps and order parameters. First of all, examination of the data points for 2-particle gaps in Fig. 11 reveal that some extrapolation is necessary, especially for small lattice spacings, to obtain the bulk limit $`N\mathrm{}`$. In order to do so, we need to fit the data with an appropriate function of $`N`$. Now results from the last section revealed that the critical behavior of this model is closely related to that of the transverse Ising model. In that model the finite lattice values away from the critical point converge exponentially to the bulk limit, modulo half-integral powers of $`N`$ (see appendix). Hence we apply a form

$$\mathrm{\Delta }_N=\left(a+\frac{b}{\sqrt{N}}+\frac{c}{N}\right)\mathrm{exp}(dN)+e$$
(70)
to the data, where $`a,b,c,d,e`$ are fit parameters. Fig. 11 shows that this form in fact fits the data very well, for all values of $`m/g`$.

As a double-check of the results of the fit, we also apply a VBS sequence extrapolation routine  to the finite-lattice sequences, which provides an independent estimate of the bulk limit. Errors for the VBS estimates were obtained by examining the columns of the VBS extrapolants, and comparing the results with different VBS parameters $`\alpha `$. Typically the VBS and the fit results agree to better than 1% accuracy, as long as we are not near the critical region where it becomes difficult to extract a reliable estimate.

Fig. 12 shows the extrapolation to the continuum limit $`x\mathrm{}`$, or $`a0`$. We see that generally there is good agreement between the VBS extrapolations and the fit (70), except towards small lattice spacings where a discrepancy opens up. An extrapolation to the continuum limit is now performed by a simple polynomial fit in powers of $`1/\sqrt{x}=ga`$. The double extrapolation therefore introduces quite large errors into our final results, compared to the original DMRG eigenvalues where the errors are at worst of order $`0.1\%`$.

#### 1 Loop energies

Fig. 13 shows our results for the “loop” energy gap $`\mathrm{\Delta }_0/g`$ for all values of $`m/g`$. We see that this gap vanishes at the same point as calculated in (63), and is zero for all $`m/g>(m/g)_c`$, as predicted by Coleman . The convergence to the bulk limit $`N\mathrm{}`$ for this case is much better than for the 2-particle gap illustrated in Fig. 11, such that away from the critical region hardly any extrapolation is needed. This gives us good accuracy for the region $`m/g<0.25`$, where the gaps are resolved to 4 figures. Near the critical region it becomes rather more difficult to extrapolate in a consistent manner.

#### 2 2-particle gap

We plot on the same figure results for the 2-particle gap $`\mathrm{\Delta }_2/g`$. Again, the gap vanishes at the critical point $`(m/g)_c`$, but on either side of the critical point there is a finite gap, and an almost linear behaviour with $`m/g`$. A linear fit through the points in the range $`m/g=[0,(m/g)_c]`$ gives

$$\frac{\mathrm{\Delta }_2}{g}0.5691.72\left(\frac{m}{g}\right),$$
(71)
which agrees very well with the prediction made in (28). Note, however, that the behaviour is not exactly linear.

The $`m/g=0`$ case is unique in the sense that we have a direct check on our results, since an analytic solution is known. For the massless case the background field has no effect on the physics of the model, since any background field is completely screened out. Hence all measurable quantities should be independent of the background field, and in particular the Schwinger boson mass

$$\frac{\mathrm{\Delta }_2}{g}=\frac{1}{\sqrt{\pi }}0.56419.$$
(72)
Our DMRG estimate gives $`\mathrm{\Delta }_2/g=0.57(1)`$, which agrees with this result within errors.

#### 3 1-particle gap

The 1-particle gap must be calculated using open boundary conditions, due to the mismatch in electric fields at either end of the chain. Fig. 14 shows the situation. A single charged particle in the middle will shift the electric field by $`\pm g`$, and hence we have differing electric fields at either end. There is a further complication in the case with open boundaries, in that applying a background field of either $`\alpha =1/2`$ or $`1/2`$ gives a different result for the ground state energy  for a finite lattice. This is due to the ‘staggered lattice’ convention we have adopted, with ‘electrons’ appearing on odd sites, and ‘positrons’ appearing on even sites. The only consistent definition of the one-particle gap for a finite lattice is therefore

$$\mathrm{\Delta }_1/g=E_1\frac{1}{2}(E_0^++E_0^{})$$
(73)
where $`E_1`$ is the 1-particle energy, $`E_0^+`$ is the ground state energy with $`\alpha =1/2`$, and $`E_0^{}`$ that for $`\alpha =1/2`$.

Apart from this complication the procedure in extrapolating to the continuum limit is the same as for the 2-particle state. Our results are shown in Fig. 13. We see that the gap vanishes for $`m/g<(m/g)_c`$, while for $`m/g>(m/g)_c`$ the 1-particle gap is very close to half the 2-particle gap. Once again, the behaviour is very nearly linear in $`m/g`$.

The pattern of eigenvalues exhibited in Fig. 13 bears an extraordinary resemblance to that of the transverse Ising model, even down to the (almost) linear behaviour with $`m/g`$. In particular, we see that the energy of the 1-particle or ‘kink’ state vanishes at the critical point, and then remains degenerate with the ground state for $`(m/g)<(m/g)_c`$. Assuming this degeneracy is exact, this indicates that a ‘kink condensate’ will form in the ground state for small mass, as discussed by Fradkin and Susskind. It also indicates the existence of a ‘dual symmetry’ in the model, additional to that discussed in Section II, which is spontaneously broken in the low-mass region, and has not been explored hitherto. There should also be a ‘dual order parameter’ associated with this symmetry. By analogy with the Ising case, one would expect the ’dual’ order parameter to be the expectation value of the kink creation/destruction operator, and the ’dual’ symmetry to correspond to inversion of this operator; but we have not verified this by explicit computation.

The other notable feature of Fig. 13 is the degeneracy (within errors) between the zero-particle gap $`\mathrm{\Delta }_0/g`$ and the 2-particle gap $`\mathrm{\Delta }_2/g`$ at small mass. Our physical pictures based on a weak-coupling representation do not give clear physical insight into this phenomenon.

All the features just discussed are of course peculiar to the special case $`\theta =\pi `$. For instance, the 1-particle gap at any other value of $`\theta `$ would be infinite.

#### 4 Order Parameters

We may also obtain estimates for the order parameters $`\mathrm{\Gamma }^5`$ and $`\mathrm{\Gamma }^\alpha `$ as functions of $`m/g`$. Our results are displayed in Fig. 15. Both order parameters are zero, within errors, for $`m/g<(m/g)_c`$. Near the critical region, particularly for $`m/g<(m/g)_c`$ it becomes quite difficult to obtain accurate estimates, but it appears that both order parameters turn over and drop abruptly to zero as the critical point is approached from above, consistent with the small exponent $`\beta =1/8`$ found in the previous section. The axial density $`\mathrm{\Gamma }^5`$ decreases steadily towards zero at large $`m/g`$, whereas $`\mathrm{\Gamma }^\alpha `$ approaches the expected value of $`1/2`$. All our results for the gaps $`\mathrm{\Delta }_0/g`$, $`\mathrm{\Delta }_1/g`$, $`\mathrm{\Delta }_2/g`$ and both order parameters are shown in Table IV for future reference.

## V RESULTS AT BACKGROUND FIELD $`\theta =0`$

The case of zero background field has been studied by many authors already, as outlined in the introduction. Our purpose here is to demonstrate how the DMRG can improve on the accuracy of existing results. The main quantities of interest are the “vector” and “scalar” state masses. The most accurate results to date are those of Sriganesh et al. , who used numerical exact diagonalisation results together with a VBS extrapolation to obtain their final estimates. The largest lattice size calculated by these authors was $`N=22`$, hence we expect to be able to do much better using our DMRG algorithm, which can go to much larger lattice sizes.

For the $`\theta =0`$ case, we find that convergence with lattice size is much more rapid than for $`\theta =\pi `$, so much so that with $`N=256`$, there is essentially no extrapolation necessary to obtain the bulk limit. Fig. 16 shows the data for $`m/g=0`$, which is the exactly soluble case. For a particular lattice spacing $`x`$, we have very good convergence with lattice size, so that for $`x=4`$ through to $`x=100`$ we have 6 figure convergence. Even for the largest lattice spacing $`x=400`$ used here, the gaps have converged to 4 figures. Our final continuum estimate is obtained by the method of linear, quadratic and cubic extrapolants in $`1/\sqrt{x}`$ as used in Ref. . An example is shown in Fig. 17. Our final estimate of $`m^{}/g=0.56419(4)`$ for this case agrees extremely well with the analytic result (72).

Table V summarises our results, which show between one and two orders of magnitude improvement in accuracy over the previous best estimates for small values of $`m/g`$. We must note, however, that for values of $`m/g>1`$ we have little or no improvement in accuracy over previous results. To explain this, first we must note that the structure of the eigenvalue function in $`x`$ shifts towards large $`x`$ for large $`m/g`$. But at large $`x`$ and large $`m/g`$ there are many “intruder” states below the vector state for finite $`N`$, artifacts of the finite lattice corresponding to states with no fermion excitations but with loops of electric flux winding around the entire lattice. This restricts the range of $`x`$ that can be used, and hence the accuracy of the calculation. In a future calculation, a way of eliminating these intruder states must be found.

A comparison of our results with previous works is shown in Fig. 18. In the low mass region, we see excellent agreement between our results and the series expansion around $`m/g=0`$. Our results are also fairly consistent with the results of Sriganesh et al., obtained by exact diagonalisation. The ’fast-moving frame’ results of Kröger and Scheu seem to be consistently a little low in this region. In the large mass region the situation is reversed, as Kröger and Scheu obtain excellent agreement with the non-relativistic expansion, while both our results and those of Sriganesh et al. seem to be slightly high. We attribute his to the problem with intruder states discussed above. The Kröger and Scheu quasi-light-cone approach appears to remain more accurate in this weak-coupling region.

Scalar state mass gaps $`m^+/g`$ may also be calculated using DMRG; here we merely demonstrate its applicability for one mass value $`m/g=1.0`$, where we obtain $`m^+/g2m/g=1.118(1)`$. This is again over an order of magnitude better that the previous best estimate of Sriganesh et al. who obtained for the same quantity 1.12(3).

## VI CONCLUSIONS

In this paper, we have demonstrated the application of the numerical DMRG technique of White and collaborators to a lattice gauge model, namely the massive Schwinger model. The results are in most cases nearly two orders of magnitude more accurate than previous calculations. Where exact diagonalisation can treat lattices up to 22 sites, DMRG gives very accurate results up to 256 sites. The long-range Coulomb interaction present in the model has proved no impediment to the DMRG technique: this is very likely connected with the fact that that the Coulomb interaction is screened, and the effective interactions are short-range.

The most interesting results were obtained for the case of background field $`\theta =\pi `$. We have performed a detailed study confirming the existence of the ‘kinks’ or ‘half-asymptotic’ particles predicted by Coleman, and have shown that there is a phase transition at $`(m/g)_c=0.3335(2)`$ belonging in the universality class of the transverse or (1+1)D Ising model. The pattern of energy eigenvalues near the critical point bears a truly remarkable resemblance to the transverse Ising model, and points to some new physical effects, beyond those discussed by Coleman.

In particular, for $`(m/g)<(m/g)_c`$ and $`\theta =\pi `$, we find that the single-kink state is degenerate with the ground state. This points to the existence of a “kink condensate”  in the ground state. It also points to the existence of a ‘dual symmetry’ and ‘dual order parameter’ in the model, which have not been discussed hitherto. We also find that the gap between the symmetric and antisymmetric ‘loop state’ combinations appears to be exactly degenerate with the 2-particle vector gap in this low-mass region.

Some calculations were also carried out at zero background field, $`\theta =0`$. The vector gap was obtained with great precision for small $`m/g`$; but at large $`m/g`$ the accuracy was spoiled a little by finite-lattice “intruder” states.

The same methods should be applicable to other lattice gauge models in (1+1) dimensions. In higher dimensions, however, the DMRG technique does not have such a large comparative advantage over other techniques. Strenuous efforts are being made to develop improved DMRG algorithms for lattice models in (2+1) D; but nothing has even been attempted in (3+1)D, as far as we are aware.

## VII ACKNOWLEDGMENTS

We would like to thank Professors Michael Creutz and Jaan Oitmaa for useful discussions on this topic, and Dr. Zheng Weihong for numerical checks of some of our results by series methods. We are grateful for computational facilities provided by the New South Wales Center for Parallel Computing, the Australian Center for Advanced Computing and Communications and the Australian Partnership for Advanced Computing. R.B. was supported by the Australian Research Council and the J. G. Russell Foundation.

## VIII APPENDIX

Here we calculate analytically the form of the finite-size scaling corrections in the 1D transverse Ising model. We follow the discussion by Hamer and Barber, hereafter referred to as I.

The Hamiltonian is taken as

$$H=\underset{m=1}{\overset{M}{}}(1\sigma _3(m))x\underset{m=1}{\overset{M}{}}\sigma _1(m)\sigma _1(m+1)$$
(74)
The model can be solved exactly by the methods of Schultz, Mattis and Lieb. The mass gap on a finite lattice of $`M`$ sites with periodic boundary conditions is

$$F(x,M)=2(1x)+2M[T_{2M}(x)T_M(x)]$$
(75)
where

$$T_M(x)=\frac{1}{M}\underset{k=0}{\overset{M1}{}}\mathrm{\Lambda }(\frac{\pi k}{M})$$
(76)
and

$$\mathrm{\Lambda }(\theta )=[(1x)^2+4x\mathrm{sin}^2\theta ]^{1/2}.$$
(77)
For fixed $`x`$, as $`M\mathrm{}`$,

$$\stackrel{~}{T}_M(x)=\frac{1}{\sqrt{x}}T_M(x)c_0+(\frac{4\pi }{x})(1x^2)^{1/2}\frac{e^{2Mw_c(x)}}{(2M)^{3/2}}$$
(78)
where

$$w_c(x)=\mathrm{sinh}^1(\frac{1x}{2\sqrt{x}})$$
(79)
(Eq. (A1.16) of ref. I) and hence

$$F(x,M)2(1x)2\pi \sqrt{\frac{2(1x^2)}{xM}}e^{2Mw_c(x)}$$
(80)
which exhibits the leading finite-size corrections, as required.


text:
# Spontaneous breaking of rotational symmetry in superconductors1footnote 11footnote 1Phys. Rev. Lett. 88, 252503 (2002)

## Abstract

We show that homogeneous superconductors with broken spin/isospin symmetry lower their energy via a transition to a novel superconducting state where the Fermi-surfaces are deformed to a quasi-ellipsoidal form at zero total momentum of Cooper pairs. In this state, the gain in the condensation energy of the pairs dominates over the loss in the kinetic energy caused by the lowest order (quadrupole) deformation of Fermi-surfaces from the spherically symmetric form. There are two energy minima in general, corresponding to the deformations of the Fermi-spheres into either prolate or oblate forms. The phase transition from spherically symmetric state to the superconducting state with broken rotational symmetry is of the first order.

In the original Bardeen-Cooper-Schrieffer (BCS) theory of bulk superconductivity, the condensate wave-function describes a quantum coherent state which is invariant under spatial and time (particle-hole) reversal transformations BCS . External perturbations which act on discrete quantum variables, like the spin of the fermions, break the particle-hole symmetry. A typical example is a metallic superconductor in a high magnetic field where the Pauli paramagnetism induces an asymmetry in the populations of the spin-up and spin-down electrons. The superconducting state is quenched via a first order phase transition once the splitting in the energy spectrum of spin-up and spin-down electrons becomes of the order of the pairing gap in the unpolarized state METALS . The crossover from the BCS to the normal state can be understood in terms of the phases space overlap between the fermionic states located at the top of their individual Fermi-surfaces. The pairing gap is maximal for the symmetric (unpolarized) state with perfectly overlapping Fermi-surfaces. As these are driven apart by the “polarizing field” the phase-space available for the pairing decreases and the gap is successively suppressed. At finite temperatures the smearing of the Fermi-surfaces increases the phase-space overlap and hence the critical field at which the superconductivity is quenched.

The superconducting state sustains larger asymmetries if the translational symmetry is broken. Larkin and Ovchinnikov and Fulde and Ferrell (LOFF) argued that the crossover from the BCS to the normal state, as the Fermi-surfaces are driven apart by the polarizing field, occurs via a spatially inhomogeneous superconducting phase LARKIN ; FULDE . The Cooper pairs carry a finite total momentum in the LOFF phase, i.e., the centers of the Fermi-spheres are shifted allowing for a partial phase-space overlap.

This paper suggests an alternative mechanism of breaking the symmetry which is based on a deformation of the spherical Fermi-surfaces, to the lowest order, into quasi-ellipsoidal form. If the total momentum of the Cooper pairs is zero, as we shall assume in the following, the deformation spontaneously breaks the rotational symmetry. The novel superconducting phase maintains its stability due to the dominance of the condensation energy of the Cooper pairs over the loss in the kinetic energy of the system caused by the deformation of the Fermi-surfaces. We shall assume that there is a single axis along which the symmetry is broken, although more complicated patterns of symmetry breaking are possible (simultaneous breaking of the rotational and translational symmetries, higher order multipole deformations of Fermi-surfaces, etc). Note that deformed or non-spherical Fermi-surfaces are common for electrons in solids; here we treat systems which are homogeneous in the normal state, i.e., any deformations of the spherical shape of the Fermi-surfaces would correspond to a metastable state.

We start with the BCS-Gorkov equations in energy-momentum representation

$`{\displaystyle \underset{\gamma }{}}\left(\begin{array}{cc}\omega E_{\alpha \gamma }& \mathrm{\Delta }_{\alpha \gamma }\\ \mathrm{\Delta }_{\alpha \gamma }^{}& \omega +E_{\alpha \gamma }\end{array}\right)\left(\begin{array}{cc}G_{\gamma \beta }& F_{\gamma \beta }\\ F_{\gamma \beta }^{}& \overline{G}_{\gamma \beta }\end{array}\right)=\delta _{\alpha \beta }\widehat{\mathrm{𝟏}},`$ (5)
where $`G_{\gamma \beta }(\omega ,𝐩)`$ and $`F_{\gamma \beta }(\omega ,𝐩)`$ refer to the full normal and anomalous retarded propagators, $`\mathrm{\Delta }_{\alpha \gamma }(\omega ,𝐩)`$ is the anomalous self-energy, and the diagonal matrix elements of the first matrix correspond to the inverse of free-single particle propagators; the Greek indeces $`\alpha ,\beta \mathrm{}=1,2`$ label the two different species, and $`\omega `$ and $`𝐩`$ refer to the particle energy and momentum. (Note that the center-of-mass momentum of particles is zero). Suppose that the rotational symmetry is broken by a deformation of the Fermi-surfaces from spherical form. The quasiparticle spectrum of species $`\alpha `$ can be parameterized, then, as

$`E_\alpha ={\displaystyle \frac{p^2}{2m_\alpha }}\mu _\alpha \left(1ϵ_\alpha \mathrm{cos}^2\theta \right),`$ (6)
where $`\mu _\alpha `$ are the chemical potentials of the particles in the undeformed state, $`\theta `$ is the angle between the particle momentum $`p`$ and the axis of symmetry breaking; the deformation of the Fermi-sphere in Eq. (6) is truncated at the lowest order non-trivial axisymmetric deformation, which is described by the $`l=2`$, $`m=0`$ term of the Legendre polynomials associated with an expansion in spherical harmonics. The constant energy surfaces of quasiparticle excitations defined by Eq. (6) represent quasi-ellipsoids of revolution with an ellipticity $`ϵ_\alpha `$. For $`ϵ_\alpha =0`$ the spectrum (6) is the true eigenstate of the unpaired, homogeneous system in the absence of external fields. (Note that the deformation described by Eq. (6) does not need to conserve the volume of the Fermi-sphere, as we impose a self-consistency condition for the total density of the system, see Eq. (13) below).

In the following we shall neglect the possible pairing among the same species ($`\mathrm{\Delta }_{\alpha \alpha }=0`$) so that only the off-diagonal elements of the anomalous self-energy matrix are non-zero. The quasiparticle excitation spectrum in the superconducting phase is determined in the standard fashion by finding the poles of the propagators in Eq. (5):

$`\omega _{1,2}=E_A\pm \sqrt{E_S^2+|\mathrm{\Delta }|^2},`$ (7)
where the symmetric and anti-symmetric parts of the spectrum (which are even and odd with respect to the time-reversal symmetry) are defined as $`E_{S,A}=(E_1\pm E_2)/2`$. The solution of Eq. (5) can be written in terms of the eigenstates (7) as

$`G_{1,2}`$ $`=`$ $`{\displaystyle \frac{u_p^2}{\omega \omega _{1,2}+i\eta }}+{\displaystyle \frac{v_p^2}{\omega \omega _{2,1}+i\eta }},`$ (8)
$`F`$ $`=`$ $`u_pv_p\left({\displaystyle \frac{1}{\omega \omega _1+i\eta }}{\displaystyle \frac{1}{\omega \omega _2+i\eta }}\right),`$ (9)
where the Bogolyubov amplitudes are

$`u_p^2={\displaystyle \frac{1}{2}}+{\displaystyle \frac{E_S}{2\sqrt{E_S^2+|\mathrm{\Delta }|^2}}},v_p^2={\displaystyle \frac{1}{2}}{\displaystyle \frac{E_S}{2\sqrt{E_S^2+|\mathrm{\Delta }|^2}}}.`$ (10)
The mean-field approximation to the anomalous self-energy yields the gap equation

$`\mathrm{\Delta }(𝐩)=2{\displaystyle \frac{d\omega ^{}d𝐩^{}}{(2\pi )^4}V(𝐩,𝐩^{})\mathrm{Im}F(\omega ^{},𝐩^{})f(\omega ^{})},`$ (11)
where $`V(𝐩,𝐩^{})`$ is the bare interaction 222For the sake of simplicity we ignore the renormalization of the pairing interaction in the particle-hole channel and assume a time-local interaction., $`f(\omega )=\left[\mathrm{exp}(\beta \omega )+1\right]^1`$is the Fermi distribution function and $`\beta `$ is the inverse temperature. The $`\omega `$ integration is straightforward in the quasiparticle approximation, since the frequency dependence of the propagator is constrained by the on-shell condition. For $`S`$-wave interactions the potential depends only on the absolute magnitude of the quasiparticle momenta. In this case the gap equation simplifies to

$`\mathrm{\Delta }(p)`$ $`=`$ $`{\displaystyle \frac{p^{}_{}{}^{}2dp^{}}{(2\pi )^2}V(p,p^{})d\mathrm{cos}\theta ^{}}`$ (12)
$`\times `$ $`{\displaystyle \frac{\mathrm{\Delta }(p^{})}{2\sqrt{E_S^2+\mathrm{\Delta }(p^{})^2}}}\left[f(\omega _1)f(\omega _2)\right].`$
Note that the deformation of the Fermi-surfaces enters the gap equation as a parameter which is determined by the minimum of the ground state energy of the superconducting phase. For strongly coupled superconductors the gap equation (12) is supplemented by the normalization condition for the net density $`\rho \rho _1+\rho _2`$ at a constant temperature. The densities of the species are given by

$`\rho _{1,2}`$ $`=`$ $`2{\displaystyle \frac{d^4p}{(2\pi )^4}\mathrm{Im}G_{1,2}(\omega ,𝐩)f(\omega )}`$ (13)
$`=`$ $`{\displaystyle \frac{d^3p}{(2\pi )^3}\left\{u_p^2f(\omega _{1,2})+v_p^2f(\omega _{2,1})\right\}},`$
where the summation is over the discrete quantum variables. The second equality follows in the quasiparticle approximation.

Next we turn to the thermodynamic properties of the superconducting state. At a fixed density and temperature the relevant thermodynamic potential is the free energy:

$`F|_{\rho ,\beta }=U\beta ^1S,`$ (14)
where $`U`$ is the internal energy and $`S`$ is the entropy. In the mean-field approximation the entropy is given by the expression

$`S`$ $`=`$ $`k_B{\displaystyle }{\displaystyle }{\displaystyle \frac{d^3p}{(2\pi )^3}}\{f(\omega _1)\mathrm{ln}f(\omega _1)+\overline{f}(\omega _1)\mathrm{ln}\overline{f}(\omega _1)`$ (15)
$`+`$ $`f(\omega _2)\mathrm{ln}f(\omega _2)+\overline{f}(\omega _2)\mathrm{ln}\overline{f}(\omega _2)\},`$
where $`\overline{f}(\omega _\pm )=[1f(\omega _\pm )]`$, $`k_B`$ is the Boltzmann constant. The internal energy, defined as the grand canonical statistical average of the Hamiltonian, is

$`U`$ $`=`$ $`{\displaystyle }{\displaystyle }{\displaystyle \frac{d^3p}{(2\pi )^3}}\{[n_1(p)E_1(p)+n_2(p)E_2(p)]`$ (16)
$`+`$ $`{\displaystyle }{\displaystyle \frac{d^3p^{}}{(2\pi )^3}}V(p,p^{})\nu (p)\nu (p^{})\},`$
where

$`n_{1,2}(p)`$ $``$ $`u_p^2f(\omega _{1,2})+v_p^2f(\omega _{2,1}),`$ (17)
$`\nu (p)`$ $``$ $`u_pv_p\left[f(\omega _1)f(\omega _2)\right].`$ (18)
The first term in Eq. (16) is the kinetic energy while the second term includes the mean field interaction among the particles in the condensate. The true ground state of the system minimizes the free energy difference $`\delta F|_{\rho ,\beta }`$ between the superconducting and normal states (the free energy in the normal state follows from Eqs. (15) and (16) when $`\mathrm{\Delta }=0`$).

The deformations of the Fermi-spheres can be described in terms of the “conformal deformation” $`ϵ=(ϵ_1+ϵ_2)/2`$ and the “relative deformation” $`\delta ϵ=(ϵ_1ϵ_2)/2`$. Then, the symmetric and anti-symmetric parts of the energy spectrum can be written as

$`E_S`$ $``$ $`{\displaystyle \frac{p^2}{2m}}\mu \left[1+ϵ\mathrm{cos}^2\theta \left(1+{\displaystyle \frac{\delta ϵ\delta \mu }{ϵ\mu }}\right)\right],`$ (19)
$`E_A`$ $``$ $`\delta \mu +\left(\mu \delta ϵ+ϵ\delta \mu \right)\mathrm{cos}^2\theta ,`$ (20)
where $`\mu =(\mu _1+\mu _2)/2`$, $`\delta \mu =(\mu _1\mu _2)/2`$ (here we ignore the difference in the masses of the spin/isospin up and down quasi-particles). Equations (12), (13) and (14) form a closed system, which determines the pairing gap and the ground state energy of a superconductor for constant density asymmetry $`\alpha =(\rho _1\rho _2)/(\rho _1+\rho _2)`$. The values of the deformation parameters $`\delta ϵ`$ and $`ϵ`$ are obtained by requiring that the free-energy attains its minimum. Note that in the weak coupling limit Eqs. (12) and (13) decouple and one may solve for $`\mathrm{\Delta }`$ as a function of $`\delta \mu `$ instead of $`\alpha `$. Apart from the fact that $`\alpha `$, rather than $`\delta \mu `$, is the measurable quantity, there is an additional reason for solving the full set of equations. The gap equation alone is symmetric under exchange $`E_AE_A`$, which implies that the solutions are symmetric under the simultaneous change of the signs of $`\delta \mu `$ and $`\delta ϵ`$. Eq. (13), however, does not have this symmetry and the solutions are distinct under the sign transformation above.

As a specific example, which illustrates the solutions above, we consider isospin-singlet (neutron-proton) pairing in nuclear matter in the $`{}_{}{}^{3}S_{1}^{}`$-$`{}_{}{}^{3}D_{1}^{}`$ channel SD\_PAIRING1 ; SD\_PAIRING2 ; SD\_PAIRING3 . The gap in the isospin symmetric case is $`\mathrm{\Delta }_{00}=12`$ MeV if we use as the bare interaction the Argonne potential and ignore the renormalization of the mass of the particles in the normal state due to interactions. The modification of the particle self-energy in nuclear medium (see for a review MUETHER ) affect the absolute magnitude of the gap and rescale its dependence on the parameters. Clearly, with these approximations, our model is schematic, however we do not expect qualitative changes when renormalization of the interaction and the bare mass are included.

The BCS solutions for the $`n`$-$`p`$ pairing has been studied for the homogeneous (translationally and rotationally) invariant state under isospin asymmetric conditions NUCL\_BCS1 ; NUCL\_BCS2 ; NUCL\_BCS3 ; NUCL\_BCS4  and the inhomogeneous state with broken translational symmetry NUCL\_LOFF ; COMMENT  (the nuclear analog of the LOFF phase. (The flavor asymmetric $`<qq>`$-condensate in high density QCD is another example of strongly coupled superconductor where the breaking of translational symmetry plays a role KRISHNA ). Consistent with the assumption $`\mathrm{\Delta }_{\alpha \alpha }=0`$ above we ignore the $`n`$-$`p`$ pairing in the $`{}_{}{}^{1}S_{0}^{}`$ channel MORTEN .

Fig. 1 shows the pairing gap as a function of the density asymmetry and the relative deformation $`\delta ϵ`$ for vanishing conformal deformation ($`ϵ=0`$). Here and below we assume a density $`\rho =0.16`$ fm-3 and a temperature $`T=3`$ MeV. The gap parameter is normalized to its value $`\mathrm{\Delta }_{00}`$ in the isospin symmetric rotationally/translationally invariant state. Although $`\alpha `$ changes in the interval $`[1;1]`$ in general, the symmetry of the equations with respect to the indexes labeling the species reduces the range of $`\alpha `$ to $`[0;1]`$. Assuming neutron excess implies that the Greek indeces equal 1 refer to neutrons and while those equal 2 refer to protons. The relative deformation obviously is not bounded and can assume both positive and negative values.

For positive values of $`\delta ϵ`$, which imply an oblate deformation for the Fermi-surface of neutrons and a prolate deformation for the Fermi-surface of protons, the solutions for the gap equation show the following features. For arbitrary constant $`\delta ϵ`$ the gap is maximal at $`\alpha =0`$ and is suppressed as the asymmetry is increased. For constant $`\alpha `$, $`\mathrm{\Delta }/\delta ϵ=0`$ corresponds to a maximum at $`\delta ϵ0`$ in the large $`\alpha `$ limit. The position of the maximum is independent of $`\alpha `$ and is located around $`\delta ϵ=0.5`$ in our model; this value also corresponds to the critical asymmetry $`\alpha _c`$ at which the superconducting state vanishes. Note that for $`\alpha `$ around $`\alpha _c`$ the gap exists only in the deformed state. For $`\alpha =0`$, Eqs. (12)-(14) are symmetric under exchange of the sign of deformation and so is the gap function. In particular, for $`\alpha =0`$, the critical deformation for positive and negative deformations coincide. For finite $`\alpha `$ the dependence of the gap on the relative deformation depends on the sign of $`\delta ϵ`$. In contrast to the positive range of $`\delta ϵ`$, where the maximum value of the gap is attained at constant $`\delta ϵ`$, for negative $`\delta ϵ`$ the maximum increases as a function of the deformation and saturates around $`\delta ϵ1`$. Quite generally, to maintain the maximal phase space overlap, the system prefers to keep the sign of $`\delta \mu `$ opposite to that of $`\delta ϵ`$. Fig. 2 shows a 2-dimensional projection of a configuration of deformed Fermi-surfaces for $`\delta ϵ=0.64`$ which minimizes the free-energy for fixed $`\alpha =0.35`$.

The difference in the free-energies of the superconducting and normal states $`\delta F`$ is shown in Fig. 3. Owing to the symmetries of the underlying equations, $`\delta F`$ is symmetric with respect to the sign change of $`\delta ϵ`$ when $`\alpha =0`$. For finite $`\alpha `$’s, a minor departure from the rotational invariant state leads to a decrease in $`\delta F`$ which develops two minima for either sign of $`\delta ϵ`$. This behavior can be traced back to the increase of the potential energy with increasing gap (cf. Fig. 1). Note that, although there are non-trivial solutions to the gap equation in the large $`\alpha `$ limit for $`\delta ϵ1`$, these solutions do not lower the energy of the system. The superconducting phase becomes unstable for $`\alpha >0.4`$ due to the increase in the kinetic energy caused by the large deformation of the Fermi-surfaces, so that $`\delta F`$ is a nearly even function of $`\delta ϵ`$. An inspection of the latent heat associated with the phase transition at finite temperatures shows that this quantity does not vanish at the crossover from the spherically symmetric to the deformed superconducting state. Consequently, the phase transition associated with the breaking of the rotational symmetry is of the first order.

To summarize, this paper suggests a novel mechanism of symmetry breaking in superconducting system with particle-hole asymmetry. The lowest order (quadrupole) deformation of the Fermi-surfaces (at zero total momentum of the Cooper pairs) increases the phase-space overlap between the Fermi-surfaces of paired quasiparticles, which is otherwise depleted by the asymmetry in the particle/hole populations. As a result, the free-energy develops minima for finite deformations, since the gain in the (negative) pairing potential energy dominates the increase in the kinetic energy caused by the deformation. Since the deformed ground state spontaneously breaks the rotational symmetry, the dynamic properties of the superconducting state with deformed Fermi-surfaces such as the sound attenuation, the infrared absorption or the Meissner effect will be anisotropic. The results above do not depend on the nature of fields inducing the asymmetry in the fermion populations nor on the nature of the pairing forces and should be applicable to a wide range of fermionic systems.