text:
# Simulation of 𝑛_𝑓=3 QCD by Hybrid Monte Carlo

## 1 Introduction

Recent lattice QCD simulations include effects of dynamical fermions. Due to the algorithmic limitation of the standard Hybrid Monte Carlo (HMC) algorithm , those simulations are limited to even numbers of degenerate flavors. In order to include dynamical effects correctly, simulations of QCD with three flavors (u,d,s quarks) are desirable. Simulations with an odd number of flavors can be performed using the R-algorithm . This algorithm, however, is not exact: it causes systematic errors of order $`\mathrm{\Delta }\tau ^2`$, where $`\mathrm{\Delta }\tau `$ is the step-size of the Molecular Dynamics evolution. A careful extrapolation to zero step-size is therefore needed to obtain exact results. Nevertheless, it is common practice to forego this extrapolation and to perform simulations with a single step-size chosen small enough that the expected systematic errors are smaller than the statistical ones. We want to point out that there is an alternative to the R-algorithm, which gives arbitrarily accurate results without any extrapolation.

Lüscher proposed a local algorithm, the so-called ”Multiboson algorithm” , in which the inverse of the fermion matrix is approximated by a suitable Chebyshev polynomial. Originally he proposed it for two flavors QCD. Boriçi and de Forcrand  noticed that the determinant of a fermion matrix can be written in a manifestly positive way using a polynomial approximation, so that one can simulate odd flavors QCD with the multiboson method. Indeed, using this method, one flavor QCD was simulated successfully . The same polynomial approximation can be applied for the HMC . Actually, in the development stage of Ref., one flavor QCD was also simulated by HMC and it was confirmed that the two algorithmically different methods — multiboson and HMC — give the same plaquette value . Here we give the formulation of the HMC algorithm with odd flavors and perform $`n_f=3`$ QCD simulations. Then we compare our results with those of the R-algorithm.

## 2 Formulation

### 2.1 $`n_f=2`$

The application of Lüscher’s idea  to $`n_f=2`$ QCD HMC was first made by the authors of Ref., and later by . The lattice QCD partition function with $`n_f=2`$ degenerate quark flavors is given by

$$Z=𝑑UdetD^2\mathrm{exp}(S_{gauge}),$$
(1)
where $`D`$ is the fermion matrix and in this study we use Wilson fermions. In the formulation of the HMC algorithm, $`detD^2`$ is treated as

$$detD^2𝑑\varphi ^{}𝑑\varphi \mathrm{exp}(\varphi ^{}D^1D^1\varphi ),$$
(2)
where the $`\gamma _5`$ hermiticity of the fermion matrix $`D`$, i.e. $`D=\gamma _5D^{}\gamma _5`$, is used.

Introducing momenta $`P`$ conjugate to the link variables $`U`$, the partition function is rewritten as

$$Z=𝑑U𝑑P\mathrm{exp}(H),$$
(3)
where the Hamiltonian $`H`$ is defined by

$$H=\frac{1}{2}P^2+S_{gauge}+\varphi ^{}D^1D^1\varphi .$$
(4)
This Hamiltonian is used for the Molecular Dynamics (MD) simulation of the standard HMC algorithm. Eq.(4) has a computational difficulty in MD simulations since one must solve $`x=D^1\varphi `$ type equations which in general take a large amount of computational time for a large fermion matrix and/or for a small quark mass.

Following Lüscher , the inverse of $`D`$ can be approximated by a polynomial:

$$1/DP_n(D)\underset{k=1}{\overset{n}{}}(DZ_k),$$
(5)
where $`Z_k`$ are the roots of the polynomial $`P_n(D)`$. We choose $`Z_k=1\mathrm{exp}(i2\pi k/(n+1))`$.

Replacing $`D^1`$ in eq.(4) by $`P_n(D)`$ we obtain an approximate Hamiltonian,

$$H_n=\frac{1}{2}P^2+S_{gauge}+\varphi ^{}P_n(D)^{}P_n(D)\varphi .$$
(6)
An advantage of using $`H_n`$ is that no solver calculation is required in the MD evolution. Instead, one needs $`n`$ multiplications by the matrix $`D`$.

$`H_n`$ does introduce some systematic errors from the polynomial approximation. For the $`n_f`$=even case, however, these errors are easily corrected by using the exact Hamiltonian of eq.(4) at the Metropolis step. This guarantees that configurations will be distributed according to the exact measure $`detD^2`$, for any polynomial $`P_n(D)`$.

However, the domain of convergence of $`P_n(D)`$ is bounded by a circle centered at $`(1,0)`$ which goes through the origin. If all eigenvalues of $`D`$ fall inside this domain, $`P_n(D)`$ converges exponentially. Otherwise, $`P_n(D)`$ does not converge, which may happen for some exceptional configurations. Our algorithm will tend to reject these configurations at the Metropolis step, leading to extremely long autocorrelation times. This domain of convergence can be changed by adopting another approximating polynomial. However, the origin must be excluded. Together with connectedness and conjugate symmetry of the spectrum, this implies that the real negative axis is always excluded from the domain of convergence for any polynomial. Configurations with real negative Dirac eigenvalues will be rejected by our polynomial algorithm.

### 2.2 $`n_f=1`$

After invention of the multiboson algorithm, Boriçi and de Forcrand  noticed that a single $`detD`$ can be treated in a manifestly positive way and an $`n_f=1`$ multiboson simulation was performed to study thermodynamics of $`n_f=1`$ QCD .

As before, the inverse of the fermion matrix $`D`$, using a polynomial of degree $`2n`$, is approximated as

$$1/D\underset{k=1}{\overset{2n}{}}(DZ_k),$$
(7)
where $`Z_k=1\mathrm{exp}(i2\pi k/(2n+1))`$. Noticing that the $`Z_k`$’s come in complex conjugate pairs, eq.(7) is rewritten as

$$1/D\underset{k=1}{\overset{n}{}}(D\overline{Z}_k)(DZ_k).$$
(8)
Using the $`\gamma _5`$ hermiticity of the fermion matrix, we find that $`det(D\overline{Z}_k)=det(DZ_k)^{}`$. Thus the determinant of $`D`$ is written as

$$det(D)det(T_n^{}(D)T_n(D))^1,$$
(9)
where $`T_n(D)_{k=1}^n(DZ_k)`$, and then we obtain

$$det(D)𝑑\varphi ^{}𝑑\varphi \mathrm{exp}(\varphi ^{}T_n^{}(D)T_n(D)\varphi ).$$
(10)
The term $`\varphi ^{}T_n^{}(D)T_n(D)\varphi `$ is manifestly positive. Then we may define the Hamiltonian of $`n_f=1`$ QCD as

$$H=\frac{1}{2}P^2+S_{gauge}+\varphi ^{}T_n^{}(D)T_n(D)\varphi .$$
(11)
With this Hamiltonian there is no difficulty to perform HMC algorithm. To improve efficiency and accuracy, one may use a polynomial of lower degree $`n`$ during the Molecular Dynamics trajectory, and a much higher degree $`mn`$ for the Metropolis step . The domain of convergence of the approximation eq.(10) is the same as for $`n_f=2`$. Exceptional configurations for which eigenvalues fall outside this domain will likewise be rejected at the Metropolis step. A further difficulty is that the sampled measure is $`det(T_m^{}(D)T_m(D))^1`$, which for exceptional configurations differs from the desired $`detD`$, increasingly so with $`m`$.

### 2.3 $`n_f=2+1`$

The partition function of $`n_f=2+1`$ QCD is given by

$$Z=𝑑Udet\stackrel{~}{D}^2detD\mathrm{exp}(S_{gauge}),$$
(12)
where the notations $`\stackrel{~}{D}`$ and $`D`$ are introduced to distinguish the two different quark masses. Using eq.(2) for $`det\stackrel{~}{D}^2`$ and eq.(10) for $`detD`$,

$$\begin{array}{c}det\stackrel{~}{D}^2detD\hfill \\ 𝑑\stackrel{~}{\varphi }^{}𝑑\stackrel{~}{\varphi }𝑑\varphi ^{}𝑑\varphi \mathrm{exp}(\stackrel{~}{\varphi }^{}\stackrel{~}{D}^1\stackrel{~}{D}^1\stackrel{~}{\varphi }H_T),\hfill \end{array}$$
(13)
where $`H_T\varphi ^{}T_n^{}(D)T_n(D)\varphi `$. We define $`n_f=2+1`$ Hamiltonian by

$$H=\frac{1}{2}P^2+S_g+\stackrel{~}{\varphi }^{}\stackrel{~}{D}^1\stackrel{~}{D}^1\stackrel{~}{\varphi }+H_T.$$
(14)
Two remarks are in order: $`(i)`$ as for $`n_f=1`$, one could use during the MD trajectory a polynomial of lower degree than for the Metropolis step; $`(ii)`$ the two bosonic fields $`\varphi `$ and $`\stackrel{~}{\varphi }`$ could be replaced by a single one, with action $`\varphi ^{}T_n^{}(D)\stackrel{~}{D}^1\stackrel{~}{D}^1T_n(D)\varphi `$. For simplicity, in this exploratory study we use two distinct bosonic fields and a single approximating polynomial.

## 3 Convergence

### 3.1 $`n_f=2`$

To see the rate of convergence of $`P_n(D)`$, we calculate the quantity $`X_n=\varphi ^{}P_n^{}(D)P_n(D)\varphi `$. In the limit $`n\mathrm{}`$, $`X_n`$ goes to $`X_{exact}\varphi ^{}D_{}^{}{}_{}{}^{1}D^1\varphi `$. First, we choose $`X_{exact}=\eta ^{}\eta `$ where $`\eta `$ is a random gaussian vector. Then the vector $`\varphi `$ is set to $`\varphi D\eta `$. The accuracy of $`X_n`$ is measured by the difference between $`X_n`$ and $`X_{exact}`$. We use a random gauge configuration. Figure 1:(top) shows $`X_n`$ versus the degree $`n`$ for different quark masses. Here the same $`\eta `$ is used for each calculation of $`X_n`$. $`X_n`$ converges to one value as $`n`$ increases, but at high degree $`n`$, $`X_n`$ diverges, which can be understood due to the rounding errors of our computer, where calculations are performed with 64-bit accuracy. Figure 1:(bottom) shows the accuracy of $`X_n`$ by $`|X_nX_{exact}|`$. Exponential convergence is seen for each quark mass, but the rate of convergence is slow for small quark masses.

### 3.2 $`n_f=1`$

We do the same analysis as for $`n_f=2`$, but for $`n_f=1`$, the value of $`X_{exact}`$ is not known. So we calculate the quantity $`X_n=\varphi ^{}T_n^{}(D)T_n(D)\varphi `$, where the vector $`\varphi `$ is a gaussian random vector, and we use a random gauge configuration. We assume that $`X_n`$ goes to a certain value in the limit of $`n\mathrm{}`$. Figure 2:(top) shows $`X_n`$ as a function of degree $`n`$. $`X_n`$ seems to converges to a certain value when the degree $`n`$ increases. At high degree $`n`$, $`X_n`$ diverges as in the case of $`n_f=2`$.

To see the rate of convergence, we calculate $`|X_nX_{max}|`$ where $`X_{max}`$ is defined by $`X_{max}=X_m`$, $`mn`$. Due to the rounding errors, we can not take $`m`$ very large. We take a maximum number $`m`$ where the rounding errors still do not appear. Figure 2:(bottom) shows $`|X_nX_{max}|`$ as a function of degree $`n`$. The dips seen in the figure are just due to the fact that at those points $`X_n=X_{max}X_m`$. The convergence seems to be exponential, but the rate of convergence is slow for small quark masses as in the $`n_f=2`$ case.

## 4 Simulations

Simulations of three flavors QCD are performed on an $`8^2\times 10\times 4`$ lattice at $`\beta =5.0`$ with $`\kappa =0.130`$ and $`0.160`$. We measure the plaquette and Polyakov loop varying the degree $`n`$ and compare them with those from the R-algorithm obtained with a step-size $`\mathrm{\Delta }\tau =0.01`$ . Figures 3 and 4:(top) show the plaquette as a function of $`n`$ at $`\kappa =0.130`$ and 0.160, respectively. Except for very small $`n`$, the results from the HMC algorithm agree with those from the R-algorithm within statistical errors. Results of the Polyakov loop are shown in Figures 3 and 4:(bottom). Except for a small discrepancy seen in Figure 3, the results from the HMC algorithm are in agreement with those from the R-algorithm. Note that convergence is not monotonic in $`n`$.

## 5 Conclusions

We formulated an odd-flavor HMC algorithm using a polynomial approximation. Simulations of three flavors QCD were performed. We found that the plaquette values are consistent with those from the R-algorithm at very small step-size. In principle the HMC algorithm is able to simulate any flavors of QCD, with arbitrary accuracy and without extrapolation \[as long as all Dirac eigenvalues are not real negative\]. However the rounding errors should be under control when we use a large lattice or/and small quark masses where one may need a polynomial of high degree $`n`$ to achieve sufficient approximation.

This work was partially supported by the Minsitry of Education, Science, Sports and Culture, Grant-in-Aid, No.11740159.


format: text

text:
# Numerical analysis of the magnetic-field-tuned superconductor-insulator transition in two dimensions

## Abstract

Ground state of the two-dimensional hard-core-boson model subjected to external magnetic field and quenched random chemical potential is studied numerically. In experiments, magnetic-field-tuned superconductor-insulator transition has already come under through investigation, whereas in computer simulation, only randomness-driven localization (with zero magnetic field) has been studied so far: The external magnetic field brings about a difficulty that the hopping amplitude becomes complex number (through the gauge twist), for which the quantum Monte-Carlo simulation fails. Here, we employ the exact diagonalization method, with which we demonstrate that the model does exhibit field-tuned localization transition at a certain critical magnetic field. At the critical point, we found that the DC conductivity is not universal, but is substantially larger than that of the randomness-driven localization transition at zero magnetic field. Our result supports recent experiment by Marković et al. reporting an increase of the critical conductivity with magnetic field strengthened.

PACS: 74.76.w Superconducting films, 71.23.An Theories and models; localized states, 75.40.Mg Numerical simulation studies, 68.35.Rh Phase transitions and critical phenomena, 75.10.Nr Spin-glass and other random models.

Keywords: magnetic-field-tuned superconductor-insulator transition, exact-diagonalization method, randomness, critical conductivity, dynamical critical exponent, finite-size-scaling method.

## 1 Introduction

Scaling argument of Abrahams, Anderson, Licciardello and Ramakrishnan  states that in two dimensions, infinitesimal amount of quenched randomness should drive itinerant extended states to localize. That is, at absolute zero temperature, the conductivity should vanish, if there exist any randomnesses. There are, however, some exceptions where the above description fails. Metal-insulator transition found in MOSFET devise , for instance, is one of recent hot topics arousing much attention. In the above-mentioned scaling theory, the randomness perturbation appears to be marginal so that some unexpected factors, for example, many-body interaction and external magnetic field, would possibly change the scenario.

Suppose that there exists an attractive interaction among particles. At the ground state,the particles would be unstable against bose condensation so that the system would fall into a superconducting state. Localization transition from the superconducting state is apparently out of the scope of the conventional localization theory, and has been studied extensively so far: In experiments, the transition is observed for metallic films  and Josephson-junction arrays . In essence, $`{}_{}{}^{4}\mathrm{He}`$ film adsorbed by porous media belongs to the same physics as well , although one cannot measure its electronic conductivity. One of the main concerns, particularly in experiments, would be the conductivity at the localization point. An argument claims that localization transition occurs at a universal condition . That is, at the transition point, irrespective of the samples examined, the electrical conductivity stays $`(2e)^2/h`$ with the electron charge $`e`$ and the Plank constant $`h`$. A great number of experiments are trying to validate this issue. From theoretical view point, the criticality itself is a matter of interest . It is to be stressed that the localization transition occurs at the ground state (zero temperature). Such ground-state criticality differs significantly from that of the finite-temperature transitions, and belongs to a peculiar universality class: From the path-integral viewpoint for the partition function, $`d`$-dimensional quantum system is regarded as a $`(d+1)`$-dimensional classical system. The system size along the imaginary-time direction is given by the inverse temperature, which diverges at zero temperature. Critical fluctuation along the imaginary-time direction (quantum fluctuation) contributes to the nature of critical phenomenon, giving rise to a new universality class. This extra contribution is characterized by the dynamical critical exponent $`z`$; a formal scaling argument on the superconductor-insulator transition will be found in the literatures . We stress that such arguments on criticality are not of mere theoretical interest, but are readily measurable in experiments; see Table 1.

In order to tune the localization transition, either magnetic field strength or amount of randomness (film thickness) has to be adjusted carefully. Because the former is far more advantageous for fine tuning, numerous experiments are devoted to the magnetic-field-tuned localization transition. In Table 1, we presented a summary of experimental results. Among them, we would like to draw reader’s attention to the recent experiment by Marković et al. , who scanned both two tuning parameters fairly systematically.

While in experiments, such the subtle issues have already come under through investigation, in computer simulation, on the contrary, it has not yet been done to incorporate the external magnetic field: So far, transitions driven solely by the random chemical potential (with zero magnetic field) have been simulated . Implementation of the magnetic field brings about a complication such that the hopping amplitude becomes complex number. Because of this, quantum Monte-Carlo simulation does not work; the inefficiency would be even more serious than that of the so-called negative sign problem. (Note that here, the magnetic field is not meant to couple to the moment of spin (as in usual spin models), but it rather couples to the charge (hopping) through the gauge twist.) Simulation of ref.  is actually dealing with both magnetic field and randomness. In the model, however, the randomness is not of random chemical potential, but is of random hopping strength; the latter may lie rather out of present scope: Note that the chemical potential term (particle number) and the superconductivity order parameter (gauge) are conjugate variables. Co-existence of these mutually conjugate terms is the source of fascination of this physics. On the contrary, it causes the aforementioned technical complications.

Here, in order to overcome the difficulty of, so to speak, ‘complex-sign problem,’ we employed the exact-diagonalization method. This method has not been used so frequently in the course of the studies on this issue. Perhaps, the limitation of the tractable system sizes had been worried over. Yet, recent extensive simulations and subsequent careful data analyses are guaranteeing its validity and usefulness . In particular, the exact diagonalization has an advantage over others that the method gives dynamical response functions  such as the AC conductivity (without resorting to the maximum-entropy method, for instance). Of course, it does manage the complex Hamiltonian elements.

The rest of this paper is organized as follows. In the next section, we explain our model Hamiltonian of the two-dimensional hard-core-boson model with magnetic field and random chemical potential. Details how we had chosen the gauge will be explicated there. In Section 3, we present numerical results. For the first time, with use of the finite-size-scaling method, we observe clear evidences of the onset of field-tuned localization transition. After determining the transition point, we evaluate experimentally accessible quantities such as the critical conductivity and the dynamical critical exponent. The last section is devoted to summary and discussions.

## 2 Model Hamiltonian — two-dimensional hard-core-boson model

In this section, we explain our model Hamiltonian and some of its physical ingredients. The model we simulated is the hard-core-boson model on the square lattice subjected to external magnetic field and quenched random chemical potential. The Hamiltonian is given by,

$$=\frac{J}{2}\underset{ij}{}\left(\mathrm{e}^{\mathrm{i}\varphi _{ij}}b_i^{}b_j+\mathrm{e}^{\mathrm{i}\varphi _{ij}}b_j^{}b_i\right)+\underset{i}{}H_i(2b_i^{}b_i1).$$
(1)
$`_{ij}`$ denotes the summation over all nearest neighbors on the square lattice of size $`L\times L`$. Periodic-boundary condition is imposed as for the horizontal ($`x`$-axis) direction. The operators $`\{b_i\}`$ ($`\{b_i^{}\}`$) are the hard-core-boson annihilation (creation) operators at site $`i`$. The gauge twist angles $`\{\varphi _{ij}\}`$ are chosen such as $`\varphi _{ij}=(i_y(L+1)/2)\pi B`$ for the horizontal $`i`$-$`j`$ pair of $`j`$ right next to $`i`$, and $`\varphi _{ij}=0`$ for the vertical $`i`$-$`j`$ pair. Because of this choice of gauge, a particle hopping clockwise around a plaquette feels, as a whole, the gauge angle of $`\pi B`$. That is, the bosons are subjected to the uniform magnetic flux $`\pi B`$ per plaquette. Hereafter, we call $`B`$ ‘magnetic field.’ (Therefore, the situation of half flux quantum per plaquette is realized at the magnetic-field strength $`B=1`$.) Quenched random chemical potentials $`\{H_i\}`$ distribute over the sector $`[\sqrt{3}\mathrm{\Delta },\sqrt{3}\mathrm{\Delta }]`$ uniformly (rectangular distribution). Therefore, the strength of randomness is controlled by the parameter $`\mathrm{\Delta }=[H_i^2]_{\mathrm{av}}`$, where $`[\mathrm{}]_{\mathrm{av}}`$ denotes the random-sample average.

Let us try to transcribe the situation realized by the Hamiltonian (1) in a manner appealing our intuition more vividly: Hard-core bosons are confined in a square lattice. Throughout this paper, we fix the particle density to be half-filled. In the absence of $`B`$ and $`\mathrm{\Delta }`$, a statistical mechanical theorem ensures that the particles are superconducting (namely, the gauge symmetry is broken spontaneously) at the ground state ; note that because of the hard-core restriction, the onset of superconductivity is a subtle issue. The random chemical potential $`\{H_i\}`$ may work so as to disturb the superconductivity. The strength of the randomness is preferably slightly less than that of the localization-transition threshold, and so the bosons are left still superconducting. Because the horizontal periodic boundary condition is imposed, the lattice is rolled up to form a cylinder. From each end of the cylinder, two bar magnets are inserted so that the cylinder is subjected to magnetic field perpendicular to the cylinder surface. This magnetic field is supposed to destroy the superconductivity so that Anderson localization may occur at a certain critical magnetic field eventually. This magnetic-field-swept transition has not yet been simulated, whereas the randomness-driven transition has already been simulated extensively . For a summary of preceding simulation results at $`B=0`$, readers may consult with Introduction of ref. .

It is meaningful to transform the boson Hamiltonian (1) into the form of quantum $`XX`$ model. With use of the mapping relations between hard-core boson and $`S=1/2`$ spin ,

$`b_i^{}`$ $`=`$ $`S_i^+`$
$`b_i`$ $`=`$ $`S_i^{},`$ (2)
the above Hamiltonian (1) is expressed in terms of the $`S=1/2`$ spin operators,

$$=\frac{J}{2}\underset{ij}{}\left(\mathrm{e}^{\mathrm{i}\varphi _{ij}}S_i^+S_j^{}+\mathrm{e}^{\mathrm{i}\varphi _{ij}}S_i^{}S_j^+\right)+2\underset{i}{}H_iS_i^z,$$
(3)
apart from a constant term. Now, we are in a position to address some notions gained from this spin representation: The onset of superconductivity, that is, the spontaneous breaking of gauge symmetry, is interpreted as the appearance of the in-plane spontaneous magnetization; that is,

$$M_{XY}^2=\underset{ij}{}\left(S_i^xS_j^x+S_i^yS_j^y\right).$$
(4)
(In the above, the spin-spin correlations over the identical sites are subtracted, because they just yield a constant term. That subtraction improves finite-size-scaling behavior.) It is proved rigorously  that the in-plane symmetry is spontaneous broken at the ground state of the two-dimensional $`XX`$ model. Therefore, it is expected that at least in the vicinity of $`B=0`$ and $`\mathrm{\Delta }=0`$, our system (1) may continue to be in a superconductor phase. The perturbations of $`B`$ and $`\mathrm{\Delta }`$ are supposed to destroy the superconductivity. In the spin language, $`B`$ introduces in-plane magnetic frustration at each plaquette; experts at the quantum Monte-Carlo technique may be convinced now why the method does not work. On the other hand, $`\mathrm{\Delta }`$ appears to be random magnetic field coupling the $`z`$-component of spin. In the next section, we will study how $`M_{XY}`$ gets disturbed by these perturbations numerically.

## 3 Numerical results

In this section, we present numerical results. We have performed exact-diagonalization calculation for the Hamiltonian (3) with system sizes up to $`L=5`$. We have fixed the particle density to be half-filled $`n(=N/L^2)=0.5`$. For those system sizes of odd $`L`$, we carried out two sets of simulations for the particle numbers of $`N=[L^2/2]`$ and $`[L^2/2]1`$ (the bracket $`[\mathrm{}]`$ denotes the Gauss notation). Thereby, we performed interpolation calculation (least-square fitting) with respect to those data so as to obtain the data at $`n=0.5`$. In the interpolation, we used the relation $`Q(n)=a(n0.5)^2+c`$, which means that the physics is symmetric with respect to $`n=0.5`$ (owing to the particle-hole symmetry). The amount of randomness $`\mathrm{\Delta }`$ is kept unchanged throughout the simulation; namely, $`\mathrm{\Delta }=0.7`$.

### 3.1 Field-tuned localization

In this subsection, we will show evidences of the superconductor-insulator transition tuned by the magnetic field. In the analysis, the spin language (3) is used. Namely, the localization transition is identified as the disappearance of the $`XY`$ (in-plane) magnetic order $`M_{XY}`$.

In Fig. 1, we plotted the square of the in-plane magnetization per site $`m^2=[M_{XY}^2]_{\mathrm{av}}/L^4`$ for $`\mathrm{\Delta }=0.7`$ and various $`B`$. The bracket $`[\mathrm{}]_{\mathrm{av}}`$ denotes the random-sample average. $`\mathrm{}`$ denotes the ground-state expectation for respective random samples. The random-sample numbers are $`1024`$, $`1024`$, $`1024`$ and $`192`$ for $`L=2`$, $`3`$, $`4`$ and $`5`$, respectively. As is noted above, for odd $`L`$, we have performed two sets of simulations. Hence, the number of random samples is effectively twice as many as that indicated above for odd $`L`$. From the plot, we see that the in-plane magnetization gets suppressed by $`B`$. However, onset of localization transition seems to be less transparent. Indication of the localization transition will be seized in the following analyses.

In Fig. 2, we plotted the Binder parameter $`U`$  for the in-plane magnetic order; that is,

$$U=1\frac{[M_{XY}^4]_{\mathrm{av}}}{3[M_{XY}^2]_{\mathrm{av}}^2}.$$
(5)
Finite-size-scaling behavior of the Binder parameter contains information whether $`M_{XY}`$ develops or not: As $`L`$ is enlarged, the Binder parameter should grow (be suppressed), provided that the order $`M_{XY}`$ is long-ranged (short ranged). At critical point, it remains scale-invariant with respect to $`L`$. That is, intersection point of the Binder-parameter curves yields the location of the transition point.

From the plot, we see that an intersection point locates at $`B0.1`$. Namely, the superconductor phase survives up to the critical magnetic field $`B0.1`$, at which the field-tuned localization transition takes place eventually. However, rather large corrections to finite-size scaling, especially for small $`L`$, prevent precise determination of the critical point.

We carried out an alternative analysis so as to complement the above observation. We calculated the gauge stiffness $`\mathrm{{\rm Y}}`$ , which is defined by the elasticity with respect to the boundary gauge twist;

$$\mathrm{{\rm Y}}=\left[\frac{^2E_\mathrm{g}(\mathrm{\Theta })}{\mathrm{\Theta }^2}|_{\mathrm{\Theta }=0}\right]_{\mathrm{av}},$$
(6)
where $`\mathrm{\Theta }`$ denotes the gauge-twist angle through the periodic boundary, and $`E_\mathrm{g}(\mathrm{\Theta })`$ is the ground-state energy under the twisted boundary condition $`\mathrm{\Theta }`$. In other words, $`\mathrm{{\rm Y}}`$ is the elasticity modulus against the introduction of a magnetic flux penetrating through the cylinder; note that our square lattice is rolled up to form a cylinder. Therefore, the stiffness should remain finite in the superconductor phase, while it vanishes in the Anderson-localization phase. In Fig. 3, we plotted the scaled stiffness $`L^2\mathrm{{\rm Y}}`$ for the same parameter range as that of Figs. 1 and 2. We observe that the stiffness is suppressed by $`B`$ as would be expected. In order to gain further information from this scaling plot, we have to consult with the formal scaling argument on $`\mathrm{{\rm Y}}`$ . According to this argument, at the localization point, the stiffness should obey the scaling relation $`\mathrm{{\rm Y}}L^{(d+z2)}`$. That is, the scaled stiffness $`L^{d+z2}\mathrm{{\rm Y}}`$ should be scale-invariant at the localization transition point. Actually, from the plot, we see that the scaled stiffness with the particular choice of $`z=2`$ exhibits scale invariance at critical point $`B0.1`$; the validity of $`z=2`$ will be confirmed in the analysis of Section 3.3. Putting together the finite-size-scaling behaviors of $`U`$ and $`\mathrm{{\rm Y}}`$, we are led to the conclusion that the field-tuned localization takes place at $`B0.1`$. In the following, we calculate experimentally accessible quantities at the point $`B=0.1`$.

### 3.2 Critical conductivity

Here, we evaluate the AC conductivity at the localization point $`B=0.1`$ determined in the above subsection. As is mentioned in Introduction, an argument claims that the DC conductivity would be a universal value of the order $`(e^{})^2/h`$ with charge of one particle $`e^{}`$ (say, Cooper pair). Various experiments have tried to validate this issue. However, the results are still remaining controversial; see Table 1. Nevertheless, it should be noted that finite conductivity itself is a novel feature lying out of the scope of the conventional localization theory in two dimensions; see Introduction.

According to the Kubo formula, the AC conductivity is expressed in terms of the current-current correlation function,

$$\sigma (\omega )=\mathrm{Re}\left[\frac{1}{\mathrm{}\omega L^2}_0^{\mathrm{}}dt\mathrm{e}^{\mathrm{i}\omega t}[J_x(t),J_x]\right]_{\mathrm{av}},$$
(7)
with the current operator $`J_x=\frac{\mathrm{i}e^{}J}{2\mathrm{}}_{j,\delta _x}\mathrm{e}^{\mathrm{i}\varphi _{j,j+\delta _x}}\delta _xa_{j+\delta _x}^{}a_j`$. This formula reduces to the resolvent form,

$$\sigma (\omega )=\mathrm{Re}\left(\frac{\mathrm{i}}{\mathrm{}\omega L^2}\left[J_x\left(\frac{\mathrm{}}{E_\mathrm{g}+\mathrm{}\omega +\mathrm{i}\eta }+\frac{\mathrm{}}{E_\mathrm{g}\mathrm{}\omega \mathrm{i}\eta }\right)J_x\right]_{\mathrm{av}}\right).$$
(8)
Hence, one is forced to calculate the inverse of the Hamiltonian matrix, which is seemingly impossible in practice. However, Gagliano and Balseiro found that the resolvent is expanded into the following continued-fraction form ,

$$f_0\left|\frac{1}{z}\right|f_0=\frac{f_0|f_0}{z\alpha _0\frac{\beta _1^2}{z\alpha _1\frac{\beta _2^2}{\mathrm{}}}},$$
(9)
with the coefficients generated by the following recursion relations,

$`|f_{i+1}`$ $`=`$ $`|f_i\alpha _i|f_i\beta _i^2|f_{i1},`$
$`\alpha _i`$ $`=`$ $`f_i||f_i/f_i|f_i,`$
$`\beta _i^2`$ $`=`$ $`f_i|f_i/f_{i1}|f_{i1}(\beta _0=0).`$ (10)
These procedures are essentially the same as those of the Lanczos tri-diagonalization. It is one of major advantages of the Lanczos diagonalization that one can calculate dynamical response functions.

In Fig. 4, we plotted the AC conductivity at the transition point ($`B=0.1`$ and $`\mathrm{\Delta }=0.7`$); we had set $`\eta =0.2`$. The data are averaged over random samples of $`1024`$, $`1024`$ and $`64`$ for $`L=3`$, $`4`$ and $`5`$, respectively. Because of finite-size energy gap above the ground state, the AC conductivity in Fig. 4 drops in the vicinity of $`\omega =0`$. AC conductivity curve forms a peak beside $`\omega =0`$, and as $`L`$ is enlarged, the peak position shifts toward the center $`\omega =0`$. Hence, for sufficiently large system sizes, the drop of $`\sigma (\omega 0)`$ may disappear, and a Drude-like peak centered at $`\omega =0`$ may emerge instead. Therefore, we regard the maximal conductivity beside $`\omega =0`$ as the DC conductivity for respective $`L`$. Since the DC conductivity for each $`L`$ exhibits finite-size correction, we need to extrapolate those finite-$`L`$ data to the value of thermodynamic limit. In Fig. 5, we plotted the finite-$`L`$ DC conductivity against $`1/L^2`$; the choice of this abscissa scale is due to the reasoning addressed in the literature . Through the least-square fit, we obtained the DC conductivity in the thermodynamic limit such as $`\sigma _\mathrm{c}=0.196\pm 0.001((e^{})^2/h)`$. Remarkably enough, we notice that this conductivity is larger than that of the randomness-driven transition $`\sigma _\mathrm{c}0.135((e^{})^2/h)`$ at $`B=0`$ , suggesting breakdown of the universality of $`\sigma _\mathrm{c}`$. We will discuss this point in the last section.

### 3.3 Dynamical critical exponent

In the subsection 3.1, we have already used the relation $`z=2`$ for the purpose of investigating scaling property of the gauge stiffness $`\mathrm{{\rm Y}}`$ (6). In this subsection, we obtain more conclusive estimate of $`z`$ with use of the Rieger-Young relation . In Fig. 6, we show the probability distribution of the logarithm of the first energy gap $`\mathrm{log}\mathrm{\Delta }E`$ over $`512`$ random samples with $`L=5`$ at the transition point ($`B=0.1`$ and $`\mathrm{\Delta }=0.7`$). According to the Rieger-Young argument, the low-energy tail of the distribution contains an information of the dynamical critical exponent. In the following, we recollect their argument briefly. Because our system is disordered, the first energy gap $`\mathrm{\Delta }E`$ distributes obeying a certain probability-distribution function $`P`$. We are considering critical phenomenon, where the physics is scale-invariant and no such particular energy scale that characterize the physics exists. Therefore, it is sensible to consider $`P(\mathrm{log}\mathrm{\Delta }E)`$ rather than $`P(\mathrm{\Delta }E)`$. Low-lying excitation may be created at a peculiar part of random sample, where local excitation costs very little energy. Therefore, the probability may be proportional to the spatial volume; namely, $`PL^d\mathrm{\Delta }E^\lambda `$ with the exponent $`\lambda `$ describing the low-energy tail. On the other hand, the finite-size-scaling theory insists that it should be of the form $`P=\stackrel{~}{P}(L/\xi _\mathrm{r})=\stackrel{~}{P}(L/\xi _\tau ^{1/z})`$ (definition of $`z`$). With use of the relation $`\mathrm{\Delta }E1/\xi _\tau `$, the distribution turns out to be a function of $`L(\mathrm{\Delta }E)^{1/z}`$. Therefore, we obtain the relation,

$$\mathrm{log}P(\mathrm{log}\mathrm{\Delta }E)=\lambda \mathrm{log}\mathrm{\Delta }E+C,$$
(11)
with $`\lambda =d/z`$.

In Fig. 6, we observe that the low-energy tail is, in fact, governed by that scaling relation with the exponent $`\lambda =1`$. In consequence, we obtain the relation $`z=d=2`$, confirming the preliminary estimate obtained in the subsection 3.1.

## 4 Summary and discussions

We have simulated the two-dimensional hard-core-boson model (1) subjected to the external magnetic field and the random chemical potential $`\mathrm{\Delta }=0.7`$. Difficulty arising from the complex hopping amplitudes $`\{\mathrm{e}^{\mathrm{i}\varphi _{ij}}\}`$, that have been preventing the application of the quantum Monte-Carlo method, is circumvented by the use of the exact-diagonalization scheme. For the first time, in computer simulation, we observed evidences of the field-tuned localization transition with the finite-size-scaling technique applied to the gauge order $`M_{XY}`$ as well as the gauge stiffness $`\mathrm{{\rm Y}}`$. After determining the location of the localization point, we carried out calculations of two experimentally accessible quantities, namely, the critical conductivity and the dynamical critical exponent. Thereby, we obtained those estimates such as $`\sigma _\mathrm{c}=0.196\pm 0.01((e^{})^2/h)`$ and $`z=2`$. First, let us make a comparison with preceding simulation results. Although no former results for $`B0`$ are available, there have been reported a number of simulation studies at $`B=0`$ . These simulation results are settling gradually to the conclusions such as $`\sigma _\mathrm{c}0.135((e^{})^2/h)`$ and $`z=2`$. Hence, contrary to common belief, we are led to conclude that the critical conductivity is not universal but increases as the external magnetic field is strengthened. As for $`z`$, on the other hand, we found that the dynamical critical exponent $`z=2`$ is retained even for $`B0`$.

Secondly, let us turn our attention to making a comparison with experiments. As is noted above, the quantities evaluated here are measurable in experiments; a summary of experimental results is presented in Table 1. Among various reports listed there, in particular, we would like to draw reader’s attention to the latest exhaustive measurement by Marković et al. , who scanned both magnetic-field strength and film thickness (randomness). Their experiment demonstrates that the critical conductivity (resistivity) does not stay universal but increases (decreases) gradually as $`B`$ is applied. Our numerical result supports this behavior . Non-universality of $`\sigma _\mathrm{c}`$ has been arousing much attention. Meanwhile, there have been proposed several ideas intending to account for that ; for instance, it was proposed that one should take into account dissipative modes (possibly un-paired normal electrons) that would give rise to extra resistivity in addition to the intrinsic one. Our result shows, however, that without resorting to such dissipative extra modes, one can account for the gradual increase of $`\sigma _\mathrm{c}`$ within the framework of the pure-boson model (1). As for $`z`$, it appears that the numerical simulation result $`z=2`$ is not in agreement with experiments. The discrepancy had already come out in the past studies of $`B=0`$. There had been proposed an attempt at altering the universality class ; It was claimed that an inclusion of long-range repulsion among particles would alter the universality class. Because in our exact-diagonalization calculation, tractable system sizes might not be sufficient to implement long-range interaction, it would be exceedingly cumbersome to validate this issue definitely.

In preliminary stage of our simulation, we had swept various parameter ranges so as to find optimal condition to observe localization transition: Note that $`\mathrm{\Delta }`$ should not exceed a threshold at which particles localize already at zero magnetic field $`B=0`$. On the other hand, $`\mathrm{\Delta }`$ should not be too weak, because superconductor-insulator transition may occur at strong magnetic field ($`B0.5`$). For such $`B`$ in close proximity to some primal fractional numbers, certain specific type of vortex-lattice structure becomes favored so that simulation data with sweeping magnetic field suffer from insystematic behaviors. Moreover, for large $`B`$, one cannot judge at all whether superconductivity is destroyed solely by the magnetic field or assistance of randomness is also significant. After sweeping various parameter ranges, we had determined the optimal randomness $`\mathrm{\Delta }=0.7`$, and thereby, we performed extensive simulations at this particular condition. It would be remained for future study to gain inclusive features for the whole region including strong magnetic field.

## Acknowledgment

Numerical calculation was performed on Alpha workstations of theoretical physics group, Okayama university, and on the vector supercomputers Hitachi SR8000/60, of the institute of solid state physics, university of Tokyo.


text:
# Analytical model for a crossover between uncorrelated and fractal behaviour of a self-repulsive chain

## Abstract

The thermodynamics of a long self-repulsive chain is studied. In $`D<4`$ dimensions it shows two distinctly different regimes, corresponding to weak and strong correlations in the system. A simple microscopic analytical model is presented which successfully describes both the regimes. The self-consistent scheme is used, in which the center of mass of a chain is fixed explicitly. This allows to take correlations into account in an indirect manner.

Problem of a crossover between the strong-correlation regime and weak correlations is of a general interest for the solid-state physics. A textbook example is here the theory of phase transitions . Not too close to the transition point Landau theory is valid, yielding the dimension-independent Landau set of critical indices. On the other hand, correlations in a critical region in a very vicinity of the transition result in different (dimension-dependent) values of critical indices, as it can be described by the renormalisation-group analysis . While the two regimes are well understood in itself, no general microscopic analytical formula describes them simultaneously.

In this paper we draw a formula of this kind, which describes the behaviour of a simpler system, namely the self-repulsive chain. The peculiarity of this system is that it does not show a phase transition, but the two regimes mentioned above are still present. We believe that an absence of a spontaneous symmetry break makes things simpler, while the physics of a crossover remains somehow similar to more complicated problems like the description of the phase transition; for instance the problem of a self-interacting chain is closely related to the theory of percolation .

Let us first define the formal (mathematical) problem. We consider paths of a random walker on a cubic lattice in $`D`$ dimensions, which returns to the starting point after $`L`$ steps (see Figure 1 for the illustration); we suppose that $`L1`$. These paths are closed polygons (maybe, self-crossing); positions of their vertexes $`x_i`$ fulfill the condition

$$x_ix_{i+1}=1;x_0=x_L.$$
(1)
We define the number of self-crosses for a given path $`X=\{x_0,x_1,\mathrm{},x_{L1}\}`$ by the formula

$$N(X)=\frac{1}{2}\underset{ij}{}\delta (x_i;x_j),$$
(2)
where $`\delta `$ is a Kronecker symbol, and introduce a statistical average as follows:

$`<f>`$ $`=`$ $`Z^1{\displaystyle \underset{X}{}}f(X)e^{\beta N(X)};`$ (3)
$`Z`$ $`=`$ $`{\displaystyle \underset{X}{}}e^{\beta N(X)}.`$ (4)
Here sums are taken over all paths $`X`$, $`f`$ is a function of $`X`$, and $`\beta `$ is a small positive parameter. In particular, we shall monitor the dispersion $`<(xu)^2>`$, where $`u=L^1_ix_i`$ is the position of the center of mass of a chain.

Physically, the above equations describe the statistics of a loop-chain with a potential energy equal to the number of its self-crosses $`N`$, at the inverse temperature $`\beta `$.

Since $`\beta 1`$, for a sufficiently small $`L`$ the statistical weight $`e^{\beta N}`$ is always close to 1, i.e. correlations in the system are negligible. Therefore we deal with a standard random-walk problem; for this case $`<(xu)^2>L`$. We shall refer to this case as an uncorrelated regime.

Let us establish a parameter which determines the magnitude of fluctuations. Since the typical volume occupied by the chain is proportional to $`L^{D/2}`$ in the uncorrelated regime, the average number of its self-crosses can be estimated as $`NL^{2D/2}`$. Therefore the dimensionless length $`\lambda =\beta ^{2/(4D)}L`$ can be introduced; the uncorrelated regime occurs at $`\beta N1`$, that corresponds to $`\lambda 1`$.

In the opposite case of a large $`\lambda `$, correlations in the chain affect its statistics strongly. The problem is similar to the well-studied case of the statistics of a self-avoiding polymer molecule . The scaling laws are dimension-dependent and can be calculated by the renormalisation-group technique. For this ”fractal” limit $`<(xu)^2>L^\varphi `$ with $`\varphi =2`$ in 1D and takes non-integer values in 2D and 3D, which are respectively close to 4/3 and 6/5. In $`D4`$ the fractal regime is virtually absent, as $`\varphi `$ takes the ”uncorrelated” value 1 (there are logarithmic corrections in $`D=4`$).

As it was declared above, the present paper is aimed to establish a formula that describes (at least qualitatively) both the uncorrelated regime and the fractal regime. The outline is as follows. First, the continuous analog (8) of the chain under study is introduced. Then, we switch to a kind of the self-consistent potential approximation (12). It yields simple formulae (16-22), that show the desired asymptotes for $`<(xu)^2>`$. The result is then compared with numerical calculations.

Let us construct a continuous analog. One can neglect the short-range correlations, as the critical behaviour in $`D<4`$ is due to the long-range ones . Let $`\tau 1`$ be a ”physically infinitesimal” length, on which correlations can be neglected; accordingly to the above estimations it should be $`\tau \beta ^{2/(4D)}`$. Since $`\beta `$ is small, in $`D<4`$ these two inequalities can be fulfilled simultaneously. Consider the paths passing through some $`x_i`$. Then the distribution function of $`x`$ at $`i+\tau `$ is given by

$$p_{i+\tau }(x;x_i)=p_0\mathrm{exp}\left(\frac{D(xx_i)^2}{2\tau }\right),$$
(5)
where $`p_0(\tau )`$ is a normalisation factor. The number of paths passing simultaneously through given $`x_i`$ and $`x_{i+\tau }`$, is obviously proportional to $`p_{i+\tau }(x_{i+\tau };x_i)`$.

Introduce also the ”partial” average number of self-crosses $`n_{ij}`$, which counts only the crosses between the ”sub-chains” of the length $`\tau `$ located near $`i`$ and $`j`$. It is defined as

$$n_{ij}(x_{i\tau }x_{j\tau },x_{i\tau }x_{(i+1)\tau },x_{j\tau }x_{(j+1)\tau })=<\underset{k,m}{}\delta (x_k;x_m)>;$$
(6)
where indices $`k`$ and $`m`$ run only within the intervals $`i\tau \mathrm{}i(\tau +1)1`$ and $`j\tau \mathrm{}j(\tau +1)1`$, respectively, whereas the averaging is carried out over the paths passing through $`x_{i\tau },x_{(i+1)\tau },x_{j\tau }`$, and $`x_{(j+1)\tau }`$.

Full path $`L`$ can be divided into $`L/\tau `$ parts of the length $`\tau `$, so that, for example, partition function $`Z`$ takes the form

$$Z=Z_0(\tau )\mathrm{exp}\left(\underset{i=0}{\overset{L/\tau }{}}\frac{Dx_{i\tau }x_{(i+1)\tau }^2}{2\tau }\frac{\beta }{2}\underset{i,j=0}{\overset{L/\tau }{}}n_{ij}\right).$$
(7)
Here the pre-exponential sum is taken over all possible sets $`\{x_0,x_\tau ,\mathrm{},x_{L\tau }\}`$.

Obviously, $`n(r,\rho _1,\rho _2)`$ peaks near $`r=0`$ and falls fast at $`r\mathrm{}`$. It can be verified  that $`_rn(r,\rho _1,\rho _2)=\tau ^2`$ at any $`\rho _1,\rho _2`$. So the continuous limit of (7) is as follows:

$$Z=Z_0[𝒟x]\mathrm{exp}\left(_0^L\frac{Dx^2(l)}{2}𝑑l\frac{\beta }{2}_0^L_0^L\delta (x(l)x(l^{}))𝑑l𝑑l^{}\right).$$
(8)
Here $`\delta `$ is the Dirac $`\delta `$-function.

Self-parallel shifts of the whole chain are irrelevant for the statistical properties under study, so we can take the path-integral only over the closed paths with the center of mass fixed at zero:

$$u=_0^Lx𝑑l=0;x(0)=x(L).$$
(9)
Continuous analog of the formula (3) for statistical averages can be written in the same manner. In particular, the distribution function of $`x`$ is equal to

$$p(x_0)=Z_p(x_0)/Z,$$
(10)
where $`Z_p(x_0)`$ is a contribution to $`Z`$ from the trajectories passing though the given point $`x_0`$ at $`l=0`$, that is $`x(0)=x_0`$. Function $`p(x_0)`$ is normalised:

$$p(x_0)𝑑x_0=1.$$
(11)
The key assumption of the present model is to approximate the second term in the exponent by the interaction with harmonic potential:

$$\frac{\beta }{2}_0^L_0^L\delta (x(l)x(l^{}))𝑑l𝑑l^{}_0^L\left(U_0+\frac{D\mathrm{\Omega }^2x^2(l)}{2}\right)𝑑l,$$
(12)
where $`\mathrm{\Omega }`$ is defined in a self-consistent way:

$$D\mathrm{\Omega }^2=\beta LD^1(^2p)|_{x_0=0}$$
(13)
(the value of $`U_0`$ may remain undefined, as it drops out of the final formulae).

After the approximation (12) is done, the distribution function can be found straightforwardly. The problem resembles Feynman’s path-integral expressions for the statistics of a quantum particle , with however an additional condition for trajectories $`x𝑑l=0`$. Let us introduce a notation

$$S_{x(l)}=_0^L\left(\frac{Dx^2}{2}U_0\frac{D\mathrm{\Omega }^2x^2}{2}\right)𝑑l$$
(14)
and define the ”minimal-action” trajectory $`\overline{x}(l;x_0)`$, which satisfy (9) and delivers a minimum of $`S`$ at given $`\overline{x}(0)=x_0`$. The conditional minimisation results in the problem of a classical oscillator in the external field

$$\frac{d^2\overline{x}}{dl^2}+\mathrm{\Omega }^2\overline{x}f=0,$$
(15)
where the value of $`f`$ should be chosen to fulfil the condition (9). The following minimal value of $`S`$ is obtained:

$`S_{min}(x_0)`$ $`=`$ $`U_0L+{\displaystyle \frac{6D\mathrm{\Lambda }(\omega )x_0^2}{L}},`$ (16)
$`\mathrm{\Lambda }(\omega )`$ $`=`$ $`{\displaystyle \frac{\omega ^2}{3}}\left({\displaystyle \frac{1}{1\omega \mathrm{ctg}\omega }}\right)`$ (17)
$`\omega `$ $`=`$ $`\mathrm{\Omega }L/2.`$ (18)
Consider now possible deviations from this trajectory: $`x(l)=\overline{x}(l)+\stackrel{~}{x}(l)`$, where $`\stackrel{~}{x}(0)=0`$ and the conditions (9) are satisfied. As the potential is harmonic it appears that $`S_x=S_{min}+S_{\stackrel{~}{x}}`$, and $`Z_p`$ can be factorised:

$$Z_p(x_0)=\stackrel{~}{Z}\mathrm{exp}(S_{min}(x_0)).$$
(19)
Here $`\stackrel{~}{Z}`$ does not depend on $`x_0`$. The last two formulae and expressions (10, 11) allow to find out the distribution function:

$$p(x_0)=\left(\frac{6D\mathrm{\Lambda }(\omega )}{2\pi L}\right)^{D/2}\mathrm{exp}\left(\frac{6D\mathrm{\Lambda }(\omega )x_0^2}{L}\right)$$
(20)
The dispersion is equal to

$$<x^2>=\frac{L}{12\mathrm{\Lambda }(\omega )}.$$
(21)
Substitution of $`^2p`$ in formula (13) gives the equation

$$\mathrm{\Lambda }(\omega )=\frac{\pi }{3D}\left(\frac{2D\omega ^2}{\pi \beta }\right)^{\frac{2}{D+2}}L^{\frac{4D}{D+2}}.$$
(22)
The left-hand and right-hand sides of this equation are schematically plotted in Figure 2 as functions of $`\omega `$. There is a single root within the range $`(0,\pi )`$. Its position is governed by $`\lambda =\beta ^{2/(4D)}L`$. At $`\lambda 1`$ we obtain $`\omega 0`$ and therefore $`\mathrm{\Lambda }(\omega )1`$, so

$$<x^2>|_{\lambda 1}=\frac{L}{12}.$$
(23)
The ”uncorrelated” dependence $`<x^2>L`$ is reproduced; moreover the factor $`1/12`$ is also correct. At the opposite limit $`\omega (\lambda 1)\pi `$, and it follows from (21-22) that

$$<x^2>|_{\lambda 1}=\frac{D}{4\pi }\left(\frac{\beta }{2\pi D}\right)^{\frac{2}{D+2}}L^{\frac{6}{D+2}}$$
(24)
In particular $`<x^2>L`$ in 4D and $`<x^2>L^2`$ in 1D. The indices in 2D and 3D are also quite reasonable. In $`D>4`$ the scheme is hardly suitable, as it is impossible to fulfill the inequality $`1\tau \beta ^{2/(4D)}`$ and the continuum analog (8) cannot be used.

It is important to discuss the following peculiarity of the approximation made. Usually, the use of mean-field schemes means that correlations like $`<x(l_1)x(l_2)>`$ at $`l_1l_2`$ are totally neglected. In our scheme however some correlations are still present, because the center of mass of the trajectory is fixed (formula (9)). Namely that circumstance results in a non-trivial dependence of $`S_{min}(L)`$ and the dimension-dependent critical index in the last formula.

We have also compared the result obtained with Monte-Carlo calculations in 1D and 2D. Figure 3 shows bi-logarithmic plots for the dispersion $`<x^2>`$ versus $`L`$ at $`\beta =10^2`$ in 1D and $`\beta =510^2`$ in 2D. Thin lines are asymptotes (23,24). An agreement is quite good even quantitatively.

The final remark is that although we have considered only a particular random-walk problem (1) with unitary steps, it is clear that the model can be applied for any random-walk self-repulsing chains. It is only required that (a) the repulsion is short-range and (b) the central limit theorem can be applied, yielding the Gaussian distribution (5). Once these requirements are fulfilled, the continuous analog of the universal form (8) can be established.

The author is grateful to P.A. Prudkovskii, A.A. Nikulin, and A.A. Fedyanin, who made valuable stylistic comments. The work was partially supported by ”Scientific Schools” program (grant 96-1596476).


text:
# String Matching in 𝑂̃(√𝑛+√𝑚) Quantum Time

## 1 Introduction

We consider the following problem: given a text $`t`$ of length $`n`$ and a pattern $`p`$ of length $`m`$, does $`p`$ occur in $`t`$? This question requires $`\mathrm{\Theta }(n+m)`$ time classically, using one of several known algorithms, e.g., the Knuth-Morris-Pratt algorithm or the Boyer-Moore algorithm. Note that the above problem is slightly different from the usual string matching problem which requires finding all occurrences of the pattern in the text.

We explore the above question on a quantum machine. Our starting point is the algorithm due to Grover\[Gr\] which searches for an element in an unordered database in $`O(\sqrt{n})`$ time. Boyer, Brassard, Hoyer and Tapp\[BBHT96\] gave a tighter analysis of Grover’s algorithm and also showed how to handle the case when the number of items searched for is unknown. Dürr and Hoyer\[DH\] used Grover’s algorithm to find the minimum in $`O(\sqrt{n})`$ time.

Finding whether the pattern matches somewhere in the text is akin to searching in an unordered database; the only issue is that checking one element of this database, i.e., a text position, for an occurrence of the pattern takes $`O(m)`$ time. In fact, this can be speeded up to $`O(\sqrt{m})`$ by viewing the act of checking whether the pattern matches at a particular text position as the act of finding a mismatch amongst $`m`$ elements; this search can be performed in $`O(\sqrt{m})`$ time with constant failure probability. This gives an overall time complexity of $`O(\sqrt{nm})`$ (actually, there are additional logarithmic factors, as will be described later).

In this paper, we show how this complexity can be improved to $`O(\sqrt{n}\mathrm{log}\sqrt{\frac{n}{m}}\mathrm{log}m+\sqrt{m}\mathrm{log}^2m)`$, by combining the above quantum search paradigm with a standard technique from parallel string matching, called Deterministic Sampling, due to Vishkin\[Vi\]. Our algorithm will work with constant failure probability. Thus, if the pattern occurs in the text, it will return some occurrence of the pattern in the text (or the leftmost occurrence, if needed) in $`O(\sqrt{n}\mathrm{log}\sqrt{\frac{n}{m}}\mathrm{log}m+\sqrt{m}\mathrm{log}^2m)`$ time, with probability which is a constant strictly more than 1/2. And if the pattern does not occur in the text then the algorithm will say so with probability which is also a constant strictly more than 1/2. The failure probability can be decreased to inverse polynomial at the expense of further logarithmic factors. Finally, note that the second component of the running time above will be due to pattern preprocessing.

## 2 Preliminaries

We use the following theorem based on Grover’s\[Gr\] database searching algorithm. The theorem itself is due to Boyer et al.\[BBHT96\].

###### Theorem 1

Given an oracle evaluating to 1 on at least $`t1`$ of the elements in an unordered database of size $`n`$, there is a quantum algorithm which returns the index of a random element on which the oracle evaluates to 1, with probability at least 3/4, in $`O(\sqrt{\frac{n}{t}})`$ time and oracle calls.

We will also need the following theorem for finding the minimum element in a database, due to Dürr and Hoyer\[DH\].

###### Theorem 2

Given a comparison oracle, there is a quantum algorithm which finds the index of the minimum element in an unordered database of size $`n`$, with probability at least 3/4, in $`O(\sqrt{n})`$ time.

We will assume a basic oracle which will compare a text and a pattern character or two pattern characters in $`O(1)`$ time. Our aim then is to use this oracle to develop suitable oracles which will enable solving the string finding problem in $`\stackrel{~}{O}(\sqrt{n}+\sqrt{m})`$ time. These new oracles could be probabilistic, i.e., they will give the correct answer with constant probability. We derive the following corollary to Theorems 1 and 2, in order to use probabilistic oracles. First, we define a probabilistic oracle formally.

Probabilistic Oracles. These are oracles which evaluate to 1 on good elements (i.e., those which are being searched for) with probability at least 3/4 and to 0 on bad elements, with probability at least 3/4. A probabilistic comparison oracle gives the correct answer with probability at least 3/4.

###### Corollary 1

Given a probabilistic oracle and a database with $`t1`$ good elements, there is a quantum algorithm which returns the index of a random good element with probability at least 3/4 in $`O(\sqrt{\frac{n}{t}}\mathrm{log}\sqrt{\frac{n}{t}})`$ time and oracle calls. Similarly, given a probabilistic comparison oracle, there is a quantum algorithm which finds the index of the minimum element with probability at least 3/4, in $`O(\sqrt{n}\mathrm{log}\sqrt{n})`$ time.

Proof. For the searching problem, each original oracle call in Theorem 1 gets replaced by $`O(\mathrm{log}\sqrt{\frac{n}{t}})`$ calls and the majority result is taken. This is to ensure that, with constant probability, none of the original $`\sqrt{\frac{n}{t}}`$ oracle calls returns a 1 on a bad element or a 0 on a good element. Similarly, for the minimum finding problem, each original oracle call in Theorem 2 gets replaced by $`O(\mathrm{log}\sqrt{n})`$ calls.

Our algorithm will use oracles whose running time will not be a constant. Therefore, the search time will be obtained by multiplying the time given by the above corollary with the time taken per oracle call.

We will require the following facts about strings.

Periodicity. A string $`p`$, $`|p|=m`$, is said to be aperiodic if any two instances of the string, one shifted to the right of the other by at most $`m/2`$, differ in some column. A string which is not aperiodic is called periodic. A periodic string $`p`$ has the form $`v^ku`$, where $`v`$ cannot be expressed as a concatenation of several instances of a smaller string, $`u`$ is a prefix of $`v`$, and $`k2`$. $`|v|`$ is said to be the period of $`p`$.

## 3  The $`\stackrel{~}{O}(\sqrt{n}\sqrt{m})`$ Time Algorithm

Consider using Grover’s algorithm in conjunction with the following natural probabilistic oracle $`f()`$ to solve the string matching problem: $`f(i)=1`$ if the pattern matches with left endpoint aligned with text position $`i`$, and $`f(i)=0`$ with probability at least 3/4, otherwise. This probabilistic oracle can be implemented so that it runs in $`\stackrel{~}{O}(\sqrt{m})`$ time as follows.

The idea is to implement the oracle $`f(i)`$ using Grover’s algorithm itself, using a deterministic oracle $`g(i,j)`$ which looks for a mismatch at location $`j`$. $`g(i,j)=1`$ if and only if the pattern with left endpoint aligned with $`t[i]`$ mismatches at location $`p[j]`$. By Theorem 1, such a location, if one exists, can be found in $`O(\sqrt{m})`$ time with probability at least 3/4. We set $`f(i)=0`$ if Grover’s algorithm with oracle $`g(i,j)`$ succeeds in finding a mismatch, and $`f(i)=1`$ otherwise. Using Corollary 1, the time taken to search using oracle $`f(i)`$ is $`O(\sqrt{n}\sqrt{m}\mathrm{log}\sqrt{n})`$ and the success probability is at least 3/4.

Next, we give the faster algorithm, first for aperiodic strings and then for periodic strings.

## 4 The $`\stackrel{~}{O}(\sqrt{n}+\sqrt{m})`$ Time Algorithm: Aperiodic Patterns

The above oracle $`f(i)`$ was too expensive to get an $`\stackrel{~}{O}(\sqrt{n}+\sqrt{m})`$ bound. A faster oracle will clearly improve the time. We do not know how to get a faster oracle directly. However, we reorganize the computation as described below and then use a faster oracle followed by a slower one, speeding up the algorithm on the whole.

We partition the text into blocks of length $`m/2`$ and use Grover’s algorithm to search for a block which contains an occurrence of the pattern. This is done using a probabilistic oracle $`h(i)`$ (to be described later), which evaluates to 1 with probability at least 3/4 if block $`i`$ has a match of the pattern (with left endpoint starting in block $`i`$), and evaluates to 0 with probability at least 3/4, otherwise. Note that by aperiodicity, the pattern can match with left endpoint at at most one text position in block $`i`$. $`h(i)`$ will take $`O(\sqrt{m}\mathrm{log}m)`$ time. It follows from Corollary 1 that the time taken for searching with the oracle $`h(i)`$ will be $`O(\sqrt{\frac{n}{m}}\sqrt{m}\mathrm{log}\sqrt{\frac{n}{m}}\mathrm{log}m=\sqrt{n}\mathrm{log}\sqrt{\frac{n}{m}}\mathrm{log}m)`$ and the success probability is a constant.

The oracle $`h(i)`$ itself will run in two steps. The first step will use deterministic sampling and will takes $`O(\sqrt{m}\mathrm{log}m)`$ time with constant success probability. This step will eliminate all but at most one of the pattern instances with left endpoint in block $`i`$. The second step will check whether this surviving instance matches the text using the $`g()`$ oracle defined in Section 3; this will take $`O(\sqrt{m})`$ time with constant success probability. We describe the two steps next.

### 4.1 Step 1: Deterministic Sampling

The oracle is based on the following theorem due to Vishkin \[Vi\].

###### Theorem 3

\[The Deterministic Sampling Theorem\] Let $`p`$ be aperiodic. Consider $`\frac{m}{2}`$ instances of $`p`$, with successive instances shifted one step to the right. Let these instances be labelled from $`1`$ to $`\frac{m}{2}`$, from left to right. Then there exists an instance $`f`$, and a set of at most $`O(\mathrm{log}m)`$ positions in $`p`$, called the deterministic sample with the following property: if all positions corresponding to the deterministic sample in instance $`f`$ of $`p`$ match the text, then none of the other instances of $`p`$ above can possibly match the text.

Proof. By aperiodicity, there is a column which contains two distinct characters and stabs all the pattern instances; pick any character in that column which is not in majority. Remove all pattern instances which do not have this character in the column being considered. Repeat $`O(\mathrm{log}m)`$ times until only one pattern instance remains. Then $`f`$ is the label of the instance which remains and the columns chosen give the deterministic sample.

Assume that a deterministic sample for the pattern has been precomputed. We will describe this precomputation later.

We now describe the first step in $`h(i)`$, where $`i`$ is a block number. We use Grover’s algorithm in conjunction with the deterministic oracle $`k(i,j)`$ which evaluates to $`1`$ if and only if the $`j`$th instance of the pattern (amongst those instances with left endpoint in block $`i`$) matches the text on its deterministic sample. Clearly, $`k(i,j)`$ takes $`O(\mathrm{log}m)`$ time. The search using $`k(i,j)`$ takes $`O(\sqrt{m}\mathrm{log}m)`$ time by Theorem 1. This search returns an instance $`j`$ of the pattern with left endpoint in block $`i`$ which has its deterministic sample matching the text (if such an instance exists), with probability at least 3/4.

### 4.2 Step 2: Direct Verification

Next, we use another application of Grover’s search to determine whether or not instance $`j`$ determined above in block $`i`$ matches the text; this is done as in Section 3 using the deterministic oracle $`g()`$. It succeeds with probability 3/4, and takes time $`O(\sqrt{m})`$.

Thus, $`h(i)`$ returns a 1 with probability at least 3/4 if block $`i`$ contains a match of the pattern, and 0 with probability at least 3/4, otherwise. The time taken to search using $`h(i)`$ is $`O(\sqrt{n}\mathrm{log}\sqrt{\frac{n}{m}}\mathrm{log}m)`$, as claimed above, and the success probability is at least 3/4.

Once a block $`i`$ is found in which $`h(i)`$ evaluates to 1, a search using oracle $`k(i,j)`$ on block $`i`$ gives the unique pattern instance $`j`$ with left endpoint in $`i`$ whose deterministic sample matches; the success probability is at least 3/4. Another search using the oracle $`g()`$ determines if this pattern instance mismatches the text; the success probability (in finding a mismatch, if any) is again at least 3/4. Thus, the total time taken is $`O(\sqrt{n}\mathrm{log}\sqrt{\frac{n}{m}}\mathrm{log}m)`$, and the probability of success is as follows.

If the pattern occurs in the text, then with probability 3/4, the search with $`h()`$ will return a block containing a match of the pattern; subsequently, with probability 3/4, the search with $`k(i,j)`$ will return a matching pattern instance, and the last search with $`g()`$ will not discover a mismatch with probability 1. Thus an occurrence of the pattern in the text will be found with $`(3/4)^2>1/2`$ probability, as claimed. And if the pattern does not occur in the text, then the last search with $`g()`$ will determine a mismatch with probability at least $`3/4>1/2`$, as required.

Finally, note that the leftmost occurrence of the pattern can be determined using the minimum finding algorithm in Corollary 1 to first find the leftmost block with $`h(i)`$ evaluating to 1, and subsequently, searching within that block as above. The time taken and the success probability are as in the previous paragraphs.

## 5 Pattern Preprocessing for Aperiodic Patterns

We show how to determine the deterministic sample in $`O(\sqrt{m}\mathrm{log}^2m)`$ time.

Determining the Deterministic Sample. Imagine $`m/2`$ copies of the pattern placed as in Theorem 3. Determining the sample will proceed in $`O(\mathrm{log}m)`$ stages. In each stage, some column and a character in that column will be identified; all surviving pattern copies which do not have this character in this column will be eliminated from future stages. This will continue until only one pattern copy remains uneliminated. Each stage will take $`O(\sqrt{m}\mathrm{log}m)`$ time and will have a constant success probability (where we count a success if the surviving pattern copies halve in cardinality).

A stage proceeds as follows. First, a column containing two distinct characters amongst the surviving pattern copies is found, with constant probability. This column will also have the property that all surviving pattern copies are stabbed by it. Two distinct characters in the above column are also identified. One of these two characters is chosen at random as the next character in the sample. Clearly, the number of stages is $`O(\mathrm{log}m)`$ with inverse polynomial (in $`m`$) failure probability because the probability of the number of surviving pattern instances halving in a stage is at least a constant.

It remains to describe how a column containing two distinct characters amongst the surviving pattern copies is found, with constant probability. Before describing this we need to mention how to find the leftmost and rightmost surviving patterns in a stage. This is done using the minimum finding algorithm of Dürr and Hoyer in conjunction with an oracle which indicates which pattern copies are consistent with the already chosen deterministic sample points; this oracle takes $`O(\mathrm{log}m)`$ time per call, giving an $`O(\sqrt{m}\mathrm{log}m)`$ time algorithm for finding the leftmost/rightmost surviving pattern copy, by Theorem 2. The success probability is at least 3/4. Now, we can describe the algorithm for finding a column with two distinct characters.

First, the leftmost and rightmost surviving pattern copies are found as above. Then a column in which these two copies differ is found using Grover’s algorithm in conjunction with a suitable oracle in $`O(\sqrt{m})`$ time; this step succeeds with probability at least 3/4. Given a column, this oracle determines whether or not the two pattern copies above differ in this column. By Theorem 1, searching for a column with two distinct characters using this oracle takes $`O(\sqrt{m})`$ time and succeeds with probability 3/4.

Thus, in time $`O(\sqrt{m}\mathrm{log}m)`$, a column containing two distinct characters and stabbing all surviving pattern copies is found, with constant probability; it is easily seen that two distinct characters in this column are also found in this process.

The total time taken in determining the deterministic sample is thus $`O(\sqrt{m}\mathrm{log}^2m)`$.

## 6 Handling Periodic Patterns

We sketch briefly the changes required to the above algorithm in order to handle periodic pattern.

For periodic patterns, the above preprocessing algorithm will not terminate with a single pattern copy but rather with several copies shifted $`|v|`$ steps to the right successively. When a stage is reached when the only surviving copies are the periodically shifted copies above, then the search for a heterogeneous column in the next $`\mathrm{\Theta }(\mathrm{log}m)`$ stages will fail. Note that for aperiodic patterns this behaviour happens with low, i.e., inverse polynomial probability.

At this point, we determine the period $`|v|`$ using two instances of the minimum finding algorithm. The first instance finds the leftmost surviving copy and the second the second leftmost; the difference of their offsets is the period. This takes $`O(\sqrt{m}\mathrm{log}m)`$ time, using the oracle which checks for consistency with the deterministic sample and also compares offsets. Given the period $`|v|`$, the following changes now need to be made to the text processing part.

Recall the oracle $`h(i)`$ from Section 4; this oracle determines whether there is a pattern instance with left endpoint in block $`i`$ which matches, first on its deterministic sample, and then on the whole. This oracle is modified as follows.

$`h(i)`$ will first determine the leftmost and the rightmost pattern instances with left endpoints in $`i`$ which match on their respective deterministic sample points; this is done using the minimum finding algorithm and takes $`O(\sqrt{m}\mathrm{log}m)`$ time with success probability at least 3/4 (see Theorem 2, this success probability can be made arbitrarily close to 1 by repeating). Let these two instance have left endpoints at text positions $`k`$ and $`l`$ respectively.

Next, $`h(i)`$ finds the longest substring (with length at most $`m`$) starting at the right boundary of text block $`i`$ which is consistent with the pattern instance starting at text position $`l`$ (and therefore consistent with the pattern instant starting at text position $`k`$ as well); this is done using the minimum finding algorithm and takes $`O(\sqrt{m})`$ time with constant failure probability. Similarly, $`h(i)`$ finds the longest substring (with length at most $`m/2`$) ending at the right boundary of text block $`i`$ which is consistent with the pattern instance starting at text position $`k`$.

Finally, using these two substrings, $`h(i)`$ can determine in $`O(1)`$ time, whether there exists a pattern instance with left endpoint in block $`i`$ which matches the text. If the length of the two substrings is less than $`m`$ then there is no such pattern instance; otherwise, all instances of the pattern which occur completely within these two substrings and starting at shifts of integer multiples of $`|v|`$ from $`k`$ are complete matches (here $`|v|`$ is the pattern period).

Thus, $`h(i)`$ determines whether or not the pattern occurs in block $`i`$ in $`O(\sqrt{m}\mathrm{log}m)`$ time, with failure probability a constant. This failure probability can be made arbitrarily close to 0 by repetition. Note that $`h(i)`$ can determine the leftmost pattern occurrence in block $`i`$ as well, if required, within the same time bounds.

The rest of the algorithm stays the same: $`h(i)`$ is used to find a block containing an occurrence of the pattern and subsequently, an occurrence of the pattern in this block is found using the above method.

## 7 Conclusions And Open Problems

We have shown how one occurrence or the leftmost occurrence of $`p`$ in $`t`$ can be found in $`\stackrel{~}{O}(\sqrt{n}+\sqrt{m})`$ time, with constant two-sided failure probability. We also note that an approximate count of the number of occurrences (within a multiplicative constant factor) can also be determined in $`\stackrel{~}{O}(\sqrt{n}+\sqrt{m})`$ using the approximate counting algorithm of Brassard, Hoyer and Tapp \[BHT\], adapted appropriately (the oracle $`h(i)`$ must now return a count of the number of matches rather than just the indication of a match). Finally, using the same algorithm, the total number of occurrences of $`p`$ can be determined in $`\stackrel{~}{O}(\sqrt{nt}+\sqrt{m})`$ time, where $`t`$ is the number of occurrences.

One open problem would be whether string matching with don’t cares can be performed in the same time bounds as above. The main challenge here to implement convolution using Fast Fourier Transforms in $`\stackrel{~}{O}(\sqrt{n})`$ time. It is not obvious how this can be accomplished.


text:
# CMB Analysis of Boomerang & Maxima & the Cosmic Parameters {Ω_{𝑡𝑜𝑡},Ω_𝑏h²,Ω_{𝑐𝑑𝑚}h²,Ω_Λ,𝑛_𝑠}

## 1. CMB Analysis of Primary Anisotropies

Experiments and Bandpowers: Anisotropies at the $`30\mu K`$ level at low multipoles revealed by COBE in 1992 were augmented at higher $`\mathrm{}`$ in some 19 other experiments, some with a comparable number of resolution elements to the 600 or so for COBE, most with many fewer. A list of these experiments to April 1999 with associated bandpowers is given in Bond, Jaffe and Knox (2000 \[BJK00\]). The anisotropy picture dramatically improved this past year, as results were announced first in summer 99 from the ground-based TOCO experiment in Chile (Miller et al. 2000), then in November 99 from Boomerang-NA, the North American test flight (Mauskopf et 1999). These two additions improved peak localization and gave evidence for $`\mathrm{\Omega }_{tot}1`$. Then in April 2000, results from the first CMB long duration balloon (LDB) flight, were announced (de Bernardis et al. 2000), followed in May 2000 by results from the night flight of Maxima (Hanany et al. 2000). Boomerang’s best resolution was $`10^{}`$, about 40 times better than that of COBE, with tens of thousands of resolution elements. Maxima had a similar resolution but covered an order of magnitude less sky. Fig. 1 shows the 150A GHz Boomerang-LDB map and the Wiener-filtered Maxima-1, to scale. The de Bernardis et al. (2000) maps at 90 and 220 GHz show the same spatial features as this 150 GHz one, with the overall intensities falling precisely on the CMB blackbody curve. The Toco, Boomerang and Maxima experiments are described elsewhere in these proceedings. They were designed to reveal the primary anisotropies of the CMB, those which can be calculated using linear perturbation theory. Fig. 1 shows the temperature power spectra for Boomerang, Maxima and prior-CMB data (Boomerang-NA+TOCO+April 99) are in good agreement. Sketching the impact of these new results on cosmic parameter estimation (Lange et al. 2000 \[Let00\], Jaffe et al. 2000 \[Jet00\]) is the goal of this paper. Space constraints preclude adequate referencing here, but these are given in the Boomerang (Let00) and Maxima+Boomerang (Jet00) parameter estimation papers (see also Bond 1996, \[B96\], for other references).

We are only at the beginning of the high precision CMB era for primary anisotropies heralded by the arrival of Boomerang and Maxima, with interferometers taking data (VSA, CBI, DASI), the single dish ACBAR about to, and new LDBs to fly in the next few years (Arkeops, Tophat, Beast/Boost), as well as Boomerang-2001 and the neo-Maxima Maxipol, both concentrating on polarization. In April 2001, NASA’s HEMT-based MAP satellite will launch, with $`12^{}`$ resolution, and in 2007, ESA’s bolometer+HEMT-based Planck satellite is scheduled for launch, with $`5^{}`$ resolution.

The CMB Analysis Pipeline: Analyzing Boomerang and other experiments involves a pipeline that takes (1) the timestream in each of the bolometer channels coming from the balloon plus information on where it is pointing and turns it into (2) spatial maps for each frequency characterized by average temperature fluctuation values in each pixel (Fig. 1) and a pixel-pixel correlation matrix characterizing the noise, from which various statistical quantities are derived, in particular (3) the temperature power spectrum as a function of multipole (Fig. 1), grouped into bands, and two band-band error matrices which together determine the full likelihood distribution of the bandpowers (Bond, Jaffe & Knox 1998 \[BJK98\], BJK00). Fundamental to the first step is the extraction of the sky signal from the noise, using the only information we have, the pointing matrix mapping a bit in time onto a pixel position on the sky. To compare the data with millions of cosmological models, as we wish to do here, the radical compression step from 2 to 3 is essential, and hinges upon an accurate representation of the likelihood surface.

There is generally another step in between (2) and (3), namely separating the multifrequency spatial maps into the physical components on the sky: the primary CMB, the thermal and kinematic Sunyaev-Zeldovich effects, the dust, synchrotron and bremsstrahlung Galactic signals, the extragalactic radio and submillimetre sources. The strong agreement among the Boomerang maps indicates that to first order we can ignore this step, but it has to be taken into account as the precision increases. The Fig. 1 map is consistent with a Gaussian distribution, thus fully characterized by just the power spectrum. Higher order (concentration) statistics (3,4-point functions, etc.) tell us of non-Gaussian aspects, necessarily expected from the Galactic foreground and extragalactic source signals, but possible even in the early Universe fluctuations. For example, though non-Gaussianity occurs only in the more baroque inflation models of quantum noise, it is a necessary outcome of defect-driven models of structure formation. (Peaks compatible with Fig. 1 do not appear in non-baroque defect models, which now appear unlikely.) Though great strides have been made in the analysis of Boomerang and Maxima, there is intense needed effort worldwide now to develop new fast algorithms to deal with the looming megapixel datasets of LDBs and the satellites (e.g., Bond et al. 1999, Szapudi et al. 2000).

## 2. Cosmic Parameters

Parameters of Structure Formation: We usually adopt the restricted set of 7 cosmological parameters used in Let00 and Jet00, $`\{\mathrm{\Omega }_\mathrm{\Lambda },\mathrm{\Omega }_k,\omega _b,\omega _{cdm},n_s,\tau _C,\sigma _8\}`$. The curvature energy is $`\mathrm{\Omega }_k1\mathrm{\Omega }_{tot}`$. The dark energy parameterized here by $`\mathrm{\Omega }_\mathrm{\Lambda }`$ could have complex dynamics associated with it, e.g., if it is the energy density of a scalar field which dominates at late times (now often termed a quintessence field, $`Q`$, with energy $`\mathrm{\Omega }_Q`$, e.g., Steinhardt 2000). One popular phenomenology is to add one more parameter, $`w_Q=p_Q/\rho _Q`$, where $`p_Q`$ and $`\rho _Q`$ are the pressure and density of the $`Q`$-field. Thus $`w_Q=1`$ and $`\mathrm{\Omega }_Q=\mathrm{\Omega }_\mathrm{\Lambda }`$ for the cosmological constant. We have also allowed $`w_Q`$ to float.

We use 2 parameters to characterize the early universe primordial power spectrum of gravitational potential fluctuations $`\mathrm{\Phi }`$, one giving the overall power spectrum amplitude $`𝒫_\mathrm{\Phi }(k_n)`$, and one defining the shape, a spectral tilt $`n_s(k_n)1+d\mathrm{ln}𝒫_\mathrm{\Phi }/d\mathrm{ln}k`$, both at some (comoving) normalization wavenumber $`k_n`$. We really need another 2, $`𝒫_{GW}(k_n)`$ and $`n_t(k_n)`$, associated with the gravitational wave component. In inflation, the amplitude ratio is related to $`n_t`$ to lowest order, with $`𝒪(n_sn_t)`$ corrections at higher order, e.g., B96. There are also useful limiting cases for the $`n_sn_t`$ relation. However, as one allows the baroqueness of the inflation models to increase, one can entertain a plethora of power spectra (with fully $`k`$-dependent $`n_s(k)`$ and $`n_t(k)`$) if one is artful enough in designing inflaton potential surfaces. As well, one can have more types of modes present, e.g., scalar isocurvature modes ($`𝒫_{is}(k_n),n_{is}(k)`$) in addition to, or in place of, the scalar curvature modes ($`𝒫_\mathrm{\Phi }(k_n),n_s(k)`$). However, our philosophy is consider minimal models first, then see how progressive relaxation of the constraints on the inflation models, at the expense of increasing baroqueness, causes the parameter errors to open up. For example, with COBE-DMR and Boomerang, we can probe the GW contribution, but the data are not powerful enough to determine much. Planck can in principle probe the gravity wave contribution reasonably well.

We use another 2 parameters to characterize the transport of the radiation through the era of photon decoupling, which is sensitive to the physical density of the various species of particles present then, $`\omega _j\mathrm{\Omega }_j\mathrm{h}^2`$. We really need 4: $`\omega _b`$ for the baryons, $`\omega _{cdm}`$ for the cold dark matter, $`\omega _{hdm}`$ for the hot dark matter (massive but light neutrinos), and $`\omega _{er}`$ for the relativistic particles present at that time (photons, very light neutrinos, and possibly weakly interacting products of late time particle decays). For simplicity, though, we restrict ourselves to the conventional 3 species of relativistic neutrinos plus photons, with $`\omega _{er}`$ therefore fixed by the CMB temperature and the relationship between the neutrino and photon temperatures determined by the extra photon entropy accompanying $`e^+e^{}`$ annihilation. Of particular importance for the pattern of the radiation is the (comoving) distance sound can have influenced by recombination (at redshift $`z_r=a_r^11`$),

$$r_s=\frac{6000}{\sqrt{3}}\mathrm{Mpc}_0^{\sqrt{a_r}}\frac{d\sqrt{a}}{(\omega _m+\omega _{er}a^1)^{1/2}(1+\omega _ba/(4\omega _\gamma /3))^{1/2}},$$
(1)
where $`\omega _\gamma =2.46\times 10^5`$ is the photon density, $`\omega _{er}=1.68\omega _\gamma `$ for 3 species of massless neutrinos and $`\omega _m\omega _{hdm}+\omega _{cdm}+\omega _b`$.

The angular diameter distance is

$`=\{d_k\mathrm{sinh}(\chi _r/d_k),\chi _r,d_k\mathrm{sin}(\chi _r/d_k)\},\mathrm{where}d_k={\displaystyle \frac{3000}{\sqrt{|\omega _k|}}}\mathrm{Mpc},`$ (2)
$`\chi _r=6000\mathrm{Mpc}{\displaystyle _{\sqrt{a_r}}^1}{\displaystyle \frac{d\sqrt{a}}{(\omega _m+\omega _Qa^{6w_Q}+\omega _ka)^{1/2}}}.`$
The 3 cases are for negative, zero and positive mean curvature, $`d_k`$ is the curvature scale and $`\chi _r`$ is the comoving distance to recombination. The location of the first acoustic peak, $`L_{Pk}0.746\pi /r_s`$ (e.g., Efstathiou and Bond 1999, hereafter EB99), depends upon $`\omega _b`$ through the sound speed as well as on $`\omega _k`$, $`\omega _\mathrm{\Lambda }`$ and $`\omega _m`$. Thus $`L_{Pk}`$ defines a functional relationship among these parameters, a degeneracy (EB99) that would be exact except for the integrated Sachs-Wolfe effect, associated with the change of $`\mathrm{\Phi }`$ with time if $`\mathrm{\Omega }_\mathrm{\Lambda }`$ or $`\mathrm{\Omega }_k`$ is nonzero.

Our 7th parameter is an astrophysical one, the Compton ”optical depth” $`\tau _C`$ from a reionization redshift $`z_{reh}`$ to the present. It lowers $`𝒞_{\mathrm{}}`$ by $`\mathrm{exp}(2\tau _C)`$ at the high $`\mathrm{}`$’s probed by Boomerang. For typical models of hierarchical structure formation, we expect $`\tau _C\mathrm{}<0.2`$. It is partly degenerate with $`\sigma _8`$ and cannot be determined at this precision by CMB data now.

The LSS also depends upon our parameter set. Here we use a set of (relatively weak) constraints on $`\mathrm{ln}\sigma _8^2`$ from cluster abundance data and on $`\mathrm{\Gamma }+(n_s1)/2`$ from galaxy clustering data (B96, Let00). $`\sigma _8^2`$ is a bandpower for density fluctuations on a scale associated with rare clusters of galaxies, $`8\mathrm{h}^1\mathrm{Mpc}`$, which we often use in place of $`𝒫_\mathrm{\Phi }(k_n)`$ or $`𝒞_{10}`$ for the amplitude parameter. The mass-density-power-spectrum-shape-parameter $`\mathrm{\Gamma }`$ depends upon $`\{\omega _m,\omega _{er},\omega _b,\mathrm{h}\}`$ and is related to the horizon scale when the energy density in relativistic particles equals that in nonrelativistic ones.

When we allow for freedom in $`\omega _{er}`$, the abundance of primordial helium, tilts of tilts ($`dn_{\{s,is,t\}}(k_n)/d\mathrm{ln}k,\mathrm{}`$) for 3 types of perturbations, the parameter count would be 17, and many more if we open up full theoretical freedom in spectral shapes. However, as we shall see, as of now only 3 or 4 combinations can be determined with 10% accuracy with the CMB. Thus choosing 7 is adequate for the present, 6 of which are discretely sampled, with generous boundaries.111The specific discrete parameter values used for the $`𝒞_{\mathrm{}}`$-database in this analysis were: ($`\mathrm{\Omega }_\mathrm{\Lambda }=`$ 0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1.0,1.1), ($`\mathrm{\Omega }_k=`$ .9,.7,.5,.3,.2,.15,.1,.05,0,-.05,-.1,-.15,-.2,-.3,-.5), ($`\tau _c=`$0, .025, .05, .075, .1, .15, .2, .3, .5) when $`w_Q=`$ –1, with slightly different ranges (and $`\mathrm{\Omega }_{tot}`$=1) when we allow $`w_Q`$ to float; ($`\omega _c=`$ .03, .06, .12, .17, .22, .27, .33, .40, .55, .8), ($`\omega _b=`$ .003125, .00625, .0125, .0175, .020, .025, .030, .035, .04, .05, .075, .10, .15, .2), ($`n_s=`$1.5, 1.45, 1.4, 1.35, 1.3, 1.25, 1.2, 1.175, 1.15, 1.125, 1.1, 1.075, 1.05, 1.025, 1.0, .975, .95, .925, .9, .875, .85, .825, .8, .775, .75, .725, .7, .65, .6, .55, .5), $`\sigma _8^2`$ was continuous, and there were 4 experimental parameters for Boomerang and Maxima (calibration and beam uncertainties), as well as other calibration parameters for some of the prior-CMB experiments. For drawing cosmological conclusions we adopt a weak prior probability on the Hubble parameter and age: we restrict $`\mathrm{h}`$ to lie in the 0.45 to 0.9 range, and the age to be above 10 Gyr.

The First Peak, $`\mathrm{\Omega }_{tot}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$: $`L_{Pk}`$ and its errors are found from the average and variance of $`\mathrm{ln}L_{Pk}`$, taken wrt the full probability function over our database described above (restricted for this exercise to the $`\tau _C=0`$ part and with the weak prior). As more CMB data were added, $`L_{Pk}`$ evolved and the errors shrunk considerably: $`240_{54}^{+70}`$ for April 99 data, $`220_{27}^{+30}`$ for TOCO+4.99 data, $`224_{21}^{+23}`$ when Boomerang-NA was added, and $`212_7^{+7}`$ when Boomerang-LDB and Maxima-1 were added to prior-CMB. The latter contrasts with $`202_7^{+7}`$ for Boomerang-LDB alone, $`226_{16}^{+17}`$ for Maxima-1 alone, and $`208_7^{+7}`$ for the combination. The numbers change a bit depending upon exactly what prior one chooses or what functional forms one averages over. For the choice here, although there is a large difference in the mean between the Maxima and Boomerang numbers, it is not unreasonable within the errors. Other ways of doing this make the discrepancy seem more statistically significant (e.g., Page 2000, these proceedings).

In Fig. 2, we show the lines of constant $`L_{Pk}/r_s`$ in the $`\mathrm{\Omega }_{tot}`$$`\mathrm{\Omega }_\mathrm{\Lambda }`$ plane, for given $`\omega _m`$ and $`\omega _b`$, using the formulas given above and discussed in more detail in EB99. The $`\pm 10`$ band around 210 corresponds to our best $`L_{Pk}`$ estimate of $`212\pm 7`$ using all current CMB data. Note that the constant $`L_{Pk}`$ lines look rather similar to the contours shown in the right panel, showing that the $`/r_s`$ degeneracy plays a large role in determining the contours. The contours hug the $`\mathrm{\Omega }_{tot}=1`$ line more closely than the allowed $`L_{Pk}`$ band does for the maximum probability values of $`\omega _m`$ and $`\omega _b`$, because of the shift in the allowed $`L_{Pk}`$ band as $`\omega _m`$ and $`\omega _b`$ vary in this plane. See also Bond et al. (2001b, \[capp2K\]) for $`L_{Pk}`$ plots in the quintessence plane, $`w_Q`$-$`\mathrm{\Omega }_Q`$, which demonstrate why $`w_Q`$ is poorly determined by CMB alone.

Marginalized Estimates of our Basic 7 Parameters: Table 1 shows there are strong detections with only the CMB data for $`\mathrm{\Omega }_{tot}`$, $`\omega _b`$ and $`n_s`$ in the minimal inflation-based 7 parameter set, and a reasonable detection of $`\omega _{cdm}`$. The ranges quoted are Bayesian 50% values and the errors are 1-sigma, obtained after projecting (marginalizing) over all other parameters. That $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is not well determined is a manifestation of the $`\mathrm{\Omega }_{tot}`$$`\mathrm{\Omega }_\mathrm{\Lambda }`$ near-degeneracy discussed above, which is broken when LSS is added because the CMB-normalized $`\sigma _8`$ is quite different for open cf. $`\mathrm{\Lambda }`$ models. Supernova at high redshift give complementary information to the CMB, but with CMB+LSS (and the inflation-based paradigm) we do not need it: the CMB+SN1 and CMB+LSS numbers are quite compatible. In our space, the Hubble parameter, $`\mathrm{h}=(_j(\mathrm{\Omega }_j\mathrm{h}^2))^{1/2}`$, and the age of the Universe, $`t_0`$, are derived functions of the $`\mathrm{\Omega }_j\mathrm{h}^2`$: representative values are given in the Table caption.

Fig. 3 shows how the parameter estimations evolved as more CMB data were added (for the weak+LSS prior). With just the COBE-DMR+LSS data, the 2-sigma contours were already localized in $`\omega _{cdm}`$. Without LSS, it took the addition of Maxima-1 before it began to localize. $`\mathrm{\Omega }_k`$ localized near zero when TOCO was added to the April 99 data, more so when Boomerang-NA was added, and much more so when Boomerang-LDB and Maxima-1 were added. Some $`n_s`$ localization occurred with just ”prior-CMB” data. $`\omega _b`$ really focussed in with Boomerang-LDB and Maxima-1, as did $`\mathrm{\Omega }_\mathrm{\Lambda }`$.

We have also considered what happens as we let $`\mathrm{\Omega }_{m\nu }/\mathrm{\Omega }_m`$, the fraction of the matter in massive neutrinos, vary from 0 to 0.3, for LSS + all of the CMB data and $`\mathrm{\Omega }_{tot}`$=1 (Bond et al. 2001a, \[$`\nu `$2K\]). Until Planck precision, the CMB data by itself will not be able to strongly discriminate this ratio. Adding HDM does have a strong impact on the CMB-normalized $`\sigma _8`$ and the shape of the density power spectrum (effective $`\mathrm{\Gamma }`$ parameter), both of which mean that when LSS is included, adding some HDM to CDM is strongly preferred in the absence of $`\mathrm{\Omega }_\mathrm{\Lambda }`$. However, though higher $`\mathrm{\Omega }_m`$ is preferred at the expense of less dark energy, significant $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is still required (see $`\nu `$2K for the evolution of the CMB+LSS 2-sigma contours in the $`(\omega _{hdm}+\omega _{cdm})`$$`\mathrm{\Omega }_\mathrm{\Lambda }`$ plane as $`\mathrm{\Omega }_{m\nu }/\mathrm{\Omega }_m`$ is varied). The $`\omega _b`$ and $`n_s`$ likelihood curves are essentially independent of $`\mathrm{\Omega }_{m\nu }/\mathrm{\Omega }_m`$.

The Future, Forecasts for Parameter Eigenmodes: We can also forecast dramatically improved precision with further analysis of Boomerang and Maxima, future LDBs, MAP and Planck. Because there are correlations among the physical variables we wish to determine, including a number of near-degeneracies beyond that for $`\mathrm{\Omega }_{tot}`$$`\mathrm{\Omega }_\mathrm{\Lambda }`$ (EB99), it is useful to disentangle them, by making combinations which diagonalize the error correlation matrix, ”parameter eigenmodes” (e.g., B96, EB99). For this exercise, we will add $`\omega _{hdm}`$ and $`n_t`$ to our parameter mix, but set $`w_Q`$=$`1`$, making 9. (The ratio $`𝒫_{GW}(k_n)/𝒫_\mathrm{\Phi }(k_n)`$ is treated as fixed by $`n_t`$, a reasonably accurate inflation theory result.) The forecast for Boomerang+DMR based on the 440 square degree patch with a single 150 GHz bolometer used in the published data is 3 out of 9 linear combinations should be determined to $`\pm 0.1`$ accuracy. This is indeed what we get in the full analysis of Let00 for CMB only. If 4 of the 6 150 GHz channels are used and the region is doubled in size, we predict 4/9 could be determined to $`\pm 0.1`$ accuracy. The Boomerang team is still working on the data to realize this promise. And if the optimistic case for all the proposed LDBs is assumed, 6/9 parameter combinations could be determined to $`\pm 0.1`$ accuracy, 2/9 to $`\pm 0.01`$ accuracy. The situation improves for the satellite experiments: for MAP, we forecast 6/9 combos to $`\pm 0.1`$ accuracy, 3/9 to $`\pm 0.01`$ accuracy; for Planck, 7/9 to $`\pm 0.1`$ accuracy, 5/9 to $`\pm 0.01`$ accuracy. While we can expect systematic errors to loom as the real arbiter of accuracy, the clear forecast is for a very rosy decade of high precision CMB cosmology that we are now fully into.