text:
# High Accuracy Many-Body Calculational Approaches for Excitations in Molecules

## Abstract

Two state-of-the-art computational approaches: quantum Monte Carlo (QMC), based on accurate total energies, and GW with exciton effects (GW-BSE), based on perturbation theory are employed to calculate ionization potentials, electron affinities, and first excited singlet and triplet energies for the silane and methane molecules. Results are in excellent agreement between these dramatically different approaches and with available experiment. The optically forbidden triplet excitation in silane is predicted to lie roughly 1 eV higher than previously reported. For methane, the impact of geometry relaxation is shown to be $``$ 2 eV for excited states. Further, in the GW-BSE method, we demonstrate that inclusion of off-diagonal matrix elements in the self-energy operator is crucial for an accurate picture.

Optical excitations play a fundamental role in technological applications such as dye chemicals, photovoltaics, laser technology, and catalysts for chemical reactions. An accurate evaluation of excited states requires more sophisticated approaches than simply a straightforward application of the mainstream ground state methods such as density functional theory and Hartree-Fock. In addition, many traditional approaches for calculating excited states can be limited by a range of factors. For example, configuration interaction with single excitations is missing correlation effects and is often inaccurate, optical transitions modeled as quasiparticle energies corrected by excitonic effects rely on an approximate description of the electron-hole interaction, and basis set correlated wave function approaches are limited by their scaling with the number of correlated electrons.

Two promising new approaches have been developed to calculate accurate excitation energies: (1) quantum Monte Carlo (QMC), based on a stochastic solution to the many-electron Schrödinger equation, can provide chemical accuracy (i.e., 0.05 eV) for the total energy while at the same time scales as N3 where N is the number of valence electrons (see, e.g., ), and (2) many-body Green’s function perturbation techniques that employ the $`GW`$ approximation  for the electron self-energy operator to calculate the single-particle excitation spectrum, followed by a solution to a Bethe-Salpeter equation  (BSE) for electron-hole excitations .

In this Letter, we demonstrate excellent agreement between QMC, GW-BSE, and available experiment for excitation energies of silane (SiH4) and methane (CH4). Calculated ionization potentials, electron affinities, and first excited singlet and triplet energies also agree very well with highly reliable standard quantum chemistry approaches; however, the major advantage of our two current approaches is numerical efficiency. Both methods require moderate computation times, and their scaling with system size is to our knowledge superior to all other quantum-chemistry approaches for computing excited states. The excellent agreement between QMC and GW-BSE is encouraging given the dramatically different methodologies involved (i.e., the former relies on accurate total energy differences while the latter is based on perturbation theory).

In our QMC approach , variational Monte Carlo is employed to find an optimized correlated many-body trial function. This trial function, $`\psi _T(R)`$, is a product of Slater determinants, $`D_n`$, and a correlation factor ,

$$\psi _T=\underset{n}{}d_nD_n^{}D_n^{}\mathrm{exp}\left[\underset{I,i<j}{}u(r_{iI},r_{jI},r_{ij})\right],$$
(1)
where $`I`$ corresponds to the ions, $`i,j`$ to the electrons and $`r_{iI},r_{jI},r_{ij}`$ to the relevant separations. Parametrization and optimization of $`u(r_{iI},r_{jI},r_{ij})`$, which represents the electron-electron and electron-electron-ion correlations, is described in Ref. . In the Slater determinant part, we employ natural orbitals rather than Hartree-Fock (HF) or density functional orbitals . To eliminate most of the remaining variational bias we use the fixed node diffusion Monte Carlo method, which is based on the property that the operator $`\mathrm{exp}(\tau H)`$, where $`H`$ is the Hamiltonian, projects out the ground state of any trial function with the same symmetry and non-zero overlap, subject to the constraint that the nodes are unchanged . All QMC results presented here are from the diffusion Monte Carlo approach.

The GW-BSE Green’s function perturbation approach consists of three successive steps. In the first step a local density approximation (LDA)  calculation is performed to obtain the electronic ground state. In the second step the LDA results are used to construct the electronic self-energy operator $`\mathrm{\Sigma }`$ (also known as the mass operator). This is done within the $`GW`$ approximation  where $`G`$ is the (LDA) single-particle Green’s function and $`W`$ is the screened Coulomb interaction based on RPA dielectric screening. Solving the equation-of-motion for the single-particle Green’s function yields the quasiparticle states (i.e., the occupied (hole) states $`v`$ and the empty (electron) states $`c`$) of the molecule. In the third step the electron-hole interaction $`vc|K^{\mathrm{eh}}|v^{}c^{}`$ between the occupied states $`v`$ and the empty states $`c`$ is evaluated. The coupled electron-hole excitations $`|S=_{vc}A_{vc}^S|vc`$ and their excitation energies $`\mathrm{\Omega }^S`$ then result from the Bethe-Salpeter equation

$$(\epsilon _c^{\mathrm{QP}}\epsilon _v^{\mathrm{QP}})A_{vc}^S+\underset{v^{}c^{}}{}vc|K^{\mathrm{eh}}|v^{}c^{}A_{v^{}c^{}}^S=\mathrm{\Omega }^SA_{vc}^S.$$
(2)
For details, we refer to Refs. .

We have also calculated excited state energies as total energy differences using LDA and the generalized gradient approximation (GGA) with the PW91 functional . In addition, the small number of correlated electrons in the present systems permits us to employ several of the most accurate quantum chemical approaches: coupled cluster with single, double, and perturbationally triple (CCSD(T)) excitations , and complete active space self consistent field (CASSCF) with an active space of 8 electrons and 9 active orbitals. For these methods, extensive basis set tests showed that energy differences are converged to better than 0.02 eV with the cc-pVQZ  basis set. All calculations presented here other than QMC and GW-BSE have been carried out using the G98 package  and the cc-pVQZ basis set.

It is important to distinguish between the two different methodologies employed here. Using total energy methods, the ionization potential (IP) and electron affinity (EA) are computed as total energy differences. In the Green’s function, the QP energies are obtained in a different manner, namely as a solution of a quasiparticle Dyson’s equation .

The single-particle HF and LDA energy spectra for SiH4 are shown in Fig. 1. While both methods are in qualitative agreement with one another, there is a large quantitative discrepancy, as expected: HF usually overestimates HOMO-LUMO gaps and eigenvalue spacings, whereas LDA typically underestimates these quantities. Note that since the first excitations are from t2 to a1 orbitals, the ground and first excited singlet states of SiH4 are of different symmetry.

It is reassuring that all of the methods (see Table I) are in good agreement with experiment for the ionization potential of SiH4. The LDA, GGA, and HF results are $``$ 0.5 eV lower than experiment, indicating a slight bias towards the ionized state; however, note the dramatic improvement when compared to their single-particle highest occupied eigenvalue. The calculations indicate a small, negative (i.e., nonphysical) electron affinity ranging from -0.1 to -0.6 eV.

Within the GW-BSE approach, the behavior of the $`a_1`$ LUMO state deserves careful discussion (see Table II). In the underlying LDA calculation, the LUMO energy is below the vacuum level (see Fig. 1) and is thus given by a localized, bound wave function. Within first-order perturbation treatment of the self-energy operator $`\mathrm{\Sigma }`$, the orbitals are kept as LDA orbitals and the LUMO energy is shifted above the vacuum level by 1.7 eV (see second column of Tab. II). Being unbound now, the QP state changes its wave function from a localized to a more extended orbital. This is described by going beyond the diagonal evaluation of $`\mathrm{\Sigma }`$ in the LDA basis in the $`GW`$ calculation and taking off-diagonal elements of $`\mathrm{\Sigma }`$ into account, which are in fact non-zero. This effect, which was not accounted for in a previous work , lowers the LUMO state by as much as 0.8 eV (see third column of Tab. II). The LUMO energy of 0.3 eV is now in good agreement with the results of the other methods. Note that the LUMO state (being unbound) converges slowly with respect to the basis, due to its delocalized nature. The low-energy electron-hole excited states, on the other hand, that are the focus of this paper, are spatially localized since the electron is attracted to the hole; these excited states are much easier to converge.

Our results for SiH4 excitation energies for the lowest-energy triplet and singlet excited states are given in Table III. Calculation of spin-singlet excitations by total energy differences poses a serious challenge to traditional approaches. For example, the first excited singlet state is described by a 2 determinant wavefunction. Accurate correlation of this state was therefore only possible using the CASSCF method. The triplet state is significantly easier to converge within traditional approaches owing to the fact that it is the lowest energy state of the given spin multiplicity. The QMC approach does not encounter the same difficulties as the traditional approaches for the singlet state since the spin configuration is given by specification of the slater determinant – singlet and triplet states are equally as simple to simulate, provided each wavefunction is orthogonal to the ground state. In the case of the excited singlet state, the DMC energy was found to be rather insensitive (i.e., $`<`$0.1 eV) to the number of determinants in the trial function, indicating that the nodes of a single determinant closely resemble those of the 2-determinant wavefunction which has the correct spatial symmetry.

The SiH4 spin singlet results are in good agreement with experiment  for the three theoretical methods listed in Table III. The optically forbidden triplet excitation is difficult to obtain from experiments. Using electron collision spectroscopy, slight spectral structures were found at around 7-8 eV . We believe that the agreement among the theoretical methods presented here strongly suggests that the triplet excitation lies roughly 1 eV higher than previously reported.

The GW-BSE results of Table III differ from the ones of Ref.  (see also Table II). The reason is again due to the off-diagonal matrix elements of the self-energy operator when expanded in the LDA orbitals, which had been neglected in Ref. . As discussed above, the full self-energy operator allows the unoccupied states to become more delocalized. Therefore the overlap between the hole and electron states, as well as the electron-hole interaction, is significantly reduced. Concomitantly, the excitation energies (column 3 in Tab. II) are not reduced compared to column 2 (as the lower LUMO energy suggests), but are in fact increased. It is interesting to note that the triplet excitation energy is more affected than the singlet, which reduces the singlet-triplet splitting from 1.1 eV to 0.7 eV within GW-BSE.

Methane and silane are isoelectronic and similar in their valence shell structure; however, the stronger potential of carbon compared to silicon causes wavefunctions and densities to be less smooth which in turn can lead to difficulties and inaccuracies within ab initio calculations. We have carried out the same calculations described above for CH4. Our results, shown in Table IV, again demonstrate strong agreement between QMC and GW-BSE.

An important issue in the discussion of molecular spectra is the interaction between electronic excitations and the molecular geometry. Despite the high dimensionality of a five-atom molecule, we are able to discuss the most important aspects of geometric relaxation in the excited states and their consequences for the total energy. Table IV illustrates the impact of relaxation from the ground state Td symmetry for each excited state molecule. Note that the optimization has been restricted to C2v symmetry. Test calculations show that further lowering the symmetry results in a total energy reduction no larger than $``$ 0.1 eV.

For all three excited states shown here, geometry relaxation yields an enormous lowering of the total energy by as much as 1.9 eV. A minor contribution (about 0.2-0.4 eV) is related to the increase of the C-H distance in the excited states. The main contribution results from breaking the $`T_d`$ symmetry of the ground state. This effect can be understood from the molecular orbital scheme (see Fig. 1). Within the ground state, all three degenerate HOMO levels are fully occupied, i.e. the ground state is non-degenerate. In the excited states, one of the three HOMO states becomes half-filled, i.e. the excited states are three-fold degenerate. The symmetry-breaking from $`T_d`$ to $`C_{2v}`$ lifts the degeneracy of the three HOMO levels. In the single-particle spectrum of Fig. 1 two HOMO levels are shifted down and one is shifted up, forming the new HOMO level of the molecule in the reduced symmetry. In the optimum geometry of the excited states this level-splitting is about 4.2 eV, although the center-of-mass of the three levels remains roughly the same. The new HOMO level, which is half-filled in all three excited states, is non-degenerate, so the excited states are non-degenerate as well.

A closer analysis of the excited-state geometry shows that the increase of the C-H bond length is relatively small (no more than 0.07 Å). The main effect results indeed from the Jahn-Teller distortion of the molecule, which changes the bond angles from 109.47 degree (tetrahedral angle) to 94 degree for the small and 118 degree for the large bond angles. This distortion increases the total energy of the ground state by 1.5 eV. In the excited states, however, this increase of the total energy is overcompensated by the splitting of the three HOMO levels described above, i.e. by the upwards shift of the new HOMO level and a correspondingly large reduction of the excitation energies. In the experimental spectrum for CH4 , a weak shoulder is observed at about 10 eV, corresponding to a vertical excitation. Compared with the onset at 8.52 eV, this indicates a relaxation gain of about 1.5 eV, supporting our calculated results.

Similar effects on the excitation energy are also found for SiH4. The gain in total energy is 0.9 eV, i.e. slightly smaller than in CH4. In the measured absorption spectrum this corresponds to the difference between the energy of maximum absorption and the low-energy onset. This difference amounts to 0.6 eV in the absorption measurements of Ref. . Our minimum-energy transition for SiH4 is 8.3 eV (in GW-BSE), in excellent agreement with the measured onset of the spectrum at 8.2 eV.

JCG acknowledges useful discussions with G.M. Rignanese and E. Chang. Support was provided by National Science Foundation Grant No. DMR-9520554 and by the Director, Office of Energy Research, Office of Basic Energy Sciences, Materials Sciences Division of the U. S. Department of Energy under Contract No. DE-AC03-76SF00098. MR acknowledges financial support by the Deutsche Forschungsgemeinschaft (Bonn, Germany) under grant No. Ro 1318/2-1. Computational resources have been provided by NCSA and by NERSC.

$``$ Present Address: Lawrence Livermore National Laboratory, 7000 East Avenue, L-415, Livermore, CA 94550.


format: text

text:
# QCD vacuum structure

## 1 Introduction

The structure of the QCD vacuum has been the subject of many lattice investigations over the years. Two phenomena have attracted most attention: chiral symmetry breaking and confinement. Instantons have been conjectured to play a key role in driving the former . Here the lattice has entered the game by trying to provide non-perturbative information on the instanton ensemble. Activity in this field has already been reviewed in  and there is not much new to add this year. This is far from implying that we have reliably obtained all the information concerning QCD instanton dynamics. Fundamental issues like the size distribution or density of instantons are still not really settled, a point that will be briefly discussed in section 3.

Instantons have also come back on stage due to the beautiful instanton-monopole link arising at finite temperature . This re-establishes some equality between these, a priori, two very different objects; instantons are made up of monopoles but monopoles can also be viewed as periodic arrays of instantons .

Monopoles bring us to the most popular scenario for confinement, the dual superconductor picture , which is beautifully at work in SUSY gauge theories . ’t Hooft’s proposal of abelian projection  has been the subject of extensive lattice studies. The strongest version of this approach, based on abelian dominance, has been criticised in several respects . Different abelian projections do not seem to be equivalent. Should one of them be preferred, it would imply that a particular U(1) is selected. The confining flux tube should then be abelian in nature and fields neutral with respect to it would be unconfined. A naive particular prediction would be no area law for adjoint Wilson loops. Although such loops are indeed screened at large distances, numerical investigations indicate the existence of a regime where they exhibit a confining behaviour with approximate Casimir scaling . Moreover,  one also observes the more ‘fundamental’ center dominance, where instead of U(1) the relevant dynamical variables are assumed to be the ones in the center of the gauge group. Note that center dominance does not a priori solve the problem of Casimir scaling (adjoint fields are blind to the center). Although in  this was intended as a criticism, indications towards the possible relevance of center vortices to confinement arose from further work on the subject . This has boosted the revival of a proposal of confinement that for very long remained asleep . It is based on the fact that center vortices, and not U(1) monopoles, are the ‘confining’ configurations. Most of the activity during the past year has been devoted to this subject. I think it would be unfair to review QCD vacuum structure without acknowledging what has captured most of the attention. I will thus present my inexpert view on vortices in section 4.

Results concerning the dual superconductor approach to confinement will not be reviewed here (for recent reviews see ). Let me only mention that there is some agreement towards the fact that abelian dominance is indeed not really the issue . Possible resolutions of the puzzles mentioned before have been proposed. An important result is that one can study condensation of monopoles by constructing a monopole creation operator . The vev of such an operator is a disorder parameter for confinement, irrespective of the abelian projection used to define it. This is considered as a strong indication in favour of dual superconductivity. Related work concerning condensation of magnetic flux will be discussed in section 4.

The review is structured as follows. I will first concentrate on the less speculative phenomena, in particular the calculation of the topological susceptibility and the $`\eta ^{}`$ mass. This will be done in section 2. Section 3 presents a few other topics related to instantons, mostly concentrating on the instanton-monopole connection at finite temperature . In section 4 results concerning center vortices will be presented.

## 2 $`\chi `$ and the $`\eta ^{}`$ mass

A lot of work has been devoted over the past years to computing the quenched topological susceptibility on the lattice. In the limit of large number of colours, it is related to the $`\eta ^{}`$ mass through the Witten-Veneziano formula ,

$$\chi _q=\frac{Q^2}{V}=\frac{f_\pi ^2}{2N_f}(m_\eta ^2+m_\eta ^{}^22m_K^2).$$
(1)
Measuring the fluctuations of topological charge on the lattice has turned out to be a difficult enterprise (for a detailed description see ). I believe, by now, we can safely say that $`\chi _q`$ has been successfully determined both for SU(2) and SU(3). Continuum extrapolations of the available lattice results give (taking $`\sqrt{\sigma }=440\mathrm{M}\mathrm{e}\mathrm{V}`$)

$`\mathrm{SU}(2)`$ $`\chi _q=(214\pm 18\mathrm{M}\mathrm{e}\mathrm{V})^4`$
$`\mathrm{SU}(3)`$ $`\chi _q=(200\pm 18\mathrm{M}\mathrm{e}\mathrm{V})^4`$
in pretty good agreement with the large $`\mathrm{N}`$ prediction $`\chi _q(180\mathrm{M}\mathrm{e}\mathrm{V})^4`$.

The situation is different for the unquenched susceptibility. In the chirally broken phase

$$\chi =\frac{f_\pi ^2m_\pi ^2}{2N_f}+O(m_\pi ^4)m_q,$$
(2)
in contrast to the behaviour in the symmetric phase where we expect $`\chi m_q^{N_f}`$. The susceptibility is hence an observable clearly exhibiting the effect of dynamical fermions. The situation as of Latt’99 did not, however, look that promising. Available results from CP-PACS and the Pisa group  failed to see any chiral behaviour and exhibited an unquenched susceptibility independent of the quark mass. Results from UKQCD were more encouraging; the expected $`m_\pi `$ dependence was indeed observed, although the value of $`f_\pi `$ extracted from the slope of the susceptibility turned out about $`20\%`$ below the physical value.

New results from UKQCD , the Pisa group  and the SESAM-T$`\chi `$L collaboration  are available. They use respectively $`N_f=2`$ clover improved, staggered and Wilson fermions. In Fig. 1 a comparison between UKQCD and SESAM-T$`\chi `$L’s data for $`\chi r_0^4`$ vs $`m_\pi ^2r_0^2`$ is presented (with the scale set by Sommer’s $`r_00.49`$ fm). Good agreement is observed. From a fit to eq. (2) keeping terms in $`m_\pi ^4`$, UKQCD quotes a value $`f_\pi =105\pm 5_{10}^{+18}`$MeV, in very good agreement with the expected physical value $`f_\pi 93`$MeV.

New results by CP-PACS also indicating the expected chiral behaviour have been presented in S. Aoki’s plenary talk at this conference.

Still, the data from the Pisa group remain puzzling. We present them in Fig. 2. The errors are rather large, hence it is difficult to judge whether the result is really inconsistent with the expected chiral behaviour. A fit independent of $`m_q`$ gives $`\chi =(163\pm 6\mathrm{M}\mathrm{e}\mathrm{V})^4`$ with a $`chi^2/dof0.37`$, while a fit with a linear homogeneous dependence in $`m_q`$ gives $`chi^2/dof0.94`$. The largest $`\beta `$ point, which seems the most problematic, might be affected by finite size effects which are rather severe for staggered fermions . Taking it out improves the quality of the linear fit to $`chi^2/dof0.3`$. It is in any case clear that improved statistics and lighter quark masses are mandatory before deciding if these data represent a problem.

To finish this section let me mention some results by SESAM-T$`\chi `$L concerning the relation of the $`\eta ^{}`$ mass to topology . Direct evaluation of the $`\eta ^{}`$ mass turns out to be very tough because it involves OZI suppressed (‘disconnected’) diagrams. Although improved estimators for these diagrams would be required, the connection with topology is already clearly exposed in . The ratio between disconnected and connected diagrams is computed in sub-ensembles of Monte Carlo configurations characterised by their topological charge. A clear correlation between this ratio and the charge is observed. In particular, the disconnected piece vanishes in the sub-ensemble with charge $`|Q|1.5`$, indicating the degeneracy of the singlet and non-singlet states in the trivial topological charge sector.

Results from CP-PACS on the $`\eta ^{}`$ mass have also been presented in S. Aoki’s plenary talk.

## 3 Instanton constituents

Looking at instantons as composite objects is not a new idea. It is indeed at work in the non-linear O(3) model in two dimensions - see for instance \- where each instanton seems to be composed of a pair of opposite Coulomb charges. It is the melting of instantons into a plasma of charges that is argued to induce the mass gap. In the deconfined phase charges remain bounded in dipoles and long range forces are screened. Similar ideas were put forward for QCD in  with fractional topological charge configurations, the merons, as fundamental objects. Merons have a singular action density and turn out analytically intractable. However, non-singular self-dual objects with fractional topological charge exist on a torus with twisted boundary conditions. They have indeed been advocated as relevant for confinement  and in some sense they can be seen as the fundamental constituents of some integer charge instantons.

Recently this constituent nature of instantons has been explicitly exposed . The most general finite temperature instanton, the so-called caloron, has been constructed. It can be considered as a periodic array of instantons. For instanton size of the order of their separation, each instanton splits into N constituent BPS monopoles. The splitting takes place when the Polyakov loop at spatial infinity is non-trivial. Let me, for the sake of simplicity, concentrate on the SU(2) case. At spatial infinity the Polyakov loop is constant and parametrised by its trace: $`2\mathrm{cos}(2\pi w)`$. The masses of the constituent monopoles are determined by $`w`$ and $`\frac{1}{2}w`$. The standard Harrington-Shepard caloron with $`w=0`$  is associated to a massive monopole at the center of the caloron plus a massless delocalised one at infinity. It is worth mentioning that the monopoles are located where the Polyakov loop has two degenerate eigenvalues . This naturally makes contact with abelian projections. In a suitable gauge, out of the N monopoles, $`\mathrm{N}1`$ are BPS static and only one is non-static. It is this one which gives topological charge to the caloron. Even more, the fermionic zero mode of the $`Q=1`$ caloron is precisely localised on the non-static monopole .

This constituent nature of calorons has already turned out to be crucial in solving a long standing controversy in SUSY Yang-Mills . Calculations of the gluino condensate based on a direct semiclassical approximation did not agree with a ‘weak’ coupling expression with matter fields added (and then decoupled) to ensure calculability. Inclusion of the constituent monopole contributions in the semiclassical expansion has brought the two expressions to an agreement.

The way to make the constituent monopoles pop-up out of lattice calorons is by enforcing a non-trivial Polyakov loop at the spatial boundary. This is most naturally achieved by introducing twisted boundary conditions . Indeed, it can be seen that the SU(2) BPS monopole with mass $`4\pi ^2/\beta `$ precisely corresponds with the $`Q=1/2`$ finite temperature fractional instanton (allowed to have non-zero magnetic charge due to the twisted b.c.). Non-triviality of the Polyakov loop can also be achieved by freezing the time links at the spatial boundary. This has been used in  to investigate the relevant configurations in finite temperature SU(2) (extracted by cooling). The trace of the boundary Polyakov loop is frozen to be zero below $`T_c`$, while it is fixed, above $`T_c`$, to the observed loop average. In the confined phase, calorons made up of two charge-$`1/2`$ monopoles dominate. Above $`T_c`$ such monopoles are still present but coming in pairs of opposite topological charge. Dominance of $`Q=1/2`$ objects for SU(2) at $`T=0`$ has also been observed by imposing magnetic twisted boundary conditions . Of course the relevant dynamical question, in particular below $`T_c`$, is whether such dominance survives the thermodynamic limit irrespective of the boundary conditions used.

The highly non-trivial behaviour of calorons arises due to the strong overlap in the periodic array of instantons . Overlap effects are also important at $`T=0`$. As shown in  the action density of overlapping instantons differs considerably from the simple addition of single instanton profiles. Consequences for the extraction of the instanton size distribution from the lattice are important. In particular, when instantons are parallel oriented in colour space, large instantons are systematically missed by instanton finders. This is not an irrelevant issue since a large instanton component has been argued to give rise to confining behaviour for the Wilson loop . It is also relevant for the instanton liquid model  since it indicates a possible failure of the 1-instanton approximation in ensembles with densities analogous to the ones obtained in lattice simulations.

A few other works have dealt with this picture of instanton constituents although from a different point of view. Ref.  presents a low-energy effective action for QCD that incorporates $`\theta `$ dependence. It can be described in terms of a Coulomb gas of fractionally charged objects, resembling the constituent monopoles described above. Also merons have come back on stage in  where regularised lattice merons and their fermionic zero modes have been obtained. Finally, ref.  studies whether instantons melt into constituents in $`CP^{(N1)}`$ models with the result that melting does not take place for $`N3`$.

Let me now mention some other instanton related works. In  first evidence for a stronger correlation between instantons and anti-instantons in dynamical configurations has been measured. This has relevance for the instanton liquid model which predicts I-A correlations in the presence of dynamical quarks . In  the low-lying mesonic spectrum has been studied in ensembles of instantons with properties as obtained from the lattice, much in the spirit of the instanton liquid model. A first attempt to reproduce the real spectrum failed due to a mixture between physical states and free lattice modes. Removing these free modes by adding a perturbative background results in an excellent agreement. Finally, the controlled cooling technique developed in  to reduce uncertainties in the analysis of the instanton content of MC configurations has been extended to SU(3) .

## 4 Center vortices

The idea that center vortices might be relevant for confinement in Yang-Mills theories is also a very old one . It is perhaps in 3 dimensions where it becomes most appealing . The 3-D vortex is a topologically stable soliton of the theory. The vortex creation operator, $`\varphi `$, is a local field whose vev signals the spontaneous symmetry breaking of the $`Z_N`$ magnetic symmetry, mapping $`\varphi e^{\frac{2\pi in}{\mathrm{N}}}\varphi `$. This allows us to postulate an effective low-energy theory in terms of such a local field. In the broken phase there are N degenerate vacua, and the Wilson loop is the creation operator of the domain wall that separates them. This wall is stable because the vacuum surrounding it is. In this picture the string tension confining quarks is related to the tension of the wall.

An extension of these ideas to 4D is not straightforward. In 4D the vortex creation operator is no longer local (for a recent discussion see ). The 4-D soliton is string-like and the natural low-energy effective theory is a theory of strings. Still, one can, following ’t Hooft , define quantised electric and magnetic $`Z_N`$ fluxes. Consider Yang-Mills theory in a 4-D torus of periods $`l_\mu `$. The gauge potential is periodic in $`x_\mu `$ up to a gauge transformation $`\mathrm{\Omega }_\mu (x)`$. Univaluedness of $`A_\mu `$ implies

$$\mathrm{\Omega }_\mu (x+l_\nu )\mathrm{\Omega }_\nu (x)=e^{\frac{2\pi in_{\mu \nu }}{\mathrm{N}}}\mathrm{\Omega }_\nu (x+l_\mu )\mathrm{\Omega }_\mu (x)$$
(3)
with twist $`n_{\mu \nu }`$ (defined mod N) allowed to be non-zero due to blindness of $`A_\mu `$ to $`Z_\mathrm{N}`$. $`m_i\frac{1}{2}ϵ_{ijk}n_{jk}`$ counts the magnetic flux in the box along direction $`i`$, while $`k_in_{0i}`$ is dual to the electric flux $`e_i`$. Electric flux along a curve $`C`$ is generated by the Wilson loop, $`W(C)`$. The ’t Hooft loop operator, $`B(C)`$, non-local in the gauge fields, creates magnetic flux along $`C`$. ’t Hooft argued that in the absence of massless particles, the vacuum should be in one of two phases parametrised by the vevs of $`W`$ and $`B`$. In the Higgs phase they show respectively perimeter and area laws. The confined phase is dual to it and the roles of $`W`$ and $`B`$ are interchanged.

The free energy of a state of given $`\stackrel{}{e}`$ and $`\stackrel{}{m}`$ is

$$e^{\beta F(\stackrel{}{e},\stackrel{}{m})}=\underset{\stackrel{}{k}}{}e^{\frac{2\pi i\stackrel{}{k}\stackrel{}{e}}{\mathrm{N}}}Z(\stackrel{}{k},\stackrel{}{m})$$
(4)
with $`Z(\stackrel{}{k},\stackrel{}{m})`$ the partition function with twisted boundary conditions specified by $`\stackrel{}{k}`$ and $`\stackrel{}{m}`$. For simplicity I have assumed that the instanton $`\theta `$ angle is zero; the $`\theta `$ dependence gives rise to peculiar phenomena which will not be discussed here.

Euclidean symmetry gives rise to an exact electric-magnetic duality between free energies of different fluxes. In the presence of a mass gap such duality implies that, in the $`\beta `$, $`l_i\mathrm{}`$ limit, some of the fluxes have to be heavy. Moreover, either all the electric or all the magnetic fluxes are light, no mixing between them takes place. Suppose it is the electric fluxes which are heavy and confined in thin strings/vortices. Assuming that at large $`\beta `$ the free energy factorizes in an electric, $`F(\stackrel{}{e})`$, plus a magnetic, $`F(\stackrel{}{m})`$, part one derives

$$F((m,0,0))=2\lambda (1\mathrm{cos}(\frac{2\pi m}{\mathrm{N}}))l_1e^{\sigma l_2l_3}$$
(5)
with $`\sigma `$ the string tension of the electric confining string. The free energy of magnetic flux decreases exponentially as we let the box become large in the plane transverse to the flux. Magnetic fluxes spread over the whole volume and condense. This behaviour parametrises the confinement phase. It is derived solely from duality and the existence of heavy electric fluxes. Indeed, duality does not tell whether it is the electric or the magnetic fluxes which condense. In the Higgs phase the roles of electric and magnetic fluxes are interchanged.

Implementing twisted boundary conditions on the lattice is rather easy. For SU(2) a twist $`n_{\mu \nu }=1`$ can be enforced by flipping the sign of all the plaquettes sitting in, for instance, the upper right corner of each $`(\mu ,\nu )`$ plane. Notice that magnetic flux is only defined modulo 2, as is the number of twisted plaquettes per plane.

This year we have seen a revival of calculations of magnetic flux free energy on the lattice, both at zero  and at finite temperature . Already in  it was proposed to use magnetic twist as a probe for phase structure. In order to compute the free energy of magnetic flux, one computes the ratio of two partition functions: $`\mathrm{exp}\{\beta F(\stackrel{}{m})\}=Z(\stackrel{}{m})/Z(\stackrel{}{0})`$. Results at zero temperature  support the exponential behaviour indicated in (5). It would be nice to check if the coefficient of the exponential decay of the free energy does indeed agree with the electric string tension. At finite temperature there are results, both for 3  and 4 dimensions , corroborating the dual behaviour of ’t Hooft and Wilson loops. There is, for non-zero $`T`$, a difference between introducing the twist in space-time or space-space planes. As shown in  space-space ’t Hooft loops show area law, corresponding to deconfinement of space-time Wilson loops. However, space-time ’t Hooft loops are screened in agreement with the observation that the spatial string tension survives above $`T_c`$. Similar results are obtained in . Previous analytic calculations of the expectation value of the ’t Hooft loop at high temperature  also support this picture.

It is worth mentioning the comparison in  between disorder operators signalling monopole and magnetic flux condensation. Both of them behave in a very similar way giving a critical temperature compatible with the standard determinations.

Another issue much more difficult to settle, is whether indeed the disordering of Wilson loops is driven by the presence of thick magnetic vortices in the vacuum as advocated in .

Smooth vortex configurations do exist and have been obtained on the lattice from cooling . For this the use of twisted boundary conditions is again essential. Here one remark is important. With twisted b.c. such that $`\stackrel{}{k}\stackrel{}{m}0`$ (mod N) the topological charge is fractional and quantised in units of $`1/\mathrm{N}`$. We have already discussed fractional charge objects in connection to monopoles. Vortices found in  also carry fractional charge.

Another nice connection between vortices and topology has been put forward in . Based on a model that describes vortices as random surfaces , topology is incorporated by providing orientation to the surface. Non-trivial topological charge comes from non-globally-orientable surfaces describable as patches of equal orientation separated by monopole lines. A prediction for the zero temperature susceptibility of $`\chi _q(T=0)=(190\pm 15\mathrm{M}\mathrm{e}\mathrm{V})^4`$ is derived, in amazingly good agreement with lattice results - see sec. 2.

How to locate thick center vortices on lattice configurations has been the subject of a big debate this year. In  an approach very similar to ’t Hooft’s abelian projection was taken (other alternative approaches will not be discussed, for a recent review see ). Center vortices are located by fixing the so-called maximal center gauge (MCG), obtained by maximising the average of $`|\mathrm{Tr}(U_\mu (x))|^2`$. Let us concentrate on the SU(2) case. Center projection consists of replacing gauge fixed links by the closest $`Z_2`$ element. Center-projected (P-)vortices correspond to coclosed sets of plaquettes taking value (-1). It is claimed that the string tension from center-projected links ($`\sigma _{Z_2}`$) agrees with the full string tension ($`\sigma _{\mathrm{SU}(2)}`$), a phenomena dubbed as center dominance. The relevance of this is, however, obscured by the fact that center dominance appears to be obtained even without gauge fixing . The physicality of P-vortices has to be judged on a different basis. Tests in  correspond to the behaviour of Wilson loops pierced by even/odd number of P-vortices, and to the scaling of the P-vortex density. The behaviour of P-vortices across the deconfinement phase transition has also been studied .

The news this year is that maximal center gauge turns out to be severely affected by lattice Gribov copies. A first indication in this direction was provided in . If, prior to fixing MCG, the configurations are driven into a smooth gauge like the Landau gauge, center dominance is lost and the density of P-vortices dramatically reduced. Worrisome is that Landau preconditioning usually gives a higher local maximum than the one from direct MCG. Further evidence comes from  where several random copies of the same configuration are made, taking from them the one that gives the higher local maximum after MCG. The number of copies is extrapolated to infinity, with the result again that a very significant part of $`\sigma _{\mathrm{SU}(2)}`$ is lost in center projection, even in the continuum limit. The debate originated about this issue (see ) seems settled in  with a very careful study of the dependence on the number of gauge copies and the results stated above.

One can do better by performing a gauge fixing free of lattice Gribov ambiguities. This is the case of the Laplacian gauge, first introduced for abelian projection  and further extended to perform center gauge fixing  (see also ). The idea is to diagonalise the adjoint Laplacian and use its two lowest eigenvectors to fix the gauge. Here instantons, monopoles and vortices arise respectively as point, 1- or 2-dimensional singularities of the gauge fixing . Indeed, in the Laplacian gauge, center dominance is recovered, although only in the continuum limit. But we have by now repeatedly said that center dominance alone is not good enough. Further investigation on the vortex content of Laplacian gauge fixed configurations is still necessary.

## 5 Closing remarks

There is still a lot of work to do to unveil the mysteries of confinement. Perhaps the relative failure of the approaches taken so far is related to our insistence on describing confinement in ‘semiclassical’ terms. We tend to bear in mind that some underlying ‘classical’ fields (be it ‘fat’ monopoles, vortices or instantons) drive the phenomena. But attempts to identify them in non-perturbative ensembles have systematically led to problems. There might be some truth in it but the key ingredient seems to be still missing.

## Acknowledgements

I thank the organisers for a very enjoyable conference. I am indebted to José Luis F. Barbón, Tony González-Arroyo and Pierre van Baal for invaluable discussions over the years. I would also like to thank Philippe de Forcrand, Adriano Di Giacomo, Jeff Greensite, Tamas Kovacs, Alvaro Montero, Carlos Pena and Nucu Stamatescu.


text:
# Classical 1D maps, quantum graphs and ensembles of unitary matrices

## I Introduction

Although graphs served as models of physical systems for a long time, the interest in properties of spectra of the corresponding quantum systems has rather short history. Quantum graphs were introduced as a “toy” model for studies of quantum chaos by Kottos and Smilansky . They observed that spectral statistics of fully connected graphs is well reproduced by random matrix theory (RMT) and may be explained using the exact trace formula. Schanz and Smilansky  proved that the spectral correlation function of simple graphs coincides with CUE expression for $`2\times 2`$ matrices. The trace formula served to calculate the spectral correlation function and to explain the deviations from the RMT predictions for star graphs . Tanner introduced directed graphs and the corresponding unitary transfer matrices  and demonstrated fast convergence of their spectral properties to RMT statistics with increasing matrix size. A model for study of the Anderson localization by means of quantum graphs was introduced by Schanz and Smilansky , while other properties of quantum graphs were recently analyzed in . Continuous dynamics on graphs was recently discussed by Barra and Gaspard , who generated Markov process by means of Poincaré surface of section.

In this work we propose another approach leading to quantum graphs. We are interested in a certain class of one dimensional (1D) piecewise linear maps with infinite family of Markov partition into equal cells. For these systems we define a family of the corresponding graphs by means of bistochastic (doubly stochastic) transfer matrices $`B`$. If the transfer matrices are *unistochastic* (i.e. there exists a unitary matrix $`U`$ such that $`B_{ij}=|U_{ij}|^2`$) the graphs may be quantized. We may thus establish a family of the unitary transfer matrices by varying the lengths of bonds of graphs and investigate their spectra. We demonstrate a link between the character of dynamics of the classical system and the spectral properties of the corresponding ensemble of unitary matrices. For a class of regular 1D systems we show that the statistics of the spectrum is locally Poissonian. If the classical system has an arbitrary small, positive KS-entropy, the spectral statistics of the corresponding ensemble of unitary matrices displays fluctuations characteristic to Circular Unitary Ensemble (CUE) . These matrices are sparse (only few nonzero elements in each row and column) and by construction differ from typical random matrices generated according to the Haar measure on $`U(N)`$.

The paper is organized as follows. In section II we introduce the class of 1D maps of our interest. Section III contains an analysis of the structure and spectral properties of quantized graphs corresponding to regular systems. In section IV we study graphs corresponding to Bernoulli shift and to some systems with arbitrary small, but positive metric entropy. Section V is devoted to study of generic chaotic systems, for which we find CUE like spectral properties of the corresponding ensembles of unitary matrices. The concluding remarks are presented in section VI, while appendix A contains the condition necessary for a bistochastic matrix to be unistochastic and the calculation of free parameters for averaging over different graphs with the same topology. In appendix B we show the construction of unitary matrices corresponding to the classical map studied in section V.

## II One dimensional dynamical systems

Consider a discrete one dimensional mapping acting on an interval $`IIR`$

$$x_{n+1}=f(x_n),$$
(1)
where $`f:II`$ is piecewise linear. Assume that $`f`$ fulfils three conditions:

1. there exists a Markov partition of the interval $`I`$ on $`M`$ equal cells $`E_i,i=1\mathrm{}M`$ and $`f`$ is linear on each cell $`E_i`$,
2. for any $`yI`$

$$\underset{xf^1(\{y\})}{}\frac{1}{f^{}(x)}=1,$$
(2)
3. the finite transfer matrix $`B`$, describing the evolution of measures uniform on the cells of Markov partition, under the action of the system $`f`$ (the Frobenius-Perron operator) is unistochastic.

Condition (i) means that the graph of the function $`f`$ (later called *diagram* of $`f`$ to avoid the confusion with graphs corresponding to the map $`f`$) is composed of straight lines (neither horizontal nor vertical) with endpoints onto $`M\times M`$ grid on the square $`I\times I`$. Some examples of functions fulfilling this condition are plotted in figure 1. If we rescale $`I`$ to $`[0,1]`$ and write $`f(x)=c_jx+b_j`$ for $`xI_j`$ where $`_jI_j=I`$, then all $`c_j`$ have to be nonzero integers and also, the numbers $`b_j`$ and $`|I_j|/|I|`$ must be rational (multiplication by a sufficiently large $`M`$ makes all these numbers integer). Any subpartition of the interval $`I`$ into $`N=Mk`$ equal cells ($`kIN`$) is also a Markov partition, so we have an infinite family of Markov partitions for our models. The Frobenius-Perron operator, reduced to densities (measures) constant on each cell $`E_j`$, may be described by a matrix $`B`$ of finite size $`M`$. Let us write the elements of $`B`$ explicitly. The system fulfilling the condition (i) is fully characterized by a set of $`M`$ linear functions

$$f_i(x)=c_ix+b_i,x[\frac{i1}{M},\frac{i}{M}).$$
(3)
By construction $`c_i`$ and $`Mb_i`$ are both integers. It is convenient to express $`c_i`$ and $`b_i`$ in terms of two integers $`l_i`$ and $`m_i`$ as

$$c_i=m_il_i,b_i=\frac{m_ii(m_il_i)}{M}.$$
(4)
It is now easily seen that the only non-zero elements of $`B`$ in the $`i`$-th row are

$$B_{ij}=\frac{1}{|m_il_i|},\mathrm{min}(m_i,l_i)+1j\mathrm{max}(m_i,l_i).$$
(5)
Closer inspection of the structure of the above formula reveals the connection between the structure (positions of the non-zero elements) of the matrix $`B`$ and the diagram of the system $`f`$ rotated clockwise by $`\pi /2`$.

The element $`B_{ij}`$ is equal to the probability of going from the cell $`E_i`$ to $`E_j`$. From (5) $`B`$ is stochastic ($`_{j=1}^NB_{ij}=1`$), all probabilities of going out of the cell $`E_i`$ sum to unity. The probability of coming to the cell $`E_j`$ from $`E_i`$ is equal to mean $`1/|f^{}(x)|`$ over $`xE_i`$ and $`f(x)E_j`$. Thus due to condition (ii) the probabilities of transitions to any fixed cell $`E_j`$ also sum to one ($`_{i=1}^NB_{ij}=1`$), so the transition matrix $`B`$ is bistochastic. It is well known that all eigenvalues of a bistochastic matrix are located in the unit circle, and the leading eigenvalue is equal to one . The condition (iii) will be used in section V.

We may compute the KS-entropy of Markov chain generated by transition matrix $`B`$

$$H_{KS}(B)=\underset{i=1}{\overset{N}{}}\stackrel{~}{p}_i\underset{j=1}{\overset{N}{}}B_{ij}\mathrm{log}B_{ij},$$
(6)
where $`\stackrel{~}{p}`$ is the normalized left eigenvector of $`B`$ corresponding to the leading eigenvalue, $`\stackrel{~}{p}B=\stackrel{~}{p}`$, $`_{i=1}^N\stackrel{~}{p}_i=1`$. In condition (i) we require that $`f`$ is linear (not only piecewise linear) on each cell $`E_j`$. Then the Markov partition on $`M`$ equal cells is a generating partition of the system, so Eq. (6) gives the dynamical entropy of the system. The transition matrix is bistochastic in our case, so all components of $`\stackrel{~}{p}`$ are nonzero ($`\stackrel{~}{p}_i=1/N`$) and $`H_{KS}=0`$ if and only if all $`B_{ij}\{0,1\}`$. It is easy to see that if $`xI:|f^{}(x)|=1`$ the system is regular. We have showed that this is also a necessary condition, provided the conditions (i) and (ii) are fulfilled, since only the systems with $`|f^{}|=1`$ have all $`B_{ij}\{0,1\}`$. Note that in general there exist regular systems with slopes larger than one, but they do not satisfy the condition (ii). An example of such a system is plotted in Fig. 1(b).

## III Regular systems

Partition of the interval $`I`$ determines the transition matrix $`B_{ij}`$ describing the symbolic dynamics of the 1D map $`f(x)`$. We introduce a directed graph with vertices corresponding to each cell $`E_i`$ of the partition and bonds corresponding to all nonzero elements of the transition matrix $`B_{ij}`$. Changing the number $`N`$ of the cells in the partition we change the dimension of the matrix $`B^{(N)}`$, so we may obtain a family of graphs corresponding to any 1D map fulfilling the condition (i).

If we consider the regular systems with $`|f^{}|=1`$ the only nonzero elements of matrices $`B_{ij}`$ are equal to 1. The matrices of this type contain only one nonzero element equal to 1 in each row and column and act as a permutation on a vector with $`N`$ elements. The corresponding graph has only one incoming and one outgoing bond at each vertex, so its valency is equal to one. It is composed from the closed loops. We say that there is an inversion at vertex $`j`$ of the graph if $`f^{}=1`$ on the corresponding cell $`E_j`$ of the partition. Dividing the interval $`I`$ into $`M`$ equal cells we find that the longest loop of the graph has at most $`M`$ vertices. If we consider a subpartition of $`I`$ into $`N=Mk`$ cells, each vertex of the graph will be replaced by $`k`$ vertices of a new graph. The loops with an even number of inversions will form $`k`$ new disconnected loops of the same period. The loops with an odd number of inversions will create $`[k/2]`$ disconnected loops with period doubled ($`[]`$ denotes the integer part of a number) and also, if $`k`$ is odd, one loop of the initial period. Fig. 2 shows an example of a regular system, for which the minimal Markov partition contains $`M=3`$ cells: $`E_1,E_2,E_3`$. Due to $`f^{}=1`$ in $`E_1`$ there exists one inversion in the 3-vertices graph corresponding to this system. For larger $`N`$ the graph is composed of $`[k/2]`$ loops with 6-vertices and one 3-vertices loop for an odd $`k`$. The classical map has one periodic orbit with period 3 and a continuous infinite family of orbits with period 6. Note similarity in the spectrum of periods of periodic orbits of the map and the graph for large odd $`k`$.

The quantization of a directed graph may be related to problem of finding a unitary matrix $`U`$, such as $`|U_{ij}|^2=B_{ij}`$ for all indices $`i,j`$ . The matrix $`U`$ may be seen as a Floquet operator governing the discrete time evolution, while $`U^n`$ represents the evolution operator at discrete time $`t=n`$. Analyzing quantum maps, corresponding to time dependent periodical systems, one studies the statistical properties of phases of the unimodular eigenvalues of the unitary matrix $`U`$ . The elements of $`U`$ may be expressed by $`U_{ij}=r_{ij}e^{iL_{ij}}`$, where the modulus are fixed by the transition matrix $`r_{ij}=\sqrt{B_{ij}}`$, and the phases $`L_{ij}`$, corresponding to lengths of the bonds of the graph, are not specified by the corresponding 1D map. We only impose that they fulfil the condition of unitarity $`U^{}U=\mathrm{𝟏}`$. Varying the phases $`L_{ij}`$ we obtain a family of graphs – an ensemble of unitary matrices. This corresponds to varying the lengths of all bonds of the graph with its topology preserved. We are interested in spectral properties of the ensemble of unitary matrices defined in this way. If $`B`$ is a permutation matrix, there is only one nonzero element in each row and each column of $`U`$, so all phases $`L_{ij}`$ appearing in matrix $`U`$ are uncorrelated.

The $`N\times N`$ unitary matrix $`U`$ corresponding to a regular 1D system is composed of $`M`$ square diagonal or antidiagonal blocks of size $`k`$. For the system presented in Fig. 2 the matrix $`U`$ has the form

$$U=\left(\begin{array}{ccc}0& 0& V_1\\ V_2& 0& 0\\ 0& V_3& 0\end{array}\right),$$
(7)
where $`M=3`$, $`(V_1)_{ml}=\delta _{m,k+1l}e^{iL_m}`$ is an antidiagonal $`k\times k`$ matrix, $`(V_2)_{ml}=\delta _{ml}e^{iL_{m+k}}`$ and $`(V_3)_{ml}=\delta _{ml}e^{iL_{m+2k}}`$ are diagonal matrices. There are $`N`$ free phases $`L_m`$ (the column index is superfluous here so we omit it). The eigenvalues of block matrix $`U`$ may be expressed by means of the spectra of its blocks $`V_i`$

$$\text{eig}(U)={}_{}{}^{M}\sqrt{\text{eig}\left(V\right)}={}_{}{}^{M}\sqrt{\text{eig}\left(\underset{i=1}{\overset{M}{}}V_i\right)},$$
(8)
where the $`M`$-th root is a $`M`$ valued function. Unitary matrix $`V=_{i=1}^MV_i`$ of size $`k`$ is diagonal or antidiagonal, depending on parity of the number of inversions in the corresponding loop in the basic $`M`$-vertex graph. The nonzero elements of $`V`$ are equal to $`e^{i_j}`$, where $`_j`$, $`j=1\mathrm{}k`$ are sums of different phases $`L_m`$, so $`_j`$ are uncorrelated. If $`V`$ is diagonal the spectrum of $`U`$ consists of $`\{e^{i(_j+2\pi m)/M}\}_{m=1\mathrm{}M}^{j=1\mathrm{}k}`$. The eigenphases of $`U`$ belonging to the same interval of the length $`2\pi /M`$ stem from different series $`_j`$, so they are uncorrelated and the spectrum has locally Poissonian statistics. If $`V`$ is antidiagonal, its spectrum contains $`\sqrt{e^{i(_j+_{k+1j})}}`$ and the spectrum of $`U`$ is composed of roots of unimodular complex numbers of order $`2M`$. The character of the spectrum remains Poissonian on the intervals of the length $`\pi /M`$.

The above argument may be generalized for the situation, in which the permutation matrix $`B`$ is reducible, composed of several permutations with different periods $`M_1\mathrm{}M_p`$. In this case the eigenphases of $`U`$ are uncorrelated in intervals of length $`\pi /\mathrm{max}M_i`$. Fig. 3 shows the commonly used spectral statistics  calculated for the system shown in Fig. 2. Note that the spectral ridigity $`\mathrm{\Delta }_3(L)`$, which measures correlations in the spectrum at the length scale $`L`$, saturates at $`L25=N/6`$ (there are $`N=150`$ eigenphases rescaled to the mean density constant and equal to $`1`$). The divisor $`6`$ is equal twice the number of cells in minimal Markov partition, since an inversion in one of the cells occurs ($`f^{}(x)=1`$ for $`xE_1`$).

## IV Simple chaotic systems

The Bernoulli shift, $`f(x)=2x\text{mod}1`$, is one of the simplest 1D maps displaying chaotic dynamics. This system fulfils the conditions (i), (ii) and (iii). The graphs corresponding to this system were analyzed by Tanner . The system, graphs and matrices are plotted in Fig. 4. Existence of a certain interval for which $`|f^{}(x)|>0`$ causes the mean valency of the graphs to become larger than 1. Some vertices of the graphs have more than one outgoing bond, and this fact is related with a positive topological entropy of the corresponding classical map. The number of such vertices grows proportionally with $`N`$ and their presence is reflected in statistical properties of spectra of the corresponding unitary matrices. For the binary graphs all vertices have two outgoing bonds and as shown by Tanner  the spectral statistics converges fast to the predictions of CUE. The traces of the unitary evolution operator $`U^n`$ may be expressed by a sum over periodic orbits of the corresponding classical graph, analogous to the semiclassical Gutzwiller trace formula . For quantum graphs the formula

$$\text{Tr}U^n=\underset{\nu 𝒫𝒪(n)}{}A_\nu e^{iL_\nu },$$
(9)
is an identity, where $`𝒫𝒪(n)`$ is the set of periodic orbits of period $`n`$, the amplitude $`A_\nu =_{(ij)\nu }\sqrt{B_{ij}}`$ and the length $`L_\nu `$ is the sum of lengths of bonds belonging to the periodic orbit $`\nu `$.

We will show that every periodic orbit of the graph corresponds one to one to a periodic orbit of the dynamical system. Assume that $`\nu =j_1j_2\mathrm{}j_n`$ is a periodic orbit of period $`n`$ of the graph, where $`j_i`$ denote the vertices belonging to $`\nu `$ and $`|f^{}|>1`$ on at least one of the cells $`E_{j_c}`$. Let the functions

$$F_{j_i}:E_{j_{i+1}}xf^1(x)E_{j_i}$$
(10)
be the local inverse of $`f`$, where we identify $`j_1`$ and $`j_{n+1}`$. The function $`F_{j_c}`$ is a contracting mapping since $`|F_{j_c}^{}|\frac{1}{2}`$. We define

$$F=F_{j_1}F_{j_2}\mathrm{}F_{j_n}$$
(11)
which is a contraction of the cell $`E_{j_1}`$ onto itself, because $`|F_{j_i}^{}|1`$ for $`ic`$. The function $`F`$ corresponds to moving along the orbit backwards. Due to the Banach theorem there exists exactly one fixed point of $`F`$ in $`E_{j_1}`$, namely $`x_\nu =F(x_\nu )`$. Hence there is a periodic orbit of period $`n`$ in the dynamical system corresponding to the orbit $`\nu `$ on the graph. The inverse is also true, so there is one to one correspondence between the orbits of the 1D chaotic system and the corresponding graphs.

The sum (9) goes over periodic orbits of the classical chaotic map. The amplitudes $`A_\nu `$, built from probabilities $`B_{ij}`$, are equal to $`_{i=0}^{n1}1/|f^{}(f^i(x_\nu ))|^{1/2}`$, so they are equal to the inverse square roots of the instabilities of the orbits of the system. Only the lengths $`L_\nu `$ are not determined yet, so we can average over the lengths of the bonds imposing the constraints due to unitarity of quantum propagator $`U`$. With growing dimension $`N`$ the spectral statistics of quantized binary graphs tends rapidly to the CUE statistics. Tanner observed it analyzing spectral form factors for graphs with rationally independent bond lengths . Fig. 5 shows the spectral statistics averaged over an ensemble of $`10^5`$ matrices of $`N=100`$, corresponding to graphs with randomly chosen bond lengths (see appendix A).

We may also analyze a system, for which $`|f^{}(x)|>1`$ only in certain regions $`I_cI`$. The KS-entropy of such a system is equal to the entropy of a system equivalent to $`f(x)`$ in $`I_c`$ times the ratio $`\text{vol}(I_c)/\text{vol}(I)`$. This allows us to construct a system with an arbitrary small, positive KS-entropy. Consider the system presented in Fig. 6 with dynamical entropy $`H_{KS}=\frac{1}{2}\mathrm{ln}2`$ (equal to the topological entropy), as an illustration of such situation. The structure of the unitary evolution operator $`U`$ is the same as the structure of the bistochastic transition matrix $`B`$ which is determined by the shape of the diagram of $`f`$. Block structure visible in Fig. 6 allows us to compute the spectra as in Eq. (8), where $`V`$ is the product of all blocks of $`U`$. If there is an even number of antidiagonal blocks, the resulting matrix $`V`$ is a product of one matrix with Poisson distribution of the spectrum and a certain number of non-diagonal matrices with a CUE-like spectrum. In ref.  we analyzed the spectra of composed ensembles consisting of products of random matrices, each pertaining to a given ensemble of unitary matrices. The CUE-like spectrum was found to be robust with respect to the multiplication by matrices of other ensembles. This is related to the fact that the Haar measure (CUE) is invariant with respect to multiplication. Following these lines of arguments, we conclude that the spectrum of $`V`$ has a CUE-like properties. The spectrum of $`U`$ is equal to the $`\stackrel{~}{M}`$-th root of the spectrum of $`V`$ ($`\stackrel{~}{M}`$ denotes the number of blocks in $`U`$)

$$e^{i\varphi _U}=e^{i(\psi _j+2l\pi )/\stackrel{~}{M}},j=1\mathrm{}k,l=0\mathrm{}\stackrel{~}{M}1,$$
(12)
where $`e^{i\psi _j}`$ is the $`j`$-th eigenvalue of $`V`$, it consists of $`\stackrel{~}{M}`$ copies of the spectrum of $`V`$ contracted $`\stackrel{~}{M}`$ times and placed one after another, so its statistical properties are locally like these of CUE. Fig. 7 shows the statistical properties of an exemplary chaotic system presented in Fig. 6, for which all spectra consist of $`\stackrel{~}{M}=2`$ identical parts of length $`\pi `$. Note that the level spacing distribution and the spectral ridigity (for small $`L`$) show good agreement with the predictions of RMT. On the other hand the spectral form factor (i.e. the Fourier transform of the two point correlation function of the spectral density ) exhibits huge fluctuations since it measures correlations in the entire spectrum and reveals the presence of 2 copies of the same spectrum. However, the average taken over a small interval of $`\tau =n/N`$ displays a CUE-like behaviour.

Local agreement with the predictions of RMT is preserved even though there exist an inversion in the map $`f`$, which implies an antidiagonal blocks in $`B^{(N)}`$. To show this we analyzed a system consisting of 4 blocks, presented in Fig. 8, while the corresponding results are shown in Fig. 9.

In the argument presented in this section we have assumed that both parts of the system (one *regular* with $`|f^{}(x)|=1`$ and the other *chaotic* with $`|f^{}(x)|>1`$) are connected, so the probability of going from one to the other is positive. If this is not true, the system splits into two subsystems, one regular and one chaotic. The spectrum of such a system will contain a mixture of eigenvalues of both subsystems, so the level spacing statistics could be described by an appropriately adopted Berry-Robnik distribution .

## V Generic systems

A generic 1D piecewise linear system fulfilling the conditions (i), (ii) and (iii) is chaotic and $`|f^{}(x)|2`$ for all $`x`$ in $`I`$. The mean valency of the corresponding graphs is also equal to or greater than 2 and the topological entropy of these graphs $`H_{\text{top}}`$ is not smaller than $`\mathrm{ln}2`$. The topological entropy of the map $`f`$ is equal to the log of the spectral radius (absolute value of the largest eigenvalue) of the *connectivity matrix* $`\stackrel{~}{B}`$ , where $`\stackrel{~}{B}_{ij}=1`$ if $`B_{ij}>0`$, and $`\stackrel{~}{B}_{ij}=0`$ otherwise. It is worth to consider, which of these graphs one may *quantize* by constructing corresponding family of unitary matrices. It is equivalent to ask, which bistochastic transition matrices $`B`$, are *unistochastic*, so there exist unitary $`U`$, such that $`|U_{ij}|^2=B_{ij}`$. It is clear that not all dynamical systems fulfilling (i) and (ii) give rise to bistochastic Frobenius-Perron matrices, which are unistochastic. As an example let us consider the system with $`M=6`$ shown in Fig. 1(c)

$$\begin{array}{cc}f_1(x)=6x,\hfill & \hfill \text{for }0x<\frac{1}{6},\\ f_{2,3}(x)=3x\frac{1}{2},\hfill & \hfill \text{for }\frac{1}{6}x<\frac{1}{2},\\ f_{4,5,6}(x)=2x1,\hfill & \hfill \text{for }\frac{1}{2}x1.\end{array}$$
(13)
The corresponding transfer matrix

$$B=\frac{1}{6}\left[\begin{array}{cccccc}1& 1& 1& 1& 1& 1\\ 2& 2& 2& 0& 0& 0\\ 0& 0& 0& 2& 2& 2\\ 3& 3& 0& 0& 0& 0\\ 0& 0& 3& 3& 0& 0\\ 0& 0& 0& 0& 3& 3\end{array}\right]$$
(14)
is clearly bistochastic but not unistochastic. Indeed one checks that the necessary conditions for unistochasticity are not fulfilled (see Appendix A Eq. (A6), $`k=3`$ and $`l=5`$). This is the smallest matrix corresponding to a dynamical system of the considered type with this property, some other may be constructed for $`M=28=14+7+4+2+1`$ or other perfect numbers. In appendix A we present a necessary condition for unistochasticity and discuss some features of the ensemble of matrices $`U`$ corresponding to the same bistochastic matrix $`B`$.

One of the maps studied by Dellnitz et al.  and showed in Fig. 10, may serve as an example of a classical map, for which transition matrices $`B`$ are unistochastic (see Appendix B). Its KS-entropy is equal to $`\frac{3}{2}\mathrm{ln}2`$, while the topological entropy equals to $`\mathrm{ln}3`$. Instead of computing the spectral radius of $`\stackrel{~}{B}^{(4)}`$ it is enough to observe that each point in $`I`$ has exactly three preimages, and the topological entropy may be computed as log of the mean number of preimages . The average should be taken with respect to the Parry measure of maximal entropy , which is not easy to specify in the general case. However, since the number of preimages is constant, its knowledge is not necessary in this case to obtain the required result $`H_{\mathrm{top}}=\mathrm{ln}3`$. Numerical investigations performed for the corresponding ensemble of $`10^5`$ matrices of size $`N=200`$ confirmed a CUE-like character of spectral statistics in the entire circle of eigenvalues as shown in Fig. 11.

Instead of averaging over an entire ensemble of unitary matrices with random phases, we may study spectral statistics for one concrete member of the ensemble. Figure 12 shows the spectral statistics for a single matrix $`U`$, in which the phases, not determined by unitarity, being set randomly (see Appendix A). Observe an agreement with CUE predictions (the form factor $`K(\tau )`$ is averaged over the window of $`\tau `$, $`\mathrm{\Delta }\tau =0.08`$). Interestingly, for this system a CUE-like behaviour was also found for a single matrix $`U_0`$ (corresponding to the map plotted in Fig. 10) with all free parameters set to zero (thus exhibiting strong correlations between phases). Matrix $`U_0`$ is real (therefore orthogonal).

Although a single member of the ensemble displays CUE-like spectrum, the averaging over free phases provides larger statistics. The form factor may be expressed as a double sum over periodic orbits

$$K(\tau )=\frac{1}{N}|\text{Tr}U^n|^2=\frac{1}{N}\underset{\nu ,\nu ^{}𝒫𝒪(n)}{}A_\nu A_\nu ^{}e^{i(L_\nu L_\nu ^{})},$$
(15)
where $`\tau =n/N`$ and the averaging goes over all allowed realizations of the graph (all bond lengths $`L_{ij}`$, for which unitarity of $`U`$ is preserved). The difference of the lengths of orbits is equal to $`L_\nu L_\nu ^{}=_{(ij)\nu }L_{ij}_{(lm)\nu ^{}}L_{lm}`$. Inserting new variables $`\stackrel{~}{L}_{ij}=L_{ij}L_{i1}L_{1j}+L_{11}`$ we have

$$L_\nu L_\nu ^{}=\underset{(ij)\nu }{}(\stackrel{~}{L}_{ij}+L_{i1}+L_{1j}L_{11})\underset{(lm)\nu ^{}}{}(\stackrel{~}{L}_{lm}+L_{l1}+L_{1m}L_{11}),$$
(16)
where the lengths $`\stackrel{~}{L}_{ij},i,j=2\mathrm{}N`$ are fixed by the condition assuring unitarity of $`U`$ (see Appendix A), while $`(2N1)`$ remaining phases $`L_{i1}`$ and $`L_{1j}`$ are drawn randomly with a uniform distribution on $`[0,2\pi )`$. Due to averaging over $`L_{i1}`$, $`L_{1j}`$ the only non vanishing terms in the sum (15) arise from periodic orbits $`\nu `$ and $`\nu ^{}`$ of the same *degeneracy class* which visit each vertex of the graph the same number of times. This observation generalizes results of Tanner for binary graphs . Therefore the calculation of $`K(\tau )`$ can be reduced to a combinatorial problem after finding the fixed lengths $`\stackrel{~}{L}_{ij}`$. The calculations of form factor for fully connected graphs  and binary graphs  shows fast convergence of the spectral statistics to RMT predictions with increasing matrix size $`N`$.

## VI Conclusions

We proposed to study a certain class of piecewise linear 1D dynamical systems, for which the Markov transition matrix (of size $`M`$) is bistochastic. It has a structure of the diagram of $`f`$, defining the classical map, rotated clockwise by $`\pi /2`$. For each subpartition of the Markov partition we construct a graph and a family of unitary matrices, of size equal to a multiplicity of $`M`$. This procedure might be thus considered as a “quantization” of a $`1D`$ classical maps, in a sense that it produces a family of corresponding unitary matrices. The above quotation marks represent the fact that no canonical action-angle quantization of classical systems with 1D phase space is possible.

We demonstrated that the spectral properties of ensembles of unitary matrices defined in this way reflect the character of the classical dynamics. For regular maps the graphs are composed of several finite disconnected loops and the spectral statistics of the corresponding unitary evolution operators is (locally) Poissonian. On the other hand, chaotic systems (with arbitrarily small, but positive metric entropy), lead to ensembles of unitary matrices with (locally) CUE-like spectral fluctuations. These unitary matrices are sparse (only a few non zero elements in each row and column), in contrast to the random matrices, typical with respect to the Haar measure on $`U(N)`$.

We have shown that the classical 1D system and the corresponding family of graphs are topologically conjugated, since each periodic orbit of the map corresponds to exactly one periodic orbit of each graph. Applying finer and finer Markov partitions to the system, introduces a natural limit in the family of the corresponding graphs. It is the limit of number of vertices tending to infinity with the set of periodic orbits preserved. This seems to be the correct limit to reach random matrix theory predictions for the spectral statistics of quantized graphs. Spectral properties of a graph may be expressed by the exact trace formula based on a sum over periodic orbits. This sum goes over all periodic orbits of the 1D system and the amplitudes are equal to inverse square roots of the instabilities of each orbit. The phases in the sum remain free parameters of the model, with constraints due to unitarity of the evolution operator. We find conditions necessary for a bistochastic transfer matrix to be unistochastic, which allow us to perform averaging over random phases (bond lengths) with unitarity preserved.

## Acknowledgments

Financial support by Polish KBN grant no 2 P03B 009 18 (P.P.) and no 2 P03B 072 19 (K.Ż. and M.K.) is gratefully acknowledged. PP is grateful to H. Schanz for an invitation and a financial support which enabled him to attend the Conference on *Random Network Models* organized in Götingen in December 2000.

## Note added in proof

After this work was completed we learnt about very recent related works on unistochastic matrices of Berkolaiko  and Tanner .

## A Unistochastic matrices

For a real bistochastic matrix $`B`$ of size $`N`$

$$B_{jk}0,\underset{j=1}{\overset{N}{}}B_{jk}=1,\underset{k=1}{\overset{N}{}}B_{jk}=1,$$
(A1)
we want to find a unitary matrix $`U`$, such that $`|U_{jk}|^2=B_{jk}`$. We can manipulate with $`N^2`$ variable phases $`L_{jk}`$ of the matrix $`U`$

$$U_{jk}=\sqrt{B_{jk}}e^{iL_{jk}}.$$
(A2)
The unitarity of $`U`$ requires that

$$U^{}U=\mathrm{𝟏}\underset{j=1}{\overset{N}{}}\sqrt{B_{kj}B_{lj}}e^{i(L_{lj}L_{kj})}=\delta _{kl}.$$
(A3)
Equations (A3) for $`k=l`$ are fulfilled due to stochasticity of $`B`$, Eq. (A1). Equations (A3) for $`k<l`$ are complex conjugated of these for $`k>l`$. Transforming the system of equations (A3)

$$U(U^{}U\mathrm{𝟏})U^{}=UU^{}(UU^{}\mathrm{𝟏})=0UU^{}\mathrm{𝟏}=0$$
(A4)
we found that $`(N1)`$ equations of the system are depending on the other ones, because the diagonal equations of the righthand side are equivalent to conditions for bistochasticity of $`B`$ – the first sum in Eq. (A1). We subtracted the unity from $`N`$ since there is one common condition stemming from each of the sums in Eq. (A1), $`_{ij=1}^NB_{ij}=N`$. The number of independent real equations in (A3) equals thus $`(N1)^2`$.

We introduce new variables $`\stackrel{~}{L}_{kj}=L_{kj}L_{k1}L_{1j}+L_{11}`$, so the equations (A3) are equivalent to

$$\underset{j=1}{\overset{N}{}}\sqrt{B_{kj}B_{lj}}e^{i(\stackrel{~}{L}_{lj}\stackrel{~}{L}_{kj})}=\delta _{kl}.$$
(A5)
Since $`\stackrel{~}{L}_{k1}=\stackrel{~}{L}_{1j}=0`$, so the $`(N1)^2`$ equations (A3) depend on $`(N1)^2`$ real independent variables $`\stackrel{~}{L}_{kj},k,j=2\mathrm{}N`$. These equations give the values of $`\stackrel{~}{L}_{kj}`$. For $`N=2`$ there exists only one phase $`\stackrel{~}{L}_{22}`$ equal to $`\pi `$ as shown by Tanner . This value does not depend on the bistochastic matrix $`B_{jk}`$. Eqs. (A5) do not depend on $`L_{k1}`$ and $`L_{1j}`$, so we may average over them, assuming the uniform distribution over the unit circle, $`P(L)=\frac{1}{2\pi }`$ for $`L[0,2\pi )`$. This is equivalent to averaging of the matrix $`U=D_1\stackrel{~}{U}D_2`$ over diagonal unitary random matrices $`D_1`$ and $`D_2`$, with a fixed matrix $`\stackrel{~}{U}`$, fulfilling $`|\stackrel{~}{U}_{jk}|^2=B_{jk}`$.

The system of equations (A3) is solvable only if

$$kl\underset{m=1\mathrm{}N}{\mathrm{max}}\sqrt{B_{km}B_{lm}}\frac{1}{2}\underset{j=1}{\overset{N}{}}\sqrt{B_{kj}B_{lj}}\text{and}\underset{m=1\mathrm{}N}{\mathrm{max}}\sqrt{B_{mk}B_{ml}}\frac{1}{2}\underset{j=1}{\overset{N}{}}\sqrt{B_{jk}B_{jl}}$$
(A6)
which is the condition for the possibility of balance the maximal term in the sums (A3) by the rest and the analog condition for the sums stemming from (A4). This necessary condition for a bistochastic matrix to be unistochastic becomes sufficient in $`3\times 3`$ case. For $`N=2`$ all bistochastic matrices are orthostochastic (and therefore also unistochastic). However for $`N3`$ several examples of bistochastic matrices, which are not unistochastic, are known , they do not satisfy the condition (A6). Statistical properties of random unistochastic matrices are studied in .

## B Construction of a family of quantum maps

The construction of quantum graphs presented in section V hinge on the possibility of finding an unitary (or orthogonal) matrix corresponding to the bistochastic matrix obtained for a given dynamical system for an arbitrary Markov subdivision of the original partition. Rather then attempting the formulation of the most general theorem concerning the existence of such unistochastic matrices, we present here a construction for a particular case of the “four legs map” presented in Fig. 10. In other words we show that the all bistochastic matrices $`B^{(N)}`$ corresponding to this classical dynamical system are unistochastic.

The original bistochastic matrix for the “four legs map” introduced in  for the initial Markov partition into four cells reads

$$B^{(4)}=\frac{1}{4}\left(\begin{array}{cccc}2& 2& 0& 0\\ 1& 1& 1& 1\\ 1& 1& 1& 1\\ 0& 0& 2& 2\end{array}\right).$$
(B1)
Let us consider the following choice of a corresponding unitary matrix

$$U^{(4)}=\frac{1}{2}\left(\begin{array}{cccc}\sqrt{2}& \sqrt{2}& 0& 0\\ 1& 1& 1& 1\\ 1& 1& 1& 1\\ 0& 0& \sqrt{2}& \sqrt{2}\end{array}\right).$$
(B2)
The bistochastic matrix for the two times finer subdivision reads

$$B^{(8)}=\frac{1}{4}\left(\begin{array}{cccccccc}2& 2& 0& 0& 0& 0& 0& 0\\ 0& 0& 2& 2& 0& 0& 0& 0\\ 1& 1& 1& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 1\\ 1& 1& 1& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 1\\ 0& 0& 0& 0& 2& 2& 0& 0\\ 0& 0& 0& 0& 0& 0& 2& 2\end{array}\right).$$
(B3)
Using our knowledge of the unitarity of $`U^{(4)}`$ we want to find a unitary matrix $`U^{(8)}`$, such that $`B_{ij}^{(8)}=|U_{ij}^{(8)}|^2`$. This can be done by replicating twice the rows of $`U^{(4)}`$ with an appropriate shift. Thus $`U^{(8)}`$ may read as follows

$$U^{(8)}=\frac{1}{2}\left(\begin{array}{cccccccc}\sqrt{2}& \sqrt{2}& 0& 0& 0& 0& 0& 0\\ 0& 0& \sqrt{2}& \sqrt{2}& 0& 0& 0& 0\\ 1& 1& 1& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 1\\ 1& 1& 1& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 1& 1& 1\\ 0& 0& 0& 0& \sqrt{2}& \sqrt{2}& 0& 0\\ 0& 0& 0& 0& 0& 0& \sqrt{2}& \sqrt{2}\end{array}\right).$$
(B4)
Note that the number of the non zero elements in $`U^{(8)}`$ is two times larger than the number of the non zero elements in $`U^{(4)}`$. By replicating each row of the matrix $`U^{(4)}`$ $`k`$ times we obtain in the similar way unitary matrices of the size $`4k`$.

The success of this construction is due to the fact that the nonzero elements of the first and fourth rows of $`U^{(4)}`$ are equal, which implies orthogonality of the following pairs of rows in $`U^{(8)}`$: (2,3), (2,5), (4,7) and (6,7). However this is not the only way of constructing the family of unitary matrices $`U^{(N)}`$. For example we found another family $`\stackrel{~}{U}^{(N)}`$ originating from

$$\stackrel{~}{U}^{(4)}=\frac{1}{2}\left(\begin{array}{cccc}\sqrt{2}& \sqrt{2}& 0& 0\\ 1& 1& 1& 1\\ 1& 1& 1& 1\\ 0& 0& \sqrt{2}& \sqrt{2}\end{array}\right),$$
(B5)
even though the procedure described above fails in this case. This fact allows us to hope that the class of classical 1D maps leading to quantum graphs (fulfilling the conditions (i), (ii) and (iii)) is not trivial.