text:
# On J. Goodman’s comment to Language Trees and Zipping

## I Introduction

In this short note we reply to J. Goodman’s comments on a paper we have recently published (“Language Trees and Zipping” PRL ).

First of all we wish to apologize with the users of electronic archives for an use of them that we judge improper. Electronic archives should be used to spread in a fast way scientific knowledge and surely they are not intended as a large audience where raising polemics whose scientific content is at least doubtful. We shall accordingly try to keep the discussion on a scientific basis which could add some technical information for the reader genuinely interested in knowing whether our method could be used and for which purposes. Nevertheless we shall just point out some inconsistencies in the Goodman’s rationale.

In a recent cond-mat submission JG-cond-mat  (anticipated by an Open letter to the editors of Physical Review Letters that we report in Appendix) J. Goodman strongly criticizes our paper on the basis of three main points:

* The paper does not relates to physics and it should not be published;
* Goodman reviews the literature to conclude that ideas like the ones discussed in our paper are generically not novel;
* Goodman provides experimental results which apparently discredit our approach and our methods.

In this Letter we shall mainly discuss in detail the last point raised by Goodman, i.e. the technical point and probably the only one deserving attention from a scientific community. As for the other two points Goodman raises let us just notice the following.

The question whether a paper unrelated to physics should or should not be published on a journal like Physical Review Letters is beyond our judging possibilities. Probably the criticism should be addressed more to the Editors of Physical Review Letters (see the Open Letter to them in Appendix). Nevertheless Physical Review Letters has a section named Interdisciplinary Physics: Biological Physics, Quantum Information, etc. where we naively thought our paper could fit. If then one considers that publication on Physical Review Letters is subject to a tough peer review probably we could not be blamed for having attempted this submission.

Goodman says:

Physics journals should not be publishing articles that have nothing to do with physics. Of course, it is completely reasonable to publish applications of physics to other fields (both because this alerts other physicists to the possibility of applying their knowledge and because those in the field of interest may have difficulty understanding the terminology or techniques). It is also completely reasonable to publish the use of non-physics techniques applied to physics in a physics journal. But this paper applies computer science techniques (gzip!) to computer science problems. This seems extremely inappropriate. One might argue that the paper discusses entropy, a concept taken from physics. But the concept was taken from physics 50 years ago, at the dawn of computer science, and there is nothing physics-specific in the use of entropy in this paper; indeed, the use is entirely in the information theory/language modeling/computer science meaning of the word.

Except from learning that the concept of entropy now belongs to Computer Science (sic !), there is one important point to stress: our method provides a general framework to analyze for instance a generic time-series: could one imagine a more physical application than the study of predictability in dynamical systems or any other phenomenon, e.g. experimental investigations, where one typically has access to the system only through a measuring device which produces a time record of a certain observable, i.e. a sequence of data?

Apart from this there is something obscure in Goodman’s comment. There are typically several reasons to criticize a paper: (a) the paper is wrong: there are evident mistakes; (b) the paper is not publishable because the results it reports have already been obtained and published somewhere else (and in this case one reports the relevant references). If the paper is correct but its relevance is limited one typically does not waste his time trying to put in a black light something that the community will anyway ignore.

From this point of view it is not clear what disturbs Goodman. Apparently the paper is correct from his point of view. Apparently he is not able to quote relevant references where a similar approach has been proposed. He is only able to say:

“The only idea in this paper that is not very well known in the field is the idea of actually using a standard compression tool to do the classification. Still, even this idea dates back (at least) to 1995, when Ken Lang and Rich Caruana tried out the idea for doing newsgroup classification. In the end though, they ended up using Naive Bayes classifiers Lang . They didn’t bother to publish the compression idea because they thought it was better viewed as an interesting thought than as a serious classification method. Still, the idea of using compress got around a bit: see an introductory tutorial by a well known practitioner in this area, Tom Mitchell Mitchell . Admittedly, however, this technique is not that widely known, because computer scientists don’t typically try anything that crude – in a couple of hours (or less), we can build from scratch tools that work better than this.” (we have only added the italic style).

Let us now come to the technical part of the paper discussing the last point raised by Goodman, namely the experimental results that apparently discredit our method. Here Goodman’s claim is as follows: “If the Benedetto, Caglioti, Loreto method fails miserably against the simplest method computer scientists could think of, then why all this noise about it?”.

## II Results on Classification by subject

In this section we report the results of the same experiments performed by J. Goodman using our method for the analysis of the newsgroup messages available from http://www.ai.mit.edu/~jrennie/20Newsgroups. (Exactly the same version used by Goodman containing 18828 messages).

It is important to stress that Goodman is criticizing our paper in a very odd way. He claims that our technique fails on an experiment (the subject classification) that was not described in our paper. In his position we would have criticized the paper on its content and not on something completely different. We are nevertheless able to perform subject classification but we have never published how the technique works in this case. It is then quite natural for us to wonder what precisely Goodman did when he says: “For zipping, I applied the algorithm of PRL ”.

On the other hand, supposing Goodman used the technique described in the paper for authorship attribution, it is hard for us to understand how it is possible to make our method working so badly. We have nevertheless accepted the challenge and analyzed the same newsgroup messages as Goodman.

The technique we use (independently of the zipper considered) is as follows:

* For each unknown texts, consider the collection of texts obtained attaching it to each one of the texts considered as reference;
* Zip the resulting files as well as all the reference files separately and compute, for each reference file the difference between the size (in bytes) of the new total file zipped and the reference file zipped;
* Write for each unknown text the ranking of all the reference texts. Once one has the entire ranking for each unknown text one can start the attribution of each unknown text to the most pertinent category. Several paths are of course possible in agreement that with the fact that any classification has some level of arbitrary. The perfect classification does not exist: the proof lies in the fact that one can always refine the list of categories ending up with one specific category for each specific item in the list.

We use the criterion to associate the unknown text to the category of the reference text that was ranked as first in its particular ranking (as described for authorship attribution in our paper). We have performed the analysis using several zippers in order to compare the results and show the robustness of the approach. We have used in particular the following data compression schemes: gzip, a modified version of gzip (denoted as bcl) which differs from gzip for the fact of simply reading the second file instead of zipping it and on the Huffman optimization long . We have also used the method proposed by Merhav and Ziv merhav-ziv  to measure the relative entropy between two sequences of characters. The results are shown in Table I for a sample of $`200`$ newsgroup messages chosen in a random way.

As a matter of fact, using only the information coming from the document ranked in the first position (and thus not exploiting the entire information of the ranking list) the three zippers display the same efficiency Goodman obtained using Naive Bayesian Methods and not, as Goodman claimed, an efficiency three times smaller. As for the computational time we agree with Goodman that data compression techniques are slower than other procedures but here we are interested in testing the efficiency of the algorithms for a large variety of problems (not only newsgroup classification). We think that the question of speed of the algorithms should be addressed once one has identified the best algorithm (from the point of view of the efficiency) for a large set of problems. It would interesting for instance to know how Naive Bayesian methods perform in the set of problems we have considered in PRL . For instance in the case of language classification, having only one short text available per language, it is not clear how a Naive Bayesian method could sort out a measure of remoteness between pairs of texts to be used after for a tree representation. In this case we can guess that the method would suffer the lack of long training texts and this could also be true for the problem of authorship recognition where the number of texts available per author is very limited.

## III Conclusions

Before concluding several remarks are in order. The results shown above demonstrate as algorithms based on data compression techniques are highly competitive. It should be stressed how Goodman did not performed the same experiments we report in PRL  but he tried to discredit our approach performing a very focused test on one particular problem: the newsgroup classification. Unfortunately in that particular case his claim was not correct and our efficiency equals that of Naive Bayesian Approaches for this specific problem. It is worth to stress how on much more controlled and clean corpora the efficiency of our method raises considerably. We refer the reader to a long paper long  where we discuss in detail all the issues simply sketched in PRL  (as well as the classification by subject problem applied to large law corpora) and we compare the results obtained with different data compression schemes: gzip, bzip2, compress, merhav-ziv, etc.

From our point of view the situation could be interpreted as follows: it is not correct to neglect some particular technique on the only basis that it is the simplest technique a computer scientist can think of. Our point of view is that the use of data compression techniques represents a very general tool in several areas and not only in the narrow field of newsgroup classification. Always remaining in the field of computational linguistics there have been several contributions showing how data compression techniques could be useful in solving different problems (an incomplete list would include bell ; teahan ; juola ; el-yaniv ; Kontoyiannis ; nevill-manning ; thaper ): language recognition, authorship recognition, language classification, classification of large corpora by subject, etc. In this spirit our paper PRL  represents only one possible example which could hopefully give a contribution to this challenging field.

It is also important to stress how data compression oriented techniques only represent one sector in the large area of language modeling for the solution of linguistic motivated problems. Also the methodologies using ”entropic” or statistical concepts are very many and used in very different contexts: n-grams, hidden Markov modeling, Bayesian approaches, context-free grammars, neural networks just to quote some existing approaches. We refer to rosenfeld  for a recent overview about statistical language modeling.

Of course the possibilities of data-compression based methods go beyond computational linguistics. More in general since these methods apply to generic sequences of characters, we are not limited to consider only texts (and linguistic applications) but we can consider other cases where the strings of characters represent a different kind of coding: time sequences, genetic sequences (DNA or proteins) etc. These features are potentially very important for fields where the human intuition can fail: DNA and protein sequences, geological time series, stock market data, medical monitoring, etc. In summary the fields of application range from time-series analysis (we refer the reader to yak  for an application to time-series generated by some specific dynamical systems), to theory of dynamical systems benci , genetic problems (here also an incomplete list would include li ; grumbach ; Loewenstern ). One of the advantages of using a data compression scheme is that in this case one can face all the problems mentioned above in the same framework.

Appendix I

An open letter to the editors of Physical Review Letters

(The letter appeared only for a few days on the site: http://research.microsoft.com/~joshuago/ and then mysteriously disappeared, substituted by a more serious comment sent to the Editors of Physical Review Letters as well as to cond-mat archives.)

```

Date: Tue, 12 Feb 2002 15:17:12 -0800
From: Joshua Goodman <joshuago@microsoft.com>
To: benedetto@mat.uniroma1.it, caglioti@mat.uniroma1.it, loreto@roma1.infn.it
Subject: Open Letter to the editors of Physical Review Letters

An open letter to the editors of Physical Review Letters

Dears Sirs,

I wish to commend you on your open-minded interdisciplinary approach

to science.   In particular, I was extremely impressed that you

published an article on Computational Linguistics and Machine Learning

in a journal ostensibly devoted to physics ("Language Trees and

Zipping" by Dario Benedetto, Emanuele Caglioti, and Vittorio Loreto,

28 January 2002, (available at

http://babbage.sissa.it/abs/cond-mat/0108530) with the following

Abstract:

   In this Letter we present a very general method for extracting

   information from a generic string of characters, e.g., a text, a DNA

   sequence, or a time series. Based on data-compression techniques, its

   key point is the computation of a suitable measure of the remoteness

   of two bodies of knowledge. We present the implementation of the

   method to linguistic motivated problems, featuring highly accurate

   results for language recognition, authorship attribution, and language

   classification.

I myself attempted to publish the following letter on Physics in

several journals, including Computational Linguistics, the Machine

Learning Journal, and Computer Speech and Language:

"Measurement of Gravitational Constant Using Household Items" by

Joshua Goodman, unpublished manuscript.

   In this Letter we present a method for measuring the gravitational

   constant ("g") using household items, such as old shoes and a

   microwave oven.   Using the timer on the microwave oven, we were able

   to measure the time for the old shoe to reach the ground from several

   different heights.   We were able to determine that "g" equals 34 feet

   per second per second, plus or minus 3.   This is within 10% of the

   current state of the art in physics.   We noticed that our curves were

   not quite parabolic, and we speculate that this is evidence of a

   wormhole in the region of the experiment.

One of the reviews I received was very positive ("I commend you on

this excellent research.   Some might argue that it is irrelevant to a

journal such as ours, but physics, and ’g’ in particular, effect us

each and every day, for instance, as we walk or pour coffee.

Furthermore, since almost all readers of this journal have had a high

school (or even college) physics class, and all of us have and use

household items such as these, it is sure to be of interest to the

readership.") Unfortunately, the other reviewers were far more

critical.   One reviewer wrote "While it is clear that physicists

should be writing papers about natural language processing and

statistical modeling, this is because they are so much smarter than we

computer scientists are.   The reverse is obviously not true."   Another

reviewer was even too lazy to read the paper, writing "As everyone

knows, it takes dozens (or even hundreds) of authors to do any good

work in physics.   I find it difficult to believe that a single

authored letter is worth reading."   And another wrote "I remember from

my coursework something about Newton’s Method and Gradient Descent.

How come neither of these are cited in this paper?"   With comments

such as those, it seems that reviewers in our field are not even

competent to review physics papers, while clearly yours know enough to

review computer science papers.

Obviously, those in my field are not nearly as broad minded or well

rounded as those in yours.   I commend you on your openness, and hope

that you will soon expand to cover other areas of interest, such as

economics, psychology, and even literature.

Sincerely yours,

Joshua Goodman

----------------------------------------------------------------------

```


format: text

text:
# Pressure dependence of the thermoelectric power of single-walled carbon nanotubes

(March 16, 2024)

## Abstract

We have measured the thermoelectric power ($`S`$) of high purity single-walled carbon nanotube mats as a function of temperature at various hydrostatic pressures up to 2.0 GPa. The thermoelectric power is positive, and it increases in a monotonic way with increasing temperature for all pressures. The low temperature $`(\mathrm{T}<40\mathrm{K})`$ linear thermoelectric power is pressure independent and is characteristic for metallic nanotubes. At higher temperatures it is enhanced and though $`S(\mathrm{T})`$ is linear again above about 100 K it has a nonzero intercept. This enhancement is strongly pressure dependent and is related to the change of the phonon population with hydrostatic pressure.

Soon after the discovery of multi-walled carbon nanotubes (MWNT)1  single–walled carbon nanotubes (SWNT) were synthesised with catalytic metal particles.2  SWNT span electronic properties from metals to variable gap semiconductors, depending on the chirality and the diameter of the rolled up graphene sheet.3  It was immediately realized that metallic tubes represent the ultimate one–dimensional conductor, and they generated considerable interest due to the possible realization of the Luttinger liquid (LL) behaviour in nature. Indeed, conductivity measurements performed on individual nanotubes have shown the canonical power–law behaviour of conductance as function of temperature and voltage.4

Early thermoelectric power (TEP) experiments performed on nanotube bundles lead to qualitatively similar results:5 -8  a metallic character due to hole carriers, a knee–feature (change of slope) around 100 K, and preparation dependent magnitude. The observed TEP was associated with intrinsic electronic behaviour, and the fact that the nanotube bundles contain both metallic and semiconducting ropes allowed a description of the unusual shape of the $`S(\mathrm{T})`$ curves in a simple two-band model, using a sufficient number of fitting parameters. It was also realized that in the high temperature range the thermoelectric power is significantly larger then expected for bundles constituted predominantly of metallic ropes.5

Since SWNTs may be contaminated by magnetic impurities coming from the catalysts (Fe, Ni, Co) needed for their growth, the high value of the TEP above 100 K questions its intrinsic origin. The influence of the magnetic impurities has clearly been experimentally demonstrated:9  in the presence of transition metal impurities the thermopower has a strong peak, just around 100 K. This peak was associated with the Kondo scattering of conduction electrons on the impurity spins. The enhancement of the TEP could be largely, albeit not completely, removed by iodine neutralization.9  A non–intrinsic, residual ”Kondo–term” may arise from magnetic impurities trapped on the body of the nanotubes, even in case of nominally pure samples.

The change of sign of the thermoelectric power upon degassing also raised the possibility of a non-intrinsic process.10  The conductivity of doped semiconducting tubes can be as high as that of the metallic ones and the overall thermopower may be determined by their contribution. Recently it was suggested that hole–doping due to oxygen adsorption explains the sign and the shape of $`S(\mathrm{T})`$, at least in case of bundles where semiconducting ropes dominate.10

In this paper we report the results of thermoelectric power experiments on high purity single–wall nanotubes under hydrostatic pressures. We show that the low–temperature linear thermoelectric power arises from metallic ropes and it is consistent with the Luttinger liquid behaviour. Our results on pressurized samples exclude the ”residual Kondo” origin of the high temperature TEP, and allow to separate the band and phonon contributions. We propose a model, in which the knee–feature is coming from the change of the population of phonons that are relevant for the electron scattering. This picture is further supported by comparison of the thermoelectric power in carbon–based metallic compounds, which have a similar local bond structure but a quite different band structure.

SWNT soot was purchased from Rice University, rich in catalytic particles. Most of these metallic blobs, covered with a graphitic shell, were eliminated using a purification method elaborated in our laboratory.11  This purification method results in open and less bundled, more dispersed nanotubes. The purified SWNTs from the suspension were deposited on a filtration membrane and a thin film of buckypaper was produced. The buckypaper was heated in vacuum to $`1200^{}\mathrm{C}`$ in order to eliminate absorbates. A thin slab of $`3.0\times 1.0\times 0.01\mathrm{mm}^3`$ was cut and mounted on a thermopower sample holder, which fits into a clamped pressure cell. Small metallic heaters installed at both ends of the sample generated the temperature gradient, which was measured with a Chromel–Constantan differential thermocouple. The pressure medium used in this study is kerosene and the maximum pressure is 2.0 GPa. The pressure was measured using a calibrated InSb pressure gauge. A similar sample arrangement and pressure cell was used in the study of the pressure dependence of the resistivity of SWNTs by Gaal et al.12

Figure 1 displays the resistance and the thermopower at ambient pressure in the $`\mathrm{T}=4300\mathrm{K}`$ region. The resistance has a negative slope in the whole temperature range, which is characteristic of an undoped sample. The low temperature increase is consistent with the signature of the LL feature.4  At low temperatures, below about 40 K, the thermopower is linear (see also the inset of Figure figure2). The well–resolved linear TEP excludes the possibility of a gapped density of states, or hopping conduction trough mid–gap impurity levels, which might be invoked to account for the low temperature upturn in the resistance. At higher temperatures the thermoelectric power is enhanced, then its slope changes around 100 K. These features are in good qualitative agreement with earlier studies on pure nanotubes.5 -8

The temperature dependence of the TEP measured at various pressures is shown in Figure 2. We found that the low–temperature linear part is not sensitive to the application of pressure, within the experimental resolution. In contrast at high temperatures the TEP is strongly pressure dependent. The distinct behaviour at low and high temperatures suggests that the thermoelectric power has different contributions in the two ranges. This observation excludes the possibility that the thermoelectric power is related to residual magnetic impurities. Pressure does not influence the sample purity or the exchange coupling (J) between the conduction electrons and the magnetic impurity, thus it could modify a possible contribution of Kondo origin solely by shifting the Kondo temperature,

$$T_K\frac{E_F}{k_B}\mathrm{exp}\left(\frac{1}{2|J|D(E_F)}\right).$$
(1)
Here $`D(\epsilon _F)`$ is the electronic density of states at the Fermi level, which might have a pressure dependence. Such a pressure dependence should, however, also be reflected in the metallic low temperature range, which is not observed. Moreover, $`T_K`$ should decrease with increasing pressure as $`D(\epsilon _F)`$ decreases. Instead, one rather observes a smearing of the shoulder due to the application of the pressure than a structure, which shifts with pressure.

The pressure study also shows that characteristic $`S(\mathrm{T})`$ shape shown in Figure 1 and published by various groups5 -8  cannot arise from oxygen adsorption, either. In the doped semiconductor picture the TEP is interpreted in terms of the particular shape of electronic density of states. Any pressure induced modification of the density of states would rescale the $`S(\mathrm{T})`$ curve in the whole temperature range, in contrast to the observed differences in the low and high temperature ranges.

The linear temperature dependence at low temperatures is characteristic of a simple metal, obeying Fermi liquid theory. Recent theoretical studies13 ; 14  of the thermal transport in the Luttinger liquid phase showed, that – though both the electrical and the thermal conduction are unconventional – the low temperature thermopower still follows Mott’s formula of the logarithmic derivative of the conductivity, $`\sigma `$ ,

$$S=\frac{\pi ^2}{3}\frac{k_B^2}{e}\left(\frac{\mathrm{ln}\sigma (\epsilon )}{\epsilon }\right)_{\epsilon =\epsilon _F}.$$
(2)
This leads to $`S=cT`$, whit a non-universal coefficient, c, which is determined by the band shape and the scattering rate. For the Luttinger liquid of spatial extension L the (2) holds only for low temperatures, while an enhancement is expected for $`k_BT>\mathrm{}v_F/gL`$ ($`v_F`$ is the Fermi velocity, while $`g`$ is the correlation parameter of the LL).14  We do not believe this is the source of the observed enhancement in our case, since for realistic values of defect free nanotube segments $`(\mathrm{L}<100\mathrm{nm})`$ the corresponding temperature is above the range covered by the experiments.

The attenuation of the shoulder in $`S(\mathrm{T})`$ with the pressure puts forward another explanation for the overall temperature dependence.8 ; 15  If inelastic scattering of electrons on phonons is of importance, then the thermoelectric power is expressed as

$$S=S_0(1+\lambda \lambda _S(\mathrm{T})),$$
(3)
where $`S_0`$ is given by (2), $`\lambda `$ is the standard mass enhancement due to renormalization by phonons and $`\lambda _S(\mathrm{T})`$ is the temperature–dependent enhancement of the thermopower. The latter quantity contains the Eliashberg function, and if the phonon density of states has a peak then $`\lambda _S(\mathrm{T})`$ will have a shoulder at the corresponding temperatures. This term gives a non–zero intercept of $`S(\mathrm{T})`$ when extrapolated from high temperatures. Furthermore, this term can fade away thus changing the slope of $`S`$ when the population of phonons, which strongly scatter the charge carriers, diminishes.

The model can also account for the observed pressure dependence. High pressure increases the force constant acting between the carbon atoms, shifting any feature in the phonon density of states to higher energies. Consequently, the population of the phonons, which strongly interact with the charge carriers, decreases and the thermopower drops. This explanation is corroborated by the fact that the temperature at which the shoulder in $`S`$ develops somewhat increases towards higher temperatures at high pressures.

The phonon mode in the case of SWNT should be the vibration of the carbon atoms in the hexagons, since similar feature in $`S`$ could be observed in other carbon–based materials, which are built up from the benzyl–ring units. Figure 3 summarizes the Seebeck coefficient of intercalated graphite,16  carbon nanohorns,17  alkali metal doped C70,18  and SWNTs. All the samples show similar break in slope in the 100 K range.

In conclusion, we have performed thermopower measurements on purified SWNT mat as a function of temperature and hydrostatic pressure. While the low temperature metallic thermopower was found to be pressure independent, the high temperature unusual TEP contribution was strongly suppressed by the application of pressure. We showed that the thermopower of our sample is due to metallic nanotubes, and we excluded non–intrinsic effects, like Kondo scattering of residual impurities or doping by adsorbed oxygen. Our results suggest that the temperature and pressure dependencies are governed by the change of population of phonons, connected with the thermal excitation of the hexagons (the building units of the SWNTs).

###### Acknowledgements.

The sample purification by Klara Hernadi and useful discussions with Walter Schirmacher are acknowledged. This work has been supported by the Swiss National Foundation for Scientific Research and the NCCR Pool ”MaNEP”.


text:
# A DETECTOR OF GRAVITATIONAL WAVES BASED ON COUPLED MICROWAVE CAVITIES

## I Introduction

In a series of papers it was studied how the effects due to the interaction between the gravitational and the electromagnetic fields could be used to detect gravitational waves bppr1 ; bppr2 . The proposed detector exploits the energy transfer induced by the gravitational wave between two levels of an electromagnetic resonator, whose frequencies $`\omega _1`$ and $`\omega _2`$ are both much larger than the angular frequency $`\mathrm{\Omega }`$ of the g.w. and satisfy the resonance conditon $`|\omega _2\omega _1|=\mathrm{\Omega }`$. 111The interaction between the g.w. and the detector is characterized by a transfer of energy and of angular momentum. Since the elicity of the g.w., i.e. the angular momentum along the direction of propagation, is $`2`$, it can induce a transition between the two levels provided their angular momenta differ by $`2`$; this can be achived by putting the two cavities at right angle or by a suitable polarization of the electromagnetic field inside the resonator. In the scheme suggested by Bernard et al. the two levels are obtained by coupling two identical high frequency cavities222Throughout this paper, we shall call resonator the whole detector made up of two coupled cavities; e.g. we shall speak about one resonator composed by two coupled spherical cavities.; the angular frequency $`\omega _1`$ is the frequency of the level symmetrical in the fields of the two cavities, and $`\omega _2`$ is that of the antisymmetrical one. The frequency difference between the symmetric and the antisymmetric level is determined by the coupling, and can be adjusted by a careful resonator design. Since the detector sensitivity is proportional to the square of the resonator quality factor, superconducting cavities should be used for maximum sensitivity.

The power transfer between the levels of a resonator made up of two pill–box cavities, mounted end–to–end and coupled by a small circular aperture in their common endwall, was checked in a series of experiments by Melissinos et al., where the perturbation of the resonator volume was induced by a piezoelectric crystal rrm1 ; rrm2 . Recently the experiment was repeated by our group with an improved experimental set–up; we obtained a sensitivity to fractional deformations of the resonator length as small as $`\delta \mathrm{}/\mathrm{}10^{20}`$ Hz-1/2 rsi .

In this paper we shall discuss the mechanism of the interaction of a gravitational wave with a detector based on two coupled resonant cavities. In previous works this issue was discussed using the concept of a dielectric tensor associated with the gravitational wave pr1 . The interaction was analyzed in the reference frame where the resonator walls were at rest even in presence of a gravitational perturbation. We shall analyze the effect in the proper reference frame attached to the detector and we shall therefore consider the interaction between the wave and the field stored inside the resonator due to the coupling of the g.w. with the mechanical structure of the detector caves .

The paper will be organized as follows: in section II the problem of finding the electromagnetic fields in a closed volume with time–varying boundary conditions is studied, and an approximate expression of the normal modes in a perturbed resonator is worked out. In section III we shall analyze the interaction of a g.w. with the mechanical structure of the detector; we shall see that the transfer of energy between a mechanical and an electromagnetic oscillation depends both on the electromagnetic field distribution inside the resonator and on the resonator geometry and mechanical properties. In section IV we shall discuss some ideas for the developement of a realistic gravitational wave detector based on spherical microwave cavities. Afterwords the coupled equations of motion for the fields in a perturbed resonator are worked out and solved. In section VIII the issue of the thermal noise of the detector’s walls is studied; other noise contribution are also considered. Finally, in the last section, the expected sensitivity of some detector configuration is shown and discussed.

## II Electromagnetic field in a resonator with perturbed boundaries

To study the mechanism of the energy transfer between the two levels of an electromagnetic resonator perturbed by a gravitational wave we shall follow classic electromagnetic theory. We shall make use of the fact that any field configuration inside the resonator can be expressed as the superposition of the electromagnetic normal modes of the given resonator slater .

If no sources are present, the electromagnetic field in vacuum is determined by the equations:

$$\stackrel{}{}\stackrel{}{E}=0$$
(1)
$$\stackrel{}{}\stackrel{}{H}=0$$
(2)
$$\stackrel{}{}\stackrel{}{E}+\mu _0\frac{\stackrel{}{H}}{t}=0$$
(3)
$$\stackrel{}{}\stackrel{}{H}ϵ_0\frac{\stackrel{}{E}}{t}=0$$
(4)
As can be easily verified, eqs. (1)–(4) are automatically satisfied if the fields satisfy the wave equations:

$$^2\stackrel{}{E}\frac{1}{c^2}\frac{^2\stackrel{}{E}}{t^2}=0$$
(5)
$$^2\stackrel{}{H}\frac{1}{c^2}\frac{^2\stackrel{}{H}}{t^2}=0$$
(6)
with $`c=(\mu _0ϵ_0)^{1/2}`$.

Let us assume that the field is contained in a resonator with perfectly conducting walls. If we impose boundary conditions on the fields, the solution of the wave equations (5)–(6) will have an infinite discrete set of normal–mode solutions orthogonal to one another and complete, in the sense that any arbitrary field in the resonator can be expressed as a sum of these normal modes with suitable amplitudes. The amplitudes of each mode can then be used to describe the field in the resonator.

By the familiar procedure of separation of variables we may assume a solution of eqs. (5)–(6) of the form:

$$\stackrel{}{E}(\stackrel{}{r},t)=\underset{n=0}{\overset{\mathrm{}}{}}_n(t)\stackrel{}{E}_n(\stackrel{}{r})$$
(7)
and

$$\stackrel{}{H}(\stackrel{}{r},t)=\underset{n=0}{\overset{\mathrm{}}{}}_n(t)\stackrel{}{H}_n(\stackrel{}{r})$$
(8)
where we have defined:

$$_n(t)\sqrt{ϵ_0}_V\stackrel{}{E}\stackrel{}{E}_n𝑑V$$
(9)
and

$$_n(t)\sqrt{\mu _0}_V\stackrel{}{H}\stackrel{}{H}_n𝑑V$$
(10)
where the integrals are performed over the resonator volume.

We require that at the walls the tangential component of $`\stackrel{}{E}`$ and the normal component of $`\stackrel{}{H}`$ vanish. With this assumption the functions $`\stackrel{}{E}_n(\stackrel{}{r})`$ and $`\stackrel{}{H}_n(\stackrel{}{r})`$ satisfy the equations:

$`k_n\stackrel{}{E}_n=\stackrel{}{}\stackrel{}{H}_n`$ (11)
$`k_n\stackrel{}{H}_n=\stackrel{}{}\stackrel{}{E}_n`$ (12)
where $`k_n=\omega _n/c`$ is the propagation constant associated with the $`n_{th}`$ mode.

It can be proved that the normal modes $`\stackrel{}{E}_n`$ and $`\stackrel{}{H}_n`$ have orthogonality properties of the form

$`{\displaystyle _V}\stackrel{}{E}_n\stackrel{}{E}_m𝑑V`$ $`=`$ $`\delta _{nm}`$ (13)
$`{\displaystyle _V}\stackrel{}{H}_n\stackrel{}{H}_m𝑑V`$ $`=`$ $`\delta _{nm}`$ (14)
For a cubical resonator the functions $`\stackrel{}{E}_n(\stackrel{}{r})`$ and $`\stackrel{}{H}_n(\stackrel{}{r})`$ are $`\mathrm{sin}(\stackrel{}{k}_n\stackrel{}{r})`$ and $`\mathrm{cos}(\stackrel{}{k}_n\stackrel{}{r})`$. For other geometries they will be other complete sets of functions. The boundary conditions and geometry determine the different modes which are distinguished by the index $`n`$. In general three numbers are needed to specify a mode; $`n`$ is an abbreviation for this set of numbers.

Let us now expand the fields $`\stackrel{}{E}`$, $`\stackrel{}{H}`$ in terms of the orthogonal functions $`\stackrel{}{E}_n`$ and $`\stackrel{}{H}_n`$; when we substitute these expansions in Maxwell’s equations and equate coefficients, so as to get the differential equations satisfied by the various coefficients, we find that equations (1) and (2) are automatically satisfied. From equations (3) and (4) we find the following equations for the expansion coefficients:

$$\frac{d_n(t)}{dt}\omega _n_n(t)=\frac{\omega _n}{𝒬_n}_n(t)$$
(15)
$$\frac{d_n(t)}{dt}+\omega _n_n(t)=0$$
(16)
We have taken into account the dissipation arising from the finite conductivity of the walls introducing the electromagnetic quality factor slater :

$$𝒬_n=\frac{\omega _n\mu _0}{R_s}\frac{_VH_n^2𝑑V}{_SH_n^2𝑑S}=\frac{G_n}{R_s}$$
(17)
$`R_s`$ is the material–dependent surface resistance of the walls, and the geometric factor $`G_n`$ of the $`n_{th}`$ mode is defined as:

$$G_n=\omega _n\mu _0\frac{_VH_n^2𝑑V}{_SH_n^2𝑑S}$$
(18)
As can readily be seen, equations (15)–(16) for the field expansion coefficients are decoupled: the modes are independent from one another and behave as simple damped harmonic oscillators.

### II.1 Perturbation of boundaries

Let us suppose that the resonator’s boundary is perturbed so that the eigenvalues and eigenfunctions of the perturbed resonator differ but little from those of the original one. The perturbation method allows to find the eigenvalues and eigenfunctions of the perturbed problem from the knowledge of the original ones. It basically consists in expanding in power series of a perturbation parameter $`\sigma `$ the new modes and frequencies goubau :

$`\stackrel{}{E}_n^{}(\stackrel{}{r})`$ $`=`$ $`\stackrel{}{E}_n(\stackrel{}{r})+\sigma \stackrel{}{e}_n(\stackrel{}{r})+𝒪(\sigma ^2)`$
$`\stackrel{}{H}_n^{}(\stackrel{}{r})`$ $`=`$ $`\stackrel{}{H}_n(\stackrel{}{r})+\sigma \stackrel{}{h}_n(\stackrel{}{r})+𝒪(\sigma ^2)`$ (19)
$`k_n^{}`$ $`=`$ $`k_n+\sigma \kappa _n+𝒪(\sigma ^2)`$
In the perturbed resonator we shall have:

$$\stackrel{}{E}(\stackrel{}{r},t)=\underset{n=0}{\overset{\mathrm{}}{}}\stackrel{~}{}_n(t)\stackrel{}{E}_n^{}(\stackrel{}{r})$$
(20)
and

$$\stackrel{}{H}(\stackrel{}{r},t)=\underset{n=0}{\overset{\mathrm{}}{}}\stackrel{~}{}_n(t)\stackrel{}{H}_n^{}(\stackrel{}{r})$$
(21)
with:

$$\stackrel{~}{}_n(t)\sqrt{ϵ_0}_V^{}\stackrel{}{E}\stackrel{}{E}_n^{}𝑑V$$
(22)
and

$$\stackrel{~}{}_n(t)\sqrt{\mu _0}_V^{}\stackrel{}{H}\stackrel{}{H}_n^{}𝑑V$$
(23)
where the integrals are now performed over the perturbed volume $`V^{}`$.

The perturbed modes will satisfy the following equations:

$`k_n^{}\stackrel{}{E}_n^{}=\stackrel{}{}\stackrel{}{H}_n^{}`$ (24)
$`k_n^{}\stackrel{}{H}_n^{}=\stackrel{}{}\stackrel{}{E}_n^{}`$ (25)
When the expansions (II.1) are inserted into the equations (24) we obtain a series of equations determining the various terms of the expansion (II.1). In the first order approximation, the perturbed fields may be written:

$$\sigma \stackrel{}{e}_n(\stackrel{}{r})=\underset{m=1}{\overset{\mathrm{}}{}}A_{nm}\stackrel{}{E}_m(\stackrel{}{r})$$
(26)
$$\sigma \stackrel{}{h}_n(\stackrel{}{r})=\frac{1}{2}𝒞_{nn}\stackrel{}{H}_n(\stackrel{}{r})+\underset{m=1}{\overset{\mathrm{}}{}}B_{nm}\stackrel{}{H}_m(\stackrel{}{r})$$
(27)
$$\sigma \kappa _n=\frac{1}{2}k_n𝒞_{nn}$$
(28)
where the sums have to be performed for $`nm`$. The expansion coefficients have the form goubau :

$$A_{nm}=\frac{k_mk_n}{k_m^2k_n^2}𝒞_{nm}$$
(29)
$$B_{nm}=\frac{k_n^2}{k_m^2k_n^2}𝒞_{nm}$$
(30)
with

$$𝒞_{nm}=_𝒱(\stackrel{}{H}_n\stackrel{}{H}_m\stackrel{}{E}_n\stackrel{}{E}_m)𝑑V$$
(31)
where $`𝒱=V^{}V`$ is the (algebraic) difference between the perturbed and the original volume.333The derivation of the eigenmodes and eigenvalues of the perturbed resonator has been made assuming a static perturbation; in the following we shall assume that it is also correct for a perturbation which has a rate of change much slower than the e.m. field characteristic frequency.

The calculation of the coupling coefficient $`𝒞_{nm}`$ depends on how the resonator is deformed by an external force. For this reason in the next section we shall briefly review the study of the mechanical behaviour of a body under the influence of an external force.

## III Analysis of the mechanical response of the detector

The interaction of a gravitational wave with the mechanical structure of the detector can be studied by means of classical, non–relativistic, linear elasticity theory landau\_elas ; lobo . If $`\stackrel{}{u}(\stackrel{}{r})`$ denotes the displacement of the mass element at point $`\stackrel{}{r}`$, relative to the centre of mass of the body in its unperturbed state, and $`\stackrel{}{f}(\stackrel{}{r},t)`$ is the volume force density which acts on the body, the displacement is the solution of the system of partial differential equations:

$$\rho \frac{^2\stackrel{}{u}}{t^2}\mu ^2\stackrel{}{u}(\lambda +\mu )\stackrel{}{}(\stackrel{}{}\stackrel{}{u})=\stackrel{}{f}(\stackrel{}{r},t)$$
(32)
with suitable boundary and initial conditions. $`\rho (\stackrel{}{r})`$ is the mass density of the body and $`\lambda `$ and $`\mu `$ are the material’s elastic Lamé coefficients. In the following we shall adopt null initial conditions:

$$\stackrel{}{u}(\stackrel{}{r},0)=\frac{\stackrel{}{u}}{t}(\stackrel{}{r},0)=0$$
(33)
The expansion theorem saulson  states that the displacement of a system in response to an applied force is equal to the superposition of the normal modes $`\stackrel{}{\xi }_\alpha (\stackrel{}{r})`$ of the system444In this section we shall label with greek indices the mechanical normal modes of the system, and with latin indices the normal modes of the electromagnetic field stored inside the system.:

$$\stackrel{}{u}(\stackrel{}{r},t)=\underset{\alpha =1}{\overset{\mathrm{}}{}}\stackrel{}{\xi }_\alpha (\stackrel{}{r})q_\alpha (t)$$
(34)
The normal modes $`\stackrel{}{\xi }_\alpha (\stackrel{}{r})`$ are the eigen–solutions to

$$\mu ^2\stackrel{}{\xi }_\alpha +(\lambda +\mu )\stackrel{}{}(\stackrel{}{}\stackrel{}{\xi }_\alpha )=\omega _\alpha ^2\rho \stackrel{}{\xi }_\alpha $$
(35)
with boundary conditions; here $`\alpha `$ is an index, or set of indices, labelling the mode of frequency $`\omega _\alpha `$. The modes are normalized so that

$$_{Vol}\stackrel{}{\xi }_\alpha (\stackrel{}{r})\stackrel{}{\xi }_\beta (\stackrel{}{r})\rho (\stackrel{}{r})𝑑V=M_\alpha \delta _{\alpha \beta }$$
(36)
where $`M_\alpha `$ is the reduced mass of the $`\alpha `$ mode. For a homogeneous system $`M_\alpha M`$, where $`M`$ is the mass of the system.

$`q_\alpha (t)`$ is the generalized coordinate of the $`\alpha `$ mode, obeying the dynamical equation of motion:

$$\ddot{q}_\alpha (t)+\frac{\omega _\alpha }{Q_\alpha }\dot{q}_\alpha (t)+\omega _\alpha ^2q_\alpha (t)=\frac{f_\alpha (t)}{M}$$
(37)
where an empirical damping term, proportional to the system velocity, has been added; $`f_\alpha (t)`$ is the generalized force, given by

$$f_\alpha (t)=_{Vol}\stackrel{}{f}(\stackrel{}{r},t)\stackrel{}{\xi }_\alpha (\stackrel{}{r})𝑑V$$
(38)
The solution of equation (37), satisfying the initial conditions (33), can be written in term of a Green function integral as courant :

$$q_\alpha (t)=\frac{1}{M\omega _\alpha }_0^tf_\alpha (t^{})\mathrm{sin}\omega _\alpha (tt^{})\mathrm{exp}\frac{tt^{}}{\tau _\alpha }dt^{}$$
(39)
with $`\tau _\alpha =2Q_\alpha /\omega _\alpha `$ is the amplitude decay time of the system.

In the frequency domain the asymptotic solution of equation (37) is easily found to be:

$$q_\alpha (\omega )=\frac{f_\alpha (\omega )/M}{\omega _\alpha ^2\omega ^2+j\omega \omega _\alpha /Q_\alpha }$$
(40)
Substituting into eq. (34), we find

$$\stackrel{}{u}(\stackrel{}{r},\omega )=\underset{\alpha =1}{\overset{\mathrm{}}{}}\frac{f_\alpha (\omega )\stackrel{}{\xi }_\alpha (\stackrel{}{r})/M}{\omega _\alpha ^2\omega ^2+j\omega \omega _\alpha /Q_\alpha }$$
(41)
It is clear from eq. (41), that if we are interested in the displacement in a narrow frequency interval $`\omega \pm \delta \omega `$, only those modes for which $`\omega _\alpha \omega `$ (and $`f_\alpha 0`$), will give a significant contribution.

### III.1 Interaction of a g.w. with the mechanical structure of the detector

An incoming gravitational wave manifests itself as a tidal force density acting on the mechanical structure of the detector. Given the expression of the gravitational force and the mechanical properties of the detector, the resulting deformation can be calculated, with the aid of the mathematical apparatus outlined in the previous section.

We are mainly interested in the evaluation of the coupling coefficient $`𝒞_{ij}(t)`$. Let us note that, for small displacements, we can write the integral over the perturbed volume as a surface integral in the form (see fig. 1):

$`𝒞_{ij}(t)={\displaystyle _𝒱}(\stackrel{}{H}_i\stackrel{}{H}_j\stackrel{}{E}_i\stackrel{}{E}_j)𝑑V`$ (42)
$`{\displaystyle _S}(\stackrel{}{H}_i\stackrel{}{H}_j\stackrel{}{E}_i\stackrel{}{E}_j)\stackrel{}{u}(t)𝑑\stackrel{}{S}`$
where the integral in the r.h.s of eq. (42) is now performed over the unperturbed detector boundary. It is worth noting that this integral can be expressed as a superposition of the mechanical normal modes of the system. Using the expansion theorem (34) we can write:

$`𝒞_{ij}(t)={\displaystyle _S}(\stackrel{}{H}_i\stackrel{}{H}_j\stackrel{}{E}_i\stackrel{}{E}_j)\stackrel{}{u}(t)𝑑\stackrel{}{S}=`$ (43)
$`{\displaystyle \underset{\alpha =1}{\overset{\mathrm{}}{}}}q_\alpha (t){\displaystyle _S}(\stackrel{}{H}_i\stackrel{}{H}_j\stackrel{}{E}_i\stackrel{}{E}_j)\stackrel{}{\xi }_\alpha 𝑑\stackrel{}{S}={\displaystyle \underset{\alpha =1}{\overset{\mathrm{}}{}}}q_\alpha (t)C_{ij}^\alpha `$
where we have defined the time–independent form factor $`C_{ij}^\alpha `$ as:555We remind that the superscript $`\alpha `$ labels the mechanical normal mode, while the subscripts $`i`$ and $`j`$ label the electromagnetic modes.

$$C_{ij}^\alpha =_S(\stackrel{}{H}_i\stackrel{}{H}_j\stackrel{}{E}_i\stackrel{}{E}_j)\stackrel{}{\xi }_\alpha 𝑑\stackrel{}{S}$$
(44)
If the external force couples strongly only to one mechanical mode of the detector (say the $`m`$), we can write the simplified expression

$$𝒞_{ij}(t)=q_m(t)C_{ij}^m$$
(45)
In summary, to have an effective coupling between the two electromagnetic modes we need:

1. that the generalized coordinate $`q_m(t)`$ is different from zero. If the system is initially at rest, this is true if and only if the generalized force, is itself different from zero, as is shown in eqs. (37) and (39);
2. that the spatial integral in eq. (44) is different from zero. In section VI we shall discuss how this depends on the symmetries of the electromagnetic field and of the perturbed volume.

## IV Detector design

In order to build a realistic detector a suitable cavity shape has to be chosen. From quite general arguments a detector based on two coupled spherical cavites looks very promising (see fig. 2).

In order to approach the interesting frequency range for g.w. detection, the mode splitting (i.e. the detection frequency) will be $`\omega _2\omega _110`$ kHz. The internal radius of the spherical cavity will be $`r100`$ mm, corresponding to a frequency of the TE011 mode $`\omega 2`$ GHz. The overall system mass and length will be $`M5`$ kg and $`L0.8`$ m. The choice of these frequencies for the resonator and mode splitting will be also useful in order to test the feasibility of a detector working at $`200`$ MHz and at a detection frequency of $`1`$ KHz.

A tuning cell, or a superconducting bellow, will be inserted in the coupling tube between the two cavities, allowing to tune the coupling strength (i.e. the detection frequency) in a narrow range around the design value.

From the point of view of the electromagnetic design the spherical cell has the highest geometrical factor, and so the highest quality factor, for a given surface resistance. For the TE011 mode of a sphere the geometric factor $`G`$ has a value $`G850\mathrm{\Omega }`$, while for a standard elliptical accelerating cavity the TM010 mode has a value of $`G250\mathrm{\Omega }`$. Looking at the best reported values of quality factor of accelerating cavities, which typically are in the range $`10^{10}`$$`10^{11}`$, we can extrapolate that the quality factor of the TE011 mode of a spherical cavity can exceed $`𝒬10^{11}`$.

From the mechanical point of view it is well know that a sphere has the highest interaction cross-section with a g.w. and that only a few mechanical modes of the sphere do interact with a gravitational perturbation (the quadrupolar ones) lobo . The mechanical design is highly simplified if the spherical geometry is used since the deformation of the sphere is given by the superposition of just one or two normal modes of vibration and thus can be easily modeled. In fact the proposed detector acts essentially as a standard g.w. resonant bar detector: the gravitational perturbation interacts with the mechanical structure of the resonator, deforming it. The e.m. field stored inside the resonator is affected by the time–varying boundary conditions and a small quantity of energy is transferred from the initially excited e.m. mode to the initially empty one, provided the g.w. frequency equals the frequency difference of the two modes. We emphasize that our detector is sensitive to the polarization of the incoming gravitational signal: once the e.m. axis has been chosen inside the resonator, a g.w with polarization axes in the direction of the field axis will drive the energy transfer between the two modes of the cavity with maximum efficiency.

Finally the spherical cells can be esily deformed in order to remove the unwanted e.m. modes degeneracy and to induce the field polarization suitable for g.w. detection. The interaction between the stored e.m. field and the time-varying boundary conditions is not trivial and depends both on how the boundary is deformed by the external perturbation and on the spatial distribution of the fields inside the resonator, as shown by the expression of the coupling coefficient $`C_{21}^m`$ (eq. (44)). It has been calculated that the optimal field spatial distribution is with the field axis of the two cavities orthogonal to each other (see fig. 3). Different spatial distributions (e.g. with the field axis along the resonators’ axis) give a smaller effect or no effect at all. A more detailed discussion of the coupling coefficient calculation is done in section VI.

## V Field equations

We shall derive the equations of motion for the fields in the coupled system from a general hamiltonian formalism louirad . The hamiltonian of the electromagnetic field inside a resonator with perfectly conducting walls can be written, in terms of the fields amplitudes:

$$H=\frac{1}{2}_V\left(ϵ_0\stackrel{}{E}\stackrel{}{E}+\mu _0\stackrel{}{H}\stackrel{}{H}\right)𝑑V$$
(46)
If we substitute eqs. (7)–(8) in eq. (46) and use the orthonormality condition, the hamiltonian for the field becomes:

$$H=\frac{1}{2}\underset{n}{}\left(_n^2+_n^2\right)$$
(47)
If we define a generalized coordinate $`X_n`$ as:

$$X_n=\frac{_n}{\omega _n}$$
(48)
and its conjugate momentum as:

$$P_n=_n$$
(49)
the hamiltonian can be written as:

$$H=\frac{1}{2}\underset{n}{}\left(P_n^2+\omega _n^2X_n^2\right)$$
(50)
which is identical to the hamiltonian of an infinite set of uncoupled harmonic oscillators. It can be easily verified that the field equations of motion can be derived from this hamiltonian by:

$`{\displaystyle \frac{H}{X_n}}`$ $`=`$ $`\dot{P}_n=\omega _n^2X_n`$
$`{\displaystyle \frac{H}{P_n}}`$ $`=`$ $`\dot{X}_n=P_n`$ (51)
which, in terms of the fields become:

$`{\displaystyle \frac{d_n}{dt}}`$ $`=`$ $`\omega _n_n`$
$`{\displaystyle \frac{d_n}{dt}}`$ $`=`$ $`\omega _n_n`$ (52)
which are identical to eqs. (15)–(16) if wall dissipation is neglected.

### V.1 Field equations for the perturbed system

In order to find the equations of motions for the fields in the perturbed resonator we shall make use of the results obtained in the previous sections. Let us consider again an external time–dependent perturbation, whose rate of change is much less than the rate of change of the fields. The normal modes in the perturbed resonator are given by eqs. (II.1)–(31). If we are looking for the fields in a frequency region where only two electromagnetic modes give a significant contribution, the hamiltonian of the system, as a function of the perturbed modes amplitudes (see eqs. (22) and (23), will be:

$$H=\frac{1}{2}\left(\stackrel{~}{}_1^2+\stackrel{~}{}_1^2+\stackrel{~}{}_2^2+\stackrel{~}{}_2^2\right)$$
(53)
If we now substitute in the above expression the expansions (II.1) we obtain the hamiltonian written in terms of the unperturbed modes amplitudes:

$`H={\displaystyle \frac{1}{2}}\left(_1^2+_1^2+_2^2+_2^2+{\displaystyle \frac{p_m^2}{M}}+M\omega _m^2q_m^2\right)+`$
$`{\displaystyle \frac{q_m}{2}}\left(C_{11}^m_1^2+C_{22}^m_2^2+2C_{12}^m_1_2\right)q_mf_m`$ (54)
where the hamiltonian of a mechanical harmonic oscillator, coupled to an external, time–dependent force, has been included.

From this hamiltonian the equations of motion for the fields can readily be obtained:

$$\frac{d_1}{dt}\omega _1_1=\frac{\omega _1}{𝒬_1}_1$$
(55)
$$\frac{d_2}{dt}\omega _2_2=\frac{\omega _2}{𝒬_2}_2$$
(56)
$$\frac{d_1}{dt}+\omega _1_1=\omega _1q_m\left(C_{11}^m_1+C_{12}^m_2\right)$$
(57)
$$\frac{d_2}{dt}+\omega _2_2=\omega _2q_m\left(C_{12}^m_1+C_{22}^m_2\right)$$
(58)
$$\frac{dq_m}{dt}\frac{p_m}{M}=\frac{\omega _m}{Q_m}q_m$$
(59)
$$\frac{dp_m}{dt}+M\omega _m^2q_m=f_mf_m^{ba}$$
(60)
where the dissipative terms have been added by hand and where the term $`f_m^{ba}`$, which describes the back–action effect of the fields on the walls is given by:

$$f_m^{ba}=\frac{1}{2}\left(C_{11}^m_1^2+C_{22}^m_2^2\right)+C_{12}_1_2$$
(61)
## VI Coupling coefficient calculation

The explicit calculation of the coupling coefficient $`C_{21}^m`$ is not trivial for an arbitrary deformation of the resonator volume and in general can be done only by numerical methods.

First let us note a general property of the coupling coefficients. Let us consider the resonant modes of the two coupled cavities. As previously noted we have the symmetric and the antisymmetric mode; for the former the electric field $`\stackrel{}{E}`$ and the magnetic field $`\stackrel{}{H}`$ in the first cavity are equal to the electric and magnetic field in the second cavity, while for the latter $`\stackrel{}{E}`$ and $`\stackrel{}{H}`$ in the first cavity are equal respectively to $`\stackrel{}{E}`$ and $`\stackrel{}{H}`$ in the second one. From this it follows that in the definition of $`C_{21}`$, the integrand expression $`(\stackrel{}{H}_2\stackrel{}{H}_1\stackrel{}{E}_2\stackrel{}{E}_1)`$ – where we remind that the subscript $`1`$ indicates the symmetric mode and $`2`$ the antisymmetric mode – is odd over the whole detector volume. For this reason if the volume perturbation, over which we perform the integration, is symmetric between the two cavities, the coupling coefficient vanishes, because the contributions to the integral coming from the two cavities, cancel each other. Otherwise, if the volume perturbation is antisymmetric (when one cavity shrinks, the other expands) the two contributions are added with the same sign, and the coupling coefficient is maximum. This general property suggests that we must find a geometrical configuration of our detector such that the volume deformation due to a g.w. is antisymmetric for the two cavities. This was pointed out already in previous works, where the argument was based on the fact that, since the g.w. carries an angular momentum equal to $`2`$, the angular momenta of the fields of the two modes should differ by $`2`$. This can be achieved by putting the two cavities at right angle or by a suitable polarization of the electromagnetic field inside the resonator.

These concepts were verified by both analytical and numerical calculations. General arguments suggested that for an ideal spherical hollow resonator, excited in the fundamental quadrupolar mechanical mode and in the TE011 electromagnetic mode, we should have666We remind that for a TE e.m. mode we have vanishing electric field on the resonator’s surface; for this reason the electric field plays no role in the coupling coefficient calculation.:

$`C_{21}^m=C_{12}^m`$ $`=`$ $`{\displaystyle _S}\left(\stackrel{}{H}_1\stackrel{}{H}_2\right)\stackrel{}{\xi }_m𝑑\stackrel{}{S}=0.4`$
$`C_{11}^m`$ $`=`$ $`{\displaystyle _S}\left(\stackrel{}{H}_1\stackrel{}{H}_1\right)\stackrel{}{\xi }_m𝑑\stackrel{}{S}=0`$ (62)
$`C_{22}^m`$ $`=`$ $`{\displaystyle _S}\left(\stackrel{}{H}_2\stackrel{}{H}_2\right)\stackrel{}{\xi }_m𝑑\stackrel{}{S}=0`$
More detailed calculations, made on a realistic model of the coupled spheres, including the central coupling cell and the e.m. input and output ports, were made by finite element methods. These calculations showed that $`C_{11}^m=C_{11}^m10^4`$, while $`C_{12}^m=C_{21}^m0.2`$.

## VII Calculation of the detector’s signal

As already pointed out in the introduction, the proposed detector exploits the energy transfer induced by the gravitational wave between two levels of an electromagnetic resonator, whose frequencies $`\omega _1`$ and $`\omega _2`$ are both much larger than the angular frequency $`\mathrm{\Omega }`$ of the g.w. and satisfy the resonance conditon $`\omega _2\omega _1\mathrm{\Omega }`$. This is an example of a frequency converter, i.e. a nonlinear device in which energy is transferred from a reference frequency to a different frequency by an external pump signal. This can be viewed as a three-bodies interaction (given by the field–wall interaction term in the hamiltonian (V.1)) which corresponds to annihilation of quanta at $`\omega _1`$ and $`\omega _m`$ and creation at $`\omega _2`$ (or vice versa). For this reason we could argue that, since for a small perturbation $`_1(t)`$ and $`_1(t)`$ will approximately be sinusoidal functions at frequency $`\omega _1`$, while $`_2(t)`$ and $`_2(t)`$ will oscillate at frequency $`\omega _2`$, in the hamiltonian only the terms varying as $`\omega _2\omega _1`$ will give a significant contribution to the interaction. The d–c terms will just give an average deformation of the detector’s walls, determining a static frequency shift of the resonant modes, while the rapidly fluctuating terms at $`\omega _2+\omega _1`$ would practically average to zero.

We shall now calculate the field that is excited by the boundary pertubation in mode $`2`$, starting from an initial condition with mode $`1`$ strongly excited in the resonator. To simplify the analysis of the system of differential equations (55)–(60) we will neglect the small perturbation, due to the external force, on the initially excited e.m. mode (mode 1), and will set $`_1A_1\mathrm{cos}(\omega _1t)`$ and $`_1A_1\mathrm{sin}(\omega _1t)`$, with constant amplitude $`A_1`$.777Actually the amplitude of mode $`1`$ is kept constant by an external rf power source. Furthermore we shall consider the coupling between two TE modes of a resonator: for these modes we have vanishing electric field on the resonator surface. Switching to the complex notation,888In the following it is understood that the physical fields are the real parts of the complex quantities. we obtain:

$$f_m^{ba}\frac{1}{2}C_{21}^m\mathrm{}(_2_1^{})$$
(63)
Finally, to further simplify our calculations, we shall choose a resonator geometry and e.m. field distribution so that $`C_{11}^m=C_{22}^m=0`$. We shall see in a following section that this choice is always possible.

With this assumptions, and taking $`𝒬_1𝒬_2𝒬`$, we can recast the coupled system of equations in the following form:

$$\ddot{}_2+\frac{\omega _2}{𝒬}\dot{}_2+\omega _2^2_2=\omega _2^2q_mC_{21}^m_1$$
(64)
$$\ddot{q}_m+\frac{\omega _m}{Q_m}\dot{q}_m+\omega _m^2q_m=\frac{f_m}{M}\frac{1}{2}\frac{C_{21}^m}{M}_2_1^{}$$
(65)
We apply the following substitutions:

$`_2(t)`$ $`=`$ $`A_2(t)\mathrm{exp}(i\omega _2t)`$
$`q_m(t)`$ $`=`$ $`Q(t)\mathrm{exp}(i\mathrm{\Omega }t)`$ (66)
$`f_m(t)`$ $`=`$ $`F(t)\mathrm{exp}(i\mathrm{\Omega }t)`$
Eqs. (64)–(65) now become:

$`\ddot{A}_2+a\dot{A}_2+bA_2`$ $`=`$ $`cQ`$
$`\ddot{Q}+d\dot{Q}+eQ`$ $`=`$ $`gA_2+F/M`$ (67)
where $`a`$, $`b`$, $`c`$, $`d`$, $`e`$ and $`g`$ are constant coefficients defined by

$`a`$ $`=`$ $`2i\left(\omega _1+\mathrm{\Omega }\right)+{\displaystyle \frac{\omega _2}{𝒬}}`$
$`b`$ $`=`$ $`\omega _2^2\left(\omega _1+\mathrm{\Omega }\right)^2+i{\displaystyle \frac{\omega _2}{𝒬}}\left(\omega _1+\mathrm{\Omega }\right)`$
$`c`$ $`=`$ $`\omega _2^2C_{21}^mA_1`$
$`d`$ $`=`$ $`2i\mathrm{\Omega }+{\displaystyle \frac{\omega _m}{Q_m}}`$ (68)
$`e`$ $`=`$ $`\omega _m^2\mathrm{\Omega }^2+i{\displaystyle \frac{\omega _m\mathrm{\Omega }}{Q_m}}`$
$`g`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{C_{21}^mA_1}{M}}`$
We can now fourier transform eqs. (VII) and solve them for $`A_2(\omega )`$. We find:

$$A_2(\omega )=2\pi \frac{cF(\omega )/M}{\left(b\omega ^2+ia\omega \right)\left(e\omega ^2+id\omega \right)gc}$$
(70)
$`_2(t)`$ is then given by:

$$_2(t)=\mathrm{exp}(i\omega _2t)_{\mathrm{}}^{\mathrm{}}\frac{cF(\omega )/M\mathrm{exp}(i\omega t)}{\left(b\omega ^2+ia\omega \right)\left(e\omega ^2+id\omega \right)gc}𝑑\omega $$
(71)
For a plane g.w travelling along the $`z`$ axis the force density, in the proper reference frame attached to the detector, has the form:

$$\stackrel{}{f}(\stackrel{}{x},t)=\frac{1}{2}\rho (\stackrel{}{x})[(\ddot{a}_{\mathrm{\hspace{0.17em}1}}^1x+\ddot{a}_{\mathrm{\hspace{0.17em}2}}^1y),(\ddot{a}_{\mathrm{\hspace{0.17em}2}}^1x\ddot{a}_{\mathrm{\hspace{0.17em}1}}^1y),(0)]$$
(72)
where $`a_j^i(t)`$, is the adimensional amplitude of the wave, and $`a_{\mathrm{\hspace{0.17em}1}}^1=a_{\mathrm{\hspace{0.17em}2}}^2`$, $`a_{\mathrm{\hspace{0.17em}2}}^1=a_{\mathrm{\hspace{0.17em}1}}^2`$. The generalized force, acting on the $`m`$ mechanical mode, then has the form

$`f_m={\displaystyle \frac{1}{2}}\ddot{a}_{\mathrm{\hspace{0.17em}1}}^1{\displaystyle _{Vol}}\left((\xi _m)_xx(\xi _m)_yy\right)\rho (\stackrel{}{x})𝑑V`$
$`{\displaystyle \frac{1}{2}}\ddot{a}_{\mathrm{\hspace{0.17em}2}}^1{\displaystyle _{Vol}}\left((\xi _m)_xy+(\xi _m)_yx\right)\rho (\stackrel{}{x})𝑑V`$ (73)
If $`a_j^i(t)`$ is given by:

$$a_j^i(t)=\left(\begin{array}{cc}h_{\mathrm{\hspace{0.17em}1}}^1\alpha (t)& h_{\mathrm{\hspace{0.17em}2}}^1\beta (t)\\ h_{\mathrm{\hspace{0.17em}2}}^1\beta (t)& h_{\mathrm{\hspace{0.17em}1}}^1\alpha (t)\end{array}\right)$$
(74)
where $`\alpha (t)`$ and $`\beta (t)`$ are sinusoidal functions of frequency $`\mathrm{\Omega }/(2\pi )`$, then $`\ddot{a}_j^i=\mathrm{\Omega }^2a_j^i`$, and if we define the effective lengths of our detector as

$`L_+`$ $`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle _{Vol}}\left((\xi _m)_xx(\xi _m)_yy\right)\rho (\stackrel{}{x})𝑑V`$
$`L_\times `$ $`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle _{Vol}}\left((\xi _m)_xy+(\xi _m)_yx\right)\rho (\stackrel{}{x})𝑑V`$ (75)
we can write

$$f_m(t)=\frac{1}{2}M\mathrm{\Omega }^2\left(L_+a_{\mathrm{\hspace{0.17em}1}}^1+L_\times a_{\mathrm{\hspace{0.17em}2}}^1\right)$$
(76)
and

$$A_2(\omega )=\frac{1}{2}\frac{c\mathrm{\Omega }^2\left(L_+h_{\mathrm{\hspace{0.17em}1}}^1+L_\times h_{\mathrm{\hspace{0.17em}2}}^1\right)}{begc}\delta (\omega \mathrm{\Omega })$$
(77)
or, making use of eq. (VII)

$$_2(t)=\frac{1}{2}\frac{\omega _2^2C_{21}^mA_1\left(L_+h_{\mathrm{\hspace{0.17em}1}}^1+L_\times h_{\mathrm{\hspace{0.17em}2}}^1\right)\mathrm{\Omega }^2\mathrm{exp}(i\omega _2t)}{\left(\omega _m^2\mathrm{\Omega }^2+i\frac{\omega _m\mathrm{\Omega }}{Q_m}\right)\left(\omega _2^2(\omega _1+\mathrm{\Omega })^2+i\frac{\omega _2(\omega _1+\mathrm{\Omega })}{𝒬}\right)\frac{(\omega _2C_{21}^mA_1)^2}{2M}}$$
(78)
The electric field amplitude the initially empty mode can be readily obtained from eq. (56); the average energy stored in mode number 2 is given by: $`U_2=(1/2)|_2|^2=(1/2)|_2|^2`$.

## VIII Noise issues

### VIII.1 Mechanical thermal noise

Thermal noise is one of the fundamental limits in the measurement of small displacements. In particular it is one of the dominant noise sources in resonant–mass detectors of g.w. and a major reason that such detectors operate at cryogenic temperatures. Since our detector exploits the coupling of the g.w. with the mechanical structure of the resonator, we have to carefully study the thermal noise contribution to our output signal.

We start again from eqs. (64)–(65) taking now the external force $`f_m(t)`$ as a stochastic force with constant power spectrum $`S_{ff}`$, given by papoulis :

$$S_{ff}=\frac{4Mk_BT\omega _m}{Q_m}$$
(79)
Making the substitutions:

$`_2(t)`$ $`=`$ $`A_2(t)\mathrm{exp}(i\omega _2t)`$
$`q_m(t)`$ $`=`$ $`Q(t)\mathrm{exp}(i(\omega _2\omega _1)t)`$ (80)
we obtain the following equations:

$`\ddot{A}_2+a\dot{A}_2+bA_2`$ $`=`$ $`cQ`$
$`\ddot{Q}+d\dot{Q}+eQ`$ $`=`$ $`gA_2+{\displaystyle \frac{f(t)}{M}}\mathrm{exp}(i(\omega _1\omega _2)t)`$ (81)
where $`a`$, $`b`$, $`c`$, $`d`$, $`e`$ and $`g`$ are defined as in eq. (VII) with the parameter $`\mathrm{\Omega }`$ replaced by the difference $`\omega _2\omega _1`$.

The first equation in (VIII.1) can be solved for $`Q(t)`$, and we are left with one equation for the variable $`A_2(t)`$:

$`{\displaystyle \frac{d^4A_2}{dt^4}}(d+a){\displaystyle \frac{d^3A_2}{dt^3}}(b+e+ad){\displaystyle \frac{d^2A_2}{dt^2}}(db+ea){\displaystyle \frac{dA_2}{dt}}+(cgeb)A_2=`$
$`c{\displaystyle \frac{f(t)}{M}}\mathrm{exp}(i(\omega _1\omega _2)t)`$ (82)
For a linear system we can immediately write the spectral density of the amplitude $`A_2`$ as papoulis :

$$S_{AA}(\omega )=|p(\omega )|^2\left(\frac{c}{M}\right)^2S_{ff}$$
(83)
with

$$p(\omega )=\omega ^4+i(a+d)\omega ^3+(e+b+da)\omega ^2i(ea+db)\omega +cgeb$$
(84)
being the fourier transform of the linear system’s impulse response function.

From the above equations we readily find the spectral densities of $`_2(t)`$ as:

$$S_{HH}(\omega )=|p(\omega \omega _2)|^2\left(\omega _2C_{21}^mA_1\right)^2\frac{4k_BT\omega _m}{MQ_m}$$
(85)
### VIII.2 Other noise sources

#### VIII.2.1 Master oscillator phase noise

To operate our device we have to feed microwave power into one resonant mode (say mode $`1`$), in order to detect the energy transfer between the full and the initially empty mode driven by the external perturbation.

To feed power into our device we shall use a voltage controlled microwave oscillator locked on mode $`1`$, at frequency $`\omega _1`$. The master oscillator phase noise is filtered through the resonator linewidth; the power spectral density has the following frequency dependence slater :

$$S_{LO}(\omega )=\frac{4\beta P_I/(\omega _1𝒬)}{(1/𝒬)^2+(\omega /\omega _1\omega _1/\omega )^2}$$
(86)
where $`P_I`$ is the power input level and $`\beta `$ is the coupling coefficient of mode $`1`$ to the output load. From the above equation we can estimate the microwave power noise spectral density at the detection frequency $`\omega _2`$:

$$S_{LO}(\omega _2)=\frac{4\beta P_I/(\omega _1𝒬)}{(1/𝒬)^2+(\omega _2/\omega _1\omega _1/\omega _2)^2}\beta \frac{P_I}{\omega _1𝒬}\left(\frac{\omega _2}{\omega _2\omega _1}\right)^2$$
(87)
This figure can be improved if the receiver discriminates the parity of the e.m. field at frequency $`\omega _2`$, i.e. if it is sensitive only to the power excited in mode number $`2`$, rejecting all contributions coming from mode number $`1`$. In this way mode $`1`$ becomes decoupled from the output load and $`\beta =0`$. The experimental set–up, based on the use of two magic–tees which accomplishes this issue is discussed in detail in rsi . Of course the mode discrimination cannot be ideal, and some power leaking from mode $`1`$ to the detector’s output will be present. Nevertheless our previous work has demonstrated that with a careful tuning of the detection electronics we can obtain $`\beta 10^{14}`$ rsi .

#### VIII.2.2 Amplifier noise

The input Johnson noise of the first amplifier in the detection electronics has to be added to the previous contributions to establish the overall noise spectral density. It can be described by the frequency independent spectral density papoulis :

$$S_{JJ}=k_BT\times 10^{(N/10)}k_BT_{eq}$$
(88)
where $`N`$ is the noise figure of the amplifier (in dB) and $`T`$ the operating temperature.

## IX Detector sensitivity

The detector sensitivity is ultimately determined by the overall effect of the various noise sources discussed in section VIII and, eventually, by several others. Of course, depending on the characteristics of the system and on the experimental set–up, different noise surces will become dominant.

We shall characterize the noise in our detector by a frequency dependent spectral density $`S_n(f)`$, with dimension Hz-1, defined as follows thorne : if a sinusoidal g.w. with known phase $`\varphi `$, known frequency $`f`$ and unknown r.m.s. amplitude $`\sqrt{2}h_0`$, impinges on the detector, and if we try to detect the wave by fourier analyzing the detector output with a bandwidth $`\mathrm{\Delta }f`$, then the amplitude signal–to–noise ratio will be:

$$\frac{S}{N}=\frac{h_0}{(S_n(f)\mathrm{\Delta }f)^{1/2}}$$
(89)
We shall also define the minimum detectable wave amplitude (at 90% C.L.) for a periodic source with known frequency and phase as:

$$h_{min}(f)=1.7\left(\frac{S_n(f)}{F}\right)^{1/2}$$
(90)
with dimension Hz-1/2. In the following calculations the average pattern function value $`F=2/5`$ has been taken maggiore .

Let us focus our attention on the system mentioned in section IV based on two spherical niobium cavities working at $`\omega _1\omega _22`$ GHz with a stored energy in the initially excited symmetric mode of $`U_110`$ J per cell. This is a small–scale system with an effective length of 0.1 m and a typical weigth of 5 kg. The lowest quadrupolar mechanical mode is at $`\omega _m4`$ kHz. In the following we shall consider an equivalent temperature of the detection electronics $`T_{eq}=30`$ K.

A possible design of the detector uses both the mechanical resonance of the structure, and the e.m. resonance. This can be accomplished if the detector is designed in order to have the mechanical mode frequency equal to the e.m. modes frequency difference $`\omega _m\omega _2\omega _1`$. In this frequency range, with reasonable values for the system parameters, the dominant noise source will be the noise coming from the brownian motion of the detector walls. The expected sensitivity of the detector for $`\omega _2\omega _1=\omega _m=4`$ kHz is shown in figure 4. In figure 5 the separate contribution of the noise sources – discussed in section VIII – to the overall noise spectral density is shown. We point out that even if in this case the dominant noise source is the walls thermal motion a lower $`T_{eq}`$ would increase the detection bandwith, as shown in figure 6.

Since our detector is based on a double resonant system (the mechanical resonator and the electromagnetic resonator) it can be operated also for frequencies $`\omega _2\omega _1\omega _m`$. At frequencies $`\omega _2\omega _11`$ kHz the master oscillator phase noise will, in general, be dominant (see sec. VIII.2.1), while at frequencies $`\omega _2\omega _110`$ kHz the noise coming from the detection electronics will dominate (at least for $`T_{eq}30`$ K), as shown in figure 7. The expected sensitivity of the detector for $`\omega _2\omega _1=10`$ kHz is shown in figure 8.

In order to work at frequencies $`\omega _2\omega _11`$ kHz a large–scale system has to be developed. A possible design could be based on two spherical cavities working at $`\omega _1\omega _2500`$ MHz, with $`\omega _2\omega _11`$ kHz. This system could have a stored energy of $`U_1800`$ J per cell, an effective length of 0.4 m and a typical weigth of 300 kg. With a rather optimistic (but not unrealistic) choice of system parameters one could obtain the sensitivity shown in figure 9. We point out that in this figure an electronics equivalent temperature of $`T_{eq}1`$ K has been used; also in this case lowering $`T_{eq}`$ corresponds to an increase of the detection bandwidth (see fig. 11).

The large–scale system could also be used at higher frequencies; in this case a good sensitivity can be achieved in a narrow detection bandwidth (see fig. 12).

## X Conclusions

A first prototype of the detector has been built and successfully tested rsi . A detector based on two coupled spherical cavities has been designed and preliminar mechanical and electromagnetic tests are being made on normal conducting prototypes. The planned timeline is as follows:

* In 2002 a bulk niobium detector (coupled spherical cavities, $`\omega =2`$ GHz, $`\omega _2\omega _1=10`$ kHz, fixed coupling) will be built at CERN;
* In 2003 a variable coupling detector will be built and tested.

In the meantime several open problems must be addressed:

* The mechanical quality factor of the detector has to be maximized in order to suppress the noise coming from the brownian motion of the detector walls. Since mechanical dissipations arise from materials intrinsic losses and from the coupling of the system to the external environment, materials with low intrinsic losses must be used for the construction of the detector and the design of a suitable suspension system has to be done carefully.
* The requirement of an high mechanical quality factor has to be matched with the requirement of high electromagnetic quality factor. This can be accomplished by the use of bulk niobium, which, at low temperatures, has low intrinsic losses both mechanical and electromagnetic, or by the use of a niobium thin film deposited on a high mechanical quality factor substrate. Both tecniques present in principle advantages and drawbacks. Several prototypes of single–cell, seamles, copper spherical cavities have been built at INFN–LNL by E. Palmieri and will be sputter–coated and tested at CERN to check the quality of noibium films deposited on spherical substrates.
* A cryogenic system with a cooling power of $`5`$ W at $`T1.8`$ K and $`P=1`$ bar has to be designed. The contribution of the cryogenic system to the noise has to be studied carefully.
* The readout electronics has to be optimized. The use of a low noise transducer, possibly based on the SQUID technology, has to be investigated.

If experimental results will be encouraging, by the end of 2003 a proposal for the construction of a g.w detector, based on superconducting rf cavities could be considered.


text:
# Excitations in time-dependent density-functional theory

Also at \]Department of Physics, Rutgers University, Frelinghuysen Road, Piscataway, NJ 08854

## Abstract

An approximate solution to the time-dependent density functional theory (TDDFT) response equations for finite systems is developed, yielding corrections to the single-pole approximation. These explain why allowed Kohn-Sham transition frequencies and oscillator strengths are usually good approximations to the true values, and why sometimes they are not. The approximation yields simple expressions for Görling-Levy perturbation theory results, and a method for estimating expectation values of the unknown exchange-correlation kernel.

Traditional density-functional theory (DFT) is a popular and efficient method for the calculation of ground-state properties of interacting many electron systems. While it has long been the method of choice for solid-state calculations JG89 , recent improvements in approximate functionals have also made it popular in quantum chemistry, because of the ability to handle large systems Kb99 . Within exact ground-state theory, the Kohn-Sham eigenvalues and orbitals are mathematical constructs, designed to reproduce the ground-state density. The one exception is that the highest occupied eigenvalue is known to equal minus the ionization potential of the system. For many years, the unoccupied Kohn-Sham levels of solids have been used to interpret the excitations in solids, despite the well-known underestimate of the band gap of insulators in LDA.

There exist several extensions of the basic formalism that allow extraction of excited-state properties. One popular approach is via time-dependent DFT (TDDFT), in which the interacting system in a time-dependent external field is mapped exactly to a non-interacting time-dependent Kohn-Sham (KS) system with the same time-dependent density RG84 . If a weak time-dependent electrical field is considered, this leads to a Dyson-like response equation for the exact susceptibility of the interacting electronic system. Poles of this susceptibility occur at the true transition frequencies PGG96 , and the strengths of these poles are related to oscillator strengths. Solution of these equations has been implemented in several quantum chemical packages, and excitation spectra of many molecules have been calculated and reported in hundreds of papers (See references in Ref. MBAG01 ). Once an accurate ground-state Kohn-Sham potential is used, transition frequencies are typically within about 0.2 eV of experiment. Oscillator strengths are usually good in these calculations (within a factor of 2), but not always. There are considerable subtleties in applying TDDFT to extract the optical response of solids 01BKBL .

The beauty of using TDDFT to extract optical spectra is that it combines moderate accuracy with inexpensive calculation, just as in the ground-state case. But even beginning from an exact ground-state KS potential, the spectrum still depends on the unknown exchange-correlation (XC) kernel, i.e., the functional derivative of the time-dependent XC potential. This is often approximated by the crude but reliable adiabatic local-density or generalized-gradient approximations. These TDDFT calculations compare very favorably with configuration-interaction singles (CIS) calculations, the only alternative that is comparable in computational cost HH99 . Higher level calculations, such as more complete CI, Bethe-Salpeter, or quantum Monte Carlo GRML01 , can be made more accurate, but cost more, limiting their use to smaller systems. For example, the TDDFT approach has recently been applied to electron-transfer problems in biological systems 01HFHH .

Although TDDFT methodology has been implemented and is being used widely, understanding of its accuracy and reliability has been slow, as well as its relation to other methods. The relation to first-order Görling-Levy perturbation theory has been shown GS99 ; BPG00 , as well as the connection with the GW approximation GS99 . The extreme case of stretched H2 has recently been studied by several authors 00CGGG ; GGBb00 ; AGR02 , although this also represents difficulties for the ground-state theory PSB95 . By using a matrix formulation of Casida C96 , the present paper shows how, when the excitations of a system are discrete, an approximate solution, that can be made arbitrarily accurate, can be used to understand and explain many trends in the results of TDDFT calculations.

To demonstrate our results, we apply them to the prototype systems of the He and Be atoms. We chose these because their exact ground-state KS potentials are known UG94 , and because their lowest allowed transitions exhibit two classes of systems we are interested in. In the He atom the 1s to 2p singlet transition is at 21.22 eV, while the KS transition, i.e., the energy difference between the 1s and 2p ground-state KS orbital energy levels, is 21.15 eV, less than 0.1 eV smaller. On the other hand, in the Be atom, the 2s to 2p level singlet transition is at 5.3 eV, but in the KS case it is at 3.6 eV. We explain below the fundamental difference between these two systems, and why the KS eigenvalues are a good approximation in the first case, but a poor approximation in the second. Furthermore, the oscillator strength for the 1s to 2p transition in the He atom is 0.27, but 0.32 in the exact KS case. Thus the oscillator strength of the He atom KS system is close to the true one, but noticeably less close than the transition frequency. Finally, the 2s to 2p oscillator strength in the Be atom is 1.37, but 2.54 in its Kohn-Sham alter ego, far closer than one would expect, given the error in transition frequencies. We explain these facts.

We denote the exact KS eigenvalues as $`ϵ_i`$ and orbitals as $`\varphi _i(𝐫)`$. Casida C96  has written the TDDFT response equations as an eigenvalue equation for the square of the transition frequencies:

$$\underset{q^{}}{}\stackrel{~}{\mathrm{\Omega }}_{qq^{}}(\omega )v_q^{}=\mathrm{\Omega }v_q,$$
(1)
where $`q`$ is a double index, representing a transition from occupied KS orbital $`i`$ to unoccupied KS orbital $`a`$, $`\omega _q=ϵ_aϵ_i`$, $`\mathrm{\Omega }=\omega ^2`$, and $`\mathrm{\Phi }_q(𝐫)=\varphi _i^{}(𝐫)\varphi _a(𝐫)`$. The matrix is

$$\stackrel{~}{\mathrm{\Omega }}_{qq^{}}=\delta _{qq^{}}\mathrm{\Omega }_q+2\sqrt{\omega _q\omega _q^{}}q|f_{\mathrm{HXC}}(\omega )|q^{},$$
(2)
where

$$q|f_{\mathrm{HXC}}(\omega )|q^{}=d^3rd^3r^{}\mathrm{\Phi }_q^{}(𝐫)f_{\mathrm{HXC}}(𝐫,𝐫^{},\omega )\mathrm{\Phi }_q^{}(𝐫^{}).$$
(3)
In this equation, $`f_{\mathrm{HXC}}`$ is the Hartree-exchange-correlation kernel, $`1/|𝐫𝐫^{}|+f_{\mathrm{XC}}(𝐫,𝐫^{},\omega )`$, where $`f_{\mathrm{XC}}`$ is the unknown XC kernel.

It has been noticed that the KS transition frequencies are often ‘good’ approximations to the true frequencies ARU98 . If the transition frequencies and oscillator strengths are expanded in powers of the coupling constant $`\lambda `$, as in GL perturbation theory GL92 , the zero-order values are the Kohn-Sham values G96 . In Eq. (3), we see that, if $`f_{\mathrm{HXC}}`$ is small, i.e., if the system is sufficiently weakly correlated, this is correct. We show below that this is not the reason why the KS values are good approximations in many systems, such as the He atom, and is especially untrue when the ground-state has a near-degeneracy, such as for the Be atom.

Our approximate solution relies on the fact that $`q|f_{\mathrm{HXC}}|q^{}`$ decays rapidly with distance from the diagonal, because the overlap of increasingly different orbitals decay by cancellation of oscillations. To zero order, we ignore all off-diagonal elements, findingVOC99  the small-matrix approximation (SMA)GKG97 :

$$\mathrm{\Omega }=\mathrm{\Omega }_q+2\omega _qq|f_{\mathrm{HXC}}(\omega _q)|q.$$
(4)
The original single-pole approximation (SPA) PGG96  can be viewed as a special case of SMA when the shift from the KS value is small:

$$\omega =\sqrt{\mathrm{\Omega }_q+2\omega _qq|f_{\mathrm{HXC}}|q}=\omega _q+q|f_{\mathrm{HXC}}(\omega _q)|q+\mathrm{}$$
(5)
In the special case of including just the exchange kernel, this result is identical to first-order GL perturbation theory GS99 .

To go beyond SPA, we use a continued-fraction method S82  for inverting a matrix with a dominant diagonal. Truncating the expansion at second order in the off-diagonal matrix elements:

$$\mathrm{\Omega }=\mathrm{\Omega }_q^{\mathrm{SMA}}+\underset{q^{}q}{}\frac{4\omega _q\omega _q^{}|q|f_{\mathrm{HXC}}(\omega _q^{\mathrm{SMA}})|q^{}|^2}{\mathrm{\Omega }_q^{\mathrm{SMA}}\mathrm{\Omega }_q^{}^{\mathrm{SMA}}}.$$
(6)
This is a key result of this paper, leading to many conclusions. First, it yields the exact GL perturbative expression to second-order in $`\lambda `$. Expanding $`f_{\mathrm{HXC}}=\lambda f_{\mathrm{HX}}+\lambda ^2f_\mathrm{C}^{[2]}+\mathrm{}`$, we find

$$\omega =\omega _q+\lambda q|f_{\mathrm{HX}}(\omega _q)|q+\lambda ^2\delta \omega _q^{(2)},$$
(7)
where the second-order shift consists of four terms:

$`\delta \omega _q^{(2)}`$ $`=`$ $`q|f_\mathrm{C}^{[2]}(\omega _q)|q+2{\displaystyle \underset{q^{}q}{}}{\displaystyle \frac{\omega _q^{}|q|f_{\mathrm{HX}}(\omega _q)|q^{}|^2}{\mathrm{\Omega }_q\mathrm{\Omega }_q^{}}}`$
$`+`$ $`q|f_{\mathrm{HX}}(\omega _q)|qq|{\displaystyle \frac{f_{\mathrm{HX}}}{\omega }}(\omega _q)|q{\displaystyle \frac{|q|f_{\mathrm{HX}}(\omega _q)|q|^2}{2\omega _q}}.`$
(8)
For the ground-state energy, the second-order correction has been identified as playing a key role in constructing accurate functionals, especially in cases of strong static correlation SPK00 . It is likely to play a similar role for excitations, and can be easily extracted from Eq. (7).

Second, we may now deduce precisely when SPA (or SMA) is valid. Defining the shift from the KS value as

$$\mathrm{\Delta }\mathrm{\Omega }_q=\mathrm{\Omega }\mathrm{\Omega }_q,$$
(9)
we rewrite Eq. (6) in the following suggestive form:

$$\mathrm{\Delta }\mathrm{\Omega }_q=\mathrm{\Delta }\mathrm{\Omega }_q^{\mathrm{SMA}}\left\{1+\underset{q^{}q}{}\frac{\mathrm{\Delta }\mathrm{\Omega }_q^{}^{\mathrm{SMA}}}{\mathrm{\Omega }_q^{\mathrm{SMA}}\mathrm{\Omega }_q^{}^{\mathrm{SMA}}}\frac{|M_{qq^{}}|^2}{M_{qq}M_{q^{}q^{}}}\right\},$$
(10)
where $`M_{qq^{}}=q|f_{\mathrm{HXC}}|q^{}`$. A simple estimate of the size of this correction can be given by assuming $`M_{qq^{}}\sqrt{M_{qq}M_{q^{}q^{}}}`$. This would be exactly true if $`f_{\mathrm{HXC}}(𝐫,𝐫^{},\omega )=f(𝐫,\omega )f(𝐫^{},\omega )`$, and is accurate to within 0.5% for exchange in the He atom. Then the SMA is valid when

$$\underset{q^{}q}{}\frac{\mathrm{\Delta }\mathrm{\Omega }_q^{}^{\mathrm{SMA}}}{\mathrm{\Omega }_q^{\mathrm{SMA}}\mathrm{\Omega }_q^{}^{\mathrm{SMA}}}1,$$
(11)
i.e., the SMA shift need only be small on the scale of the separation between transition frequencies. Thus, even when the corrections to KS transition frequencies are large, SMA can remain valid if the poles are well-separated.

In the special case when the SMA correction is small compared to a KS transition itself, this result simplifies to

$$\underset{q^{}q}{}\frac{q|f_{\mathrm{HXC}}|q^{}}{\omega _q\omega _q^{}}1.$$
(12)
For example, the allowed transitions from the ground-state of the He atom are listed in Table 1. Comparing KS transition frequencies with physical ones, we find them good to within less then 0.1 eV. This implies all matrix elements $`q|f_{\mathrm{HXC}}|q^{}`$ are small, and SPA is valid. Calculations within ALDA for this case 00CGGG  show no difference between full solution of the response equations and the SPA result. The column marked SMA lists SMA results with our best estimate of $`f_{\mathrm{HXC}}`$ for this case, a hybrid of exact exchange with ALDA antiparallel correlation BPG00 . For the transition to 2p, the exact exchange result is 21.37 eV PGB00 , showing that there must be substantial cancellation by correlation effects. For this system and others like it, GL perturbation theory converges slowly, while our expansion converges rapidly.

Our other prototype is the excitations from the ground state of the Be atom. For the 2s $``$ 2p transition, the expectation value of $`f_{\mathrm{HXC}}`$ is relatively large. We expect SMA to work quite well for that transition, but less well for others, since the one strong transition contributes to the correction in Eq. (10), especially having a small denominator. This is born out by the frequency results in Table 1. Within ALDA 00CGGG , the SMA transition is at 5.27 eV, but the full calculation is 5.08 eV.

In the SMA, in which off-diagonal elements are neglected, the eigenvalues in Eq. (1) remain unit vectors, and the oscillator strengths retain their KS values. When we include the change due to the off-diagonal elements to leading order, we find:

$$f=\frac{2}{3}\left(\omega _q\stackrel{}{\mu }_q^2+\underset{q^{}q}{}\frac{4q|f_{\mathrm{HXC}}|q^{}\omega _q\omega _q^{}\stackrel{}{\mu }_q\stackrel{}{\mu }_q^{}}{\mathrm{\Omega }_q^{\mathrm{SMA}}\mathrm{\Omega }_q^{}^{\mathrm{SMA}}}\right).$$
(13)
Here $`\stackrel{}{\mu }_q`$ denotes the KS dipole matrix element. Our first conclusion from this important result is that it contains the exact GL expression for oscillator strengths to first order in $`\lambda `$, by ignoring the correlation kernel. The sum is rapidly converging, as the matrix element decays rapidly with principal quantum number. Using $`f_\mathrm{X}`$ for the He atom, the corrections of Eq. (13), summed over only bound-bound transitions, reduce the oscillator strength by 11%, whereas the exact answer is 15% lower than the KS value. The remaining reduction is due to either correlation effects or transitions to the continuum.

If we estimate off-diagonal elements with geometric means of diagonal elements, we find

$$ff_q\left(1+2\underset{q^{}q}{}\frac{\sqrt{\mathrm{\Delta }\mathrm{\Omega }_q^{\mathrm{SMA}}\mathrm{\Delta }\mathrm{\Omega }_q^{}^{\mathrm{SMA}}}}{\mathrm{\Omega }_q^{\mathrm{SMA}}\mathrm{\Omega }_q^{}^{\mathrm{SMA}}}\sqrt{\frac{f_q^{}}{f_q}}\right).$$
(14)
The effect of off-diagonal matrix elements is to mix various KS oscillator strengths. For the dominant transition , if Eq.(11) is satisfied for excitation energies, it is also satisfied for oscillator strengths. The correction to an oscillator strength is first-order in the off-diagonal matrix element, as opposed to the second-order correction to the SMA transition-frequency shift. Thus, fractional corrections to KS oscillator strengths will generally be larger than those to SMA shifts. So even with large SMA shifts, the associated oscillator strengths can be good.

The oscillator strengths of the He and Be atoms confirm our previous conclusions. For the well-separated transitions (He atom), the KS oscillator strengths are close to the true oscillator strengths, but not as close as the transition frequencies. The deviations estimated in Eq. (14) are consistent with those of the transition frequencies. Similarly, in the case of the Be atom, we see that the 2s$``$2p KS oscillator strength is good to within a factor of 2, because the corrections due to other transitions are quite mild. On the other hand, the higher transitions have KS oscillator strengths that are an order of magnitude different from the true ones, because of the huge corrections due to the first transition.

Ours is not an accurate numerical solution of the TDDFT response equations, but is rather a method for understanding results. For example00CGGG , the He oscillator strengths within ALDA are good to within 5%. The present work shows that this reflects the accuracy of $`q|f_{\mathrm{HXC}}^{\mathrm{ALDA}}|q^{}`$ for these transitions. Our approximation scheme handles well-separated poles. As the system-size grows, e.g., even for small clusters, excitations become more dense on the $`\omega `$-axis, and kernel corrections can cause levels to reorder, etc. We believe that these methods can still provide estimates of expected trends, but this remains to be tested. It is unclear what happens in the thermodynamic limit, without better approximations to the XC kernel01BKBL . For finite systems, the Tamm-Dancoff approximation neglects all matrix elements between one forward and one backward transition (reducing the size of matrices by a factor of 2), and has been found to often yield excellent results compared to a full solutionHH99 . However, this approximation violates the Thomas-Reiche-Kuhn sum rule, whereas our expansion satisfies it order-by-order.

We have focussed exclusively on the allowed transitions, because we need exact results on both transition frequencies and oscillator strengths for our analysis. To analyze the limitations of SMA for singlet-triplet splittingsVOC99  would require knowledge of the exact oscillator strengths of spin-flipping transitions.

So far, we have focussed on approximate solutions of the TDDFT response equations for the exact ground-state KS potential. In practice, this potential is approximated. Local-density and generalized-gradient approximations have potentials that are too shallow, so that Rydberg states are not bound. This can be corrected by some addition of the correct asymptotic behavior CS00 , or by use of an orbital-dependent functional GKKG99 ; IHB99 ; G99 , whose derivative yields an accurate potential at large distances. From accurate calculations on the He and Be atoms, we find the principal effects of using either exchange-only or LDA-SIC potentials to be a shift in the orbital energies, numerically identical to the error in the ionization potential. But the KS dipole matrix elements are extremely accurate in these approximate potentials, so that the dominant error in oscillator strengths comes from the errors in eigenvalues. Thus current technology allows accurate calculation of KS oscillator strengths.

We conclude with an observation that should be useful for development of approximations to $`f_{\mathrm{XC}}`$. In cases where poles are well-separated and SMA is valid, as can be determined by comparing oscillator strengths, the difference between KS transition frequencies and physical ones yields an exact expectation value of $`f_{\mathrm{HXC}}`$, to all orders in coupling-constant $`\lambda `$. This would provide an invaluable benchmark for testing approximate XC kernels, similar to the widespread use exact XC potentials have enjoyed in the ground-state caseUG94 .

We thank C. J. Umrigar for providing us with his exact Kohn-Sham potentials for the He and Be atoms. H. Appel acknowledges the support of a fellowship of the Studienstiftung des deutschen Volkes. This work was supported by the National Science Foundation under grant no. CHE-9875091, and also partially by the Petroleum Research Fund.