text:
# Quintessence and Spontaneous Leptogenesis

## Abstract

We propose in this paper a scenario of spontaneous baryogenesis in cosmological models of Quintessence by introducing a derivative coupling of the Quintessence scalar $`Q`$ to the baryon current $`J_B^\mu `$ or the current of the baryon number minus lepton number $`J_{BL}^\mu `$. We find that with a dimension-5 operator $`_\mu QJ_{BL}^\mu `$ suppressed by the Planck mass $`M_{pl}`$ or the Grand Unification Scale $`M_{GUT}`$, baryon number asymmetry $`n_B/s10^{10}`$ can be naturally explained via leptogenesis. We study also the isocurvature baryon number fluctuation generated in our model.

PACS numbers: 98.80.Cq

Evidence is increasing that Universe is spatially flat and accelerating at present time. The simplest account of the cosmic acceleration seems to be a remnant small cosmological constant $`\mathrm{\Lambda }`$, however many physicists are attracted by the idea that a new form of matter is causing the cosmic acceleration. This new form of matter is not clustered gravitationally on the scale of the galaxy clusters and has been usually called dark energy or Quintessence.

For Quintessence, typically one assumes the existence of a scalar field $`Q`$ (or multi scalar fields) which is taken to be homogeneous in space and via its kinetic and potential energy contributions to the energy-momentum tensor tuned to provide an equation of state leading to accelerated expansion, beginning to dominate the matter content of the Universe today. Models of Quintessence can have interesting cosmological properties including tracking behavior.

Generally minimal or non-minimal couplings of Quintessence with gravitation are considered in the literature. Being a dynamical field, Quintessence is expected to have interactions with the ordinary matters, however as argued in Refs. for most of cases the couplings are strongly constrained. But there are exceptions. For example, Carroll has considered an interaction of form $`QF_{\mu \nu }\stackrel{~}{F}^{\mu \nu }`$ with $`F_{\mu \nu }`$ being the electromagnetic field strength tensor and studied its implication on the rotation of the plane of polarization of light coming from distant sources. In this paper we introduce a type of interaction of Quintessence with the matter, which in terms of an effective lagrangian is given by

$`_{eff}={\displaystyle \frac{c}{M}}_\mu QJ^\mu ,`$ (1)
where $`M`$ is the cut-off scale which for example could be the Planck mass $`M_{pl}`$ or the scale of Grand Unification Theory $`M_{GUT}`$, and $`c`$ is the coupling constant which characterizes the strength of Quintessence interacting with the ordinary matter in the Standard Model of the electroweak theory. Specifically we take in this paper $`J^\mu `$ to be the baryon current $`J_B^\mu `$ or the current of baryon number minus lepton number $`J_{BL}^\mu `$, and study their implications on baryogenesis. The lagrangian in Eq.(1) involves derivative and obeys the symmetry $`QQ+\mathrm{𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡}`$ , so the Quintessence potential will not be modified by the quantum corrections.

The mechanism of generating the baryon number asymmetry in this paper follows closely the spontaneous baryogenesis. The term in Eq.(1), when $`\dot{Q}`$ is non-zero during the evolution of spatial flat Friedmann-Robertson-Walker Universe, violates CPT invariance and generates an effective chemical potential $`\mu _b`$ for baryons, $`i.e.`$,

$`{\displaystyle \frac{c}{M}}_\mu QJ_B^\mu c{\displaystyle \frac{\dot{Q}}{M}}n_B=c{\displaystyle \frac{\dot{Q}}{M}}(n_bn_{\overline{b}}),`$ (2)
$`\mu _b=c{\displaystyle \frac{\dot{Q}}{M}}=\mu _{\overline{b}}.`$ (3)
In thermal equilibrium the baryon number asymmetry is given by (when $`Tm_b`$)

$`n_B={\displaystyle \frac{g_bT^3}{6}}({\displaystyle \frac{\mu _b}{T}}+𝒪({\displaystyle \frac{\mu _b}{T}})^3)c{\displaystyle \frac{g_b\dot{Q}T^2}{6M}},`$ (4)
where $`g_b`$ counts the internal degree of freedom of the baryon. Using the familiar expression for entropy density

$$s=\frac{2\pi ^2}{45}g_{}T^3,$$
(5)
we arrive at the final expression for the baryon to entropy ratio

$$n_B/s\frac{15c}{4\pi ^2}\frac{g_b\dot{Q}}{g_{}MT}.$$
(6)
$`\dot{Q}`$ in Eq.(6) can be obtained by solving the equation of motion of Quintessence given below

$$\ddot{Q}+3H\dot{Q}+V^{^{}}(Q)=\frac{c}{M}(\dot{n}_B+3Hn_B),$$
(7)
where H is the Hubble constant and $`V(Q)`$ is the potential of Quintessence field. For the radiation dominated era the Hubble constant is

$$H=\frac{1}{2t}=1.66g_{}^{1/2}\frac{T^2}{M_{pl}}.$$
(8)
The right-handed side of Eq.(7) is about $`\frac{cg_b}{6}\frac{T^2}{M^2}(\ddot{Q}+H\dot{Q})`$ and can be neglected unless in the very early universe when the temperature $`T`$ is close to the cut-off scale.

In this paper we consider a model which has the tracking property, where the potential has a modified exponential form,

$$V(Q)=f(Q)e^{\frac{\lambda }{m_{pl}}Q}.$$
(9)
In Eq.(9), $`f(Q)`$ is a function of Q which is required to change more slowly than $`e^{\frac{\lambda }{m_{pl}}Q}`$ in the regime of tracking. For the detailed discussions on the specific form of $`f(Q)`$ and model properties we refer to Ref.. Here we simply summarize the salient features of this model. For $`\lambda ^2>3(1+w_B)`$ with $`w_B`$ being the ratio of the pressure to the energy density of the background fluid, this model has an attractor solution,

$`\mathrm{\Omega }_Q`$ $``$ $`{\displaystyle \frac{\rho _Q}{\rho _Q+\rho _B}}={\displaystyle \frac{3}{\lambda ^2}}(1+w_B),`$ (10)
$`w_Q`$ $`=`$ $`w_B.`$ (11)
During radiation dominated era, $`\rho _B=\frac{\pi ^2}{30}g_{}T^4`$ and $`w_B=1/3`$, so we have

$`\rho _Q`$ $``$ $`{\displaystyle \frac{2\pi ^2}{15(\lambda ^24)}}g_{}T^4,`$ (12)
$`\dot{Q}`$ $``$ $`{\displaystyle \frac{4\pi }{3\sqrt{10(\lambda ^24)}}}g_{}^{1/2}T^2.`$ (13)
Thus we obtain the ratio of baryon number to entropy, which is given by

$`{\displaystyle \frac{n_B}{s}}{\displaystyle \frac{5c}{\pi \sqrt{10(\lambda ^24)}}}g_bg_{}^{1/2}{\displaystyle \frac{T}{M}}.`$ (14)
Taking $`\lambda ^2>20`$  constrained by nucleosynthesis, and $`g_{}𝒪(100),g_b𝒪(1)`$ we have

$`{\displaystyle \frac{n_B}{s}}0.01c{\displaystyle \frac{T}{M}}.`$ (15)
In the scenario of spontaneous baryogenesis, the baryon asymmetry is generated in thermal equilibrium. This requires baryon number violating interactions occur rapidly ($`\mathrm{\Gamma }_{\mathit{}}>H`$). However if the B-violating interactions keep in equilibrium until $`\dot{Q}0`$, the final baryon asymmetry will be zero. Denoting the epoch when the B-violating interactions freeze out by $`T_D`$, i.e. $`\mathrm{\Gamma }_{\mathit{}}(T_D)H(T_D)`$, the final baryon number asymmetry obtained is

$`{\displaystyle \frac{n_B}{s}}|_{T_D}0.01c{\displaystyle \frac{T_D}{M}}.`$ (16)
For $`T_D`$ around the Grand Unification scale $`10^{16}GeV`$, Eq.(16) shows that it is quite easy to have $`n_B/s10^{10}`$ required by nucleosynthesis and measurement of Cosmic Microwave Background anisotropy, for both $`M=M_{pl}`$ and $`M=M_{GUT}`$. However, if the B-violating interactions of Grand Unified Theory conserve $`BL`$, the asymmetry generated will be erased by the electroweak Sphaleron . In this case Sphaleron $`B+L`$ interactions will make $`T_D`$ as low as around 100 GeV and $`n_B/s`$ generated will be too small unless $`c`$ is unreasonably large.

Now we discuss the mechanism of leptogenesis  in our model by identifying $`J^\mu `$ with $`J_{BL}^\mu `$ in Eq.(1). Doing the calculation with the same procedure as above for $`J^\mu =J_B^\mu `$ we have the final asymmetry of the baryon number minus lepton number

$`{\displaystyle \frac{n_{BL}}{s}}0.01c{\displaystyle \frac{T_D}{M}}.`$ (17)
The asymmetry $`n_{BL}`$ in (17) will be converted to baryon number asymmetry when electroweak Sphaleron $`B+L`$ interaction is in thermal equilibrium which happens for temperature in the range of $`10^2\mathrm{GeV}10^{12}\mathrm{GeV}`$. $`T_D`$ in (17) is the temperature below which the $`BL`$ violating interactions freeze out.

In the Standard Model of electroweak theory, $`BL`$ symmetry is exactly conserved, however it could be violated by the (majorana) neutrino masses via higher dimensional effective operators. For instance, we take a dimension-5 operator,

$`_{\mathit{}}={\displaystyle \frac{2}{f}}l_Ll_L\varphi \varphi +h.c.,`$ (18)
where $`f`$ is a scale of new physics beyond the Standard Model which generates the $`BL`$ violations, $`l_L`$ and $`\varphi `$ are the left-handed lepton and Higgs doublets respectively. When the Higgs field gets vacuum expectation value $`<\varphi >v`$, the left-handed neutrino receives a majorana mass $`m_\nu \frac{v^2}{f}`$.

In the early universe the lepton number violating rate induced by the interaction in (18) is

$`\mathrm{\Gamma }_{\mathit{}}0.04{\displaystyle \frac{T^3}{f^2}}.`$ (19)
This rate at $`T_D`$ is required to be slower than the Universe expansion rate $`1.66g_{}^{1/2}T^2/M_{pl}`$, which will put a lower limit on the scale $`f`$, consequently a upper limit on the neutrino mass $`m_\nu `$ .

To get the upper limit on $`m_\nu `$, we try to minimize the value of $`f`$. This can be done by minimizing the freezing out temperature $`T_D`$. From Eq.(14), one can see that if requiring $`c4\pi `$ so that the effective lagrangian (Eq.(1)) is sensible for the perturbative calculation, $`T_D`$ must be larger than $`10^{10}`$ GeV for $`M=M_{pl}`$, and $`10^7`$GeV for $`M=M_{GUT}10^{16}`$GeV. As a result of it, $`f`$ should be larger than about $`10^{13}`$ GeV and the neutrino mass should be smaller than $`4`$ eV for $`M=M_{pl}`$. In the case that $`M=M_{GUT},f>4.9\times 10^{11}`$GeV and $`m_\nu <0.1`$keV. And for both cases, $`T_D`$ are much smaller than the cut-off scales. This makes our approximation of neglecting the r.h.s of Eq.(6) valid. We have also done a numerical calculation with result which supports for it.

So far what we have considered is the case that Quintessence is homogeneous in space. Since Quintessence is a scalar field one expects its fluctuations in space-time. In the following we study the isocurvature baryon number fluctuations generated in our model.

Since the baryon number asymmetry is proportional to $`\dot{Q}`$ and in the tracking regime $`\dot{Q}\sqrt{V(Q)}`$, we have

$`{\displaystyle \frac{n_B}{s}}\sqrt{V(Q)},`$ (20)
and

$`\delta ({\displaystyle \frac{n_B}{s}}){\displaystyle \frac{V^{^{}}}{2\sqrt{V}}}\delta Q,`$ (21)
$`{\displaystyle \frac{\delta (n_B/s)}{n_B/s}}={\displaystyle \frac{V^{^{}}}{2V}}\delta Q.`$ (22)
To get $`\delta Q`$, we are solving the equation of motion given below,

$`\ddot{\delta Q_k}+3H\dot{\delta Q_k}+{\displaystyle \frac{k^2}{a^2}}\delta Q_k+V{}_{}{}^{\prime \prime }\delta Q_k=0,`$ (23)
where $`\delta Q_k`$ is the Fourier transform of $`\delta Q`$. To have an analytical formula for $`\delta Q`$ we closely follow  by writing $`V^{^{\prime \prime }}`$ as

$`V^{^{\prime \prime }}={\displaystyle \frac{3}{2}}H\dot{c}_Q^2+{\displaystyle \frac{9}{4}}(1c_Q^2)(w_B+c_Q^2+2)H^2,`$ (24)
where $`c_Q^2`$ is the sound speed of Quintessence defined by

$$c_Q^2\frac{\dot{p}_Q}{\dot{\rho }_Q}=w_Q\frac{\dot{w}_Q}{3H(1+w_Q)}.$$
(25)
Consider the case that the sound speed is approximately constant, the first term of the rhs. of Eq.(24) can be dropped out, then the solution to Eq.(23) is:

$$\delta Q_k\eta ^{1/2}[C_1J_\nu (k\eta )+C_2J_\nu (k\eta )],$$
(26)
where $`\eta `$ is the conformal time ($`d\eta \frac{dt}{a}`$), and $`a\eta T^1`$ during the radiation dominated era. In Eq.(26), $`J_{\pm \nu }(k\eta )`$ are the first kind of Bessel functions, $`C_1`$, $`C_2`$ are constants which are time independent, $`\nu ^2`$ and $`\beta `$ are,

$`\nu ^2{\displaystyle \frac{1}{4}}\beta ,`$ (27)
$`\beta {\displaystyle \frac{9}{4}}(1c_Q^2)(w_B+c_Q^2+2).`$ (28)
The modes of fluctuations which are interesting in cosmology are those outside the Hubble radius in very early universe, i.e. $`k\eta 1`$. In such a limit the Bessel function can be approximated by

$$J_\nu (k\eta )\frac{(k\eta )^\nu }{2^\nu \mathrm{\Gamma }(1+\nu )},(k\eta 0).$$
(29)
On the initial value of $`\delta Q_k`$ denoted by $`\delta Q_{ki}`$, we assume the primordial fluctuation of the Quintessence scalar is generated during inflation. Being almost massless $`\delta Q_{ki}`$ is simply

$$|\delta Q_{ki}|=\frac{H_{in}}{\sqrt{2k^3}},$$
(30)
where $`H_{in}`$ is the Hubble constant during inflation. Furthermore the kinetic energy of the Quintessence will be diluted by inflation. Thus the Quintessence energy density remains potential energy dominated after inflation until the time $`\eta _t`$ when it goes into the tracking regime. This period of time is called the potential phase in the Ref.. During this period $`\rho _Q`$ is nearly constant and $`c_Q^2w_B2`$ , hence $`\beta 0`$ and $`V^{^{\prime \prime }}0`$. We can see that $`\delta Q_k`$ is constant in the potential phase. But in the regime of tracking ($`\eta \eta _t`$),

$`c_Q^2`$ $`=`$ $`w_Q,`$ (31)
$`\beta `$ $`=`$ $`{\displaystyle \frac{9}{4}}(1w_Q)(w_B+w_Q+2)={\displaystyle \frac{9}{4}}(1w_Q)({\displaystyle \frac{7}{3}}+w_Q).`$ (32)
Numerically $`1<w_Q1/3`$, so we have $`\beta >1/4`$. We obtain that

$`\nu ^2<0,`$ (33)
$`\nu =i|\nu |=i\sqrt{\beta {\displaystyle \frac{1}{4}}}`$ (34)
$`J_{\pm \nu }(k\eta ){\displaystyle \frac{\mathrm{exp}(\pm i|\nu |\mathrm{ln}\frac{k\eta }{2})}{\mathrm{\Gamma }(1\pm i|\nu |)}}.`$ (35)
Substituting $`J_{\pm \nu }(k\eta )`$ above and $`\delta Q_{ki}`$ (Eq.(30)) into (26), we will get $`\delta Q_k`$.

We have also numerically solved Eq.(23) for the Quintessence model we use in this paper. In Fig.1 we plot $`\delta Q_k`$ as a function of red-shift $`z`$. One can see that $`\delta Q_k`$ is constant in the potential phase i.e. $`\delta Q_{ki}=\delta Q_{kt}`$, then oscillates with amplitude which decays like $`\eta ^{1/2}`$ during tracking. This confirms our analytical analysis above.

Denoting the $`\delta Q_k`$ at the temperature $`T_D`$ when the baryon number asymmetry freezes in by $`\delta Q_{kD}`$, we have the ratio of $`|\delta Q_{kD}|`$ to the $`|\delta Q_{ki}|`$

$`|{\displaystyle \frac{\delta Q_{kD}}{\delta Q_{ki}}}|({\displaystyle \frac{\eta _t}{\eta _D}})^{1/2}=({\displaystyle \frac{T_D}{T_t}})^{1/2},`$ (36)
where $`T_t`$ is the temperature when $`Q`$ transits from the potential phase to the tracking regime. Making use of the initial value of $`\delta Q_k`$ generated during inflation, $`\delta Q_{ki}=\frac{H_{in}}{\sqrt{2k^3}}`$, we have

$`|\delta Q_{kD}|`$ $``$ $`{\displaystyle \frac{H_{in}}{\sqrt{2k^3}}}({\displaystyle \frac{T_D}{T_t}})^{1/2}.`$ (37)
Following , we calculate the RMS fluctuation of baryon number per logarithmic interval in our model which is given by

$$(\frac{\delta n_B}{n_B})_k\frac{H_{in}}{4\pi }|V^{^{}}(Q_D)/V(Q_D)|(\frac{T_D}{T_t})^{1/2}.$$
(38)
For the Quintessence model with potential in Eq.(9),

$`({\displaystyle \frac{\delta n_B}{n_B}})_k`$ $``$ $`{\displaystyle \frac{\lambda }{4\pi }}{\displaystyle \frac{H_{in}}{M_{pl}}}({\displaystyle \frac{T_D}{T_t}})^{1/2}`$ (39)
$``$ $`10^5({\displaystyle \frac{T_D}{T_t}})^{1/2},`$ (40)
which will be smaller than $`10^5`$ since $`T_t>T_D`$. To get the result above we have used the constraints on $`\frac{H_{in}}{M_{pl}}10^5`$ from the limit of the detecting of primordial gravitational waves predicted by inflation models. In addition we have also taken $`\lambda ^2𝒪(10)`$ to satisfy the constraints from the nucleosynthesis. The isocurvature baryon number fluctuation in (39) is small and is not conflict with the Boomerang and Maxima-1 data on the anisotropy of cosmic microwave background.

The isocurvature baryon number fluctuation in the spontaneous baryogenesis has been studied by Turner, Cohen and Kaplan , where they find the amplitude sizable. In their model the potential of the ilion field $`\theta `$ is $`m^2\theta ^2`$. The produced baryon asymmetry is linearly proportional to the initial value of the ilion field $`\theta _0`$. Therefore quantum fluctuations in the initial $`\theta _0`$ during inflation will directly induce isocurvature baryon number fluctuations. In our model the Quintessence has the tracking behavior, and the produced baryon number asymmetry is almost independent of the initial conditions of Quintessence field. Consequently a small isocurvature perturbation is expectable.

The tracking behavior of the Quintessence, however plays an important role in our discussions. Because of the tracking, $`\dot{Q}`$ can be large enough and the baryon number asymmetry required can be generated. For models without tracking property, such as the model with potential $`V(Q)=\mathrm{\Lambda }^4(1\pm \mathrm{cos}\frac{Q}{M})`$, we have checked and found that $`\dot{Q}`$ is too small to account for $`n_B/s10^{10}`$.

Our mechanism for leptogenesis, however does put a lower limit on the reheating temperature $`T_R`$. It requires $`T_R`$ be higher than $`T_t`$ and $`T_D`$. For cut-off $`M`$ of Eq.(1) in the range of $`10^{16}10^{19}`$GeV, $`T_D`$ is about $`10^710^{10}`$GeV, which will not put a strong constraint on models of inflation and reheating.

In summary we have introduced in this paper an effective interaction of the Quintessence scalar to the ordinary matter and studied the possibilities of baryogenesis. Our model can naturally explain the baryon number asymmetry of our Universe $`n_B/s10^{10}`$ via spontaneous leptogenesis. Our scenario provides a unified description for the present accelerated expansion and baryon number asymmetry generation of our Universe by the dynamics of Quintessence.

Acknowledgments: We thank David Kaplan for correspondence on spontaneous baryogenesis and Robert Brandenberger, Fabio Finelli, Z-H. Lin for discussions. This work is supported in part by National Natural Science Foundation of China under Grant No. 10047004 and by Ministry of Science and Technology of China under Grant No. NKBRSF G19990754.


format: text

text:
# A Balanced Truncation Primer

## I Introduction

Theoretical physics endeavors to produce models that predict observed physical phenomena. Though sometimes the challenge is to develop a mathematical language describing the system of interest, often (especially in the study of complex systems) one can write down the exact dynamics and gain little insight — the resulting expressions are too cumbersome, too messy, or too ill-conditioned to be useful. By exploiting symmetries and other characteristics of the particular system, one may find a simpler equivalent description of the dynamics. This description may be further simplified by using approximation techniques, e.g., asymptotic limits and small-parameter expansions. Much of the art of the field lies in finding and choosing ad hoc methods for deriving these simpler models; however, more systematic methods are clearly desirable.

Control theorists have developed a variety of model reduction techniques that systematically produce simple models of complex systems. These notes will describe balanced truncation Dullerud and Paganini (2000); Zhou et al. (1996), a model reduction technique for linear systems which is readily available in a variety of formats (e.g. MATLAB). Balanced truncation has recently been applied in physics contexts Rahn et al. (2001); Sznaier et al. (2001), and the results suggest it will prove a useful tool for treating large systems in both classical and quantum settings.

Much of this discussion comes directly from Dullerud and Paganini Dullerud and Paganini (2000). We will present the important concepts and results necessary to apply balanced truncation, omitting both the proofs and the algorithms. We refer the reader to Dullerud and Paganini (2000) and Zhou et al. (1996) for a more complete mathematical discussion, and to MATLAB toolboxes and their documentation alg  for the computational methods.

In section II we will describe the input-output paradigm of control theory, and introduce state-space models, the class of systems treatable by balanced truncation. Given an arbitrary state-space model, we will characterize the smallest state-space model with identical input-output characteristics in section III, and in section IV we will show how balanced truncation is used to find smaller models with controlled approximation errors.

## II Input-Output Maps and State-Space Models

In many physics settings one is more concerned with the macroscopic behavior of a large system, and less concerned with the system’s microscopic details. As an admittedly contrived example, consider a pendulum in a plane at whose free end is a tank partially filled with a classical fluid (see Fig. 1). Suppose that at time $`t=0`$ the system is at rest in its stable equilibrium, and our only method of disturbing the system is to exert a time-varying torque $`\tau (t)`$ at the pivot. Suppose further that we are concerned only with the time evolution of pendulum’s angle $`\theta (t)`$ — not with the distribution of the fluid, given by some high-dimensional variable $`\mathrm{\Phi }(t)`$.

An exact model including the full fluid state $`\mathrm{\Phi }`$ would give the system’s exact response to a driving torque, but would be quite impractical. As we only wish to describe the angular response to the driving torque (a mapping from one degree of freedom to another) one suspects a lower-dimensional model might suffice. For example, one might try treating the fluid as a point mass attached to the pendulum by a non-linear spring.

Systems of this sort are naturally phrased in a control theory language. In the typical control scenario, a time-varying input $`u(t)`$ drives a system with state $`x(t)`$ giving output signal $`y(t)`$, and the system dynamics have the form

$$\begin{array}{ccc}\hfill \dot{x}(t)& =& f(x(t),u(t))\hfill \\ \hfill y(t)& =& h(x(t),u(t)).\hfill \end{array}$$
(1)
Such systems are depicted by a block diagram as shown in Fig. 2. In the example of the pendulum, the system’s state is given by $`x=(\theta ,\dot{\theta },\mathrm{\Phi },\dot{\mathrm{\Phi }})`$ so as to describe the evolution with the first-order dynamics (1). The input to the system is $`u(t)=\tau (t)`$, and the system’s output is $`y(t)=\theta (t)`$.

Together with some initial condition (typically $`x(0)=0`$), the functions $`f`$ and $`h`$ in (1) define an input-output map $`\mathrm{\Psi }`$ taking $`u`$ to $`y`$. Since it is this relation with which we are concerned, rather than the system’s internal dynamics, a theoretical model for the system will suffice if it describes $`\mathrm{\Psi }`$. Given some system of the form (1), model reduction aims to produce simpler models (i.e. models with a lower-dimensional state $`x`$) that approximate the original input-output map.

We will consider models of the form

$$\begin{array}{ccc}\hfill \dot{x}& =& Ax+Bu\hfill \\ \hfill y& =& Cx+Du\hfill \end{array}$$
(2)
where $`x^n`$, $`u^m`$, and $`y^p`$, and $`A`$, $`B`$, $`C`$ and $`D`$ are time-independent real matrices of sizes $`n\times n`$, $`n\times m`$, $`p\times n`$ and $`p\times m`$ respectively. (This entire discussion also holds for complex-valued systems.) Such linear models are called state-space models, and the order of a model is $`n`$, the dimension of the state $`x`$. For compactness, the model with matrices $`A`$, $`B`$, $`C`$ and $`D`$ is denoted by

$$\left(\begin{array}{cc}A& B\\ & \\ C& D\end{array}\right).$$
(3)
(This notation should not be confused with the $`(n+p)\times (n+m)`$ matrix with blocks $`A`$, $`B`$, $`C`$ and $`D`$.) Often the output will only depend on the system state, i.e. $`D=0`$.

Consider the dynamics (2) under a change-of-basis $`z=Tx`$, where $`T`$ is invertible but need not be unitary. The dynamics may then be written as

$$\begin{array}{ccc}\hfill \dot{z}& =& TAT^1z+TBu\hfill \\ \hfill y& =& CT^1z+Du.\hfill \end{array}$$
(4)
Thus changing the basis of the state space defines a mapping on state-space models given by

$$\left(\begin{array}{cc}A& B\\ & \\ C& D\end{array}\right)\left(\begin{array}{cc}TAT^1& TB\\ & \\ CT^1& D\end{array}\right).$$
(5)
(Note that $`D`$ is unchanged.) We will call such maps similarity transformations. Because these transformations are merely a rewriting of the system dynamics, the input-output map remains the same.

As a given state-space model is not a unique description of an input-output map $`\mathrm{\Psi }`$, in the next section we ask: what is the lowest-order model with the same $`\mathrm{\Psi }`$ as the given model? After finding a lowest-order exact model, we will use balanced truncation to find lower-order models approximating $`\mathrm{\Psi }`$. This procedure is summarized in Fig. 3.

## III Lowest-Order Exact Models

To find state-space models of lowest order with the exact $`\mathrm{\Psi }`$ of a given model, we will split the problem into two parts: controllability and observability. We will then combine these ideas to find minimal realizations of the system.

### III.1 Controllability

Assume the system is initially in the state $`x(0)=0`$; given a time $`\tau >0`$, the controllable states $`x_f`$ are those for which there is an input signal $`u(t)`$ yielding $`x(\tau )=x_f`$. The dynamics (2) can be integrated to yield

$$x(\tau )=_0^\tau e^{A(\tau t)}Bu(t)𝑑t,$$
(6)
which gives a linear map $`\mathrm{\Psi }_c^{(\tau )}:ux(\tau )`$. The controllable states form the image of this map. Since the map is linear, the image is a subspace of the state-space $`^n`$, called the controllable subspace. We denote this subspace by $`𝒞_{AB}`$, as controllability only depends on the matrices $`A`$ and $`B`$.

Using properties of the matrix exponential it can be shown that the controllable subspace is the image of the controllability matrix:

$$𝒞_{AB}=\mathrm{Image}[BABA^2B\mathrm{}A^{n1}B].$$
(7)
Thus we see that $`C_{AB}`$ is independent of $`\tau `$. It can also be shown that the controllable subspace is the image of the controllability gramian, an $`n\times n`$ matrix given by

$$X_C=_0^\tau e^{At}BB^{}e^{A^{}t}𝑑t,$$
(8)
and the orthogonal subspace of uncontrollable states is given by the controllability gramian’s kernel. If $`C_{AB}=^n`$ (i.e. there exists an input signal $`u`$ to prepare any state $`x(\tau )`$) then we say that $`(A,B)`$ is a controllable pair.

It can be shown that given any $`A`$ and $`B`$, we can find a similarity transformation such that the transformed matrices have the block structure

$$\stackrel{~}{A}=TAT^1=\left[\begin{array}{cc}\stackrel{~}{A}_{11}& \stackrel{~}{A}_{12}\\ 0& \stackrel{~}{A}_{22}\end{array}\right]\stackrel{~}{B}=TB=\left[\begin{array}{c}\stackrel{~}{B}_1\\ 0\end{array}\right]$$
(9)
with $`(\stackrel{~}{A}_{11},\stackrel{~}{B}_1)`$ a controllable pair. Writing the state vector as $`x=(x_1,x_2)`$ corresponding to this block structure, we have

$$\begin{array}{ccccccc}\hfill \dot{x}_1& =& \stackrel{~}{A}_{11}x_1& +& \stackrel{~}{A}_{12}x_2& +& \stackrel{~}{B}_1u\hfill \\ \hfill \dot{x}_2& =& & & \stackrel{~}{A}_{22}x_2.& & \end{array}$$
(10)
Because $`x(0)=0`$, these dynamics yield $`x_2(t)=0`$ for all time. Thus the dynamics for $`x_1`$ reduce to

$$\begin{array}{ccc}\hfill \dot{x}_1& =& \stackrel{~}{A}_{11}x_1+\stackrel{~}{B}_1u.\hfill \end{array}$$
(11)
Because $`(\stackrel{~}{A}_{11},\stackrel{~}{B}_1)`$ is a controllable pair, we may choose an input $`u`$ to prepare any state $`x_1(\tau )`$, and thus the transformed controllable subspace $`𝒞_{\stackrel{~}{A}\stackrel{~}{B}}`$ is given by the states of the form $`(x_1,0)`$. The orthogonal subspace, given by states of the form $`(0,x_2)`$, is irrelevant to the input-output map since no input can affect these states.

### III.2 Observability

We now consider another problem with a similar structure. Suppose the system is in some initial state $`x(0)=x_0`$ and $`u=0`$ for all time. Based on the output $`y(t)`$ for $`0t\tau `$, can we uniquely identify $`x_0`$? Integrating the dynamics (2) yields

$$y(t)=Ce^{At}x_0,$$
(12)
which gives a linear map $`\mathrm{\Psi }_o^{(\tau )}:x_0y`$ (where by $`y`$ we mean the output signal $`y(t)`$ for $`0t\tau `$). Suppose two initial states $`x_0`$ and $`x_1`$ give the same $`y`$. As $`\mathrm{\Psi }_o^{(\tau )}`$ is linear, the initial state $`x_0x_1`$ must give $`y=0`$. We call initial states giving output $`y=0`$ unobservable, since any unobservable state may be added to any other initial state without changing the output.

The unobservable states form the kernel of $`\mathrm{\Psi }_o^{(\tau )}`$; as the map is linear these states form a subspace, called the unobservable subspace and denoted by $`𝒩_{CA}`$. It can be shown that the unobservable subspace is given by the kernel of the observability matrix:

$$𝒩_{CA}=\mathrm{ker}\left[\begin{array}{c}C\\ CA\\ \mathrm{}\\ CA^{n1}\end{array}\right].$$
(13)
Thus $`𝒩_{CA}`$ is independent of $`\tau `$. It can also be shown that $`𝒩_{CA}`$ is given by the kernel of the observability gramian

$$Y_O=_0^\tau e^{A^{}t}C^{}Ce^{At}𝑑t,$$
(14)
and the orthogonal subspace of observable states is given by the observability gramian’s image. If $`𝒩_{CA}=\{\stackrel{}{0}\}`$, then the entire space is observable, and we say that $`(C,A)`$ is an observable pair. As the observable states are given by the image of (14) and the controllable states are given by the image of (8), $`(C,A)`$ is an observable pair if and only if $`(A^{},C^{})`$ is a controllable pair.

It can be shown that given any $`C`$ and $`A`$, we can find a similarity transformation such that the transformed matrices have the block structure

$$\begin{array}{ccccc}\hfill \stackrel{~}{A}& =& TAT^1& =& \left[\begin{array}{cc}\stackrel{~}{A}_{11}& 0\\ \stackrel{~}{A}_{21}& \stackrel{~}{A}_{22}\end{array}\right]\hfill \\ & & & & \\ \hfill \stackrel{~}{C}& =& CT^1& =& \left[\begin{array}{cc}\stackrel{~}{C}_1& \mathrm{\hspace{0.33em}0}\end{array}\right]\hfill \end{array}$$
(15)
with $`(\stackrel{~}{C}_1,\stackrel{~}{A}_{11})`$ an observable pair. Writing the state vector as $`x=(x_1,x_2)`$ corresponding to this block structure, we have

$$\begin{array}{cccccc}\hfill \dot{x}_1& =& \stackrel{~}{A}_{11}x_1& & & \\ \hfill \dot{x}_2& =& \stackrel{~}{A}_{21}x_1& +& \stackrel{~}{A}_{22}x_2& \\ \hfill y& =& \stackrel{~}{C}_1x_1.& & & \end{array}$$
(16)
Thus the time evolution of $`x_1`$ is never affected by $`x_2`$, and the output signal $`y`$ only depends on $`x_1`$. Because $`(\stackrel{~}{C}_1,\stackrel{~}{A}_{11})`$ is an observable pair, we can uniquely identify an initial state $`x(0)=(x_1,0)`$ based on $`y`$. The transformed unobservable subspace $`𝒩_{\stackrel{~}{C}\stackrel{~}{A}}`$ is given by the states of the form $`(0,x_2)`$, and is irrelevant to the input-output map since no output can be affected by these states.

### III.3 Minimal Realizations

The notions of controllability and observability give us a means of deciding whether a state affects the system’s input-output map: if a state is unobservable, it does not affect the output, and if a state is uncontrollable, it is unaffected by the input. Only those states that are both controllable and observable are of relevance.

We say that a state-space model given by matrices $`(A,B,C,D)`$ is a minimal realization if no lower-order model gives the same input-output map $`\mathrm{\Psi }`$. The intuition above can be made precise as follows: it can be shown that a model is a minimal realization if and only if all states are both controllable and observable, i.e. $`(A,B)`$ is a controllable pair and $`(C,A)`$ is an observable pair.

Given a state-space model given by $`A`$, $`B`$, $`C`$ and $`D`$ we may find a minimal realization by isolating only those dimensions which are both controllable and observable. To do so, we perform a Kalman decomposition, which simultaneously performs the transformations (9) and (15) as follows. For any state-space model there exists a similarity transformation such that the matrices of transformed model have the block structure

$$\left(\begin{array}{cc}TAT^1& TB\\ & \\ CT^1& D\end{array}\right)=\left(\begin{array}{ccccc}\stackrel{~}{A}_{11}& 0& \stackrel{~}{A}_{13}& 0& \stackrel{~}{B}_1\\ \stackrel{~}{A}_{21}& \stackrel{~}{A}_{22}& \stackrel{~}{A}_{23}& \stackrel{~}{A}_{24}& \stackrel{~}{B}_2\\ 0& 0& \stackrel{~}{A}_{33}& 0& 0\\ 0& 0& \stackrel{~}{A}_{43}& \stackrel{~}{A}_{44}& 0\\ & & & & \\ \stackrel{~}{C}_1& 0& \stackrel{~}{C}_3& 0& D\end{array}\right)$$
(17)
and, writing the state vector as $`x=(x_1,x_2,x_3,x_4)`$ corresponding to the block structure, the controllable states are of the form $`(x_1,x_2,0,0)`$ and the observable states are of the form $`(x_1,0,x_3,0)`$. Thus only states of the form $`(x_1,0,0,0)`$ are both controllable and observable.

Eliminating all but the states $`(x_1,0,0,0)`$ yields the state-space model

$$\left(\begin{array}{cc}\stackrel{~}{A}_{11}& \stackrel{~}{B}_1\\ & \\ \stackrel{~}{C}_1& D\end{array}\right).$$
(18)
Since we have only eliminated states irrelevant to the input-output map, this reduced system has the exact same $`\mathrm{\Psi }`$ as the original system. Further, since states of the form $`(x_1,0,0,0)`$ were both controllable and observable, $`(\stackrel{~}{A}_{11},\stackrel{~}{B}_1)`$ is a controllable pair and $`(\stackrel{~}{C}_1,\stackrel{~}{A}_{11})`$ is an observable pair; thus this model is a minimal realization.

Sometimes the uncontrollable and unobservable states are obvious from the form of a physical system’s dynamics, e.g. some degrees of freedom are uncoupled, or some symmetry can be exploited. In other circumstances, the physics may not make it clear which states can be eliminated without affecting the input-output map. Especially in the latter situation, an algorithmic method such as the Kalman decomposition (available in MATLAB) can be quite advantageous. Nonetheless, one feels intuitively that one should always be able to find a minimal realization analytically if one is sufficiently clever. However, we do not expect in general that standard analytic methods will be useful in seeking lower-order approximations, as we will do in the next section.

## IV Lower-Order Approximate Models

To apply balanced truncation, we will first assume that we have reduced the model to a minimal realization. We make the further assumption that the resulting system is exponentially stable, i.e. all eigenvalues of $`A`$ have strictly negative real part. (Various methods exist to extend these methods to unstable systems, e.g., Sznaier et al. (2001), but we assume stability to prove the standard result.)

### IV.1 Quantifying Observability and Controllability

In the previous section we distinguished between states that were observable and unobservable; we now wish to quantify the observability of the observable states. To do so, consider the output signal for $`t0`$ that results from the initial state $`x(0)=x_0`$ when there is no input $`u`$. This signal is given by (12). Using the $`_2`$ norm for the signal $`y`$, we have

$$\begin{array}{ccc}\hfill y^2& =& _0^{\mathrm{}}y^{}(t)y(t)𝑑t\hfill \\ & =& x_0^{}\left(_0^{\mathrm{}}e^{A^{}t}C^{}Ce^{At}𝑑t\right)x_0\hfill \\ & =& x_0Y_Ox_0\hfill \end{array}$$
(19)
where $`Y_O`$ is the observability gramian given by (14) with $`\tau \mathrm{}`$. Recall that the kernel of $`Y_O`$ is independent of $`\tau `$, so choosing $`\tau \mathrm{}`$ will not change which states are observable and which states are unobservable. (System stability is required to ensure convergence of the integral in this limit.) Scaling the norm-squared of the output by the norm-squared of the initial state yields

$$\frac{y^2}{x_0^2}=\frac{x_0^{}Y_Ox_0}{x_0^{}x_0},$$
(20)
which quantifies the observability of the states in the direction of $`x_0`$.

From (19) we see that the observability gramian is Hermitian and positive semidefinite. As we have a minimal realization, all states $`x_0`$ are observable ($`x_0Y_Ox_0>0`$), and so $`Y_O`$ is strictly positive definite. Thus in geometric terms analogous to the moment-of-inertia tensor, $`Y_O`$ defines an “observability ellipsoid” in the state space, with the longest principal axes along the most observable directions (see Fig. 4). A similarity transformation given by $`T`$ transforms the observability gramian by $`Y_O(T^1)^{}Y_OT^1`$. As $`T`$ need not be unitary, this transformation may rescale the ellipsoid’s axes as well as rotate them.

We quantify controllability in a similar fashion. Suppose that at a time well in the past ($`\tau \mathrm{}`$) the system is in the state $`x=0`$, and some input $`u(t)`$ drives the system for $`t0`$, yielding a final state $`x(0)=x_0`$. As we have a minimal realization, all states $`x_0`$ are controllable and therefore can be prepared in this fashion. For each state $`x_0`$ there is a minimum signal size $`u_{\mathrm{opt}}^2=_{\mathrm{}}^0u(t)^{}u(t)𝑑t`$ required to yield $`x(0)=x_0`$. The smaller this minimum signal, the more sensitive this state is to the input signal. Thus states with a smaller $`u_{\mathrm{opt}}^2`$ are said to be more controllable.

Consider the controllability gramian $`X_C`$ given by (8) with $`\tau \mathrm{}`$. Just as $`Y_O`$, $`X_C`$ is Hermitian and positive semidefinite, and because we have assumed a minimal realization, $`X_C`$ is strictly positive definite and therefore invertible. It can be shown that

$$u_{\mathrm{opt}}^2=x_0^{}X_C^1x_0.$$
(21)
It follows that

$$\left(\frac{u_{\mathrm{opt}}^2}{x_0^2}\right)^1=\frac{x_0^{}X_Cx_0}{x_0^{}x_0},$$
(22)
which quantifies the controllability of the states in the direction of $`x_0`$. Just as with $`Y_O`$, $`X_C`$ defines a “controllability ellipsoid” in state space, with the longest principal axes along the most controllable directions. A similarity transformation given by $`T`$ transforms the controllability gramian by $`X_CTX_CT^{}`$.

We have thus found that the gramians give us a useful measure of a state’s observability and controllability. One might ask whether the observability and controllability matrices of (13) and (7) could serve a similar purpose, but in fact they are only useful for determining whether a given state is observable/controllable. One might also ask if $`\tau \mathrm{}`$ is necessary — we could choose finite $`\tau `$, and quantify controllability and observability on a finite time horizon. However, finite $`\tau `$ does not lead to the error bound in the approximations of the next section, which is the main result of these notes.

### IV.2 Balanced Truncation

With the above quantification of observability and controllability, one might be tempted to prescribe some algorithm like eliminating the least observable or least controllable dimensions in the state space to yield a lower-order approximate model. However, such an approach would not necessarily be successful. Suppose, for example, that the least observable states were in the direction of the unit vector $`\widehat{x}`$, but that states in this direction were extremely controllable. Thus a small signal $`u`$ might lead to the internal state $`x=\lambda \widehat{x}`$ with $`\lambda `$ large. Though this state is the least observable, $`\lambda `$ might be sufficiently large that the resulting output signal is non-negligible.

Instead, we wish to use observability and controllability to yield a single measure of a state’s importance to the input-output map $`\mathrm{\Psi }`$. It can be shown that given two positive definite square matrices $`X_C`$ and $`Y_O`$ of the same size, there exists an invertible $`T`$ such that

$$TX_CT^{}=(T^1)^{}Y_OT^1=\mathrm{\Sigma }$$
(23)
where $`\mathrm{\Sigma }`$ is diagonal with positive real diagonal entries. Such a similarity transformation is called a balancing transformation. Geometrically, balancing transforms the observability and controllability ellipsoids so that they are identical and their principal axes lie on the coordinate axes of the state space (see Fig. 5).

The resulting $`\mathrm{\Sigma }`$ is unique up to permutation of the diagonal elements, so we may choose $`T`$ yielding the transformed gramians

$$\stackrel{~}{X}_C=\stackrel{~}{Y}_O=\mathrm{\Sigma }=\left[\begin{array}{cccc}h_1& & & \\ & h_2& & \\ & & \mathrm{}& \\ & & & h_n\end{array}\right]$$
(24)
with $`h_1h_2\mathrm{}h_n>0`$. The $`h_i`$ are called Hankel Singular Values (HSVs); in this transformed system $`h_i`$ is the quantitative measure of both the observability and controllability of the unit basis vector $`\widehat{e}_i`$. Thus the basis vectors have been sorted in order of relevance to the input-output map.

Once the system is balanced, we may truncate the state-space dimensions with low HSVs to yield lower-order approximate models. Intuitively, the smaller the HSVs corresponding to truncated dimensions, the better the approximation. We will now make this idea precise.

Let the balanced state-space model with sorted HSVs be given by matrices $`\stackrel{~}{A}`$, $`\stackrel{~}{B}`$, $`\stackrel{~}{C}`$ and $`D`$. Choosing some $`r`$ such that $`h_r`$ is strictly greater than $`h_{r+1}`$, we write the state vector $`x=(x_1,x_2)`$ where $`x_1`$ gives the first $`r`$ coordinates, and $`x_2`$ gives the last $`nr`$ coordinates. We then write the state-space model in the corresponding block structure

$$\left(\begin{array}{cc}\stackrel{~}{A}& \stackrel{~}{B}\\ & \\ \stackrel{~}{C}& D\end{array}\right)=\left(\begin{array}{ccc}\stackrel{~}{A}_{11}& \stackrel{~}{A}_{21}& \stackrel{~}{B}_1\\ \stackrel{~}{A}_{12}& \stackrel{~}{A}_{22}& \stackrel{~}{B}_2\\ & & \\ \stackrel{~}{C}_1& \stackrel{~}{C}_2& D\end{array}\right)$$
(25)
and truncate the $`nr`$ least significant dimensions of the model, yielding the order $`r`$ model

$$\left(\begin{array}{cc}\stackrel{~}{A}_{11}& \stackrel{~}{B}_1\\ & \\ \stackrel{~}{C}_2& D\end{array}\right).$$
(26)
This procedure is called balanced truncation.

We now give bounds on the resulting approximation error. Let $`u(t)`$ be some input signal on $`t(\mathrm{},\mathrm{})`$ that is finite with respect to the $`_2`$ norm

$$z^2=_{\mathrm{}}^{\mathrm{}}z(t)^{}z(t)𝑑t.$$
(27)
Let $`y(t)`$ and $`\stackrel{~}{y}(t)`$ be the resulting output signals of the original and truncated systems respectively. Stability of the original system (25) and the strict inequality $`h_r>h_{r+1}`$ guarantees stability of the truncated system (26); since both systems are stable, the output signals $`y`$ and $`\stackrel{~}{y}`$ are also finite with respect to the norm (27), as is the error signal $`\stackrel{~}{y}y`$. Now, letting $`h_1^{\mathrm{tr}}>h_2^{\mathrm{tr}}>\mathrm{}>h_k^{\mathrm{tr}}`$ be the distinct HSVs of the truncated $`nr`$ dimensions, it can be shown that

$$h_1^{\mathrm{tr}}\underset{u}{\mathrm{max}}\frac{\stackrel{~}{y}y}{u}2\underset{j=1}{\overset{k}{}}h_j^{\mathrm{tr}}$$
(28)
Thus we have both an upper and lower bound on the worst-case errors resulting from balanced truncation. Given the Hankel Singular Values of a system, we may truncate to a system of desired order and bound the resulting error, or we may choose the order of truncation based on the maximum tolerable error.

The advantages of balanced truncation as method of approximating the input-output map are two-fold. First, the error is bounded as above, in contrast to more traditional approximations where errors estimated by the order of some small parameter, e.g. $`O(ϵ^2)`$. Second, the algorithmic process is efficient. Balanced truncation does not necessarily yield the optimal order-$`r`$ approximation in the sense of the error bound (28), but finding the optimal approximation may be computationally difficult.

For single-input single-output systems (i.e. $`u`$, $`y`$), the error bound may also be described in terms of frequency responses. Given a sinusoidal input, a state-space model will give an output with the same frequency and a frequency-dependent amplitude and phase shift:

$$u(t)=\mathrm{sin}(\omega t)y(t)=A_\omega \mathrm{sin}(\omega t+\varphi _\omega ).$$
(29)
Let $`A_\omega `$ and $`\varphi _\omega `$ be the amplitude and phase shifts for the exact system, and let $`\stackrel{~}{A}_\omega `$ and $`\stackrel{~}{\varphi }_\omega `$ be the shifts for the approximate system. Writing the exact and approximate shifts in complex form, we have

$$h_1^{\mathrm{tr}}\underset{\omega }{\mathrm{max}}\left|\stackrel{~}{A}_\omega e^{i\stackrel{~}{\varphi }_\omega }A_\omega e^{i\varphi _\omega }\right|2\underset{j=1}{\overset{k}{}}h_j^{\mathrm{tr}}.$$
(30)
Thus even though sinusoidal inputs and outputs are unbounded in the $`_2`$ norm (27), the error is controlled in this fashion.

## V Conclusion

According to Dullerud and Paganini (2000), balanced realizations first appeared in the control literature in 1981, and the proof of the error bound on truncated models first appeared in 1984. Until recently, however, the physics community has made little use of this powerful tool. We believe that balanced truncation will be of particular value when building simulations of and theoretical models for the evolution of macroscopic quantities in large complex systems, and it is hoped that these notes will be helpful in such studies.


text:
# 1 Introduction

## 1 Introduction

Quarks are among the most fundamental building blocks of nature. Unlike the previous elementary constituents, quarks are found to possess a remarkable property, confinement, i.e., we cannot isolate quarks by dissociating hadrons. Historically, fundamental theories of nature have emerged through reactions of isolated elements to weak external perturbations. Quark confinement, however, has the consequence that we cannot directly measure the fundamental properties of quarks by experiment.

Fortunately, quarks are found to possess another remarkable property, asymptotic freedom, whereby interactions among quarks become weak at high energies. This enabled us to identify quantum chromodynamics (QCD) as the fundamental theory of quarks. QCD has been quite successful in explaining high energy processes of particles.

On the other hand, because the coupling parameter of QCD becomes large at low energies, basic properties of hadrons, such as mass spectrum, spatial sizes, decay constants, etc., are not calculable by conventional analytic methods. In other words, we cannot reconstruct hadrons from quarks based on the dynamics of QCD yet. In order to answer the fundamental question of whether QCD is correct also at low energies, we have to carry out non-perturbative calculations of hadrons directly from the first principles of QCD.

Calculation of low-energy hadron properties is also required for determination of fundamental parameters of nature. Due to quark confinement, quark masses and coupling parameters have to be inferred indirectly from a comparison of experimental results for hadron masses etc. with a theoretical calculation of them as functions of the fundamental parameters of QCD. Also for electroweak interactions, although the core part of the reaction can be reliably calculated in perturbation theory, uncertainties in QCD corrections to the quark currents propagate to the determination of the Cabibbo-Kobayashi-Maskawa parameters.

Numerical simulation of quarks based on the lattice formulation of QCD is currently the only method to calculate low-energy non-perturbative properties of hadrons. Recently, predictions with big impact are beginning to be produced through development of dedicated parallel computers . In particular, a big breakthrough was achieved by the CP-PACS computer developed at the University of Tsukuba . In this paper, I summarize major results from CP-PACS, focusing on the most fundamental topics of QCD, i.e., the light hadron spectrum and light quark masses.

In Sec. 2, I introduce the lattice formulation of QCD and explain why dedicated computers have been developed by physicists. Results for the light hadron spectrum and light quark mass are presented in Secs. 3 and 4. Conclusions are given in Sec. 5.

## 2 Lattice QCD and dedicated parallel computers

We formulate QCD on a four-dimensional hypercubic lattice with a finite lattice spacing $`a`$ and finite lattice size $`Na`$. The real world is defined by the limit of vanishing lattice spacing $`a0`$ keeping the lattice size $`Na`$ finite or sufficiently large (the continuum limit). Before taking this limit, the theory is finite and mathematically well-defined. Therefore, we can apply various non-perturbative techniques. Here, we perform numerical simulations on finite lattices to calculate hadrons.

To accomplish this, the most time-consuming part of the calculation is the inversion of quark propagation kernels. Numerically, this is an inversion of large sparse complex matrices. For the case of Wilson-type lattice quark actions, the typical size of the matrix is $`12V\times 12V`$, where $`12=3\times 4`$ is the freedom of color and spin, and $`V=N^4`$ is the lattice volume. Because the condition number is inversely proportional to the quark mass, the inversion is numerically more intensive when we decrease the quark mass. Even with the latest supercomputers, it is difficult to simulate the light $`u,d`$ quarks directly. Therefore, in addition to the continuum extrapolation, we need to extrapolate the results at typically around the $`s`$ quark mass to the physical $`u,d`$ quark mass point. This procedure is called “the chiral extrapolation”.

In summary, in order to extract a prediction for the real world, we have to perform continuum and chiral extrapolations. To get a precise and reliable result, it is essential to have good control of these extrapolations. This requires a large-scale systematic calculation, and thus huge computer power. The power of vector supercomputers in 1980’s was not sufficient, which motivated several groups of lattice physicists to construct parallel computers dedicated to lattice calculations .

The historical development of computer speed is shown in Fig. 2. The frontier of speed has been advanced by parallel computers since around 1990. Dedicated machines developed by physicists contributed much to this trend. At the University of Tsukuba, we have developed two machines, QCDPAX (1990)  and CP-PACS (1996) . CP-PACS is a parallel computer achieving peak performance of 614.4 GFLOPS with 2048 single-processor nodes. Since the first power-on in 1996, intensive lattice QCD calculations have been made on CP-PACS. In the following, I focus on the studies of hadron spectrum and light quark masses.

## 3 Light hadron mass spectrum

Precise calculation of the hadron mass spectrum directly from the first principles of QCD is one of the main goals of lattice QCD. Because the computer power required is enormous, we study the issue step by step as follows:

: Calculate in the quenched approximation, in which the effects of dynamical pair creation and annihilation of quarks are neglected.

: Include the dynamical $`u,d`$ quarks in a degenerate approximation, while heavier quarks are treated in the quenched approximation (two-flavor full QCD).

: Include the dynamical $`s`$ quark (2+1 flavor full QCD).

: Introduce the $`u,d`$ mass difference, dynamical $`c`$ quark, etc.

### 3.1 Quenched studies

With the quenched approximation, we can reduce the computer time by a factor of several hundred preserving the basic properties of QCD (confinement, asymptotic freedom, and spontaneous breakdown of chiral symmetry). The effects of the approximation are expected to be about 10% in the spectrum.

The first studies in this approximation had already been made in early 80’s. The issue turned out to be quite tough and computationally demanding. Actually, it took about ten years until the first systematic study, performing all the extrapolations, was attempted in 1993 . From this study, the quenched light hadron spectrum was found to be consistent with experiment within the errors of about 10%. However, the quality of data was insufficient to test the accuracy of numerical extrapolations, and also the final errors were too large to resolve the quenching artifact.

The status was significantly improved by CP-PACS . From an intensive computation involving about 100 times more floating-point calculations than the previous studies, it became possible to perform systematic tests on the quality of the extrapolations. Errors in the final results for the light hadron spectrum are now confidently estimated to be less than about 1% for mesons and about 3% for baryons. These errors include statistical and all systematic errors except for those from the quenched approximation itself.

The quenched light hadron spectrum is shown in Fig. 2. Experimental values of meson masses — $`M_\pi `$, $`M_\rho `$, and either $`M_K`$ or $`M_\varphi `$ — are used as inputs to fix the lattice spacing $`a`$ (or the strong coupling constant), the average $`u,d`$ quark mass $`m_{ud}`$, and the $`s`$ quark mass $`m_s`$. Masses for other hadrons are among the predictions of QCD. From Fig. 2, we find that the global pattern of the spectrum is well reproduced. At the same time, we clearly see that the quenched spectrum deviates from experiment by about 10% (7 standard deviations) for the worst case: Results from the $`M_K`$-input and $`M_\varphi `$-input are discrepant by, and the hyperfine splitting between $`K`$ and $`K^{}`$ mesons is smaller than experiment by about 10%. Decuplet baryon splittings are also small. Because all other systematic errors are well controlled, we identify the discrepancy as quenching artifact.

### 3.2 Two-flavor full QCD

From the quenched simulation we find that, to calculate hadronic quantities with precision better than 10%, we have to incorporate dynamical quarks. A naive extension of the quenched simulation to full QCD is difficult because several hundred times more computer time is required. A partial solution is given by improvement of the lattice theory, whereby continuum properties are realized on coarser lattices. From a preparatory study , we find that the combination of renormalization-group improved glue action and clover-improved Wilson quark action is effective in removing major lattice artifacts. This reduces the computer time by a factor of about ten.

Although more might still be hoped for, we can start the first systematic studies of two-flavor full QCD on the CP-PACS performing both chiral and continuum extrapolations . To identify dynamical quark effects clearly, we carried out another quenched simulation using the same improved action.

Figure 4 shows the continuum extrapolation of vector meson masses. We find that the two quenched results (open and filled squares) lead to universal values in the continuum limit $`a=0`$. They deviate from the experimental values (diamonds) as noted in the previous subsection. On the other hand, the full QCD results (circles) extrapolate to values much closer to experiment. Accordingly, we find no big discrepancies between the results from $`M_K`$ and $`M_\varphi `$ inputs. This means that the quenching artifacts in the spectrum are mostly removed by introducing dynamical $`u,d`$ quarks.

The whole light hadron spectrum can now be approximately reproduced from QCD by adjusting just three parameters; the strong coupling constant, $`m_{ud}`$ and $`m_s`$. This provides us with a strong confirmation that QCD is correct also at low energies. Remaining small deviations from experiment may be explained by the quenched approximation of the $`s`$ quark. In order to confirm this, however, uncertainties from the chiral and continuum extrapolations should be reduced to the level of our quenched study. This goal is reserved for future work.

## 4 Light quark masses

Adjustment of $`m_{ud}`$ and $`m_s`$ in the calculation of the hadron spectrum provides us with the most direct determination of quark masses from QCD. Fig. 4 shows $`m_s`$ in the $`\overline{\mathrm{MS}}`$ scheme at $`\mu =2`$ GeV as functions of the lattice spacing $`a`$. Results from the quenched improved action, which are consistent with the quenched standard action in the continuum limit, are omitted for clarity.

On the lattice, there are several alternative definitions for the quark mass, with the difference being $`O\left(a\right)`$. In Fig. 4, they are denoted as VWI and AWI (vector and axial-vector Ward-Takahashi identity). See  for details. As shown in Fig. 4, different definitions lead to different values of $`m_q`$ at finite lattice spacings. This has been a big source of uncertainty in the previous calculations. From Fig. 4, we see that they converge to a universal value in the continuum limit, as theoretically expected. This confirms the quality of our calculations.

The quenched value for $`m_s`$, however, differs by about 20% between $`M_K`$-input and $`M_\varphi `$-input. This is equivalent to the quenched artifact discussed in Sec. 3.1. We note that the quenched artifact is larger in the light quark masses than in the hadron spectrum. When we turn on the dynamical $`u,d`$ quarks, the discrepancy between the inputs disappears within our errors, in accord with the fact that the full QCD reproduces the hadron spectrum better (Sec. 3.2).

The light quark masses in the continuum limit are summarized in Table 1. The errors include our estimates for systematic errors from chiral and continuum extrapolations and renormalization factors. We note that the masses from two-flavor QCD are 20–30% smaller than those from the quenched QCD. In particular, our $`s`$ quark mass in $`N_f=2`$ QCD is about 90 MeV, which is significantly smaller than the value $`150`$ MeV often used in phenomenology. Our results are, however, consistent with recent estimates $`m_s=83`$–130 MeV and $`m_{ud}=3.4`$–5.3 MeV from QCD sum rules  and $`m_s/m_{ud}=24.4\pm 1.5`$ from one-loop chiral perturbation theory .

## 5 Conclusions

An intensive calculation of the light hadron spectrum in the quenched approximation has revealed discrepancies of about 10% from experiment. With two flavors of dynamical $`u,d`$ quarks, the first systematic study performing both the continuum and chiral extrapolations has shown that the discrepancies are mostly removed, providing us with a strong confirmation that QCD is the correct theory of quarks at low energies as well.

Precision calculation of hadron masses enables us to determine fundamental parameters of quarks directly from the first principles of QCD. We find that the dynamical quark effect is as large as 20–30% in light quark masses. Noticeable and sizable dynamical quark effects are observed also in B meson decay constants , the equation of state at high temperatures , and the topological structure of the QCD vacuum . At the same time, new types of lattice fermions have begun to be applied to study hadronic matrix elements relevant to the $`\mathrm{\Delta }I=1/2`$ rule and the CP violation parameters .

Because a shift in fundamental parameters can have significant implications to phenomenological studies of the standard model, it is urgent to evaluate dynamical quark effects in other hadronic quantities as well. The influence of dynamical $`s`$ quark should also be studied. Studies in these directions are under way to clarify the precise structure of the standard model.

The studies presented in this report have been performed by the CP-PACS Collaboration. I thank the members of the Collaboration for discussions and comments. This work is in part supported by the Grants-in-Aid of Ministry of Education, Science and Culture (Nos. 12304011 and 13640260) and JSPS Research for Future Program.


text:
# 1 The integrated densities (see text for definition) for various ground state or isomeric nuclear configurations at constant density and 𝑍/𝐴 ratio. The size of the box was equal to 𝑑=26𝑓𝑚 - subfigures a) and b), 𝑑=23.4𝑓𝑚 - subfigures c) and d), 𝑑=20.8𝑓𝑚 -subfigures e) and f). The spherical nuclear phase is shown in the subfigures a) (scc) and d) (bcc), rod–like phase – subfigures b) and c), slab–like phase – subfigure e), bubble–like phase: subfigure f). Note different lattice types for the same nuclear shape.

## Abstract

In the present study we investigate the static properties of nuclei in the inner crust of neutron stars. Using the Hartree-Fock method in coordinate space, together with the semiclassical approximation, we examine the patterns of phase transitions.

The properties of the outer parts of neutron stars are important for the understanding of several observational issues. Despite the fact that the crust is typically less than 10% of the radius of the star one expects that such phenomena as: thermal X-ray emission from the stellar surface, X-ray burst sources, or the sudden speed-ups in the rotation rate of some neutron stars may be related to the structure of the crust.

In the inner crust of the star, the increasing pressure and density force nuclei to loose some of their neutrons. Hence deeper in the star, the nuclei are immersed in a neutron gas. Since nuclei in such an environment cannot be explored experimentally, our understanding of matter under these conditions is based on modelling. An agreement has been reached in the literature concerning the existence of the following chain of phase transitions as the density increases: nuclei $``$ rods $``$ plates $``$ tubes $``$ bubbles $``$ uniform nuclear matter . The appearance of different phases is attributed to the interplay between Coulomb and surface energies. Most published works are based on the minimization of some density functional in a single Wigner–Seitz cell and neglect the shell effects. The Hartree–Fock calculations performed so far are limited to Wigner–Seitz cells with geometries reflecting the nuclear shape and not the lattice structure . Hence the shell effects associated with the scattering of unbound neutrons on nuclear inhomogeneities were not properly taken into account.

There is no well established terminology for the energy corrections we are considering here, even though the problem has been addressed to some extent by other authors. In the case of finite systems, the energy difference between the true binding energy and the liquid drop energy of a given system is typically referred to as shell correction energy. In field theory a somewhat similar energy appears, due to various fluctuation induced effects and it is generically referred to as the Casimir energy :

$$E_{Casimir}=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}𝑑\epsilon \epsilon [g(\epsilon ,𝐥)g_0(\epsilon )],$$
(1)
where $`g_0(\epsilon )`$ is the density of states per unit volume for the fields in the absence of any objects, $`g(\epsilon ,𝐥)`$ is the density of states per unit volume in the presence of some “foreign” objects, such as plates, spheres, etc., and $`𝐥`$ is an ensemble of geometrical parameters describing these objects and their relative geometrical arrangement. A similar formula can be written for neutron matter energy

$$E_{nm}=_{\mathrm{}}^\mu 𝑑\epsilon \epsilon g(\epsilon ,𝐥)_{\mathrm{}}^{\mu _0}𝑑\epsilon \epsilon g_0(\epsilon ,𝐥),$$
(2)
with the notable difference in the upper integration limit. In the above equation $`g_0(\epsilon ,𝐥)`$ stands for the Thomas–Fermi or liquid drop density of states of the inhomogeneous phase and $`g(\epsilon ,𝐥)`$ is the true quantum density of states in the presence of inhomogeneities. The parameters: $`\mu `$ and $`\mu _0`$ are determined from the requirement that the system has a given average density:

$$\rho =_{\mathrm{}}^\mu 𝑑\epsilon g(\epsilon ,𝐥)=_{\mathrm{}}^{\mu _0}𝑑\epsilon g_0(\epsilon ,𝐥).$$
(3)
Since in infinite matter the presence of various inhomogeneities does not lead to the formation of discrete levels, the effects we shall consider here arise from the scattering states, which is in complete analogy with the procedure for computing the Casimir energy.

We shall only consider here spherical bubble–like, rod–like and plate–like phases. One can distinguish two types of “bubbles”: i) nuclei–like structures embedded in a neutron gas and ii) void–like structures. By voids we mean the regions in which the nuclear density is significantly lower than in the surrounding space. In the first case i), the single particle wave functions can be roughly separated into two classes, those localized mostly inside the nuclei–like structures and those which are completely delocalized. A fermion in a delocalized state will spend some time inside the “nuclei” too, but since the potential experienced by a nucleon is deeper there, the local momentum is larger and thus the relative time and relative probability to find a nucleon in this region is smaller . One can approximately replace then the “nuclei” by an effective repulsive potential of roughly the same shape. In the case of a “bubble”, when the probability to find a nucleon inside a “bubble” is reduced, again such an approximation appears reasonable. The same holds true for rod–like or plate–like structures. There are of course a number of “resonant” delocalized states, whose amplitudes behave in an opposite manner. However, the number of such “resonant” states is small and brings only small corrections. In all these phases the shell effects depend on the structure and stability of periodic orbits in the system .

In the semiclassical approximation it can be shown that the leading contribution to the shell energy (which may be thought of as the interaction energy), for two obstacles being either spherical ($`i=0`$), rod–like ($`i=1`$) or slab–like ($`i=2`$) nuclei located at a distance $`d`$ is given by :

$$E_{shell}\frac{\mathrm{}^2L^iR^{2i}}{8m_n}\left(\frac{3}{\pi }\right)^{\frac{2+i}{6}}\frac{(\rho _n^{out})^{\frac{2+i}{6}}}{d^{\frac{6i}{2}}}\mathrm{cos}\left(2k_F^ndi\frac{\pi }{4}\right),$$
(4)
where $`m_n`$ is the neutron mass, $`k_F^n`$ and $`\rho _n^{out}`$ denote the Fermi momentum and density of unbound neutrons, respectively111 In the Refs. the factor “$`1/\pi `$” has been missed in the expression for the shell energy of rod–like phase. However it does not change the conclusions.. We have assumed that rods and slabs are parallel to each other. In this equation, $`L`$ defines the length of the obstacle, and $`R`$ its radius (in the case of a slab it is defined as half of its width). It turns out that the semiclassical estimations yield energy corrections of the same order of magnitude as the liquid drop energy differences between the various phases of the inhomogeneous nuclear matter . It has lead to the hypothesis that the inner crust may have a quite complicated structure, probably even completely disordered .

Unfortunately the semiclassical approach has some drawbacks. First, it assumes that the obstacles are inpenetrable scatterers which may overestimate the amplitude of the shell energy. Second, the method lacks of the mutual interplay between the shell energy and the liquid drop energy. Such a coupling is quite important since it allows a part of the shell energy to be “absorbed” into deformation. Hence a microscopic treatment of the problem is needed where both effects, i.e. the one coming from the liquid drop energy (where mainly the surface and Coulomb terms are involved) and the shell energy term are treated on the same footing. Therefore we have applied the Hartree-Fock (HF) approach in the coordinate space with a Skyrme effective interaction (see  for details). The analysis of coexistence and stability of different nuclear phases requires construction of an adiabatic path between different configurations. It has lead us to abandon the Wigner-Seitz approximation for the Coulomb interaction, solving instead the Poisson equation for the electric potential within the box with periodic boundary conditions, the total charge of the box being zero. In our approach the electrons are assumed to form a uniform relativistic gas where screening effects are neglected.

Summarizing, we minimize an energy functional for the neutron-proton-electron matter within the box. The energy functional describes both the liquid drop part of the energy and the shell effects. Namely:

$$E_{tot}(\rho _p,\rho _n,\rho _e)=E_{TF}(\rho _p,\rho _n,\rho _e)+E_{shell}^{in}(\rho _p^{in},\rho _n^{in})+E_{shell}^{out}(\rho _n^{out},\rho _p^{out}),$$
(5)
where $`\rho _n,\rho _p,\rho _e`$ denote the density of neutrons, protons and electrons, respectively. The charge neutrality condition implies $`\rho _e=\rho _p`$. The total nuclear density $`\rho =\rho _n+\rho _p`$. In the above expression $`E_{TF}`$ represents the Thomas-Fermi energy of the system obtained by constructing the energy functional for a Skyrme interaction. We have used the SLy4 force  which give a correct description of neutron matter at low densities .

The shell energy consists of two parts. The one denoted by $`E_{shell}^{in}`$ is determined by the density of neutrons $`\rho _n^{in}`$ and protons $`\rho _p^{in}`$ inside the nucleus, whereas $`E_{shell}^{out}`$ is a function of unbound neutrons $`\rho _n^{out}`$ and protons $`\rho _p^{out}`$. Since $`\rho _p^{out}=0`$ for $`\rho <0.08`$ fm-3 , one can neglegt the contribution to the energy due to the unbound protons.

The HF calculations confirm the importance of the shell effects associated with unbound neutrons . In particular, they show that the shell effects associated with unbound neutrons may reverse the phase transition order predicted by the liquid drop based approaches. Moreover, the number of phase transitions increases since the same phase may appear for various density ranges. Our results suggest also that the number of different phases present in the crust is larger since one has to take into account various lattice geometries which may compete. This important fact was disregarded in earlier investigations since in liquid drop based approaches the system favours only one lattice type for a given nuclear shape. However it is no longer the case when one considers the neutron shell energy which is very sensitive to the spatial order in the system.

In Fig.1, we show several nuclear configurations that can be formed at subnuclear densities. One notices the presence of different types of crystal lattices for the same nuclear shapes. The integrated densities are defined by:

$`\rho _i^z(x,y)`$ $`=`$ $`{\displaystyle _0^d}\rho _i(x,y,z)𝑑z,`$ (6)
$`\rho _i^y(x,z)`$ $`=`$ $`{\displaystyle _0^d}\rho _i(x,y,z)𝑑y,`$ (7)
where $`i=n,p`$ and $`d`$ denotes the size of the box.

Once a phase is formed, there is a positional order maintained by the Coulomb repulsion between spherical nuclei, rods or slabs. The Coulomb energy is a smooth function of the void displacement but not the shell energy. Since several orbits contribute to the shell effects (except for the slab–like phase), the displacement of a single bubble–like or rod–like void from its equilibrium position in the lattice will give rise to interference effects .

The main conclusion of the present work is that the amplitude of the shell effects is comparable to the energy differences between the various phases determined in simpler liquid drop type models. Our results suggest that the inhomogeneous phase has a complicated structure, maybe completely disordered, with several types of coexisting shapes. The few results highlighted here shows the richness of these systems where both static and dynamic properties are challenging to describe.

ACKNOWLEDGMENTS

This research was supported in part by the Polish Committee for Scientific Research (KBN) and the Wallonie/Brussels-Poland integrated action program. Authors would like to thank Hubert Flocard for providing us with the numerical code solving the Poisson equation. We also appreciate many discussions with A. Wirzba, Y. Yu, S.A. Chin and H. Forbert.


text:
# Spatial field correlation, the building block of mesoscopic fluctuations

## Abstract

The absence of self averaging in mesoscopic systems is a consequence of long-range intensity correlation. Microwave measurements suggest and diagrammatic calculations confirm that the correlation function of the normalized intensity with displacement of the source and detector, $`\mathrm{\Delta }R`$ and $`\mathrm{\Delta }r`$, respectively, can be expressed as the sum of three terms, with distinctive spatial dependences. Each term involves only the sum or the product of the square of the field correlation function, $`FF_E^2`$. The leading-order term is the product, the next term is proportional to the sum. The third term is proportional to $`[F(\mathrm{\Delta }R)F(\mathrm{\Delta }r)+[F(\mathrm{\Delta }R)+F(\mathrm{\Delta }r)]+1]`$.

Short-range correlation in waves transmitted through random media is manifest in the intensity speckle pattern. The leading contribution, $`C_1`$, to the cumulant correlation function $`C`$ of intensity normalized to its ensemble average on the output surface of the sample is given by the square of the field correlation function, $`C_1=F_E^2`$. Neglecting internal reflection from the surface, its dependence upon displacement $`\mathrm{\Delta }r`$ on the output surface is given by $`C_1(\mathrm{\Delta }r)=F_E^2(\mathrm{\Delta }r)F(\mathrm{\Delta }r)=(\mathrm{sin}k\mathrm{\Delta }r/k\mathrm{\Delta }r)^2\mathrm{exp}(\mathrm{\Delta }r/\mathrm{}_s)`$, where $`k`$ is the wave vector and $`\mathrm{}_s`$ is the scattering mean free path Shapiro . This term dominates intensity fluctuations. Defining a correlation length, $`\delta r`$, as the first zero of $`C_1`$ gives $`\delta r=\pi /k=\lambda /2`$. The intensity is correlated far beyond $`\delta r`$ as a result of scattering within the medium Stephen ; Feng88 ; mello ; GenackPRL90 ; kk ; dutch  so that intensity values in remote speckle spots are not statistically independent. This gives rise to two additional contributions to $`C`$, which can therefore be expressed as, $`C=C_1+C_2+C_3`$ Feng88 ; Shapiro99 , and leads to greatly enhanced mesoscopic fluctuations AltshulerAltshuler . The $`C_2`$ term produces an enhancement in total transmission fluctuations over that given by the field factorization approximation by a factor of $`L/\mathrm{}`$, where $`\mathrm{}`$ is the transport mean free path and $`L`$ is the sample length Stephen ; dutch . The $`C_3`$ term is the source of universal conductance fluctuations, which are enhanced by a factor of $`(L/\mathrm{})^2`$ AltshulerAltshuler ; Laibowitz . The magnitude of $`C_1`$ at a point is unity as a result of normalization, whereas the magnitudes of $`C_2`$ and $`C_3`$ are expansions in $`1/g`$ with leading terms of order $`1/g`$ and $`1/g^2`$ Feng88 , respectively, where $`g`$ is the dimensionless conductance. Since the onset of localization is at $`g=1`$ gang4 , the two terms beyond the field factorization approximation for $`C`$, reflect the approach to localization.

In this Letter, we use microwave measurements and diagrammatic calculations to show that each of the contributions to $`C`$ may be expressed in terms of the square of the field correlation function with regard to displacements of the source, $`\mathrm{\Delta }R`$ and detector, $`\mathrm{\Delta }r`$. The $`C_1`$ term is $`F(\mathrm{\Delta }R)F(\mathrm{\Delta }r)`$, the $`C_2`$ term is proportional to $`[F(\mathrm{\Delta }R)+F(\mathrm{\Delta }r)]`$, while the $`C_3`$ term is proportional to $`[F(\mathrm{\Delta }R)F(\mathrm{\Delta }r)+[F(\mathrm{\Delta }R)+F(\mathrm{\Delta }r)]+1]`$. When the intensity correlation is considered at a shifted frequency, $`\mathrm{\Delta }\nu `$, the full correlation function remains a sum of three terms each being a product of the corresponding terms in C, thus, $`C_i=A_i(\mathrm{\Delta }\nu )C_i(\mathrm{\Delta }r,\mathrm{\Delta }R)(i=1,2,3)`$. Absorption alters the magnitudes of $`C_2`$ and $`C_3`$, but it does not change the spatial structure of these terms.

Initial measurements of angular intensity correlation, carried out in the far field of weakly scattering media, gave $`C`$, which was essentially equal to $`C_1`$ Freund88 ; Li94 . Recently, measurements of the spatial correlation of the field on the sample surface have yielded the $`C_1`$ contribution directly PRE00 . Measurements of intensity correlation between points on the sample surface and the interior of the sample on a scale greater than the wavelength have allowed the observation of $`C_2`$ GenackPRL90 . In addition measurements have been made of the spectral correlation of the $`C_1`$ Genack87  and $`C_2`$ terms Garcia89 ; Ad90 ; GGPS  and of the temporal correlation of the $`C_1`$ Maret87 , $`C_2`$ FSPRB97  and $`C_3`$ FSPRL98  terms in colloidal samples. However, the variation of $`C`$ on a subwavelength scale, as well as the structure of correlation with displacement of both the source and detector, have not been reported previously. The present measurements and calculations allow us to discern the structure of intensity correlation and to relate it to the correlation of the underlying field.

Measurements are made in a disordered dielectric sample contained within a reflecting tube. Radiation of frequency $`\nu `$ emitted by a source antenna at $`\stackrel{}{R}`$ at one end of the tube is detected at a point $`\stackrel{}{r}`$ at the other end. We denote the intensity at $`\stackrel{}{r}`$ due to a source at $`\stackrel{}{R}`$ by $`I_\nu (\stackrel{}{r},\stackrel{}{R})`$ and consider the normalized cumulant correlation function,

$$C(\mathrm{\Delta }r,\mathrm{\Delta }R)\delta I_\nu (\stackrel{}{r},\stackrel{}{R})\delta I_\nu ^{}(\stackrel{}{r}^{},\stackrel{}{R}^{})/I_\nu (\stackrel{}{r},\stackrel{}{R})I_\nu ^{}(\stackrel{}{r}^{},\stackrel{}{R}^{}),$$
(1)
where $`\delta I_\nu `$ is the deviation of the intensity from its ensemble average value, $`\mathrm{\Delta }r=|\stackrel{}{r}\stackrel{}{r}^{}|`$ and $`\mathrm{\Delta }R=|\stackrel{}{R}\stackrel{}{R}^{}|`$ are the displacements across the output and input surfaces, respectively, and $`\mathrm{}`$ denotes the average over an ensemble of random realizations. The leading contribution to $`C`$ obtained by factorizing the fields is Stephen ; Shapiro ; PRE00 ; FreundEliyahou :

$$C_1(\mathrm{\Delta }r,\mathrm{\Delta }R)=|E_\nu (\stackrel{}{r},\stackrel{}{R})E_\nu ^{}^{}(\stackrel{}{r}^{},\stackrel{}{R}^{})|^2/I_\nu (\stackrel{}{r},\stackrel{}{R})I_\nu ^{}(\stackrel{}{r}^{},\stackrel{}{R}^{}).$$
(2)
The samples studied are random mixtures of 1.27-cm-diameter polystyrene spheres at a volume filling fraction of 0.52. They are contained inside a 100 cm long copper tube with a diameter of 7.5 cm capped with thin Plexiglas face plates. Field and intensity spectra are obtained in the frequency range 16.8 - 17.8 GHz with a step size of 625 kHz using a vector network analyzer. Measurements are taken in an ensemble of 690 random samples by rotating the tube to create new configurations of the spheres after each set of spectra are taken. Field spectra are obtained in each sample realization by translating an antenna detector along a line to each of 50 locations separated by 1.06 mm on the output surface for each of two fixed antenna sources at the incident surface, which are separated by $`\mathrm{\Delta }R=3\pm 0.1\mathrm{cm}d`$. Field spectra are also taken by translating the source along a line for each of two fixed antenna detectors at the output surface, separated by $`\mathrm{\Delta }r=d`$. At the separation $`d`$, the short-range term has largely decayed. The antennas are aligned perpendicular to the line of translation. Intensity spectra are obtained by squaring the field spectra.

The spatial variations of $`C(\mathrm{\Delta }r,0)`$ and $`C_1(\mathrm{\Delta }r,0)`$ are shown in Fig. 1. The $`C_1`$ contribution is directly obtained by squaring the field correlation function, shown in the inset of Fig. 1. Subtracting $`C_1`$ from the full intensity correlation function gives the difference $`[CC_1]`$, shown in Fig. 2 for a single source ($`\mathrm{\Delta }R=0`$) and for two sources separated by $`\mathrm{\Delta }R=d`$. This difference includes the terms beyond the field factorization approximation.

Measurements of $`C(0,\mathrm{\Delta }R)`$ and $`C_1(0,\mathrm{\Delta }R)`$, i.e. for a fixed detector and a scanning source, have also been performed. Within experimental error, $`C(\mathrm{\Delta }r,0)`$ and $`C(0,\mathrm{\Delta }R)`$ were found to be identical functions of their respective arguments. The same is true for $`C_1`$. Similarly, we find that plots of $`[CC_1]`$ versus $`\mathrm{\Delta }R`$ with $`\mathrm{\Delta }r=0`$ and $`\mathrm{\Delta }r=d`$ are nearly the same as those shown in Fig. 2. Thus $`\mathrm{\Delta }R`$ and $`\mathrm{\Delta }r`$ can be interchanged as required by reciprocity.

Measurements of $`C_1(\mathrm{\Delta }r,\mathrm{\Delta }R)`$ for $`\mathrm{\Delta }R=0`$ and $`\mathrm{\Delta }R=d`$ are presented in Fig. 3. Within the noise level of $`10^4`$, the two functions have the same variation with $`\mathrm{\Delta }r`$, $`C_1(\mathrm{\Delta }r,d)=2\times 10^3C_1(\mathrm{\Delta }r,0)`$. This numerical factor, $`C_1(0,d)`$, is equal to the value of $`C_1(d,0)`$ within the uncertainty in $`d`$. This result, taken together with the aforementioned symmetry with respect to interchanging $`\mathrm{\Delta }r`$ and $`\mathrm{\Delta }R`$, suggests that $`C_1`$ can be written as the product of two identical functions, $`C_1(\mathrm{\Delta }r,\mathrm{\Delta }R)=F(\mathrm{\Delta }R)F(\mathrm{\Delta }r)`$.

We now examine $`[CC_1](\mathrm{\Delta }r,\mathrm{\Delta }R)`$, which is do minated by $`C_2`$ in our sample. This function is seen in Fig. 2 to fall to nearly one half its value when either $`\mathrm{\Delta }R`$ or $`\mathrm{\Delta }r`$ increase beyond $`\delta r`$ when there is no displacement of the other variable. This shows that $`C_2`$ is given by the addition of two equal terms. The comparison of the short-range variation of $`[CC_1]`$ with $`C_1`$ in Fig. 2 suggests that the additive form factors are identical to $`F`$, so that the dominant contribution to $`CC_1`$ is proportional to $`[F(\mathrm{\Delta }R)+F(\mathrm{\Delta }r)]`$. This would imply that $`C_2(\mathrm{\Delta }r)`$ approach a constant value for $`\mathrm{\Delta }r>\delta r`$, whereas the measurement of $`[CC_1](\mathrm{\Delta }r,d)`$ is seen in Fig. 2 to fall slightly with increasing displacement. This could be the consequence of a slight departure from a quasi-1D geometry at the output face of the sample. There, the average intensity is slightly larger at the center than at the edges since the wave can spread beyond the cross section of the tube. Notwithstanding this effect, our experimental results suggest that both $`C_1`$ and $`C_2`$ can be expressed in terms of a single form factor $`F(x)`$, where $`x`$ stands for either $`\mathrm{\Delta }r`$ or $`\mathrm{\Delta }R`$. $`C_1`$ and $`C_2`$ contain, respectively, the product and the sum of two form factors. In Fig. 2, the correlation function $`[CC_1](\mathrm{\Delta }r,d)`$ is seen to lie above the dotted curve, which is proportional to $`C_1`$. This suggests a constant contribution to $`C`$. For $`\mathrm{\Delta }r>30`$ mm, the correlation function becomes negative, but here the noise becomes larger than the signal because of the reduced number of pairs of points with increasing $`\mathrm{\Delta }r`$. Such a long-range correlation for large values of $`\mathrm{\Delta }R`$ may be part of $`C_3`$.

The structure of the joint spatial and frequency dependences of $`C_1`$ and $`C_2`$, is obtained from measurements of the correlation functions $`C_1(\mathrm{\Delta }\nu ,\mathrm{\Delta }r)`$ and $`[CC_1](\mathrm{\Delta }\nu ,\mathrm{\Delta }r)`$ for $`\mathrm{\Delta }R=0`$, shown in Figs. 4a and 4b, respectively. The semilog representations in Fig. 4 show that, within the limits set by the noise level, $`C_i`$ have the same frequency dependence for any $`\mathrm{\Delta }r`$, while $`C_i`$ have the same spatial dependence for any $`\mathrm{\Delta }\nu `$ for $`i=1,2`$. Thus their spatial and spectral variations for a single source are given by $`C_i(\mathrm{\Delta }\nu ,\mathrm{\Delta }r)=A_i(\mathrm{\Delta }\nu )C_i(\mathrm{\Delta }r)`$. The noise level found in $`C_1`$ is low compared to that in $`C_2`$ because the field correlation function, $`F_E`$, is computed and then squared to obtain $`C_1`$, giving a signal to noise ratio which is the square of that for the field correlation function. The form of the intensity correlation function suggested by experiment is borne out in diagrammatic calculations, which are briefly summarized below.

The three terms in $`C`$ may be represented diagrammatically. The diagram corresponding to the $`C_1`$ term describes two non-interacting diffusons attached to pairs, $`GG^{}`$, of averaged Green’s functions dutch ; PniniBook . This diagram factorizes into a product of two field functions:

$$E_\nu (\stackrel{}{r},\stackrel{}{R})E_\nu ^{}(\stackrel{}{r}^{},\stackrel{}{R}^{})=d^3r_1d^3r_2G_\nu (\stackrel{}{r},\stackrel{}{r}_1)G_\nu ^{}^{}(\stackrel{}{r}^{},\stackrel{}{r}_1)T_{\nu \nu ^{}}(\stackrel{}{r}_1,\stackrel{}{r}_2)G_\nu (\stackrel{}{r}_2,\stackrel{}{R})G_\nu ^{}^{}(\stackrel{}{r}_2,\stackrel{}{R}^{}),$$
(3)
where $`T_{\nu \nu ^{}}(\stackrel{}{r}_1,\stackrel{}{r}_2)`$ denotes the diffusion ladder and integration is performed over $`\stackrel{}{r}_1`$, $`\stackrel{}{r}_2`$ inside the tube. For the quasi one-dimensional geometry, the diffuson is independent of its transverse coordinates, whereas the Green’s functions decay rapidly on a scale of the mean free path. Taking $`\stackrel{}{R}=\stackrel{}{R}^{}`$ and $`\nu =\nu ^{}`$, we obtain

$$E_\nu (\stackrel{}{r},\stackrel{}{R})E_\nu (\stackrel{}{r}^{},\stackrel{}{R})=\left(\frac{4\pi }{\mathrm{}}\right)d^3r_1G_\nu (\stackrel{}{r},\stackrel{}{r}_1)G_\nu ^{}(\stackrel{}{r}^{},\stackrel{}{r}_1)I_\nu (\stackrel{}{r}_1,\stackrel{}{R})F_E(\mathrm{\Delta }r),$$
(4)
where $`I_\nu (\stackrel{}{r},\stackrel{}{R})`$ is the intensity at $`\stackrel{}{r}`$, normalized to its average value at the output face of the tube. Thus, $`C_1(\mathrm{\Delta }r,\mathrm{\Delta }R)=F_E^2(\mathrm{\Delta }r)F_E^2(\mathrm{\Delta }R)F(\mathrm{\Delta }R)F(\mathrm{\Delta }r)`$.

The diagrams corresponding to the $`C_2`$ and $`C_3`$ terms describe two incoming and two outgoing diffusons which interact in the bulk of the medium. In these diagrams, each pair $`GG^{}`$ of external Green’s functions contributes a spatial form factor $`F_E`$ as in Eq. (4). These give

$`C`$ $`=`$ $`C_1+C_2+C_3`$ (5)
$`=`$ $`A_1(\mathrm{\Delta }\nu ,\alpha )F(\mathrm{\Delta }R)F(\mathrm{\Delta }r)+{\displaystyle \frac{2}{3g}}A_2(\mathrm{\Delta }\nu ,\alpha )\left[F(\mathrm{\Delta }R)+F(\mathrm{\Delta }r)\right]`$
$`+`$ $`{\displaystyle \frac{2}{15g^2}}A_3(\mathrm{\Delta }\nu ,\alpha )\left[1+F(\mathrm{\Delta }R)+F(\mathrm{\Delta }r)+F(\mathrm{\Delta }R)F(\mathrm{\Delta }r)\right],`$
where $`g`$ is the average conductance of the sample. The coefficients $`A_i(i=1,2,3)`$ depend on the absorption coefficient $`\alpha `$ and the frequency shift $`\mathrm{\Delta }\nu `$. The structure of Eq. (5) is similar to that of the correlation in transmission, obtained in the multichannel formalism Feng88 ; mello ; dutch . In the present case, however, all terms are described by a single spatial form factor. The field factorization term, $`A_1(\mathrm{\Delta }\nu =0,\alpha )`$ is unity and independent of absorption by definition, while $`A_2(\mathrm{\Delta }\nu =0,\alpha )`$ is given by kk ; ps ; brewer :

$$A_2(\mathrm{\Delta }\nu =0,\alpha )=\frac{3}{16\alpha }\left[\frac{\mathrm{sinh}2\alpha 2\alpha (2\mathrm{cosh}2\alpha )}{\mathrm{sinh}^2\alpha }\right].$$
(6)
The coefficient $`A_3(\alpha )`$ depends weakly on $`\alpha `$ and its limiting values are $`A_3(0)=1,A_3(\mathrm{})=15/16`$.

For our samples, $`g7`$ and $`\alpha 3`$ PRE00 . From the measurement of $`[CC_1]`$ at $`\mathrm{\Delta }r=0`$, $`\frac{2}{3g}A_2=0.076`$ with $`A_2=0.87`$. This is in agreement with Eq. (6) and calculated corrections due to localization effects kk ; GGPS ; ps ; brewer ; PniniBook . Using the measured $`C_1(\mathrm{\Delta }r)`$ as the functional form $`F(\mathrm{\Delta }r)`$, following Eq. (5) and neglecting $`C_3`$, we obtain a good fit of the spatial structure of the intensity correlation function, as shown in Fig. 1.

These considerations have applications to present efforts to enhance the capacity of wireless communication by utilizing multiple antennas to detect the multiply scattered field Moustakas ; PhysToday . Antenna separation should be larger than $`\delta r`$ and the number of statistically independent antennas equals the inverse of the degree of long-range intensity correlation, $`3g/2`$.

In conclusion, we have uncovered the connection between the field and intensity correlation functions in the spatial structure of the three contributions to $`C`$. In contrast to the case of angular correlation Feng88 , intensity correlation can be expressed in terms of a single form factor obtained from the field correlation function. We have demonstrated the multiplicative character of $`C_1`$ and the additive character of $`C_2`$. Calculations predict a mixed character for $`C_3`$, which includes a multiplicative, an additive, and a constant term of equal amplitude. We observe the infinite-range component of $`C_3`$ in the residual correlation when both the source and detector are displaced by more than the correlation length. Determining the proper breakup of $`C`$ into its components is of particular importance when considering simultaneous variations in space, time and frequency. Each term is a product of the corresponding $`C_1`$, $`C_2`$ and $`C_3`$ correlation function in the appropriate variables.

We thank A.A. Chabanov for valuable discussions. Support from the National Science Foundation (DMR 9973959), Army Research Office (DAAD 190010362), the United States-Israel Binational Science Foundation (BSF) and the Groupement de Recherche PRIMA are gratefully acknowledged.