text:
# Monte Carlo reconstruction of the inflationary potential

## I Introduction

Given a detailed model of physics at the Grand Unified Theory (GUT) scale that supports inflation one can deduce many of the overall properties of the universe from first principles. However, new strides in precision cosmology are likely to present us with the inverse problem: given an accurate measurement of the primordial perturbation spectrum, can we determine the particle physics that drove inflation? This question leads to “reconstruction” – rebuilding a portion of the inflationary potential from observations – whose apotheosis is the work of Lidsey et al. lidsey95 . The conclusion of this effort was that while a limited reconstruction may be possible, it is overly optimistic to hope to deduce more than the rough form of a small piece of the inflaton potential on the basis of astrophysical observations alone.

Despite this apparent setback, the prospect of being able to probe ultra-high energy particle physics via astrophysical measurements remains tempting. In this paper we present a new approach to the task, Monte Carlo reconstruction. Unlike previous workers, we tackle the problem stochastically, by generating large numbers of inflationary models and their associated potentials and identifying those that lead to a spectrum of primordial perturbations which fits inside a specified window of parameter space. We sacrifice the prospect of unambiguously reconstructing the inflationary potential, but instead map out the ensemble of viable potentials. We can then assess the extent to which these potentials resemble one another, and thus the degree of ambiguity in the “reconstruction”.

The starting point for Monte Carlo reconstruction is the Hubble slow roll expansion liddle94 , which leads to a set of first order “flow” equations hoffman00 ; kinney02  that describes the inflationary dynamics. The Hubble slow roll expansion can be systematically continued to arbitrary order but, as we explain below, the truncated expansion actually yields an exact solution of the full equations, drawn from a subset of the overall solution space. We select a “universe” by assigning random values to each term in the truncated slow roll expansion, and integrate the flow equations forward to the end of inflation (or far enough into the future to show that the solution approaches an inflationary attractor). We can then determine the epoch that corresponds to cosmological structure formation and thus the predicted parameters of the primordial spectrum, allowing us to select the models which fit into a pre-determined window of parameter space.

As with any other reconstruction program, we face limitations dictated by the physics of inflation itself. Since the field rolls slowly, it traverses only a small portion of the potential during the epoch in which the primordial perturbations are laid down. Consequently, even if all the potentials we recover for a specific set of cosmological parameters overlap in the region where the cosmologically relevant portion of the spectrum is laid down, they will still diverge outside this region. However, because we reconstruct the entire potential, we can estimate the degree to which the reconstructed potentials diverge from one another.

In the following section we summarize the inflationary dynamics and the flow equations, and in Section 3 apply these to spell out the process of Monte Carlo reconstruction in detail. We present the results of a variety of simulations in Section 4, and test how accurately Monte Carlo reconstruction can recover a known inflationary potential. As expected, a fully consistent reconstruction requires the unambiguous detection of tensor modes in the CMB, but even in the absence of a detectable tensor signal we can still recover some information about the inflationary potential. If we assume the error bars on the scalar spectral index and tensor-to-scalar ratio that are expected for the Planck mission, the underlying potential cannot be unambiguously recovered. However, if we decrease these error bars (perhaps by adding information from large scale structure surveys, or better measurements of CMB polarization) by a factor of 5, then we can start to recover a $`\varphi ^n`$ potential. Our conclusions are summarized in Section 5.

## II Inflation and flow

We consider inflation driven by a single homogeneous scalar field $`\varphi `$ (the inflaton) with potential $`V(\varphi )`$ and equation of motion

$$\ddot{\varphi }+3H\dot{\varphi }+V^{}(\varphi )=0,$$
(1)
where $`H(\dot{a}/a)`$ is the Hubble parameter. We assume a spatially flat, homogeneous and isotropic universe where the inflaton field is the only contribution to the energy-momentum tensor, so that the Einstein field equations have the familiar form

$$H^2=\left(\frac{\dot{a}}{a}\right)^2=\frac{8\pi }{3m_{\mathrm{Pl}}^2}\left[V(\varphi )+\frac{1}{2}\dot{\varphi }^2\right],$$
(2)
and

$$\frac{\ddot{a}}{a}=\frac{8\pi }{3m_{\mathrm{Pl}}^2}\left[V(\varphi )\dot{\varphi }^2\right].$$
(3)
Here $`m_{\mathrm{Pl}}=G^{1/2}10^{19}\mathrm{GeV}`$ is the Planck mass. These background equations, along with the equation of motion (1), form a coupled set of differential equations that describe the evolution of the universe. The limit $`\dot{\varphi }=0`$ corresponds to a de Sitter universe, with the scale factor increasing exponentially in time

$`H=\sqrt{\left({\displaystyle \frac{8\pi }{3m_{\mathrm{Pl}}^2}}\right)V(\varphi )}=\mathrm{const},`$ (4)
$`ae^{Ht}.`$ (5)
In general $`H`$ is not exactly constant, but varies as the field $`\varphi `$ evolves along the potential $`V(\varphi )`$. A powerful way of describing the dynamics of an inflationary universe with a varying field (and non-trivial potential) is to express the Hubble parameter as a function of the field $`\varphi `$, $`H=H(\varphi )`$, which is consistent provided $`\varphi `$ is monotonic in time. The equations of motion become grishchuk88 ; muslimov90 ; salopek90 ; lidsey95

$`\dot{\varphi }`$ $`=`$ $`{\displaystyle \frac{m_{\mathrm{Pl}}^2}{4\pi }}H^{}(\varphi ),`$ (6)
$`\left[H^{}(\varphi )\right]^2{\displaystyle \frac{12\pi }{m_{\mathrm{Pl}}^2}}H^2(\varphi )`$ $`=`$ $`{\displaystyle \frac{32\pi ^2}{m_{\mathrm{Pl}}^4}}V(\varphi ).`$ (7)
These are completely equivalent to the second-order equation of motion (1). The second of the above equations is referred to as the Hamilton-Jacobi equation, and can be written in the useful form

$$H^2(\varphi )\left[1\frac{1}{3}ϵ(\varphi )\right]=\left(\frac{8\pi }{3m_{\mathrm{Pl}}^2}\right)V(\varphi ),$$
(8)
where $`ϵ`$ is defined to be

$$ϵ\frac{m_{\mathrm{Pl}}^2}{4\pi }\left(\frac{H^{}(\varphi )}{H(\varphi )}\right)^2.$$
(9)
The physical meaning of $`ϵ`$ can be seen by expressing Eq. (3) as

$$\left(\frac{\ddot{a}}{a}\right)=H^2(\varphi )\left[1ϵ(\varphi )\right],$$
(10)
so that the condition for inflation $`(\ddot{a}/a)>0`$ is given by $`ϵ<1`$. The scale factor is given by

$$ae^N=\mathrm{exp}\left[_{t_0}^tH𝑑t\right],$$
(11)
where the number of e-folds $`N`$ is

$$N_t^{t_e}H𝑑t=_\varphi ^{\varphi _e}\frac{H}{\dot{\varphi }}𝑑\varphi =\frac{2\sqrt{\pi }}{m_{\mathrm{Pl}}}_{\varphi _e}^\varphi \frac{d\varphi }{\sqrt{ϵ(\varphi )}}.$$
(12)
It is convenient to use $`N`$ as the measure of time during inflation. We take $`t_e`$ and $`\varphi _e`$ to be the time and field value at end of inflation. Therefore $`N`$ is defined as the number of e-folds before the end of inflation and increases as one goes backward in time, $`dt>0dN<0`$. The sign convention for $`\sqrt{ϵ}`$ must be applied carefully, and we take it to have the same sign as $`H^{}(\varphi )`$:

$$\sqrt{ϵ}+\frac{m_{\mathrm{PL}}}{2\sqrt{\pi }}\frac{H^{}}{H}.$$
(13)
Liddle, Parsons and Barrow use this as the starting point for an infinite hierarchy of slow roll parameters liddle94 :

$`\sigma `$ $``$ $`{\displaystyle \frac{m_{\mathrm{Pl}}}{\pi }}\left[{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{H^{\prime \prime }}{H}}\right)\left({\displaystyle \frac{H^{}}{H}}\right)^2\right],`$ (14)
$`{}_{}{}^{\mathrm{}}\lambda _{\mathrm{H}}^{}`$ $``$ $`\left({\displaystyle \frac{m_{\mathrm{Pl}}^2}{4\pi }}\right)^{\mathrm{}}{\displaystyle \frac{\left(H^{}\right)^\mathrm{}1}{H^{\mathrm{}}}}{\displaystyle \frac{d^{(\mathrm{}+1)}H}{d\varphi ^{(\mathrm{}+1)}}}.`$ (15)
The evolution of these parameters during inflation is determined by a set of “flow” equations hoffman00 ; kinney02 ,

$`{\displaystyle \frac{dϵ}{dN}}`$ $`=`$ $`ϵ\left(\sigma +2ϵ\right),`$ (16)
$`{\displaystyle \frac{d\sigma }{dN}}`$ $`=`$ $`5ϵ\sigma 12ϵ^2+2\left({}_{}{}^{2}\lambda _{\mathrm{H}}^{}\right),`$ (17)
$`{\displaystyle \frac{d\left({}_{}{}^{\mathrm{}}\lambda _{\mathrm{H}}^{}\right)}{dN}}`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}1}{2}}\sigma +\left(\mathrm{}2\right)ϵ\right]\left({}_{}{}^{\mathrm{}}\lambda _{\mathrm{H}}^{}\right)+{}_{}{}^{\mathrm{}+1}\lambda _{\mathrm{H}}^{}.`$ (18)
Here the time variable is the number of e-folds $`N`$, where

$$\frac{d}{dN}=\frac{d}{d\mathrm{log}a}=\frac{m_{\mathrm{Pl}}}{2\sqrt{\pi }}\sqrt{ϵ}\frac{d}{d\varphi }.$$
(19)
The derivative of a slow roll parameter at a given order is higher order in slow roll. Taken to infinite order, this set of equations completely specifies the cosmological evolution, up to the normalization of the Hubble parameter $`H`$. When truncated at finite order, by assuming that the $`{}_{}{}^{\mathrm{}}\lambda _{\mathrm{H}}^{}`$ are all zero above some fixed value of $`\mathrm{}`$, the solution of the flow equations still yields an exact solution of the background equations, albeit one that is drawn from a special subset of the overall solution space.

## III Monte Carlo reconstruction

The inflationary dynamics are fully specified by the values of the slow roll parameters $`[ϵ,\sigma ,{}_{}{}^{\mathrm{}}\lambda _{\mathrm{H}}^{}]`$ at a fixed time, which serve as initial conditions for the flow equations, truncated at order $`M`$ in the slow roll expansion. We now show that this information determines the inflationary potential $`V(\varphi )`$, up to a constant multiplier. The starting point is the Hamilton-Jacobi equation,

$$V(\varphi )=\left(\frac{3m_{\mathrm{Pl}}^2}{8\pi }\right)H^2(\varphi )\left[1\frac{1}{3}ϵ(\varphi )\right].$$
(20)
We have $`ϵ(N)`$ trivially from the flow equations. In order to calculate the potential, we need to determine $`H(N)`$ and $`\varphi (N)`$. With $`ϵ`$ known, $`H(N)`$ can be determined by inverting the definition of $`ϵ`$, Eq. (9):

$$\frac{1}{H}\frac{dH}{dN}=ϵ.$$
(21)
This can be viewed as the lowest-order member of the system of flow equations (18). Similarly, $`\varphi (N)`$ follows from the first Hamilton-Jacobi equation (7):

$$\frac{d\varphi }{dN}=\frac{m_{\mathrm{PL}}}{2\sqrt{\pi }}\sqrt{ϵ}.$$
(22)
Using these equations and Eq. (20), the form of the potential can then be fully reconstructed from the numerical solution for $`ϵ(N)`$.111Expressing the potential in this way is certainly not new. This is the basis of “functional reconstruction” method of Refs. hodges90 ; copeland93 , and well as the “Stewart-Lyth inverse problem” of Ref. beato00 . The only necessary observational input is the normalization of the Hubble parameter $`H`$, which enters the above equations as an integration constant. Here we use the simple condition that the density fluctuation amplitude (as determined by a first-order slow roll expression) be of order $`10^5`$,

$$\frac{\delta \rho }{\rho }\frac{H}{2\pi m_{\mathrm{Pl}}\sqrt{ϵ}}=10^5.$$
(23)
A more sophisticated treatment would perform a full normalization to the COBE CMB data bunn94 ; stompor95 . The value of the field, $`\varphi `$, also contains an arbitrary, additive constant.

Given a solution to the flow equations it is straightforward to determine the observational predictions for that model, even before the potential is computed. For instance, measurements of CMB anisotropies dodelson97 ; kinney98  may determine the tensor/scalar ratio $`r`$, the spectral index $`n`$, and the “running” of the spectral index $`dn/d\mathrm{log}k`$, and we focus on these quantities here. To lowest order, the relationship between the slow roll parameters and the observables is especially simple: $`r=ϵ`$, $`n1=\sigma `$, and $`dn/d\mathrm{log}k=0`$. To second order in slow roll, the observables are given by liddle94 ; stewart93 ,

$$r=ϵ\left[1C\left(\sigma +2ϵ\right)\right],$$
(24)
for the tensor/scalar ratio, and

$$n1=\sigma \left(53C\right)ϵ^2\frac{1}{4}\left(35C\right)\sigma ϵ+\frac{1}{2}\left(3C\right)\left({}_{}{}^{2}\lambda _{\mathrm{H}}^{}\right)$$
(25)
for the spectral index. The constant $`C4(\mathrm{ln}2+\gamma )`$, where $`\gamma 0.577`$ is Euler’s constant. Derivatives with respect to wavenumber $`k`$ can be expressed in terms of derivatives with respect to $`N`$ as liddle95

$$\frac{d}{dN}=\left(1ϵ\right)\frac{d}{d\mathrm{ln}k},$$
(26)
The scale dependence of $`n`$ is then given by the simple expression

$$\frac{dn}{d\mathrm{ln}k}=\left(\frac{1}{1ϵ}\right)\frac{dn}{dN},$$
(27)
which can be evaluated to third order in slow roll by using Eq. (25) and the flow equations. The final result following the evaluation of a particular path in $`M`$ dimensional “slow roll space” is a point in “observable parameter space”, i.e. $`(r,n,dn/d\mathrm{log}k)`$, corresponding to the observational prediction for that particular model. This process can be repeated for a large number of models, and used to study the attractor behavior of the inflationary dynamics. In fact, the models cluster strongly in the observable parameter space kinney02 .

We are now in a position to ask the key question: given an observational constraint, what dynamics – and underlying inflationary potential – are compatible with that constraint? We proceed by generating a large number of “universes” with random values (in a sense that is made explicit below) for the slow roll parameters. Given an allowed region in the observable parameter space, centered on specified values of $`r`$,$`n`$, and $`dn/d\mathrm{log}k`$ with a width given by the error bars on these quantities, we can construct an ensemble of inflationary potentials consistent with with the specified “window” in parameter space. This Monte Carlo approach to reconstruction differs from previously used techniques lidsey95 ; grivell00 ; escalante02  in that it generates an ensemble of potentials consistent with a given constraint, rather than attempting to produce a parameterized form for the potential and then constrain those parameters. The advantage of Monte Carlo reconstruction is that by evaluating the flow equations to high enough order we need only very mild a priori assumptions about the form of the potential. In addition, by probing the parameter space in a uniform way, we obtain what is in some sense a “fair sample” of the potentials consistent with a given observational bound. Even so, we do not have a metric on the space of initial conditions, so it is not possible to derive statistical inferences from the ensemble of potentials.

The condition for the end of inflation is that $`ϵ=1`$. Integrating the flow equations forward in time will yield two possible outcomes. One possibility is that the condition $`ϵ=1`$ may be satisfied for some finite value of $`N`$, which defines the end of inflation. We identify this point as $`N=0`$ so that the primordial fluctuations are actually generated when $`N50`$. Alternatively, the solution can evolve toward an inflationary attractor with $`r=0`$ and $`n>1`$, in which case inflation never stops.222See Ref. kinney02  for a detailed discussion of the fixed-point structure of the slow roll space. In reality, inflation must stop at some point, presumably via some sort of phase transition, such as the “hybrid” inflation mechanism linde91 ; linde94 ; copeland94 . Here we make the simplifying assumption that the observables for such models are the values at the late-time attractor.

To summarize, the algorithm for Monte Carlo reconstruction is as follows:

1. Specify a “window” of parameter space: e.g. central values for $`n1`$, $`r`$ or $`dn/d\mathrm{ln}k`$ and their associated error bars.
2. Select a random point in slow roll space, $`[ϵ,\eta ,{}_{}{}^{\mathrm{}}\lambda _{\mathrm{H}}^{}]`$, truncated at order $`M`$ in the slow roll expansion.
3. Evolve forward in time ($`dN<0`$) until either (a) inflation ends ($`ϵ>1`$), or (b) the evolution reaches a late-time fixed point ($`ϵ={}_{}{}^{\mathrm{}}\lambda _{\mathrm{H}}^{}=0,\sigma =\mathrm{const}.`$)
4. If the evolution reaches a late-time fixed point, calculate the observables $`r`$, $`n1`$, and $`dn/d\mathrm{ln}k`$ at this point.
5. If inflation ends, evaluate the flow equations backward $`N`$ e-folds from the end of inflation. Calculate the observable parameters at this point.
6. If the observable parameters lie within the specified window of parameter space, compute the potential and add this model to the ensemble of “reconstructed” potentials.
7. Repeat steps 2 through 6 until the desired number of models have been found.

In principle it is possible to carry out Monte Carlo reconstruction with no assumptions about the convergence of the hierarchy of slow roll parameters. In practice, the flow equations (18) must be truncated at some finite order and evaluated numerically. Moreover, for any given path in the parameter space, we do not know a priori the correct number of e-folds $`N`$ at which to evaluate the observables, since this depends on details such as the energy density during inflation and the reheat temperature lidsey95 . Consequently, after truncating to order $`M`$ in slow roll, we select our models’ parameters randomly from the following uniform distributions:

$`N`$ $`=`$ $`[40,70]`$ (28)
$`ϵ`$ $`=`$ $`[0,0.8]`$ (29)
$`\sigma `$ $`=`$ $`[0.5,0.5]`$ (30)
$`{}_{}{}^{2}\lambda _{\mathrm{H}}^{}`$ $`=`$ $`[0.05,0.05]`$ (31)
$`{}_{}{}^{3}\lambda _{\mathrm{H}}^{}`$ $`=`$ $`[0.025,0.025],`$ (32)
$`\mathrm{}`$ (33)
$`{}_{}{}^{M+1}\lambda _{\mathrm{H}}^{}`$ $`=`$ $`0.`$ (34)
and so forth, reducing the width of the range by factor of five for each higher order in slow roll. The series is closed to order $`M`$ by taking $`{}_{}{}^{M+1}\lambda _{\mathrm{H}}^{}=0`$. We use $`M=5`$ for the calculations in this paper. The exact choice of ranges for the initial parameters does not have a large influence on the result of the Monte Carlo process, as long as they are chosen such that the slow roll hierarchy is convergent. As noted above, the form of the flow equations (18) ensures that the derivative of $`{}_{}{}^{\mathrm{}}\lambda _{\mathrm{H}}^{}`$ depends only on parameters of order $`\mathrm{}`$ and $`\mathrm{}+1`$, so this truncation still leads to an exact evaluation of the flow equations to infinite order. We are selecting a finite subset out of an infinite number of possibilities for initial conditions, but the background evolution for a given model is evaluated exactly.333As usual in all analyses of inflation we ignore the back reaction of quantum fluctuations on the background evolution.

The result of integrating the flow equations for a particular model is a “path” in the slow roll space parameterized by the number of e-folds $`N`$. Figures 1-3 show examples of paths plotted in the $`\sigma ϵ`$ plane for different assumptions about the central values for the parameters $`r`$, $`n`$, and $`dn/d\mathrm{log}k`$. Also plotted is the corresponding reconstructed potential for each path. Surprisingly, while quite complex behavior is possible for the slow roll parameters, this behavior is only weakly reflected in the shape of the potential itself.

## IV Monte Carlo Reconstruction in Practice

We now describe two concrete applications of this method. First, we select central values in various regions of the observable parameter space with the error bars expected from Planck, and reconstruct the inflationary potential based on this “synthetic” data. We find that Planck will allow us to determine the qualitative form of the potential, but is insufficiently precise for well-constrained reconstruction. This is consistent with the conclusions of other studies lidsey95 ; grivell00 . Second, we choose a region of parameter space centered on the values for $`n`$, $`r`$ and $`dn/d\mathrm{ln}k`$ for a specified potential – in this case $`V(\varphi )\varphi ^4`$ – and apply our reconstruction algorithm with different sized error bars, thus determining the observational precision that is needed in order to reliably reconstruct a known potential.

Figure 4 shows the results of two different “reconstructions”: for $`r=0.02`$, $`n=0.95`$ and $`dn/d\mathrm{log}k=0`$, and for $`r=0`$, $`n=.93`$ and $`dn/d\mathrm{log}k=0`$. In both cases the errors bars are those anticipated for the Planck mission’s measurement of the spectrum kinney98 ; copeland97 , namely $`\delta r0.01`$, $`\delta n0.01`$ and $`\delta dn/d\mathrm{log}k0.01`$. The first case has a tensor contribution that would be resolved by Planck, and the normalization of the potential is relatively tightly constrained, $`V10^{11}m_{\mathrm{Pl}}^4`$. The second reconstruction assumes that the tensor amplitude is not detectable by Planck, with a central value $`r0.0`$. In this case, the normalization of the potential is very poorly constrained, although the shape of the potential is consistent with standard “small-field” models dodelson97 . The clear conclusion from these plots, and the other cases that we have examined, is that Planck may be able to determine the qualitative features of the potential, but it will not permit the quantitative reconstruction of the inflationary potential on its own.

It should not come as any particular surprise that Monte Carlo reconstruction does not yield a tight constraint on the possible form of the inflationary potential. However, the question this immediately raises is how much more accurately – relative to the precision expected from Planck – would we need to measure the spectrum in order to be able to make a quantitative statement about the functional form of the potential. We tackle this issue by choosing values of $`n`$, $`r`$ and $`dn/d\mathrm{log}k`$ consistent with a known potential, in this case the canonical $`\lambda \varphi ^4`$ model, which is typical of “large-field” models and has a tensor fluctuation amplitude observable by Planck. With this assumption, the central values for the spectral index and tensor-to-scalar ratio are $`n=0.943`$ and $`r=0.02`$ (evaluated 50 e-folding before the end of inflation). For $`V\varphi ^4`$, $`dn/d\mathrm{log}k`$ is close to zero, and we take the central value to be precisely zero.

Using the above choices for the central values, we have performed Monte Carlo reconstruction with Planck-sized errors bars and with $`\delta n=0.002`$, $`\delta r=0.002`$ and $`\delta dn/d\mathrm{log}k=0.005`$. The second choice leads to a box in the $`(n,r)`$ plane that is 25 times smaller than the bound predicted for Planck and thus represents a substantial improvement in experimental precision. We divided the $`dn/d\mathrm{log}k`$ bound by 2 rather than 5 solely for computational convenience. The second reconstruction involved computing approximately 90 million “trial” inflationary models in order to filter out 100 solutions that passed through the specified window in parameter space. This took several days of CPU time on a fast desktop machine, and would have taken even longer if we had applied the same scaling to $`\delta dn/d\mathrm{log}k`$ as we did to $`\delta n`$ and $`\delta r`$.

If we plot graphs to analogous to those in Fig. 4 for these two reconstructions we find that, unsurprisingly, the overall shape of both sets of reconstructed potentials are consistent, but that there is less spread in the set with the smaller error bars. However, to determine whether any sort of quantitative reconstruction is possible, we adopt the prior that the potentials are proportional to $`(\varphi \varphi _0)^m`$, where $`\varphi _0`$ is an unknown offset, and then perform a least-squares fit to determine the best-fit value of $`m`$. In practice, it is sufficient to normalize the height of the reconstructed potential to precisely unity at $`\varphi =0`$ and then fit it to the following functional form:

$$(1c\varphi )^m$$
(35)
This choice forces the fitted (and rescaled) potential to be unity when $`\varphi =0`$. Adding normalization as a third parameter to the fit does not have a significant impact on the computed values of $`m`$.

The results for both simulations are shown in Fig. 5. The error bars expected from the Planck mission do not allow one to conclude that $`V\varphi ^4`$, but the tighter bounds on $`n`$ and $`r`$ do rule out values of $`m`$ markedly different from 4.444In reality, if the computed value of $`m`$ is significantly different from 4, what we actually learn is that the potential cannot be well described by Eq. (35) – since setting $`m`$ to a value significantly different from 4 produces values of $`n`$ and $`r`$ well outside our assumed range.

We thus tentatively conclude that while Planck cannot measure the perturbation spectrum accurately enough to put even mild quantitative constraints on the potential, increasing the accuracy with which both $`n`$ and $`r`$ can be recovered by a factor of 5 would allow one to start making meaningful reconstructions of the potential.

## V Discussion

We have presented a novel method for reconstructing the inflationary potential, based on the powerful flow equation approach to inflationary dynamics. This method, which we dub Monte Carlo reconstruction, differs from other prescriptions for reconstruction in that it is stochastic in nature and involves extremely weak prior assumptions about the form of the potential. A stochastic approach to reconstruction has distinct advantages compared with, for example, methods which expand the potential in a Taylor series and fit the coefficients of the Taylor expansion to observable parameters. In particular, the method can be applied equally well to both poor data sets and to high-quality data: parameters such as the tensor spectral index do not enter the reconstruction in any direct way, making the method very simple and robust. Also, Monte Carlo reconstruction naturally produces shapes for the potential outside the region directly constrained by observation. This follows from our principal assumption, single field inflation, and the attractor behavior of the inflationary dynamics which tends to ensure that the reconstructed potentials overlap outside the region in which we have direct observational input.

We want to make clear what Monte Carlo reconstruction is not: it is a stochastic method, but not a statistical one. We do not have a metric on the space of initial conditions. Consequently, we cannot use the “density” of models in any particular parameter space to infer the relative likelihood of one parameter region over another. Plots of models in either the space of possible potentials or the space of exponent $`m`$ (from $`\varphi ^m`$, as in Fig. 5) are properly interpreted as exclusion plots, indicating which regions are either consistent or inconsistent with the data. However, we cannot determine the relative likelihood of different initial points in slow roll space without an understanding of inflationary initial conditions. Many other reconstruction attempts truncate the slow roll expansion at second or third order, but the approach here can be extended to arbitrary order in slow roll, and we have checked that our results do not depend on the specific level at which we truncate slow roll. However, any slow roll ansatz effectively rules out potentials with small features where the higher order derivatives of the potential and $`H(\varphi )`$ are large wang97 . If the feature is located in the region of the potential that corresponds to the primordial spectrum, strict limits can be placed on the size and slope of the feature adams01 . However, if the feature is outside this region then we will obviously not be able to reconstruct it.

The most straightforward application of Monte Carlo reconstruction is to simply generate an ensemble of potentials consistent with some observational constraint. This simple and robust procedure can be applied to data sets (for example the forthcoming data from the MAP satellite) for which direct functional reconstruction of the potential will likely be impossible. While this method suffers from the same fundamental limitations as any other method, it will allow us to at least answer qualitative questions about the form of the inflationary potential: for example, is the potential convex or concave? Higher quality data will allow for more quantitative constraints on the potential. Detection of a nonzero tensor component will give information about the normalization of the potential, to within an order of magnitude or two with the accuracy projected for Planck. Even absent a detection of tensor modes, it may be possible to reach conclusions about the shape of the potential, if not its normalization. More quantitative information can be gleaned from more accurate data sets. For instance, an improvement in parameter resolution by a factor of five or so over Planck will make information about the exponent $`\varphi ^m`$ of the potential available in at least a rough sense. Other possible applications would be plotting the results of the Monte Carlo reconstruction in the space of derivatives of the potential, $`V,V^{},V^{\prime \prime }`$, and so on, for comparison with Refs. hansen01 ; caprini02 .

Monte Carlo reconstruction is easily generalizable. For example, one need not use the slow roll approximation to calculate the power spectrum associated with a particular choice of “slow roll” parameters generated by the flow equations. (Note that the slow roll expansion is completely distinct from the slow roll approximation. The solutions we generate for the background evolution are exact.) For calculating the fluctuation spectrum associated with a particular path in the parameter space, one could equally well apply the method of uniform approximations introduced by Habib, et al. habib02 . It is also straightforward to solve for the perturbation spectrum by numerically solving the exact equation for quantum modes in the inflationary spacetime: as in the calculation of the potential itself, all the information required to do so is contained in the solution to the flow equations. Thus the method is not only exact in principle, but can be made so in practice if necessary. In addition, the same techniques could be applied to a set of flow equations based around an expansion other than the Hubble slow roll expansion of Liddle et al., such as the expansion used in Ref. schwarz01 . Doing so would be useful to investigate the dependence of the reconstructed potentials on the details of the truncation scheme for the flow equations. We expect this dependence to be small.

## Acknowledgments

RE and WHK are supported by ISCAP and the Columbia University Academic Quality Fund. ISCAP gratefully acknowledges the generous support of the Ohrstrom Foundation. We would like to thank Ed Copeland for useful discussions. RE thanks the Aspen Center for Physics, where part of this work was conducted. Some of the computational time used for these calculations were kindly provided by the High Energy Theory group at Brown University.


format: text

text:
# Dynamic Exchange Coupling in Magnetic Bilayers

## Abstract

A long-ranged dynamic interaction between ferromagnetic films separated by normal-metal spacers is reported, which is communicated by nonequilibrium spin currents. It is measured by ferromagnetic resonance (FMR) and explained by an adiabatic spin-pump theory. In FMR the spin-pump mechanism of spatially separated magnetic moments leads to an appreciable increase in the FMR line width when the resonance fields are well apart, and results in a dramatic line-width narrowing when the FMR fields approach each other.

The giant magnetoresistance Baibich et al. (1988) accompanying realignment of magnetic configurations in metallic multilayers by an external magnetic field is routinely employed in magnetic read heads and is essential for high-density nonvolatile magnetic random-access memories. These typically consist of ferromagnetic/normal/ferromagnetic (F/N/F) metal hybrid structures, i.e., magnetic bilayers which are an essential building block of the so called spin valves. The static Ruderman-Kittel-Kasuya-Yosida (RKKY) interlayer exchange between ferromagnets in magnetic multilayers Grünberg et al. (1986) is suppressed in these devices by a sufficiently thick nonmagnetic spacer N or a tunnel barrier. The interest of the community shifts increasingly from the static to the dynamic properties of the magnetization Back et al. (1999). This is partly motivated by curiosity, partly by the fact that the magnetization switching characteristics in memory devices is a real technological issue. A good grasp of the fundamental physics of the magnetization dynamics becomes of essential importance to sustain the exponential growth of device performance factors.

In this Letter we study the largely unexplored dynamics of magnetic bilayers in a regime when there is no discernible static interaction between the magnetization vectors. Surprisingly, the magnetizations still turn out to be coupled, which we explain by emission and absorption of nonequilibrium spin currents. Under special conditions the two magnetizations are resonantly coupled by spin currents and carry out a synchronous motion, quite analogous to two connected pendulums. This dynamic interaction is an entirely new concept and physically very different from the static RKKY coupling. E.g., the former does not oscillate as a function of thickness and its range is exponentially limited by the spin-flip relaxation length of spacer layers and algebraically by the elastic mean free path. This coupling can have profound effects on magnetic relaxation and switching behavior in hybrid structures and devices.

The unit vector $`𝐦=𝐌/M`$ of the magnetization $`𝐌(t)`$ of a ferromagnet changes its direction in the presence of a noncollinear magnetic field. The motion of $`𝐦`$ in a single domain is described by the Landau-Lifshitz-Gilbert (LLG) equation

$$\frac{d𝐦}{dt}=\gamma 𝐦\times 𝐇_{\text{eff}}+\alpha 𝐦\times \frac{d𝐦}{dt},$$
(1)
with $`\gamma `$ being the absolute value of the gyromagnetic ratio. The first term on the right-hand side represents the torque induced by the effective magnetic field $`𝐇_{\text{eff}}=F/𝐌`$, where the free-energy functional $`F[𝐌]`$ consists of the Zeeman energy, magnetic anisotropies, and exchange interactions Heinrich and Cochran (1993). The second term in Eq. (1) is the Gilbert damping torque which governs the relaxation towards equilibrium. The intrinsic damping in bulk metallic ferromagnets, $`\alpha ^{(0)}`$, typically 0.002-0.025, appears to be governed by spin-orbit interactions Kunes and Kambersky (2002) in the 3d transition metals. The magnetization vector can be forced into a resonant precession motion by microwave stimulation. This ferromagnetic resonance (FMR) is measured via the absorption of microwave power using a small rf field at a frequency $`\omega `$ polarized perpendicular to the static magnetic moment as a function of the applied dc magnetic field, see the right inset in Fig. 1. The absorption is given by the imaginary part of the susceptibility $`\chi ^{\prime \prime }`$ of the rf magnetization component along the rf driving field. This FMR signal has a Lorentzian line shape with a width $`\mathrm{\Delta }H=(2/\sqrt{3})\alpha \omega /\gamma `$ when defined by the inflection points (i.e., the extrema of $`d\chi ^{\prime \prime }/dH`$), see the left inset in Fig. 1.

When two or more ferromagnets are in electrical contact via nonmagnetic metal layers, interesting new effects occur. Transport of spins accompanying an applied electric current driven through a magnetic multilayer causes a torque on the magnetizations Slonczewski (1996), which at sufficiently high current densities leads to spontaneous magnetization-precession and switching phenomena Myers et al. (1999). Even in the absence of an applied charge current, spins are injected into the normal metal by a ferromagnet with moving magnetization. This causes additional magnetic damping, provided that the spin-flip relaxation rate of normal metal is high Tserkovnyak et al. (2002a). The present Letter focuses on the discovery of novel dynamic effects in F1/N/F2 structures in the limit when the spin-flip scattering in N is weak. Let us first sketch the basic physics. A precessing magnetization $`𝐦_i`$ “pumps” a spin current $`𝐈_{si}^{\text{pump}}𝐦_i`$ into the normal metal Tserkovnyak et al. (2002a). We focus on weakly excited magnetic bilayers close to the parallel alignment, so that $`𝐈_{si}^{\text{pump}}𝐦_j`$ for arbitrary $`i,j=1,2`$. The spin momentum perpendicular to the magnetization direction cannot penetrate a ferromagnetic film beyond the (transverse) spin-coherence length, $`\lambda _{\text{sc}}=\pi /|k_\text{F}^{}k_\text{F}^{}|`$, which is determined by the spin-dependent Fermi wave vectors $`k_\text{F}^,`$ and is smaller than a nanometer for 3d metals Stiles and Zangwill (2002). A transverse spin current ejected by one ferromagnet can therefore be absorbed at the interface to the neighboring ferromagnet, thereby exerting a torque $`𝝉`$. Each magnet thus acts as a spin sink which can dissipate the transverse spin current ejected by the other layer.

The theoretical basis of this picture is the adiabatic spin-pumping mechanism Tserkovnyak et al. (2002a) and magnetoelectronic circuit theory Brataas et al. (2000). N is assumed thick enough to suppress any RKKY Grünberg et al. (1986), pin-hole Bobo et al. (1999), and magnetostatic (Néel-type) Rücker et al. (1995) interactions. We consider ultrathin films with a constant magnetization vector across the film thickness Heinrich and Cochran (1993), which are nonetheless thicker than $`\lambda _{\text{sc}}`$ and, therefore, completely absorb transverse spin currents. In the experiments described below, N is thinner than the electron mean free path, so that the electron motion inside the spacer is ballistic. Precessing $`𝐦_i`$ pumps spin angular momentum at the rate Tserkovnyak et al. (2002a)

$$𝐈_{si}^{\text{pump}}=\frac{\mathrm{}}{4\pi }g^{}𝐦_i\times \frac{d𝐦_i}{dt},$$
(2)
where $`g^{}`$ is the dimensionless “mixing” conductance Brataas et al. (2000) of the F/N interfaces, which can be obtained via ab initio calculations of the scattering matrix Xia et al. (2002) or measured via the angular magnetoresistance of spin valves Bauer et al. (2002) as well as FMR line widths of F/N and F/N/F magnetic structures Tserkovnyak et al. (2002a); Mizukami et al. (2001a); Urban et al. (2001). Note that $`g^{}`$ must be renormalized for the intermetallic interfaces considered here Bauer et al. (2002). We assume identical Fi/N interfaces with real-valued $`g^{}`$, as suggested by calculations for various F/N combinations Xia et al. (2002). When the spacer is not ballistic, its diffuse resistance can simply be absorbed into the value of $`g^{}`$, which should then be interpreted as the mixing conductance of an F/N interface in series with the half of the spacer. When, furthermore, the spacer is thicker than the spin-diffusion length, the spin-pumping exchange between the magnetic layers becomes exponentially suppressed with the spacer thickness Tserkovnyak et al. (2002a).

Alloy disorder at the interfaces scrambles the distribution function. Disregarding spin-flip scattering in the normal metal, an incoming spin current on one side leaves the normal-metal node by equal outgoing spin currents to the right and left Bauer et al. (2002). (As the interfacial scrambling is only partial and the spacer is ballistic, the last statement should not be taken literally, but as an effective theory which is valid after renormalizing the interfacial conductance parameters.) On typical FMR time scales, this process occurs practically instantaneously. The net spin torque at one interface is therefore just the difference of the pumped spin currents divided by two: $`𝝉_1=(𝐈_{s2}^{\text{pump}}𝐈_{s1}^{\text{pump}})/2=𝝉_2`$. When one ferromagnet is stationary, see the left drawing in Fig. 2, the dynamics of the other film, Fi, is governed by the LLG equation with a damping parameter $`\alpha _i=\alpha _i^{(0)}+\alpha _i^{}`$ enhanced with respect to the intrinsic value $`\alpha _i^{(0)}`$ by $`\alpha _i^{}=\gamma \mathrm{}g^{}/(8\pi \mu _i)`$, where $`\mu _i`$ is the total magnetic moment of Fi. Since $`\mu _i`$ scales linearly with the volume of the ferromagnet and $`g^{}`$ scales with its interface area, $`\alpha _i^{}`$ is inversely proportional to the film thickness.

When both magnetizations are allowed to precess, see the right drawing in Fig. 2, the LLG equation expanded to include the spin torque reads

$`{\displaystyle \frac{d𝐦_i}{dt}}`$ $`=`$ $`\gamma 𝐦_i\times 𝐇_{\text{eff}}^i+\alpha _i^{(0)}𝐦_i\times {\displaystyle \frac{d𝐦_i}{dt}}`$ (3)
$`+\alpha _i^{}\left[𝐦_i\times {\displaystyle \frac{d𝐦_i}{dt}}𝐦_j\times {\displaystyle \frac{d𝐦_j}{dt}}\right],`$
where $`j=1(2)`$ if $`i=2(1)`$. As a simple example, consider a system in the parallel configuration, $`𝐦_1^{(0)}=𝐦_2^{(0)}`$, with matched resonance conditions. In addition, let us assume the resonance precession is circular. If we linearize Eq. (3) in terms of small deviations $`𝐮_i(t)=𝐦_i(t)𝐦_i^{(0)}`$ of the magnetization direction $`𝐦_i`$ from its equilibrium value $`𝐦_i^{(0)}`$, we find that the average magnetization deviation $`𝐮=(𝐮_1\mu _1+𝐮_2\mu _2)/(\mu _1+\mu _2)`$ is damped with the intrinsic Gilbert parameter $`\alpha ^{(0)}`$, whereas the difference $`\mathrm{\Delta }𝐮=𝐮_1𝐮_2`$ relaxes with enhanced damping constant $`\alpha =\alpha ^{(0)}+\alpha _1^{}+\alpha _2^{}`$.

Measuring the spin torques requires independent control of the precessional motion of the two F layers, with FMR absorption line widths of isolated films dominated by the intrinsic Gilbert damping. Both conditions were met by high-quality crystalline Fe(001) films grown on 4x6 reconstructed GaAs(001) substrates by Molecular Beam Epitaxy Urban et al. (2001); Enders et al. (2001). Fe(001) films were deposited at room temperature (RT) from a thermal source at a base pressure of less then $`2\times 10^{10}`$ Torr and the deposition rate was $`1`$ ML/min. For the experiments discussed below, single Fe ultrathin films with thicknesses $`d_F=`$ 11,16,21,31 ML were grown directly on GaAs(001) and covered by a 20 ML protective Au(001) cap layer. The magnetic anisotropies as measured by FMR are described by a constant bulk term and an interface contribution inversely proportional to $`d_F`$. The Fe ultrathin films grown on GaAs(001) and covered by gold have magnetic properties nearly identical to those in bulk Fe, modified only by sharply defined interface anisotropies. The in-plane uniaxial anisotropy arises from electron hybridization between the As dangling bonds and the iron interface atoms. These Fe films were then regrown as one element of a magnetic bilayer structure and in the following referred to as F1 layers. They were separated from a thick Fe layer, F2, of 40 ML thickness by a 40 ML Au spacer. The magnetic bilayers were covered by 20 ML of protective Au(001). The complete structures are therefore GaAs/Fe(8,11,16,21,31)/40Au/40Fe/20Au(001), where the integers represent the number of MLs. The electron mean free path in thick films of gold is 38 nm Enders et al. (2001) and, consequently, the spin transport even in the 40 ML (8 nm) Au spacer is purely ballistic. The interface magnetic anisotropies allowed us to separate the FMR fields of the two Fe layers with resonance-field differences that can exceed 5 times the FMR line widths, see Fig. 1. Hence, the FMR measurements for F1 in double layers can be carried out with a nearly static F2.

The FMR line width of F1 increases in the presence of F2. The difference $`\mathrm{\Delta }H^{}`$ in the FMR line widths between the magnetic bilayer and single-layer structures is nearly inversely proportional to the thin-film thickness $`d_F`$ Urban et al. (2001), proving that $`\mathrm{\Delta }H^{}`$ originates at the F1/N interface. Secondly, $`\mathrm{\Delta }H^{}`$ is linearly dependent on microwave frequency for both the in-plane (the saturation magnetization parallel to the film surface) and perpendicular (the saturation magnetization perpendicular to the film surface) configurations, strongly implying that the additional contribution to the FMR line width can be described strictly as an interface Gilbert damping Urban et al. (2001). At the FMR, the film precessions are driven by an applied rf field. When the resonance fields are different, one layer (say F1) is at resonance with maximum precessional amplitude while the other layer (F2) is off resonance with small precessional amplitude, see Fig. 2. The spin-pump current for F1 reaches its maximum while F2 does not emit a significant spin current at all. F2 acts as a spin sink causing the nonlocal damping for F1. The N/F2 interface provides a “spin-momentum brake” for the F1 magnetization. The corresponding additional Gilbert parameter $`\alpha ^{}`$ for a 16 ML Fe is significant, being similar in magnitude to the intrinsic Gilbert damping in isolated Fe films, $`\alpha ^{(0)}=0.0044`$.

These assertions can be tested by employing the in-plane uniaxial anisotropy in F1 to intentionally tune the resonance fields for F1 and F2 into a crossover which is shown in the shaded area of Fig. 1. When the resonance fields are identical, $`H_1=H_2`$, the rf magnetization components of F1 and F2 are parallel to each other, see the right drawing in Fig. 2. The total spin currents across the F1/N and N/F2 interfaces therefore vanish resulting in zero excess damping for F1 and F2, see Eq. (3), which is experimentally verified, as shown in Fig. 3. For a theoretical analysis, we solved Eq. (3) and determined the total FMR signal as a function of the difference between the resonance fields $`H_2H_1`$. The theoretical predictions are compared with measurements in Fig. 3. The remarkable good agreement between the experimental results and theoretical predictions provides strong evidence that the dynamic exchange coupling not only contributes to the damping but leads to a new collective behavior of magnetic hybrid structures.

We have additionally carried out our measurements on samples with Au spacer thickness between 14 and 100 monolayers. The weak dependence of the FMR response on the spacer thickness fully supports our picture of the long-ranged dynamic interaction.

In conclusion, we found decisive experimental and theoretical evidence for a new type of exchange interaction between ferromagnetic films coupled via normal metals. In contrast to the well-known oscillatory exchange interaction in the ground state, this coupling is dynamic in nature and long ranged. Precessing magnetizations feel each other through the spacer by exchanging nonequilibrium spin currents. When the resonance frequencies of the ferromagnetic banks differ, their motion remains asynchronous and net spin currents persist. However, when the ferromagnets have identical resonance frequencies, the coupling quickly synchronizes their motion and equalizes the spin currents. Since these currents flow in opposite directions, the net flow across both F1/N and N/F2 interfaces vanishes in this case. The lifetime of the arising collective motion is limited only by the intrinsic local damping. These effects can be well demonstrated in FMR measurements.

We are grateful to B. I. Halperin for stimulating discussions. This work was supported in part by the NEDO International Joint Research Grant Program “Nano-magnetoelectronics”, NSF Grant DMR 02-33773, and the FOM.


text:
# Violation of Ioffe-Regel condition but saturation of resistivity of the high 𝑇_𝑐 cuprates

## Abstract

We demonstrate that the resistivity data of a number of high $`T_c`$ cuprates, in particular La2-xSrxCuO4, are consistent with resistivity saturation, although the Ioffe-Regel condition is strongly violated. By using the f-sum rule together with calculations of the kinetic energy in the $`tJ`$ model, we show that the saturation resistivity is unusually large. This is related to the strong reduction of the kinetic energy due to strong correlation effects. The fulfilment of the Ioffe-Regel condition for conventional transition metal compounds is found to be somewhat accidental.

Some strongly correlated metals, in particular many high $`T_c`$ cuprates, show an exceptionally large resistivity, $`\rho `$, at high temperatures $`T`$. Thus $`\rho `$ can reach values of several m$`\mathrm{\Omega }`$cm. This is in strong contrast to almost all other metals. Typically, a metal with a very large resistivity shows resistivity saturation. Thus when $`\rho `$ reaches values of the order 0.1 m$`\mathrm{\Omega }`$cm, the slope of $`\rho (T)`$ is typically reduced substantially, although not necessarily to zero. This happens when the apparent mean free path $`l`$ becomes comparable to the separation, $`d`$, of two atoms, the Ioffe-Regel condition. This kind of behavior has been observed for many metals and it used to be considered a universal behavior.

The metals mentioned above and the alkali-doped fullerenes are apparent exceptions to this saturation behavior, with the cuprates forming a particular challenge. The deviation from the “universal” behavior can be illustrated by considering Ioffe-Regel condition for La2-xSrxCuO4. We assume a “large” Fermi surface of cylindrical shape, containing $`1x`$ electrons (or $`1+x`$ holes). Assuming that $`l=a`$, where $`a`$ is the lattice parameter in the CuO2 plane, we obtain the saturation resistivity

$$\rho _{\mathrm{sat}}=\frac{0.7}{\sqrt{1\pm x}}\mathrm{m}\mathrm{\Omega }\mathrm{cm}.$$
(1)
If we instead assume a “small” Fermi surface, $`\sqrt{1\pm x}`$ in the denominator is replaced by $`\sqrt{x}`$. In either case, the experimental resistivity for small $`x`$ is much larger than the saturation resistivity above, supporting the conclusion that there is no saturation in these systems.

Here we demonstrate that the experimental data for the strongly correlated high $`T_c`$ cuprates are, nevertheless, consistent with saturation. The saturation resistivity can, however, be much larger than the Ioffe-Regel value. This is not entirely surprising, since the Ioffe-Regel condition is based on a semiclassical picture, which is not valid when $`ld`$. The good agreement between the Ioffe-Regel condition and the saturation for weakly correlated transition metal compounds is therefore somewhat accidental. The large saturation resistivity is shown to be due to a large reduction of the kinetic energy, due to strong correlation effects. Since the alkali-doped fullerenes appear to lack saturation, this puts these compounds in a special class, different from the cuprates.

To discuss resistivity saturation, we use the f-sum rule in a form appropriate for the models discussed here (only nearest neighbor hopping and no on-site matrix elements of the current operator)

$$\frac{2}{\pi }_0^{\mathrm{}}\sigma (\omega )𝑑\omega =\frac{1}{2}\frac{d^2e^2}{N\mathrm{\Omega }\mathrm{}^2}T_K,$$
(2)
where $`\sigma (\omega )`$ is the optical conductivity, $`d`$ is the separation of the sites, $`\mathrm{\Omega }`$ is the volume of a unit cell, $`N`$ is the number of unit cells and $`T_K`$ is the expectation value of the kinetic energy operator. The prefactor $`1/2`$ refers to the two-dimensional case and is replaced by $`1/3`$ for the three-dimensional case. To obtain an approximate upper limit to the resistivity, we assume that the (Drude) peak at $`\omega =0`$ has been smeared out and that $`\sigma (\omega )`$ is a smooth function. The removal of the $`\omega =0`$ peak may be due to any scattering mechanism, electron-phonon, electron-electron or impurity scattering and, depending on the system, it may or may not happen at large values of $`T`$. We emphasize that our theory below does not depend on any particular scattering mechanism. We furthermore assume that $`\sigma (\omega )=0`$ for $`|\omega |>W`$, where $`W`$ is the one-particle band width. This should be a good approximation, since $`|\omega |>W`$ would involve multiple electron-hole pair excitations and have a small weight. Explicit calculations for a model of transition metal compounds and for the $`tJ`$ model confirm this. If $`\sigma (\omega )`$ had a box shape, $`\sigma (\omega )\sigma (0)`$ for $`0|\omega |W`$, the integral over $`\sigma (\omega )`$ would be $`\sigma (0)W`$ and we could write

$$\sigma (\omega =0)=\frac{\gamma \mathrm{}}{W}_0^{\mathrm{}}\sigma (\omega )𝑑\omega ,$$
(3)
with $`\gamma =1`$. For a more realistic shape of $`\sigma (\omega )`$, where $`\sigma (\omega )`$ has a maximum for $`\omega =0`$, it follows that $`\gamma >1`$. Following the considerations for a model of transition metal compounds with electron-phonon scattering, we may set $`\gamma =1.9`$. In the $`tJ`$ model, considering electron-electron scattering, we find that $`\sigma (\omega )`$ is more concentrated to small $`\omega `$, giving a larger value of $`\gamma `$. If the electrons would tend to localize, $`\gamma `$ could be much smaller than unity. This should, however, not happen at the large values of $`T`$ considered here. Measurements of $`\sigma (\omega )`$ for La2-xSrxCuO4 at large $`T`$ show a similar behavior as assumed in Eq. (3), apart from some sharp phonon structures at finite $`\omega `$.

We first consider noninteracting electrons in a band with the orbital degeneracy $`n`$. The scattering is assumed to be due to the electron-phonon interaction. For $`TW`$ and neglecting the influence of the electron-phonon coupling on the kinetic energy,

$$T_K=2n_{W/2}^\mu \epsilon N(\epsilon )𝑑\epsilon 2n\alpha WN,$$
(4)
where $`N(\epsilon )`$ is the density of states per orbital and spin, $`\mu `$ is the chemical potential and $`\alpha `$ depends on the shape $`N(\epsilon )`$ and the filling. For semi-elliptical $`N(\epsilon )`$ and filling 0.4, $`\alpha =0.10`$. This result is relatively independent of the specific shape of $`N(\epsilon )`$, and using, e.g., a constant $`N(\epsilon )`$ leads to $`\alpha =0.12`$. For a three-dimensional system we then obtain

$$\sigma (0)=\frac{\pi n}{3}\alpha \gamma \frac{d^3}{\mathrm{\Omega }}\frac{e^2}{\mathrm{}d},$$
(5)
which provides an approximate upper limit to the resistivity for $`TW`$. Saturation happens if $`\rho (T)`$ grows so rapidly that this approximate upper limit is approached for experimentally accessible values of $`T`$. This happens for some transition metal compounds, with the A15 compounds, such as Nb3Sb, being particularly striking examples. The important orbitals in Nb3Sb are the Nb $`d`$-orbitals, leading to $`n=5`$. Considering the A15 lattice of Nb3Sb, we obtain the saturation resistivity 0.14 m$`\mathrm{\Omega }`$cm, which compares well with the experimental saturation resistivity of about 0.15 m$`\mathrm{\Omega }`$cm. This illustrates the usefulness of this approach for obtaining a saturation resistivity.

We next consider the strongly correlated cuprates. To calculate the kinetic energy $`T_K`$, we need an appropriate model for these systems. The transport properties should mainly be determined by the Cu-O antibonding band of Cu $`x^2y^2`$ and O $`2p`$ character. Thus we consider one orbital per site and the orbital degeneracy $`n=1`$. These orbitals are put on a two-dimensional square lattice, and we consider a a layered system with the lattice parameters ($`a`$,$`a`$,$`2c`$). We further assume that the Coulomb interaction $`U`$ between two electrons on the same site is so large that double occupancy of a site can be neglected. This can be described by the $`tJ`$ model. The properties of this model have been studied extensively by Jaklic and Prelovsek. The system has the hole doping $`x`$ $`(>0)`$. We emphasize that the $`tJ`$ model is only used to estimate $`T_K`$. We do not exclude other scattering mechanisms; we just neglect their contribution to the kinetic energy.

For simplicity, we first consider $`J=0`$. By expanding in $`1/T`$, we obtain the large $`T`$ limit of the kinetic energy as

$$T_K=\frac{\mathrm{Tr}T_Ke^{H/T}}{\mathrm{Tr}e^{H/T}}\frac{1}{T}[\frac{\mathrm{Tr}T_K^2}{\mathrm{Tr}1}]_{\mathrm{no}\mathrm{doubble}\mathrm{occ}},$$
(6)
where the right hand site is evaluated for states without double occupancy. Since the expression contains $`T_K^2`$, this involves a hole hopping to a neighboring site and back. In the $`tJ`$ model this is only possible if the neighboring site has no holes. The probability for having a hole or no hole on a given site is $`x`$ and $`(1x)`$, respectively. The probability for hopping between two nearest neighbor sites is then approximately $`x(1x)`$. Summing over the four nearest neighbors, this gives the large $`T`$ limit

$$T_K=4Nt^2\frac{x(1x)}{T}.$$
(7)
Using similar arguments, we can obtain a rough estimate for the $`T=0`$ kinetic energy

$$[T_K]_{T=0}=4tx(1x)N.$$
(8)
This result is shown by the full curve in Fig. 1. A similar result, but with a smaller coefficient, can be obtained variationally, by assuming that the holes behave as spinless fermions in a ferromagnetic background (broken curve in Fig. 1). These results can be compared with results from an exact diagonalization of a $`4\times 4`$ $`tJ`$ model with $`J=0`$ (crosses in Fig. 1). The figure illustrates the approximate $`x(1x)`$ dependence of the kinetic energy. In the following we use the $`tJ`$ model with $`J/t=0.3`$ and $`t=0.4`$ eV. We solve the model, using the finite $`T`$ exact diagonalization technique of Jaklic and Prelovsek. We consider a finite cluster of $`4\times 4`$ atoms and use the grand canonical ensemble. The resulting numerical results for the kinetic energy for small values of $`T`$ and $`x`$ can be approximately expressed in the form of Eq. (8), but with the prefactor 4 reduced to about 3.4. For the values of $`T`$ of interest here, the $`T`$ dependence of $`T_K`$ is weak, and it is neglected in the following. If the system has a strong electron-phonon interaction, the kinetic energy should be somewhat further reduced.

Combining Eqs. (2, 8), we predict that the right hand-side of Eq. (2) should be proportional to $`x`$, for small $`x`$. The effective number of carriers per $`\mathrm{Cu}`$ atom $`N_{eff}(x)`$ can be estimated integrating the optical conductivity up to $`\omega _c=1.2`$ eV, as suggested by Yamada et al.. This also includes some contribution from interband transitions, and Yamada et al. therefore subtracted $`N_{eff}(0)`$, which should contain only interband transitions. They found that $`N_{eff}(x)N_{eff}(0)`$ is indeed approximately proportional to $`x`$, up to $`x=0.12`$. Based on the absolute numbers by Uchida et al., we find that our factor of proportionality in Eqs. (2, 8) agrees with experiment to within about 20-30 $`\%`$. Due to the problems of removing the interband transitions, this should be a good agreement.

We consider conduction in the ab-plane of a cuprate and assume that the $`\omega =0`$ peak is gone. Insertion of Eq. (8), but with the prefactor 3.4, in Eqs. (2, 3) gives an approximate upper limit to the resistivity for small and intermediate values of $`T`$

$$\rho =\frac{0.4}{x(1x)}\mathrm{m}\mathrm{\Omega }\mathrm{cm},$$
(9)
where we have used distance $`c=6.4`$ Å between the CuO2 planes. Inserting the calculated kinetic energy for the $`tJ`$ model would lead to a somewhat larger resistivity for small $`x`$ and it would introduce a weak $`T`$ dependence. The saturation resistivity in Eq. (9) is much larger than Eq. (1) for small $`x`$. It illustrates that one should expect saturation at a substantially larger resistivity than predicted by the Ioffe-Regel condition. The saturation resistivity in Eq. (9) is also much larger than the result derived for a model of weakly correlated transition metal compounds

$$\rho \frac{0.2d}{n}\frac{0.5}{n},$$
(10)
where $`d`$ is expressed in Å. This is partly due to the degeneracy being just $`n=1`$ for the $`tJ`$ model but $`n=5`$ for our model of the weakly correlated transition metal compounds. Furthermore, the strong correlation drastically reduces the kinetic energy in the cuprates, which shows up as the factor $`x(1x)`$ in the denominator of Eq. (9). Finally, the large value of the lattice parameter $`c`$ also increases $`\rho `$ for the cuprates.

We now consider the $`T`$ dependent resistivity for different doping. The experimental results of Takagi et al. are shown in Fig. 2. For $`T300`$ K, more recent results have been obtained by Ando et al.. Their resistivities are qualitatively similar but generally smaller than the results in Fig. 2. The results by Takagi et al. are compared with the saturation resistivity in Eq. (9). For all values of $`x`$, the experimental resistivity is smaller than the saturation resistivity predicted here. The experimental results do therefore not demonstrate absence of saturation. Actually, for $`x=0.04`$ and $`x=0.07`$ the data show signs of saturation for values not much smaller than the expected saturation resistivities, but at much larger values than the saturation resistivity predicted by the Ioffe-Regel criterion (Eq. (1)) (indicated by an arrow). This shows that the Ioffe-Regel criterion is invalid for these systems. It is interesting to plot $`x\rho (T)`$. In such a plot the saturation resistivity becomes almost independent of $`x`$. Indeed the curves of $`x\rho (T)`$ for $`x=0.04`$ and $`x=0.07`$ fall almost on top of each other, suggesting that close to saturation there is a scaling of $`\rho (T)`$ by $`1/x`$, as implied by the arguments above. For $`x=0.15`$ and $`x=0.34`$ one would have to study much higher (and experimentally unaccessible) values of $`T`$ to determine whether or not there is saturation.

La2-xSrxCuO4 has a structural transition from an orthorhombic to a tetragonal phase. This transition happens, however, at a somewhat lower temperature than the break in the slopes of the curves for $`x=0.04`$ and $`x=0.07`$ in Fig. 2.

We have also considered the resistivity of YBa2Cu3O6+δ, Nd2-xCexCuO4-y, Bi2Sr2Ca1-xYxCu2O8+y, Tl2Ba2CuO6+δ and HgBa2Can-1CunO2n+2+x. In all these cases we find that the experimental resistivity does not exceed the expected saturation resistivity. In the case of Bi2Sr2Ca1-xYxCu2O8+y saturation has been reported. This provides further support for our theory.

Above we have provided an upper limit to the resistivity based on the f-sum rule and the calculated kinetic energy and without specifying a scattering mechanism. We have also calculated the resistivity of the $`tJ`$ model directly, although the results are not very accurate for the small clusters considered here. For $`x=0.15`$, the calculated resistivity is of the right order of magnitude, while for $`x=0.07`$ and $`x=0.04`$ it is substantially smaller than the experimental results in Fig. 2. This suggests that there is an additional scattering mechanism, beyond the electron-electron scattering included explicitly in the $`tJ`$ model, which becomes important for $`x=0.04`$ and $`x=0.07`$. This additional mechanism apparently makes the resistivity increase so rapidly for small $`T`$ that the approximate upper limit (9) is approached and the resistivity shows sign of saturation. The estimate of the upper limit (9) is nevertheless correct, unless the additional scattering mechanism influences the kinetic energy in a substantial way.

Above we have discussed the situation when $`T`$ is much smaller than the band width, which applies for values of $`T`$ which can be reached experimentally. The kinetic energy $`T_K`$ then has a weak $`T`$ dependence, leading to an upper limit for the resistivity which is only weakly $`T`$ dependent. As can be seen from Eq. (7), however, the $`T`$ dependence becomes strong for very large $`T`$ ($`>W/4`$). Inserting Eq. (7) into Eqs. (2,3) then leads to an upper limit for $`\rho (T)`$ which grows linearly with $`T`$. A similar result would be found even for noninteracting electrons scattered by impurities, i.e., for a $`T`$-independent mechanism, where the $`T`$ dependence can then be traced back to the $`T`$ dependence of the Fermi functions. This illustrates that “saturation” does not necessarily imply that the resistivity becomes $`T`$ independent, just that its slope is reduced substantially at large $`T`$. The effects discussed in this paragraph, however, refer to very large values of $`T`$, of the order of 10000 K. They are therefore only of academic interest, unless the band width is very small.

To summarize, we have shown that the experimental resistivity data for a number of cuprates are consistent with resistivity saturation. The saturation resistivity is, however, much larger than for, e.g., A15 compounds, due to the strong suppression of the kinetic energy in the strongly correlated cuprates. The Ioffe-Regel condition is therefore violated for a many of these systems.

We would like to thank P. Horsch and B. Keimer for useful discussions, and the Max-Planck-Forschungspreis for financial support.


text:
# Random matrix theory, the exceptional Lie groups, and 𝐿-functions.

## 1 Introduction

Most work on the connection between random matrix theory and $`L`$-functions has concentrated on random matrices chosen from ensembles related to the classical compact groups. Montgomery , Rudnick and Sarnak , and Bogomolny and Keating  calculated the correlation functions of the zeros of the Riemann zeta-function, scaled to have unit mean spacing, in the limit as $`T`$, the extent of the averaging range up the critical line, tends to infinity. Their results suggest that these correlation functions coincide with ones relating to the eigenvalues of unitary matrices in the limit as the matrix size, $`N`$, tends to infinity. In the latter case the average is defined with respect to the uniform invariant (Haar) measure on the unitary group $`U(N)`$; that is, with respect to the Circular Unitary Ensemble (CUE) of Random Matrix Theory (RMT). There is extensive numerical evidence in support of this connection , which is expected to extend to the zeros of any given principal $`L`$-function. Katz and Sarnak  conjectured that the distributions of low-lying zeros in families of $`L`$-functions are the same as those of the eigenvalues of matrices from the various classical compact groups (e.g. the orthogonal group $`O(N)`$ and the symplectic group $`USp(2N)`$, as well as $`U(N)`$), the particular group in question being determined by the symmetry of the family. This is also supported by numerical evidence .

It was suggested by Keating and Snaith  that the value distribution of a given principle $`L`$-function on its critical line coincides, in the limit as $`T\mathrm{}`$, with the value distribution of the characteristic polynomials of random unitary matrices, defined again by an average with respect to Haar measure for $`U(N)`$, in the limit as $`N\mathrm{}`$. The random matrix value distribution was calculated in  by expressing the group average in terms of an integral over the matrix eigenvalues, using a formula for the measure due to Weyl , and then relating the resulting ($`N`$-dimensional) integral to one evaluated by Selberg. This idea was later extended in line with the Katz-Sarnak philosophy to relate the value distribution within a given family of $`L`$-functions at the centre of the critical strip to the value distribution of the characteristic polynomials associated with elements of the appropriate classical compact group in the $`N\mathrm{}`$ limit . Again, the random matrix calculations were performed using Weyl’s integration formula and the Selberg integral. One interesting feature of these calculations is that in all cases the logarithm of the characteristic polynomial, normalized appropriately (by $`\mathrm{log}N`$), satisfies a central limit theorem in the limit $`N\mathrm{}`$. This is in agreement with a theorem of Selberg which states that the logarithm of the Riemann zeta function, normalized appropriately (by $`\mathrm{log}\mathrm{log}T`$), also satisfies a central limit theorem in the limit $`T\mathrm{}`$. For further related developments see .

It is in this context that we now ask whether there is a connection between $`L`$-functions and random matrices from the non-classical groups. A particularly interesting class of these groups, closely related to the classical groups, is that of the exceptional Lie groups. Our purpose in this note is to point out that a number of the key constructions which serve to provide the link with random matrix theory in the classical case have analogues for the exceptional Lie groups. We illustrate this by computing the moments and value distribution of the characteristic polynomials of matrices associated with the 7- and 14-dimensional representations of one particular exceptional Lie group, $`G_2`$. The methods employed are again the appropriate Weyl integration formula and generalizations of the Selberg integral conjectured by Macdonald (and known as Macdonald constant term identities) , proven for $`G_2`$ by Zeilberger  and Habsieger  and by Opdam  in the general case. These methods extend to all of the other exceptional Lie groups. We then go on to describe a one-parameter family of $`L`$-functions over a function field, whose value distribution coincides with that of the characteristic polynomials associated with the 7-dimensional representation of $`G_2`$ in the limit as the size of the finite field grows (this was proved by Katz ). The link with finite fields is natural, because $`N`$ is fixed for the exceptional groups and the $`L`$-functions in question (whose zeros correspond to eigenvalues) are polynomials.

This note is organized as follows. In Section 2 we review properties of $`G_2`$ necessary for our random matrix calculations. These calculations are performed in Section 3. The $`L`$-functions associated with $`G_2`$ are constructed in Section 4. In Section 5 we conclude with a brief discussion of the generalization to the other exceptional Lie groups.

## 2 Preliminaries about $`G_2`$

As pointed out in the Introduction, the exceptional Lie groups are closely related to the classical matrix groups. One particularly natural way of seeing this relationship is via their Lie algebras. An important class of Lie algebras, because they form building blocks of more general algebras, is the class of complex semi-simple Lie algebras. This class allows a complete categorization and is elegantly summarized in the possible Dynkin diagrams which encapsulate the allowed root systems; the root systems describe the structure constants of the Lie algebra (standard references for this material include ). The result of the analysis is that the structure of possible root systems is highly constrained. Indeed the only possibilities fall into four infinite families, $`a_n`$, $`b_n`$, $`c_n`$ and $`d_n`$ plus five exceptional cases, $`g_2`$, $`f_4`$, $`e_6`$, $`e_7`$ and $`e_8`$. Each complex Lie algebra has a compact real form and this real form is the Lie algebra of a compact group; the compact real forms of $`a_n`$, $`b_n`$, $`c_n`$ and $`d_n`$ are the Lie algebras of $`SU(n+1)`$, $`SO(2n+1)`$, $`Sp(2n)`$ and $`SO(2n)`$ respectively. The five exceptional cases are the Lie algebras of the groups $`G_2`$, $`F_4`$, $`E_6`$, $`E_7`$ and $`E_8`$. What is significant for our purposes is the many of the constructions which exist for the classical groups have analogues for the exceptional groups.

Here, for concreteness, we will mainly focus on the smallest exceptional case, namely $`G_2`$, for which the Lie algebra is the compact real form of $`g_2`$. The group $`G_2`$ is the automorphism group of the octonions, and has an embedding into $`SO(7)`$.

The group is 14-dimensional and has rank 2. The six positive roots may be taken to be

$`\alpha _1`$ $`=`$ $`\left(\begin{array}{c}1\\ 0\end{array}\right);\alpha _2=\left(\begin{array}{c}3/2\\ \sqrt{3}/2\end{array}\right);`$ (5)
$`\alpha _3`$ $`=`$ $`\alpha _1+\alpha _2=\left(\begin{array}{c}1/2\\ \sqrt{3}/2\end{array}\right);\alpha _4=2\alpha _1+\alpha _2=\left(\begin{array}{c}1/2\\ \sqrt{3}/2\end{array}\right);`$ (10)
$`\alpha _5`$ $`=`$ $`3\alpha _1+\alpha _2=\left(\begin{array}{c}3/2\\ \sqrt{3}/2\end{array}\right);\alpha _6=3\alpha _1+2\alpha _2=\left(\begin{array}{c}0\\ \sqrt{3}\end{array}\right).`$ (15)
The complete set of roots is $`R=\{\pm \alpha _i\},i=1\mathrm{}6`$. The set of short roots is $`R_S=\{\pm \alpha _1,\pm \alpha _3,\pm \alpha _4\}`$; the set of long roots is $`R_L=\{\pm \alpha _2,\pm \alpha _5,\pm \alpha _6\}`$. The Weyl group $`W`$ of $`G_2`$ is a dihedral group with $`12`$ elements.

Let $`T`$ be a maximal torus of $`G_2`$, which is isomorphic to a product of two circles $`S^1\times S^1`$. Every element of $`G_2`$ is conjugate in $`G_2`$ to an element of $`T`$, which is unique up to conjugation by the Weyl group.

Weyl’s integration formula reads: If $`d\mu _{\mathrm{inv}}(g)`$ is the Haar probability measure on $`G_2`$, $`dt`$ is the Haar probability measure on $`T`$ and $`F`$ is a continuous function on $`G_2`$, invariant under conjugation, then

$`{\displaystyle _{G_2}}F(g)𝑑\mu _{\mathrm{inv}}(g)={\displaystyle \frac{1}{12}}{\displaystyle _T}F(t)|\mathrm{\Delta }(t)|^2𝑑t`$ (16)
where

$`\mathrm{\Delta }(t)={\displaystyle \underset{\sigma W}{}}(det\sigma )t^{\sigma (\delta )}=t^\delta {\displaystyle \underset{\alpha >0}{}}(1t^\alpha )`$ (17)
(the equality is Weyl’s denominator formula), where $`\delta `$ is half the sum of the positive roots:

$`\delta ={\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha >0}{}}\alpha =5\alpha _1+3\alpha _2.`$ (18)
and $`\sigma (\delta )`$ means the vector obtained from $`\delta `$ by the Weyl group element $`\sigma `$. If we parametrize a particular maximal torus by $`t=(t_1,t_2)`$, then $`t^\alpha `$ is an expression of the form $`t_1^{\alpha .e_1}t_2^{\alpha .e_2}`$ where $`e_1`$ and $`e_2`$ are (two-component) vectors; the vectors $`e_1`$ and $`e_2`$ and range of the parameters $`t_1`$ and $`t_2`$ depend on the choice of maximal torus. We derive our results below without needing to make an explicit choice for this torus.

For any pair of integers $`[n_1,n_2]`$ there is an irreducible representation $`\rho _{[n_1,n_2]}`$ which has highest weight $`\lambda _{[n_1,n_2]}=n_1\omega _1+n_2\omega _2`$ where $`\omega _1=\alpha _4`$ and $`\omega _2=\alpha _6`$ are the fundamental weights.

The character $`\chi _\lambda `$ of the representation evaluated at the group element $`tT`$ is

$`\chi _\lambda (t)=\mathrm{Tr}\left[\rho _\lambda (t)\right]={\displaystyle \underset{\mu }{}}d_\mu t^\mu ,`$ (19)
where $`\mu `$ are the weights of the representation, $`d_\mu `$ is the multiplicity of the weight $`\mu `$.

The character $`\chi _\lambda (t)`$ and dimension $`d_\lambda `$ of the representation $`\rho _\lambda `$ are given by Weyl’s formulae:

$`\chi _\lambda (t)`$ $`=`$ $`{\displaystyle \frac{\underset{\sigma W}{}(det\sigma )t^{\sigma (\lambda +\delta )}}{_{\sigma W}(det\sigma )t^{\sigma (\delta )}}};`$ (20)
$`d_\lambda `$ $`=`$ $`{\displaystyle \frac{\underset{\alpha >0}{}(\lambda +\delta ).\alpha }{_{\alpha >0}\delta .\alpha }}.`$ (21)
The sums in (20) are over elements $`\sigma `$ in the Weyl group. With these conventions, the orthogonality relation for characters is

$`{\displaystyle \frac{1}{12}}{\displaystyle _T}|\mathrm{\Delta }(t)|^2\chi _{[n_1,n_2]}(t)\overline{\chi _{[m_1,m_2]}(t)}𝑑t=\delta _{n_1,m_1}\delta _{n_2,m_2}.`$ (22)
We will be particularly interested in the fundamental representations $`[1,0]`$ (induced from the embedding of $`G_2`$ into $`SO(7)`$), and $`[0,1]`$ (the adjoint representation), which have characters

$`\chi _{[1,0]}(t)`$ $`=`$ $`1+{\displaystyle \underset{\alpha R_S}{}}t^\alpha ;`$ (23)
$`\chi _{[0,1]}(t)`$ $`=`$ $`2+{\displaystyle \underset{\alpha R}{}}t^\alpha ;`$ (24)
and dimensions

$`d[1,0]=7;d[0,1]=14.`$ (25)
## 3 Characteristic polynomials

We will focus on the characteristic polynomials

$`Z(U_\rho ,\theta ):=det(1U_\rho e^{i\theta })`$ (26)
of (unitary) matrices $`U_\rho `$ coming from a given representation $`\rho `$ of the group. The group elements which these matrices represent can be thought of as being chosen randomly from the group with respect to the uniform invariant (Haar) measure. We will calculate explicit expressions for the averages (over the group, with respect to Haar measure) of $`|Z|^s`$ for complex numbers $`s`$ (see  for analogous calculation relating to $`U(N)`$, $`O(N)`$, and $`USp(2N)`$).

First let us consider the modulus of $`Z`$. We wish to calculate

$`<|Z(U_\rho ,\theta )|^s>_{G_2}`$ $`=`$ $`{\displaystyle |Z(U_\rho (g),\theta )|^s𝑑\mu _{\mathrm{inv}}(g)}`$ (27)
$`=`$ $`{\displaystyle |det(1U_\rho (g)e^{i\theta })|^s𝑑\mu _{\mathrm{inv}}(g)}.`$
Since the integrand is a class function, this integral reduces to an integral over the maximal torus $`T`$:

$`<|Z(U_\rho ,\theta )|^s>_{G_2}={\displaystyle \frac{1}{12}}{\displaystyle _T}|\mathrm{\Delta }(t)|^2|det(1U_\rho (t)e^{i\theta })|^s𝑑t.`$ (28)
### 3.1 The seven-dimensional representation

This representation is induced by the embedding of $`G_2`$ as a subgroup of $`SO(7)`$. From (23), we can calculate that

$`Z(U_{[1,0]},\theta )=det(1U_{[1,0]}(t)e^{i\theta })=(1e^{i\theta }){\displaystyle \underset{\alpha R_S}{}}(1t^\alpha e^{i\theta }).`$ (29)
We see that $`Z`$ has a zero at $`\theta =0`$, as is the case for the characteristic polynomial of any odd-dimensional orthogonal matrix. Let us define $`\widehat{Z}`$ as

$`\widehat{Z}(U_{[1,0]},\theta )=(1e^{i\theta })^1Z(U_{[1,0]},\theta ).`$ (30)
We now present a formula for $`<|\widehat{Z}(U_{[1,0]},\theta )|^s>`$ at $`\theta =0`$.

We have that

$`<|\widehat{Z}(U_{[1,0]},0)|^s>_{G_2}`$ $`=`$ $`{\displaystyle \frac{1}{12}}{\displaystyle _T}|\mathrm{\Delta }(t)|^2|{\displaystyle \underset{\alpha R_S}{}}(1t^\alpha )|^s𝑑t`$ (31)
$`=`$ $`{\displaystyle \frac{1}{12}}{\displaystyle _T}|\mathrm{\Delta }(t)|^2{\displaystyle \underset{\alpha R_S}{}}(1t^\alpha )^sdt,`$
but

$`|\mathrm{\Delta }(t)|^2={\displaystyle \underset{\alpha R}{}}(1t^\alpha ),`$ (32)
and so

$`<|\widehat{Z}(U_{[1,0]},0)|^s>_{G_2}={\displaystyle \frac{1}{12}}{\displaystyle _T}{\displaystyle \underset{\alpha R}{}}(1t^\alpha )^{k_\alpha }dt,`$ (33)
where

$`k_\alpha =\{\begin{array}{cc}s+1\hfill & \text{if }\alpha R_S\text{,}\hfill \\ 1\hfill & \text{if }\alpha R_L\hfill \end{array}.`$ (34)
Consider first the case that $`s`$ is an integer. The value of the integral (33) is then the constant term in the expression

$`{\displaystyle \frac{1}{12}}{\displaystyle \underset{\alpha R}{}}(1t^\alpha )^{k_\alpha }.`$ (35)
The value of this constant term is in turn one of Macdonald’s celebrated constant term conjectures , proved for $`G_2`$ by Zeilberger  and Habsieger  (see Opdam  for a uniform proof):

$`{\displaystyle \frac{(3k_S+3k_L)!(2k_S)!(2k_L)!(3k_L)!}{12(2k_S+3k_L)!(k_S+2k_L)!(k_S+k_L)!(k_S)!((k_L)!)^2}},`$ (36)
where $`k_S`$ (resp. $`k_L`$) is the value of $`k_\alpha `$ for the short (resp. long) roots.

Thus for the representation $`[1,0]`$, $`k_S=s+1`$ and $`k_L=1`$, and so for $`s`$ a positive integer or zero,

$`<|\widehat{Z}(U_{[1,0]},0)|^s>_{G_2}={\displaystyle \frac{(3s+6)!(2s+2)!}{(2s+5)!(s+3)!(s+2)!(s+1)!}}.`$ (37)
It follows from Carlson’s theorem (see ) that then

$`<|\widehat{Z}(U_{[1,0]},0)|^s>_{G_2}={\displaystyle \frac{\mathrm{\Gamma }(3s+7)\mathrm{\Gamma }(2s+3)}{\mathrm{\Gamma }(2s+6)\mathrm{\Gamma }(s+4)\mathrm{\Gamma }(s+3)\mathrm{\Gamma }(s+2)}}`$ (38)
for $`\mathrm{Re}s>3/2`$. To see this, note that one may deduce directly from (33) that $`2^{6s}<|\widehat{Z}(U_{[1,0]},0)|^s>_{G_2}`$ is bounded when $`\mathrm{Re}s>0`$, and from Stirling’s formula that

$$2^{6s}\frac{\mathrm{\Gamma }(3s+7)\mathrm{\Gamma }(2s+3)}{\mathrm{\Gamma }(2s+6)\mathrm{\Gamma }(s+4)\mathrm{\Gamma }(s+3)\mathrm{\Gamma }(s+2)}$$
(39)
is also bounded in the same half-plane. The function

$$2^{6s}\left(<|\widehat{Z}(U_{[1,0]},0)|^s>_{G_2}\frac{\mathrm{\Gamma }(3s+7)\mathrm{\Gamma }(2s+3)}{\mathrm{\Gamma }(2s+6)\mathrm{\Gamma }(s+4)\mathrm{\Gamma }(s+3)\mathrm{\Gamma }(s+2)}\right)$$
(40)
is therefore regular and bounded in $`\mathrm{Re}s>0`$, and vanishes when $`s`$ is a non-negative integer. Carlson’s theorem therefore implies (38).

Now that (38) has been established, we may immediately write down expressions for the probability density functions associated with the value distributions of $`\mathrm{log}|\widehat{Z}(U_{[1,0]},0)|`$,

$`P_1(x)={\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{\Gamma }(3iy+7)\mathrm{\Gamma }(2iy+3)}{\mathrm{\Gamma }(2iy+6)\mathrm{\Gamma }(iy+4)\mathrm{\Gamma }(iy+3)\mathrm{\Gamma }(iy+2)}}\mathrm{e}^{iyx}𝑑y`$ (41)
and $`|\widehat{Z}(U_{[1,0]},0)|`$,

$`P_2(x)={\displaystyle \frac{1}{2\pi ix}}{\displaystyle _{ci\mathrm{}}^{c+i\mathrm{}}}{\displaystyle \frac{\mathrm{\Gamma }(3s+7)\mathrm{\Gamma }(2s+3)}{\mathrm{\Gamma }(2s+6)\mathrm{\Gamma }(s+4)\mathrm{\Gamma }(s+3)\mathrm{\Gamma }(s+2)}}x^s𝑑s`$ (42)
for any $`c>0`$. One can easily deduce asymptotic properties of the probability density functions from these integrals; see, for example .

We note finally that $`\widehat{Z}`$ is real and positive at $`\theta =0`$ and so its phase there is zero.

### 3.2 The fourteen-dimensional representation

In this case the determinant $`Z`$ has a double zero at $`\theta =0`$ (corresponding to the twice repeated weight 0, see (24)). Thus we define

$`\widehat{Z}(U_{[0,1]},\theta )=(1e^{i\theta })^2det(1U_{[0,1]}e^{i\theta }).`$ (43)
Similar calculations to the previous case show that for $`s`$ a positive integer or zero, $`<|\widehat{Z}(U_{[0,1]},0)|^s>_{G_2}`$ is given by the constant term in

$`{\displaystyle \frac{1}{12}}{\displaystyle \underset{\alpha R}{}}(1t^\alpha )^{k_\alpha },`$ (44)
where in this case,

$`k_S=k_L=s+1.`$ (45)
Thus for $`s`$ a non-negative integer the Macdonald identity quoted above implies that

$`<|\widehat{Z}(U_{[0,1]},0)|^s>_{G_2}={\displaystyle \frac{(6s+6)!(2s+2)!}{12(5s+5)!((s+1)!)^3}}.`$ (46)
Once again, Carlson’s theorem may be applied (in this case after multiplication by $`2^{12s}`$) to show that

$`<|\widehat{Z}(U_{[0,1]},0)|^s>_{G_2}={\displaystyle \frac{\mathrm{\Gamma }(6s+7)\mathrm{\Gamma }(2s+3)}{12\mathrm{\Gamma }(5s+6)\left[\mathrm{\Gamma }(s+2)\right]^3}}`$ (47)
for $`\mathrm{Re}s>7/6`$. This can then be used to write down expressions for the probability density functions associated with the value distributions of $`\mathrm{log}|\widehat{Z}(U_{[0,1]},0)|`$ and $`|\widehat{Z}(U_{[0,1]},0)|`$, as in the previous section.

Also, as in the case of the representation $`[1,0]`$, the phase of $`\widehat{Z}`$ is zero.

### 3.3 Other representations and $`\theta 0`$

For other representations of $`G_2`$ and for $`\theta 0`$, the integrand is not of the form

$`{\displaystyle \underset{\alpha R}{}}(1t^\alpha )^{k_\alpha },`$ (48)
and, as far as we are aware, closed form expressions for these integrals are not known.

## 4 Value distribution problems over function fields

Our purpose now is to outline the analogy between number fields and function fields over a finite field in the context of the value distribution of zeta- and $`L`$-functions for these cases and the predicted behaviour in terms of Random Matrix Theory.

It was proved by Selberg that the logarithm of the Riemann zeta function on the critical line has a Gaussian value distribution , and the same is true for all $`L`$-functions  under suitable assumptions. Selberg also investigated the “$`q`$-analogue” of these results for the value distribution of the family of Dirichlet $`L`$-functions  at a point on the critical line and there too obtained a Gaussian value distribution. Precisely, for $`q`$ prime we have $`q2`$ primitive characters $`\chi `$ modulo $`q`$, and for fixed $`t`$ consider the $`q2`$ numbers

$$\frac{argL(\frac{1}{2}+it,\chi )}{\sqrt{\frac{1}{2}\mathrm{log}\mathrm{log}q}}$$

($`\chi `$ varies over all primitive/nonprincipal characters modulo $`q`$). Then as $`q\mathrm{}`$, these numbers are distributed as a standard Gaussian.

Our purpose in this section is two-fold: first it is to point out that there are corresponding results for various families of $`L`$-functions over function fields, with the rôle of the Gaussian being replaced by various distributions from Random Matrix Theory; and second to construct an example for which the appropriate random-matrix distribution is the $`G_2`$ result calculated above.

Likewise, there are similar results for the moments of the $`L`$-functions, which in the number field setting are mostly conjectural , but in the function field setting can sometimes be proved.

### 4.1 Zeta functions

Let $`k`$ be a finite field of cardinality $`q`$ and $`X/k`$ a (smooth, geometrically connected, proper) curve defined over $`k`$. The zeta function of $`X/k`$ is given by the series

$$Z(X;T)=\mathrm{exp}(\underset{n=1}{\overset{\mathrm{}}{}}N_n\frac{T^n}{n}),$$

where $`N_n=\mathrm{\#}X(k_n)`$ is the number of points of $`X`$ over the field $`k_n`$, the extension of $`k`$ of degree $`n`$. The series is absolutely convergent for $`|T|<1`$.

Trivial example: take $`X=𝐏^1`$, the projective line. The number of points of $`𝐏^1`$ over the finite field $`k_n`$ is $`\mathrm{\#}k_n+1=q^n+1`$ and so

$$Z(𝐏^1;T)=\frac{1}{(1qT)(1T)}.$$

This zeta function has an Euler product

$$Z(X;T)=\underset{p}{}(1T^{\mathrm{deg}p})^1,|T|<1,$$

where $`p`$ runs over all closed points of $`X`$. (In the example of $`𝐏^1`$, the closed points $`p`$ correspond to irreducible monic polynomials $`p(x)k[x]`$ with the addition of the “point at infinity”.)

It turns out that $`Z(X;T)`$ is a rational function of $`T`$, of the form

$$Z(X;T)=\frac{P(X;T)}{(1T)(1qT)}$$

with $`P(X;T)1+T𝐙[T]`$ a monic integer polynomial of degree $`2g`$, $`g`$ being the genus of the curve $`X`$, which we can write as $`P(X;T)=_{j=1}^{2g}(1\alpha _jT)`$. The inverse roots $`\alpha _j`$ are thus algebraic integers. Further, there is a functional equation $`Tq/T`$:

$$Z(X;\frac{1}{qT})=q^{1g}T^{22g}Z(X;T).$$

If we set $`T=q^s`$ then the functional equation translates into $`s1s`$.

The “Riemann Hypothesis for curves over a finite field” (proved in the general case by A. Weil) is that all the inverse roots $`\alpha _j`$ of $`P(X;T)`$ have absolute value $`\sqrt{q}`$, that is as a function of the variable $`s`$ all zeros are on the line $`\mathrm{Re}s=1/2`$.

What is especially important for our purpose is that the polynomial $`P(X;T)`$ is the characteristic polynomial of a matrix: there is a unique conjugacy class $`\mathrm{\Theta }_XUSp(2g)`$ in the unitary symplectic group such that $`P(X;T)=det(Iq^{1/2}T\mathrm{\Theta }_X)`$.

### 4.2 Families of curves

Now consider a “family” of curves $`X/k`$. In order to study the behaviour of $`P(X;T)`$ as $`X`$ varies, it suffices to understand the distribution of the conjugacy classes $`\mathrm{\Theta }_X`$. In several cases, it is known that as $`q\mathrm{}`$ these become equidistributed in $`USp(2g)`$ (with respect to Haar measure).

For instance, this is the case for the family $`_g`$ of all $`k`$-isomorphism classes of (smooth, geometrically connected, proper) curves of given genus $`g`$ \[15, Theorem 10.7.15\].

This allows one to compute arithmetic quantities such as the moments of $`P(X;T)`$ as $`X`$ varies in $`_g(k)`$ by using the corresponding (non-arithmetic) computation in Random Matrix Theory for $`USp(2g)`$. Thus one finds that for $`T`$ fixed, say $`q^1T=1`$, one has

$$\underset{q\mathrm{}}{lim}\frac{1}{\mathrm{\#}_g(k)}\underset{X_g(k)}{}P(X,q^{1/2})^s=_{USp(2g)}det(IA)^sd_{Haar}(A).$$

The moments of the characteristic polynomial in $`USp(2g)`$ were computed in  and are given by

$$_{USp(2g)}det(IA)^sd_{Haar}(A)=2^{2gs}\underset{j=1}{\overset{g}{}}\frac{\mathrm{\Gamma }(1+g+j)\mathrm{\Gamma }(1/2+s+j)}{\mathrm{\Gamma }(1/2+j)\mathrm{\Gamma }(1+s+g+j)}.$$

The probability density functions for the value distributions associated with the polynomial and its logarithm may then be written as integrals, as in (3.1) and (42) . In the case of the logarithm of the characteristic polynomial, the limit distribution when $`g\mathrm{}`$ is a Gaussian.

### 4.3 $`L`$-functions attached to exponential sums

We consider one-variable exponential sums constructed as follows: let $`k`$ be a finite field with $`q`$ elements as above, $`f(x)`$ and $`h(x)k(x)`$ be rational functions, $`\psi `$ a nontrivial additive character of $`k`$ (e.g. for $`k=𝐙/p𝐙`$ take $`\psi (x)=\mathrm{exp}(2\pi iax/p)`$, $`0a𝐙/p𝐙`$), and $`\chi `$ a multiplicative character of $`k^\times `$. Set

$$S(\psi ,\chi ;f,h;q)=\underset{x}{}\psi (f(x))\chi (h(x)),$$

the sum running over all $`xk`$ which are not poles of $`f,h`$ and such that $`h(x)0`$. For the finite extension $`k_n`$ of degree $`n`$ of $`k`$, we get nontrivial characters $`\psi _n=\psi \mathrm{Tr}_{k_n/k}`$ and $`\chi _n=\psi \mathrm{N}_{k_n/k}`$ by composing with the trace and norm maps. Correspondingly we get exponential sums for $`k_n`$

$$S_n(\chi ,\psi ;f,h):=S(\chi _n,\psi _n;f,h;q^n).$$

The $`L`$-function is defined as

$$L(S,T)=\mathrm{exp}(\underset{n=1}{\overset{\mathrm{}}{}}S_n(\chi ,\psi ;f,h)\frac{T^n}{n}).$$

These have an Euler product decomposition and are rational functions of $`T`$.

In many cases of interest to us, it turns out that $`L(T)`$ is in fact a polynomial of the form $`det(Iq^{1/2}\mathrm{\Theta }_S)`$ with $`\mathrm{\Theta }_S`$ a unitary matrix.

### 4.4  Gauss sums

Given a nontrivial additive character $`\psi `$ of $`k`$ and a nontrivial multiplicative character $`\chi `$ of $`k^\times `$, one defines the Gauss sum $`g(\chi ,\psi )`$ by

$$g(\chi ,\psi )=\underset{x0}{}\chi (x)\psi (x).$$

Correspondingly we get Gauss sums for $`k_n`$

$$g_n(\chi ,\psi ):=g(\chi _n,\psi _n).$$

To compute the $`L`$-function

$$L(g(\chi ,\psi ),T)=\mathrm{exp}(\underset{n=1}{\overset{\mathrm{}}{}}g_n(\chi ,\psi )\frac{T^n}{n})$$

one can use the Hasse-Davenport relations:

$$g_n(\chi ,\psi )=(g(\chi ,\psi ))^n.$$

These give

$$L(g(\chi ,\psi ),T)=1+Tg(\chi ,\psi ).$$

As is well known, $`|g(\chi ,\psi )|=\sqrt{q}`$. Thus we may write $`g(\chi ,\psi )=\sqrt{q}e^{i\theta _\chi }`$ (we omit the dependence on $`\psi `$ which is of a trivial nature). Setting $`s=1/2+i\theta /\mathrm{log}q`$, $`T=q^s=q^{1/2}e^{i\theta }`$, we find

$$L(g(\chi ,\psi ),T)=1+e^{i(\theta _\chi \theta )}.$$

The $`q2`$ angles $`\{\theta _\chi :\chi \chi _0\}`$ are uniformly distributed in $`[0,2\pi )`$. This is easy to see from Deligne’s estimate on hyper-Kloosterman sums, see \[12, section 1.3.3\]. Thus the moments of $`L(g(\chi ,\psi ),T)`$ and of its logarithm, averaged over $`\chi `$ and taken as $`q\mathrm{}`$, are the same as those for the function $`1+e^{i\theta }`$.

### 4.5 Kloosterman sums (i)

These are the sums

$$\mathrm{Kl}(a,p)=\underset{x_1x_2=amodp}{}\mathrm{exp}\frac{2\pi i}{p}(x_1+x_2).$$

More generally for a finite field $`k`$ with $`q`$ elements, take a nontrivial additive character $`\psi `$ and $`a0`$, and set

$$\mathrm{Kl}(a,q)=\underset{x_1x_2=a}{}\psi (x_1+x_2).$$

This sum is real (replace $`xx`$), and as Weil proved satisfies

$$|\mathrm{Kl}(a,q)|2\sqrt{q}.$$

The associated $`L`$-function is a polynomial of degree $`2`$

$$L(\mathrm{Kl}(a,q),T)=1+\mathrm{Kl}(a,q)T+qT^2.$$

It is of the form $`det(Iq^{1/2}\mathrm{\Theta }_a)`$ with $`\mathrm{\Theta }_aSU(2)`$. It was shown by Katz  that as $`q\mathrm{}`$, the $`q1`$ conjugacy classes $`\{\mathrm{\Theta }_a:ak^\times \}`$ become equidistributed in $`SU(2)`$ with respect to Haar measure. This implies we can compute the value distribution of $`L`$ and $`\mathrm{log}L`$ via RMT on $`SU(2)`$.

### 4.6 Kloosterman sums (ii)

We next look at an example of exponential sums in several variables: hyper-Kloosterman sums are $`n`$-variable sums ($`n2`$) generalizing the previous example, given by

$$\mathrm{Kl}_n(a,q)=\underset{x_1x_2\mathrm{}x_n=a}{}\psi (x_1+x_2+\mathrm{}+x_n).$$

Replacing $`\psi `$ by $`\psi \mathrm{Tr}_{k_m/k}`$ gives the sum $`\mathrm{Kl}_n(a,q^m)`$.

The associated $`L`$-function is defined as

$$L_a(T):=\mathrm{exp}\left((1)^n\underset{m=1}{\overset{\mathrm{}}{}}\mathrm{Kl}_n(a,q^m)\frac{T^m}{m}\right).$$

This $`L`$-function is a polynomial of degree $`n`$. It was proved by Katz  that it can be written as $`det(Iq^{(n1)/2}T\mathrm{\Theta }_n(a,q))`$, with $`\mathrm{\Theta }_n(a,q)K_n`$, where $`K_n`$ is the compact group $`USp(2n)`$, $`n`$ even, $`SU(n)`$, $`n`$, $`q`$ odd, $`SO(n)`$, $`q`$ even, $`n7`$ odd, and $`G_2`$, $`q`$ even , $`n=7`$. Moreover, as $`a`$ varies through all nonzero elements of $`k`$, the $`q1`$ conjugacy classes $`\mathrm{\Theta }_n(a,q)`$ of $`K_n`$ become equidistributed there as $`q\mathrm{}`$ while keeping the type of $`K_n`$ fixed. For instance, taking $`n=7`$ and $`q=2^r`$, $`r\mathrm{}`$ gives $`2^r1`$ conjugacy classes $`\{\mathrm{\Theta }_7(a,2^r):0a𝐅_{2^r}\}`$ which become equidistributed in $`G_2`$ as $`r\mathrm{}`$.

### 4.7 An exponential sum associated to $`G_2`$

Let $`p`$ be a prime, $`p17`$, $`k=𝐙/p𝐙`$, $`\chi _{(2)}`$ the unique quadratic character (Legendre symbol) of $`k^\times `$, and $`\psi `$ a nontrivial additive character of $`k`$, that is $`\psi (x)=e^{2\pi iax/p}`$ for some $`ak^\times `$. Consider for $`tk^\times `$ the exponential sum

$$𝒦𝒯(t)=\underset{xk^\times }{}\chi _{(2)}(x)\psi (x^7+tx).$$

These sums were studied by Nick Katz and the results below are due to him .

Note that $`\overline{𝒦𝒯(t)}=\chi _{(2)}(1)𝒦𝒯(t)`$ and so $`𝒦𝒯(t)`$ is real if $`\chi _{(2)}(1)=1`$, that is if $`p=1mod4`$, and imaginary if $`\chi _{(2)}(1)=1`$, i.e. if $`p=3mod4`$.

In view of the transformation properties under complex conjugation, we divide the exponential sum $`𝒦𝒯(t)`$ by the quadratic Gauss sum $`g(\chi _{(2)})`$ to get a real number. Furthermore, there is a (unique) choice of sign $`ϵ_p=\pm 1`$ so that<sup>1</sup><sup>1</sup>1At the time of writing we do not know how to determine $`ϵ_p`$.

$$𝒦𝒯^{}(t)=ϵ_p\frac{𝒦𝒯(t)}{g(\chi _{(2)})}$$

is minus the trace of a matrix $`\mathrm{\Theta }_tSO(7)`$: $`𝒦𝒯(t)=\mathrm{Tr}\mathrm{\Theta }_t`$. Moreover, this matrix turns out to lie in $`G_2`$.

The associated $`L`$-function is a polynomial of degree $`7`$, which is a characteristic polynomial of the element $`\mathrm{\Theta }_t`$ of $`G_2`$,

$$L(𝒦𝒯^{}(t),T)=det(I\mathrm{\Theta }_tT).$$

As $`t`$ varies in $`(𝐙/p/Z)^\times `$, these $`p1`$ conjugacy classes $`\mathrm{\Theta }_t`$ become equidistributed in $`G_2`$ as $`p\mathrm{}`$. Thus the value distribution of $`L(𝒦𝒯^{}(t),T)`$ at fixed $`T`$ is computed by RMT for $`G_2`$.

## 5 Other Lie groups

The random matrix calculations reported here for $`G_2`$ generalize straightforwardly to the other exceptional Lie groups. In each case one has a Weyl integration formula, which allows the moments of the characteristic polynomials associated with representations of the group to be written as integrals over the Cartan subgroup, and a Macdonald identity, which enables the integrals to be evaluated as ratios of $`\mathrm{\Gamma }`$-functions. This prompts the question as to whether families of finite-field $`L`$-functions can be constructed whose value distributions are given by each of the other exceptional Lie groups.

As a final remark, we note that in  Macdonald gives constant term formulae for affine root systems (which are related to Kac-Moody algebras). This suggests the intriguing possibility of extending the ideas described in this paper to families of random matrices arising from the representations of the associated infinite dimensional groups.

## Acknowledgments

We are grateful to Nick Katz for a very helpful correspondence. JPK is also grateful to the American Institute of Mathematics for generous hospitality during the completion of this manuscript.