text:
# Vortex lattice formation in a rotating Bose-Einstein condensate

## Abstract

We study the dynamics of vortex lattice formation of a rotating trapped Bose-Einstein condensate by numerically solving the two-dimensional Gross-Pitaevskii equation, and find that the condensate undergoes elliptic deformation, followed by unstable surface-mode excitations before forming a quantized vortex lattice. The origin of the peculiar surface-mode excitations is identified to be phase fluctuations at the low-density surface regime. The obtained dependence of a distortion parameter on time and that on the driving frequency agree with the recent experiments by Madison et al. \[Phys. Rev. Lett. 86, 4443 (2001)\].

Quantized vortices have long been studied in superfluid 4He as the topological defects characteristic of superfluidity Donnelly ; Pismen . However, the relatively high density and strong repulsive interaction complicate the theoretical treatments of the Bose-Einstein condensed liquid, and the healing length of the atomic scale makes the experimental visualization of the quantized vortices difficult. The recent achievement of Bose-Einstein condensation in trapped alkali-metal atomic gases at ultra low temperatures has stimulated intense experimental and theoretical activity. The atomic Bose-Einstein condensates(BECs) have the weak interaction because they are dilute gases, thus being free of the above difficulties that superfluid 4He is subject to. Quantized vortices in the atomic BECs have recently been created experimentally by Matthews et al. JILA , Madison et al. ENS1  and Abo-Shaeer et al. MIT .

By rotating an asymmetric trapping potential, Madison et al. at ENS succeeded in forming a quantized vortex in 87Rb BEC for a stirring frequency that exceeds a critical value ENS1 . Vortex lattices were obtained for higher frequencies. The ENS group subsequently observed that vortex nucleation occurs via a dynamical instability of the condensate ENS2 . For a given modulation amplitude and stirring frequency, the steady state of the condensate was distorted to an elliptic cloud, stationary in the rotating frame, as predicted by Recati et al. Recati . An intrinsic dynamical instability Sinha  of the steady state transformed the elliptic state into a more axisymmetric state with vortices. However, the origin of that instability and how it leads to the formation of vortex lattices remains to be investigated

The ENS group found that the minimum rotation frequency $`\mathrm{\Omega }_{\mathrm{nuc}}`$ at which one vortex appears is 0.65$`\omega _{}`$, where $`\omega _{}`$ is the transverse oscillation frequency of the cigar-shaped trapping potential, independent of the number of atoms or the longitudinal frequency $`\omega _z`$ ENS3 . There is a discrepancy between the observations and the theoretical considerations based on the stationary solution of the Gross-Pitaevskii equation(GPE) Isoshima ; cri ; $`\mathrm{\Omega }_{\mathrm{nuc}}`$ is significantly larger than its equilibrium estimates.

The present paper addresses these issues by numerically solving the GPE that governs the time evolution of the order parameter $`\psi (𝐫,t)`$:

$`(i\gamma )\mathrm{}{\displaystyle \frac{\psi }{t}}=\left[{\displaystyle \frac{\mathrm{}^2}{2m}}^2+V_{\mathrm{tr}}+g|\psi |^2\mu \mathrm{\Omega }L_z\right]\psi .`$ (1)
Here $`g=4\pi \mathrm{}^2a/m`$ is the coupling constant, proportional to the 87Rb scattering length $`a`$5.77 nm. The high anisotropy of the cigar-shaped potential used in the ENS experiments ($`\omega _{}/\omega _z14`$) may permit the two-dimensional analysis. We thus focus on the two-dimensional dynamics of Eq. (1) by assuming the trapping potential

$$V_{\mathrm{tr}}(𝐫)=\frac{1}{2}m\omega _{}^2\{(1+ϵ_x)x^2+(1+ϵ_y)y^2\},$$
(2)
where $`\omega _{}=2\pi \times 219`$ Hz, and the parameters $`ϵ_x=0.03`$ and $`ϵ_y=0.09`$ describe small deviations of the trap from the axisymmetry, corresponding to the ENS experiments ENS1 . The centrifugal term $`\mathrm{\Omega }L_z=i\mathrm{}\mathrm{\Omega }(x_yy_x)`$ appears in a system rotating about the $`z`$ axis at a frequency $`\mathrm{\Omega }`$. An important characteristic parameter of the two-dimensional system is $`C=8\pi Na/L`$, with the total number $`N`$ of the condensate atoms and $`L`$ the size of the system along the $`z`$ axis.

The term with $`\gamma `$ in Eq. (1) introduces the dissipation. Although the detailed mechanism of the dissipation is yet to be understood, we include this term in the GPE because of the following reasons. First, collective damped oscillations of the condensate have been observed by Jin et al. and Mewes et al. MIT2 , which shows the presence of some dissipative mechanisms and is consistent with the solution of the GPE with the dissipative term Choi . Second, even if the trapping potential is rotated fast enough, vortex lattices will never be formed without dissipation, because vortex lattice correspond to local minimum of the total energy in the configuration space rota . In other words, the observation of the vortex lattices ENS1  implies the presence of dissipation. Following Ref. Choi , we use $`\gamma =0.03`$ throughout this work.

In this paper, we focus on two kinds of responses of the condensate, namely, the response to the sudden turn on of the rotation of the potential, and that to the slow turn on. The numerical calculations are performed with the implicit Crank-Nicolson method. The chemical potential $`\mu `$ is adjusted at all times so as to preserve the number of condensate atoms.

We first prepare an equilibrium condensate with $`C=1400`$ trapped in a stationary potential. Figure 1 shows the typical dynamics of the condensate density $`|\psi (𝐫,t)|^2`$ after the potential begins to rotate suddenly with $`\mathrm{\Omega }=0.7\omega _{}`$ HP . The condensate is elongated along the $`x`$ axis because of the small anisotropy of $`V_{\mathrm{tr}}`$ \[Eq. (2)\], and the elliptic cloud oscillates. Then, the boundary surface of the condensate becomes unstable, exciting the surface waves, which propagate along the surface. The excitations are likely to occur on the surface whose curvature is low, i.e., parallel to the longer axis of the ellipse. The ripples on the surface develop into the vortex cores around which superflow circulates. Subject to the dissipative vortex dynamics, some vortices enter the condensate, forming a vortex lattice. As the vortex lattice is being formed, the axial symmetry of the condensate is recovered by transferring angular momentum into quantized vortices.

This peculiar dynamics is understood by investigating the phase of $`\psi (𝐫,t)`$ as shown in Fig. 2. There are some lines where the phase changes discontinuously from black to white, which corresponds to the branch cuts between the phases 0 and 2$`\pi `$. Their ends represent phase defects, i.e. vortices. As soon as the rotation starts, some defects begin to enter. When the defects are on the outskirts of the condensate where the amplitude $`|\psi (𝐫,t)|`$ is almost negligible, they neither contribute to the energy nor the angular momentum of the system. These defects come into the boundary surface of the condensate within which the amplitude grows up. Then the defects compete with each other and induce the above surface waves due to interference. There the selection of the defects starts, because their further invasion into the condensate costs the energy and the angular momentum. As is well known in the study of rotating superfluid helium rota , the rotating drive pulls vortices into the rotation axis, while repulsive interaction tends to push them apart; this competition yields a vortex lattice whose vortex density depends on the rotation frequency. In our case, some vortices enter the condensate and form a lattice dependent on $`\mathrm{\Omega }`$, while excessive vortices are repelled and escape to the outside. Remarkably, the phase profile of Fig. 2(b) reveals that the repelled vortices also form a lattice on the outskirts of the condensate. Since they cannot be seen in the corresponding density profile of Fig. 1(f), they may be called “ghost” vortices.

The distortion of the condensate to an elliptic cloud was theoretically studied by Recati et al. Recati . Assuming the quadrupolar velocity field $`𝐯(𝐫)=\alpha (xy)`$, they obtained the distortion parameter $`\alpha =\mathrm{\Omega }(R_x^2R_y^2)/(R_x^2+R_y^2)`$ in the steady states as a function of $`\mathrm{\Omega }`$, where $`R_x`$ and $`R_y`$ are the sizes of the condensate along each direction. Madison et al. observed that, after the rotation of $`\mathrm{\Omega }=0.7\omega _{}`$ starts suddenly, $`\alpha `$ oscillates during a few hundred milliseconds and then falls abruptly to a value below 0.1 when vortices enter the condensate from its boundary surface ENS2 . Figure 3 shows the oscillation of $`\alpha `$ and the increase of the angular momentum $`\mathrm{}_z/\mathrm{}`$ per atom in our dynamics of Fig. 1. This figure closely resembles Fig. 3 of Ref. ENS2 , and our scenario is thus consistent with the experimental results.

The dispersion relation of the surface waves Khawaja  of a rotating condensate is obtained as follows. The substitution of the Madelung transformation $`\psi =|\psi |e^{i\theta }`$ into Eq. (1) yields

$`{\displaystyle \frac{|\psi |^2}{t}}+\left[|\psi |^2(𝐯\mathrm{\Omega }\times 𝐫)\right]=0,`$ (3a)
$`m{\displaystyle \frac{𝐯}{t}}=\left\{\delta \mu +{\displaystyle \frac{m}{2}}\left(𝐯\mathrm{\Omega }\times 𝐫\right)^2{\displaystyle \frac{m}{2}}\mathrm{\Omega }^2r^2\right\},`$ (3b)
where $`𝐯=(\mathrm{}/m)\theta `$ and $`\delta \mu =V_{\mathrm{tr}}+g|\psi |^2\frac{\mathrm{}^2}{2m|\psi |}^2|\psi |\mu `$. By linearizing Eq. (3), we obtain the dispersion relation of the surface waves $`\omega ^2=R(\omega _{}^2\mathrm{\Omega }^2)k`$ with the condensate size $`R`$. The dependence of $`\omega `$ on $`R`$ shows that the surface with lower curvature is excited more easily, which is clearly seen in Fig. 1. The group velocity $`d\omega /dk=\sqrt{R(\omega _{}^2\mathrm{\Omega }^2)}/(2\sqrt{k})`$ agrees with the propagation velocity of the surface wave in our simulation. As discussed by Al Khawaja et al. Khawaja , the surface waves are connected with the low energy excitations studied by Stringari and Dalfovo et al. Stringari , and Isoshima and Machida Isoshima .

We have examined the critical frequency at which vortices can enter with $`C=420`$ corresponding to the experimental condition ENS1 , and found that only when $`\mathrm{\Omega }`$ is larger than $`\mathrm{\Omega }_{c1}0.57\omega _{}`$ can vortices enter the condensate and form a lattice. This critical frequency is closer to the observed value $`\mathrm{\Omega }_{\mathrm{nuc}}=0.65\omega _{}`$, than the values obtained in previous literature. The critical frequencies have been studied from thermodynamic or stability arguments Isoshima ; cri ; 3D , but they are generally much smaller than the observed value. It should be noted that these critical frequencies only give the necessary condition which enables a vortex to exist stably at the center of the trap. Actually, vortices should be nucleated at the boundaries and come into the condensate; the condition of $`\mathrm{\Omega }`$, which realizes such nonlinear dynamics, may be generally different from that obtained from the stability arguments. Isoshima and Machida examined the local instability of nonvortex states toward vortex states within the Bogoliubov theory Isoshima . They note that this local instability could give higher critical frequencies than the global instability, which comes from comparing their total energy, for example, about $`0.6\omega _{}`$ for the ENS experimental condition Machida . The stability analysis cannot answer what happens when the system goes over the linear region. Our paper reveals the importance of the nonlinear dynamics beyond the stability analysis. The detailed studies of the dependence of the number of vortices on $`\mathrm{\Omega }`$, that of $`\mathrm{\Omega }_{c1}`$ on the number of condensate atoms Isoshima  will be reported elsewhere.

Finally, we study the evolution of the condensate in a potential rotating with a time-dependent frequency $`\mathrm{\Omega }(t)`$. The experiments were made by Madison et al. ENS2 . The condensate at every moment was in a stationary state with an elliptic cloud obtained by Recati et al. under the Thomas-Fermi limit Recati , since the time dependence was slow enough. One of the important observations was that, for an ascending ramp, the condensate follows branch I of the nonvortex state until $`\mathrm{\Omega }_c0.75\omega _{}`$. Beyond this critical frequency, the condensate changes from branch I to the more axisymmetric state, nucleating vortices. This critical frequency agrees with the frequency at which branch I becomes dynamically unstable Sinha ; however, the detailed nature of the bifurcation has remained to be clarified. We study the response of the condensate of $`C=420`$ for the ascending frequency $`\mathrm{\Omega }(t)=ct`$ with $`c=2\pi \times 400`$ Hz/s; they are about the same values as those used in the ENS experiment. The obtained results are consistent with the experimental ones, as shown in Fig. 4. As $`\mathrm{\Omega }`$ increases, the condensate becomes gradually elliptic following branch I, when the nucleated “ghost” vortices come to the surface of the condensate. When $`\mathrm{\Omega }`$ is raised to the critical value $`\mathrm{\Omega }_{c2}0.75\omega _{}`$, the vortices jostle with each other and excite the surface waves, reflecting the dynamical instability. Then, similar to the first set of calculations, some vortices enter the condensate, form a lattice and increase the angular momentum, while the condensate recovers the axisymmetry apart from branch I. The small discrepancy of the trace of our calculation from branch I in Fig. 4 may be attributed to a deviation from the Thomas-Fermi limit; the similar discrepancy is reported by the Madison et al. too ENS2 .

It is necessary to compare two critical frequencies $`\mathrm{\Omega }_{c1}0.57\omega _{}`$ and $`\mathrm{\Omega }_{c2}0.75\omega _{}`$. What causes the difference? In the first set of calculations where the potential begins to rotate abruptly with a frequency $`\mathrm{\Omega }`$, the condensate first deviates much from the stationary state $`\alpha (\mathrm{\Omega })`$ with the quadrupolar velocity field. As shown in Fig. 3, it seemingly relaxes to that state through the damped oscillation of the large amplitude, but it changes on the way to another state with a vortex lattice which has lower energy than the nonvortex state $`\alpha (\mathrm{\Omega })`$. The critical frequency $`\mathrm{\Omega }_{c1}`$ means that the state with vortices can be realized dynamically from the nonvortex state only when $`\mathrm{\Omega }>\mathrm{\Omega }_{c1}`$. In the second calculation, as $`\mathrm{\Omega }`$ increases slowly, the condensate follows quasistatically the nonvortex state $`\alpha (\mathrm{\Omega })`$ of branch I, even after it becomes metastable for $`\mathrm{\Omega }>\mathrm{\Omega }_{c1}`$. Only when $`\mathrm{\Omega }`$ exceeds $`\mathrm{\Omega }_{c2}`$, the condensate can change from this state to the vortex state, reflecting the dynamical instability Sinha .

Recently García-Ripoll et al. and Feder et al. study this system by three-dimensional analysis of the GPE 3D , taking account of the bending of vortex lines in the cigar-shaped potential. A bent vortex moves with the self-induced velocity proportional to its curvature Donnelly , which is lacking in the two-dimensional dynamics. However, due to the high anisotropy of the cigar-shaped potential, vortices are aligned almost straight along the symmetry axis. It is therefore unlikely that the bending effect alters the results of our work qualitatively and significantly.

In conclusion, we have studied the dynamics of vortex lattice formation of a rotating trapped BEC by numerically solving the two-dimensional GPE, and obtained the following physical picture. When the trapping potential begins to rotate at a sufficiently fast frequency, the condensate is distorted to an elliptic shape and oscillates. Then, its boundary surface becomes unstable, exciting surface waves. The origin of these ripples is identified to be the interference of violent phase fluctuations that occur on the outskirts of the condensate; and some of the surface ripples develop into the vortex cores which then form a vortex lattice. The critical frequency for this process is found to be $`\mathrm{\Omega }_{c1}0.57\omega _{}`$ for the condition at the ENS group. On the other hand, if the rotation frequency is raised slowly from zero, the condensate becomes elliptic and follows the nonvortex elliptic state with the quadrupolar flow Recati . When $`\mathrm{\Omega }`$ exceeds $`\mathrm{\Omega }_{c2}0.75\omega _{}`$, the condensate deviates from the state, nucleating vortices via the dynamical instability Sinha . The whole dynamics is concerned with the “ghost” vortices, which are the phase defects outside the condensate. The obtained results are consistent with the experimental results of the ENS group ENS1 ; ENS2 .

The authors thank K. Machida for useful discussions. M.U. acknowledges support by a Grant-in-Aid for Scientific Research (Grant No. 11216204) by the Ministry of Education, Science, Sports, and Culture of Japan, and by the Toray Science Foundation.


format: text

text:
# DESY 01-019 April 2001 False Vacuum Decay after Inflation

## Abstract

Inflation is terminated by a non-equilibrium process which finally leads to a thermal state. We study the onset of this transition in a class of hybrid inflation models. The exponential growth of tachyonic modes leads to decoherence and spinodal decomposition. We compute the decoherence time, the spinodal time, the size of the formed domains and the homogeneous classical fields within a single domain.

An inflationary phase in the early history of the universe can explain its present flatness and homogeneity as well as the anisotropy of the cosmic microwave background . Inflation has to terminate in some non-equilibrium process which eventually leads to a thermal state. This is elegantly realized in hybrid models of inflation where the time evolution of the inflaton field $`\phi `$ controls the transition from a false vacuum state to the true vacuum . Hybrid inflation arises naturally in supersymmetric theories with spontaneous symmetry breaking . The dynamics of the inflaton field during the inflationary phase is then determined either by supergravity corrections , by quantum corrections  or by the effects of the supersymmetry breaking sector of the theory .

The phase transition at the end of inflation is of crucial importance since it determines the initial conditions for the further evolution towards the thermal state. Important quantities, which have to be obtained, are the reheating temperature and the initial abundances of weakly coupled particles, such as gravitinos, which are not in thermal equilibrium. In hybrid inflation the decay of the false vacuum takes place via spinodal decomposition , similar to a second-order phase transition. Such a decay has already been discussed in the literature in the context of new inflation  and also for scalar $`\varphi ^4`$-theory in Minkowsky space .

In the following we shall study the onset of this transition in a class of hybrid inflation models , where the evolution is controlled by the slow-roll motion of the inflaton field. This will allow us to follow in detail the breakup of the homogeneous inflationary phase into domains and to compute the decoherence time, the spinodal time, the average size of the formed domains and the homogeneous classical fields inside a single domain.

False vacuum inflation

Hybrid inflation can be based on the superpotential

$`W=W_G+W_S=\lambda T\left(M_G^2\mathrm{\Sigma }^2\right)+M_S^2\left(\beta +S\right).`$ (1)
Here $`W_G`$ describes the supersymmetric breaking of a global symmetry, whereas the Polonyi superpotential $`W_S`$ accounts for supersymmetry breaking . The discrete $`Z_2`$ symmetry of the model leads to the well known domain wall problem . We consider this case only for simplicity; the extension to a continuous $`U(1)`$ symmetry, where no such problems exist, is straightforward. The mass $`M_S`$ is the scale of supersymmetry breaking, associated with the gravitino mass $`m_{3/2}=M_S^2/M_\mathrm{P}\mathrm{exp}(2\sqrt{3})`$, and the constant $`\beta (2\sqrt{3})M_\mathrm{P}`$, where $`M_\mathrm{P}=(8\pi G_N)^{1/2}=2.4\times 10^{18}\text{GeV}`$ is the Planck mass. Successful inflation requires $`M_S^2\lambda M_G^2`$. A rotation of the fields $`T`$ and $`S`$ by the small angle $`\xi M_S^2/(\lambda M_G^2)`$ changes the superpotential (1) to the standard form of an O’Raifeartaigh model .

In the case of global supersymmetry one obtains from the superpotential (1) the scalar potential

$`V_0=\lambda ^2|M_G^2\mathrm{\Sigma }^2|^2+4\lambda ^2|T|^2|\mathrm{\Sigma }|^2+M_S^4,`$ (2)
with the ground state

$`T=0,|\mathrm{\Sigma }|=M_G.`$ (3)
For $`|T|^2>M_G/2`$, $`V_0`$ has a local minimum at $`\mathrm{\Sigma }=0`$ where it is flat in $`S`$ and $`T`$.

Supergravity corrections are included in the potential

$$V=e^{K/M_\mathrm{P}^2}\left[\left|\frac{W}{z_i}+\frac{z_i^{}W}{M_\mathrm{P}^2}\right|^23\frac{|W|^2}{M_\mathrm{P}^2}\right],$$
(4)
where $`z_i=T,\mathrm{\Sigma },S`$, $`i=1\mathrm{}3`$. We choose $`K`$ to be the canonical Kähler potential, $`K=_i|z_i|^2`$, which is singled out in $`N=2`$ supersymmetry . Up to corrections $`𝒪(1/M_\mathrm{P}^n)`$, the expectation values of the three fields $`T`$, $`\mathrm{\Sigma }`$ and $`S`$ in the true vacuum are

$$T=0,|\mathrm{\Sigma }|=M_G,S=(2\sqrt{3})M_\mathrm{P}.$$
(5)
The supergravity corrections distinguish between the real and the imaginary part of $`T(\phi +i\chi )/\sqrt{2}`$. For $`\phi >\phi _c=M_G`$, $`V`$ has again a local minimum at $`\mathrm{\Sigma }=(\sigma +i\rho )/\sqrt{2}=0`$. The energy density is then $`V\lambda ^2M_G^4`$, which leads to an inflationary phase with Hubble parameter $`H=\lambda M_G^2/(\sqrt{3}M_\mathrm{P})`$. Here the potential is almost flat in $`\phi `$ which plays the role of the inflaton. The mass of $`\rho `$ is large, $`m_\rho 2\lambda \phi >2\lambda M_G`$. The minimum with respect to $`S`$ and $`\chi `$ lies at $`S=\chi =0`$, but the corresponding masses are small, $`m_S=𝒪(H)`$ and $`m_\chi =𝒪(H\phi /M_\mathrm{P})`$.

In this paper we are interested in the decay of the false vacuum driving inflation, which takes place near the critical point $`\phi _c=M_G`$. This phase transition is determined by the dynamics of the lightest scalar fields, the inflaton $`\phi `$ and $`\sigma `$, which becomes massless at $`\phi _c`$. The scalar potential for $`\phi `$ and $`\sigma `$ reads, up to corrections $`𝒪(1/M_\mathrm{P}^2)`$,

$`V(\phi ,\sigma )=\lambda ^2\left(M_G^2{\displaystyle \frac{\sigma ^2}{2}}\right)^2+\lambda ^2\phi ^2\sigma ^2+\kappa \phi .`$ (6)
The linear term in the potential is $`𝒪(1/M_\mathrm{P})`$, with

$`\kappa =2^{3/2}(2\sqrt{3})\lambda {\displaystyle \frac{M_S^2M_G^2}{M_\mathrm{P}}}0.6\lambda m_{3/2}M_G^2.`$ (7)
For a certain range of parameters the time evolution of the inflaton field during inflation is determined by the linear term in the potential . Representative values are $`M_S1.6\times 10^{10}\text{GeV}`$, which implies $`m_{3/2}130\text{GeV}`$, $`\lambda 1.3\times 10^5`$ and $`M_G10^{16}\text{GeV}`$. This yields $`H3.1\times 10^8\text{GeV}`$ and $`\kappa (4.6\times 10^9\text{GeV})^3`$. From the full supergravity potential (4) one obtains for the slow-roll conditions of the inflaton $`\phi `$, as long as $`\sigma ^2`$ is sufficiently small,

$`ϵ`$ $`=`$ $`{\displaystyle \frac{M_\mathrm{P}^2}{2}}\left({\displaystyle \frac{V^{}}{V}}\right)^2={\displaystyle \frac{4\xi ^2\beta ^2}{M_\mathrm{P}^2}}\left(1{\displaystyle \frac{\phi ^3}{2\sqrt{2}\xi \beta M_\mathrm{P}^2}}+\mathrm{}\right)1,`$ (8)
$`\eta `$ $`=`$ $`M_\mathrm{P}^2{\displaystyle \frac{V^{\prime \prime }}{V}}={\displaystyle \frac{3}{2}}{\displaystyle \frac{\phi ^2}{M_\mathrm{P}^2}}+\mathrm{}1.`$ (9)
These conditions are clearly satisfied for values of $`\phi M_\mathrm{P}`$ and in particular below the critical value $`\phi _c`$. The solution of the equation of motion for the homogeneous classical field reads,

$`\phi _{cl}(t)=M_G{\displaystyle \frac{\kappa }{3H}}t,`$ (10)
where we have fixed the origin of time by $`\phi (0)=\phi _c=M_G`$. The transition from the inflationary phase to the thermal phase starts at $`t=0`$.

Fluctuations near the critical point

At $`\phi =\phi _c`$ the mass term of the $`\sigma `$-field changes sign. This is analogous to a second-order phase transition with $`\phi (t)`$ playing the role of a time dependent temperature of the thermal bath. In flat de Sitter space with metric

$`ds^2=dt^2e^{2Ht}d\stackrel{}{x}^2,`$ (11)
the quantum fields $`\phi (t,\stackrel{}{x})`$ and $`\sigma (t,\stackrel{}{x})`$ can be expanded into plane waves which are momentum eigenstates,

$`\phi (t,\stackrel{}{x})`$ $`=`$ $`\phi _{cl}(t)+e^{\frac{3}{2}Ht}{\displaystyle \frac{d^3k}{(2\pi )^{3/2}}\left(a_\phi (\stackrel{}{k})\phi _k(t)e^{i\stackrel{}{k}\stackrel{}{x}}+a_\phi ^{}(\stackrel{}{k})\phi _k^{}(t)e^{i\stackrel{}{k}\stackrel{}{x}}\right)},`$ (12)
$`\sigma (t,\stackrel{}{x})`$ $`=`$ $`e^{\frac{3}{2}Ht}{\displaystyle \frac{d^3k}{(2\pi )^{3/2}}\left(a_\sigma (\stackrel{}{k})\sigma _k(t)e^{i\stackrel{}{k}\stackrel{}{x}}+a_\sigma ^{}(\stackrel{}{k})\sigma _k^{}(t)e^{i\stackrel{}{k}\stackrel{}{x}}\right)}.`$ (13)
Note the absence of $`\sigma _{cl}(t)`$, the homogeneous part of $`\sigma (t,\stackrel{}{x})`$. This is due to the initial condition $`\sigma _{cl}(0)=0`$ at the end of inflation. The annihilation and creation operators satisfy canonical commutation relations,

$`[a_\phi (\stackrel{}{k}),a_\phi ^{}(\stackrel{}{k^{}})]=[a_\sigma (\stackrel{}{k}),a_\sigma ^{}(\stackrel{}{k^{}})]=\delta ^3(\stackrel{}{k}\stackrel{}{k^{}}),`$ (14)
and the canonically conjugate operators to the fields $`\phi `$ and $`\sigma `$ are

$`\pi _\phi (t,\stackrel{}{x})=e^{3Ht}\dot{\phi }(t,\stackrel{}{x}),\pi _\sigma (t,\stackrel{}{x})=e^{3Ht}\dot{\sigma }(t,\stackrel{}{x}).`$ (15)
The assumption of a constant Hubble parameter is well justified during inflation and the onset of the subsequent phase transition.

Near the critical point $`\phi =\phi _c`$ the potential for the quantum fields $`\phi `$ and $`\sigma `$ takes the form,

$`V(\phi ,\sigma )={\displaystyle \frac{1}{2}}D^3t\sigma ^2+\mathrm{}.`$ (16)
where $`D^3=4\lambda ^2\kappa M_G/(3H)`$. For the parameters given above one has $`D=3.4H`$. Note, that this equation also holds in other inflationary models where the inflaton is slowly rolling near $`\phi =\phi _c`$, since $`D^3t`$ is the first term in an expansion of the $`\sigma `$-mass around $`t=0`$.

The mode functions $`\phi _k(t)`$ and $`\sigma _k(t)`$ satisfy the equations of motion,

$`\ddot{\phi _k}+\left(k^2e^{2Ht}{\displaystyle \frac{9}{4}}H^2\right)\phi _k=0,`$ (17)
$`\ddot{\sigma _k}+\left(k^2e^{2Ht}{\displaystyle \frac{9}{4}}H^2D^3t\right)\sigma _k=0.`$ (18)
The solution of eq. (17) is well known from studies of density fluctuations during inflation. It is given by a Hankel function

$`\phi _k(t)={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\pi }{H}}\right)^{1/2}\text{H}_{3/2}^{(1)}(z),`$ (19)
where $`z=k\mathrm{exp}(Ht)/H`$ and

$`\text{H}_{3/2}^{(1)}(z)=\left({\displaystyle \frac{2}{\pi z}}\right)^{1/2}e^{iz}{\displaystyle \frac{z+i}{z}}.`$ (20)
The solution of eq. (18) has to be determined numerically. However, asymptotically, for $`tt_k`$ and $`tt_k`$, where

$`k^2e^{2Ht_k}{\displaystyle \frac{9}{4}}H^2=|D^3t_k|,`$ (21)
useful approximate analytical solutions can be found111For momenta $`3H/2<k<k_{max}(D^3/(2He))^{1/2}3H`$ eq. (21) has more than one solution for $`t_k`$. We choose the smallest value for $`t_k`$ to construct an approximate solution. For a more detailed discussion, see .. At large negative times we again have the Hankel function analogous to (19) as solution. This is the de Sitter-space initial condition.

For $`tt_k`$ eq. (18) reduces to

$`\ddot{\sigma _k}D^3t\sigma _k=0.`$ (22)
which describes a scalar field in Minkowski space with mass squared $`\mu ^2(t)=D^3t`$, decreasing with time. The solutions are given by the Airy functions Ai and Bi . An approximate analytic solution of eq. (18) can be obtained by matching Airy and Hankel functions at $`t=t_k`$,

$`\sigma _k(t)=C_A(k)\text{Ai}(z)+C_B(k)\text{Bi}(z),`$ (23)
where $`z=Dt`$. For soft momenta, $`k<H`$, one finds

$`C_A=i\sqrt{{\displaystyle \frac{\pi }{2D}}},C_B=\sqrt{{\displaystyle \frac{\pi }{2D}}}.`$ (24)
For $`z=Dt1`$ the behaviour of $`\phi _k(t)`$ is dominated by the asymptotic growth of $`\text{Bi}(z)`$,

$`\text{Bi}(z){\displaystyle \frac{1}{\sqrt{\pi }}}z^{1/4}\mathrm{exp}\left({\displaystyle \frac{2}{3}}z^{3/2}\right).`$ (25)
The physical meaning of such tachyonic instabilities was studied in detail by Guth and Pi for fixed negative mass squared, $`\mu ^2(t)=const`$, in the context of new inflation . More recently, the dynamics of spinodal decomposition has been investigated by several groups .

The situation considered here is similar to the previously discussed cases with respect to the occurrence of spinodal decomposition. However, a crucial difference is that the negative curvature of the potential is neither constant nor suddenly ‘turned on’. The change from positive to negative curvature is, on the contrary, controlled by the slow-roll of the inflaton field across the critical point $`\phi =\phi _c`$.

An important quantity is the growth of the total fluctuation, the variance

$`\sigma ^2(t)_r=\sigma ^2(t)\sigma ^2(0),`$ (26)
where

$`\sigma ^2(t)={\displaystyle \frac{d^3k}{(2\pi )^3}e^{3Ht}|\sigma _k(t)|^2}{\displaystyle 𝑑kf_k(t)}.`$ (27)
Results for the spectral function $`f_k(t)`$ are shown in Fig. 1. The expected rapid growth with time is clearly visible. For small $`k`$ one has $`f_kk^2`$, since $`|C_B(k)|const`$. For larger momenta, $`k>k_{}3H`$, also the time $`t_k`$ is large, which leads to a delayed growth of these modes and a rapid fall-off of $`f_k(t)`$ at large $`k`$. For very large momenta, $`kk_{}`$, one sees the behaviour of a massless field, i.e. $`f_kk`$, which is redshifted with increasing time due to the Hubble expansion.

The growth of the variance $`\sigma ^2(t)_r`$ is shown in Fig. 2. It clearly exhibits the characteristic feature of our model, a faster than exponential growth with time. Eventually, this growth is terminated by the back reaction of the produced fields. The spinodal time at which the formation of domains is completed can be estimated by comparing the variance with the global minimum of the potential,

$`\sigma ^2(t_{sp})_r2M_G^2.`$ (28)
From Fig. 2 one reads off $`t_{sp}3/H`$.

Tachyonic modes $`\sigma _k`$ become classical at the decoherence time $`t_{dec}`$. At later times the system is characterized by a classical probability distribution $`f(\sigma _k,t)`$ for the amplitude $`\sigma _k`$. The homogeneous mode of the system $`\sigma _{cl}(t)=\sigma (t,\stackrel{}{x})`$ remains zero after spinodal decomposition .

The decoherence time $`t_{dec}`$ is determined by the requirement that the tachyonic modes become classical, i.e., the product $`|\sigma _k(t)\pi _{\sigma k}(t)|`$ has to be larger than $`1/2`$, the minimal value obtained for an oscillating mode . For large times, $`Dt1`$, the product is determined by the exponentially growing Airy function $`\text{Bi}(z)`$,

$`|\sigma _k(t)\pi _{\sigma k}(t)|{\displaystyle \frac{1}{\pi }}D|C_B(k)|^2\mathrm{exp}\left({\displaystyle \frac{4}{3}}(Dt)^{3/2}\right).`$ (29)
For soft modes $`k<H`$, the condition $`|\sigma _k(t_{dec})\pi _{\sigma k}(t_{dec})|=R_{dec}1`$ yields the decoherence time

$`t_{dec}{\displaystyle \frac{1}{D}}\left({\displaystyle \frac{3}{4}}\mathrm{ln}(2R_{dec})\right)^{2/3}.`$ (30)
For $`R_{dec}=100`$ this implies $`t_{dec}3/D1/H`$. Results for the decoherence time obtained from the numerical integration are shown in Fig. 3.

From eq. (10) and from Fig. 3 one obtains the values of the inflaton field $`\phi `$ and the variance $`\sigma ^2^{1/2}`$ at the decoherence time $`t_{dec}3/D`$,

$`M_G\phi _{cl}(t_{dec})10^{12}\text{GeV},\sigma ^2(t_{dec})_r^{1/2}10^9\text{GeV}.`$ (31)
These values can serve as initial conditions for the further time evolution of the system . Note, that they are different for the two fields and much larger than the frequently assumed initial condition $`H/(2\pi )`$.

Formation of domains

The rapid growth of tachyonic modes leads to a breakup of the homogeneous inflationary phase into domains. This spinodal decomposition is conveniently described by means of fields averaged over a volume $`V_l=l^3`$ characterized by some smearing function . For simplicity we shall use just a momentum cutoff ($`k_l=2\pi /l`$),

$`\sigma _l(t,\stackrel{}{x})=e^{\frac{3}{2}Ht}{\displaystyle _{k<k_l}}{\displaystyle \frac{d^3k}{(2\pi )^{3/2}}}\left(a_\sigma (\stackrel{}{k})\sigma _k(t)e^{i\stackrel{}{k}\stackrel{}{x}}+a_\sigma ^{}(\stackrel{}{k})\sigma _k^{}(t)e^{i\stackrel{}{k}\stackrel{}{x}}\right).`$ (32)
One then obtains for the correlation function of the smeared operator ($`x=|\stackrel{}{x}|`$),

$`\sigma _l(t,\stackrel{}{x})\sigma _l(t,\stackrel{}{0})={\displaystyle \frac{1}{2\pi ^2x}}e^{3Ht}{\displaystyle _{k<k_l}}𝑑kk|\sigma _k(t)|^2\mathrm{sin}kx.`$ (33)
At time $`t`$ the physical distance and the physical momentum are given by $`y=x\mathrm{exp}(Ht)`$ and $`p=k\mathrm{exp}(Ht)`$, respectively.

An analytic estimate for the correlation function can be obtained based on the matching of Hankel and Airy functions described above. Consider first the case of short times, $`Dt<1`$, and short distances, $`1/k_l<x1/H`$. Here one finds,

$`\sigma _l(t,\stackrel{}{x})\sigma _l(t,\stackrel{}{0}){\displaystyle \frac{1}{4\pi ^2y^2}}\left(\mathrm{cos}(k_lx)1\right).`$ (34)
This is the result expected for a massless scalar field in Minkowsky space. The wavelength of the oscillation is given by the ultraviolet cutoff, $`l=2\pi /k_l`$.

For large times, $`Dt1`$, the correlation function is dominated by the tachyonic modes with $`k<k_{}`$. For distances $`x<1/k_{}`$ the deviation from homogeneity is small,

$`\sigma _l(t,\stackrel{}{x})\sigma _l(t,\stackrel{}{0})\sigma ^2(t)_r\left(1+𝒪((k_{}x)^2)\right).`$ (35)
For large distances, $`x>1/k_{}`$, the correlation function rapidly falls off.

Numerical results for the correlation function are shown in Fig. 4 for two different values of the ultraviolet cut-off $`k_l`$. For a large cut-off, $`k_l=10H`$, one reads off the correlation length $`l_{corr}1/(3H)1/k_{}`$, which corresponds to the average size of a domain. For cut-offs below $`1/l_{corr}`$, e.g. $`k_l=H`$, the expected decrease of the correlation function is observed. The energy density stored in spatial variations of the $`\sigma `$-field is roughly given by $`\rho _{grad}M_G^2/l_{corr}^2`$, which is much smaller than the total energy density $`V\lambda ^2M_G^4`$.

So far we have ignored all back reaction effects. They will terminate the slow-roll motion of the inflaton field, damp the growth of the tachyonic modes and modify the spinodal time $`t_{sp}`$. However, we do not expect a change of our qualitative picture with small decoherence time, $`t_{dec}<t_{sp}`$, and asymmetric classical homogeneous fields $`\phi `$ and $`\sigma `$ larger than $`H/(2\pi )`$. The determination of the properties of the final thermal state, in particular the reheating temperature and the energy fraction of topological defects, requires further work.

Thanks to the slow-roll motion of the inflaton field $`\phi `$ across the critical point $`\phi _c`$ we have been able to follow in detail the breakup of the homogeneous inflationary phase into domains. The tachyonic modes of the field $`\sigma `$, which become classical before the spinodal decomposition ends, lead to an average domain size O(1/H). Within a single domain they appear as a homogeneous classical field.


text:
# Spectral gap of doubly stochastic matrices generated from equidistributed unitary matrices

## 1 Introduction

Recently it was proposed by Kottos and Smilansky  that quantum graphs can serve as a good “toy” model for quantum chaos. Numerical simulations confirmed that for some classes of graphs, including the fully connected graphs, the Bohigas-Giannoni-Schmit Conjecture  holds true. By this we mean that various statistical functions of the spectrum of the graphs converge to the corresponding functions obtained in Random Matrix Theory (RMT)  in the limit as the number of vertices is taken to infinity.

On the graphs, the quantum evolution operator assumes the form of a finite unitary matrix $`𝐔`$. In  it was argued that the classical counterpart of the quantum system is the Markov chain on the directed bonds of the graph with the matrix $`𝐌`$ of the transition probabilities given by

$$M_{b,b^{}}=|U_{b,b^{}}|^2.$$
(1)
It was conjectured that there is a link between ergodicity of the Markov chain and the convergence of the quantum statistics to RMT results.

However, in  it was shown that the statistics for a special family of graphs, the quantum star graphs, converge to a non-RMT limit as the size of the graphs increases. However “classically” the star graphs have nice properties: the corresponding Markov chain possesses a unique attracting equilibrium distribution which corresponds to the eigenvalue 1 of $`𝐌`$. Thus ergodicity of each graph does not immediately imply RMT-like behaviour.

It is well-known that a good quantitative measure of ergodicity of a Markov chain is its spectral gap. Given the spectrum $`\{1,\lambda _2,\lambda _3,\mathrm{}\lambda _N\}`$ of $`𝐌`$, the spectral gap is defined as

$$g=1\underset{i=2,\mathrm{}N}{\mathrm{max}}|\lambda _i|.$$
(2)
The spectral gap indicates the speed of convergence of any initial distribution to the equilibrium one: the greater $`g`$ is, the faster is the convergence. Tanner  suggested that if for a sequence of graphs the spectral gap of the corresponding Markov chains is uniformly bounded away from zero, then the spectral statistics will converge to their RMT form. This seems to be a good indicator of the conformance to RMT: while for each star graph the gap is non-zero, it decreases to zero as we increase the size of the graph. And in the numerically verified cases of RMT-like behaviour the gap stays bounded away from zero.

To investigate the connection between random unitary matrices and the ergodicity of their “classical” analogues, we ask the following question: given a large random unitary matrix, find the probability distribution of the spectral gap of the corresponding Markov chain.

## 2 Preliminaries

Since the most natural ensemble of the unitary matrices is the Circular Unitary Ensemble (CUE), we shall restrict our attention to it, although the generalisation of the problem to other ensembles is straightforward.

###### Definition 1.

CUE$`(N)`$ is defined as the ensemble of all unitary $`N\times N`$ matrices endowed with the probability measure that is invariant under every automorphism

$$𝐔\mathrm{𝐕𝐔𝐖},$$
(3)
where $`𝐕`$ and $`𝐖`$ are any two $`N\times N`$ unitary matrices.

###### Definition 2.

An entrywise nonnegative $`N\times N`$ matrix $`𝐌`$ is called doubly stochastic if

$$\underset{i=1}{\overset{N}{}}M_{i,j}=1j\text{ and }\underset{j=1}{\overset{N}{}}M_{i,j}=1i.$$
(4)
The set of all such $`N\times N`$ matrices we denote by DS$`(N)`$.

If $`𝐔`$ is unitary then the matrix $`𝐌`$ defined by (1) is doubly stochastic. For a Markov chain it means that the uniform measure is invariant, that is $`(1/N,\mathrm{},1/N)`$ is a left eigenvector of $`𝐌`$ with the eigenvalue 1.

###### Remark 1.

It is easy to check that if $`𝐀,𝐁\text{DS}(N)`$ then $`\mathrm{𝐀𝐁}\text{DS}(N)`$, thus DS$`(N)`$ forms a semigroup.

Now we define the function $`S:\text{CUE}(N)\text{DS}(N)`$ to map $`𝐔`$ to $`𝐌`$ according to (1). It is clear that $`S`$ induces a probability measure on the semigroup DS$`(N)`$. Thus our question is

> what is the probability distribution of the spectral gap $`g`$ if the matrices $`M`$ are selected randomly with the probability induced by the correspondence $`S`$.

We denote the gap probability density function by $`f_N(g)`$ and the corresponding cumulative distribution by $`F_N(g)=_0^gf_N(g^{})𝑑g^{}`$. For convenience we also denote by $`\lambda _2(𝐌)`$ the second largest (by modulus) eigenvalue of the doubly stochastic matrix $`𝐌`$,

$$|\lambda _2(𝐌)|=\underset{i=2,\mathrm{}N}{\mathrm{max}}|\lambda _i|$$
(5)
and, as a consequence, $`g=1|\lambda _2(𝐌)|`$.

## 3 Numerics and a conjecture

First of all we are going to present the results of some numerical simulations. The algorithm we used is simple: random matrices from CUE$`(N)`$ are generated using the Hurwitz parametrisation  (see also ), the spectrum of the corresponding doubly stochastic matrix $`M`$ is computed and the gap is calculated according to the definition (2).

The cumulative distribution functions obtained for $`N=2`$, $`3`$, $`5`$, $`10`$, $`20`$, $`40`$, $`80`$ for sample sizes of $`10^5`$ are shown on Figure 1. From the plot we can see that the typical gap tends to 1 as we increase $`N`$. Thus it is natural to conjecture that the gap probability density function $`f_N(g)`$ converges to $`\delta (g1)`$ in the sense of distributions as $`N\mathrm{}`$ or, in other words,

$$\mathrm{Pr}\left\{0ga\right\}0\text{as }N\mathrm{}a<1.$$
(6)
Further evidence in support of this guess is presented by Figure 2, where the mean $`𝖤|\lambda _2|`$ and the standard deviation $`\sigma `$ of the second largest eigenvalue are plotted as functions of $`N`$. A good fit to the mean is given by $`𝖤|\lambda _2|1/\sqrt{N}`$ (which is supported by heuristic observations presented in Section 4). Note that since $`0|\lambda _2|1`$ and

$$𝖤|\lambda _2|a\mathrm{Pr}\left\{a|\lambda _2|1\right\}+0\mathrm{Pr}\left\{0|\lambda _2|a\right\}=a\mathrm{Pr}\left\{0ga\right\},$$
(7)
to verify (6) it is enough to show that $`𝖤|\lambda _2|0`$ as $`N\mathrm{}`$.

## 4 Verification of the conjecture

### 4.1 $`N=2`$ case

In the case of $`N=2`$ it is easy to calculate the gap probability explicitly. For completeness we include the derivation although it is obvious from Figure 1 that the distribution for $`N=2`$ is not typical. A $`2\times 2`$ unitary matrix can be parametrized as

$$𝐔=e^{i\alpha }\left(\begin{array}{cc}e^{i\psi }\mathrm{cos}\varphi & e^{i\chi }\mathrm{sin}\varphi \\ e^{i\chi }\mathrm{sin}\varphi & e^{i\psi }\mathrm{cos}\varphi \end{array}\right),$$
(8)
with

$$0\alpha 2\pi ,0\chi 2\pi ,0\psi 2\pi ,0\varphi \pi /2.$$
(9)
The Haar (uniform) measure on the CUE(2) can be expressed in terms of the parameters $`\alpha `$, $`\chi `$, $`\psi `$ and $`\varphi `$ as $`(2\pi )^3d\alpha d\chi d\psi d\mathrm{sin}^2\varphi `$. The corresponding doubly stochastic matrix is

$$𝐌=\left(\begin{array}{cc}\mathrm{cos}^2\varphi & \mathrm{sin}^2\varphi \\ \mathrm{sin}^2\varphi & \mathrm{cos}^2\varphi \end{array}\right)$$
(10)
with eigenvalues 1 and $`\mathrm{cos}^2\varphi \mathrm{sin}^2\varphi `$. Thus the gap is equal to $`2\mathrm{min}(1\mathrm{sin}^2\varphi ,\mathrm{sin}^2\varphi )`$ and its distribution is uniform on the interval $`[0,1]`$ (which agrees with Figure 1).

### 4.2 Mean of the second largest eigenvalue

###### Theorem 1.

Let $`\lambda _2(𝐌)`$ be the second largest (by modulus) eigenvalue of the random doubly stochastic matrix $`𝐌`$ drawn from $`\mathrm{DS}(N)`$ with the probability induced by the correspondence $`S:\mathrm{CUE}(N)\mathrm{DS}(N)`$. Then

$$𝖤|\lambda _2|0\text{as }N\mathrm{}.$$
(11)
As a consequence, the gap distribution $`f_N(g)`$ converges to $`\delta (g1)`$ in the sense of distributions as $`N\mathrm{}`$.

###### Proof.

The eigenvalues of $`𝐌\text{DS}(N)`$ are, in general, complex. It would be more convenient to deal with nonnegative real eigenvalues. The trick is to consider the matrix $`𝐀=𝐌^T𝐌`$ which, due to the semigroup property is also doubly stochastic. In addition, it is symmetric (all eigenvalues are real) and positive definite (i.e. all eigenvalues are nonnegative). To relate the eigenvalues of the matrices $`𝐀`$ and $`𝐌`$ we write for the second largest eigenvalue of $`𝐀`$,

$$\lambda _2(𝐀)=\underset{|x|=1,(x,e)=0}{\mathrm{max}}(𝐀x,x)(𝐀y,y)=(𝐌y,𝐌y)=|\lambda _2(𝐌)|^2,$$
(12)
where $`e=(1,\mathrm{},1)^T`$ is the eigenvector of both $`𝐀`$ and $`𝐌`$ with the eigenvalue 1, and $`y`$ is the eigenvector of $`𝐌`$ corresponding to $`\lambda _2(𝐌)`$. From (12) and the inequality $`(𝖤\xi )^2<𝖤(\xi ^2)`$ we infer that to prove the Theorem it is enough to show that the mean of $`|\lambda _2(𝐀)|`$ converges to zero as $`N\mathrm{}`$.

The simplest way to estimate the second largest eigenvalue of $`𝐀`$ is to compute its trace. Since $`𝐀`$ is positively defined, we have

$$\mathrm{Tr}𝐀^n1+\lambda _2(𝐀)^n.$$
(13)
Thus if we can find such $`n`$ that $`𝖤(\mathrm{Tr}𝐀^n)1`$, it will imply that $`𝖤(\lambda _2(𝐀))0`$. It turns out that it is enough to take $`n=2`$. However let us consider $`n=1`$ too.

$$\mathrm{Tr}𝐀=\underset{i=1}{\overset{N}{}}(𝐌^T𝐌)_{i,i}=\underset{i,j=1}{\overset{N}{}}(𝐌^T)_{i,j}(𝐌)_{j,i}=\underset{i,j=1}{\overset{N}{}}M_{j,i}^2=\underset{i,j=1}{\overset{N}{}}|U_{i,j}|^4.$$
(14)
And calculating the mean

$$𝖤(\mathrm{Tr}𝐀)=\underset{i,j=1}{\overset{N}{}}𝖤|U_{i,j}|^4=N^2𝖤|U_{1,1}|^4,$$
(15)
where due to the invariance of measure of CUE$`(N)`$ the means of the different matrix elements are equal. To calculate $`𝖤|U_{1,1}|^4`$ one can either integrate over the measure of CUE$`(N)`$, as was done in , or apply various invariance considerations . The result is $`𝖤|U_{1,1}|^4=2/(N^2+N)`$ and therefore $`𝖤(\mathrm{Tr}𝐀)2`$ as $`N\mathrm{}`$. This agrees with the numerical observation that $`𝖤|\lambda _2(𝐌)|1/\sqrt{N}`$: the eigenvalues of $`𝐀`$ are then of order $`1/N`$ and there are $`N1`$ of them (not counting the 1).

In the case $`n=2`$ we have

$`\mathrm{Tr}𝐀^2={\displaystyle \underset{i,j,k,l=1}{\overset{N}{}}}(𝐌^T)_{i,j}(𝐌)_{j,k}(𝐌^T)_{k,l}(𝐌)_{l,m}={\displaystyle \underset{i,j,k,l=1}{\overset{N}{}}}|U_{j,i}U_{j,k}U_{l,i}U_{l,k}|^2`$
$`={\displaystyle \underset{ik,jl}{}}|U_{j,i}U_{j,k}U_{l,i}U_{l,k}|^2+{\displaystyle \underset{i=k,jl}{}}|U_{j,i}U_{l,i}|^4`$
$`+{\displaystyle \underset{ik,j=l}{}}|U_{j,i}U_{j,k}|^4+{\displaystyle \underset{i=k,j=l}{}}|U_{j,i}|^8.`$ (16)
Applying the averaging we again find that due to the invariance of the measure all contributions from the first sum are the same and there are $`(N^2N)^2`$ of them; there are $`2N^2(N1)`$ contributions from the second and the third sum (the contributions are equal) and $`N^2`$ contributions from the last. Counting the number of terms in each sum, we write

$`𝖤\mathrm{Tr}𝐀^2=(N^2N)^2𝖤|U_{1,1}U_{1,2}U_{2,1}U_{2,2}|^2`$
$`+2N^2(N1)𝖤|U_{1,1}U_{1,2}|^4+N^2𝖤|U_{1,1}|^8,`$ (17)
where the averages can be calculated using ,

$`𝖤|U_{1,1}U_{1,2}U_{2,1}U_{2,2}|^2`$ $`=`$ $`{\displaystyle \frac{N^2+N+2}{N^2(N^21)(N+2)(N+3)}},`$ (18)
$`𝖤|U_{1,1}U_{1,2}|^4`$ $`=`$ $`{\displaystyle \frac{4}{N(N+1)(N+2)(N+3)}},`$ (19)
$`𝖤|U_{1,1}|^8`$ $`=`$ $`{\displaystyle \frac{24}{N(N+1)(N+2)(N+3)}}.`$ (20)
Bringing everything together, we obtain

$$𝖤\mathrm{Tr}𝐀^2=\frac{N^3+8N^2+17N2}{(N+1)(N+2)(N+3)}=1+\mathrm{O}\left(\frac{1}{N}\right),$$
(21)
which effectively finishes the proof. ∎

## 5 Shape of the distribution of $`|\lambda _2(𝐌)|`$

It is now clear that the distributions $`f_N(|\lambda _2|)`$ have singular limit. The natural question to ask, then, is what is the shape of $`f_N(|\lambda _2|)`$ for finite, but large, $`N`$ and whether there are sequences of normalizing coefficients $`a_N`$ and $`b_N`$ such that the random variables $`(|\lambda _2|a_N)/b_N`$ have continuous limiting distribution. While we are unable to give a definite answer, we can make a guess. Since $`|\lambda _2|=\mathrm{max}_{j=2,\mathrm{},N}|\lambda _j|`$, it is not unnatural to conjecture that for large $`N`$ the distribution of the second largest eigenvalue $`|\lambda _2|`$ is well approximated by the Generalized Extreme Value distribution (GEV) , specified by its cumulative distribution function

$$G(x)=\mathrm{exp}\left\{\left[1+\xi \left(\frac{xa}{b}\right)\right]^{1/\xi }\right\},$$
(22)
for some ($`N`$-dependent) $`\xi `$, $`a`$ and $`b`$: the GEV distribution is used to model maxima of sequences of random variables.

To check our conjecture we fit the numerical data for $`N=80`$ to the form (22) using an Splus routine by Coles . The routine minimizes the log-likelihood function

$$l(a,b,\xi )=\underset{i=1}{\overset{k}{}}\left\{\mathrm{ln}b\left[1+\xi \left(\frac{x_ia}{b}\right)\right]^{1+1/\xi }\left[1+\xi \left(\frac{xa}{b}\right)\right]^{1/\xi }\right\},$$
(23)
where $`x_i`$ is the sequence of observations of $`|\lambda _2|`$. The result of the routine is plotted as the probability density function $`G^{}(x)`$ for fitted $`\xi =0.07`$, $`a`$ and $`b`$ against the histogram of the numerical data for $`|\lambda _2|`$. One can see that the agreement is very good. The probability and quantile plots, which we do not present here, also show good agreement. As $`N`$ is taken to infinity, the values of $`a`$ and $`b`$ decay to zero, so that the limit of the distribution is the step function, and $`\xi `$ saturates at a nonzero value.

## 6 Conclusions

We have shown that if we take element-wise modulus squared of a large unitary matrix $`𝐔`$, the resulting doubly stochastic matrix $`𝐌`$ will have large spectral gap: with high probability all $`N1`$ “free” eigenvalues will have small modulus. In particular, it means that the Markov chain corresponding to $`𝐌`$ will quickly relax from any initial distribution into the uniform stationary state.

While our observation is by no means a proof of Tanner’s  conjecture that the families of graphs with large gap in their “classical” spectrum will adhere to RMT predictions, it is another powerful argument in its favour. We also note that the unitary evolution matrix of a general graph is sparse due to the topological restrictions . Thus, although arbitrary $`𝐔`$ can be considered as the scattering matrix for an $`N`$-star graph, the generalisation of the gap distribution to graphs with less trivial topology still remains to be studied.

The results for the mean of the traces of $`𝐌^T𝐌`$ and $`(𝐌^T𝐌)^2`$, obtained in the course of our proof, completely agree with the numerical observation that the mean gap increases with $`N`$ like $`1\mathrm{O}\left(N^{1/2}\right)`$. We also presented numerical evidence in favour of the claim that for finite $`N`$ the distribution of the second largest eigenvalue is well approximated by the generalized extreme value distribution.

The author is very grateful to M. Grinfeld, B. Gutkin, J.P. Keating, U. Smilansky, N. Snaith, M.G. Stepanov, S. Subba Rao Tata, G. Tanner and K. Życzkowski. Without their suggestions and encouragement this work would not be possible. The author was supported by the Israel Science Foundation, a Minerva grant, and the Minerva Center for Nonlinear Physics.

## References


text:
# Doping dependence of the electrical and thermal transport properties of CaB6

## I Introduction

Recent experiments on hexaborides with divalent metal cations of the alkaline-earth series, namely Ca, Sr and Ba, have revealed some unusual physical properties of these materials. In particular, an itinerant type of weak ferromagnetic order, stable up to temperatures between 600 and 900 K, has been observed in CaB6 and related alloys in a very narrow range of electron doping.  The electronic properties of alkaline-earth hexaborides place these materials close to a metal-insulator transition.  This has been confirmed theoretically by recent calculations of the electronic structure of divalent hexaborides which have shown that, with the exception of a small region around the X point of the cubic Brillouin zone, the valence- and the conduction bands are separated by a large gap of several eV.  The slight band overlap is reduced if the inter-octahedron B-B bond distances increase, finally resulting in the opening of a gap over the entire momentum space. It has also been predicted that the transport properties of SrB6 may be strongly dependent on doping, with changes from a $`p`$\- to an $`n`$-type metal induced by slight shifts of the Fermi energy around the zero doping value.  These observations have led to the speculation that the ground state of undoped divalent hexaborides may be characterized by a Bose condensate of bound electron-hole pairs, or excitons.  Weak ferromagnetism may then develop because of a spontaneous time-reversal symmetry breaking via doping. More recent calculations of the single-particle excitation spectrum of CaB6, however, have led to the claim that this material is not a semimetal but a semiconductor with a minimum band gap of 0.8 eV.  If this possibility is considered, the experimental observations, which suggest that binary hexaborides are close to a metal-insulator transition, would indicate the presence of a defect band of itinerant charge carriers. In an attempt to resolve the uncertainties with respect to the correct description of this material, we have studied the influence of the chemical composition on the physical properties of CaB6. The electrical resistivity, magnetoresistance, low-temperature specific heat and the optical conductivity of stoichiometric and doped CaB6 samples have been studied and the results have been presented in a previous publication.  Below, we present the results of measurements of the thermal conductivity and the thermopower. While data on the electrical conductivity and the electronic contribution to the thermal conductivity provide information on the density of electronic states $`N(E)`$ close to the Fermi energy $`E_\mathrm{F}`$, the thermopower may serve to establish the energy derivative of $`N(E)`$, i.e., an additional important detail of the electronic excitation spectrum. One of the samples was prepared in such a way as to obtain material with a close to stoichiometric composition, denoted as CaB6. The second sample, which we denote as Ca1-δB6, contained a small number of vacancies on the calcium sites, inadvertently introduced during the flux-growth procedure and leading to a certain degree of self doping which is difficult to control. This paper is organized as follows. After a brief description of the sample preparation and the experimental methods used in this investigation in section II, we present, in section III, the results of our measurements and their analysis. The conclusions are presented in section IV.

## II Samples and experimental methods

Binary hexaborides can be synthesized in a narrow range of composition, with the tendency to be boron-rich with concentrations of metal vacancies up to several percent. Stoichiometric CaB6 crystals can be obtained close to the border of the metal-rich phase boundary. For our experiments the single crystals were grown by a slow-cooling procedure in aluminum flux  starting with a nominal ratio CaB3 for the stoichiometric CaB6 sample, and with CaB12 for the metal-deficient Ca1-δB6 sample. The crystals were removed from the flux by leaching in a concentrated sodium hydroxide solution. Subsequent etching with HNO3 was intended to remove possible surface contaminations. The samples on which our transport measurements were made were of prism-type shape with approximate overall dimensions of $`4.2\times 0.5\times 0.45\mathrm{mm}^3`$ for CaB6 and $`4.7\times 0.45\times 0.4\mathrm{mm}^3`$ for Ca1-δB6.

The thermoelectric power $`S`$ and the thermal conductivity $`\kappa `$ of both samples were measured simultaneously by means of a standard steady-state heat-flow technique. A commercial 4He gas-flow cryostat was used for cooling the sample holder. At one end of a prism-shaped sample, the thermal contact to a copper heat sink was achieved by using high-conductance silver epoxy. The sample heater, consisting of a $`100\mathrm{\Omega }`$ ruthenium-oxide chip resistor, was attached to the other end of the prism by the same method. Joule heating caused by heater currents of the order of a few mA provided the necessary heat-flow and hence a thermal gradient along the crystals. In order to measure both the temperature difference $`\mathrm{\Delta }T`$ between two contacts mounted perpendicularly to the heat flow as well as the thermoelectric voltage of our samples with respect to chromel, we used two pairs of calibrated 0.025 mm Au-Fe (0.07 at. %) versus chromel thermocouples.  In order to minimize spurious thermal voltages, the thermocouple leads were connected to uninterrupted copper wires reaching three home-built low-noise voltage amplifiers, mounted directly on top of the cryostat insert.

## III Experimental results and Analysis

### A Thermoelectric power

The thermoelectric power $`S(T)`$ of both samples is shown in Fig. 1 on linear scales. The negative sign of $`S`$ clearly demonstrates the dominant $`n`$-type character of the investigated crystals. The rather large values of $`S`$, if compared to the typical values observed for the thermoelectric power of common metals, are an indication for the low itinerant charge-carrier concentration in these materials, compatible with the very low values of the electrical conductivity reported previously.  The overall temperature dependence of $`S`$ is obviously non-linear, with features that are reminiscent of band-structure effects affecting the thermopower of elements and crystalline alloys.  Nevertheless, below $`T^{}=20`$ and 40 K for CaB6 and Ca1-δB6, respectively, the curves show an approximately linear temperature dependence. A negative thermoelectric power varying linearly with $`T`$ is usually identified as the diffusion thermopower $`S`$ of metals, where the Fermi-energy $`E_\mathrm{F}`$ is located within the conduction band. In the free-electron approximation,

$$S=\frac{\pi ^2k_\mathrm{B}^2T}{3|e|E_\mathrm{F}/s},$$
(1)
with $`E_\mathrm{F}`$ as the Fermi-energy measured from the bottom of the conduction band, and $`s`$ a factor, typically of the order of unity, describing the energy dependence of the scattering time $`\tau E^{s3/2}`$ for the scattering mechanism dominating the low temperature behaviour.  From the average slope of $`S(T)`$ below $`T^{}`$, we may thus calculate a rescaled Fermi energy $`E_\mathrm{F}/s`$. These values are $`E_\mathrm{F}/s=4.5`$ and 8.7 meV for CaB6 and Ca1-δB6, respectively. Considering the entire temperature dependence of $`S`$, in particular the abrupt changes of slopes of $`S(T)`$ around $`T^{}`$, suggests that ultimately, $`S(T)`$ cannot be explained by simply taking into account electronic states in the conduction band alone. In the following, we try to identify the possible reasons for the measured departure from linearity in the temperature dependence of $`S`$.

With the common assumption that the current densities $`𝐣`$ and $`𝐰`$ for electrical and thermal transport, respectively, respond linearly to their driving forces, we can write

$`\begin{array}{cc}\hfill 𝐣& =e^2L_0𝐄{\displaystyle \frac{e}{T}}L_1_\mathrm{r}𝐓,\hfill \\ \hfill 𝐰& =eL_1𝐄{\displaystyle \frac{1}{T}}L_2_\mathrm{r}𝐓,\hfill \end{array}`$ (2)
where $`𝐄`$ is the external electric field and $`_\mathrm{r}𝐓`$ is the temperature gradient. In Eq. 2, the transport coefficients $`L_\mathrm{i}`$ (i=0,1,2) explicitly fulfill the Onsager relations.  For $`𝐣=0`$, we obtain $`𝐰=\kappa _{\mathrm{el}}_\mathrm{r}T`$ and $`𝐄=S_\mathrm{r}T`$, whereas for $`_\mathrm{r}T=0`$, we get $`𝐣=\sigma 𝐄`$. These are the equations defining the electronic contribution to the thermal conductivity $`\kappa _{\mathrm{el}}`$, the thermoelectric power $`S`$, and the electrical conductivity $`\sigma `$, such that

$`\kappa _{\mathrm{el}}`$ $`={\displaystyle \frac{1}{T}}\left(L_2{\displaystyle \frac{L_1^2}{L_0}}\right),`$ (3a)
$`S`$ $`={\displaystyle \frac{1}{|e|T}}{\displaystyle \frac{L_1}{L_0}},`$ (3b)
$`\sigma `$ $`=e^2L_0.`$ (3c)
The integrals $`L_\mathrm{i}`$ (i=0,1,2), are defined as

$$L_\mathrm{i}=_{\mathrm{}}^+\mathrm{}\stackrel{~}{\sigma }(E)(EE_\mathrm{F})^i\frac{f_0}{E}𝑑E,$$
(4)
with $`\stackrel{~}{\sigma }(E)`$ representing the conductivity spectrum, including all system-dependent features. By partial integration of Eq. 4 and by using the Sommerfeld expansion scheme, we get

$$S=\frac{\pi ^2k_\mathrm{B}^2T}{3|e|}\left(\frac{\mathrm{ln}\stackrel{~}{\sigma }(E)}{E}\right)_{E=E_\mathrm{F}},$$
(5)
which approximates Eq. 3b very well at temperatures $`k_\mathrm{B}TE_\mathrm{F}`$. Since we intend to interpret our data over the entire covered temperature range, we will not use the low-temperature approximation represented by Eq. 5, in our analysis. For our purposes the general form of $`S`$ shown in Eq. 3b is more suitable.

Since $`\stackrel{~}{\sigma }(E)>0`$, the coefficient $`L_0`$ is strictly positive and the sign of $`S`$ is determined by the sign of the integral $`L_1`$. By inspecting the integrand of $`L_1`$, represented in Eq. 4, and also by considering Eq. 3b, it may be concluded that the states with energies higher than the chemical potential, i.e., $`E>E_\mathrm{F}`$, provide a negative contribution to the thermoelectric power, whereas those states which are located below $`E_\mathrm{F}`$ contribute with the opposite sign. At low temperatures, where only the states very close to $`E_\mathrm{F}`$ contribute to the electronic transport, the measured $`S`$ values are negative. This leads us to conclude that $`\stackrel{~}{\sigma }(E)`$ must increase monotonously across $`E_\mathrm{F}`$, so that the states contributing negatively to $`S(T)`$ in $`L_1`$ acquire a larger weight than those located below the Fermi-energy, which contribute positively to $`S(T)`$. In order to reproduce the observed reduction of the slope of $`S(T)`$ at $`T^{}`$ in the calculation, however, we are forced to introduce a “feature” in $`\stackrel{~}{\sigma }(E)`$ centered at an energy $`E^{}`$ such that $`|E^{}E_\mathrm{F}|k_\mathrm{B}T^{}`$. Since the absolute value of $`S/T`$ decreases above $`T^{}`$, the feature in $`\stackrel{~}{\sigma }(E)`$ must either reduce the negative contribution to $`S(T)`$ given by the states above $`E_\mathrm{F}`$, or else enhance the positive contribution to $`S(T)`$ of the states below $`E_\mathrm{F}`$. In what follows, we assume that in our case the second possibility is realized and we add an additional part to the $`\stackrel{~}{\sigma }(E)`$ spectrum, decreasing linearly with increasing energy below $`E_\mathrm{F}`$. The complete chosen $`\stackrel{~}{\sigma }(E)`$ spectrum, corresponding to the Eqs. 7 and 10 below, is shown in Fig. 2.

Based on this preliminary analysis and the model outlined above, we have calculated the electrical transport properties of CaB6 and Ca1-δB6. The Fermi energy is fixed within a parabolic conduction band. Close to the bottom of this band we position an additional band, which may be interpreted as a defect band. With these assumptions, the conductivity spectrum may be written as

$$\stackrel{~}{\sigma }(E)=\stackrel{~}{\sigma }^n(E)+\stackrel{~}{\sigma }^p(E).$$
(6)
For an isotropic parabolic conduction band, the Boltzmann equation in the relaxation time approximation leads to

$$\stackrel{~}{\sigma }^n(E)=\{\begin{array}{cc}0\hfill & (E<E_\mathrm{n})\hfill \\ \frac{1}{3\pi ^2m_\mathrm{n}}\tau (E)\left(\frac{2m_\mathrm{n}}{\mathrm{}^2}\right)^{3/2}(EE_\mathrm{n})^{3/2}\hfill & (E>E_\mathrm{n}),\hfill \end{array}$$
(7)
where the bottom of the conduction band has been fixed to $`E_\mathrm{n}=0\mathrm{meV}`$. By setting $`\stackrel{~}{\sigma }(E)=\stackrel{~}{\sigma }^n(E)`$ in Eq. 5, it is easy to recover Eq. 1. For the effective mass of the conduction electrons we used $`m_n=0.28m_0`$ with $`m_0`$ as the free electron mass. We only considered two scattering mechanisms of the charge carriers in the conduction band, i.e., scattering by acoustic lattice vibrations, approximated by a rate

$$\tau _{\mathrm{lv}}^1=\alpha E^{1/2}T,$$
(8)
and scattering by ionized impurities, with a rate

$$\tau _{\mathrm{ii}}^1=\beta E^{3/2},$$
(9)
where $`\alpha `$ and $`\beta `$ are two constants to be determined by the fitting procedure. If the defect states which are forming a band below the bottom of the conduction band coexist with localized states, the latter may act, when ionized, as scattering centers for the electrons. This justifies considering a term $`\tau _{\mathrm{ii}}^1`$ in the scattering rate of the electrons. As we will show below, the existence of localized donor states may also be compatible with our thermal conductivity data. Another possibility is to attribute the rate $`\tau _{\mathrm{ii}}^1`$ to the vacancies at the Ca2+-ion sites which would act as negatively charged scattering centers. With the common assumption that the two scattering mechanisms do not interfere with each other, we may use Matthiessen’s rule in the form $`\tau =(\tau _{\mathrm{lv}}^1+\tau _{\mathrm{ii}}^1)^1`$ for the electrons’ average scattering time.

For the defect band, we postulate

$$\stackrel{~}{\sigma }^p(E)=\{\begin{array}{cc}\gamma (E_{\mathrm{me}}E)\hfill & (E<E_{\mathrm{me}})\hfill \\ 0\hfill & (E>E_{\mathrm{me}}),\hfill \end{array}$$
(10)
with a mobility edge $`E_{\mathrm{me}}`$ and $`\gamma `$ a constant. Since the quantities $`L_0`$ and $`L_1`$ which dictate the temperature dependence of the thermopower, appear in the ratio $`L_1/L_0`$, $`S`$ depends only on the relative magnitude of the constants $`\alpha `$, $`\beta `$, and $`\gamma `$. This allows some reduction of the number of fitting parameters. In Fig. 1, the solid curves represent the results of the calculation of $`S(T)`$. By considering the simplicity of the model and the small number of free parameters ($`\alpha /\gamma ,\beta /\gamma ,E_\mathrm{F},\text{and }E_{\mathrm{me}}`$), the agreement with the experimental data is quite remarkable. The Fermi energies which emerge from the fitting procedure amount to 14 and 28 meV for CaB6 and Ca1-δB6, respectively, corresponding to 0 K charge-carrier densities $`n_\mathrm{c}`$ of $`1.1\times 10^{24}\mathrm{m}^3`$ and $`3.2\times 10^{24}\mathrm{m}^3`$, or to $`0.8\times 10^4`$ and $`2.3\times 10^4`$ charge carriers per unit cell, respectively. Preliminary results of low-temperature Hall effect measurements performed on a CaB6 sample with the same resistivity as our vacancy-doped sample lead, when interpreted with a single-band model, to a charge carrier density of $`4.2\times 10^{24}\mathrm{m}^3`$ (Ref. ), hence compatible with our results. The value of the mobility edge $`E_{\mathrm{me}}`$, introduced to reproduce the kink in $`S(T)`$ at the observed temperatures is approximately $`+8`$ meV for both samples. We note, however, that this value strongly depends on the choice of the energy exponent of the low-temperature scattering term. In this sense we have to admit that this result depends on the employed model.

Our calculation of the thermoelectric power relies essentially on the model of the conductivity spectrum (modulo a factor), which we may now test by calculating other electronic transport properties such as the electrical resistivity $`\rho (T)`$ and the thermal conductivity $`\kappa _{\mathrm{el}}(T)`$ and compare the results with the experimental data. In Fig. 3, we have plotted the rescaled electrical resistivity of CaB6 calculated with Eq. 3c, together with the experimental data taken on two samples of the same batch. The agreement between the calculated curve and the experimental data is acceptable, especially if one considers primarily the temperature dependence.

By using Eqs. 3c, 4, and 7, combined with the measured resistivity curves shown in Fig. 3, we conclude that the mean free path $`l_\mathrm{F}=\mathrm{}k_\mathrm{F}\tau (E_\mathrm{F})/m_\mathrm{n}`$ of the conduction electrons in CaB6 is less than 100 Å. The same length is obtained, obviously enough, by employing the Drude formula $`\sigma =ne^2\tau /m_n`$ and by using $`n=n_\mathrm{c}`$, the value cited above for the charge carrier concentration of CaB6. Since the electrons’ wave vector at the Fermi-energy is of the order of $`10^2`$ Å-1, it seems fair to assume that the Ioffe’s criterion for localization ($`k_\mathrm{F}l_\mathrm{F}<1`$) is close to be fulfilled in CaB6.

### B Thermal transport

In Fig. 4, we show the measured thermal conductivities $`\kappa _{\mathrm{meas}}`$ of CaB6 and Ca1-δB6 on logarithmic scales between 6 and 300 K. In the further analysis we assume that only conduction-electron and lattice contributions to $`\kappa _{\mathrm{meas}}`$ need to be considered.

A measure for the electronic contribution to the thermal conductivity is given by the quantity $`L^{}=\kappa _{\mathrm{el}}/\sigma T`$, which can be calculated from the Eqs. 3a and 3c. Just like the thermoelectric power $`S`$, $`L^{}`$ depends only on the relative magnitude of the constants $`\alpha `$, $`\beta `$, and $`\gamma `$. In a degenerate one-band model, the value

$$\kappa _{\mathrm{el}}/\sigma T=L_0=\frac{\pi ^2}{3}\left(\frac{k_\mathrm{B}}{e}\right)^2$$
(11)
is denoted as the Lorenz number. The ratio $`L^{}`$ can be much higher in semiconductors than in metals, due to the so-called ambipolar diffusion, i.e., the diffusion of electron-hole pairs. In Fig. 5 we plotted the temperature dependence of $`L^{}`$, calculated by using the fitting parameters provided by our analysis of $`S(T)`$, together with the value $`L_0=2.45\times 10^8(\mathrm{V}/\mathrm{K})^2`$. As expected, $`\kappa _{\mathrm{el}}/\sigma T`$ is equal to the Lorenz number at very low temperatures. This is so because only the electrons of the conduction band contribute to the transport. The influence of the ambipolar diffusion to the electronic thermal conductivities of CaB6 and Ca1-δB6 is manifest at intermediate temperatures. Nevertheless, this effect is still quite small if compared to enhancements of $`L^{}`$ by two orders of magnitude observed in intrinsic semiconductors.

The electronic contributions to the thermal conductivities of CaB6 and Ca1-δB6, shown in Fig. 4 as solid lines, have beeen calculated by using the Wiedemann-Franz relation $`\kappa _{\mathrm{el}}=L^{}\sigma T`$, with $`L^{}`$ calculated from the Eqs. 3a and 3c as shown in Fig. 5. For the electrical conductivity $`\sigma `$, we used the values which have typically been measured on close-to-stoichiometric and vacancy-doped samples, respectively.  For both hexaborides, the estimated $`\kappa _{\mathrm{el}}`$ is at least two orders of magnitude smaller than $`\kappa _{\mathrm{meas}}`$ at all temperatures covered in this study. Since the experimental uncertainty in our measurements of the thermal conductance is of the order of $`0.5\%`$, the procedure for estimating $`\kappa _{\mathrm{el}}`$ is not critical and may be considered as sufficiently accurate for the validity of the following considerations.

The lattice contributions $`\kappa _{\mathrm{ph}}=\kappa _{\mathrm{meas}}\kappa _{\mathrm{el}}`$ to the thermal conductivities of both samples are shown in Fig. 6, revealing monotonic increases with $`T`$ up to 18 K and 23 K for CaB6 and Ca1-δB6, respectively. At these temperatures, $`\kappa _{\mathrm{ph}}`$ passes over maxima of 124 and 74 $`\mathrm{Wm}^1\mathrm{K}^1`$, respectively, and decreases with varying slopes up to room temperature. We note that at low temperatures $`\kappa _{\mathrm{ph}}`$ of Ca1-δB6 attains lower values than $`\kappa _{\mathrm{ph}}`$ for CaB6. This seems reasonable because of an expected enhanced scattering of the phonons by defects in the Ca-deficient material. The crossing of the curves at a temperature close to 30 K is a rather surprising experimental observation which indicates a more efficient high-frequency scattering mechanism of the phonons in the structurally more perfect sample.

The following quantitative analysis of our data is based on a relaxation-time approximation in the context of the Debye model. The Debye approximation is valid at temperatures which are low enough with respect to the Debye temperature $`\mathrm{\Theta }_\mathrm{D}`$, where we may assume that the acoustic branches $`\omega (k)`$ of the vibrational modes are well approximated by a linear $`k`$ dependence. Since the temperature range covered in our experiment is considerably lower than the Debye temperature $`\mathrm{\Theta }_\mathrm{D}=783`$ K,  this approximation is justified to some extent. With this simplification, the thermal conductivity may be calculated from

$$\kappa _{\mathrm{ph}}^{\mathrm{fit}}=3nk_\mathrm{B}\left(\frac{T}{\mathrm{\Theta }_\mathrm{D}}\right)^3v_{\mathrm{ph}}^2_0^{\mathrm{\Theta }_\mathrm{D}/T}\frac{x^4e^x}{(e^x1)^2}\tau (x,T)𝑑x,$$
(12)
where $`\tau (x,T)`$ is the mean time between collisions of a phonon, $`n`$ the number density of atoms in the crystal, and $`x=\mathrm{}\omega /k_\mathrm{B}T`$. By limiting the total density of allowed degrees of freedom to $`3n`$, it is possible to calculate the cut-off frequency $`\omega _\mathrm{D}`$, which is related to the slope of $`\omega (k)`$. This cut-off frequency defines the Debye temperature $`\mathrm{\Theta }_\mathrm{D}`$ of the material which is, in this sense, a measure of the phonon velocity and is usually extracted from low-temperature specific-heat or elastic-constant measurements. For the Debye temperature $`\mathrm{\Theta }_\mathrm{D}`$, we used the above cited value which was extracted from specific-heat measurements on stoichiometric CaB6 Since upon introducing a very small number of calcium vacancies, no significant changes in the lattice contribution to the specific heat are expected, we used the same Debye temperature for CaB6 and Ca1-δB6. This argument is supported by the observation that there are no significant differences between the low-temperature specific heats of CaB6 and Ca1-xLaxB6 with x=0.005.  The number density of atoms in CaB6 can be calculated by using the lattice constant $`a=4.146`$ Åreported in Ref. . From the relation

$$n=1/6\pi ^2(k_\mathrm{B}\mathrm{\Theta }_\mathrm{D}/\mathrm{}v_{\mathrm{ph}})^3,$$
(13)
we may also calculate the average velocity of the vibrational modes for which we obtain $`v_{\mathrm{ph}}=5722\mathrm{m}/\mathrm{sec}`$. This value is compatible with the average sound velocity $`v_{\mathrm{ph}}=4200\mathrm{m}/\mathrm{sec}`$, valid for EuB6  when scaled with the mass-density factor $`\sqrt{\rho _{\mathrm{EuB}_6}/\rho _{\mathrm{CaB}_6}}`$.

For evaluating the integral in Eq. 12, we approximate the phonon relaxation rate by

$$\tau ^1(x,T)=(\tau _{\mathrm{Cas}}^1+\tau _{\mathrm{disl}}^1+\tau _{\mathrm{Rayl}}^1+\tau _{\mathrm{res}}^1)(x,T).$$
(14)
The different terms on the right hand side of Eq. 14 represent the scattering of phonons at grain boundaries, by the strain fields surrounding dislocations, by point defects and via resonant scattering, respectively. The resulting fits are shown in Fig. 6 as solid lines. In what follows, we discuss each term appearing in Eq. 14 separately. The values of the free parameters involved in our fitting procedure (see below) are given in Table I and the resulting total scattering rates of the phonons are shown as a function of frequency in Fig. 7 for $`T=10`$ and 100 K.

The temperature and frequency independent scattering length $`l_{\mathrm{Cas}}`$ in the first term in Eq. 14, $`\tau _{\mathrm{Cas}}^1=v_{\mathrm{ph}}/l_{\mathrm{Cas}}`$, is the so-called Casimir lenght, where the phonon mean free path is limited by the sample dimensions.  The values of 0.303 mm and 0.083 mm for CaB6 and Ca1-δB6, respectively, are of the order of magnitude of the smallest sample dimensions, resulting from the variation of the sample cross section perpendicular to the direction of the thermal current.

The second term in Eq. 14, $`\tau _{\mathrm{disl}}^1=A\omega `$, is usually attributed to phonon-electron scattering but, as pointed out in Ref. , can also be caused by phonon scattering due to strain fields surrounding dislocations. As may be seen in Table I, this term is slightly larger for CaB6 than for Ca1-δB6. In view of the higher density of charge carriers in Ca1-δB6 derived in section III A, this rules out that mainly phonon-electron scattering is responsible for this term. Phonon-electron scattering may also be excluded by inspecting the relevant scattering rate which can be written as

$$\tau _{\mathrm{cond}}^1=\frac{m_n^2C_\mathrm{D}^2}{2\pi \rho _{\mathrm{CaB}_6}v_{\mathrm{ph}}\mathrm{}^3}\omega ,$$
(15)
with $`C_\mathrm{D}\frac{2}{3}E_\mathrm{F}`$ as the deformation potential. The fit values for $`A`$ given in Table I result in Fermi-energies of the order of 1-2 eV, i.e., values which exceed $`E_\mathrm{F}`$ evaluated in section III A by two orders of magnitude. For phonon scattering at dislocations we have

$$\tau _{\mathrm{disl}}^1=0.06n_\mathrm{d}(b\gamma _\mathrm{G})^2\omega ,$$
(16)
with $`n_\mathrm{d}`$ as the dislocation density, $`\gamma _\mathrm{G}`$ the Grüneisen parameter, and b the Burger vector. By using the lattice constant for $`b4`$ Åand $`\gamma _\mathrm{G}1`$, we obtain $`n_\mathrm{d}2\times 10^{15}\mathrm{m}^2`$, a plausible result.

The enhanced presence of scattering centers at the vacancy sites in Ca1-δB6 is unequivocally reflected in the third term of Eq. 14. Rayleigh scattering by point-like mass defects leads, in the continuum approximation, to a scattering rate term $`\tau _{\mathrm{Rayl}}^1=B\omega ^4`$, with

$$B=\frac{n_\mathrm{p}a^6}{4\pi v_{\mathrm{ph}}^3}\left(\frac{\delta M}{M}\right)^2.$$
(17)
In Eq. 17, $`n_\mathrm{p}`$ is the number density of point defects in Ca1-δB6, $`a=4.146`$ Åis its lattice constant, $`M`$ is the mass of a unit cell, and $`\delta M`$ is the supplementary or missing mass due to a point defect in a sphere with the radius of the point defect itself. Since the extra mass term appears squared in Eq. 17, it is not possible to decide whether $`\delta M`$ is positive or negative. By assuming that each point defect in Ca1-δB6 provides one free electron to the ensemble of conduction electrons and that for stoichiometric CaB6$`n_\mathrm{c}=0`$, we have $`n_\mathrm{p}=3.2\times 10^{24}\mathrm{m}^3`$, i.e., the same value as the electron density calculated in section III A. Again, by using for $`B`$ the value which we obtained from the fitting procedure, we note that $`\delta M/M`$ is of order unity. In view of the rather crude approximations, this result is encouraging, as it seems to indicate that the concentration of itinerant electrons in CaB6 is equivalent to the concentration of impurities (or vacancies) with approximately the same mass as the mass of the unit cells. The scattering rate term $`\tau _{\mathrm{Rayl}}^1`$ is only of significance for fitting the $`\kappa _{\mathrm{ph}}`$ data of the vacancy-doped Ca1-δB6 sample. Our attempts to evaluate the free parameter $`B`$ for CaB6 resulted in a vanishing contribution of $`\tau _{\mathrm{Rayl}}^1`$ to the total scattering rate of the stoichiometric hexaboride. Considering the estimated error range of a few times $`10^{44}`$ sec3 for $`B`$, this is at least consistent with a reduced amount of defects in CaB6.

The last term appearing in Eq. 14 needs a detailed examination. Resonant scattering of phonons has already been invoked in the analysis of the thermal conductivity data of EuB6 in a previous investigation.  The appearance of $`\tau _{\mathrm{res}}^1`$ in the scattering rate of EuB6 was intepreted as being due to a resonant scattering of low-momentum vibrational states in the dispersionless part of the $`\omega (k)`$ branches with acoustic phonons of constant energy that seem to occupy an extended region of the Brillouin zone at higher momenta. Indeed, a rapid flattening of the LA branches at an energy of approximately 12 meV has been observed in inelastic neutron scattering experiments performed on the trivalent hexaborides XB6 (X = Ce, La, Sm).  Because of similar atypical temperature dependencies of $`\kappa _{\mathrm{ph}}`$ above 50 K, we follow the arguments given in Ref.  and introduce a resonant scattering term of the form

$$\tau _{\mathrm{res}}^1=C\frac{\omega ^4}{(\omega _{\mathrm{res}}^2\omega ^2)^2}F(T),$$
(18)
with

$$F(T)=\frac{1}{\mathrm{exp}(\mathrm{}\omega _{\mathrm{res}}/k_\mathrm{B}T)1}$$
(19)
as the bosonic population factor. In our intepretation of the data, the crossing of the thermal conductivity curves at $`T=30`$ K has to be traced back to a much stronger resonant scattering rate in CaB6. As may be seen in Table I, the fit value $`C`$ is more than one order of magnitude larger for CaB6 than for Ca1-δB6.

Although, by taking into account the resonant-scattering term, we achieve a fairly accurate description of our thermal conductivity data, it seems rather unlikely that the large difference in the temperature dependence of $`\kappa _{\mathrm{ph}}`$ may be accounted for by considering phonon-phonon scattering processes alone, especially if we consider the only small difference in chemical composition between the two materials.

In order to find a more plausible explanation for the experimentally observed crossing of $`\kappa _{\mathrm{ph}}(T)`$, we note that a strong reduction of the lattice thermal conductivity of semiconductors induced by defect donor states has been reported in a large number of publications.  The intensity of the phonon scattering processes invoking such states depends strongly on the electronic excitation spectra induced by these defects. In particular, donor states in cubic semiconductors, leading to the existence of localized singlet-triplet centers just below the conduction band, are also expected to provide the source for a strong resonant phonon scattering.  Physical realizations of such systems have been found in the valley-orbit-split ground state of the donors As, Sb, P, and Li in Ge. The singlet-triplet valley-orbit splitting of the fourfold degenerate ground state of As in Ge, for example, is 4.23 meV, a value which is of the same order of magnitude as the resonant energies of 10.85 and 9.38 meV, quoted in Table I for CaB6 and Ca1-δB6, respectively. If the density of defects in the crystal is such that the average distance between the donor states is of the order of the extension of the defect wave function, characterized by a Bohr radius of typically 40 Å, we expect a broadening of the degenerate energy-eigenvalues, finally resulting in the formation of a band of defect states. The typical extension of the defect wave function should be compared with the expected mean distance between the electrons in CaB6, estimated by calculating the third root of the inverse of the electron density derived in section III A. This leads to values of 70 Å and 100 Å for Ca1-δB6 and CaB6, respectively. If the electron densities at 0 K derived in section III A are due to the existence of donor states in CaB6, we are thus close to concentrations for which a band is expected to form. A schematic representation of this is shown in Fig. 8. Because of the statistical distribution of distances between donor states, it seems possible to realize a situation where both, localized donor states and defect-band itinerant states, coexist. Since stoichiometric CaB6 is closer to perfection than the vacancy-doped sample, the observation of a stronger high-frequency scattering mechanism in CaB6, as mentioned above, might appear as rather surprising. However, as we noted above, an increasing distance between the defects due to a lower concentration of defects would narrow the tight-binding bands in CaB6, eventually resulting in a larger fraction of localized donor states acting as resonant scattering centers for the phonons. Hence, if this scheme is valid, the higher resonant scattering rate in CaB6 makes sense.

All our attempts to cast the scattering of phonons at localized donor states into a reasonable equation for its relaxation rate and to obtain a fair description of the high temperature features of $`\kappa _{\mathrm{ph}}`$ for CaB6 failed. This failure may either be due to a breakdown of the Debye model at higher frequencies, or may indicate that localized donor states in CaB6 do not scatter phonons. With respect to the first possibility it may well be that using the Debye model at the frequencies around the resonance is not justified. A strong resonant scattering mechanism may, in some cases (see, for example, Ref. ), modify the unperturbed dispersion relations by opening a gap in the phonon spectrum at the resonant frequency, thus turning the Debye scheme into a very poor approximation at the frequencies around the resonance. This may also be the reason for the poor quality of the fit to the CaB6 data, for which the resonant-scattering is quite significant. Since the dispersion relations of the lattice vibrational modes in CaB6 have not been measured yet, we cannot offer a more detailed discussion. In case of deviations of the acoustic branches of CaB6 from a non-dispersive $`k`$ dependence at the frequencies which are relevant for our investigation, the functional form of $`\tau _{\mathrm{res}}^1`$ would certainly change.

Concluding this section, we recall that our thermal conductivity data show a strong and somewhat unexpected dependence on calcium vacancy-doping at elevated temperatures, which we interpret as being due to a resonant scattering of the vibrational modes. For a more realistic computation of the thermal transport properties at the high end of the temperature regime covered in our investigation, more information concerning the phonon dispersion relations of CaB6 is necessary, however.

## IV Summary and conclusion

We have measured the thermoelectric power and the thermal conductivity of close to stoichiometric and vacancy-doped CaB6 between 5 and 300 K. The high negative values of the thermoelectric power $`S(T)`$ of both materials indicate a low concentration of itinerant $`n`$-type charge-carriers. We have achieved a reasonable interpretation of $`S(T)`$ by using a relaxation-time approximation of Boltzmann’s equation, by considering a sizeable gap between valence and conduction band states, and by postulating the existence of an additional term in the conductivity spectrum which may be related to the existence of a defect band in proximity to the lower conduction band edge. A Debye-type relaxation-time approximation is successful in describing the thermal conductivity data across the entire temperature range covered in the experiments. The usual terms contributing to the total scattering rate of the vibrational states agree well with our expectations. A rather unusual temperature dependence of $`\kappa _{\mathrm{ph}}`$ is observed above approximately 30 K. A resonant type of scattering seems to influence the mean free path of the phonons at these temperatures. Quite surprisingly, its limiting effect on the time period between two vibrational-mode collisions is much stronger in the material closer to stoichiometry. Although postulating a resonant-type of phonon-phonon scattering leads to a satisfactory description of our experimental observation, an alternative possible interpretation which considers the existence of localized donor states is suggested. Since the lattice dispersion relations of CaB6 have not yet been investigated experimentally, we cannot offer any rigorous statements about the physical nature of this resonant type of scattering of the vibrational modes.

## V Acknowledgements

We thank R. Monnier and M. E. Zhitomirsky for stimulating discussions. This work was financially supported by the Schweizerische Nationalfonds zur Förderung der Wissenschaftlichen Forschung.