text:
# Almost any quantum spin system with short-range interactions can support toric codes

## 1 Introduction

Error correction is a key ingredient in any good recipe for quantum information processing (see  for an overview and further references). Many ingenious schemes have been invented to that effect. A particularly interesting approach has been suggested by Kitaev and colleagues in a series of beautiful papers , namely the possibility of implementing quantum error correction on the physical level.

The idea is quite simple. In , Kitaev considers a specific Hamiltonian $`H_\mathrm{\Lambda }`$ for a quantum spin-1/2 system on a finite lattice $`\mathrm{\Lambda }^2`$ with periodic boundary conditions (the spins are taken to reside on the bonds of $`\mathrm{\Lambda }`$). For any finite $`\mathrm{\Lambda }`$, the ground state of $`H_\mathrm{\Lambda }`$ is fourfold degenerate; information is stored in the corresponding four-dimensional eigenspace. Addition of a small local perturbation of a certain kind modifies the Hamiltonian to $`H_\mathrm{\Lambda }(ϵ)`$, where $`ϵ`$ is the perturbation strength. The effect of the perturbation is to introduce an energy splitting between the degenerate ground-state levels of the unperturbed Hamiltonian. Kitaev then argues that there exists a constant $`ϵ_0`$ so that, at low temperatures and for all $`\left|ϵ\right|ϵ_0`$, the energy splitting for sufficiently large $`\mathrm{\Lambda }`$ is given by $`\mathrm{exp}(c\sqrt{N_\mathrm{\Lambda }})`$, where $`N_\mathrm{\Lambda }`$ is the number of bonds in $`\mathrm{\Lambda }`$. In other words, in the thermodynamic limit $`\mathrm{\Lambda }^2`$, the ground state is still fourfold degenerate, and any sufficiently weak perturbation is ”washed out” by the system itself.

The four-body interactions comprising Kitaev’s Hamiltonian were originally considered by him in Ref.  as generators for a family of stabilizer codes , which he termed ”toric codes.” The remarkable feature of toric codes is the fact that, despite their apparent nonoptimality (in the sense of ), they require only local operations for their implementation and can correct any number of errors (by making the lattice large enough). The bulk of Kitaev’s analysis of toric codes was concerned with their properties as ”conventional” quantum error-correcting codes  that require active intervention through frequent measurements and other external processing. The issue of constructing ”self-correcting” quantum spin systems on the basis of toric codes has been taken up again only very recently by Dennis, Kitaev, Landahl, and Preskill . Their approach, however, is centered around topological features of toric codes and delves deep into such subjects as nonabelian gauge theory .

On the other hand, the very idea of physical error correction is so tantalizing, both practically and conceptually, that one cannot help but wonder: how generic are phenomena of this kind? In this Letter I show that a few results in statistical mechanics of quantum spin systems point towards the conclusion that physical error correction is fairly common in such systems, under quite reasonable conditions. These preliminary findings are, in my opinion, suggestive and interesting enough, and might spur further research into this topic both in the quantum information community and in the statistical mechanics community. A detailed analysis is forthcoming .

## 2 Laying out the ingredients

First of all, let us agree on the ingredients necessary for analysis of a self-correcting quantum spin system. Let $`\mathrm{\Lambda }^\nu `$ be a finite lattice. We assume here that $`\nu 2`$ is an integer. Let $`_0`$ be the $`(2S+1)`$-dimensional Hilbert space of a single particle of spin $`S`$. Spins are situated on lattice sites $`\mathrm{}\mathrm{\Lambda }`$ (in Kitaev’s construction, spins were located on lattice bonds). In order to retain a superficial analogy with stabilizer codes, we will assume that the unperturbed Hamiltonian $`H_\mathrm{\Lambda }`$ is classical, i.e., the interactions comprising it generate an abelian subalgebra of the algebra $`(_\mathrm{\Lambda })`$ of all linear operators on $`_\mathrm{\Lambda }:=_\mathrm{}\mathrm{\Lambda }_{\mathrm{}}`$, where $`_{\mathrm{}}`$ is an isomorphic copy of $`_0`$. That is,

$$H_\mathrm{\Lambda }=\underset{M\mathrm{\Lambda }}{}\mathrm{\Phi }_M,$$
(1)
where each $`\mathrm{\Phi }_M`$ is a Hermitian operator on $`_M:=_\mathrm{}M_{\mathrm{}}`$, and $`[\mathrm{\Phi }_M,\mathrm{\Phi }_N]=0`$. We assume periodic boundary conditions (that is, the lattice $`\mathrm{\Lambda }`$ is drawn on the torus $`^\nu /`$). We let $`\{|\underset{¯}{\sigma }\}`$ be the orthonormal basis of $`_\mathrm{\Lambda }`$ in which $`H_\mathrm{\Lambda }`$ is diagonal; the basis vectors are labelled by classical spin configurations, $`\underset{¯}{\sigma }=\{\sigma (\mathrm{})\}_\mathrm{}\mathrm{\Lambda }`$ with $`\sigma (\mathrm{})\{S,S+1,\mathrm{},S1,S\}`$. We also assume that the smallest eigenvalue of $`H_\mathrm{\Lambda }`$ is equal to zero, and that its geometric multiplicity is $`m2`$. We denote the corresponding eigenspace by $`_\mathrm{\Lambda }^g`$.

The effect of errors is modelled by introducing an off-diagonal perturbation term to the Hamiltonian:

$$H_\mathrm{\Lambda }(ϵ):=H_\mathrm{\Lambda }+ϵP,$$
(2)
where $`ϵ`$ is a positive constant and $`P`$ is a Hermitian operator whose exact form is, for the moment, left unspecified. Addition of the $`ϵP`$ term will perturb the eigenvalues of $`H_\mathrm{\Lambda }`$, resulting in energy splitting between orthogonal ground states of the original (unperturbed) Hamiltonian. Consequently, we define

$$\mathrm{\Delta }E_\mathrm{\Lambda }(ϵ):=\underset{|\underset{¯}{\sigma }_\mathrm{\Lambda }^g}{\mathrm{max}}\underset{¯}{\sigma }|H_\mathrm{\Lambda }(ϵ)\underset{¯}{\sigma }.$$
(3)
The basic idea behind a self-correcting quantum spin system thus boils down to the following. Information is stored in the ground-state eigenspace of the unperturbed Hamiltonian $`H_\mathrm{\Lambda }`$. The multiplicity $`m`$ is, obviously, dictated by the desired storage capacity: when $`m=2^k`$, our ”ground-state memory cell” will hold $`k`$ qubits. Errors will cause some of the information to leak out into excited states. In order for error correction to take place, the system should be able to recover its ground state from sufficiently weak perturbations at sufficiently low temperature (the fact that we have to work with low temperatures is clear since we deal with the ground state). That is, we hope that there exists a threshold value $`ϵ_0`$ such that

$$\underset{\mathrm{\Lambda }^\nu }{lim}\mathrm{\Delta }E_\mathrm{\Lambda }(ϵ)=0,ϵϵ_0.$$
(4)
This condition, however, is necessary but not sufficient for error correction. It may happen that the $`m`$-fold degeneracy of the ground state does not survive in the thermodynamic limit, and we still lose some storage capacity. We therefore require that the ground state of the perturbed Hamiltonian remains $`m`$-fold degenerate for all $`ϵϵ_0`$ in the thermodynamic limit. In the next section we elaborate further on these requirements for self-correction and show that they are quite easy to fulfill in a wide variety of quantum spin systems. Before we proceed, we must remark that these are essentially ”bare-bones” desiderata; more elaborate performance measures will also be investigated elsewhere .

## 3 Putting it together

The main question is: which restrictions ensue on the unperturbed Hamiltonian $`H_\mathrm{\Lambda }`$ and on the perturbation $`P`$? It turns out that this question can be answered using the same methods that are employed for constructing low-temperature phase diagrams for classical spin systems with quantum perturbations . Thus the Hamiltonian $`H_\mathrm{\Lambda }`$ can be comprised by $`n`$-spin interactions (for fixed $`n`$) that satisfy the Peierls condition : the energy ”cost” of a local perturbation $`\omega ^{}`$ of a translationally invariant ground state $`\omega `$ is of the order of the surface area of the region that encloses the part of the lattice on which $`\omega `$ and $`\omega ^{}`$ differ. Additionally, the unperturbed Hamiltonian is assumed to have a spectral gap $`g>0`$ (i.e., its first nonzero eigenvalue $`g`$). Admissible perturbations are formed by sums of translates of an arbitrary Hermitian operator $`P_0`$, whose support (the set of sites on which the action of $`P_0`$ is nontrivial) is finite and encloses the origin of the lattice. Thus

$$P=\underset{\mathrm{}\mathrm{\Lambda }}{}P_{\mathrm{}},$$
(5)
where $`P_{\mathrm{}}=\gamma (\mathrm{})P_0`$ with $`\gamma (\mathrm{})`$ being the automorphism induced by the translation of the lattice $`\mathrm{\Lambda }`$ that maps the origin $`0`$ to the site $`\mathrm{}`$ and respects the periodic boundary conditions. Also, both the unperturbed and the perturbed Hamiltonians are assumed to be invariant under unitary transformations induced by a symmetry group acting transitively on the set $`\{|\underset{¯}{\sigma }\}_\mathrm{\Lambda }^g`$.

Assuming these conditions are satisfied, we invoke the following theorem : there exists a constant $`ϵ_0`$ such that, for all $`ϵϵ_0`$, the perturbed spin system has $`m`$ translationally invariant ground states in the thermodynamic limit. Furthermore, if the $`m`$ translationally invariant ground states of the unperturbed Hamiltonian are invariant under some additional symmetries, these invariance properties carry over to the ground states of the perturbed system. The threshold value $`ϵ_0`$ of the parameter $`ϵ`$ is determined by developing a low-temperature expansion  of the perturbed partition function using a modified Lie-Trotter product formula,

$$\mathrm{tr}e^{\beta H_\mathrm{\Lambda }(ϵ)}=\underset{N\mathrm{}}{lim}\mathrm{tr}\left\{\left[\left(1\mathrm{I}\frac{ϵP}{N}\right)e^{(1/N)H_\mathrm{\Lambda }}\right]^{N\beta }\right\}.$$
(6)
The trace is expanded in the basis $`\{|\underset{¯}{\sigma }\}`$, thus allowing for combinatorial analysis on a ”space-time” grid, where the space axis is labelled by lattice sites $`\mathrm{}`$ and the time axis is labelled by values $`0,(1/N)\beta ,(2/N)\beta ,\mathrm{},\beta `$. The perturbation theory is controlled by a suitable coarse-graining of the time axis, which then allows to determine the threshold value $`ϵ_0`$ that will render the contributions of the perturbation terms $`P_{\mathrm{}}`$ sufficiently small. (See  for details; we only remark here that translation invariance requirements can be lifted, the perturbation theory still goes through .) Similar space-time analysis of the error rate has been described heuristically by Dennis et al. .

Another important issue is the following: while the infinite-volume ground state of the perturbed system retains the degeneracy of the original (unperturbed) ground state, the degeneracy may be lost in the case of finite lattice size. This phenomenon is referred to as ”obscured symmetry breaking” , whereby low-lying eigenstates of the finite-system Hamiltonian converge to additional ground states in the thermodynamic limit. In this case we will have, for any finite $`\mathrm{\Lambda }`$,

$$\mathrm{\Delta }E_\mathrm{\Lambda }(ϵ)>0.$$
(7)
It is therefore important to obtain an estimate of the convergence rate in (4); namely, given some $`\delta >0`$, find $`N_0`$ such that

$$\mathrm{\Delta }E_\mathrm{\Lambda }(ϵ)<\delta ,\left|\mathrm{\Lambda }\right|N_0.$$
(8)
Knowing the convergence rate allows us to appraise the resources needed to implement error correction with the desired accuracy $`\delta `$. In this respect, an estimate of the form

$$\mathrm{\Delta }E_\mathrm{\Lambda }(ϵ)=e^{c\left|\mathrm{\Lambda }\right|},ϵϵ_0,$$
(9)
where the constant $`c`$ depends on $`ϵ`$, would be ideal — an exponential gain in error-correction accuracy could then be achieved with polynomial resources. This exponential convergence rate is, in fact, one of the most attractive features of Kitaev’s construction in Ref. .

On the other hand, the rate at which $`\mathrm{\Delta }E_\mathrm{\Lambda }(ϵ)`$ converges to zero is determined by the unperturbed Hamiltonian $`H_\mathrm{\Lambda }`$, the perturbation $`P`$, and the perturbation strength $`ϵ`$. It is therefore important to know what we can expect in a generic setting. Obviously, the exponential behavior, as in Eq. (9), is optimal, but it may as well turn out that the particular implementation (e.g., with a different Hamiltonian) does not allow for it. We can, however, hope for a slower (but still quite decent) convergence rate. According to a theorem of Horsch and von der Linden , certain quantum spin systems possess low-lying eigenstates of the finite-lattice Hamiltonian with

$$\mathrm{\Delta }E_\mathrm{\Lambda }(ϵ)=c/\left|\mathrm{\Lambda }\right|.$$
(10)
The conditions for this to hold are the following. There has to exist an order observable $`O_\mathrm{\Lambda }`$ of the form

$$O_\mathrm{\Lambda }=\underset{\mathrm{}\mathrm{\Lambda }}{}O_{\mathrm{}},$$
(11)
where each $`O_{\mathrm{}}`$ is a Hermitian operator such that $`[O_{\mathrm{}},O_{\mathrm{}^{}}]=0`$. Furthermore, for any interaction term $`\mathrm{\Phi }_{M\mathrm{\Lambda }}`$ in the perturbed Hamiltonian $`H_\mathrm{\Lambda }(ϵ)`$ (these also include the perturbation terms), we will have $`[\mathrm{\Phi }_M,O_{\mathrm{}}]=0`$ unless $`\mathrm{}M`$. The operators $`\mathrm{\Phi }_M`$ and $`O_{\mathrm{}}`$ are required to be uniformly bounded (in $`M`$ and $`\mathrm{}`$ respectively), and the cardinality of the support set $`M`$ must not exceed some fixed constant $`C`$ (the latter condition has also to be fulfilled for the perturbation theory described above to converge). Finally, if $`|\phi `$ is an eigenstate of $`H_\mathrm{\Lambda }(ϵ)`$, then we must have $`\phi |O_\mathrm{\Lambda }|\phi =0`$, but $`\phi |O_\mathrm{\Lambda }^2|\phi \zeta \left|\mathrm{\Lambda }\right|^2`$ (here the constant $`\zeta `$ depends on $`O_{\mathrm{}}`$). The latter conditions are taken as manifestations of ”obscured symmetry breaking” . Examples of systems for which the Horsch-von der Linden theorem holds include  the Ising model in the transverse magnetic field or the Heisenberg antiferromagnet with a Néel order.

## 4 Conclusion

Where does it all take us? It appears, from the discussion in the previous section, that any quantum spin system, whose Hamiltonian is formed by mutually commuting $`n`$-body interactions that satisfy the Peierls condition, can recover from sufficiently weak quantum (i.e., off-diagonal) perturbations at low temperatures. The admissible perturbations can be either finite-range , or exponentially decaying . Under these (quite general) conditions, it follows from rigorous perturbation theory for quantum spin systems that there exists a critical perturbation strength $`ϵ_0`$, such that for all $`ϵ<ϵ_0`$ the degeneracy and symmetry properties of the ground state of the original (unperturbed) system survive in the thermodynamic limit. Furthermore, even if ground-state degeneracy is removed by perturbation of the finite-size system, the effect of the error (perturbation) is effectively ”washed out” in the thermodynamic limit, as the low-lying excited states of the perturbed system converge to additional ground states.

However, the systems we have considered in this Letter are assumed to have classical Hamiltonians and discrete symmetries. What about truly quantum Hamiltonians and continuous symmetry (e.g., the quantum Heisenberg model)? The situation here is not so easy. For instance, it is apparent from our discussion that, in order to be self-correcting, the perturbed system must exhibit an order-disorder transition as the parameter $`ϵ`$ is tuned: in the ”ordered phase,” error correction is possible; in the ”disordered phase,” occurrence of errors results in irrevocable loss of information. (This has already been noted by Dennis et al. .) Since we require the ground state of the perturbed Hamiltonian to exhibit the same degeneracy as the corresponding state of the unperturbed Hamiltonian, it makes sense to talk about spontaneously broken symmetry in the ordered phase (i.e., when $`ϵ<ϵ_0`$). However, according to the so-called Goldstone theorem , symmetry cannot be broken in a system with continuous symmetry and a gap. It would certainly be worthwile to explore physical error correction in systems with continuous symmetries as well, but the models in which it can work will not be as easy to find.

## Acknowledgments

The author has benefitted from a few discussions with V. Moroz. This work was supported by the U.S. Army Research Office through MURI grant DAAD19-00-1-0177.


format: text

text:
# References

Integrable Structure behind WDVV Equations

Henrik Aratyn
Department of Physics,

University of Illinois at Chicago,

845 W. Taylor St.,

Chicago, IL 60607-7059

e-mail: aratyn@uic.edu

Johan van de Leur
Mathematical Institute,

University of Utrecht,

P.O. Box 80010, 3508 TA Utrecht,

The Netherlands

e-mail: vdleur@math.uu.nl

Abstract

An integrable structure behind Witten–Dijkgraaf–Verlinde–Verlinde (WDVV) equations is identified with reduction of a Riemann-Hilbert problem for a homogeneous $`\widehat{GL}(N,\text{})`$ loop group. Reduction requires the dressing matrices to be fixed points of a loop group automorphism of order two resulting in a sub-hierarchy of $`\widehat{gl}(N,\text{})`$ hierarchy containing only odd symmetry flows. The model possesses Virasoro symmetry and imposing Virasoro constraints ensures homogeneity property of the Darboux-Egoroff structure. Dressing matrices of the reduced model provide solutions of the WDVV equations.

1. Introduction

Massive topological field theories can be classified locally by the Darboux–Egoroff system of differential equations:

$$\frac{}{u_k}\beta _{ij}=\beta _{ik}\beta _{kj},ikj,\underset{k=1}{\overset{N}{}}\frac{}{u_k}\beta _{ij}=0,ij,$$
(1)
$$\beta _{ij}=\beta _{ji},$$
(2)
where in addition one also assumes that

$$\underset{k=1}{\overset{N}{}}u_k\frac{}{u_k}\beta _{ij}=\beta _{ij}.$$
(3)
This was observed by B. Dubrovin in , . The $`\beta _{ij}`$, $`1i,jN`$ with $`ij`$ are called the “rotation coefficients”. To be more precise, these Darboux–Egoroff equations are the compatibility conditions of the following linear system of differential equations depending on a spectral parameter $`\lambda `$:

$`{\displaystyle \frac{}{u_i}}A_{jk}(u,\lambda )`$ $`=`$ $`\beta _{ji}(u)A_{ik}(u,\lambda );ij=1,\mathrm{},N`$ (4)
$`{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{}{u_j}}A_{ik}(u,\lambda )`$ $`=`$ $`\lambda A_{ik}(u,\lambda );i=1,\mathrm{},N`$ (5)
for each $`k=1,\mathrm{},N`$. Let $`V(u)=(\beta _{ij}(u)(u_ju_i))_{1i,jn}`$ and assuming that this matrix is diagonalizable, one can construct a local semisimple Frobenius manifold, which provides an algebraic formulation of a massive topological field theory as follows. Find an invertible power series

$$A(u,\lambda )=A_0(u)+A_1(u)\lambda +A_2(u)\lambda ^2+\mathrm{},$$
(6)
which is a solution of this system, then

$$\underset{k=1}{\overset{N}{}}u_k\frac{}{u_k}A_0(u)=V(u)A_0(u).$$
(7)
It is then straightforward to check that $`A_0(u)^1V(u)A_0(u)`$ is a constant matrix, hence there exists an invertible complex matrix $`S`$ such that

$$\mu =\underset{i=1}{\overset{N}{}}\mu _iE_{ii}=S^1A_0(u)^1V(u)A_0(u)S$$
(8)
is a diagonal matrix. Now let $`A(u)=A_0(u)S`$ with $`A(u)=(a_{ij}(u))_{1i,jN}`$, then

$$\underset{k=1}{\overset{N}{}}u_k\frac{}{u_k}A(u)=A(u)\mu .$$
(9)
Dubrovin then showed  that there exists on the domain $`u_iu_j`$ and $`a_{11}a_{21}\mathrm{}a_{N1}0`$ a local semisimple Frobenius manifold with scaling dimensions $`\mu _\alpha \mu _1`$ (see also ). The diagonal metric is given by

$$ds^2=\underset{i=1}{\overset{n}{}}a_{i1}^2(u)(du_i)^2.$$
(10)
The flat coordinates $`x^\alpha `$ are determined by

$$\underset{\beta }{}\eta _{\alpha \beta }\frac{x^\beta (u)}{u_i}=a_{i1}(u)a_{i\alpha }(u),$$
(11)
where

$$\eta _{\alpha \beta }=\underset{i=1}{\overset{n}{}}\alpha _{i\alpha }(u)\alpha _{i\beta }(u).$$
(12)
Finally, the structure constants

$$c_{\alpha \beta \gamma }(x(u))=\frac{^3F(x(u))}{x^\alpha x^\beta x^\gamma }=\underset{i=1}{\overset{n}{}}\frac{a_{i\alpha }(u)a_{i\beta }(u)a_{i\gamma }(u)}{a_{i1}(u)}$$
(13)
are the third order derivatives of the WDVV prepotential $`F(x)`$.

In this paper, we use a Riemann–Hilbert problem for the loop group of $`GL(N)`$ to obtain a similar structure. The symmetry conditions (2) for the rotation coefficients are obtained by using the Cartan involution on the loop group. The model possesses Virasoro symmetry and imposing Virasoro constraints ensures the homogeneity property (3) of the Darboux-Egoroff structure.

2. Symmetry flows, Riemann-Hilbert problem and the $`\tau `$ function

We will introduce the $`G=\widehat{GL}(N)`$ symmetry flows in terms of (extended) Riemann-Hilbert problem involving two subgroups of the Lie loop group $`G`$ defined as :

$$G_{}=\left\{gG\right|g(\lambda )=1+\underset{i<0}{}g^{(i)}\},G_+=\left\{gG\right|g(\lambda )=\underset{i0}{}g^{(i)}\},$$
(14)
here $`g_i`$ has grading $`i`$ with respect to a homogeneous gradation defined by a derivation $`\lambda d/d\lambda `$. It also holds that $`G_+G_{}=I`$. Let the loop algebra corresponding to $`G`$ be $`\widehat{𝒢}=\widehat{gl}(N)`$. This algebra splits into the direct sum $`\widehat{𝒢}=\widehat{𝒢}_+\widehat{𝒢}_{}`$ with respect to the homogeneous gradation, where $`\widehat{𝒢}_\pm `$ are Lie algebras associated with the subgroups $`G_\pm `$.

Let the “un-dressed” wave (matrix) function be

$$\mathrm{\Psi }_0(𝐮,\lambda )=\mathrm{exp}\left(\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}E_{jj}^{(n)}u_j^{(n)}\right)=\underset{j=1}{\overset{N}{}}E_{jj}e^{_{n=1}^{\mathrm{}}\lambda ^nu_j^{(n)}}$$
(15)
We now define an extended Riemann-Hilbert factorization problem for the homogeneous gradation :

$$\mathrm{\Psi }_0(𝐮,\lambda )g(\lambda )=\mathrm{\Theta }^1(𝐮,\lambda )M(𝐮,\lambda )$$
(16)
where $`g:S^1G_{}G_+`$ while $`\mathrm{\Theta }G_{}`$, $`MG_+`$. We use the multi-time notation with $`(𝐮)=(𝐮_\mathrm{𝟏},\mathrm{},𝐮_𝐍)`$ to denote $`N`$ multi-flows $`𝐮_𝐣`$. Each argument $`𝐮_𝐣`$, $`j=1,\mathrm{},N`$ is an abbreviated notation for the multi-flows $`u_j^{(n)}`$ with $`n`$ running from $`1`$ to $`\mathrm{}`$.

For the symmetry flows one derives from (16) expressions:

$`{\displaystyle \frac{}{u_j^{(n)}}}\mathrm{\Theta }(𝐮,\lambda )`$ $`=`$ $`\left(\mathrm{\Theta }E_{jj}^{(n)}\mathrm{\Theta }^1\right)_{}\mathrm{\Theta }(𝐮,\lambda )`$ (17)
$`{\displaystyle \frac{}{u_j^{(n)}}}M(𝐮,\lambda )`$ $`=`$ $`\left(\mathrm{\Theta }E_{jj}^{(n)}\mathrm{\Theta }^1\right)_+M(𝐮,\lambda )`$ (18)
where $`(\mathrm{})_\pm `$ denote the projections onto $`\widehat{𝒢}_\pm `$. Note, that $`_{j=1}^N\mathrm{\Theta }/u_j^{(n)}=0`$ for any $`n>0`$.

In the following we will denote for brevity $`u_j=u_j^{(1)}`$ and $`_j=/u_j^{(1)}`$ for $`j=1,\mathrm{},N`$.

It is natural to associate to the Riemann-Hilbert factorization a one-form :

$$𝒥=\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}\mathrm{Res}_\lambda \left(tr\left(\mathrm{\Theta }^1\frac{d\mathrm{\Theta }}{d\lambda }E_{jj}^{(n)}\right)\right)du_j^{(n)}=\mathrm{Res}_\lambda \left(tr\left(\mathrm{\Theta }^1\frac{d\mathrm{\Theta }}{d\lambda }\mathrm{\Psi }_0^1d\mathrm{\Psi }_0\right)\right)$$
(19)
which due to equation (17) is closed : $`d𝒥=0`$ with $`d=_{j=1}^N_{n=1}^{\mathrm{}}du_j^{(n)}\frac{}{u_j^{(n)}}`$. Hence, locally, $`𝒥`$ can be written as $`d`$ of a zero-form being some function of the dynamical variables of the model. Conventionally, this function is denoted as the logarithm of the $`\tau `$-function:

$$𝒥=d\mathrm{log}\tau =\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}\frac{\mathrm{log}\tau }{u_j^{(n)}}du_j^{(n)}$$
(20)
This equation can effectively be integrated to yield an explicit expression for the $`\tau `$-function. To do this it is required to extend the formalism from the loop groups to the centrally extended groups. As in case of the AKNS hierarchy with the $`\tau `$-function defined in  we embed the framework in the centrally extended group with elements denoted as $`\widehat{g}`$ and propose :

$$\tau \left(𝐮\right)=0|\widehat{\mathrm{\Psi }_0g}|0$$
(21)
where $`|0`$ is the highest weight state such that $`X_0|0=0`$ and $`0|X_0=0`$ and furthermore $`0|\widehat{k}|0=1`$, where $`\widehat{k}`$ denotes the direction associated to the central element within the affine Kac-Moody algebra.

Let $`\theta ^{(1)}`$ be the grade $`1`$ term of expansion:

$$\mathrm{\Theta }=1+\theta ^{(1)}\lambda ^1+\theta ^{(2)}\lambda ^2+\mathrm{}.$$
(22)
From equations (19) and (20) (for $`n=1`$) one finds :

$$(\theta ^{(1)})_{ii}=_i\mathrm{log}\tau .$$
(23)
Introduce now the “rotation cefficients” $`\beta _{ij}`$ with $`ij`$ as the off-diagonal elements of $`\theta ^{(1)}`$ :

$$\beta _{ij}=(\theta ^{(1)})_{ij},ij=1,\mathrm{},N$$
(24)
For $`\theta ^{(1)}`$ matrix we find from (17) :

$$_j\theta ^{(1)}=[E_{jj},\theta ^{(2)}]+[\theta ^{(1)},E_{jj}]\theta ^{(1)}$$
(25)
from which it follows that :

$$_j(\theta ^{(1)})_{ik}=\left(\theta ^{(1)}\right)_{ij}(\theta ^{(1)})_{jk},ji,jk.$$
(26)
Accordingly, the “rotation coefficients” $`\beta _{ij}`$ given by (24) satisfy (1).

As found in , the matrix elements of the dressing matrix $`\mathrm{\Theta }=(\theta _{ij})_{1i,jN}`$ are given by:

$$\theta _{ij}(𝐮,\lambda )=\frac{1}{\tau (𝐮)}\times \{\begin{array}{cc}\frac{1}{\lambda }\beta _{ij}(𝐮_𝐣[\lambda ^1])\tau (𝐮_𝐣[\lambda ^1])\hfill & \hfill ij\\ \tau (𝐮_𝐢[\lambda ^1])\hfill & \hfill i=j\end{array}$$
(27)
where the multi-flow arguments $`𝐮_𝐣`$ were shifted according to $`(𝐮_𝐣[\lambda ^1])=(u_j^{(1)}1/\lambda ,u_j^{(2)}1/2\lambda ^2,\mathrm{})`$.

The matrix coefficients of the wave matrix $`\mathrm{\Psi }`$:

$$\mathrm{\Psi }(𝐮,\lambda )=\mathrm{\Theta }(𝐮,\lambda )\mathrm{\Psi }_0(𝐮,\lambda )$$
(28)
and the matrix coefficients of $`M(𝐮,\lambda )=\mathrm{\Psi }(𝐮,\lambda )g(\lambda )`$ satisfy the linear system (4)-(5), which are the compatibility conditions of equations (1). This is the Darboux-Egoroff system provided we also impose the symmetry condition (2) on the rotation coefficients. To do this we reduce the integrable hierarchy associated with the flow equation equations (17),(18) using the automorphism $`\sigma `$ :

$$\sigma \left(X(\lambda )\right)=\left(\left(X(\lambda )\right)^T\right)^1,XG=\widehat{GL}(N)$$
(29)
which leaves the evolution equations (17)-(18) invariant for odd flows only (labeled by $`n`$ being an odd integer). Accordingly, we define the integrable sub-hierarchy by constraining the dressing matrices $`\mathrm{\Theta }(𝐮,\lambda )`$ and $`M(𝐮,\lambda )`$ to the fixed points of the loop group automorphism $`\sigma `$ (29) :

$$\mathrm{\Theta }^1(𝐮,\lambda )=\mathrm{\Theta }^T(𝐮,\lambda ),M^1(𝐮,\lambda )=M^T(𝐮,\lambda )$$
(30)
with $`\mathrm{\Theta }(𝐮,\lambda )`$ and $`M(𝐮,\lambda )`$ depending only on odd coordinates $`𝐮`$: ( $`u_j^{(2k+1)},k=0,1,\mathrm{}`$). The odd flows of the reduced sub-hierarchy preserve the conditions (30). The fixed points of the automorphism $`\sigma `$ form a subgroup of $`G=\widehat{GL}(N)`$, called a twisted loop group of $`GL(N)`$ (see e.g. ). For the first term $`\theta ^{(1)}`$ of the expansion (22) the constraint (30) implies that $`\theta ^{(1)}=\theta ^{(1)T}`$. Accordingly, the “rotation coefficients” $`\beta _{ij}`$ associated to the subhierarchy via (24) are symmetric and give rise to the Darboux-Egoroff metrics.

3. Condition of conformal invariance and homogeneity

We propose here the Virasoro constraints in the setting of the Riemann-Hilbert problem approach to the integrable models. In this setting the Virasoro constraints are seen to arise as symmetries of the hierarchy.

The action of the Virasoro symmetry flows on the dressing matrices is given by:

$`\delta _k^V\mathrm{\Theta }(𝐮,\lambda )`$ $`=`$ $`(\mathrm{\Theta }L_k\mathrm{\Theta }^1)_{}\mathrm{\Theta }=(\mathrm{\Theta }l_k\mathrm{\Theta }^1)_{}\mathrm{\Theta }{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}{\displaystyle \frac{\mathrm{\Theta }}{u_j^{(n+k)}}}`$ (31)
$`\delta _k^VM(𝐮,\lambda )`$ $`=`$ $`(\mathrm{\Theta }L_k\mathrm{\Theta }^1)_+M=(\mathrm{\Theta }l_k\mathrm{\Theta }^1)_+M{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}{\displaystyle \frac{M}{u_j^{(n+k)}}}`$ (32)
where

$`L_k`$ $`=`$ $`\lambda ^{k+1}{\displaystyle \frac{d}{d\lambda }}+{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}E_{jj}^{(k+n)},k1`$
$`=`$ $`\lambda ^{k+1}{\displaystyle \frac{d}{d\lambda }}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nU^{(n)}\lambda ^{k+n},`$
$`U^{(n)}`$ $`=`$ $`\mathrm{diag}(u_1^{(n)},\mathrm{},u_N^{(n)})={\displaystyle \underset{i=1}{\overset{N}{}}}u_i^{(n)}E_{ii}`$
and $`l_k=\lambda ^{k+1}d/d\lambda `$. Both $`l_k`$ and $`L_k`$ for $`k1`$ generate a subalgebra of the Witt (centerless Virasoro) algebra:

$$[l_k,l_r]=(kr)l_{k+r},[L_k,L_r]=(kr)L_{k+r}$$
(34)
The signs in eqs. (31) and (32) have been chosen in such a way that the transformation $`\delta _k^V`$ also satisfies the Virasoro algebra $`[\delta _k^V,\delta _r^V]=(kr)\delta _{k+r}^V`$ when applied on the dressing matrices $`\mathrm{\Theta }`$ and $`M`$. Furthermore, the transformations $`\delta _k^V`$ commute with the flows $`/u_j^{(n)}`$:

$$[\delta _k^V,\frac{}{u_j^{(n)}}]\mathrm{\Theta }(𝐮,\lambda )=0$$
(35)
Accordingly, the Virasoro transformations from eqs. (31) and (32) constitute additional symmetries of the integrable model defined by the Riemann-Hilbert factorization problem. One can show that only the even Virasoro flows will preserve the conditions (30) and define additional symmetries of the reduced sub-hierarchy associated with the twisted loop group of $`GL(N)`$. We can therefore consistently impose the constraints $`\delta _k^V\mathrm{\Theta }(𝐮,\lambda )=0`$ with even $`k`$. For $`k=0`$ the equation (31) specializes to:

$$\delta _0^V\mathrm{\Theta }(𝐮,\lambda )=\lambda \frac{d}{d\lambda }\mathrm{\Theta }\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}\frac{\mathrm{\Theta }}{u_j^{(n)}}$$
(36)
as follows from identity $`(\mathrm{\Theta }\lambda d\mathrm{\Theta }^1/d\lambda )_{}\mathrm{\Theta }=\lambda d\mathrm{\Theta }/d\lambda `$. Hence for $`k=0`$ the constraint $`\delta _k^V\mathrm{\Theta }(𝐮,\lambda )=0`$ takes a form

$$\left(\lambda \frac{d}{d\lambda }E_0\right)\mathrm{\Theta }(𝐮,\lambda )=0$$
(37)
where we have introduced the Euler vector field $`E_0`$:

$$E_0=\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}\frac{}{u_j^{(n)}}$$
(38)
The homogeneity condition (37) is also a condition for the conformal Frobenius manifold.

By projecting in equation (36) on the terms with $`1`$ grade we obtain :

$$\delta _0^V\theta ^{(1)}=\theta ^{(1)}\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}\frac{\theta ^{(1)}}{u_j^{(n)}}$$
(39)
where $`\theta ^{(1)}`$ is the grade $`1`$ term of expansion: (22).

Hence, the condition (37) implies that

$$E_0\theta ^{(1)}=\theta ^{(1)}$$
(40)
Plugging relations (23) and (24) into equation (39) we obtain :

$`\delta _0^V_i\mathrm{log}\tau `$ $`=`$ $`_i\delta _0^V\mathrm{log}\tau =_i{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}{\displaystyle \frac{\mathrm{log}\tau }{u_j^{(n)}}},i=1,\mathrm{},N`$ (41)
$`\delta _0^V\beta _{kl}`$ $`=`$ $`\beta _{kl}{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}{\displaystyle \frac{\beta _{kl}}{u_j^{(n)}}}`$ (42)
4. Virasoro flows on the $`\tau `$ function

This subsection studies action of the general Virasoro transformation on the $`\tau `$-function of the underlying integrable system.

Recall an expression for derivatives of $`\mathrm{log}\tau `$ given in (19) and (20). Applying $`_i`$ on both sides of eq. (20) we obtain an identity:

$$\frac{^2\mathrm{log}\tau }{u_j^{(n)}u_i}=\mathrm{Res}_\lambda \left(tr\left(\mathrm{\Theta }E_{jj}^{(n)}\mathrm{\Theta }^1E_{ii}\right)\right)i,j=1,\mathrm{},N,n>0$$
(43)
Combining eqs. (31) and (19)-(20) we find

$$\delta _k^V\frac{\mathrm{log}\tau }{u_j^{(n)}}=\mathrm{Res}_\lambda \left(tr\left(\mathrm{\Theta }E_{jj}^{(n)}\mathrm{\Theta }^1\frac{d}{d\lambda }(\mathrm{\Theta }L_k\mathrm{\Theta }^1)_{}\right)\right)$$
(44)
Since $`\lambda d(\lambda ^1)/d\lambda =\lambda ^1`$ the last equation becomes

$$\delta _k^V_i\mathrm{log}\tau =\mathrm{Res}_\lambda \left(tr\left(E_{ii}\mathrm{\Theta }l_k\mathrm{\Theta }^1\right)\right)\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}\mathrm{Res}_\lambda \left(tr\left(\mathrm{\Theta }E_{jj}^{(k+n)}\mathrm{\Theta }^1E_{ii}\right)\right)$$
(45)
Using the identity (43) we rewrite the last equation as

$$\delta _k^V_i\mathrm{log}\tau =\mathrm{Res}_\lambda \left(tr\left(E_{ii}\mathrm{\Theta }l_k\mathrm{\Theta }^1\right)\right)\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}_i\frac{\mathrm{log}\tau }{u_j^{(k+n)}}$$
(46)
The cases $`k=0,1`$ were calculated and eq. (45) gives :

$$_i\delta _k^V\mathrm{log}\tau =_i\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}\frac{\mathrm{log}\tau }{u_j^{(n+k)}},k=0,1$$
(47)
and

$$_i\delta _1^V\mathrm{log}\tau =_i\underset{j=1}{\overset{N}{}}\underset{n=2}{\overset{\mathrm{}}{}}nu_j^{(n)}\frac{\mathrm{log}\tau }{u_j^{(n1)}}.$$
(48)
Only $`k=0`$ case applies to the reduced case.

From eqs. (31) and (32) we find :

$`\delta _k^V\left(\mathrm{\Theta }^1M\right)`$ $`=`$ $`\mathrm{\Theta }^1(\mathrm{\Theta }L_k\mathrm{\Theta }^1)_{}M\mathrm{\Theta }^1(\mathrm{\Theta }L_k\mathrm{\Theta }^1)_{}M=L_k\left(\mathrm{\Theta }^1M\right)`$ (49)
$`=`$ $`\lambda ^{k+1}{\displaystyle \frac{d}{d\lambda }}\left(\mathrm{\Theta }^1M\right){\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}E_{jj}^{(n+k)}\left(\mathrm{\Theta }^1M\right)`$
Recalling the factorization problem (16) we find

$$\delta _k^V\left(\mathrm{\Theta }^1M\right)=\mathrm{\Psi }_0(𝐮,\lambda )\lambda ^{k+1}\frac{dg}{d\lambda }$$
(50)
since

$$\lambda ^{k+1}\frac{d}{d\lambda }\mathrm{\Psi }_0(𝐮,\lambda )=\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}E_{jj}^{(n+k)}\mathrm{\Psi }_0(𝐮,\lambda )$$
(51)
Recall, the definition (21) of the $`\tau `$-function $`\tau _g\left(𝐮\right)=0|\widehat{\mathrm{\Psi }_0g}|0`$, rewritten here in the form which explicitly shows dependence of $`\tau `$ on the loop group element $`g(\lambda )`$.

It now appears that

$$\delta _k^V\tau _g=0|\lambda ^{k+1}\widehat{\mathrm{\Psi }_0\frac{dg}{d\lambda }}|0=\tau _{\lambda ^{k+1}dg/d\lambda }$$
(52)
So under infinitesimal $`S^1`$ diffeomorpism $`gg+\lambda ^{k+1}dg/d\lambda `$ the $`\tau `$-function transforms as

$$\tau _g\tau _{g+\lambda ^{k+1}dg/d\lambda }=\tau _g+\delta _k^V\tau _g$$
(53)
5. Virasoro symmetry and dressing

Let

$$\delta _j^{(n)}=\frac{}{u_j^{(n)}}E_{jj}^{(n)},\delta _j=\delta _j^{(1)}=_j\lambda E_{jj},j=1,\mathrm{},N$$
(54)
be a set of commuting operators :

$$[\delta _j^{(n)},\delta _i^{(k)}]=0,$$
(55)
The above operators and the Virasoro operators $`L_k`$ from (S0.Ex1) commute with each other:

$$[\delta _j^{(n)},L_k]=0,i,j=1,\mathrm{},N$$
(56)
Both operators (54) and (S0.Ex1) annihilate the bare wave function $`\mathrm{\Psi }_0`$ from (15) :

$$L_k\mathrm{\Psi }_0=0,\delta _j^{(n)}\mathrm{\Psi }_0=0,j=1,\mathrm{},N.$$
(57)
Similarly the operators:

$$𝒟_j^{(n)}=\mathrm{\Theta }\delta _j^{(n)}\mathrm{\Theta }^1,_k=\mathrm{\Theta }L_k\mathrm{\Theta }^1,$$
(58)
satisfy the commutation relations (56) and (34) and annihilate the wave (matrix) function $`\mathrm{\Psi }(𝐮,\lambda )`$ defined in (28) :

$$_k\mathrm{\Psi }=0,𝒟_j^{(n)}\mathrm{\Psi }=0,j=1,\mathrm{},N.$$
(59)
Since

$$\mathrm{\Theta }_j\mathrm{\Theta }^1=_j+\mathrm{\Theta }(_j\mathrm{\Theta }^1)=_j+\left(\mathrm{\Theta }\lambda E_{jj}\mathrm{\Theta }^1\right)_{}$$
(60)
we obtain

$$𝒟_j=𝒟_j^{(1)}=\mathrm{\Theta }\left(_j\lambda E_{jj}\right)\mathrm{\Theta }^1=_j\left(\mathrm{\Theta }\lambda E_{jj}\mathrm{\Theta }^1\right)_+=_j\lambda E_{jj}V_j$$
(61)
where

$$V_j[\theta ^{(1)},E_{jj}],(V_j)_{kl}=\left(\delta _{lj}\delta _{kj}\right)\beta _{kl}$$
(62)
for $`\theta ^{(1)}`$ defined in (22).

Consider

$$\lambda \frac{d\mathrm{\Psi }}{d\lambda }=\lambda \frac{d\mathrm{\Theta }}{d\lambda }\mathrm{\Psi }_0+\mathrm{\Theta }\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}E_{jj}^{(n)}\mathrm{\Psi }_0$$
(63)
Recall, now from (36) that the condition $`\delta _{k=0}^V\mathrm{\Theta }(𝐮,\lambda )=0`$ (37) implies

$$\lambda \frac{d\mathrm{\Theta }}{d\lambda }=(\mathrm{\Theta }\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}E_{jj}^{(n)}\mathrm{\Theta }^1)_{}\mathrm{\Theta }$$
(64)
and accordingly

$$\lambda \frac{d\mathrm{\Psi }}{d\lambda }=(\mathrm{\Theta }\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}E_{jj}^{(n)}\mathrm{\Theta }^1)_+\mathrm{\Psi }$$
(65)
Let us now study the operator $`_{k=1}=\mathrm{\Theta }L_1\mathrm{\Theta }^1`$ from (58) for $`\mathrm{\Theta }`$ satisfying the homogeneity condition (64).

$`_1`$ $`=`$ $`{\displaystyle \frac{d}{d\lambda }}+{\displaystyle \frac{d\mathrm{\Theta }}{d\lambda }}\mathrm{\Theta }^1+{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}\mathrm{\Theta }E_{jj}^{(n1)}\mathrm{\Theta }^1`$
$`=`$ $`{\displaystyle \frac{d}{d\lambda }}\lambda ^1{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}(\mathrm{\Theta }E_{jj}^{(n)}\mathrm{\Theta }^1)_{}+{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}\mathrm{\Theta }E_{jj}^{(n1)}\mathrm{\Theta }^1`$
$`=`$ $`{\displaystyle \frac{d}{d\lambda }}+\lambda ^1{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}nu_j^{(n)}(\mathrm{\Theta }E_{jj}^{(n)}\mathrm{\Theta }^1)_+`$
and for $`u_j^{(n)}=0,n>1`$ one obtains

$$_1=\frac{d}{d\lambda }+U+\lambda ^1[\theta ^{(1)},U]=\frac{d}{d\lambda }+U+\lambda ^1V$$
(67)
The components of the matrix $`V=[\theta ^{(1)},U]`$ are

$$V_{ij}=(u_ju_i)\theta _{ij}^{(1)}=(u_ju_i)\beta _{ij},i,j=1,\mathrm{},N$$
(68)
It follows from the construction that the wave (matrix) function $`\mathrm{\Psi }(𝐮,\lambda )`$ defined in (28) is annihilated by $`_1`$ and $`𝒟_j`$

$$_1\mathrm{\Psi }(𝐮,\lambda )=0\frac{d\mathrm{\Psi }}{d\lambda }=(U+\lambda ^1V)\mathrm{\Psi }$$
(69)
and

$$𝒟_j\mathrm{\Psi }(𝐮,\lambda )=0\frac{\mathrm{\Psi }}{u_j}=(\lambda E_{jj}+V_j)\mathrm{\Psi }.$$
(70)
Compatibility of the above equations amounts to $`[𝒟_j,_1]=0`$, which follows by dressing of $`[\delta _j,L_1]=0`$. From commutation relations $`[𝒟_j,_1]=0`$ and $`[𝒟_i,𝒟_j]=0`$ one obtains :

$`_jV`$ $`=`$ $`[V_j,V]`$ (71)
$`_jV_i`$ $`=`$ $`_iV_j+[V_j,V_i]`$ (72)
$`[V,E_{jj}]`$ $`=`$ $`[V_j,U]`$ (73)
6. Isomonodromic Tau function $`\tau _𝐈`$

We now consider the coefficients $`_j\mathrm{log}\tau `$ of the one-form $`𝒥`$ from (19) with the conformal constraint (64) imposed. Recall, that this constraint follows from the homogeneity condition (37) and in view of relations (41) and (42) it holds in the reduced case that $`E_0\tau =\mu \tau `$ and $`E_0\beta _{kl}=\beta _{kl}`$ where $`E_0=_{j=1}^Nu_j/u_j`$ and $`\mu `$ is independent of variables $`u_j,j=1,\mathrm{},N`$.

Plugging relation (64) into the defining equations (19) and (20) we obtain for $`n=1`$:

$$_j\mathrm{log}\tau =\mathrm{Res}_\lambda \left(tr\left(E_{jj}(\mathrm{\Theta }\underset{j=1}{\overset{N}{}}\underset{n=1}{\overset{\mathrm{}}{}}nu_j^{(n)}E_{jj}^{(n)}\mathrm{\Theta }^1)_{}\right)\right),j=1,\mathrm{},N$$
(74)
In the rest of the paper we set all $`u_i^{(k)}=0`$ for $`k>1`$ leaving the model defined only in terms of the canonical coordinates $`u_j,j=1,\mathrm{},N`$ and assume from now on that $`\sigma (g(\lambda ))=g(\lambda )`$ and that (30) holds.

The formula (74) simplifies significantly and takes the form of

$$_j\mathrm{log}\tau =\mathrm{Res}_\lambda \left(tr\left(E_{jj}(\mathrm{\Theta }\lambda U\mathrm{\Theta }^1)_{}\right)\right)=tr\left(E_{jj}\left([\theta ^{(2)},U]+\frac{1}{2}[\theta ^{(1)},[\theta ^{(1)},U]]\right)\right)$$
(75)
Using invariance of the trace $`tr(A[B,C])=tr([A,B]C)`$ we obtain expression:

$$_j\mathrm{log}\tau =\frac{1}{2}tr\left(V_jV\right)$$
(76)
with $`V_j,V`$ matrices defined in (62) and (68) and related through the commutation relation (73).

The equation (71) can be cast into the Hamiltonian form (see also )

$$_jV=\{V,H_j\}$$
(77)
with respect to the standard Poisson bracket $`\{,\}`$ on $`so(N)`$ and with the Hamiltonians

$$H_i=\frac{1}{2}\underset{ji}{}\frac{V_{ij}^2}{u_iu_j}.$$
(78)
The commutation relation (73) in components gives a relation

$$(u_ku_j)(V_i)_{kj}=\left(\delta _{ik}\delta _{ij}\right)V_{kj},kj=1,\mathrm{},N$$
(79)
which can be used to prove that

$$\frac{1}{2}tr\left(V_iV\right)=\underset{ji}{}\frac{V_{ij}V_{ji}}{u_iu_j}=\underset{ji}{}\frac{V_{ij}^2}{u_iu_j}=2H_i$$
(80)
From equations (71) and (72) it follows that $`dtr_{i=1}^N\left(V_iV\right)du_i=0`$ and the one-form $`_{i=1}^NH_idu_i`$ is closed. Locally, this one-form can be represented by the derivative of the logarithm of the so-called isomonodromy tau-function $`\tau _I`$

$$d\mathrm{log}\tau _I=\underset{i=1}{\overset{N}{}}H_idu_i=\frac{1}{2}d\mathrm{log}\tau .$$

which implies that upto a multiplicative constant :

$$\tau _I=1/\sqrt{\tau }.$$
(81)
This formula was also obtained in a different way in .

7. Frobenius structure

Expand, the dressing matrix $`M(u,\lambda )`$ according to the grading :

$$M(u,\lambda )=M_0(u)+M_1(u)\lambda +M_2(u)\lambda ^2+\mathrm{}$$
(82)
Then, according to the condition (30) the zero grade term $`M_0(u)=M(u,0)=(m_{ij}(u))_{1i,jN}`$ is an orthogonal matrix $`M_0^T(u)=M_0^1(u)`$ which as seen from eqs. (18) and (32) also satisfies the flow equations :

$$_jM_0(u)=[\theta ^{(1)}(u),E_{jj}]M_0(u)=V_j(u)M_0(u),\delta _0^VM_0(u)=V(u)M_0(u)$$
(83)
and hence also

$$\underset{j=1}{\overset{N}{}}u_j_jM_0(u)=V(u)M_0(u).$$
(84)
Define

$$𝒱=M_0^1VM_0,$$
(85)
it is easy to see that $`𝒱`$ satisfies:

$$_j𝒱=0,𝒱^T=𝒱,\delta _0^V𝒱=0,j=1,\mathrm{},N$$
(86)
and accordingly $`𝒱`$ is a constant. We can bring $`𝒱`$ into Jordan normal form $`\mu =\widehat{\mu }+R_0`$, with $`\widehat{\mu }`$ its semisimple part and $`R_0`$ its nilpotent part, i.e. there exists a complex invertible matrix $`S`$, such that

$$𝒱=S\mu S^1,$$
(87)
Now define

$$M(u)=M_0(u)S=(m_{ij}(u))_{1i,jN},$$
(88)
then using (83) we obtain:

$$_jM(u)=[\theta ^{(1)}(u),E_{jj}]M(u)=V_j(u)M(u),\delta _0^VM(u)=M(u)\mu $$
(89)
and hence also

$$\underset{j=1}{\overset{N}{}}u_j_jM(u)=M(u)\mu .$$
(90)
Now assume that $`𝒱`$ is diagonalizable, i.e. the matrix $`V(u)`$ has no nilpotent part, then

$$\mu =\underset{i=1}{\overset{N}{}}\mu _iE_{ii},\text{with}\mu _i=\mu _{N+1i}$$
(91)
and the columns of $`M(u)`$ are eigenvectors for the Euler vectorfield $`E_0`$. Hence, on the domain $`u_iu_j`$ and $`m_{11}m_{21}\mathrm{}m_{N1}0`$ we have a local semisimple Frobenius manifold (see e.g. ) with Lamé coefficients

$$h_i=m_{i1}.$$
(92)
Define,

$$\eta =(\eta _{\alpha \beta })_{1\alpha ,\beta N}=M^TM=S^TS,\text{and denote}\eta ^1=(\eta ^{\alpha \beta })_{1\alpha ,\beta N},$$
(93)
then $`\mu \eta +\eta \mu =0`$ and the derivatives with respect to the flat coordinates $`x^\alpha ,\alpha =1,\mathrm{},N`$ are given by

$$\frac{}{x^\alpha }=\underset{i=1}{\overset{N}{}}\frac{m_{i\alpha }}{m_{i1}}\frac{}{u_i}$$
(94)
and

$$\frac{}{u_i}=\underset{\alpha ,\beta =1}{\overset{N}{}}\eta ^{\alpha \beta }m_{i1}m_{i\beta }\frac{}{x^\alpha }.$$
(95)
The metric equals

$$ds^2=\underset{i=1}{\overset{N}{}}h_i^2(du_i)^2=\underset{\alpha ,\beta =1}{\overset{N}{}}\eta _{\alpha \beta }dx^\alpha dx^\beta ,\text{with}\eta _{\alpha \beta }=\underset{i=1}{\overset{N}{}}m_{i\alpha }m_{i\beta }.$$
(96)
and the structure constants are given by

$$c_{\alpha \beta \gamma }=\underset{i=1}{\overset{N}{}}\frac{m_{i\alpha }m_{i\beta }m_{i\gamma }}{m_{i1}}.$$
(97)
Since we have constructed everything in such a way that the columns of $`M(u)`$ are eigenvectors, our Frobenius manifold has scaling dimensions $`\mu _\alpha \mu _1`$.

Now, multiply the wave (matrix) function $`\mathrm{\Psi }(u,\lambda )`$, defined in (28), from the left by

$$M^1(u)=S^1M_0^1(u)=\eta ^1S^TM_0^T(u)=\eta ^1M^T(u)$$
(98)
and from the right by $`S`$ to obtain the wave function $`𝒴(u,\lambda )`$:

$$𝒴(u,\lambda )=M^1(u)\mathrm{\Psi }(u,\lambda )S=S^1M_0^1(u)M(u,\lambda )g^1(\lambda )S.$$
(99)
$`𝒴(u,\lambda )`$ satisfies, as in  :

$`𝒴(u,\lambda )^T\eta 𝒴(u,\lambda )`$ $`=`$ $`\eta `$ (100)
$`{\displaystyle \frac{𝒴(u,\lambda )}{u_i}}`$ $`=`$ $`\lambda _i(u)𝒴(u,\lambda )`$ (101)
$`{\displaystyle \frac{𝒴(u,\lambda )}{x^\alpha }}`$ $`=`$ $`\lambda C_\alpha (u)𝒴(u,\lambda )`$ (102)
$`{\displaystyle \frac{𝒴(u,\lambda )}{\lambda }}`$ $`=`$ $`\left(𝒰(u)+\lambda ^1\mu \right)𝒴(u,\lambda )`$ (103)
where

$`_i(u)`$ $`=`$ $`M^1(u)E_{ii}M(u)`$ (104)
$`𝒰(u)`$ $`=`$ $`M^1(u)UM(u)`$ (105)
$`\mu `$ $`=`$ $`M^1(u)V(u)M(u)=S^1𝒱S`$ (106)
and

$$C_\alpha (u)=\underset{\beta ,\gamma =1}{\overset{N}{}}c_{\alpha \beta }^\gamma E_{\beta \gamma }=\underset{\beta ,\gamma ,\delta =1}{\overset{N}{}}c_{\alpha \beta \delta }\eta ^{\delta \gamma }E_{\beta \gamma }.$$
(107)
For general $`V(u)`$, which may contain a nilpotent part, we have $`\mu ^T\eta +\eta \mu =0`$. Define again $`𝒴(u,\lambda )`$ by (99) then also

$$(u,\lambda )=M^1(u)M(u,\lambda )S=I+O(\lambda )$$
(108)
such that $`𝒴(u,\lambda )=(u,\lambda )S^1g^1(\lambda )S`$ and

$$(u,\lambda )^T\eta (u,\lambda )=\eta .$$
(109)
From equation (103) one obtains :

$$g^1\frac{g}{\lambda }=S^1(\frac{}{\lambda }𝒰\frac{\mu }{\lambda })S^1$$
(110)
Assume, that $`g=g_{}g_+`$ then

$$g_{}^1\frac{g_{}}{\lambda }+\frac{g_+}{\lambda }g_+^1=g_+S^1(\frac{}{\lambda }𝒰\frac{\mu }{\lambda })S^1g_+^1$$
(111)
Define $`G=g_+S^1S^1`$ then

$$g_{}^1\frac{g_{}}{\lambda }=\frac{1}{\lambda }G_0S\mu S^1G_0^1=\frac{1}{\lambda }G_0𝒱G_0^1$$
(112)
with matrix $`G_0=(g_+S^1S^1)_0=(g_+)_0`$ being orthogonal for a fixed point $`g`$ of the automorphism $`\sigma `$.

Acknowledgements
H.A. is partially supported by NSF (PHY-9820663).


text:
# Study of the SU(2) High Density System with Wilson Fermion

## 1 Introduction

High density state of the strongly interacting matter is a very attracting problem in the connection to the recent experiments of the Relativistic Heavy Ion Collider at BNL . Physics of extremely high density state is also important for the understanding the astrophysical objects such as neutron star. Usually in the statistical physics, high density state is investigated with use of the ground canonical ensemble with chemical potential which is introduced as a Lagrange Multiplier for the conserved charge. However, as is well known, chemical potential makes action of the QCD complex and causes a problem in the numerical convergence of simulation based on the monte calro method .

SU(2) gauge theory is the simplest non-abelian gauge theory and holds the pseud-real property which keeps action real even with chemical potential. Though the real QCD is not SU(2) but SU(3), SU(2) gauge theory is expected to play an role as a QCD-like theory which informs us important features of the non-abelian field theory . SU(2) gauge theory as a residual interaction in the Color Super Conducting state is also one of the promising possibility and the high density state of SU(2) gauge theory can be realized .

We calculate thermodynamical quantities as functions of chemical potential. Ground canonical ensemble states are generated through update based on the Woodbery formula . Iwasaki improved action is adopted for the gauge action. Though the improved action seems to stabilize the simulation better than the naive Wilson action , our simulation still shows instability at around $`\mu =1.0`$. However, the susceptibility of the Polyakov line $`\frac{L}{\mu }`$ shows peak at around $`\mu =0.7`$ which suggests the deconfinement phase transition. Therefore, our simulation can correspond to the state around across the phase transition.

## 2 Thermodynamical quantities

With use of Wilson fermions and Iwasaki improved action we calculate thermodynamical quantities, such as Polyakov line, gluon energy density and baryon number density with $`4^3\times 8`$ lattice. The expectation value of the baryon number density is given by,

$$<n>=\frac{1}{\beta V_s}\frac{}{\mu }\mathrm{log}Z$$
(1)
where $`V_s`$ is the spatial volume $`N_x\times N_y\times N_z`$. Energy density, $`\epsilon `$, is given as,

$$\epsilon =\frac{1}{V_s}\left(\frac{}{\beta }+\frac{\mu }{\beta }\frac{}{\mu }\right)\mathrm{log}Z$$
(2)
The derivative of the partition function is composed of two parts,

$$(\mathrm{log}Z)^{}=\frac{1}{Z}𝒟U𝒟\overline{\psi }𝒟\psi (S_G^{}S_F^{})e^{(S_GS_F)},$$
(3)
where $`S_G`$ and $`S_F`$ are gluon action and fermion action, respectively. We denote the contribution of the gluon action part (first term in the r.h.s. of eq.(3)) by gluon energy density. Figure 1 displays Polyakov line, gluon energy density and baryon number density as functions of $`\mu `$ where $`\beta =0.7`$ and $`\kappa =0.150`$, respectively. All quantities start to have non-zero value at about $`\mu =0.4`$ and rise up with chemical potential in the region $`0.4<\mu <0.8`$. None of them is order parameter of the phase transition in the exact sense, and since the lattice size is small, no sharp change is seen. However, growing up of the these quantities indicates that quarks and gluons become free from the confinement force at finite chemical potential about $`0.4<\mu <0.8`$.

Figure 2 shows susceptibility of the Polyakov line as the functions of the chemical potential $`\mu `$. The peak appears at around $`\mu =0.7`$ which suggests that our simulation crosses the deconfinement phase transition. However, peak point seems to depend on the hopping parameter $`\kappa `$. Larger hopping parameter seems to give lower phase transition $`\mu _C`$.

## 3 Propagators

Figure 3 shows our preliminary results of the propagator of mesons and diquarks with $`\kappa =0.175`$. In the lower figure, at $`\mu =0.0`$, pion and scalar diquark(b5b), and a0 and pseudscalar diquark(b1b) are degenerate each other, as expected. However, such degeneration of meson and diquark disappears with finite chemical potential. Because diquark possesses definite baryon number, the effect of the chemical potential on diquark and anti-diquark are in the opposite sign. Therefore, the symmetry of the propagator in time direction and anti-time direction is broken and asymmetry in nt appears(b1b and b5b in the lower figure of Fig.3). On the other hand, mesons, which have vanishing baryon number, keep symmetric property in nt direction . According to figure 2, $`\mu =0.8`$ with $`\kappa =0.175`$ is in the region upper side than the peak point.

## 4 Concluding Remarks

SU(2) finite density state is investigated with Wilson fermion and Iwasaki improved action. Peak of the susceptibility appears at around $`\mu =0.7`$, but it seems to depend on the $`\kappa `$. It means physical size of the lattice spacing will change with $`\kappa `$ and careful treatment will be required in the procedure of the chiral limit. The propagators of the mesons and diquark state (Baryon in SU(2)) have also discussed . Detailed analyses for both meson and diquark will appear somewhere else.

Acknowledgment

One of the authors(C.N.) would like to acknowledge the financial support by the Soryushi Shogakukai. This work is supported by Grant-in-Aide for Scientific Research by Monbu-Kagaku-sho, Japan (No.11440080 and No. 12554008). Simulations were performed on SR8000 at IMC, Hiroshima University, SX5 at RCNP, Osaka university, SR8000 at KEK and VPP5000 at Science Information Processing Center, Tsukuba university.


text:
# 1 Introduction

## 1 Introduction

Observations show significant abundance fluctuations along and between lines of sight in the interstellar medium (ISM), suggesting the ISM is far from homogeneous. For example, observations using the Interstellar Medium Absorption Profile Spectrograph (IMAPS) show local variations in the D/H ratio of about a factor of three (Jenkins et al. 1999; Sonneborn et al. 2000). Deuterium is mainly produced during Big Bang nucleosynthesis111However, Mullan & Linsky (1999) have proposed that D can be formed in stellar flares from M dwarf stars. and only destroyed in pre-main sequence stars. Therefore, the D in the ISM traces the fraction of gas not processed through stars, and the detection of local variations in the D/H ratio indicates incomplete mixing of processed material into the ISM. However, such an inhomogeneity picture contradicts the almost homogeneous distribution seen in O/H, which suggest that the ISM mixes completely on timescales no longer than that for chemistry reach an equilibrium.

These observations pose the question of how chemical elements such as D can avoid being mixed well enough to suppress local variations in the D/H ratio. Furthermore, the ISM is regulated by supernova explosions that form well structured and explosive flows in addition to diffuse turbulence that acts on the smaller scales. It is therefore, necessary to understand the mixing process in such a medium by carrying out direct simulations of the evolution of the ISM polluted with inhomogeneities of different length scales.

Section 2 deals with the description of the three-dimensional model used in this study, as well as, with the setup and evolution of the simulations. In section 3 we present the results of these simulations and discuss them in the light of the classical diffusion theory. Finally, in section 4 we present a summary of the main results.

## 2 Supernova-Driven ISM Model

### 2.1 Mixing Simulations

The model includes a fixed gravitational field provided by the stars in the disk, and an ideal-gas equation of state. Radiative cooling is treated assuming collisional ionization equilibrium, using a piecewise cooling function following that presented in Figure 2 of Dalgarno & McCray (1972) with an ionization fraction of 0.1 at temperatures below $`10^4`$ K and a temperature cutoff at 10 K. Background heating due to starlight varies with $`z`$ as described in Wolfire et al. (1995), and at $`z=0`$ it is chosen to balance radiative cooling at 8000 K.

Type Ib, Ic, and II SNe are set up at the beginning of their Sedov phases, with radii determined by their progenitor masses, as described in detail by Avillez (2000). Sixty percent of the SNe are set up within associations, and the rest are set up at random sites. OB associations are set up in a layer with a scale height of 46 pc (from the midplane) following the distribution of the molecular gas in the Galaxy, while the isolated SNe are set with an exponential distribution with a scale height of 90 pc. The rates of occurrence of SNe in the Galaxy are normalized to the volume under study, with the Galactic rate for types Ib and Ic in the Galaxy taken to be $`2\times 10^3`$ yr-1, while for type II it is $`1.2\times 10^2`$ yr-1 (Cappellaro et al. 1997). The total rate of these SNe in the Galaxy is $`\tau _{\mathrm{gal}}=1.4\times 10^2`$ yr-1, or one SN every 71 yr.

The interstellar gas is initially distributed in a smooth disk with the vertical distribution of Dickey & Lockman (1990). In addition, an exponential profile representing the $`z`$distribution of the warm ionized gas with a scale-height of 1 kpc in the Galaxy as described in Reynolds (1987) is used.

We model inhomogeneities in the ISM with a tracer field distributed in a checkerboard pattern of empty and full squares. We followed the evolution of models with squares having length $`L`$ of 25, 50 and 500 pc in the $`x`$ and $`y`$directions (parallel to the Galactic plane). Slab symmetry of the tracer field is assumed along the $`z`$-direction (perpendicular to the Galactic plane). In the current work we report on simulations using five SN rates: $`\sigma /\sigma _{\mathrm{gal}}=1`$, 10, 20, 30 and 50. The simulations are run until 400 Myr for $`\sigma /\sigma _{\mathrm{gal}}=1`$ and 200 Myr for the remaining supernova rates. This evolution time is enough for the system to reach steady state. A summary of these runs is presented in Table 1.

The computational domain contains a section of the Galaxy with an area of 1 kpc2 and vertical extension from -10 to 10 kpc. The innermost edge lies 8.5 kpc from the Galactic centre. The computational grid has a resolution of 10 pc, except in the layer between -250 and 250 pc, where three levels of adaptive mesh refinement are used, yielding a finest resolution of 1.25 pc. These simulations use the piecewise-parabolic method of Colella & Woodward (1984), a third-order scheme implemented in a dimensionally-split (Strange 1968) manner in combination with the AMR scheme of Berger & Colella (1989), with the subgridding scheme of Bell et al. (1994).

### 2.2 Global Evolution of the ISM

Once disrupted by the SN explosions, which start occurring at time $`t=0`$, the disk never returns to its initial state. Instead, regardless of the initial vertical distribution of the disk gas, a thin disk of cold gas forms in the Galactic plane, and, above and below, a thick, inhomogeneous gas disk forms. The code does not explicitly follow ionization states, but we can trace gas with temperature $`T10^4`$ K and scale height of 180 pc, which we designate Hi, and gas with $`10^4T10^5`$ K and scale height of 1 kpc, which we designate Hii. These distributions reproduce those described in Dickey & Lockman (1990) and Reynolds (1987), respectively. The upper parts of the thick Hii disk form the disk-halo interface, where a large scale fountain is set up by hot ionized gas escaping in a turbulent convective flow (see Avilez 2000 and Avillez & Mac Low 2001).

## 3 Evolution of ISM Inhomogeneities

The tracer field is a scalar whose mass flux follows that of the density and whose value is unity in the regions of the flow where it is present and zero elsewhere in the flow. By definition the tracer field should keep its value of 1, and therefore its average value $`<C>`$ on a checkerboard should always be 0.5 regardless of how it is advected. Numerical schemes introduce numerical diffusion leading sooner or later to a change of the value of the tracer field. However, its average value over the grid should remain constant. The simulations show that as the disk gas evolves to steady state $`<C>`$ approaches 0.5, with variations of no more than 2.4% occurring for a supernova rate of 30 times that in the Galaxy.

We define complete mixing to be when the maximum and minimum of the tracer field become numerically indistinguishable from each other and $`<C>`$, i.e., $`C_{\mathrm{max}}=C_{\mathrm{min}}=<C>0.5`$ ($`\mathrm{\Delta }C=C_{\mathrm{max}}C_{\mathrm{min}}<0.01`$). Table 1 presents the time scales for complete mixing in checkerboards with 50 and 500 pc scale-lengths for resolutions of 5, 2.5 and 1.25 pc, while Figures 1 and 2 show the time variation of $`C_{\mathrm{max}}`$ for a resolution of 5 pc and the variation of the mixing time with supernova rate for different finer resolutions, respectively.

Figure 1 shows that there is a steady decrease of $`C_{\mathrm{max}}`$ until complete mixing occurs. Futhermore, $`C_{\mathrm{max}}`$ and $`C_{\mathrm{min}}`$ take almost the same amount of time to vary from a difference $`\mathrm{\Delta }C=1`$ to $`0.05`$ as the time from $`\mathrm{\Delta }C=0.05`$ to $`0.01`$, when we define complete mixing to occur. This indicates that smaller scale structures dominate the mixing process after the gas has been completely mixed on larger scales.

Figures 1 and 2 show that the time taken for complete mixing to occur decreases with increasing supernova rate, that is, increased energy input into the ISM. However, the rate of increase of mixing reaches a saturation level for $`\sigma /\sigma _{Galaxy}>20`$, beyond which the increase in the supernova rate does not accelerate mixing in the ISM. For SN rates greater than 20 times the Galactic rate we see that the mixing time scales differ by only a few Myr (see Table 1), with complete mixing occurring at $`107`$ Myr for $`\mathrm{\Delta }x=5`$ pc. This results from the fact that for $`\sigma /\sigma _{Galaxy}20`$, after the first 50 Myr of evolution a large fraction fo the disk gas has a temperature greater than $`10^5`$ K and at 100 Myr of evolution all the gas is hotter than $`10^5`$ K. As a consequence of the high SN rate the disk becomes hotter and cooling becomes inefficient. At this high temperature the sound speed is correspondingly high, so the gas mixes within a 40-50 Myr period (Roy & Kunth 1995) leading to complete mixing at some 110 Myr.

The time taken for complete mixing to occur in a checkerboard with 50 pc squares varies from $`300`$ Myr for $`\sigma /\sigma _{Galaxy}=1`$ to some 120 Myr for $`\sigma /\sigma _{Galaxy}=20`$ and 107 for $`\sigma /\sigma _{Galaxy}>30`$ for $`\mathrm{\Delta }x=5`$ pc, while for $`\mathrm{\Delta }x=1.25`$ pc mixing varies from 347 Myr for $`\sigma /\sigma _{Galaxy}=1`$ to 150 Myr for $`\sigma /\sigma _{Galaxy}=20`$. On 500 pc squares complete mixing is surprisingly only delayed by around 20 Myr (Table 1). The similarity between these time scales indicates that the mixing process is almost independent of the length-scale of the inhomogeneities, and therefore, of the initial setup adopted for the tracer field. Rather, it is regulated by supernovae and their energy input into the ISM.

According to classical mixing length theory the mixing time scales of these boards should follow the relation $`\tau \mathrm{}^2`$, where $`\mathrm{}`$ is the length of the squares. Therefore, the mixing time of inhomogeneities with a scale length of 500 pc should be 25 times longer than the mixing time for 50 pc squares. The computed similarity in the mixing time scales (Table 1) for the different inhomogeneity sizes shows that classical mixing length theory does not describe mixing in the ISM well. This seems to be because much of the ISM is dominated by nearly laminar flows inside of SNRs and superbubbles rather than by diffuse turbulence. The latter does act on smaller scales and in cooler gas, in which classical mixing length theory may apply better.

## 4 Summary

In this paper we try to shed light on the time scales that it takes to erase inhomogeneities in the ISM and how this process depends on the inhomogeneity length scale. We find: (1) The time scale to erase inhomogeneities in the ISM is nearly independent of their length scale on scales above 10 pc or so; (2) Mixing time scales decrease with increasing supernova rate until all the gas is hot and the mixing rate saturates. This occurs when $`\sigma /\sigma _{\text{Galaxy}}20`$, when the mixing time scale is $`107120`$ Myr for $`\mathrm{\Delta }x=5`$ pc and 150 Myr for $`\mathrm{\Delta }x=1.25`$ pc; (3) Even after some $`80\%`$ of the gas has been mixed, the simulations show that the ISM remains clumpy for at least twice as long if we examine scales smaller than a kpc. Even if the rate of SNe is increased to ten times the Galactic rate, this clumpy structure does not disappear over a time of a few tens of Myr.

This work shows that inhomogeneities that are present in the ISM, assuming that no further inhomogeneities are introduced into the system, take a very long time to be erased. The time for complete mixing (greater than 107 Myr, which occurs for $`\sigma /\sigma _{\text{Galaxy}}=50`$ and $`\mathrm{\Delta }x=5`$ pc) is already longer that any time scale for chemical equilibrium to occur, and therefore, the ISM will show a clumpy distribution. The local variations occur because, even in the presence of numerical diffusion, gas does not mix quickly between hot and cold regions.

Acknowledgements.  The authors acknowledge useful discussions with D. York, and E. Jenkins, and with R. Ferlet, who also emphasized the importance of working on this problem. This work was supported by an NSF CAREER grant (AST 99-85392).


text:
# REFERENCES

Superconducting energy gaps, low temperature specific heat, and quasiparticle spectra of MgB2
Hyoung Joon Choi, David Roundy,∗† Hong Sun, Marvin L. Cohen,∗† & Steven G. Louie∗†
Department of Physics, University of California at Berkeley, Berkeley, CA 94720, USA.
Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA.

Magnesium diboride is remarkable not only for its unusually high transition temperature of 39K for a $`sp`$-bonded metal and its reduced isotope-effect exponent, but also for its unusual superconducting gap behavior. Some earlier experiments seem consistent with conventional single-gap superconductivity, but recent measurements indicate the existence of more than one gap. For example, specific heat measurements show an anomalous large bump at low temperature, inconsistent with single-gap BCS theory. Moreover, there is a significant variation in the measured values for the gaps. Here, we report first-principles calculations of the $`k`$\- and $`T`$-dependent superconducting gap $`\mathrm{\Delta }(\stackrel{}{k},T)`$ in MgB2 and its manifestation in various measured quantities. Because its Fermi surface has disconnected sheets with very different electron-phonon coupling strengths, our calculations show that near $`T=0`$, the values of $`\mathrm{\Delta }(\stackrel{}{k})`$ cluster into two groups: large values ($``$ 6.5 to 7.5 meV) on the strongly coupled sheets and small values ($``$ 1 to 3 meV) on the weakly coupled sheets. The calculated gap, quasiparticle density of states, and specific heat and their temperature dependences are in agreement with the recent measurements which support that MgB2 is a multiple gap superconductor. In fact, theory predicts four prominent values for the gap at low $`T`$.

Most theoretical work on MgB2 has been focused on the electronic structure, the phonon structure, and the electron-phonon interaction. It is shown that the electron-phonon interaction in MgB2 varies strongly on the Fermi surface. Also, recent two-band model and first-principles studies of the transition temperature support a multigap scenario. However, no detailed quantitative calculations have been presented for the temperature and $`k`$-space dependence of the superconducting gap. The calculational framework used here is based on the strong coupling formalism of superconductivity established by Eliashberg. The Eliashberg formalism is a more general case of the original formulation of the BCS theory for phonon-mediated pairing, taking into account strong electron-phonon coupling. When applying the strong coupling theory to MgB2, we include the momentum dependency of the anisotropic Eliashberg equation and also take into account the anharmonic effect on phonon frequencies. Material-dependent parameters in the equation are obtained by ab initio pseudopotential density functional calculations. This approach has already been shown to describe the superconducting transition of MgB2 very accurately; the obtained specific heat mass enhancement factor $`\lambda =0.61`$, the transition temperature $`T_c=39K`$, and isotope-effect exponent $`\alpha _B=0.31`$ are all in excellent agreement with experiments. Here, we compute the gap function $`\mathrm{\Delta }(\stackrel{}{k},\omega )`$ and the renormalization factor $`Z(\stackrel{}{k},\omega )`$ from the fully anisotropic Eliashberg equation at various temperatures. This allows us to address all the superconducting properties except those related to the presence of an applied magnetic field.

Figure 1 shows the calculated superconducting gap $`\mathrm{\Delta }(\stackrel{}{k})`$ on the Fermi surface at 4 K. The Fermi surface of MgB2 consists of four sheets: two-dimensional (2D) light hole and heavy hole sheets forming coaxial cylinders along $`\mathrm{\Gamma }`$ to $`A`$, a three-dimensional (3D) hole sheet connecting regions near $`K`$ and $`M`$, and a 3D electron sheet connecting regions near $`H`$ and $`L`$. The 2D light and heavy hole sheets are derived from $`\sigma `$-antibonding states of boron $`p_{x,y}`$ orbitals, while the 3D hole and electron sheets are derived from $`\pi `$-bonding and antibonding states of boron $`p_z`$ orbitals, respectively. The calculated density of states at the Fermi energy is 0.115 states/eV$``$atom$``$spin; 44 % of which comes from the 2D $`\sigma `$ cylindrical sheets and the rest comes from the 3D $`\pi `$ sheets. The superconducting gap is nonzero everywhere on the Fermi surface. Thus, the symmetry of the superconducting gap is s-wave, but the size of the gap changes on the Fermi surface. The gap values cluster into two groups. The largest gap is on the 2D light hole $`\sigma `$ cylindrical sheet (shown red in Fig. 1), where the average gap is 7.2 meV with variations of less than 0.1 meV. On the 2D heavy hole $`\sigma `$ cylindrical sheet (shown orange in Fig. 1), the superconducting gap ranges from 6.4 meV to 6.8 meV, having an average of 6.6 meV, with maximum value near $`\mathrm{\Gamma }`$ and minimum value near A. The average of the gap values on the two 2D $`\sigma `$ cylindrical sheets is 6.8 meV. The superconducting gap is significantly smaller and more spread out on the 3D $`\pi `$ sheets (shown green and blue in Fig. 1), ranging from 1.2 meV to 3.7 meV. Out of the two 3D $`\pi `$ sheets, the hole sheet shows a slightly larger gap with a peak at 2.1 meV and the electron sheet shows a somewhat smaller gap with a peak at 1.5 meV. The average of the gap values on the two 3D $`\pi `$ sheets is 1.8 meV. This very large variation in the gap value originates from the differences in electron-phonon coupling strengths on the different sheets of the Fermi surface. Our ab initio calculations demonstrate that the electron-phonon coupling strength is 4 to 5 times stronger on the cylindrical $`\sigma `$ sheets (exceeding 2.5 for some states) than on the $`\pi `$ sheets. Our result is consistent with the recent experiments reporting two gaps because the superconducting gap in MgB2 can be grossly viewed as consisting of a large gap of $``$ 6.8 meV on the strongly coupled 2D $`\sigma `$-cylindrical sheets and a small gap of $``$ 1.8 meV on the weakly coupled 3D $`\pi `$ sheets. However, this simplication is limited. As seen in Fig. 1, the large and small gaps both have substantial variation in value: 6.4 to 7.2 meV on the 2D $`\sigma `$ sheets and 1.2 to 3.7 meV on the 3D $`\pi `$ sheets. The experimental interpreted gaps have a range of 1.5 to 3.5 meV for the small gap and a range of 5.5 to 8 meV for the large gap.

Figure 2 shows the temperature dependence of the superconducting energy gap on the Fermi surface from 4K to 38K which may be probed in tunneling, optical, and specific heat measurements. The vertical blue curves in Fig. 2 present the distribution of the superconducting gap values at various temperatures, and the red lines are curves of the form $`\mathrm{\Delta }(T)=\mathrm{\Delta }(0)\sqrt{1(T/T_c)^p}`$ fitted separately to the averaged values of the calculated gap on the 2D $`\sigma `$ sheets and those on the 3D $`\pi `$ sheets. Our calculated results give $`\mathrm{\Delta }(0)=6.8`$ meV ($`2\mathrm{\Delta }(0)/k_BT_c=4.0`$) for the averaged gap at $`T=0`$ on the 2D $`\sigma `$ sheets, and $`\mathrm{\Delta }(0)=1.8`$ meV ($`2\mathrm{\Delta }(0)/k_BT_c=1.06`$) for that on the 3D $`\pi `$ sheets. The theory shows that the superconducting gap opens at $`T_c`$ throughout the Fermi surface, and that the gap widens much faster on the 2D $`\sigma `$ sheets than on the 3D $`\pi `$ sheets at $`T{}_{}{}^{<}T_{c}^{}`$.

From the calculated gap function $`\mathrm{\Delta }(\stackrel{}{k},\omega )`$, the quasiparticle density of states is given by

$$N(\omega )/N(0)=\text{Re}\frac{\omega +i\mathrm{\Gamma }}{\sqrt{(\omega +i\mathrm{\Gamma })^2\mathrm{\Delta }(\stackrel{}{k},\omega )^2}},$$
(1)
where $`\mathrm{}`$ indicates an average over a surface of constant $`\omega `$. Figure 3 depicts the theoretical quasiparticle density of states calculated with an assumed finite lifetime $`\mathrm{\Gamma }`$ of 0.1 meV. The quasiparticle density of states at 4K shows three discernable peaks. One peak is at 2.2 meV, which is slightly greater than the dominant size of the gap (2.1 meV) on the 3D $`\pi `$ hole sheet. The other two peaks are at 6.8 and 7.2 meV, which are the maximal sizes of the superconducting gap on the 2D heavy and light hole $`\pi `$ sheets, respectively. The quasiparticle density of states can be deduced from tunneling experiments and various spectroscopic measurements, but a direct quantitative comparison requires knowledge of various physical parameters involved in a specific experiment.

A fundamental measurement which can provide a probe of the superconducting quasiparticle spectrum is the temperature-dependent specific heat. Figure 4 shows the calculated $`T`$-dependent specific heat of MgB2 through the superconducting transition. For the normal-state specific heat $`C_N=\gamma T`$, we obtain $`\gamma `$ = 2.62 mJ/mol$``$K2, which agrees well with experimental values 2.6 mJ/mol$``$K2 and 2.7 mJ/mol$``$K2. For the superconducting-state specific heat $`C_S`$, we first calculate the free energy difference ($`F_SF_N`$) of the superconducting and normal states and then obtain the specific heat difference by

$$C_SC_N=T\frac{d^2}{dT^2}(F_SF_N).$$
(2)
The measured low $`T`$ specific heat shows substantial magnitude and a large bump at about 10K. This anomalous behavior is inconsistent with a 1-gap BCS model. However, our results are in excellent agreement with experiment. The origin of the bump in our calculated curve is the existence of low energy excitations above the small gap on the weakly coupled 3D $`\pi `$ sheets of the Fermi surface. The overall shape and magnitude of the calculated specific heat curve agrees very well with experiments, especially below 30K, and shows a jump of $`\mathrm{\Delta }C/\gamma T_c=1.0`$ at $`T_c`$ which is within the range of the measured values.

This work was supported by National Science Foundation Grant No. DMR00-87088, and by the Director, Office of Science, Office of Basic Energy Sciences of the U. S. Department of Energy under Contract DE-AC03-76F0098. Computational resources have been provided by the National Science Foundation at the National Center for Supercomputing Applications and by the National Energy Research Scientific Computing Center. Authors also acknowledge financial support from the Miller Institute (H.J.C.) and from the Berkeley Scholar Program funded by the Tang Family Foundation (H.S.).

Correspondence should be addressed to S.G.L. (e-mail: sglouie@uclink.berkeley.edu).


text:
# Covariant model for proton-proton bremsstrahlung: comparison with high-precision data

## A Introduction

The proton-proton bremsstrahlung ($`pp\gamma `$) process has attracted significant attention over the years, both experimentally and theoretically. In recent years, several new experiments have been performed , inspiring many new theoretical investigations . In particular, a number of microscopic models have been developed to describe the $`pp\gamma `$ process. Examples are the potential model of Nakayama et al and the covariant model of Martinus et al. The theoretical predictions of these models could be compared to the $`pp\gamma `$ cross sections and analyzing powers for the TRIUMF experiment at 280 MeV . The agreement of theory with these TRIUMF data is rather good, especially for the cross sections, provided that the experimental cross sections are renormalized . Some outstanding discrepancies occur, however, for certain asymmetric proton angles.

More recently, the first high-precision data from the KVI experiment at 190 MeV became available , and many more data are forthcoming . When comparing these data with theory, a pronounced and undisputable discrepancy between theory and experiment was observed in specific kinematic regions . The size of the discrepancy between theory and experiment is disturbing, since what primarily enters are the two-nucleon ($`NN`$) interaction and the electromagnetic coupling of the photon to the $`NN`$ system, both of which are believed to be accurately known at this energy. The high precision of the new KVI data allows one, in principle, to study smaller effects, like those arising from negative-energy states, the $`\mathrm{\Delta }`$-isobar, and meson-exchange currents. It is therefore important to identify the possible reasons for the discrepancies.

In this paper, we compare the covariant model of Ref.  with the KVI data available so far and analyze the discrepancies. We demonstrate that the dominant contribution to $`pp\gamma `$ for the specific problematic kinematic regions results when the $`NN`$ interaction is evaluated at energies below about 15 MeV, and that, in fact, at least a major part of the problem resides in the low-energy behavior of the $`NN`$ interaction models used.

This paper is organized as follows. First, we review briefly the covariant model for $`pp\gamma `$ and compare it to the KVI data in the problematic kinematic regions. Next, we demonstrate the sensitivity of the bremsstrahlung cross sections in these regions to the properties of the low-energy $`NN`$ interaction. We then show that the problems can be significantly alleviated accordingly by improving the description of the interaction at low energies. We end by summarizing our conclusions.

## B Covariant model for bremsstrahlung

We first give a short review of the microscopic model of Martinus et al. for $`pp\gamma `$ . In this covariant model the $`NN`$ $`T`$-matrix is obtained by solving the Bethe-Salpeter equation for the two-nucleon system  in the equal-time approximation with the one-boson exchange (OBE) kernel of Ref. . This OBE model for the $`NN`$ interaction was, at that time, fitted to the Virginia Tech $`np`$ partial-wave solution , by adjusting the meson-nucleon coupling constants and the form factor parameters, and a reasonable agreement was obtained.

This covariant $`NN`$ $`T`$-matrix enters the model for $`pp\gamma `$, in which a number of contributions can be distinguished. The most important ones in the energy regime we consider here are the “nucleonic” contributions, consisting of single-scattering terms, i.e. photon emission off the external proton legs (see Fig. 1, a-b), and the contribution commonly known as rescattering (see Fig. 1, c). The model is relativistic covariant and therefore negative-energy states are included in a natural way. The relevance of these negative-energy states is small for energies around 200 MeV. The reason  is that the contributions from the single-scattering diagrams where the intermediate nucleons are in a negative-energy state are canceled by similar contributions coming from the rescattering diagram. This cancellation holds for terms up to order $`𝒪(q)`$, where $`q`$ is the photon momentum, and is a consequence of the soft-photon theorem .

Contributions from the $`\mathrm{\Delta }`$-isobar and from magnetic meson-exchange currents, containing in particular the $`\omega \pi \gamma `$ and $`\rho \pi \gamma `$ decay graphs, are also taken into account. These two-body current terms are included in a perturbative way, since they are small in general. The coupling constants of the photon to the various mesons were determined from the radiative decay widths of the vector mesons. As one can observe from Fig. 2 below, at an energy of 190 MeV the contribution of the two-body currents is small. These terms, however, increase in size with energy and can become appreciable around the pion-production threshold at 280 MeV and above.

This covariant bremsstrahlung model  is theoretically well founded and many of its ingredients have been tested in other calculations such as those for electron scattering on the deuteron . Therefore, one did not expect major discrepancies with new experimental $`pp\gamma `$ data at energies below 280 MeV.

## C Comparison with the KVI data

In the KVI experiment, $`pp\gamma `$ cross sections and analyzing powers were measured for 190 MeV incoming proton energy, with the scattered protons detected at small forward angles, and with the photon emitted in the backward hemisphere . A typical example of the KVI data  in comparison with theory is shown is Fig. 2. The data are plotted for two different kinematic situations, with $`\theta _1`$ and $`\theta _2`$ the fixed angles of the outgoing protons in the laboratory frame, and as function of the polar angle of emission of the photon, $`\theta _\gamma `$. The theoretical predictions are from the model as published in Ref. .

For the asymmetric proton angles $`\theta _1=8^{}`$, $`\theta _2=16^{}`$ (upper left panel) the cross section shows a large discrepancy between theory and experiment for values of $`\theta _\gamma `$ corresponding to the backward peak in the cross section. For the symmetric proton angles $`\theta _1=\theta _2=16^{}`$ (upper right panel), on the other hand, the cross section shows a much better agreement between theory and data. The contribution of the two-body currents is seen to be minor, and thus the discrepancy for $`\theta _1=8^{}`$, $`\theta _2=16^{}`$ is unlikely to come from this source. In the following, we focus therefore on the nucleonic contribution as the cause of the problem.

Taking a closer look at the two kinematics presented in Fig. 2 reveals that the one which poses problems ($`\theta _1=8^{}`$, $`\theta _2=16^{}`$) is dominated by the contribution from the $`{}_{}{}^{1}S_{0}^{}`$ wave in the $`NN`$ $`T`$-matrix, while for the second one for which the agreement is much better ($`\theta _1=\theta _2=16^{}`$) the $`{}_{}{}^{3}P`$ waves are as important as the $`{}_{}{}^{1}S_{0}^{}`$ wave. The contribution of the rescattering diagram is relatively small, suggesting that the problem resides already at the level of the single-scattering diagrams.

In the bremsstrahlung calculation the $`NN`$ $`T`$-matrix is evaluated at different energies for each of the diagrams shown in Fig. 1. The energy is lowest for the diagrams where the photon is emitted by one of the incoming protons. The value of this energy is shown in Fig. 3 as a function of $`\theta _\gamma `$. For $`\theta _1=8^{}`$, $`\theta _2=16^{}`$, the $`T`$-matrix is evaluated at 7 and 11 MeV kinetic energy for the cases corresponding to the peaks in the bremsstrahlung cross sections, i.e. at about $`\theta _\gamma =20^{}`$ and $`140^{}`$, respectively. For the minimum in the cross section, around $`\theta _\gamma =75^{}`$, the $`T`$-matrix is evaluated at 24 MeV when the photon is emitted by one of the incoming protons. In contrast, one observes that for $`\theta _1=\theta _2=16^{}`$, the $`T`$-matrix is evaluated at more than 20 MeV for $`\theta _\gamma =160^{}`$ and at more than 40 MeV for the minimum around $`\theta _\gamma =75^{}`$; see the upper panel of Fig. 3. Thus, we conclude that theory and experiment agree well for those situations where the $`T`$-matrix is evaluated above about 15 MeV, while the discrepancies occur for the cases corresponding to a lower energy.

In order to demonstrate the sensitivity of the $`pp\gamma `$ cross section to the low-energy $`NN`$ interaction, we plot in the lower panel of Fig. 3 the cross section for $`\theta _1=8^{}`$, $`\theta _2=16^{}`$, where we only include the $`{}_{}{}^{1}S_{0}^{}`$ wave, for three values of the scattering length. Here we simply changed, for illustrative purpose only, the value of $`g_\epsilon ^2`$, the coupling constant of the “$`\epsilon `$-meson” in the OBE model.

In Table 1 we list for the two cases, $`\theta _1=8^{}`$, $`\theta _2=16^{}`$ and $`\theta _1=\theta _2=16^{}`$, the cross sections calculated by including in the $`NN`$ $`T`$-matrix only one of the important partial waves, viz. the $`{}_{}{}^{1}S_{0}^{}`$ wave and the $`{}_{}{}^{3}P`$ waves. This is done for the value of $`\theta _\gamma `$ corresponding to the backward peak or plateau in the cross section, i.e. $`\theta _\gamma =140^{}`$ in the first case and $`\theta _\gamma =160^{}`$ in the second case, respectively. In both cases the first number listed gives the cross sections with the $`T`$-matrix as used in Ref.  (the second number gives the corresponding value in the fit described below). It is seen that the contribution of the $`{}_{}{}^{1}S_{0}^{}`$ wave is significantly larger for the case $`\theta _1=8^{}`$, $`\theta _2=16^{}`$ compared to the case $`\theta _1=\theta _2=16^{}`$. By far the largest contribution to this $`{}_{}{}^{1}S_{0}^{}`$ cross section comes from radiation off the incoming proton legs. The $`{}_{}{}^{3}P`$ wave cross sections arise mainly from radiation off the outgoing proton legs. For $`\theta _1=8^{}`$, $`\theta _2=16^{}`$ emission from the incoming proton legs dominates the bremsstrahlung cross section, while for $`\theta _1=\theta _2=16^{}`$ emission from initial and final proton legs are comparable. The total cross sections listed in the bottom row of Table 1 include also the rescattering contribution.

At the time of its construction, the OBE model of Ref.  was fitted to the $`np`$ phase-shift parameters of the Virginia Tech group, where the interest of the fit was at the higher energies. In view of our findings with the new high-precision $`pp\gamma `$ data, it becomes important to improve the description at low energies. In order to investigate whether at least the major part of the discrepancy between theory and the KVI data can be removed, we have performed a preliminary refit with an emphasis of obtaining better agreement with the Nijmegen $`pp`$ partial-wave analysis of Ref. , at energies starting from 10 MeV up to 215 MeV.

In Fig. 4 we present the same calculation as in Fig. 2, but now for the $`NN`$ model after the preliminary new fit. The description of the $`pp\gamma `$ cross sections in the region of the backward peak for $`\theta _1=8^{}`$, $`\theta _2=16^{}`$ has improved substantially. This is mainly due to the improvement in the description of the $`{}_{}{}^{1}S_{0}^{}`$ wave at low energies. Also the description of the $`{}_{}{}^{3}P`$ partial waves has improved. In Table 1 we now again list for the two kinematical situations the contribution to the $`pp\gamma `$ cross section of the different important partial waves. The new results should be compared to the values obtained with the model as published in Ref. , and also listed in Table 1.

A preliminary investigation suggests that for those regions in phase space where similar problems occurred as demonstrated here for the case $`\theta _1=8^{}`$, $`\theta _2=16^{}`$, at least a large part of the problem is due to an inadequate description of the low-energy behavior of the $`NN`$ $`T`$-matrix. Detailed results will be presented in Ref. .

Finally, it is worth pointing out that the sensitivity of the $`pp\gamma `$ cross sections to the low-energy $`NN`$ interaction raises the question about the comparison between microscopic bremsstrahlung models of the type discussed here and soft-photon descriptions. The soft-photon amplitude of Ref.  gives a remarkably good description of the experimental $`pp\gamma `$ cross sections, including the data from KVI. In constructing the soft-photon amplitude, a specific choice is made for the on-shell points used in the calculation. For the kinematic situations discussed here, this choice does not correspond to the $`NN`$ interaction at low energies, and therefore, the relation to microscopic models is not obvious.

## D Conclusions

In conclusion, the $`pp\gamma `$ cross section at 190 MeV varies strongly as a function of the angles of the protons. This can be traced back to an increasing contribution of the $`{}_{}{}^{1}S_{0}^{}`$ wave for the regions in phase space where the $`NN`$ interaction is probed at low energies. For $`\theta _1=8^{}`$, $`\theta _2=16^{}`$, the $`{}_{}{}^{1}S_{0}^{}`$ wave dominates strongly at the peak of the cross section. In contrast, for $`\theta _1=\theta _2=16^{}`$, the $`{}_{}{}^{3}P`$ waves are as important as the $`{}_{}{}^{1}S_{0}^{}`$ wave.

The present study shows that the predicted cross sections in certain kinematic regions are sensitive to the $`NN`$ interaction at low energy. In all cases, the analyzing powers are affected less by changes in the interaction. We have in particular not included in these calculations the Coulomb interaction, which plays an important role in the $`pp`$ system at low energies. The results of such an analysis will be given in Ref. . Only a limited set of KVI data has been published so far . With the complete $`pp\gamma `$ data set , one may hope to be able to test in a quantitative way the validity of the various microscopic bremsstrahlung models.

## Acknowledgments

We acknowledge discussions with our experimental colleagues at KVI. The research of R.G.E.T. was made possible by a fellowship of the Royal Netherlands Academy of Arts and Sciences. We thank M.C.M. Rentmeester for his help in setting up the codes for fitting the $`NN`$ model.


text:
# Double parton scattering processes in 𝑝𝑝 collisions and the scale factors.

## 1 Double Parton Collisions

Multiple parton interaction processes, where different pairs of partons have an hard scattering in the same hadronic collision, become experimentally feasible at high energies because of the growing flux of partons. The hard partonic interactions occur independently in different points inside the overlap region of the two colliding hadrons and are separated by a distance of the order of the hadron radius. These processes are characterized by two very different scales: the hadronic dimension and the inverse of the hard scale in the partonic subprocesses. Therefore the interference effects betweeen the partonic interactions can be neglected and the expression of the cross sections for these processes have a simple probabilistic interpretation.
In the case of double parton collisions, for the process $`ppAB`$, the cross section is

$`\sigma _D(AB)`$ $`=`$ $`{\displaystyle \frac{m}{2}}{\displaystyle 𝑑x_1𝑑x_1^{}𝑑x_2𝑑x_2^{}d^2b}`$ (1)
$`\times `$ $`\mathrm{\Gamma }(x_1,x_2;b)\widehat{\sigma }^A(x_1,x_1^{})\widehat{\sigma }^B(x_2,x_2^{})\mathrm{\Gamma }^{}(x_1^{},x_2^{};b)`$
The non-perturbative input in this expression is the two-body parton distribution $`\mathrm{\Gamma }(x_1,x_2;b)`$ where $`x_1`$ and $`x_2`$ are the momentum fractions and $`b`$ is the distance of the two partons in transverce space. $`\widehat{\sigma }^A`$ and $`\widehat{\sigma }^B`$ are the partonic cross sections for the production of the final states $`A`$ and $`B`$ respectively and $`m=1(m=2)`$ when $`A`$ and $`B`$ are indistinguishable (distinguishable) processes.

If one assumes that the the momentum fractions are not correlated and the dependence on $`b`$ can be factorized, the inclusive double parton scattering cross section for the two parton processes $`A`$ and $`B`$ reduces to the simplest expression

$$\sigma _D=m\frac{\sigma _A\sigma _B}{2\sigma _{eff}}$$
(2)
where $`\sigma _A`$ and $`\sigma _B`$ are the inclusive single scattering cross sections for producing the final states $`A`$ and $`B`$ respectively. The information on the structure of the hadron in transverse space is summarized in the value of the scale factor $`\sigma _{eff}`$ which is universal i.e. does not depend on the considered process. The experimental value quoted by CDF is $`\sigma _{eff}=14.5\pm 1.7_{2.3}^{+1.7}`$ mb .

Using this simple factorized expression one find that the double parton collisions give important effects in several processes at the LHC, for example in the production of $`(W,Z)b\overline{b}`$, $`W+jets`$, four $`jets`$, $`b\overline{b}b\overline{b}`$, etc . An important case is represented by the production of $`Wb\overline{b}`$ in a double parton scattering, which is a large source of background in the Higgs boson production mechanism $`ppW(Hb\overline{b})`$ .

## 2 Scale Factors

The scale factor $`\sigma _{eff}`$ is related to the transverse size of the region where the hard interactions take place. More precisely, this parameter is given by the overlap in transverse space of the matter distribution of the two interacting hadrons. The measured value of the scale factor is sizably smaller than the naive expectations and it is an indication of important correlations effects in the hadron structure . In fact, while correlations in the momentum fractions will affect the two-body parton distributions only for large values of $`x`$, the correlations in the transverse parton coordinates are always present due to the binding force. In a recent paper  the small value of the effective cross section was explained by a simple correlated model of the partonic population in the proton. In this model the distributions in transverse space of gluons and sea quarks are correlated with the configuration taken by the valence quarks, in such a way that extended configurations of the valence give rise to a large number of partons while compact configurations give rise to a smaller number, as one expects for the growth of the energy of the gluonic field with the distance between the valence quarks. On the other hand in the model the momentum fractions of the initial parton pairs are uncorrelated. One of the main features is that there are two different kinds of partons as far as the distributions in transverse space are concerned: the valence quarks and the sea quarks and gluons. As a consequence the scale factor in a reaction involving valence quarks will be different than the scale factor in a reaction involving sea quarks and gluons .

The inclusive double parton scattering cross section is therefore expressed as

$`\sigma _D(A,B)`$ $`=`$ $`{\displaystyle \frac{m}{2}}{\displaystyle \underset{ijkl}{}}{\displaystyle 𝑑x_1𝑑x_1^{}𝑑x_2𝑑x_2^{}d^2b}`$ (3)
$`\times `$ $`\mathrm{\Gamma }_{ij}(x_1,x_2;b)\widehat{\sigma }_{ik}^A(x_1,x_1^{})\widehat{\sigma }_{jl}^B(x_2,x_2^{})\mathrm{\Gamma }_{kl}^{}(x_1^{},x_2^{};b)`$
where the indices $`i`$ and $`j`$ refer to the different kinds of partons, while dependence of $`\mathrm{\Gamma }`$ on $`x_1`$, $`x_2`$ and $`b`$ is factorized: $`\mathrm{\Gamma }_{ij}(x_1,x_2;b)=G_i(x_1)G_j(x_2)F_j^i(b)`$, and $`G_i(x)`$ are the usual one-body parton distributions.

The model introduces the geometrical coefficients

$$\mathrm{\Theta }_{kl}^{ij}=d^2bF_k^i(b)F_{l}^{j}{}_{}{}^{}(b),$$
(4)
wich have the dimension of the inverse of a cross section, so that the double scattering cross section is given by

$$\sigma _D(A,B)=\frac{m}{2}\underset{ijkl}{}\mathrm{\Theta }_{kl}^{ij}\sigma _{ij}(A)\sigma _{kl}(B)$$
(5)
where $`\sigma _{ij}(A)`$ is the inclusive cross section for the process $`ijA`$.

By using the equation

$$\underset{ijkl}{}\mathrm{\Theta }_{kl}^{ij}\sigma _{ij}(A)\sigma _{kl}(B)=\frac{\sigma (A)\sigma (B)}{\sigma _{eff}}$$
(6)
we can therefore obtain the scale factors, as a function of the c.m. energy, in various reactions where the rate of double parton scatterings is large enough to measure the corresponding scale factor. We have considered the following production processes: equal sign $`WW`$ , $`Wjj`$, $`Wb\overline{b}`$, 4-jets and $`b\overline{b}b\overline{b}`$.

Our results are plotted in Fig.1, and show strong difference between the scale factors in the channels with and without $`W`$ bosons, as a consequence of the different role played by the valence in the two cases. The scale factor in reactions involving valence quarks is in fact larger than that in reactions involving gluons and sea quarks. This is due to the fact that in the model the average transverse distance of a pair of valence quarks is larger than the average transverse distance of a pair of gluons or sea quarks. The overall qualitative indication of our calculations is that the scale factors are process and energy dependent and that they are the physical observables wich can provide a first quantitative information on the correlations in transverse space of the initial state partons.

## Acknowledgments

This work was partially supported by the Italian Ministry of University and of Scientific and Technological Research by means of the Fondi per la Ricerca scientifica-Università di Trieste.