text:
# Spectral scrambling in Coulomb-blockade quantum dots

## Abstract

We study the fluctuations of an energy level as a function of the number of electrons $`m`$ added to a Coulomb-blockade quantum dot. A microscopic calculation in the limit of Koopmans’ theorem predicts that the standard deviation of these fluctuations behaves as $`\sqrt{m}`$ in the absence of surface charge but is linear in $`m`$ when the effect of a surface charge in a finite geometry is included. The microscopic results are compared to a parametric random-matrix approach. We estimate the number of electrons it takes to scramble the spectrum completely in terms of the interaction strength, the dimensionless Thouless conductance, and the symmetry class.

The simplest model for describing a quantum dot in the Coulomb-blockade regime is the constant interaction (CI) model , in which the electrons occupy single-particle levels and the Coulomb interaction is taken as an average electrostatic energy that depends only on the number of electrons. When the single-particle dynamics in the dot are chaotic (or diffusive), one can apply random matrix theory (RMT) to describe the statistical properties of the single-particle wave functions and energies within an energy window whose width is the Thouless energy. RMT was successful in describing the conductance peak height distributions and their sensitivity to time-reversal symmetry . However, other measured statistics, most notably the peak spacing distribution , indicated that it is necessary to include interaction effects beyond the simple CI model . The best way to take into account interactions while retaining a single-particle framework is the Hartree-Fock (HF) approximation, which has been used to explain some of the observed features of the peak spacing statistics .

The HF single-particle wave functions and energies are calculated self-consistently and therefore can change as electrons are added to the dot. In the statistical regime (i.e., in a chaotic or diffusive dot), this phenomenon is called scrambling. Scrambling was observed in the decay of correlations between the $`m`$-th excited state in the dot and the ground-state of a dot with an additional $`m`$ electrons , and indirectly through the saturation of the peak-to-peak correlations as a function of temperature . A phenomenological way to describe scrambling is to consider a discrete set of Hamiltonians (corresponding to the different number of electrons in the dot) that are random but have the correct symmetries. Such a set is known as a discrete Gaussian process (GP), and can be embedded in a continuous GP, i.e., random matrices that depend on a continuous parameter . This parametric approach leads to a nearly Gaussian peak spacing distribution , and explains the saturation of the number of correlated peaks versus temperature . While the parametric approach is appealing in its simplicity, it is not clear how well it describes features obtained in a microscopic approach. For example, in the parametric approach, the fluctuation standard deviation (FSD) of a typical level upon the addition of $`m`$ electrons into the dot is proportional to $`m`$. One can also evaluate the FSD of a level in the microscopic approach through the fluctuations of the diagonal interaction matrix elements (in the limit where the single-particle wave functions do not change with the addition of electrons, hence the off-diagonal matrix elements can be neglected). When only the bulk screened Coulomb interaction is accounted, we find that the associated FSD behaves like $`\sqrt{m}`$, unlike the standard parametric approach. However, in the finite geometry of a quantum dot, the level fluctuation may be dominated by the variation of the mean-field potential arising from surface charge. We shall show that in this case the FSD is linear in $`m`$ (for $`1mg`$, where $`g`$ is the dimensionless Thouless conductance), in overall agreement with parametric RMT. The main results of the microscopic approach are summarized in Eqs. (12), (19) and (20) in the absence of surface charge and Eqs. (25), (30) and (31) in the presence of surface charge. Extrapolating our results to larger values of $`m`$, we also estimate the dependence of the number of added electrons $`m_c`$ required for complete scrambling on $`g`$ and the interaction strength both in the absence and presence of surface charge (Eqs. (21) and (33), respectively).

Scrambling implies variation of both eigenstates and energy levels. However, in this work we approach the problem from the limit where the wave functions do not change (Koopmans’ theorem ), and we only estimate the FSD of energy levels. In the parametric approach this limit is valid when the change in the parameter upon the addition of $`m`$ electrons is small compared to the mean parametric distance between avoided crossings. Complete scrambling occurs when the energy fluctuations become comparable to $`\mathrm{\Delta }`$ (the mean-level spacing in the middle of the band).

We first discuss the parametric approach. The variation of the single-particle energies and eigenfunctions (e.g., in a mean-field approximation) with the addition of electrons into the dot is described by a parametric variation of the Hamiltonian . We denote by $`H(x_𝒩)`$ the effective single-particle Hamiltonian of the dot with $`𝒩`$ electrons. We assume that the dot is either diffusive or ballistic with chaotic dynamics, and that the statistical properties of $`H(x_𝒩)`$ are not modified by the interactions. Restricting ourselves to the universal regime, i.e., to $`g`$ levels in the vicinity of the Fermi energy , we assume that $`H(x_𝒩)`$ belongs to one of the Gaussian ensembles of random matrices whose symmetry class $`\beta `$ is independent of $`x_𝒩`$. The sequence of Hamiltonians $`H(x_𝒩)`$ forms a discrete GP that can be embedded in a continuous GP $`H(x)`$ . A simple GP is given by  $`H(x)=\mathrm{cos}xH_1+\mathrm{sin}xH_2`$, where $`H_1`$ and $`H_2`$ are $`N\times N`$ uncorrelated random matrices chosen from the Gaussian ensemble of symmetry class $`\beta `$: $`P(H)e^{\frac{\beta }{2a^2}\mathrm{Tr}H^2}`$. We choose $`a=(2/N)^{1/2}\zeta `$ so that the average level density (a semicircle) has a constant band width of $`2a\sqrt{2N}=4\zeta `$ and the mean-level spacing in the middle of the spectrum is $`\mathrm{\Delta }=\pi a/\sqrt{2N}=\pi \zeta /N`$. The average distance between avoided crossings is given by the inverse of the rms level velocity $`\delta x_c=\mathrm{\Delta }[\overline{(ϵ_\alpha /x)^2}]^{1/2}=\pi (\beta /2N)^{1/2}`$, where $`ϵ_\alpha `$ is an energy level. This distance is larger for the GUE by a factor of $`\sqrt{2}`$ compared with the GOE. Two-point parametric correlators become universal when the parameter $`x`$ is measured in units of $`\delta x_c`$; i.e., as a function of a scaled parameter $`\overline{x}x/\delta x_c`$.

The energy levels $`ϵ_\alpha `$ scramble as the parameter $`\overline{x}`$ changes: $`\delta ϵ_\alpha =ϵ_\alpha (\overline{x}+\delta \overline{x})ϵ_\alpha (\overline{x})`$. The variance of $`\delta ϵ_\alpha `$ is estimated in the limit $`\delta \overline{x}1`$ to be $`\sigma ^2(\delta ϵ_\alpha )=\mathrm{\Delta }^2(\delta \overline{x})^2`$ using first order perturbation theory in $`\delta \overline{x}`$ (i.e., ignoring the change of the single-particle wave function as $`\overline{x}\overline{x}+\delta \overline{x}`$).

In the parametric approach it is assumed that $`\overline{x}`$ changes by $`\delta \overline{x}_1`$ upon the addition of one electron into the dot (independently of the number of electrons $`𝒩`$). Thus for the addition of $`m`$ electrons $`\delta \overline{x}_m=m\delta \overline{x}_1`$, and as long as $`\delta \overline{x}_m1`$, we can still use first order perturbation theory to find

$$\sigma (\delta ϵ_\alpha ^{(m)})=\mathrm{\Delta }m\delta \overline{x}_1=m\sigma (\delta ϵ_\alpha ^{(1)}),$$
(1)
where $`\delta ϵ_\alpha ^{(m)}ϵ_\alpha ^{(𝒩+m)}ϵ_\alpha ^{(𝒩)}`$ ($`ϵ_\alpha ^{(𝒩)}`$ is the energy of level $`\alpha `$ in a dot with $`𝒩`$ electrons).

To relate the parametric approach to a microscopic mean-field approach, we describe the effect that adding electrons has on a particular single-particle HF level. The Hamiltonian of the dot is

$$H=\underset{ij}{}h_{ij}^{(0)}a_i^{}a_j+\frac{1}{4}\underset{ijkl}{}v_{ijkl}^Aa_i^{}a_j^{}a_la_k,$$
(2)
where $`h^{(0)}`$ is the one-body part (accounting for the disorder or the chaotic confining potential) and $`v_{ij;kl}^Aij|v|klij|v|lk`$ are the antisymmetrized matrix elements of the two-body interaction. The HF Hamiltonian in the self-consistent basis for $`𝒩`$ electrons is given by $`h_{\alpha \gamma }^{(𝒩)}=ϵ_{\alpha \gamma }^{(0)}+_{\delta =1}^𝒩v_{\alpha \delta ;\gamma \delta }^A`$ (where the sum is over the lowest $`𝒩`$ occupied levels), and the eigenvalues of $`h_{\alpha \gamma }^{(𝒩)}`$ are the HF single-particle energies $`ϵ_\alpha ^{(𝒩)}`$. The diagonal part of the Hamiltonian (2) in the HF basis is

$$H_{\mathrm{diagonal}}=\underset{\alpha }{}ϵ_{\alpha \alpha }^{(0)}\widehat{n}_\alpha +\frac{1}{2}\underset{\alpha ,\gamma }{}v_{\alpha \gamma }^A\widehat{n}_\alpha \widehat{n}_\gamma ,$$
(3)
where $`\widehat{n}_\alpha `$ is the number operator for single-particle state $`\alpha `$ and $`v_{\alpha \gamma }^Av_{\alpha \gamma ;\alpha \gamma }^A`$.

In Koopmans’ limit, the single-particle wave functions do not change when an electron is added to the dot. Consequently, the change in the single-particle energy $`ϵ_\alpha `$ when an electron is added to the $`𝒩+1`$-st level is given by a diagonal interaction matrix element

$$\delta ϵ_\alpha ^{(1)}ϵ_\alpha ^{(𝒩+1)}ϵ_\alpha ^{(𝒩)}v_{\alpha 𝒩+1}^A.$$
(4)
Assuming the HF wave functions satisfy similar statistics as in the non-interacting case, we can estimate the variance of the matrix element in (4). The two-body interaction used in the HF approximation is often an effective interaction, e.g., an RPA screened interaction . Alternatively, one can employ a short-range dressed interaction  $`v(𝒓𝒓^{})=\lambda \mathrm{\Delta }V\delta (𝒓𝒓^{})`$, where $`\lambda `$ is the interaction strength and $`V`$ is the system’s volume. A diagonal interaction matrix element is then given by $`v_{\alpha \gamma }=\lambda \mathrm{\Delta }V𝑑𝒓\psi _\alpha ^{}(𝒓)\psi _\gamma ^{}(𝒓)\psi _\alpha (𝒓)\psi _\gamma (𝒓)`$, and its variance

$`\sigma ^2(v_{\alpha \gamma })`$ $`=\lambda ^2\mathrm{\Delta }^2V^2{\displaystyle }d𝒓_1{\displaystyle }d𝒓_2[\psi _\alpha ^{}(𝒓_1)\psi _\gamma ^{}(𝒓_1)\psi _\alpha (𝒓_1)\psi _\gamma (𝒓_1)\psi _\alpha (𝒓_2)\psi _\gamma (𝒓_2)\psi _\alpha ^{}(𝒓_2)\psi _\gamma ^{}(𝒓_2)`$ (6)
$`\psi _\alpha ^{}(𝒓_1)\psi _\gamma ^{}(𝒓_1)\psi _\alpha (𝒓_1)\psi _\gamma (𝒓_1)\psi _\alpha (𝒓_2)\psi _\gamma (𝒓_2)\psi _\alpha ^{}(𝒓_2)\psi _\gamma ^{}(𝒓_2)].`$
Only the connected part of the ensemble average over the product of eight wave functions contributes to the r.h.s. of (6). In the following we restrict our calculations to the GUE symmetry. All possible pairwise and quadruplet-wise contractions should be taken into account, and they can be regrouped to give

$`\sigma ^2(v_{\alpha \gamma })`$ $`=\lambda ^2\mathrm{\Delta }^2V^2{\displaystyle }d𝒓_1{\displaystyle }d𝒓_2\{[\psi _\alpha ^{}(𝒓_1)\psi _\alpha (𝒓_1)\psi _\alpha (𝒓_2)\psi _\alpha ^{}(𝒓_2)\psi _\gamma ^{}(𝒓_1)\psi _\gamma (𝒓_1)\psi _\gamma (𝒓_2)\psi _\gamma ^{}(𝒓_2){\displaystyle \frac{1}{V^4}}]`$ (7)
$`+`$ $`\left[\psi _\alpha ^{}(𝒓_1)\psi _\alpha (𝒓_1)\psi _\gamma (𝒓_2)\psi _\gamma ^{}(𝒓_2)\psi _\gamma ^{}(𝒓_1)\psi _\gamma (𝒓_1)\psi _\alpha (𝒓_2)\psi _\alpha ^{}(𝒓_2){\displaystyle \frac{1}{V^4}}\right]`$ (8)
$`+`$ $`\psi _\alpha ^{}(𝒓_1)\psi _\alpha (𝒓_2)\psi _\gamma (𝒓_1)\psi _\gamma ^{}(𝒓_2)\psi _\gamma ^{}(𝒓_1)\psi _\gamma (𝒓_2)\psi _\alpha (𝒓_1)\psi _\alpha ^{}(𝒓_2)\},`$ (9)
where the $`1/V^4`$ factors are subtracted to avoid double counting of pairwise contractions (here $`\mathrm{}`$ denotes an ensemble average and not just the connected part).

Wave-function correlations appearing in Eq. (7) were calculated in Ref. . Including contributions of order $`1/g^2`$ and assuming that the distance between levels $`\alpha `$ and $`\gamma `$ is less than the Thouless energy, we have

$`\sigma ^2(v_{\alpha \gamma })`$ $`={\displaystyle \frac{\lambda ^2\mathrm{\Delta }^2}{V^2}}{\displaystyle }d𝒓_1{\displaystyle }d𝒓_2\{[1+k_d(r)(1+\mathrm{\Pi }(𝒓_1,𝒓_2))+\mathrm{\Pi }(𝒓_1,𝒓_2)+b(𝒓_1,𝒓_2)]^21`$ (10)
$`+`$ $`[1+k_d(r)\mathrm{\Pi }(𝒓_1,𝒓_2)+c(𝒓_1,𝒓_2)]^21+[k_d(r)+\mathrm{\Pi }(𝒓_1,𝒓_2)]^2\}.`$ (11)
The function $`k_d(r)`$ ($`𝒓=𝒓_1𝒓_2`$) is given by $`k_d(r)(\pi \nu )^2\mathrm{Im}G^R(𝒓)^2`$ (where $`G^R`$ is the retarded Green function and $`\nu `$ the electron density of states per unit volume) and describes the short-range correlations. In 2D $`k_d(r)=\mathrm{exp}(r/\mathrm{})J_0^2(k_Fr)`$ (where $`\mathrm{}`$ is the mean free path and $`J_0`$ is a Bessel function), while in 3D $`k_d(r)=\mathrm{exp}(r/\mathrm{})(\mathrm{sin}k_Fr/k_Fr)^2`$. $`\mathrm{\Pi }(𝒓_1,𝒓_2)=(\pi \nu )^1_{𝑸0}\varphi _𝑸(𝒓_1)\varphi _𝑸(𝒓_2)/D𝑸^2`$ is the long-range diffuson propagator where $`\varphi _𝑸`$ is the eigenfunction of the diffusion operator corresponding to the eigenvalue $`D𝑸^2`$ ($`D`$ is the diffusion constant). The functions $`b,c`$ in Eq. (10) are of order $`1/g^2`$ ($`c\mathrm{\Pi }^2(𝒓_1,𝒓_2)/2`$ for $`\mathrm{}rL`$  where $`L`$ is the linear size of the dot). The main contributions to (10) arise from the long-range terms $`\mathrm{\Pi }^2`$ ($`𝑑𝒓_1\mathrm{\Pi }(𝒓_1,𝒓_2)=0`$ and does not contribute). Using the normalization conditions $`𝑑𝒓_1𝑑𝒓_2\psi _\alpha ^{}(𝒓_1)\psi _\alpha (𝒓_1)\psi _\alpha (𝒓_2)\psi _\alpha ^{}(𝒓_2)=1`$ and $`𝑑𝒓_1𝑑𝒓_2\psi _\alpha ^{}(𝒓_1)\psi _\alpha (𝒓_1)\psi _\gamma (𝒓_2)\psi _\gamma ^{}(𝒓_2)=1`$, it follows that  $`𝑑𝒓_1𝑑𝒓_2\left[k_d(r)\left(1+\mathrm{\Pi }(𝒓_1,𝒓_2)\right)+b(𝒓_1,𝒓_2)\right]=0`$ and $`𝑑𝒓_1𝑑𝒓_2\left[1+k_d(r)\mathrm{\Pi }(𝒓_1,𝒓_2)+c(𝒓_1,𝒓_2)\right]=0`$, respectively. Thus the functions $`b`$ and $`c`$ do not contribute to order $`1/g^2`$ and

$$\sigma ^2(v_{\alpha \gamma })=2\kappa \left(\lambda \mathrm{\Delta }/g\right)^2,$$
(12)
where $`g=2\pi \nu DL^{d2}`$ is the dimensionless Thouless conductance. For a cube in $`d`$ dimensions $`\kappa =\frac{4}{\pi ^4}_{𝒏0}1/|𝒏|^4`$, where the summation is over $`𝒏=(n_1,\mathrm{},n_d)`$ with $`n_i`$ non-negative integers \[$`𝒏(0,\mathrm{},0)`$\].

We next consider the addition of $`m`$ electrons into the dot. Extending Koopmans’ limit to $`m`$ electrons, i.e., assuming the single-particle wave functions do not change with the addition of $`m`$ electrons, we can express the change in a single-particle level $`\alpha `$ as

$$\delta ϵ_\alpha ^{(m)}ϵ_\alpha ^{(𝒩+m)}ϵ_\alpha ^{(𝒩)}\underset{i=1}{\overset{m}{}}v_{\alpha 𝒩+i}^A.$$
(13)
For $`m<g`$, the variance of (13) is given by $`\sigma ^2(\delta ϵ_\alpha ^{(m)})=m\sigma ^2(v_{\alpha \gamma })+m(m1)(\overline{v_{\alpha \gamma }v_{\alpha \delta }}\overline{v_{\alpha \gamma }}\overline{v_{\alpha \delta }})`$ (since the variances and covariances are independent of the particular wave functions). In analogy to Eq. (7), the covariance of $`v_{\alpha \gamma }`$ and $`v_{\alpha \delta }`$ is given by

$`\overline{v_{\alpha \gamma }v_{\alpha \delta }}`$ $``$ $`\overline{v_{\alpha \gamma }}\overline{v_{\alpha \delta }}=`$ (14)
$`\lambda ^2`$ $`\mathrm{\Delta }^2`$ $`V^2{\displaystyle }d𝒓_1{\displaystyle }d𝒓_2\{[\psi _\alpha ^{}(𝒓_1)\psi _\alpha (𝒓_1)\psi _\alpha (𝒓_2)\psi _\alpha ^{}(𝒓_2)\psi _\gamma ^{}(𝒓_1)\psi _\gamma (𝒓_1)\psi _\delta (𝒓_2)\psi _\delta ^{}(𝒓_2){\displaystyle \frac{1}{V^4}}]`$ (16)
$`+[\psi _\alpha ^{}(𝒓_1)\psi _\alpha (𝒓_1)\psi _\delta (𝒓_2)\psi _\delta ^{}(𝒓_2)\psi _\gamma ^{}(𝒓_1)\psi _\gamma (𝒓_1)\psi _\alpha (𝒓_2)\psi _\alpha ^{}(𝒓_2){\displaystyle \frac{1}{V^4}}]\},`$
and using the known wave-function correlations

$`\overline{v_{\alpha \gamma }v_{\alpha \delta }}`$ $``$ $`\overline{v_{\alpha \gamma }}\overline{v_{\alpha \delta }}={\displaystyle \frac{\lambda ^2\mathrm{\Delta }^2}{V^2}}{\displaystyle }d𝒓_1{\displaystyle }d𝒓_2\{[1+k_d(r)(1+\mathrm{\Pi }(𝒓_1,𝒓_2))+\mathrm{\Pi }(𝒓_1,𝒓_2)+b(𝒓_1,𝒓_2)]`$ (17)
$`\times `$ $`[1+k_d(r)\mathrm{\Pi }(𝒓_1,𝒓_2)+c(𝒓_1,𝒓_2)]1+[1+k_d(r)\mathrm{\Pi }(𝒓_1,𝒓_2)+c(𝒓_1,𝒓_2)]^21\}.`$ (18)
Terms such as $`k_d^2\mathrm{\Pi }`$, $`k_d\mathrm{\Pi }^2`$ and $`k_dc`$ in Eq. (17) contribute $`(\mathrm{}/L)^2\lambda ^2\mathrm{\Delta }^2/g^3`$ and behave like $`\lambda ^2\mathrm{\Delta }^2/g^3`$ in the ballistic limit $`\mathrm{}L`$. Other terms contribute $`1/g^4`$. Thus to leading order in $`1/g^2`$

$$\overline{v_{\alpha \gamma }v_{\alpha \delta }}\overline{v_{\alpha \gamma }}\overline{v_{\alpha \delta }}0.$$
(19)
Using Eqs. (13), (12) and (19) we find

$$\sigma ^2(\delta ϵ_\alpha ^{(m)})=m\sigma ^2(\delta ϵ_\alpha ^{(1)})=2m\kappa (\lambda \mathrm{\Delta }/g)^2.$$
(20)
We conclude that the statistical fluctuations accumulated by sequentially adding electrons into the dot are added independently, namely, the FSD of an energy level behaves like $`\sqrt{m}`$, contrary to the parametric approach where it is linear in $`m`$ (see Eq. (1)).

A complete scrambling of the spectrum corresponds to a number of electrons $`m_c`$ for which $`\sigma (\delta ϵ_\alpha ^{(m_c)})\mathrm{\Delta }`$. Using (20) we find

$$m_c\left(\frac{g}{\lambda }\right)^2.$$
(21)
We note that since we have neglected wave-function scrambling due to off-diagonal interaction matrix elements, our estimate in Eq. (21) should only be regarded as an upper bound.

Next we consider the effect due to excess negative charge on the boundaries of the dot. The effective interaction due to such charge is  $`v(𝒓,𝒓^{})=V_\kappa (𝒓)+V_\kappa (𝒓^{})`$ where $`V_\kappa (𝒓)`$ describes the variation of the mean-field potential upon the addition of one electron. For a disk of radius $`R`$ in 2D we have $`V_\kappa (𝒓)=(e^2/2\kappa R)(R^2r^2)^{1/2}`$ where $`\kappa `$ is the inverse screening length. The related interaction matrix element is now $`v_{\alpha \gamma }=𝑑𝒓|\psi _\alpha (𝒓)|^2V_\kappa (𝒓)+𝑑𝒓|\psi _\gamma (𝒓)|^2V_\kappa (𝒓)`$. Using the wavefunction correlations to order $`1/g^2`$, its variance is

$`\sigma ^2(v_{\alpha \gamma })`$ $`=`$ $`{\displaystyle \frac{2}{V^2}}{\displaystyle 𝑑𝒓_1𝑑𝒓_2\left[\left(V^2|\psi _\alpha (𝒓_1)|^2|\psi _\alpha (𝒓_2)|^21\right)+\left(V^2|\psi _\alpha (𝒓_1)|^2|\psi _\beta (𝒓_2)|^21\right)\right]}`$ (24)
$`\times V_\kappa (𝒓_1)V_\kappa (𝒓_2)={\displaystyle \frac{4}{\beta V^2}}{\displaystyle }d𝒓_1{\displaystyle }d𝒓_2\{k_d(r)[1+(\beta +2)\mathrm{\Pi }(𝒓_1,𝒓_2)/\beta ]`$
$`+\mathrm{\Pi }(𝒓_1,𝒓_2)+b(𝒓_1,𝒓_2)+c(𝒓_1,𝒓_2)\}V_\kappa (𝒓_1)V_\kappa (𝒓_2).`$
The leading order contribution to (24) comes from $`\mathrm{\Pi }(𝒓_1,𝒓_2)`$. Using the expansion of $`\mathrm{\Pi }(𝒓_1,𝒓_2)`$ in terms of the diffusion modes $`\varphi _𝑸`$, the respective contribution is

$$\sigma ^2(v_{\alpha \gamma })=\frac{4}{\beta V^2}\frac{1}{\pi \nu }\underset{𝑸0}{}\frac{\left[𝑑𝒓\varphi _𝑸(𝒓)V_\kappa (𝒓)\right]^2}{D𝑸^2}=a\frac{\lambda ^2\mathrm{\Delta }^2}{\beta g}.$$
(25)
For a disk in 2D the coefficient $`a=(2\pi )^1_𝒒\left[𝑑𝒓\varphi _𝒒(𝒓)(1r^2)^{1/2}\right]^2/q^2`$, where $`\varphi _𝒒`$ are the diffusion modes in a disk of radius $`R=1`$ ($`𝑸=R^1𝒒`$). The short-range term $`k_d`$ in (24) contributes a term $`(\mathrm{}/L)^2(\mathrm{\Delta }^2/\beta g)`$ that for a diffusive dot is parametrically small compared with (25), but in the ballistic limit its contribution becomes comparable to (25). The surface charge contribution to the variance (25) dominates the bulk contribution (12). The bulk-surface cross-correlation term is given by $`(2\lambda \mathrm{\Delta }/V^2)𝑑𝒓_1𝑑𝒓_2\left[V^2|\psi _\alpha (𝒓_1)|^2|\psi _\alpha (𝒓_2)|^21\right]V_\kappa (𝒓_2)`$. The dominating contribution due to $`k_d`$ ($`(\mathrm{}/L)^2\mathrm{\Delta }^2/g`$) is suppressed in a diffusive dot but is comparable to (25) in the ballistic limit.

Similarly, the covariance of $`v_{\alpha \gamma }`$ and $`v_{\alpha \delta }`$ (due to surface charge) is evaluated to be

$`\overline{v_{\alpha \gamma }v_{\alpha \delta }}\overline{v_{\alpha \gamma }}\overline{v_{\alpha \delta }}`$ (26)
$`={\displaystyle \frac{1}{V^2}}{\displaystyle 𝑑𝒓_1𝑑𝒓_2\left[\left(V^2|\psi _\alpha (𝒓_1)|^2|\psi _\alpha (𝒓_2)|^21\right)+3\left(V^2|\psi _\alpha (𝒓_1)|^2|\psi _\beta (𝒓_2)|^21\right)\right]V_\kappa (𝒓_1)V_\kappa (𝒓_2)}`$ (27)
$`={\displaystyle \frac{2}{\beta V^2}}{\displaystyle 𝑑𝒓_1𝑑𝒓_2\left[k_d(r)\left(1+\frac{3\beta +2}{\beta }\mathrm{\Pi }(𝒓_1,𝒓_2)\right)+\mathrm{\Pi }(𝒓_1,𝒓_2)+b(𝒓_1,𝒓_2)+3c(𝒓_1,𝒓_2)\right]}`$ (28)
$`\times V_\kappa (𝒓_1)V_\kappa (𝒓_2).`$ (29)
The leading order contributions for both a diffusive and ballistic dots are similar to the ones contributing to the variance, and comparing (26) with (24), we find

$$\overline{v_{\alpha \gamma }v_{\alpha \delta }}\overline{v_{\alpha \gamma }}\overline{v_{\alpha \delta }}\frac{1}{2}\sigma ^2(v_{\alpha \gamma }).$$
(30)
We therefore obtain the following relation

$$\sigma ^2(\delta ϵ_\alpha ^{(m)})=\frac{1}{2}m(m+1)\sigma ^2(\delta ϵ_\alpha ^{(1)})=\frac{1}{2}m(m+1)a\frac{\mathrm{\Delta }^2}{\beta g}.$$
(31)
Thus, for $`m1`$, the FSD of a level is linear in the number of added electrons $`m`$, in agreement with the result (1) of the parametric approach. The scrambling parameter $`\delta \overline{x}_1`$ can be determined by comparing the variance (1) in the parametric approach with the corresponding variance (31) in the microscopic approach

$$\delta \overline{x}_1=\left(\frac{a}{2\beta g}\right)^{1/2}.$$
(32)
A complete scrambling of the spectrum corresponds to a parametric change of one avoided crossing, i.e., $`\delta \overline{x}_m1`$. Using (32), we find

$$m_c\frac{1}{\delta \overline{x}_1}(\beta g)^{1/2}.$$
(33)
When an interaction scale $`\lambda `$ is introduced (e.g., screening due to external gates), we find $`m_c(\beta g)^{1/2}/\lambda `$.

In a recent experiment  the standard deviation $`\sigma _p(m)\sigma (V_{𝒩+m}V_𝒩)`$, where $`V_𝒩`$ is the position in gate voltage for the $`𝒩`$-th peak, was measured as a function of the separation $`m`$ in peak number. In Koopmans’ limit $`V_{𝒩+m}V_𝒩ϵ_{𝒩+m}^{(𝒩+m)}ϵ_𝒩^{(𝒩)}=\mathrm{\Delta }_m^{(𝒩+m)}+\delta ϵ_𝒩^{(m)}`$, where $`\mathrm{\Delta }_m^{(𝒩+m)}ϵ_{𝒩+m}^{(𝒩+m)}ϵ_𝒩^{(𝒩+m)}`$ is the distance between levels separated by $`m`$ consecutive spacings in a dot with a fixed number of electrons ($`𝒩+m`$). The latter quantity is unrelated to scrambling and its fluctuation for $`m<g`$ is completely determined from RMT: $`\sigma (\mathrm{\Delta }_m^{(𝒩+m)})\mathrm{\Delta }\mathrm{ln}m`$. In principle one can extract $`\sigma (\delta ϵ_𝒩^{(m)})`$ from the measured $`\sigma _p(m)`$. In practice this is difficult to do with the available experimental statistics since for $`m<m_c`$, $`\sigma (\delta ϵ_𝒩^{(m)})`$ is smaller than $`\sigma (\mathrm{\Delta }_m^{(𝒩+m)})`$, and $`\sigma _p(m)`$ is dominated by the $`\mathrm{ln}m`$ term of RMT. Nevertheless, $`\sigma _p(m)`$, when plotted versus $`\mathrm{ln}m`$, appears to follow a straight-line behavior with a small concave deviation , in qualitative agreement with an additional contribution due to scrambling.

In conclusion, we have studied spectral scrambling when several electrons are added to a Coulomb-blockade quantum dot. In the absence of surface charge the fluctuation standard deviation of a level upon the addition of $`m`$ electrons behaves like $`\sqrt{m}`$, but when surface charge is present this standard deviation is linear in $`m`$. The latter result is the one obtained in the conventional parametric random matrix approach. Our analysis was carried out within the framework of the HF-Koopmans picture and was limited to energy fluctuations. A microscopic description of wave functions scrambling is an outstanding problem.

We acknowledge H. Baranger, C. M. Marcus and A. D. Mirlin for useful discussions. We also thank S. Altman for a Yale-Weizmann cooperation fellowship to support the visit of Y.G. to Yale, where part of this work was completed. This work was supported in part by the Department of Energy grant No. DE-FG-0291-ER-40608, by the U.S.-Israel Binational Science Foundation, by the DIP Foundation, by the ISF of the Israel Academy of Science and Humanities, and the Minerva Foundation.


format: text

text:
# Recovering physical parameters from galaxy spectra using MOPED

## 1 Introduction

Most of the information about the physical properties of galaxies comes from their electromagnetic spectrum. It is therefore of paramount importance to be able to extract as much physical information as possible from it. In principle, it is straightforward to determine physical parameters from an individual galaxy spectrum. The method consists of building synthetic stellar population models which cover a large enough range in the parameter space and then use a merit function (typically $`\chi ^2`$ if the pixels are not correlated) to evaluate which suite of parameters fits the observed spectrum best. There are two obvious limitations of the above method: first, the number of parameters that govern the spectrum of a galaxy may be very large and thus difficult to explore fully; secondly, in the case of ongoing large redshifts surveys which will provide us with about a million galaxy spectra, it will be computationally very expensive (and possibly intractable for redshift surveys like the 2dF and SDSS) to apply a brute-force $`\chi ^2`$ analysis to each individual spectrum which itself may contain of the order of $`10^3`$ data points.

A less obvious route to tackle the high computational requirement is to compress the original data set, giving more weight to those pixels in the spectrum that carry most information about a given parameter. In this paper we show how this can be done in an optimal way. It is worth remembering that data compression is commonly applied to galaxy spectra, either by the instrument, through the use of photometric filters, or in the interpretation, by concentrating on specific spectral features and ignoring others. Not surprisingly, this empirical data compression is not optimal since it is ad hoc. For example, the photometric $`B`$ filter alone is not optimal to recover the age of a galaxy. On the other hand, more sophisticated and non-empirical methods have been proposed for extracting information from galaxy spectra, some of them as old as the Johnson’s filter system. Many of these are based on Principal Component Analysis (PCA) or wavelet decomposition (???????????). PCA projects galaxy spectra onto a small number of orthogonal components. The weighting of each component corresponds to its relative importance in the spectrum. However while these components appear to correlate reasonably well with physical properties of galaxies, their interpretation is difficult since they do not have known, specific physical properties – they can be amalgams of different properties. To interpret these components, we have to return to model spectra and compare them with the components (?). This is a disadvantage of PCA since one important goal of the analysis is to study the evolution of the physical properties which dramatically affect galaxy spectra, such as the age, metallicity, star formation history or dust content. It is important to recognise that PCA can play an important role if there is no underlying model for how the data should behave. If such a model exists, then one can do better by using projections of the data which are designed to give the parameters of the model as accurately as possible.

An optimal parameter–extraction method, which we term MOPED111The MOPED algorithm has a patent pending (Multiple Optimised Parameter Estimation and Data compression) was developed in ?. The purpose of this paper is to apply the MOPED method to a specific set of observations, to estimate their physical parameters and thus demonstrate the usefulness of the approach. We also want to demonstrate with a specific example the massive speed up factor that the method provides, making computationally-intensive problems into much more accessible ones. The outline of the paper is at follows: in section II we briefly describe the method used; in section III we discuss the problems involved in extracting physical parameters from galaxy spectra: how to choose the best parametrisation for the star formation history; in section IV we derive physical parameters for the Kennicutt atlas and present the best fitting models in an appendix.

## 2 The method

The MOPED algorithm was presented in detail in ?. Here we briefly recall the main ingredients and steps of the method.

The main idea of the method is that, in practice, some of the data may tell us very little about the parameters we are trying to estimate, either through being very noisy, or through having no sensitivity to the parameters. So in principle, we may be able to throw away some data without losing much information about the parameters. It is obvious that simply throwing away some of the data is not in general optimal; it will usually lose information. On the other hand, by constructing linear combinations of the data we might do better and then we can throw away the linear combinations which tell us least. In fact one can do much better than this. Providing the noise has certain properties, one can reduce the size of the dataset down to a handful of linear combinations – one for each parameter – which contain as much information as the entire dataset. It is by no means obvious that this can be done.

Given a set of data x (in our case the spectrum of a galaxy) which includes a signal part $`𝝁`$ and noise $`𝐧`$, i.e. $`𝐱=𝝁+𝐧`$, the idea then is to find weighting vectors $`𝐛_m`$ such that $`y_m𝐛_m^t𝐱`$ contain as much information as possible about the parameters (star formation rates, metallicity etc.). These numbers $`y_m`$ are then used as the data set in a likelihood analysis.

In ? an optimal and lossless method was found to calculate $`𝐛_m`$ for multiple parameters (as is the case with galaxy spectra). The definition of lossless here is that the Fisher matrix at the maximum likelihood point (see ?) is the same whether we use the full dataset or the compressed version. The Fisher matrix gives a good estimate of the errors on the parameters, provided the likelihood surface is well described by a multivariate Gaussian near the peak. We find that the method is lossless provided that the noise is independent of the parameters. This is not exactly true for galaxy spectra, owing to the presence of a shot noise component from the source photons. However, the increase in parameter errors is very small in this case (see ?). The weights required are

$$𝐛_1=\frac{𝐂^1𝝁_{,1}}{\sqrt{𝝁_{,1}^t𝐂^1𝝁_{,1}}}$$
(1)
and

$$𝐛_m=\frac{𝐂^1𝝁_{,m}\underset{q=1}{\overset{m1}{}}(𝝁_{,m}^t𝐛_q)𝐛_q}{\sqrt{𝝁_{,m}𝐂^1𝝁_{,m}_{q=1}^{m1}(𝝁_{,m}^t𝐛_q)^2}}.$$
(2)
where a comma denotes the partial derivative with respect to the parameter $`m`$ and $`C`$ is the covariance matrix with components $`C_{ij}=n_in_j`$. $`m`$ runs from 1 to the number of parameters $`M`$, and $`i`$ and $`j`$ from 1 to the size of the dataset (the number of flux measurements in a spectrum). To compute the weight vectors requires an initial guess of the parameters. We term this the fiducial model.

The dataset $`\{y_m\}`$ is orthonormal: i.e. the $`y_m`$ are uncorrelated, and of unit variance. The new likelihood is easy to compute (the $`y_m`$ have mean $`y_m=𝐛_m^t𝝁`$), namely:

$$\mathrm{ln}(\theta _\alpha )=\mathrm{constant}\underset{m=1}{\overset{M}{}}\frac{(y_my_m)^2}{2}.$$
(3)
Further details are given in ?.

It is important to note that if the covariance matrix is known for a large dataset (e.g. a large galaxy redshift survey) or it does not change significantly from spectrum to spectrum, then the $`y_m`$ need be computed only once for the whole dataset, thus with massive speed up factors in computing the likelihood as will be shown in sections 3 and 4. Note that the $`y_m`$ are only orthonormal if the fiducial model coincides with the correct one. In practice one finds that the recovered parameters are almost completely independent of the choice of fiducial model, but one can iterate if desired to improve the solution.

### 2.1 Estimating errors in the recovered parameters

If we have $`M`$ parameters, we can get the marginal error on each of the parameters from the Fisher matrix

$$F_{ij}\frac{^2\mathrm{ln}}{\theta _i\theta _j}$$
(4)
where $`i`$ and $`j`$ run from 1 to $`M`$ (Note the size of the dataset appears nowhere here - it is only relevant in computing $`\mathrm{ln}`$). The averaging is over many realisations.

The conditional error on $`\theta _i`$ (fixing all other parameters) is just $`(F_{ii})^{1/2}`$ (no summation). This is not too relevant as we wish to estimate all the parameters simultaneously. More relevant then is the marginal error, which is

$$\sigma _i=\sqrt{(F^1)_{ii}}.$$
(5)
So the procedure to follow is to estimate $`F_{ij}`$ from the likelihood surface near the peak, invert the $`M\times M`$ Fisher matrix and use the diagonal components of the inverse matrix to assign marginal errors. All this assumes that the likelihood surface is well-approximated by a multivariate Gaussian, not just at the peak, but also far enough down the likelihood hill (so $`\mathrm{ln}`$ drops by around unity).

#### Fisher matrix estimation

Assuming that the maximisation finds the maximum precisely, at position $`𝜽_0`$ (an $`M`$dimensional vector in parameter space), with value $`\mathrm{ln}_0`$, then a Taylor expansion of $`\mathrm{ln}`$ around the maximum gives

$$\mathrm{ln}(𝜽_0+\mathrm{\Delta }𝜽)\mathrm{ln}_0+\frac{1}{2}\frac{^2\mathrm{ln}}{\theta _i\theta _j}\mathrm{\Delta }\theta _i\mathrm{\Delta }\theta _j$$
(6)
where the summation convention is assumed. It will not apply for the rest of this section.

We estimate the diagonal components of the second derivatives by keeping all components of $`𝜽`$ constant apart from a single $`\theta _i`$.

$`F_{ii}`$ $``$
$`{\displaystyle \frac{1}{\mathrm{\Delta }\theta _i^2}}\left[\mathrm{ln}(𝜽_0+\mathrm{\Delta }\theta _i𝐞_i)+\mathrm{ln}(𝜽_0\mathrm{\Delta }\theta _i𝐞_i)2\mathrm{ln}_0\right].`$
Similarly, the off-diagonal terms are estimated from

$`F_{ij}`$ $``$ $`{\displaystyle \frac{1}{2\mathrm{\Delta }\theta _i\mathrm{\Delta }\theta _j}}[`$
$`\mathrm{ln}(𝜽_0+\mathrm{\Delta }\theta _i𝐞_i+\mathrm{\Delta }\theta _j𝐞_j)+`$
$`\mathrm{ln}(𝜽_0\mathrm{\Delta }\theta _i𝐞_i\mathrm{\Delta }\theta _j𝐞_j)`$
$`\mathrm{ln}(𝜽_0\mathrm{\Delta }\theta _i𝐞_i+\mathrm{\Delta }\theta _j𝐞_j)`$
$`\mathrm{ln}(𝜽_0+\mathrm{\Delta }\theta _i𝐞_i\mathrm{\Delta }\theta _j𝐞_j)].`$
The above procedure is not computationally expensive - for $`M=17`$ we need $`17\times 16/2`$ off-diagonal terms, and $`17`$ diagonal terms. This requires about 600 likelihood evaluations. The marginal errors are estimated from $`\sigma _i1/\sqrt{(F^1)_{ii}}`$.

## 3 Determining physical parameters from galaxy spectra

### 3.1 The problem

Our aim is to determine the star formation history, metallicity and dust evolution of the stellar population of a galaxy from its spectrum. We wish to do this for galaxy spectra that typically contain thousands of data points. Additionally, current large surveys contain on the order of $`10^6`$ galaxies. It is computationally very expensive to compute a “brute-force” likelihood analysis with all the data points. For example, for a spectrum with $`10^3`$ points for which we want to determine 17 parameters, it takes about 1.5 hours of CPU time on an Athlon 1Ghz machine to find the minimum of the likelihood surface. With $`10^6`$ spectra, it would take about 170 years, and if the theoretical models change, it would take as long again. Clearly, the brute force approach is too slow by a factor of several hundred. As we show below, we can estimate 17 parameters of a galaxy with MOPED in about 20 seconds.

### 3.2 Choosing the optimal parametrisation

Which parametrisation should we choose in order to determine the star formation history (SFH) of a galaxy? The common procedure in the literature is to assume that the SFH proceeded as a decaying exponential law with one parameter, while more sophisticated modifications allow for the presence of a burst (or two) at a given time. Ideally one would prefer to avoid any assumptions about the form of the SFH. In fact, star formation in galaxies takes place in giant molecular clouds which are short lived (about $`10^7`$ years) due to the explosion of SN. It is therefore clear that we could parametrise the SFH in a rather model–independent way by dividing time into bins of width smaller than $`10^7`$ years and using the height of each bin as the parameter. Of course, for galaxies as old as $`10^{10}`$ years this approach becomes useless, since one can hardly deal with over a thousand parameters. We therefore choose a different approach. We divide time into a fixed number of wider bins, and choose their positions and widths according to the following scheme. We consider bursts of star formation at the beginning and end of each time bin (at a fixed metallicity), and require that the maximum fractional difference in the final spectrum is the same for each bin. The maximum difference depends on how many bins we wish to consider. To illustrate: for a galaxy whose age is 14 Gyr with a maximum error in the height of each bin of a factor of 2 (100% error), the optimal bin boundaries will be at 0.03, 0.23, 0.42, 0.66, 0.92, 2.56, 6.33 and 14 Gyr (considering fixed solar metallicity). Further exploration with different metallicities showed us that the bin boundaries differ by a small amount. In fact the bin boundaries are very close to equally spaced bins in logarithmic space – which for the above case are: 0.02, 0.05, 0.13, 0.33, 0.84, 2.15, 5.49 and 14.0 Gyr. For convenience, we choose equally-spaced logarithmic bin boundaries in what follows. Therefore, we are faced now with 8 parameters to determine the SFH. Note that we have not made any specific assumption about the actual shape of the SFH. Within this framework metallicity is extremely easy to parametrise since to each bin we simply assign an extra parameter which is the metallicity of the stars formed in that bin.

Another parameter to consider is the dust content of the galaxy. The importance of dust is being recognised (e.g. ?), but modelling is difficult and currently much less sophisticated than the stellar population modelling. Since the purpose of this paper is only to illustrate the usefulness of the method, we will not develop very sophisticated models for the process of dust emission and absorption. We use the ? parametrisation as sufficient to describe the major effect of dust absorption on the integrated light of galaxies. The Calzetti model depends only on one parameter: the amount of dust in the galaxy. We will therefore use this single parameter to describe the global effect of dust. Note that as more realistic dust parametrisations are proposed (e.g. ??), the search parameter space increases and data compression methods become increasingly necessary.

Emission lines are not naturally predicted in the spectrophotometric models. They could be included at the expense of a larger parameter space to explore. An alternative approach, which we follow here, is not to include the emission lines in the analysis at all. Thus we will not expect to fit the emission lines, but only the continuum and absorption lines. Similarly, we could include one (or more) parameters characterising the velocity dispersion of the galaxy. With the coarse spectral resolution of the models, we have chosen not to do this here.

### 3.3 Recovery of parameters of simulated spectra

In what follows, and for convenience, we will use the set of synthetic stellar population models developed in ?. We emphasize that the choice of models is not crucial at any point in our argument. Furthermore, any suit of models can be used in principle. It is worth noting that given the low spectral resolution of the models (20 Å) we are rather limited in our capacity to extract information from the spectrum. Also, and since we are aiming (in this particular paper) to extract parameters from the ? atlas, we will limit the wavelength coverage between 3800 and 7000 Å. Note that this is very restrictive since with a larger spectral coverage and, more importantly, better spectral resolution of the models, we would be able to extract a larger number of parameters with a smaller error (see section 5).

We are now in a position to test our method. We can do this by building synthetic models with known star formation histories, metallicities and dust and then try the parameter recovery. Fig. 1 shows the model spectrum (solid black line) which has been constructed out of five bins with values: 0.4, 0.1, 0.0, 0.0 and 0.5 with fixed solar metallicity (no dust) in each bin and to which artificial Poisson noise has been added ($`S/N=10`$). The dashed line shows the best recovered model with bin values corresponding to: 0.41, 0.10, 0.0, 0.0 and 0.49. The agreement is remarkable. Fig. 2 shows another example for constant metallicity, again the efficiency of the method is excellent.

Encouraged by the above results we now consider a more general case where we allow the metallicity of each bin to take an arbitrary value. We also include dust as a free parameter. The first synthetic test case has got artificial Poisson noise added with $`S/N=10`$. In this case the recovery of all 17 parameters is not optimal. The reason for this is that with such a low $`S/N`$ most of the information that is contained in the absorption features has disappeared; we are left only with a continuum. Given the fact of our limited wavelength range in this example (3800 to 7000 Å), it is almost impossible to break the “famous” age–metallicity degeneracy with only the continuum shape. This is in stark contrast with the first test case where the metallicity and dust where kept fixed, then it was possible with only the continuum to derive the heights in the bins. Therefore, we now increase the $`S/N`$ until we are able to recover the 17 parameters. This happens for $`S/N=20`$ and is illustrated in Fig. 3. It can be seen that the parameter recovery is very good.

Given the success of the above tests we decided to apply the method to a real sample, but before doing it we turn our attention to the physical interpretation of the eigenvectors $`𝐛`$.

### 3.4 Eigenvectors or where is the information?

We now concentrate our attention on the physical meaning of the eigenvectors computed from eq. (1) and (2). As shown in ? the choice of the fiducial model is not important. One can always iterate, but this appears to be quite unnecessary. Figures 4 and 5 show the SFH and metallicity eigenvectors respectively as a function of wavelength in the range 3600 to 7000 Å for a model in which we have parametrised the SFH by using 8 bins. Figure 4 concerns the SFH eigenvectors for bins corresponding to : 0.02, 0.05, 0.13, 0.33, 0.84, 2.15, 5.49 and 14.0 Gyr (from top to bottom), while the 8 panels on Figure 5 are the corresponding metallicity eigenvectors. Not surprisingly, for young ages, the age eigenvector likes to weight most those pixels bluewards of 4500 Å, while those redwards of 4500 Å get very little weight (note that the eigenvector can always be multiplied by $`1`$). As the population ages most of the information about the SFH becomes more distributed among wavelength, thus making it very difficult to design a narrow band filter which would capture most of the information. A similar situation occurs for the metallicity bins (see Figure 5, for young ages most of the information is concentrated on the bins bluewards of 5000 Å while for older ages the information becomes more equally distributed among different absorption features. Note that for recovering metallicity information the absorption features are more relevant than the continuum.

## 4 Physical parameters for the spectra in the Kennicutt atlas

In this section we apply the method to the Kennicutt atlas (?). Although the Kennicutt sample is certainly small, about 50 galaxies, the S/N of the spectra is very high. Furthermore, the largest redshift survey available at the moment, the 2dF, suffers from serious complications in calibrating photometrically the continuum of the spectra (?), thus, for the time being, we concentrate in illustrating the method on the small but significant atlas of observed spectra provided by ?. We have re-binned the data to the same spectral resolution as the models, with the new flux calculated as the mean flux unit wavelength in the bin. It is worth keeping in mind that the spectrophotometric calibration of the ? atlas has an error of 10% and that our models do not include emission lines in the spectra. It should be stressed that despite the excellent fits obtained to the ? atlas, one should not over–interpret the results since the 10% uncertainty in the spectrophotometric calibration affects the parameter determination. As stated above in the paper the physical parameters we are aiming at deriving are: the amount of star formation in the time bins chosen – here we chose 8, the metallicity in each bin and the global dust content, using the ? formula for the dust extinction. This gives a total of 17 parameters for each galaxy. This is the most effective way to parametrise the SF since it does not depend on any previous knowledge about its shape, i.e. we are not assuming that it is a declining exponential or similar. Note that if we use the full dataset instead of the compressed data, we get the same results (although it takes much longer (see section 5)).

Results are presented in the next 12 figures, where the best fit to the spectrum is drawn in the left panel and the corresponding SFH is plotted with horizontal error bars denoting the width of the bins and vertical error bars for the marginal uncertainty in the height of the bin. The derived metallicity for each bin is also labeled in the right panel as well as the dust with corresponding errors.

The first result is somewhat surprising. We find that all galaxies classified as E/S0 (NGC3245, NGC3379, NGC3516, NGC4472, NGC4648, NGC4889, NGC5866, NGC6052) show a reasonable level of recent star formation activity although in all cases more than 50% of the stellar mass is older than 5 Gyr. NGC3245 (S0) shows recent star formation below 7% of the total mass as does NGC6052 (S0). NGC3379, NGC3516, NGC4472, NGC4648, NGC4889 and NGC5866 all show a significant ( $`>`$ 30% of the total mass) amount of star formation in the last 3 Gyr (Note NGC3516 has a Seyfert nucleus). Note also that the metallicities derived are very reasonable since the final mass–weighted metallicity for the above galaxies is about 1.5 times the solar value222due to the fact that the ? sample is nearby, the slit of the spectrograph can only sample the inner regions of the galaxy, where the metallicity is higher.. Note that these findings are in excellent agreement with the recent findings about the nature of spheroids by e.g. ????. Also significant is that the method is able to recognize current star formation, despite the fact of not having any emission lines included in the model. For example, for galaxies with slight Hα (like NGC1357, NGC3147, NGC3227, NGC3921 and NGC4750) the method always shows a significant level of recent star formation in the most recent bins, which agrees with the star formation rate that would be estimated using the Hα line itself (NB NGC3227 is complicated by the presence of a Seyfert component). As expected, galaxies classified as Sa to Sc show more recent bursts of star formation, but note that star formation activity does not correlate with the Hubble type. In other words, star formation history in galaxies does not proceed as a single exponential decaying law that changes according to the Hubble type in a monotonic manner but it is much more like a sequence of burst events. We now turn our attention to the metallicity evolution of the stars. The general trend is that high-metallicity stars are formed at an early epoch, the typical values being over the solar value. On the other hand the youngest bursts tend to have slightly sub-solar metallicities (about half the solar value). Note that this trend is in good agreement with predictions from infall models (e.g. ?). These models are motivated by the fact that if all the gas was available for consumption into stars from the very beginning (closed–box model), then a large number of low metallicity stars should be observed now. What is observed is that the metallicities of stars in galaxies (both disks and ellipticals) are gaussianly distributed. An obvious solution to this problem is not to allow all the gas to be acreted into the galaxy at the beginning but at a slower pace. As a result of this the late forming generations of stars will have a reservoir of fresh low metal gas which will decrease the metallicity that would be predicted by a closed–box model, which increases monotonically with time. Inspection of the figures in the appendix shows that this is not the case, but that late forming generations have a lower metallicity. Note also that where the star formation is effectively zero, the metallicity is not recovered accurately, as one might expect. We note also that dust values seem to be very reasonable.

## 5 Discussion and conclusions

We have presented a fast and efficient algorithm (MOPED) to recover physical parameters from galaxy spectra. The algorithm is based on a data representation technique which discards most of the data set that does not contain any information with respect to a given parameter. By doing so we are able to speed up the parameter determination by a large factor: typically of the order of $`N/M=`$ number of data points/number of parameters. We note here that if the noise matrix is non-diagonal, the speed-up is much larger $`N^3/M`$ (see Gupta & Heavens, in preparation). We have also applied the method to the ? atlas and derive star formation and metallicity histories for the NGC galaxies in the survey. We have found that most spheroids, albeit having at least 50% stars formed at ages older than 5 Gyr, have had recent star formation episodes. The quality of the spectral models does not allow us yet to determine if those are the consequence of recent mergers.

The method presented in section 2.1 is good at estimating the local errors around the maximum. Unfortunately, the thing that usually happens is that the error estimation is mostly dominated by local maxima that are in distant regions of the parameter space, this is our case. In order to explore this we record all the maxima that the method finds starting from different initial random guesses (in this case 5000) 333The space parameter is so large that one cannot compute the whole likelihood surface, instead we used a conjugate gradient method to find the maximum of the likelihood. In order to test that we had found the absolute maximum we repeated the process from new random initial guesses. We found that 400 of these guesses suffice to find the global maximum..We then explore if some of these solutions are also allowed in the $`\chi ^2`$ sense. We have done this for all the spectra in the ? atlas. Our findings are as follows: for systems with only absorption features the number of degenerate models is almost negligible, thus the best fitting model (which is plotted in the appendix) corresponds quite closely to an absolute maximum. On the other hand, for those spectra with strong emission lines degeneracies are more severe, to such extent that for NGC1569, NGC3310, NGC3690, NGC4194, NGC4449, NGC4485, NGC4631, NGC4670, NGC4775, NGC5996, NGC6052 and NGC6240 other acceptable fits that differ significantly from the best fit presented in the appendix exist. For example, an acceptable fit is also one where most of the star formation takes place in the youngest bins and the metallicity is about solar. This of course would be cured if we were able to have more detailed absorption features in the spectrum and especially more resolution, since this would constrain some of the models that now appear as plausible. Indeed, a larger wavelength coverage would also help. In any case, it is worth reminding that this degeneracy is generic of the present problem we are dealing with, and that using MOPED only facilitates the search in a large parameter space for such degeneracies and to understand the physics behind the problem.

The speed up factors obtained in the present analysis are quite remarkable. On a 1Ghz Athlon machine it takes much less than a second for an initial starting guess to find a (possibly local) maximum of the likelihood. For 400 initial guesses it takes about 20 seconds of CPU time to find the global maximum. If one used the whole data set the same calculation would take about an hour and a half. Extrapolating these numbers one can see that for analysing a big dataset like the 2dF or SDSS, the time required using all of the data would be about 170 years, longer if higher-resolution models are used. This is obviously too long to be achieved, especially if one takes into account that the theoretical model may change, necessitating a full repeat of the analysis. With MOPED, it would take a few months on a single workstation, or much less by exploiting the trivially parallel nature of the problem. Finally, the rapid algorithm here has allowed us to explore more fully the likelihood surface, and we find that for galaxies with significant very recent star formation there is often more than one acceptable parameter fit. There is clearly a need for theoretical stellar atmospheric models with spectral resolution better than 1Å if we are to extract all the rich information that large redshift surveys contain in their galaxy spectra.

Acknowledgments

We are grateful to Karl Glazebrook, Ofer Lahav and Marc Kamionkowski for useful discussions.

## Appendix I

The following 12 figures show the best fitting model (thick solid line) for each spectrum (thin solid line) in the ? atlas (left panel) and the plot of the star formation history for the 8 age bins (right panel). The SFH is plotted with horizontal error bars denoting the width of the bins and vertical error bars for the marginal uncertainty in the height of the bins. Also shown for each model, is the value of the metallicity derived for each bin in solar units (left to right correspond to youngest to oldest bin) and the overall amount of extinction according to the ? model. For those bins with zero star formation the metallicity is irrelevant and set to zero (and zero error). Note that there are five galaxies with Seyfert nuclei: NGC3227, 3516, 5548, 6764, and 7469.


text:
# Scattering theory for p-forms on hyperbolic real space

## 1. Introduction

Scattering theory for the wave equation has been developed by Lax and Phillips essentially for functions, at first in euclidean setting (see  and the references therein), and more recently, after the papers of Faddeev and Pavlov (see ) and the construction of a Fourier analysis for functions on symmetric spaces due to Helgason (see ), in non-euclidean setting. Particular attention has been devoted in , , ,  to the wave equation for functions on the hyperbolic space $`H^{n+1}`$ and on quotients of $`H^{n+1}`$ by groups of isometries.

The purpose of this paper is to extend the Lax-Phillips theory to the case of the wave equation for $`p`$-forms, or more generally for sections of vector bundles, on complete, non compact Riemannian manifolds. The basic motivation is connected, on one hand, with the spectral properties of the Laplace-Beltrami operator on differential forms, and hence in some way with the “geometry” of the manifold, and on the other side with the validity of Huygens’ principle for the wave equation, and therefore with the “analysis” that can be developed on the manifold.

A general approach to these problems on a general manifold turns out to be extremely difficult to deal with, for the lack of information about the behaviour of Laplace-Beltrami operator on a generic complete non compact Riemannian manifold, and of analytic tools mimicking the Radon and Fourier transforms.

These facts motivate our choice to investigate the problem at first on the hyperbolic space $`H^{n+1}`$. In this case, the spectrum of the Laplace-Beltrami operator on $`p`$-forms on $`H^{n+1}`$ has been explicitely computed by Donnelly (see ). Moreover in recent years (see ) it has been developed a Radon transform for $`p`$-forms on $`H^{n+1}`$, providing an essential tool to the study of the wave equation.

In section 2, we recall some essential facts about $`H^{n+1}`$ and its spectral properties. In section 3, we prove that it is not possible to develop a scattering theory in the Lax-Phillips sense for the wave equation for all $`p`$-forms on $`H^{n+1}`$. In section 4, we show that, as a consequence, Huygens’ principle for the wave equation in this context does not hold. In section 5 we show that since Huygens’ principle does hold for the wave equation for coclosed $`p`$-forms when $`n=2p`$, a scattering theory in the Lax-Phillips sense can be developed if we impose some restrictions on the class of forms considered. This result implies some conclusions about the spectrum of the Laplace-Beltrami operator on coclosed $`p`$-forms on $`H^{n+1}`$, when $`n=2p`$.

This paper summarizes some results of my PhD Thesis, compiled under the direction of Prof. Edoardo Vesentini, to whom I am deeply grateful for his constant and careful support.

## 2. Preliminaries

We will denote by $`H^{n+1}`$ the real hyperbolic space of dimension $`n+1`$, either in the half-space model

$$H^{n+1}=\{(x_1,\mathrm{},x_{n+1})R^{n+1}x_{n+1}>0\}$$

endowed with the Riemannian metric

$$ds^2=\frac{dx_1^2+dx_2^2+\mathrm{}+dx_{n+1}^2}{x_{n+1}^2},$$

or in the Poincaré model

$$H^{n+1}=\{(x_1,\mathrm{},x_{n+1})^{n+1}|x|^2<1\}$$

with the Riemannian metric

$$d\sigma ^2(x)=\left(\frac{2}{1|x|^2}\right)^2(dx_1^2+\mathrm{}+dx_{n+1}^2),$$

where $`|x|^2=x_1^2+x_2^2+\mathrm{}+x_{n+1}^2`$.

We will denote by $`\mathrm{\Lambda }^p(H^{n+1})`$ the space of smooth $`p`$-forms on $`H^{n+1}`$, and by $`\mathrm{\Lambda }_c^p(H^{n+1})`$ the space of smooth compactly supported $`p`$-forms on $`H^{n+1}`$. As usual, $`L_p^2(H^{n+1})`$ will be the Hilbert space of square integrable $`p`$-forms on $`H^{n+1}`$, that is the completion of $`\mathrm{\Lambda }_c^p(H^{n+1})`$ with respect to the $`L^2`$ norm induced by the Riemannian metric:

$$\omega _{L^2}^2=_{H^{n+1}}\omega \omega .$$

Since $`H^{n+1}`$ is a complete Riemannian manifold, the Laplace-Beltrami operator $`\mathrm{\Delta }_p:\mathrm{\Lambda }_c^p(H^{n+1})\mathrm{\Lambda }_c^p(H^{n+1})`$ is an essentially selfadjoint operator in $`L_p^2(H^{n+1})`$. We will denote its closure again by the same symbol $`\mathrm{\Delta }_p`$.

H. Donnelly in  has explicitely computed the spectrum of the Laplace-Beltrami operator on $`p`$-forms on $`H^{n+1}`$. He has proved that $`\mathrm{\Delta }_p`$ has point spectrum if and only if $`n+1=2p`$. In this case, the only eigenvalue of $`\mathrm{\Delta }_p`$ is zero, occurring with infinite multiplicity, and the continuous part of the operator is unitarily equivalent to the operator

$$L^2(^+,dr,MM)L^2(^+,dr,MM)$$

$$f_1(x)f_2(x)\left(\frac{1}{4}+x^2\right)f_1(x)\left(\frac{1}{4}+x^2\right)f_2(x),$$

where $`M`$ is a Hilbert space of countably infinite dimension and $`dr`$ is the Lebesgue measure; hence, if $`n+1=2p`$, the absolutely continuous spectrum of $`\mathrm{\Delta }_p`$ is $`\sigma _{ac}(\mathrm{\Delta }_p)=[\frac{1}{4},+\mathrm{})`$, with constant multiplicity.

If $`p\frac{n+1}{2}`$, the point spectrum is empty, and $`\mathrm{\Delta }_p`$ is unitarily equivalent to the operator

$$L^2(^+,dr,MM)L^2(^+,dr,MM)$$

$$f_1(x)f_2(x)\left((\frac{n2p}{2})^2+x^2\right)f_1(x)\left((\frac{n2p+2}{2})^2+x^2\right)f_2(x).$$

Hence $`\sigma (\mathrm{\Delta }_p)=\sigma _{ac}(\mathrm{\Delta }_p)=[\text{min}\{(\frac{n2p}{2})^2,(\frac{n2p+2}{2})^2\},+\mathrm{})`$, and the multiplicity of the spectrum is not constant.

In the following we will suppose $`p\frac{n+1}{2}`$.

## 3. General case

Consider the Cauchy problem for the wave equation for $`p`$-forms on $`H^{n+1}`$:

(3.1)
$$\{\begin{array}{c}\mathrm{\Delta }_p\omega +_{tt}^2\omega =0\hfill \\ \omega (0)=\omega _1\hfill \\ _t\omega (0)=\omega _2,\hfill \end{array}$$

with $`\omega _1,\omega _2`$ smooth and compactly supported. For every $`t0`$

(3.2)
$$<\mathrm{\Delta }_p\omega (t),\omega (t)>+<_t\omega (t),_t\omega (t)>=<\mathrm{\Delta }_p\omega _1,\omega _1>+<\omega _2,\omega _2>,$$

because

$$_t(<\mathrm{\Delta }_p\omega (t),\omega (t)>+<_t\omega (t),_t\omega (t)>)=$$

$$=<\mathrm{\Delta }_p(_t\omega ),\omega >+<\mathrm{\Delta }_p\omega ,_t\omega >+$$

$$+<_{tt}^2\omega ,_t\omega >+<_t\omega ,_{tt}^2\omega >=0.$$

Hence, the energy preserved by the system is given by

$$(\omega _1,\omega _2)_H^2=<\mathrm{\Delta }_p\omega _1,\omega _1>+<\omega _2,\omega _2>0,$$

which defines indeed a norm on $`\mathrm{\Lambda }_c^p(H^{n+1})`$, because $`(\omega _1,\omega _2)_H=0`$ implies:

$$\{\begin{array}{c}d\omega _1=0\hfill \\ \delta \omega _1=0\hfill \\ \omega _2=0,\hfill \end{array}$$

that is

$$\{\begin{array}{c}\mathrm{\Delta }_p\omega _1=0\hfill \\ \omega _2=0.\hfill \end{array}$$

Since $`p\frac{n+1}{2}`$ and therefore $`0`$ is not in the point spectrum, $`\mathrm{\Delta }_p\omega _1=0`$ implies $`\omega _1=0`$, showing in conclusion that

$$(\omega _1,\omega _2)_H=0(\omega _1,\omega _2)=(0,0).$$

It is therefore meaningful to take the completion $`H`$ of the space of smooth compactly supported initial data with respect to this latter norm, which is induced by the scalar product

(3.3)
$$<(\omega _1,\omega _2),(\stackrel{~}{\omega }_1,\stackrel{~}{\omega }_2)>=<\mathrm{\Delta }_p\omega _1,\stackrel{~}{\omega }_1>+<\omega _2,\stackrel{~}{\omega }_2>.$$

The operator $`A`$ defined on

$$D(A)=\mathrm{\Lambda }_c^p(H^{n+1})\mathrm{\Lambda }_c^p(H^{n+1})$$

by

(3.4)
$$A=\left(\begin{array}{cc}0& I\\ \mathrm{\Delta }_p& 0\end{array}\right)$$

is essentially skewselfadjoint in $`H`$. Hence its closure $`\overline{A}`$ is the generator of a one-parameter group of unitary operators $`U(t)`$ on $`H`$ which maps every initial data into the value at time $`t`$ of the corresponding solution of the wave equation for $`p`$-forms on $`H^{n+1}`$.

The operator $`i\overline{A}`$ is essentially selfadjoint and can be shown to be unitarily equivalent to the operator

$$\left(\begin{array}{cc}\sqrt{\mathrm{\Delta }_p}& 0\\ 0& \sqrt{\mathrm{\Delta }_p}\end{array}\right)$$

on $`L_p^2(H^{n+1})L_p^2(H^{n+1})`$ (see ).

Now, thanks to Donnelly’s result, $`\sqrt{\mathrm{\Delta }_p}`$ is unitarily equivalent to the direct sum of multiplication operators

$$L^2(^+,dr,MM)L^2(^+,dr,MM)$$

$$f_1(x)f_2(x)\left(\sqrt{(\frac{n2p}{2})^2+x^2}\right)f_1(x)\left(\sqrt{(\frac{n2p+2}{2})^2+x^2}\right)f_2(x).$$

As a consequence, if $`n+12p`$ the spectrum of $`i\overline{A}`$ is purely absolutely continuous, equal to

$$(\mathrm{},\mathrm{max}\{\left|\frac{n2p}{2}\right|,\left|\frac{n2p+2}{2}\right|\}]$$

$$$$

$$[\mathrm{min}\{\left|\frac{n2p}{2}\right|,\left|\frac{n2p+2}{2}\right|\},+\mathrm{})$$

and does not have constant multiplicity.

Now, if it would be possible to develop a scattering theory in the Lax-Phillips sense for the wave equation for $`p`$-forms on $`H^{n+1}`$, there would exist an incoming or an outgoing translation representation of $`(H,U(t)=e^{\overline{A}t})`$, that is there would exist a unitary operator

$$R:HL^2(,N),$$

such that

$$RU(t)R^1=T_0(t),$$

where $`N`$ is an auxiliary Hilbert space and $`T_0(t)`$ denotes the translation of $`t`$ units to the right. In turn, this would imply that $`i\overline{A}`$ has a purely absolutely continuous spectrum on the whole real axis, and moreover that the spectrum of $`i\overline{A}`$ has constant multiplicity. As a matter of fact, if there would exist an incoming or an outgoing translation representation, $`U(t)`$ on $`H`$ would be unitarily equivalent to $`T_0(t)`$ on $`L^2(,N)`$, hence its generator $`\overline{A}`$ would be unitarily equivalent to the operator $`\frac{d}{dx}`$ on $`L^2(,N)`$. As a consequence, $`i\overline{A}`$ would be unitarily equivalent to $`i\frac{d}{dx}`$, and hence, through a Fourier transform, to the multiplication operator

$$L^2(,N)L^2(,N)$$

$$f(x)xf(x),$$

whose spectrum is purely absolutely continuous, covers the whole real line and has constant multiplicity.

Since we have shown that the multiplicity of the spectrum of $`i\overline{A}`$ is not constant, we can conclude that there can be no incoming or outgoing translation representation for $`(H,U(t))`$. Hence:

###### Theorem 3.1.

It is not possible to develop a scattering theory following the Lax-Phillips approach for all $`p`$-forms on $`H^{n+1}`$.

Only in the cases $`n=2p`$ or $`n=2p2`$ one can hope that, restricting in some way the class of forms, it might be possible to develop a scattering theory for the wave equation in the Lax-Phillips sense.

## 4. Huygens’ principle

Consider the Cauchy problem (3.1) for the wave equation for $`p`$-forms, with smooth initial data $`\omega _1`$, $`\omega _2`$, and fix a point $`x_0H^{n+1}`$. We will say that Huygens’ principle holds if the following conditions are fulfilled:

1. if the initial data are zero in the ball $`B_{H^{n+1}}(x_0,R)`$ then the solution $`\omega (x,t)`$ at time $`t`$ of the initial value problem is zero in the ball $`B_{H^{n+1}}(x_0,Rt)`$;
2. if the initial data vanish outside the ball $`B_{H^{n+1}}(x_0,R)`$, then the solution $`\omega (x,t)`$ at time $`t`$ of the initial value problem vanishes in the ball $`B_{H^{n+1}}(x_0,|t|R)`$.

The spectral properties of the Laplace-Beltrami operator can be used in order to prove that Huygens’ principle fails in the case of the wave equation for $`p`$-forms on $`H^{n+1}`$, when $`p\frac{n+1}{2}`$.

Consider the Cauchy problem (3.1) for the wave equation for $`p`$-forms on $`H^{n+1}`$, with initial data $`\omega _1,\omega _2`$ smooth and compactly supported.

As noticed before, the operator $`A`$ defined on $`\mathrm{\Lambda }_c^p(H^{n+1})\mathrm{\Lambda }_c^p(H^{n+1})`$ by (3.4) is essentially skewselfadjoint in $`H`$ endowed with the energy norm, and its closure $`\overline{A}`$ is the generator of a one-parameter group of unitary operators $`U(t)`$ on $`H`$ which describes the solution of the wave equation for $`p`$-forms on $`H^{n+1}`$.

Now, let $`D_+`$ be the subspace of $`H`$ of those smooth initial data $`(\omega _1,\omega _2)`$ such that for every $`t>0`$

$$\{\begin{array}{c}(d\omega _1)(x,t)=0\hfill \\ (\delta \omega _1)(x,t)=0\hfill \\ \omega _2(x,t)=0\hfill \end{array}$$

for $`d_{H^{n+1}}(x,x_0)<t`$, and

$$(\omega _1(x,t),\omega _2(x,t))=(U(t)(\omega _1,\omega _2))(x).$$

The closure $`\overline{D_+}`$ of $`D_+`$ is a closed subspace of $`H`$ such that

$$U(t)\overline{D_+}\overline{D_+}$$

for every $`t>0`$. To show that

(4.1)
$$\underset{t}{}U(t)\overline{D_+}=\{\underset{¯}{0}\}$$

note that, if $`(\omega ,\eta )U(t)\overline{D_+}`$, for $`t>0`$, there exists a sequence $`\{(\omega _n,\eta _n)\}D_+`$ such that

$$U(t)(\omega _n,\eta _n)(\omega ,\eta )$$

in $`H`$, that is:

$$\{\begin{array}{c}d\omega _nd\omega \hfill \\ \delta \omega _n\delta \omega \hfill \\ \eta _n\eta \hfill \end{array}$$

Now,

$$\{\begin{array}{c}(d\omega _n)(x,t)=0\hfill \\ (\delta \omega _n)(x,t)=0\hfill \\ \eta _n(x,t)=0\hfill \end{array}$$

when $`d_{H^{n+1}}(x,x_0)<t`$, and this implies:

$$\{\begin{array}{c}(d\omega )(x,t)=0\hfill \\ (\delta \omega )(x,t)=0\hfill \\ \eta (x,t)=0\hfill \end{array}$$

when $`d_{H^{n+1}}(x,x_0)<t`$. Therefore, if $`(\omega ,\eta )_tU(t)\overline{D_+}`$, then

$$\{\begin{array}{c}(d\omega )(x,t)=0\hfill \\ (\delta \omega )(x,t)=0\hfill \\ \eta (x,t)=0\hfill \end{array}$$

when $`d_{H^{n+1}}(x,x_0)<t`$ for every $`t>0`$, and this implies $`(\omega ,\eta )=\underset{¯}{0}`$ in $`H`$, i.e. that 4.1 holds.

Now, if Huygens’ principle would hold, then

$$\overline{\underset{t}{}U(t)\overline{D_+}}=H.$$

In fact, if $`(\omega ,\eta )=0`$ for $`d_{H^{n+1}}(x_0,x)>R`$, then, by Huygens’ principle,

$$U(t)(\omega ,\eta )=0\text{for}d_{H^{n+1}}(x_0,x)<tR$$

for every $`t>0`$. In particular,

$$U(t+R)(\omega ,\eta )=0\text{for}d_{H^{n+1}}(x_0,x)<t+RR=t$$

for every $`t>0`$. This in turn would imply

$$U(R)(\omega ,\eta )D_+.$$

Thus, $`U(R)D_+`$ would contain the set of all initial data with support in $`B_{H^{n+1}}(x_0,R)`$, and hence

$$\underset{t}{}U(t)D_+\underset{t}{}U(t)\overline{D_+}$$

would contain all compactly supported data, which are dense in $`H`$.

Therefore, if Huygens’ principle would hold, $`D_+`$ would be an outgoing subspace for $`U(t)`$, and, thanks to Sinai’s theorem, there would exist an outgoing translation representation for $`H`$, $`U(t)`$. Now, this contradicts the properties of the spectrum of $`\mathrm{\Delta }_p`$. In conclusion:

###### Theorem 4.1.

Huygens’ principle does not hold for the wave equation for $`p`$-forms on $`H^{n+1}`$, when $`p\frac{n+1}{2}`$.

## 5. Coclosed case

In the case $`p=\frac{n}{2}`$, the spectrum of the Laplace-Beltrami operator is $`[0,+\mathrm{})`$, with non constant multiplicity. This fact makes it conceivable that, under some restrictions on the forms considered, a scattering theory in the Lax-Phillips sense could be developed for the wave equation.

This is indeed the case, if the class of coclosed $`p`$-forms, with $`n=2p`$, in $`H^{n+1}`$ is considered.

Indeed, it has been shown in  that Huygens’ principle holds for the wave equation for coclosed $`p`$-forms on $`H^{n+1}`$, when $`n=2p`$.

As a consequence:

###### Theorem 5.1.

Incoming and outgoing subspaces exist for the wave equation for coclosed $`p`$-forms on $`H^{n+1}`$, when $`p=\frac{n}{2}`$.

Proof. Consider the space of coclosed smooth, compactly supported $`p`$-forms in $`H^{n+1}`$, with $`n=2p`$, and the Cauchy problem (3.1) for the wave equation with coclosed smooth initial data $`\omega _1`$, $`\omega _2`$. The condition

$$\delta \omega _1(x)=\delta \omega _2(x)=0$$

is equivalent to

$$t\delta \omega (t,x)=\delta \omega _t(t,x)=0,$$

where $`\omega (x,t)`$ is the solution of the Cauchy problem corresponding to the initial data $`\omega _1,\omega _2`$.

Indeed,if the initial data $`\omega _1,\omega _2`$ are coclosed and $`\omega (x,t)`$ denotes the solution of the corresponding Cauchy problem, then $`\delta \omega (x,t)`$ is the solution of the Cauchy problem for the wave equation for $`(p1)`$-forms, with initial data $`\delta \omega _1,\delta \omega _2`$. Now, thanks to the uniqueness of the solution of the Cauchy problem for the wave equation for $`(p1)`$-forms, if $`\delta \omega _1=0,\delta \omega _2=0`$ then

$$\delta \omega (x,t)=0$$

for every $`t`$.

Consider now the Cauchy problem for the wave equation in the space of coclosed forms. The energy left invariant by the system is expressed by the norm

$$(\omega ,\eta )_{\stackrel{~}{H}}^2=<\mathrm{\Delta }_p\omega ,\omega >+\eta ^2=d\omega ^2+\eta ^20,$$

which is induced by the scalar product (3.3) restricted to the set of coclosed forms. The Hilbert space $`\stackrel{~}{H}`$, completion of the space $`D`$ of all smooth, compactly supported coclosed data with respect to this norm, is a closed subspace of $`H`$.

We will denote by $`\stackrel{~}{A}`$ the restriction of the operator $`A`$, expressed by (3.4), to $`D`$. $`\stackrel{~}{A}`$ maps $`D`$ into $`D`$, and is essentially skewselfadjoint in $`\stackrel{~}{H}`$. Consequently, it is the generator of a one-parameter group of unitary operators $`\stackrel{~}{U}(t)`$, which describes the evolution of the system, mapping initial data to the corresponding solution at time $`t`$.

Let now $`\stackrel{~}{D}_+`$ be the subspace of $`\stackrel{~}{H}`$ of those initial data $`(\omega _1,\omega _2)`$ such that for every $`t>0`$

$$\{\begin{array}{c}(d\omega _1)(x,t)=0\hfill \\ \omega _2(x,t)=0\hfill \end{array}$$

for $`d_{H^{n+1}}(x,x_0)<t`$, where $`x_0`$ is a fixed point of $`H^{n+1}`$ and

$$(\omega _1(x,t),\omega _2(x,t))=(\stackrel{~}{U}(t)(\omega _1,\omega _2))(x).$$

Its closure $`\overline{\stackrel{~}{D}}_+`$ in $`\stackrel{~}{H}`$ is a closed subspace of $`\stackrel{~}{H}`$ such that

$$\stackrel{~}{U}(t)\overline{\stackrel{~}{D}_+}\overline{\stackrel{~}{D}_+}.$$

Moreover, it satisfies

$$\underset{t}{}\stackrel{~}{U}(t)\overline{\stackrel{~}{D}_+}=\{\underset{¯}{0}\}.$$

Finally, the validity of Huygens’ principle implies that

$$\overline{\underset{t}{}\stackrel{~}{U}(t)\overline{\stackrel{~}{D}_+}}=\stackrel{~}{H},$$

because, if Huygens’principle holds, then, for every $`R>0`$, $`\stackrel{~}{U}(R)\stackrel{~}{D}_+`$ contains every smooth compactly supported coclosed data with support in $`B(x_0,R)`$. Hence, $`\overline{\stackrel{~}{D}_+}`$ is an outgoing subspace for $`\stackrel{~}{U}(t)`$ in $`\stackrel{~}{H}`$. Analogously, the subspace $`\overline{\stackrel{~}{D}_{}}`$ of $`\stackrel{~}{H}`$, consisting of those initial data $`(\omega _1,\omega _2)`$ such that for every $`t<0`$

$$\{\begin{array}{c}(d\omega _1)(x,t)=0\hfill \\ \omega _2(x,t)=0\hfill \end{array}$$

for $`d_{H^{n+1}}(x,x_0)<t`$, can be shown to be an incoming subspace for $`\stackrel{~}{U}(t)`$ in $`\stackrel{~}{H}`$.∎

Theorem 5.1 implies that the spectrum of $`i\overline{\stackrel{~}{A}}`$ is purely absolutely continuous, with constant multiplicity, and equal to $`(\mathrm{},+\mathrm{})`$. As a consequence:

###### Theorem 5.2.

When $`n=2p`$, the spectrum of the Laplace-Beltrami operator on coclosed $`p`$-forms on $`H^{n+1}`$ is purely absolutely continuous, with constant multiplicity, equal to $`[0,+\mathrm{})`$.