format: text

text:
# Symmetry breaking and the random-phase approximation in small quantum dots

## I Introduction

The random-phase approximation (RPA) is a very successful many-body theory, widely used to describe Fermi systems in general. It was originally proposed by Bohm and Pines to describe the response function of the electron gas and, since then, it has been rederived in many different ways and contexts Pines . Important developments of this approximation for finite systems, including the establishment of a formal scheme for both excited and ground states, were obtained many years ago in nuclear physics Rowe . Unfortunately, great technical difficulties have hindered for a long time the solution of the general RPA equations, imposing the recourse to severe simplifying assumptions. In recent times, however, a renewed interest in the RPA as a systematic theory of ground state correlations in atomic nuclei has appeared (see for instance Ref. Nuclear, ).

In the electron-gas theory the RPA has played a benchmark role in the development of accurate exchange-correlation energy functionals, specially in the high density limit Langr . When turning to finite electronic systems, RPA calculations are generally more involved, due to the lack of translational symmetry. Nevertheless, extensive applications to metallic clusters (see, e.g., the review by Brack Brack ) and recent studies of molecules have been performed Furche .

A convenient starting point to treat Fermion systems is, in many cases, a mean field description like the Hartree-Fock (HF) approach. Selfconsistency between the mean field and the single-particle orbitals and total energy minimization are the basic conditions at this level. It may happen that the selfconsistent solution breaks one of the symmetries of the exact Hamiltonian, a well known phenomenon called spontaneous symmetry breaking. It has been predicted that the RPA should improve the ground state energy. More importantly, it should restore symmetries which could be broken at the mean field level. It is our aim in this work to show with a realistic calculation for interacting electrons that this is indeed the case. We have chosen as an appropriate system a two-dimensional parabolic quantum dot, partly because of the interest attracted by semiconductor nanostructures as a paradigm of artificial quantum systems, and, also, because these systems are believed to be good candidates where electron localization or formation of Wigner molecules should be relevant in some regimes (see Ref. Reimann,  for a recent review). We remind that Wigner molecules in quantum dots are associated with the breaking of rotational symmetry of the mean field.

The RPA is, in principle, formulated to address the excited vibrational states, or vibrons for short, of a Fermion system by means of a superposition of particle-hole and hole-particle transitions on an a priori unknown correlated ground state. A basic assumption consists in treating each particle-hole pair as a boson, i.e., the so-called quasi-boson approximation Rowe ; RS ; BR86 . Thus defined, the RPA spectrum describes the intrinsic vibrations of the system and, imposing the condition of excitation vacuum, also yields an approximation for the correlated ground state. In cases of mean field broken symmetry a remarkable property of the RPA is that it separates out the collective excitation associated with each broken symmetry as a zero-energy or spurious RPA mode, with an inertial parameter which is also fixed by the approximation. As we shall show below, these spurious RPA modes play a central role in the restoration of symmetries broken at the mean field level. The physical excitations of the system consist, therefore, of vibrons and collective modes associated with spurious solutions, both kinds being described within RPA.

It should be noted that the restoration of broken symmetries can also be attained via projection techniques (see Refs. RS ; BR86 ). Examples demonstrating their use for the case of quantum dots have been recently presented by Yannouleas and Landman Yan02 . In this paper, however, we shall focuss our analysis on the RPA description of the ground state and its connection to the excited states.

The electronic islands formed at the GaAs-AlGaAs interface by using appropriate electrodes constitute a remarkable example of controllable quantum system. Many properties can be described, in a first approach, by using a model of electrons confined by an external potential well of parabolic type and restricting the electronic motion to the interface plane. These model quantum dots have been intensively studied using a variety of methods, including semiclassical, mean field and density-functional approaches as well as exact methods attempting a direct solution of the Schrödinger equation Reimann . Within the latter category we find the diagonalization in a basis, the coordinate-space and the Monte-Carlo solutions. In the exact methods one normally imposes from the beginning the Hamiltonian symmetries to the solution and, therefore, the symmetry breaking is not possible (unless degeneracy is present). This striking difference with the mean-field solution does not necessarily imply that the latter one is an artifact. Indeed, the RPA provides a consistent scheme for the corresponding physical interpretation. We attempt in this work a clarification of this issue by comparing in a test case, such as small quantum dots, exact, mean-field and RPA results.

The paper is organized as follows: Section II briefly introduces the HF mean field approach and discusses the appearance of symmetry-broken solutions. Section III is devoted to the RPA excitations, introducing the quasi-boson approximation and discussing the RPA spurious mode connected with a global rotation when the electrons are localized in 2-, 6- and 12-electron dots. In Sec. IV we discuss the RPA ground state focussing on the correlation energy and symmetry restoration, comparing the RPA results with exact values. The conclusions are finally drawn in Sec. V.

## II The Mean-Field solution

We consider a system of $`N`$ electrons whose motion is restricted to the $`xy`$ plane, where a parabolic potential induces the formation of an electron island. The full Hamiltonian reads

$$=\underset{i=1}{\overset{N}{}}\left[\frac{𝐩^2}{2m}+\frac{1}{2}m\omega _0^2r^2\right]_i+\underset{i>j=1}{\overset{N}{}}\frac{e^2}{\kappa r_{ij}},$$
(1)
where $`\kappa `$ and $`m`$ are the dielectric constant and electron’s effective mass; $`\omega _0`$ is the external confinement frequency and we have used polar coordinates ($`r^2=x^2+y^2`$). By introducing the length $`\mathrm{}_0=(\mathrm{}/m\omega _0)^{1/2}`$ and energy $`E_0=\mathrm{}^2/m\mathrm{}_0^2=\mathrm{}\omega _0`$ units we may rewrite Hamiltonian (1) in terms of a single adimensional parameter

$$R_W=\frac{e^2/(\kappa \mathrm{}_0)}{\mathrm{}\omega _0}$$
(2)
as

$$=\underset{i=1}{\overset{N}{}}\left[\frac{1}{2}^2+\frac{1}{2}r^2\right]_i+R_W\underset{i>j=1}{\overset{N}{}}\frac{1}{r_{ij}}.$$
(3)
The $`R_W`$ parameter was introduced in Ref. rw,  as a measure of the relative importance of electron-electron interaction to confinement potential strength. In Ref. YL1,  it was used as a measure to define the different phases in quantum dots and, in particular, a formation of Wigner molecules. In our unit system it is also fulfilled that $`\mathrm{}=m=1`$.

### II.1 The Hartree-Fock equations

Introducing the basis of oscillator and spin states that diagonalizes the first term in Eq. (3), $`\{|a\eta ;a=1,\mathrm{}𝒩;\eta =,\}`$, where $`a`$ labels the orbital part and $`\eta `$ the spin, an arbitrary single-particle orbital $`|i`$ is expanded as

$$|i=\underset{a\eta }{}B_{a\eta }^{(i)}|a\eta .$$
(4)
The oscillator states $`|a`$ are characterized by radial ($`n_a`$) and angular momentum ($`m_a`$) quantum numbers –see, e.g., Ref. Hawrylak, . We shall assume that each orbital $`i`$ has non zero components only for a given spin orientation $`\eta _i=`$ or $`\eta _i=`$, i.e., $`B_{a\eta }^{(i)}=\delta _{\eta \eta _i}B_{a\eta _i}^{(i)}`$. This condition imposes good spin symmetry on the single-particle orbitals and it is the only symmetry that we shall keep, leaving totally unspecified the remaining spatial symmetries. Note that the Slater determinant built with these single-particle orbitals will be an eigenstate of the total $`S_z`$ operator but not, in general, of $`𝐒^2`$ not2 . In the chosen basis, the HF equations are written as a system of nonlinear eigenvalue equations for the matrix of $`B`$ coefficients. The corresponding equation for the $`a`$-th component of orbital $`i`$ reads

$`\epsilon _a^{(0)}B_{a\eta _i}^{(i)}`$ $`+`$ $`{\displaystyle \underset{c=1}{\overset{𝒩}{}}}B_{c\eta _i}^{(i)}[{\displaystyle \underset{bd=1}{\overset{𝒩}{}}}V_{abcd}\left({\displaystyle \underset{k=1}{\overset{N}{}}}B_{b\eta _k}^{(k)}B_{d\eta _k}^{(k)}\right)`$ (5)
$`{\displaystyle \underset{bd=1}{\overset{𝒩}{}}}V_{abdc}\left({\displaystyle \underset{k=1}{\overset{N}{}}}\delta _{\eta _i\eta _k}B_{b\eta _k}^{(k)}B_{d\eta _k}^{(k)}\right)]`$
$`=`$ $`\epsilon _iB_{a\eta _i}^{(i)}.`$
In these equations $`\epsilon _a^{(0)}`$ is the energy of the oscillator state $`|a`$ and, analogously, $`V_{abcd}`$ and $`V_{abdc}`$ are the matrix elements of the two-body interaction defined above ($`R_W/r_{12}`$) with the corresponding oscillator states. We can take advantage of the fact that these matrix elements can be obtained analytically in terms of the radial ($`n_a`$) and angular momentum ($`m_a`$) quantum numbers for each basis state $`|a`$.

The first and second terms within the square brackets of Eq. (5) are the well-known direct and exchange contributions to the HF problem, respectively not1 . We note that the direct term is independent of spin, while the exchange yields a different contribution for $`\eta _i=`$ and $`\eta _i=`$. Using the selfconsistency approach, both direct and exchange terms are determined from the ’preceding’ iteration and then the eigenvalue problem can be separately solved for each spin subset. Successive repetitions of the process lead eventually to the converged solution. It is worth to stress that this formulation of the HF equations is particularly well suited to numerical treatment, since the use of analytical expressionsChak  for $`V_{abcd}`$ minimizes numerical errors. The only remaining source of numerical error is, in fact, the precission of the eigenvalue solver. Nevertheless, a necessary approximation of the method is the truncation of the basis size $`𝒩`$ and one must check the stability of the results with this parameter.

### II.2 Some HF results

The oscillator basis is defined by including all states below a given cut-off $`E_c`$, i.e., $`\epsilon _a^{(0)}<E_c`$. The value of $`E_c`$ has been chosen large enough to satisfy the above mentioned convergence criterion. In particular, for $`N=2`$ and 6 it sufficed to include the first 55 and 153 oscillator states, respectively. For $`N=12`$ we used up to 190 states which gave good HF convergence. With the purpose to present in the following sections a systematic comparison of HF, RPA and exact results with varying $`R_W`$ we have taken as test case the $`N=2`$ quantum dot computing its properties for $`R_W=0.5`$ to 2.5 in steps of 0.5. Besides, some selected cases for $`N=6`$ and 12 have also been calculated.

The HF density

$$\rho _{\mathrm{HF}}(𝐫)=\underset{i,\mathrm{𝑜𝑐𝑐}}{}|\phi _i(𝐫)|^2,$$
(6)
where $`\phi _i(𝐫)𝐫|i`$ refers to the HF wave functions (4), are shown in Fig. 1 for the $`N=2`$ system. We note that between $`R_W=1`$ and 1.5 the system changes from a circularly symmetric solution to a spontaneously broken one, where the two electrons localize in opposite positions in the mean-field frame. These intrinsic localization becomes more and more conspicuous as the $`R_W`$ parameter is increased; an evident manifestation that if electron repulsion is strong enough the favoured solution consists of particles located as far as possible from each other.

In Fig. 2 we display the density for the $`N=6`$ and 12 quantum dots, using $`R_W=1.89`$. This value has been chosen after Pederiva et al. Pederiva  who, in turn, adjusted it to the experiments by Tarucha et al.Tarucha . Both cases show a clear symmetry breaking; the $`N=6`$ having the electrons localized on a ring while the 12-electron dot has a central electron dimer surrounded by a ring with 10 electrons. Regarding the HF total energies, we leave the discussion for the following sections where it will be compared with RPA and exact results.

## III RPA excitations

This section presents a brief discussion of the main points of the RPA for excitations, which are an essential ingredient for the ground state description in the same theory.

### III.1 The Quasi-Boson approximation

It is assumed that the system’s excitations are created by operators of the type

$$O_\lambda ^{}=\underset{mi}{}\left(X_{mi}^{(\lambda )}a_m^{}a_iY_{mi}^{(\lambda )}a_i^{}a_m\right),$$
(7)
whose action on the as yet unknown ground state $`|0`$ yields the excited vibrational states $`|\lambda `$ ($`\lambda >0`$); i.e.,

$$|\lambda =O_\lambda ^{}|0.$$
(8)
Hereafter, we shall refer to the electronic vibrational states as vibrons. Using the standard notation we refer to particle states (unoccupied HF orbitals) by subscripts $`m`$, $`n`$ and to hole states (below the HF Fermi level) by $`i`$, $`j`$. We shall denote the HF Slater determinant for $`N`$ electrons as $`|v`$ since it will act as our particle-hole vacuum. The coefficients $`X^{(\lambda )}`$ and $`Y^{(\lambda )}`$ in Eq. (7) are a set of complex amplitudes characterizing each particular excitation.

The quasi-boson approximation (QBA) amounts to treat each particle-hole pair ($`mi`$) as an elementary boson, thus fulfilling

$$[a_i^{}a_m,a_n^{}a_j][b_{mi},b_{nj}^{}]=\delta _{mn}\delta _{ij},$$
(9)
where the first equality introduces the boson operators $`b_{mi}^{}=a_m^{}a_i`$. The result (9) is exact in the average, i.e., when considered as expectation value with the uncorrelated vacuum $`|v`$ as $`v|[a_i^{}a_m,a_n^{}a_j]|v=\delta _{mn}\delta _{ij}`$; see Refs. Rowe, ; RS, ; BR86, . In this approximation, the creation $`O_\lambda ^{}`$ and annihilation $`O_\lambda `$ operators obey bosonic commutation relations

$$[O_\lambda ,O_\lambda ^{}^{}]=\delta _{\lambda \lambda ^{}}.$$
(10)
At the RPA level the eigenstates and eigenvalues are obtained by simplifying the exact equations of motion

$$[H,O_\lambda ^{}]=\mathrm{}\omega _\lambda O_\lambda ^{}$$
(11)
for the vibron operators $`O_\lambda ^{}`$ and excitation energies $`\mathrm{}\omega _\lambda `$ within the particle-hole space. With the aid of the QBA the particle-hole part of the commutator $`[H,O_\lambda ^{}]`$ can be derived Rowe ; RS ; BR86  and Eq. (11) then transforms to a generalized eigenvalue problem, well known in the literature. We stress that henceforth all commutators are to be understood in the above QBA. The RPA equations read

$$\left(\begin{array}{cc}A\hfill & B\hfill \\ B^{}\hfill & A^{}\hfill \end{array}\right)\left(\begin{array}{c}\hfill X\\ \hfill Y\end{array}\right)=\mathrm{}\omega _\lambda \left(\begin{array}{c}\hfill X\\ \hfill Y\end{array}\right),$$
(12)
where the submatrices $`A`$ and $`B`$ are given by

$`A_{minj}`$ $`=`$ $`\epsilon _{mi}\delta _{mn}\delta _{ij}+V_{mjin}V_{mjni}`$
$`B_{minj}`$ $`=`$ $`V_{mnij}V_{mnji}.`$ (13)
Here, $`\epsilon _{mi}`$ is the energy of each bosonic pair, given in terms of the HF eigenvalues as $`\epsilon _{mi}=\epsilon _m\epsilon _i`$.

The matrix elements in Eqs. (III.1) may be computed in the oscillator basis using expressions similar to those of the HF eigenvalue problem. For instance,

$`V_{mjin}`$ $`=`$ $`\delta _{\eta _m\eta _i}\delta _{\eta _n\eta _j}\times `$ (14)
$`{\displaystyle \underset{abcd=1}{\overset{𝒩}{}}}V_{abcd}B_{a\eta _m}^{(m)}B_{b\eta _j}^{(j)}B_{c\eta _i}^{(i)}B_{d\eta _n}^{(n)}.`$
In practice, the numerical solution of Eq. (12) requires the diagonalization of a rather large matrix, whose dimension depends on the number of hole $`N_h`$ and particle $`N_p`$ states. To be consistent with our imposed spin symmetry we restrict to bosonic pairs where each member of the pair has the same spin.

### III.2 Spurious modes

It was proved by Thouless Th  that, when the HF solution corresponds to a minimum in the energy surface, the RPA Eq. (12) provides only real frequencies $`\omega _\lambda `$. Furthemore, if the full Hamiltonian has any symmetry which is broken by the mean field there is a corresponding RPA mode with zero frequency (a spurious mode) which is orthogonal to the other vibrational (nonzero energy) modes. In other words, the generators of symmetries broken on the mean field level create eigenstates with zero energy in RPA.

As shown below, we shall be concerned with rotational symmetry for which the infinitesimal generator is the orbital angular momentum $`L_z`$. To treat the spurious mode related to rotation it is convenient to introduce the canonical conjugate operators $`L_z`$ and an angle operator $`\mathrm{\Phi }`$; the latter defined by the following relations Th ; MW69

$`[H,L_z]`$ $`=`$ $`0,`$
$`[H,C]`$ $`=`$ $`\mathrm{}L_z`$
$`[L_z,C]`$ $`=`$ $`\mathrm{}𝒥.`$ (15)
where $`C=i𝒥\mathrm{\Phi }`$ is an anti-Hermitian operator and $`𝒥`$ is the moment of inertia.

When the above Eqs. are implemented in the RPA they lead to the so-called Thouless-Valatin moment of inertiaTV62  $`𝒥_0`$ and to the two RPA vectors

$`L_z`$ $``$ $`{\displaystyle \underset{mi}{}}\left(\mathrm{}_{mi}^{(z)}b_{mi}^{}+\mathrm{}_{mi}^{(z)}b_{mi}\right),`$
$`C`$ $``$ $`{\displaystyle \underset{mi}{}}\left(c_{mi}b_{mi}^{}c_{mi}^{}b_{mi}\right).`$ (16)
The coefficients of the RPA $`L_z`$ operator are directly given by the single-particle HF matrix elements. In contrast, for the $`C`$ operator one needs to solve the linear system of equations

$$\left(\begin{array}{cc}A\hfill & B\hfill \\ B^{}\hfill & A^{}\hfill \end{array}\right)\left(\begin{array}{c}c\\ c^{}\end{array}\right)=\mathrm{}\left(\begin{array}{c}\mathrm{}^{(z)}\\ \mathrm{}^{(z)}\end{array}\right).$$
(17)
Once these two sets of coefficients are determined, the Thouless-Valatin moment of inertia $`𝒥_0`$ may be calculated as

$$𝒥_0=\underset{mi}{}\left(\mathrm{}_{mi}^{(z)}c_{mi}+\mathrm{}_{mi}^{(z)}c_{mi}^{}\right).$$
(18)
The operators of Eqs. (16), together with the RPA vibrons $`O_\lambda ^{}`$, $`O_\lambda `$, form a complete set for any operator linear in the bosons $`b_{mi}^{}`$ and $`b_{mi}`$. From the point of view of its Fermionic character this could represent an arbitrary one-particle-one-hole operator of the type

$$F=\underset{mi}{}f_{mi}b_{mi}^{}+f_{im}b_{mi}.$$
(19)
The corresponding expansion in terms of RPA excitations (for a single spurious mode) reads

$`F`$ $`=`$ $`{\displaystyle \underset{\omega _\lambda >0}{}}\left([O_\lambda ,F]O_\lambda ^{}+[F,O_\lambda ^{}]O_\lambda \right)`$ (20)
$`+`$ $`{\displaystyle \frac{1}{𝒥_0}}\left([F,C]L_z+[L_z,F]C\right).`$
When applied to the boson operators themselves, the above result implies that we can expand each $`b_{mi}^{}`$ and $`b_{mi}`$ in terms of the vibrons ($`O_\lambda ^{}`$, $`O_\lambda `$) and spurious modes ($`L_z`$,$`C`$). We shall return to this point in more detail in the next section.

### III.3 Strength function and sum rules

The strength function for a general Hermitian single-particle operator $`F`$ reads

$$S(E)=\underset{\lambda }{}\delta (E\mathrm{}\omega _\lambda )\left|\lambda |F|0\right|^2.$$
(21)
Important quantities related to the strength function are its energy moments, usually known as sum-rules. The $`k`$-th moment $`S_k`$ is given by

$$S_k=\underset{\lambda }{}(\mathrm{}\omega _\lambda )^k\left|\lambda |F|0\right|^2.$$
(22)
It is well known that the RPA fulfills several exact sum rules (see Ref. Lippa,  for recent applications of sum rules to quantum dots). In particular, for the dipole operators $`F=_ix_i`$ or $`F=_iy_i`$ one has $`S_1=\mathrm{}^2N/(2m)`$. More importantly, the RPA also fulfills the Kohn theorem for parabolic confinement (see, e.g., Ref. Chak, ) stating that the only allowed dipole excitation is the center of mass rigid motion at an energy $`\mathrm{}\omega _0`$. These exact results constitute good checks of selfconsistency between HF and RPA, as well as of convergence with space size, in numerical calculations.

### III.4 Calculated spectra

Figure 3 displays the eigenvalues $`\mathrm{}\omega _\lambda `$ for two cases of the $`N=2`$ quantum dot as a function of the number of particle states included in the RPA equations. The first case ($`R_W=1`$) is circularly symmetric at the mean field level while the second one ($`R_W=2`$) has rotational broken symmetry (cf. Fig. 1). The evolution of the eigenvalue set with $`N_p`$ shows a remarkable difference in these two cases. While all the eigenvalues of the circular dot stabilize for a high enough value of $`N_p`$, the symmetry broken solution exhibits one state whose energy keeps decreasing as the number of particle states is increased. It clearly corresponds to the appearance of an RPA spurious mode connected with the broken rotational symmetry, which in the limit $`N_p\mathrm{}`$ should lie at zero energy. Though stable, the first excited state of the circular case lies rather low in energy, a hint of the proximity in $`R_W`$ value to the broken symmetry transition point. Actually, for $`R_W=0.5`$ the first vibron settles down at a higher energy, $`E0.7E_0`$.

Figure 4 shows the corresponding spectra for the $`N=6`$ and 12 quantum dots. Both HF solutions correspond to broken symmetry cases (cf. Fig. 2) and, therefore, their RPA spectra should display spurious solutions. This is clearly the situation for $`N=6`$, with a well separated low energy mode. For $`N=12`$, however, the situation is less convincing and we actually obtain two modes which seem to separate from the higher vibron states. We stress that for the 12 electron dot the strong increase in number of particle-hole pairs prevents us from extending the calculation to higher $`N_p`$ values. Nevertheless, a tentative interpretation may be given taking into account the mean field structure of this dot shown in Fig. 2. One would expect a real spurious mode, corresponding to the rigid rotation of the full system, and a soft (low energy) mode corresponding to the rotation of the inner electron dimer with respect to the outer ring of electrons. These two modes could thus explain the behaviour of the two lower RPA eigenvalues in the right panel of Fig. 4.

We analyze next the strength function results. Figure 5 displays some representative cases, for the dot with six electrons. The upper panels show the electric ($`x`$) and magnetic dipole ($`\mathrm{}_z`$) strength functions computed within the HF method, while the lower ones are the analogous results within the RPA. In HF all excitations are above the gap between the Fermi level and the first unoccupied orbital (the HOMO-LUMO gap) which for $`N=6`$ takes the value $`1.8E_0`$. In RPA, however, one obtains states below this gap. Actually, the strength function for the $`\mathrm{}_z`$ operator concentrates almost exclusively on the first RPA eigenvalue, at a very small energy; a clear manifestation of rotational spurious mode discussed above. For the dipole operator there is again a single mode, this time at an energy $`1E_0`$, and the total sum rule $`S1`$ is $`3E_0\mathrm{}_0^2`$. As also discussed above, this proves that our results fulfill Kohn theorem with an excellent accuracy.

The Thouless-Valatin moments of inertia obtained for the quantum dots with spurious modes are summarized in Tab. I. These inertial parameters define physical rotational bands, characterized by the energies $`\mathrm{}^2M^2/(2𝒥_0)`$ with $`M=0,1,\mathrm{}`$ not3 . We have also checked numerically that the Thoules-Valatin moment of inertia coincides, to a good accuracy, with the value obtained from a constrained HF calculation for $`=\lambda L_z`$ as

$$𝒥_0=\frac{d^2}{d\lambda ^2}=\frac{dL_z}{d\lambda }.$$
(23)
We stress that the equivalence between the two moments of inertia can be fulfilled if and only if a self-consistent HF minimum solution is found.

## IV The ground state in the RPA

The RPA ground state $`|0`$ is defined from the condition that it is the vacuum for all vibrons

$$O_\lambda |0=0.$$
(24)
In this section we shall discuss the following aspects of the RPA ground state: i) total energy; ii) restoration of the HF broken symmetries.

### IV.1 Correlation energy

The total energy in the RPA ($`E_{\mathrm{RPA}}`$) can be split in a mean field contribution ($`E_{\mathrm{HF}}`$) and a correction ($`\mathrm{\Delta }_{\mathrm{RPA}}`$), as

$$E_{\mathrm{RPA}}=E_{\mathrm{HF}}\mathrm{\Delta }_{\mathrm{RPA}}.$$
(25)
Note that the correction $`\mathrm{\Delta }_{\mathrm{RPA}}`$, with a minus sign, gives the standard correlation energy of the system. A decrease in ground state energy ($`\mathrm{\Delta }_{\mathrm{RPA}}>0`$) will imply an improvement with respect to the mean field theory. However, since RPA is not based on energy minimization, it is not bound to fulfill the variational principle and, therefore, $`E_{\mathrm{RPA}}`$ could be even lower than the exact energy. In other words, the correlation energy could be overestimated in RPA, as it has been suggested in the literature BR86 .

The RPA correction reads RS ; BR86

$$\mathrm{\Delta }_{\mathrm{RPA}}=\frac{1}{2}\left(\mathrm{Tr}A\underset{\omega _\lambda >0}{}\mathrm{}\omega _\lambda \right).$$
(26)
The above Eq. (26) includes the contribution from the vibrons at positive frequencies $`\omega _\lambda `$ and, also, from the spurious mode associated with a rotational excitation. Note that Eq. (26) is the result of a partial cancellation between two large terms. In practice, this may cause $`\mathrm{\Delta }_{\mathrm{RPA}}`$ to converge rather slowly with space dimension, i.e., a high number $`N_p`$ of HF particle orbitals may be needed (see below).

### IV.2 Numerical results for ground-state energies

We have analyzed the RPA ground state for the same dots whose excited states have been discussed in the preceding sections. Figure 6 displays the evolution of the ground state energy $`E_{\mathrm{RPA}}`$ with number of particle states $`N_p`$. The convergence with $`N_p`$ depends on each specific case: while the $`N=2`$, $`R_W=1`$ energy converges fast enough, the $`N=12`$ result is clearly still increasing for the maximum $`N_p`$ used in the numerics. This is not surprising, since one naturally expects that the larger $`N_h`$ the larger the space requirements for RPA convergence. A similar conclusion may be drawn regarding the intensity of the interaction, as given by the parameter $`R_W`$. In view of the obtained results, $`E_{\mathrm{RPA}}`$ must be extrapolated, in general, from the range of the computationally feasible values. We have empirically found that a 4-parameter function, $`ab/(1+cN_p)^{1/d}`$, leads to quite reasonable fits (see Fig. 6).

To assess the quality of $`E_{\mathrm{RPA}}`$ we compare next with available exact energies, i.e., with ab initio solutions of the many-body Schrödinger equation. In fact, the studied cases were chosen in part to allow for this comparison. For the $`N=2`$ dot the exact solution can be obtained by introducing center of mass and relative coordinates. Only the relative coordinate problem needs to be solved numerically and, because of the circular symmetry, it reduces to an straightforward 1d radial problem. We have solved this equation using a uniform radial grid and standard integration methods (see also Refs. YL2, ; Sim,  for recent discussions of the 2-electron dot exact solution as well as Ref. Yan02  for the application of the projection methods to this system). For the larger dots we shall compare with the Monte Carlo results of Pederiva et al. Pederiva .

Figure 7 presents a systematics for the $`N=2`$ dots with different values of the interaction-confinement ratio. As anticipated in Sec. IV-A, $`E_{\mathrm{RPA}}`$ lies below the exact value, a manifestation of an overestimated correlation energy in RPA. Nevertheless, for most of the results the RPA clearly improves the HF energy, since $`E_{\mathrm{RPA}}`$ is closer to the exact value than $`E_{\mathrm{HF}}`$. This is not the case, however, when $`R_W=1`$; the exact result lying almost halfway of RPA and HF in this instance. For $`N=6`$ and 12 the results in Fig. 6 increase towards the exact values. The extrapolation in these systems leads to $`\mathrm{\Delta }_{\mathrm{RPA}}(6)=1.27E_0`$ and $`\mathrm{\Delta }_{\mathrm{RPA}}(12)=3.21E_0`$; to be compared with the diffusion Monte Carlo values $`\mathrm{\Delta }_{\mathrm{MC}}(6)=1.13E_0`$ and $`\mathrm{\Delta }_{\mathrm{MC}}(12)=2.64E_0`$. The RPA overestimation of correlation energy is thus of $`12`$% and $`21`$% for 6 and 12, respectively. The larger difference for the 12-electron dot could be partly attributed to the extrapolation procedure being based on a relatively small space. The deviations with respect to the exact values for $`N=2`$ are not surprising since RPA and HF approaches are many-body theories aiming at an accurate description of large enough systems.

### IV.3 Occupation numbers and expectation value

In order to construct the RPA ground state we use Eq. (24) and the condition

$$L_z|0=0,$$
(27)
that insures the rotational invariance of the ground state. One seeks solutions of the form Rowe ; RS ; BR86

$$|0=N_0e^S|v$$
(28)
with the $`S`$ operator involving the creation of two bosons,

$$S=\frac{1}{2}\underset{minj}{}Z_{minj}b_{mi}^{}b_{nj}^{}.$$
(29)
In Eqs. (28) and (29) $`N_0`$ is a normalization constant and the matrix $`Z_{minj}`$ is, in general, complex and symmetric in the boson indexes, i.e., $`Z_{minj}=Z_{njmi}`$.

Using the general indentity

$$Fe^S=e^S\left(F+[F,S]+\frac{1}{2}[[F,S],S]+\mathrm{}\right)$$
(30)
for the operators $`FO_\lambda `$ and $`FL_z`$; as well as Eq. (24), Eq. (27) and the QBA one finds

$$\left(F+[F,S]\right)|v=0.$$
(31)
For the case of non-broken symmetry the above requirement leads to the following equations for the $`Z_{minj}`$ coefficients

$$Y_{mi}^{(\lambda )}=\underset{nj}{}Z_{minj}X_{nj}^{(\lambda )}.$$
(32)
When the mean-field breaks rotational symmetry we have to complement the above system of equations for vibrons, which are reduced by one equation, with the additional condition for the spurious mode

$$\mathrm{}_{mi}^{(z)}=\underset{nj}{}Z_{minj}\mathrm{}_{nj}^{(z)}.$$
(33)
Our interest is focused on the occupation numbers of the HF orbitals in the correlated ground state, $`0|a_m^{}a_m|0`$ and $`0|a_i^{}a_i|0`$. Using

$`[a_m^{}a_m,b_{nj}^{}]`$ $`=`$ $`\delta _{mn}b_{mj}^{},`$
$`[a_i^{}a_i,b_{nj}^{}]`$ $`=`$ $`\delta _{ij}b_{ni}^{}`$ (34)
one obtains

$`0|a_m^{}a_m|0`$ $`=`$ $`{\displaystyle \underset{inj}{}}Z_{minj}0|b_{mi}^{}b_{nj}^{}|0`$
$`0|a_i^{}a_i|0`$ $`=`$ $`1{\displaystyle \underset{mnj}{}}Z_{minj}0|b_{mi}^{}b_{nj}^{}|0.`$ (35)
By substituting in the above expressions the expansion of the boson operators $`b_{mi}^{}`$ in terms of vibrons and spurious modes, Eq. (20), we finally obtain

$`0|a_m^{}a_m|0`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i\lambda }{}}|Y_{mi}^\lambda |^2+0|\mathrm{\Phi }^2|0{\displaystyle \underset{i}{}}|\mathrm{}_{mi}^z|^2`$ (36)
$`0|a_i^{}a_i|0`$ $`=`$ $`1{\displaystyle \frac{1}{2}}{\displaystyle \underset{m\lambda }{}}|Y_{mi}^\lambda |^2`$ (37)
$``$ $`0|\mathrm{\Phi }^2|0{\displaystyle \underset{m}{}}|\mathrm{}_{mi}^z|^2`$
These equations generalize the result known in the literature Rowe  by introducing an additional term related to the canonical variables of the spurious mode $`\{L_z,\mathrm{\Phi }\}`$. It should be noted that the factor $`1/2`$ in Eqs. (36) and (37) is introduced following Ref. 14, , where the occupation numbers were calculated using Fermionic anticommutator rules, without referring to the QBA.

The above discussion can be easily extended to obtain the expectation value of any 1-body operator such as, e.g., the particle density

$$\widehat{\rho }(𝐫)=\underset{\alpha \beta }{}\rho _{\alpha \beta }(𝐫)a_\alpha ^{}a_\beta ,$$
(38)
where indexes $`\alpha `$ and $`\beta `$ run over all the HF set of orbitals and $`\rho _{\alpha \beta }(𝐫)=\phi _\alpha ^{}(𝐫)\phi _\beta (𝐫)`$ (with $`\phi _\alpha `$ the HF wave functions). Omitting for the sake of presentation the spatial dependence, one can write the RPA ground state density as

$$0|\widehat{\rho }|0=\underset{i}{}\rho _{ii}\mathrm{\Delta }\rho _h+\mathrm{\Delta }\rho _p$$
(39)
with

$`\mathrm{\Delta }\rho _h`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\lambda ijm}{}}\rho _{ij}Y_{mi}^{(\lambda )}Y_{mj}^{(\lambda )}`$ (40)
$`+`$ $`0|\mathrm{\Phi }^2|0{\displaystyle \underset{ijm}{}}\rho _{ij}\mathrm{}_{mj}^{(z)}\mathrm{}_{mi}^{(z)}`$
$`\mathrm{\Delta }\rho _p`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\lambda mnj}{}}\rho _{mn}Y_{mj}^{(\lambda )}Y_{nj}^{(\lambda )}`$ (41)
$`+`$ $`0|\mathrm{\Phi }^2|0{\displaystyle \underset{mnj}{}}\rho _{mn}\mathrm{}_{mj}^{(z)}\mathrm{}_{nj}^{(z)}`$
The occupation numbers Eqs. (36) and (37), and density expectation value Eq. (39) depend on the matrix element $`0|\mathrm{\Phi }^2|0`$. Involving the two-body operator $`\mathrm{\Phi }^2`$ this matrix element is not determined within the QBA MW69 . Therefore, we propose the following procedure to fix its value. Treating the unknown matrix element as a parameter ($`\gamma 0|\mathrm{\Phi }^2|0`$) we introduce the function

$$F(\gamma )=𝑑𝐫\left[\rho (\gamma ,𝐫)\stackrel{~}{\rho }(\gamma ,r)\right]^2$$
(42)
where $`\rho (\gamma ,𝐫)`$ refers to the RPA density Eq. (39) and

$$\stackrel{~}{\rho }(\gamma ,r)=\frac{1}{2\pi }_0^{2\pi }𝑑\theta \rho (\gamma ,𝐫)$$
(43)
is the corrresponding angular average. Since $`|0`$ is rotationally invariant, it is natural to require the physical value $`\gamma =\gamma _0`$ to be a minimum of $`F(\gamma )`$. Therefore, it fulfills the condition

$$\frac{dF(\gamma )}{d\gamma }|_{\gamma =\gamma _0}=0.$$
(44)
Resolving this equation we obtain the unknown quantity as

$$\gamma _0=\frac{A}{B},$$
(45)
where

$`A`$ $`=`$ $`{\displaystyle 𝑑𝐫\underset{k=0}{\overset{2}{}}\underset{m=1}{\overset{2}{}}(1)^{k+m}(a_k\stackrel{~}{a}_k)(b_m\stackrel{~}{b}_m)},`$
$`B`$ $`=`$ $`{\displaystyle 𝑑𝐫\left(\underset{m=1}{\overset{2}{}}(1)^m(b_m\stackrel{~}{b}_m)\right)^2},`$ (46)
and we have defined the $`a`$ and $`b`$ coefficients according to Eqs. (40) and (41) as

$`\mathrm{\Delta }\rho _h`$ $``$ $`a_1+\gamma b_1`$
$`\mathrm{\Delta }\rho _p`$ $``$ $`a_2+\gamma b_2`$
$`a_0`$ $``$ $`{\displaystyle \underset{i}{}}\rho _{ii}.`$ (47)
The tilde coefficients are defined through analogous equations with the circularly averaged densities.

### IV.4 Numerical expectation values

The above discussed Fig. 7 contains also the numerical values of the hole state occupation for $`N=2`$ –there is only one occuppied orbital (spin degenerate) in this case. The occupation is close to 1 for weak interaction (low $`R_W`$) and it decreases in an important way as the interaction strength is increased. This proofs that the RPA ground state deviates from the HF one when interaction becomes strong and correlations allow particles and holes to be found above and below the Fermi level, respectively. The occupations for the other two dots analyzed in this work are shown in Fig. 8. In these cases we have two (for $`N=6`$) and four (for $`N=12`$) different hole levels with occupations varying with their respective energies. In general, the RPA occupation numbers show a smoother distribution than the HF ones, the maximum variations occurring nearby the Fermi level.

The RPA density is compared with the exact one in Fig. 9, for the 2-electron dot. For all $`R_W`$ parameter values the RPA density is circularly symmetric; fulfilling the symmetry restoration discussed above. Therefore, in Fig. 9 we focus exclusively on the radial dependence. While the dot edge is well reproduced, a conspicuous feature is the underestimation of the central density by the RPA, specially at large $`R_W`$ values. Our result is in qualitative agreement with that of Reinhard Reinhard  for electrons in jellium spheres (where a circularly symmetric mean field was imposed). As shown by Reinhard the RPA correctly describes the low $`q`$ components (large $`r`$’s) of the form factor $`F(q)=𝑑r\mathrm{exp}(iqr)\rho (r)`$, but it fails for the large $`q`$ contributions (low $`r`$’s). The overcorrection of both central density and correlation energy seem thus to be peculiarities of the RPA.

Analogous results for the 6- and 12-electron dots are presented in Fig. 10. An excellent restoration of the circular symmetry is obtained for the $`N=6`$ dot –we remind the reader that the corresponding HF density was shown in Fig. 2. For $`N=12`$ the RPA density, though more circular than the HF one, still has some residual deformation. This can be surely attributed to incompleteness of the RPA space considered in the numerical calculation. In fact, for this dot the spurious mode separation (Fig. 4), the ground state energy (Fig. 6), and the density (Fig. 10) are all indicating that convergence with number of particle states is rather slow and, therefore, very difficult to be achieved in a numerical calculation. We emphasize that circular symmetry can be restored within the RPA only if the contribution from the spurious mode is taken into account.

## V Summary and outlook

The RPA is a general framework where mean field theory can be improved. Focussing on two-dimensional quantum dots we have reviewed the description of excited and ground states within the RPA. Some of the results known in the literature have been generalized to include broken symmetries at the mean field level. Quite importantly, the RPA provides a physical interpretation of this broken-symmetry phenomenon: mean-field solutions lacking one symmetry of the Hamiltonian represent, according to the RPA, true internal deformations of the system having an associated collective motion at vanishing energy (the spurious RPA mode). The ground state energy, occupation numbers, as well as general expectation values with 1-body operators, including the spurious-mode contribution, have been obtained. We proposed the procedure to restore the rotational symmetry broken at the mean field, which can be extended for other symmetry breaking cases.

Small quantum dots have provided us a good scenario to assess the above properties of the RPA by numerical calculations. For large enough values of the interaction-confinement ratio $`R_W`$ the HF mean field breaks circular symmetry; the electrons being localized in specific geometric distributions. In these cases we obtain an spurious RPA mode related to the global rotation, with its associated moment of inertia. The corresponding rotational spectra may be associated with the rotation of a Wigner molecule Koon . A systematics with $`R_W`$ for the 2-electron dot has been presented, while for 6- and 12-electron dots the calculations have been performed with an $`R_W`$ value suggested by the experiments.

When comparing with exact results, the RPA generally overestimates the correlation energy and the central dot densities by an amount which depends on the electron number and the Hamiltonian parameters. Another general behaviour of RPA is the slow convergence with space dimensions, requiring in practical applications a very high number of empty HF orbitals. For systems with many electrons this renders the numerical calculation quite a formidable task. In spite of these difficulties the RPA corrections essentially improve the HF energies with respect to exact values. It should be noted that the QBA only approximately satisfies the Pauli principle (see Eq.(9)). To solve this deficiency of the QBA could be one way to improve the description of the ground state. On the other hand, a more systematic treatment of the ground state energies would require, probably, alternative many body techniques.

The exact spin symmetry has been imposed on the single-particle states in this work, both at the HF and RPA levels. An interesting extension is, therefore, to relax this constraint allowing the HF orbitals to be general two-component spinors. The allowed symmetry breaking would eventually lead to ground states consisting of vectorial spin textures, lacking a single quantization axis. This type of states have been predicted for quantum dots in large magnetic fields and the question arises what collective spurious modes will be associated with them. This and related issues are left for a future work.

###### Acknowledgements.

This work was supported by Grant No. BFM2002-03241 from DGI (Spain). R. G. N. gratefully acknowledges support from the Ramón y Cajal programme (Spain).


text:
# Interference Commensurate Oscillations in Q1D Conductors

(February 24 2003, Submitted to Physical Review Letters)

## Abstract

We suggest an analytical theory to describe angular magnetic oscillations recently discovered in quasi-one-dimensional conductor (TMTSF)2PF6 \[see Phys. Rev. B 57, 7423 (1998)\] and define the positions of the oscillation minima. The origin of these oscillations is related to interference effects resulting from an interplay of quasi-periodic and periodic (”commensurate”) electron trajectories in an inclined magnetic field. We reproduce via calculations existing experimental data and predict some novel effects.

preprint: Lebed-Naughton

Quasi-one-dimensional (Q1D) organic conductors (TMTSF)2X (X = PF6, ClO4, etc.) demonstrate a variety of unique properties in a magnetic field in their superconducting \[1-3\], field-induced spin-density-wave \[1-3\], and metallic \[4-22,1-3\] phases. In the metallic phase \[4-22,1-3\], the Q1D electron spectra \[1-3\] of these compounds,

$$ϵ^\pm (𝐩)=\pm v_F(p_xp_F)+2t_b\mathrm{cos}(p_yb^{})+2t_c\mathrm{cos}(p_zc^{}),v_Fp_Ft_bt_c,$$
(1)
are characterized by two open sheets of Fermi surface (FS). Therefore, traditional magnetic oscillations (related to Landau level quantization)  cannot exist in these materials. \[Here $`+()`$ stands for the right (left) sheet of the Fermi surface, $`v_F`$ and $`p_F`$ are the Fermi velocity and Fermi momentum, respectively; $`t_b`$ and $`t_c`$ are the overlapping integrals between electron wave functions; $`\mathrm{}1`$\].

Surprisingly, the metallic phases of (TMTSF)2X materials exhibit a number of unconventional magnetic oscillations related to an open Q1D FS (1). Among them are ”magic angles” (MA) \[4-9,1-3\], the first angular oscillations with a clear Fermi-liquid (FL) physical meaning - Danner-Kang-Chaikin (DKC) oscillations , the ”third angular effect” (TAE) \[11-19\], and rich angular oscillations recently discovered by Lee and Naughton  in (TMTSF)2PF6 and by Yoshino et al.  in (DMET)2I3. We call the latter ”interference commensurate” (IC) oscillations which reflects their physical meaning revealed in this Letter. Note that, despite all these observations of different magnetic oscillations \[4-19\], a question on the existence or not of FL behavior in Q1D metallic phases of (TMTSF)2X conductors is still controversial \[3,9,20-22\].

Numerical solutions  of the Boltzmann kinetic equation for Q1D metal (1) in an inclined magnetic field,

$$𝐇=H(\mathrm{cos}\theta \mathrm{cos}\varphi ,\mathrm{cos}\theta \mathrm{sin}\varphi ,\mathrm{sin}\theta ),$$
(2)
have given a hint on the FL nature of the IC oscillations. Nevertheless, due to very complex behavior of these oscillations, their physical meaning has not been revealed and their properties have not been described in detail. In particular, a hypothesis  that minima of resistivity $`\rho _{zz}(𝐇)`$ correspond to some ”commensurate” directions of a magnetic field,

$$\mathrm{sin}\varphi =N(\frac{b^{}}{c^{}})\mathrm{tan}\theta ,$$
(3)
(where N is an integer) has not been theoretically proven.

The goals of our Letter are: 1) to derive an analytical expression for $`\rho _{zz}(𝐇)`$ and to define the positions of its minima; 2) to reveal a quantum interference nature of the IC oscillations, 3) to compare our results with experiment , and 4) to predict some novel qualitative effects. In particular, we demonstrate that the origin of the IC oscillations is related to special ”commensurate” electron trajectories in a magnetic field (different from the MA trajectories ), where an average electron velocity along the $`𝐳(𝐜^{})`$-axis is non-zero, $`<v_z(t)>_t0`$ (see Fig.1). Note that an importance of this condition for low-dimensional conductors was pointed out by DKC , Osada et al. , and by Kartsovnik and Yakovenko et al.  in different context.

Below, we consider the most general electron trajectories, characterized by two angles, $`\theta `$ and $`\varphi `$ \[see Eq. (2)\]. We show that they correspond to two novel types of angular magnetic oscillations: ”commensurate” oscillations and some ”generalized DKC” oscillations. In particular, we demonstrate that $`\rho _{zz}(𝐇)`$ in this case is defined by quantum interference effects related to periodic and quasi-periodic electron motion along the open Q1D FS (1). For small values of $`\theta `$ \[see Eq.(2)\], these effects occur between some narrow areas on the Q1D FS (1) \[i.e., ”effective stripes” (ES) parallel to $`𝐩_z`$-axis \], with their positions being dependent on the magnetic field orientation (see Fig.1c). These unique features demonstrate an unusual physical meaning of the IC oscillations. We also show that the ”commensurate” directions of a magnetic field (3) indeed correspond to minima in $`\rho _{zz}(𝐇)`$ at large enough $`\theta `$, whereas at smaller $`\theta `$ ($`|\theta |5^o`$ in our case) $`\rho _{zz}(𝐇)`$ minima occur for only even or only odd values of N in Eq. (3), depending on the particular value of $`\theta `$.

Let us discuss the physical meaning of the IC oscillations by analyzing quasi-classical electron trajectories in an inclined magnetic field (2). As usual, the trajectories are solutions of the equations of motion :

$$d𝐩/dt=(e/c)[𝐯(𝐩)\times 𝐇],𝐯(𝐩)=dϵ(𝐩)/d𝐩.$$
(4)
For Q1D electron spectrum (1), Eqs.(4) can be rewritten as follows :

$$d(p_yb^{})/dt=\omega _b(\theta ),d(p_zc^{})/dt=\omega _c(\theta ,\varphi )\omega _c^{}(\theta ,\varphi )\mathrm{sin}(p_yb^{}),$$
(5)
where

$`\omega _b(\theta )=ev_FHb^{}\mathrm{sin}\theta /c,\omega _c(\theta ,\varphi )=ev_FHc^{}\mathrm{cos}\theta \mathrm{sin}\varphi /c,`$
$`\omega _c^{}(\theta ,\varphi )=(v_y^o/v_F)(ev_FHc^{}/c)\mathrm{cos}\theta \mathrm{cos}\varphi ,v_y^o=2t_bb^{}.`$ (6)
From Eqs.(5), one defines electron trajectories in a reciprocal plane $`(p_y,p_z)`$:

$$p_zc^{}=p_z^oc^{}+\frac{\omega _c(\theta ,\varphi )}{\omega _b(\theta )}(p_yp_y^o)b^{}+\frac{\omega _c^{}(\theta ,\varphi )}{\omega _b(\theta )}[\mathrm{cos}(p_yb^{})\mathrm{cos}(p_y^ob^{})],$$
(7)
with a velocity component along the $`𝐳`$-axis being

$`v_z(p_y)=2t_cc^{}\mathrm{sin}(`$ $`p_z^oc^{}+{\displaystyle \frac{\omega _c(\theta ,\varphi )}{\omega _b(\theta )}}(p_yp_y^o)b^{}`$ (8)
$`+{\displaystyle \frac{\omega _c^{}(\theta ,\varphi )}{\omega _b(\theta )}}[\mathrm{cos}(p_yb^{})\mathrm{cos}(p_y^ob^{})]),p_yt.`$
Note that electron trajectories in a reciprocal plane $`(p_y,p_z)`$ become periodic  for ”commensurate” directions of a magnetic field (3) (see Figs.1a,b). Therefore, an average velocity over electron path,

$$<v_z(t)>_t\mathrm{sin}[p_z^oc^{}\frac{\omega _c(\theta ,\varphi )}{\omega _b(\theta )}p_y^ob^{}\frac{\omega _c^{}(\theta ,\varphi )}{\omega _b(\theta )}\mathrm{cos}(p_y^ob^{})]J_N[\frac{\omega _c^{}(\theta ,\varphi )}{\omega _b(\theta )}],$$
(9)
(which is zero for all ”non-commensurate” trajectories), for ”commensurate” orbits (3,7), becomes zero only at zero values of the N-order Bessel functions:

$$J_N[\frac{\omega _c^{}(\theta ,\varphi )}{\omega _b(\theta )}]J_N[\frac{2t_bc^{}\mathrm{cos}\varphi }{v_F\mathrm{tan}\theta }]=0.$$
(10)
From general metals theory , it is known that $`<v_z(t)>_t0`$ results in saturating behavior of the conductivity $`\sigma _{zz}(\theta ,\varphi )`$ for high magnetic fields, whereas, at $`<v_z(t)>_t=0`$, $`\sigma _{zz}(\theta ,\varphi )`$ decreases with increasing field. Therefore, there appear maxima in $`\sigma _{zz}(\theta ,\varphi )`$ \[i.e., minima in $`\rho _{zz}(\theta ,\varphi )`$\] at ”commensurate” angles (3) (see Fig.2) which are a novel type of angular magnetic oscillations. In a similar way, Bessel functions zeros (10) lead to the appearance of maxima in $`\rho _{zz}(𝐇)`$ (see Fig.2) which are a generalization of the DKC oscillations corresponding to $`N=0`$ in Eq. (10). In other words, the IC oscillations are characterized by minima in $`\rho _{zz}(𝐇)`$ at ”commensurate” angles (3) which are modulated by oscillatory Bessel functions (10).

Let us discuss electron motion along the $`𝐳`$-axis in real space, $`z(t)=_0^tv_z(t^{^{}})𝑑t^{^{}}`$. Here, contributions to $`z(t)`$ from different Brillouin zones \[see Eq.(8)\] are in-phase for ”commensurate” trajectories (3) (see Fig.1c). These interference effects correspond to maxima of conductivity $`\sigma _{zz}(\theta ,\varphi )`$ \[i.e., minima of resistivity $`\rho _{zz}(\theta ,\varphi )`$\] (see Fig. 2). For $`\omega _c^{}(\theta ,\varphi )/\omega _b(\theta )1`$ (i.e., for small $`\varphi `$), the integral $`z(t)=_0^tv_z(t^{^{}})𝑑t^{^{}}`$ can be evaluated by means of the stationary-phase method and is determined by interference effects between some ES located near the points $`p_yb^{}=\mathrm{arcsin}[\omega _c(\theta ,\varphi )/\omega _c^{}(\theta ,\varphi )]+2\pi n`$ and $`p_yb^{}=\pi \mathrm{arcsin}[\omega _c^{}(\theta ,\varphi )/\omega _c(\theta ,\varphi )]+2\pi n`$ on Q1D FS (1) (see Fig.1c), where $`n`$ is an integer.

To develop an analytical theory of the IC oscillations and to further demonstrate their quantum interference nature, we make use of the quasi-classical approximation for electron motion , where a magnetic field is introduced by the Peierls substitution , $`𝐩𝐏(e/c)𝐀`$. We choose a vector-potential of magnetic field (2), in the following form

$$𝐀=(0,x\mathrm{sin}\theta ,x\mathrm{cos}\theta \mathrm{sin}\varphi +y\mathrm{cos}\theta \mathrm{cos}\varphi )H.$$
(11)
In this case, electron wave functions in a mixed $`(p_y,x)`$-representation , $`\mathrm{\Psi }_ϵ^\pm (p_y,x)=\mathrm{exp}(\pm ip_Fx)\psi _ϵ^\pm (p_y,x)`$, are solutions of the Schrodinger equations ,

$$(iv_F\frac{d}{dx}+2t_b\mathrm{cos}[p_yb^{}\frac{\omega _b(\theta )x}{v_F}])\psi _ϵ^\pm (p_y,x)=ϵ\psi _ϵ^\pm (p_y,x),$$
(12)
and can be expressed as

$$\psi _ϵ^\pm (p_y,x)=\mathrm{exp}(\pm i\frac{ϵ}{v_F}x)\mathrm{exp}[\pm \frac{2it_b}{\omega _b(\theta )}(\mathrm{sin}[p_yb^{}\frac{\omega _b(\theta )x}{v_F}]\mathrm{sin}[p_yb^{}])].$$
(13)
Note that the $`𝐳`$-component of the quasi-classical velocity operator, $`\widehat{𝐯}(𝐩)=d\widehat{ϵ}(𝐩)/d𝐩`$ , in a gauge (11) is equal to

$$\widehat{v}_z(p_z,x,y)=2t_cc^{}\mathrm{sin}[p_zc^{}\frac{\omega _c(\theta ,\varphi )x}{v_F}\frac{\omega _c^{}(\theta ,\varphi )y}{v_y^0}],y=i\frac{d}{p_y}.$$
(14)
It is possible to show that wave functions (13) are eigenfunctions of velocity operator (14) with their eigenvalues being:

$`\widehat{v}_z(p_z,x,y)\psi _ϵ^\pm (p_y,x)=2t_cc^{}\mathrm{sin}[`$ $`p_zc^{}{\displaystyle \frac{\omega _c(\theta ,\varphi )x}{v_F}}\pm {\displaystyle \frac{\omega _c^{}(\theta ,\varphi )}{\omega _b(\theta )}}`$
$`\times `$ $`(\mathrm{cos}[p_yb^{}{\displaystyle \frac{\omega _b(\theta )x}{v_F}}]\mathrm{cos}[p_yb^{}])]\psi _ϵ^\pm (p_y,x).`$ (15)
Since wave functions (13) and their matrix elements (15) are known, $`\sigma _{zz}(H,\theta ,\varphi )`$ can be evaluated by means of the Kubo formalism. After straightforward but rather routine calculations, we obtain:

$`\sigma _{zz}(H,\theta ,\varphi )`$ $`{\displaystyle _0^{2\pi }}d(p_yb^{}){\displaystyle _0^{\mathrm{}}}dx\mathrm{exp}({\displaystyle \frac{x}{v_F\tau }})`$ (16)
$`\times \mathrm{cos}({\displaystyle \frac{\omega _c(\theta ,\varphi )x}{v_F}}+{\displaystyle \frac{\omega _c^{}(\theta ,\varphi )}{\omega _b(\theta )}}[\mathrm{cos}(p_yb^{}{\displaystyle \frac{\omega _b(\theta )x}{v_F}})\mathrm{cos}(p_yb^{})]),`$
where $`\tau `$ is an electron relaxation time. Double integral (16) can be rewritten in a simple form:

$$\sigma _{zz}(H,\theta ,\varphi )=\sigma _{zz}(0)\underset{N=\mathrm{}}{\overset{+\mathrm{}}{}}(\frac{J_N^2[\omega _c^{}(\theta ,\varphi )/\omega _b(\theta )]}{1+\tau ^2[\omega _c(\theta ,\varphi )N\omega _b(\theta )]^2}),$$
(17)
where for a Q1D conductor (1):

$$\rho _{zz}(H,\theta ,\varphi )1/\sigma _{zz}(H,\theta ,\varphi ).$$
(18)
To summarize, Eqs.(17) and (18) provide analytical expressions for experimentally measured resistivity $`\rho _{zz}(H,\theta ,\varphi )`$ . As is seen from Eq. (17), $`\sigma _{zz}(H,\theta ,\varphi )`$ possesses maxima \[i.e., $`\rho _{zz}(H,\theta ,\varphi )`$) possesses minima\] at $`\omega _c(\theta ,\varphi )=N\omega _b(\theta )`$ if $`\omega _c(\theta ,\varphi ),\omega _b(\theta )1/\tau `$. This coincides with the ”commensurability” condition (3). It is important to note that this theory, based on Eqs.(17),(18), predicts no angular oscillations  at MA directions of a magnetic field  (i.e., at $`\varphi =\pi /2`$) since $`\omega _c^{}=0`$ in Eq.(17). Therefore, the previous interpretation \[15-19\] of the IC oscillations as a simple combination of the MA effects and the DKC oscillations is too oversimplified. We also stress that the integration of the quasi-periodic function in Eq.(16) corresponds to interference effects in matrix elements of velocity operator (15) and, thus, directly demonstrates the interference nature of the IC oscillations. Indeed, two different periods in Eq.(16) become commensurate at $`\omega _c(\theta ,\varphi )=N\omega _b(\theta ,\varphi )`$ \[i.e., for ”commensurate” trajectories (3)\] and, therefore, the integral (16) is increased due to these interference effects.

Let us consider Eq.(17) at small enough angles, $`\varphi \pi /2`$ and $`\theta 2t_bc^{}/v_F`$. In this case, one can make use of an asymptotic expression for the Bessel functions in Eqs.(10),(17): $`J_N(2t_bc^{}/v_F\mathrm{tan}\theta )\mathrm{cos}(2t_bc^{}/v_F\mathrm{tan}\theta \pi /4\pi N/2)`$. Therefore, depending on the value of parameter $`2t_bc^{}/v_F\mathrm{tan}\theta `$, the Bessel functions of even orders are bigger than those of odd orders or vise versa. At $`\omega _c\tau ,\omega _b\tau 1`$, this results in the appearance of minima of $`\rho _{zz}(\varphi ,\theta )`$ in Eqs.(17),(18) only for all even or only for all odd values of the integer $`N`$. We call this phenomenon ”even-odd” angular resonance (see Fig.2b). It is important that Eq.(17) demonstrates also another kind of angular oscillations (i.e., ”extended DKC” oscillations) related to zeros of the Bessel functions (10). By analyzing Eqs.(16)-(18), it is possible to show that $`\rho _{zz}(H,\theta ,\varphi )`$ is characterized by an unusual linear behavior for ”non-commensurate” directions of a magnetic field and small $`\theta 2t_bb^{}/v_F`$,

$$\rho _{zz}(H,\theta ,\varphi )|H|,$$
(19)
whereas, for commensurate directions (3), $`\rho _{zz}(H,\theta ,\varphi )`$ saturates with increasing magnetic field. The latter effect is different from the linear magnetoresistance predicted in Ref. .

In Fig.2, we compare the experimental data  with Eqs. (17),(18) using the same values of parameters, $`t_a/t_b=8.5`$ and $`\omega _c(\theta =0,\varphi =\pi /2,H=1T)\tau =15`$, for all three theoretical curves. These curves not only demonstrate qualitative but quantitative agreement between theory (17),(18) and experiment , in a broad region of magnetic field orientations, $`\varphi 20^o`$ . Note that, for $`\theta =3^o`$, $`\rho _{zz}(H,\theta ,\varphi )`$ minima appear both theoretically and experimentally only for odd integers $`N`$ in Eq.(3) \[see Fig.(2)\] which is an agreement with the ”odd” angular resonance effect discussed above. The agreement is far better than what has been achieved with numerical solutions of the kinetic equation \[15-19\], which do not account in full for the interference nature of the IC oscillations. We have initiated more detailed experiments to try to confirm other effects predicted in the Letter such as ”even” angular resonance, ”generalized DKC” oscillations (10), and linear magnetoresistance (19).

Let us discuss peculiarities of the interference effects at small angles $`\theta 2t_bb^{}/v_F`$ in a clean limit where $`\omega _c(\theta ,\varphi )\tau ,\omega _b(\theta )\tau \mathrm{}`$. In this case, the integral (16) can be evaluated by means of the stationary-phase method and is determined by the close proximity of the following two series of ”effective points” (EP) in $`(p_y,x)`$-space:

$`\omega _b(\theta )z_n/v_F=\mathrm{arcsin}[\omega _c(\theta ,\varphi )/\omega _c^{}(\theta ,\varphi )]+2\pi n,z=xv_Fp_yb^{}/\omega _b(\theta ),`$ (20)
$`\omega _b(\theta )z_n/v_F=\pi \mathrm{arcsin}[\omega _c(\theta ,\varphi )/\omega _c^{}(\theta ,\varphi )]+2\pi n,z=xv_Fp_yb^{}/\omega _b(\theta ),`$ (21)
where $`n`$ is an integer. In quantum mechanical language, these EP (20),(21) correspond to ES on the Q1D FS (1) discussed above. The contributions to (16) from the EP (20),(21) are characterized by non-zero phase shifts with their values being dependent on magnetic field orientation (2). Note that these phase shifts are integer values of $`2\pi `$ for different $`n`$ in Eqs.(20),(21) only at $`\omega _c(\theta ,\varphi )=N\omega _b(\theta ,\varphi )`$. Therefore, at small angles $`\theta `$, the IC oscillations (3),(10) can be interpreted in terms of interference effects between these two infinite series of electron waves (20),(21). It is easy to show that, in our picture, different integers $`n`$ in Eqs.(20),(21) correspond to different Brillouin zones for a quasi-classical electron motion picture in the extended Brillouin zone.

This work was supported in part by National Science Foundation, grant number DMR-0076331, the Department of Energy, grant number DE-FG02-02ER63404, and by the INTAS grants numbers 2001-2212 and 2001-0791. One of us (AGL) is thankful to N.N. Bagmet and E.V. Brusse for useful discussions.


text:
# 1 Introduction

## 1 Introduction

$`1+1`$ dimensional integrable models are theoretical laboratories where one can study selected problems in field theory. These include the question whether the continuum limit of the lattice regularized version of a field theory coincides with the theory constructed by continuum techniques. One can study here structural questions related to asymptotic freedom, topological excitations etc., and, in the case that the theories agree, the nature of lattice artifacts.

The form factor bootstrap method  of constructive field theory has been used to construct many integrable models. The method consists of three main steps. The starting point is the exact scattering matrix that is calculated (sometimes conjectured) by the bootstrap method . Next the form factors (matrix elements of local operators between scattering states) are calculated by solving the form factor (FF) axioms . Finally correlation functions are constructed by inserting a complete set of intermediate states between the local operators.

The program has so far only been carried out completely for the Ising model . The S–matrix is known for many integrable models and if there are no internal degrees of freedom also the many–particle form factors are easily constructed as a multiple product of the basic two–particle form factor (which is known again for most models). On the other hand, for models with internal degrees of freedom the solution of the FF axioms, matrix functional difference equations in this case, is not known in general. Nevertheless, the form factors of some models belonging to this class, the Sine–Gordon model, the chiral Gross–Neveu model, and the O(3) and O(4) nonlinear $`\sigma `$–models, are available . Finally (with the exception of the Ising model) the $`r`$–particle contributions to the correlation functions can only be computed numerically and therefore this sum has to be truncated after the first few terms.

Among models with internal (isospin) degrees of freedom, it is the O(3) model, where one can go furthest. This is because in this special case the many–particle form factors (after removing a simple factor) are polynomials in the particle rapidities which can be determined recursively . Fortunately it is also a model which is of particular interest since it has many properties akin to Yang–Mills theory in 4 dimensions in that it exhibits asymptotic freedom and has instanton solutions.

The question mentioned above concerning the equivalence between the FF construction and that from the lattice regularization has been addressed in ref. . The two constructions are both non–perturbative and in both cases it is the low energy properties that are most easily accessible. However in the bootstrap approach it is simpler to obtain reliable results at higher energies. For example for the Fourier transform of the 2–point correlation function of Noether currents, the contributions of higher $`r`$–particle states become significant at typical energies rapidly increasing with $`r`$. Of course, despite the fact that in practice renormalized perturbation theoretic computations often seem to be quite accurate down to unexpectedly low energies, the property of asymptotic freedom is a statement pertaining to asymptotically high energies, and to establish this property in the bootstrap framework it is necessary to have rigorous knowledge on the contributions of $`r`$–particle states with arbitrarily large $`r`$. Although the latter has not yet been achieved, based on the exact form factors up to 6 particles, Max Niedermaier and one of the present authors (J.B.) presented convincing evidence for remarkable scaling properties of the higher intermediate particle contributions , which are very powerful because they make it possible to include dominant contributions from all states in the form factor expansion of correlation functions. Furthermore the scaling hypothesis is completely compatible with asymptotic freedom.

Although the lattice regularization is an important non–perturbative method (and practically the only one available in 4d QCD), there are many physical phenomena which are inaccessible in this approach for example nuclear structure functions at small Bjorken $`x`$. In the framework of the bootstrap approach in 2d however, we hope that properties of the structure functions at small $`x`$ can be extracted, an insight which would hopefully also be relevant for QCD. This hope would be realized if we could find evidence for scaling properties of the structure functions similar to those already observed for the 2–point function spectral densities.

With this goal we have initiated a project to compute the structure functions in the O(3) sigma model, and with the known form factors we have computed contributions from intermediate states up to 5 particles. To establish their scaling properties however it is helpful to extend the list of known form factors to higher particle number. In this paper we discuss the calculation of the 7–particle form factors, which although in principle trivial, is technically challenging because for $`r=7`$ one has to deal with quite large polynomials. We will report on the structure functions in a future publication. Here we restrict attention to the contribution of the 7–particle intermediate states to the spin–spin correlation function, which yields an additional test of the scaling hypothesis. We also exhibit simple bounds on the square of the known form factor polynomials, which if they could be generalized to an arbitrary number of particles might be useful for analytic estimates to establish general properties.

## 2 Some basic definitions

The central object of attention in this paper is the 7–particle form factor of the spin operator of the O(3) model. This is the $`r=7`$ case of

$$0|S^a(0)|a_1,\theta _1;\mathrm{};a_r,\theta _r=\frac{2}{\sqrt{\pi }}f_{a_1\mathrm{}a_r}^a(\theta _1,\mathrm{},\theta _r),$$
(1)
where the particles are labelled by their isospin indices and rapidities and the spin operator $`S^a`$ also carries an isospin index. The normalization factor $`2/\sqrt{\pi }`$ is introduced here for later convenience. The spin operator itself is normalized by

$$0|S^a(0)|a_1,\theta _1=\delta ^{aa_1},b,\chi |a,\theta =4\pi \delta ^{ab}\delta (\chi \theta ),$$
(2)
which means that the Minkowski propagator has unit residue at $`p^2=M^2`$, where $`M`$ is the physical mass of the O(3) particle.

Introducing the squared form factors

$$\delta ^{ab}F^{(r)}(u)=\underset{a_1\mathrm{}a_r}{}f_{a_1\mathrm{}a_r}^a(\theta _1,\mathrm{},\theta _r)^{}f_{a_1\mathrm{}a_r}^b(\theta _1,\mathrm{},\theta _r)$$
(3)
the spectral density is given by

$$\rho ^{\mathrm{spin}}(\mu )=\frac{4}{\pi }\underset{k=0}{\overset{\mathrm{}}{}}\rho ^{(2k+1)}(\mu ),$$
(4)
where

$$\rho ^{(r)}(\mu )=_0^{\mathrm{}}\frac{\mathrm{d}u_1\mathrm{}\mathrm{d}u_{r1}}{(4\pi )^{r1}}F^{(r)}(u)\delta (\mu M^{(r)}(u)).$$
(5)
Here $`u_j=\theta _j\theta _{j+1}`$ are rapidity differences and $`M^{(r)}(u)`$ is the $`r`$–particle invariant mass. The Fourier transform of the correlation function is represented as the Stieltjes transform of the spectral density:

$$I^{\mathrm{spin}}(p^2)=_0^{\mathrm{}}d\mu \frac{\rho ^{\mathrm{spin}}(\mu )}{p^2+\mu ^2}.$$
(6)
We will compare (the truncation of) (6) to the results of perturbation theory. The 2–loop order perturbative result is

$$p^2I^{\mathrm{spin}}(p^2)=\lambda _1\left\{\frac{1}{\alpha (p)}+(2+\xi _0)+(2+\xi _0)\alpha (p)+O\left(\alpha ^2(p)\right)\right\}.$$
(7)
Here the running coupling function $`\alpha (p)`$ is the solution of

$$\frac{1}{\alpha (p)}+\mathrm{ln}\alpha (p)=\mathrm{ln}\frac{p}{M}$$
(8)
and the parameter $`\xi _0`$ gives the connection between the perturbative mass parameter $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}`$ and the exact mass gap $`M`$. In the O(3) model their ratio is known exactly :

$$\xi _0=\mathrm{ln}\frac{M}{\mathrm{\Lambda }_{\overline{\mathrm{MS}}}}=\mathrm{ln}81\mathrm{\hspace{0.17em}\hspace{0.17em}1.07944}.$$
(9)
The overall constant $`\lambda _1`$ cannot be calculated in perturbation theory, but has been determined exactly  using the scaling hypothesis (see also Sect. 5):

$$\lambda _1=\frac{4}{3\pi ^2}.$$
(10)
Instead of the physical form factor (1) it is convenient to consider the reduced form factors $`g_{a_1\mathrm{}a_r}^a`$ defined by

$$f_{a_1\mathrm{}a_r}^a(\theta _1,\mathrm{},\theta _r)=\frac{\pi ^{\frac{3r}{2}1}}{2}\mathrm{\Psi }(\theta _1,\mathrm{},\theta _r)g_{a_1\mathrm{}a_r}^a(\theta _1,\mathrm{},\theta _r).$$
(11)
Here

$$\mathrm{\Psi }(\theta _1,\mathrm{},\theta _r)=\underset{i<j}{}\psi (\theta _i\theta _j),$$
(12)
with

$$\psi (\theta )=\frac{\theta i\pi }{\theta (2\pi i\theta )}\mathrm{tanh}^2\frac{\theta }{2}.$$
(13)
We note that (13) is the basic 2–particle form factor and the product (12) would be the complete $`r`$–particle form factor if there were no isospin degrees of freedom.

The advantage of using the reduced form factors is that they are polynomials and can (in principle) be recursively calculated . The simplest nontrivial case is $`r=3`$:

$$g_{a_1a_2a_3}^a(\theta )=\delta ^{aa_3}\delta ^{a_1a_2}(\theta _2\theta _1)+\delta ^{aa_2}\delta ^{a_1a_3}(\theta _1\theta _32i\pi )+\delta ^{aa_1}\delta ^{a_2a_3}(\theta _3\theta _2).$$
(14)
The higher reduced form factors very rapidly become complicated because of the large number of isospin components and rapidity variables.

The quantity entering the spectral densities and two–point functions is the absolute square of the form factors, summed over the internal symmetry indices. For the reduced form factors the corresponding quantities are

$$G^{(r)}(\theta _1,\mathrm{},\theta _r)=\frac{1}{3}\underset{aa_1\mathrm{}a_r}{}|g_{a_1\mathrm{}a_r}^a(\theta _1,\mathrm{},\theta _r)|^2.$$
(15)
For the $`r=3`$ example we get

$$G^{(3)}(\theta _1,\theta _2,\theta _3)=2\left[(\theta _1\theta _2)^2+(\theta _1\theta _3)^2+(\theta _2\theta _3)^2\right]+12\pi ^2.$$
(16)
## 3 The 7–particle form factors

The 7–particle reduced form factors

$$g_{a_1\mathrm{}a_7}^a(\theta _1,\mathrm{},\theta _7)$$
(17)
are polynomials of degree 15 in the rapidity differences. To count the independent isospin components we can obviously fix the operator index $`a=3`$ since the rest of the components are related by simple isospin transformations. Then we get 274 non–vanishing components not counting terms related by the $`12`$ isospin symmetry twice. These components can be calculated from the known 5–particle form factors using the inhomogeneous FF axioms . The resulting polynomials contain typically $`10^4`$ terms. Actually it is not necessary to do this calculation for all the 274 components since we can find 4 components such that the rest can be obtained from them by using the homogeneous FF axioms (corresponding to permutations of the rapidity variables). After having computed all terms we checked once more that all the FF axioms are satisfied by (17) and also that the coefficients of the leading terms in the last variable (those proportional to $`\theta _7^5`$) are proportional to the 6–particle reduced form factors as they should .

After the calculation of the form factors the next step is to compute the sum of the absolute square of the components since this is needed in the calculation of the 7–particle contribution to the spin–spin correlation function and other physical quantities. This calculation is rather challenging because at intermediate steps one has to handle quite large polynomials (of degree 30 in 7 variables consisting of about $`410^5`$ terms). The final result can be simplified enormously because, being a polynomial symmetric under permutations of the rapidity variables, it can be expanded in terms of the basic symmetric polynomials $`\sigma _1,\sigma _2,\mathrm{},\sigma _7`$ (defined in Eq. (A.3) of ref. ). Written in this way, the final result consists of only 3214 terms and is reduced to a manageable size. 111Which is however still too large to usefully reproduce here.

To be sure that we have obtained the correct result for the square $`G^{(7)}(\theta _1,\mathrm{},\theta _7)`$ we have performed a number of checks. First we checked that the leading terms in the last variable (the coefficient of $`\theta _7^{10}`$) are twice the analogous 6–particle square $`G^{(6)}(\theta _1,\mathrm{},\theta _6)`$ . Next we checked that the overall leading terms (the terms of total degree 30) agree with those of a simple symmetric polynomial $`P^{(7)}(t)`$ . This is defined by (for $`r`$ particles)

$$P^{(r)}(t)=\underset{\mathrm{perms}}{}P_0^{(r)}(t),$$
(18)
where

$$P_0^{(r)}(t)=\underset{ji>1}{}\left[(\theta _i\theta _j)^2+t\pi ^2\right].$$
(19)
In (18) the sum extends over the $`r!`$ permutations of the rapidity variables $`\theta _1,\mathrm{},\theta _r`$ and in (19) the correction terms depending on the parameter $`t`$ do not affect the overall leading terms, they are included here for later convenience. Finally, we verified that $`G^{(7)}(\theta _1,\mathrm{},\theta _7)`$ vanishes at the points

$$\theta _1\theta _2=\theta _2\theta _3=i\pi ,\theta _j=\mathrm{arbitrary}j>3.$$
(20)
The fact that $`G^{(7)}`$ passed all these nontrivial tests gives us some confidence in the correctness of our results.

## 4 The spin–spin correlation function

Having calculated the square of the 7–particle form factors the corresponding contributions to the spin–spin correlation function and the spectral density can now be computed straightforwardly. The only difficulty is that to get correct numerical results one has to use high precision arithmetic in order to avoid huge rounding errors in calculating $`G^{(7)}`$. As explained above, to put it into a manageable form, we expressed it in terms of the basic symmetric polynomials $`\sigma _1,\mathrm{},\sigma _7`$. Although it is a sum of absolute squares and thus it is obviously positive, it is not manifestly positive in this form. Actually we experienced large cancellations between positive and negative contributions. After rescaling all rapidities by $`\pi `$ all coefficients of the polynomial become integers and we can illustrate the large degree of cancellation by considering the following integer rapidities:

$$\theta _1=12,\theta _2=11,\theta _3=10,\theta _4=9,\theta _5=8,\theta _6=7,\theta _7=0.$$
(21)
We can exactly calculate the sum of positive and negative contributions in this case and we get

$`\mathrm{positive}\mathrm{part}`$ $`:`$ $`26512984286926120356867283491794696305689824`$
$`\mathrm{negative}\mathrm{part}`$ $`:`$ $`26512984286926120356855853186993819121689824`$
$`\mathrm{total}`$ $`:`$ $`11430304800877184000000`$
showing a 21–digit cancellation here!

We have computed the 7–particle contribution to the spin–spin correlation function using the VEGAS integration routine, and a subroutine computing $`G^{(7)}`$ invoking quartic precision (32 digit) arithmetic. The results, compared with the prediction of two–loop perturbation theory, are shown in Figure 1. Note that here we have no free parameters at our disposal in the perturbative calculation as would be the case in most other models. As discussed above, we know the exact relation between the perturbative $`\mathrm{\Lambda }`$ parameter and the particle mass $`M`$ and also the absolute normalization of the perturbative curve is available in this model. The form factor results are in very good agreement with perturbation theory in the expected (high) energy range. The small deviation for energies $`p/M10^4`$ can be accounted for by the contribution of $`r>7`$ intermediate particles. To illustrate the fact that this good agreement is quite nontrivial, in Figure 1 we also show the perturbative result with the overall factor changed arbitrarily by 5%.

## 5 Scaling

The 7–particle results also corroborate the scaling hypothesis  for the spectral densities. To study this aspect we introduce the modified $`r`$–particle spectral density depending on the logarithmic variable $`x`$ by the definition

$$\mu \rho ^{(r)}(\mu )=R^{(r)}(x),\mu =Me^x.$$
(22)
These are defined for all $`r2`$, where the cases with $`r`$ even refer to the spectral densities of the two–point function of the isospin Noether current. The graph of this function is a bell–shaped curve starting as zero at $`x=\mathrm{ln}r`$, reaching its maximum $`M^{(r)}`$ at $`x=\xi ^{(r)}`$ and then slowly decreasing for larger $`x`$. Let us introduce the rescaled spectral density $`Y^{(r)}`$ by

$$Y^{(r)}(z)=\frac{1}{M^{(r)}}R^{(r)}\left(\xi ^{(r)}z\right).$$
(23)
It has been found (based on the study of up to 6 particles) that the shape of the rescaled spectral density and the parameters $`\xi ^{(r)},M^{(r)}`$ satisfy self–similarity

$$\underset{r\mathrm{}}{lim}Y^{(r)}(z)Y(z)$$
(24)
with universal shape function $`Y(z)`$ and asymptotic scaling,

$$M^{(r)}M_{}r^\gamma ,\xi ^{(r)}\xi _{}r^{1+\alpha }$$
(25)
for large $`r`$, with some coefficients $`M_{},\xi _{}`$ and exponents $`\gamma ,\alpha `$. Whereas the properties of the form factors are consistent with $`\gamma =1`$ only, the other exponent can only be determined numerically with result $`\alpha =0.27`$. For those readers unfamiliar with this model we point out the amusing fact that the self–similarity holds for all $`r2`$ and thus interrelates the spectral functions of two different isovector operators, the spin field and the conserved vector current. Similarly (25) is equally valid for even and odd $`r`$ values, provided we use the normalization introduced in (1).

Self–similarity continues to be satisfied very well as demonstrated in Figure 2. To test asymptotic scaling, we used the fitted numerical values based on results for up to 6 particles to “predict” for the $`r=7`$ case

$$M_{\mathrm{sc}}^{(7)}=0.03188\mathrm{and}\xi _{\mathrm{sc}}^{(7)}=17.73,$$
(26)
and compared it to the values directly determined from our 7–particle results

$$M^{(7)}=0.03189\mathrm{and}\xi ^{(7)}=17.77.$$
(27)
Similarly for the integrals

$$c^{(r)}=_M^{\mathrm{}}d\mu \rho ^{(r)}(\mu )=_0^{\mathrm{}}dxR^{(r)}(x)h^{(r)}=_0^{\mathrm{}}\frac{\mathrm{d}x}{x}R^{(r)}(x)$$
(28)
we predict

$$c_{\mathrm{sc}}^{(7)}=1.464\mathrm{and}h_{\mathrm{sc}}^{(7)}=0.04879,$$
(29)
and actually get

$$c^{(7)}=1.46(1)\mathrm{and}h^{(7)}=0.0488(1).$$
(30)
## 6 Ultra–positivity

Because of the large cancellation between positive and negative terms in the representation of $`G^{(r)}`$ in terms of symmetric polynomials, it is natural to ask if there exists an alternative representation that is manifestly positive. It is indeed possible to find such representations. If we arrange the rapidities in decreasing order (which is always possible since the polynomial is symmetric under permutations) then all $`u_j=\theta _j\theta _{j+1}`$ rapidity differences are positive and we found (for all available particle numbers $`3r7`$) that if we expand $`G^{(r)}`$ in terms of these differences then all coefficients are positive. Moreover, it is possible to find upper and lower bounds both of which are of the simple form (18) such that, expanded in terms of the rapidity differences $`u_j`$, each and every term in the expansion of $`G^{(r)}`$ has smaller coefficient than the corresponding one of the upper bound $`P^{(r)}(t_u^{(r)})`$ and similarly larger than the coefficient of the corresponding term of the lower bound $`P^{(r)}(t_l^{(r)})`$.

The value of the parameter characterizing the lower bound is uniformly $`t_l^{(r)}=2`$ for all cases $`3r7`$, whereas the $`t_u^{(r)}`$ giving the lowest upper bound are

$$t_u^{(3)}=2,t_u^{(4)}=9/4,t_u^{(5)}=2.434,t_u^{(6)}=2.576,t_u^{(7)}=2.688.$$
(31)
The last three values are numerically approximate 222e.g. the exact value for the case $`r=5`$ is $`t_u^{(5)}=208^{1/6}`$ and come from the requirement that the constant term of $`P^{(r)}(t_u^{(r)})`$ be larger than that of $`G^{(r)}`$. This is clearly necessary, but at least in the available cases for $`r5`$ it is also sufficient 333The case $`r=4`$ is an exception since the requirement in this case yields a value $`(34/3)^{1/3}=2.246`$, which is not quite sufficient..

The integrals over the spin and current spectral densities are given as sums of the integrals defined in (28):

$`C^{\mathrm{spin}}`$ $`=`$ $`{\displaystyle _M^{\mathrm{}}}d\mu \rho ^{\mathrm{spin}}(\mu )={\displaystyle \frac{4}{\pi }}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}c^{(2k+1)},`$
$`C^{\mathrm{curr}}`$ $`=`$ $`{\displaystyle _M^{\mathrm{}}}d\mu \rho ^{\mathrm{curr}}(\mu )={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}c^{(2k+2)}.`$ (32)
An outstanding question is whether these are finite or infinite in the FF construction. Certainly renormalized perturbation theory predicts that $`C^{\mathrm{spin}}=\mathrm{}=C^{\mathrm{curr}}`$. Also the validity of the scaling hypothesis requires this because in this scenario $`c^{(r)}`$ grows as $`r^\alpha `$; the numerical evidence thereof is shown in Figure 3.

If the bounds we found above for $`3r7`$ continue to be true for all $`r`$, then the simple structure of (18) allows to study the structure of the correlation function analytically. The existence of an upper limit of simple form may help proving the existence of the correlation function whereas the existence of the lower bound may facilitate the construction of a proof that $`C^{\mathrm{spin}}=\mathrm{}=C^{\mathrm{curr}}`$ independent of the validity of the scaling hypothesis.

We are not going to discuss these questions further in this paper, but to illustrate the usefulness of the existence of the simple lower bound we have calculated the integrals $`c^{(r)}(2)`$. These are defined analogously to the ones in (28) using $`P^{(r)}(2)`$ instead of $`G^{(r)}`$. The simplicity of the integrand allows us to calculate these integrals numerically quite effortlessly 444This is because one can represent $`c^{(r)}(2)`$ as an integral with the integrand just involving the polynomial $`P_0^{(r)}(2)`$ in Eq. 19 (instead all its permutations), at the cost of extending the limits of integration over the $`u_i`$ from $`\mathrm{}`$ to $`+\mathrm{}`$. up to $`r=11`$. The results of the numerical integration are given in Table 1.

Although the integrals $`c^{(r)}(2)`$ are (apparently) decreasing with $`r`$, it is possible that they are decreasing slow enough to make the series (32) diverge. Indeed, as seen in Figure 3, $`rc^{(r)}(2)`$ is monotonically increasing for all $`r`$ evaluated up to now.

Acknowledgements

J. B. wishes to thank the Max–Planck–Institut für Physik in Munich where most of this work has been done, for hospitality. This investigation was supported in part by the Hungarian National Science Fund OTKA (under T030099 and T034299).