text:
# Transient QED effects in absorbing dielectrics

## I Introduction

The rate and the spectral and spatial characteristics of the spontaneous decay of an atom depend on the properties of the atom and of the radiation field, and on the interaction between them. The radiation field changes by the presence of other matter . One can try and manipulate the emission properties once the influence of this medium is understood.

In quantum optics of linear dielectrics, one tries to describe the material medium in an effective way with the help of the classical dielectric function $`\epsilon (𝐫,\omega )`$, which in general is a complex function of both position and frequency and in this full generality describes the propagation and loss of light at each point in the dielectric. Sometimes it is possible to neglect the spatial variations (including local field effects), dispersion and losses altogether. The spontaneous emission rate of an atom in such a simple dielectric is the refractive index $`n`$ of the medium times the rate $`\mathrm{\Gamma }_0`$ in vacuum .

The situation becomes more complicated when material dispersion has to be taken into account . Since the Kramers-Kronig relations tell that dispersion and loss always come together (be it not always at the same frequencies), one should like to include losses as well in order to describe all frequencies in one theory. The damped polariton model  provides us with such a microscopic theory. From that theory it was shown that the radiative spontaneous emission rate equals $`\mathrm{\Gamma }_0`$ times the real part of the refractive index at the transition frequency .

The quantum mechanical treatment of dissipative systems is more complicated than the classical one, because of the extra requirement that equal-time commutation relations do not change over time . Based on the damped-polariton model and on the fluctuation-dissipation theorem, phenomenological quantization theories were constructed that meet these requirements. In these theories, the dielectric function is an input function and the Maxwell field operators satisfy quantum Langevin equations with both loss and quantum noise terms . With the use of a Green-function approach, the phenomenological quantization theories have been generalized to inhomogeneous dielectrics, first for multilayer systems and later for general $`\epsilon (𝐫,\omega )`$ . Field commutation relations turn out only to depend on the analytical properties of the Green function. However, the calculation of spontaneous emission inside such a medium would involve the actual computation of the Green function, which for general $`\epsilon (𝐫,\omega )`$ is not easy.

A special case of the former theories is the quantum optical description of inhomogeneous systems at frequencies where both dispersion and losses can be neglected. Then a description in terms of modes is possible, where the mode functions are harmonic solutions of the classical wave equation featuring a position-dependent dielectric “constant” $`\epsilon (𝐫)`$ . This encompasses the now theoretically and experimentally very active research area of the so-called photonic crystals , where a periodic modulation of the refractive index at the scale of the wavelength of light can drastically modify the mode structure compared to vacuum. By increasing the refractive-index contrast, even a photonic bandgap can open up, giving rise to a frequency interval for which waves cannot travel in the crystal in any direction, so that spontaneous emission would be inhibited completely. Until now, such a bandgap has not been found conclusively in the optical regime . It has been proposed to look for frequencies close to material resonances, where refractive indices can be quite substantially higher or lower than 1 .

Interesting new effects have been predicted for bandgap systems, such as photon-atom bound states and non-exponential spontaneous decay at the edges of the gap . A current debate is whether the Weisskopf-Wigner approximation can be used in the calculation of spontaneous emission near an edge of a photonic bandgap. This question seems to depend strongly on the analytic or singular behavior of the density of states at the edges of the gap, which has recently been calculated for face-centered cubic and diamond-like crystal structures . If near the edge of the bandgap a large part of the modes has a cavity-like structure, producing nonzero dwell-times near the emitting atom, then an emitted photon has a nonzero probability of being reabsorbed, which would give Rabi-like oscillations of the atomic population that are missed in the Weisskopf-Wigner approximation.

Non-exponential decay can also be caused by the interference of possible decay-channels: for short times after the excitation of the atom, a larger frequency interval of the medium states plays a part in the decay process than for later times. Ultimately, only the refractive index at the atomic transition frequency plays a role, all in concordance with the energy-time uncertainty relation. This interference process already happens for spontaneous emission in vacuum. However, when the medium has a strong jump in the density of states around the atomic transition frequency, the interference effect will change substantially.

To separate the latter cause of non-exponential decay from the former, it is interesting to consider the spontaneous emission inside homogeneous lossy dielectrics with strong and narrow material resonances, where the density of states can also change very rapidly. Here all states correspond to simple plane wave modes, so that real reabsorption processes do not play a role. In this article, we use the damped-polariton model formulated by Huttner and Barnett  to study the interference effects of spontaneous emission. If absorption is neglected in the damped-polariton model, then we are left with the Hopfield model of a dielectric , which has a frequency bandgap inside which the refractive index is purely imaginary. The analogy between this polariton band gap system and photonic crystals was drawn in .

The organization of the paper is as follows: in section II we introduce the theory and solve its equations of motion using Laplace transformations. In section III we show that the consistency of our solutions depends on the validity of a number of velocity sum rules, which are then proved. In section IV, we find that for long times all field operators can be expressed in terms of the initial bath operators, and we give an interpretation of the result. We also show how to relate the result to phenomenological quantization theories. Before we can discuss transient effects of spontaneous emission in section VI, we discuss in section V the Lorentz oscillator model and the point scattering model. We show how both these models can be found from the damped-polariton theory by choosing a suitable coupling to the bath. The paper ends with a discussion of the results and with conclusions in section VII.

## II The model and solutions of the equations of motion

The damped-polariton theory describes the interaction of light with an absorbing homogeneous medium. The coupling of the matter to a frequency continuum is the cause of the light absorption. The continuum could be a phonon bath or something else, but for the moment that is not specified: it is a collection of harmonic oscillators with a frequency-dependent coupling to the matter fields. Since the medium is homogeneous, the dynamics can be separated into a transverse and a longitudinal part. In this article we concentrate on the transverse excitations as described by the following Hamiltonian :

$$H=H_{\mathrm{em}}+H_{\mathrm{mat}}+H_{\mathrm{bath}}+H_{\mathrm{int}},$$
(1)
with

$`H_{\mathrm{em}}`$ $`=`$ $`{\displaystyle \text{d}^3k\mathrm{}\stackrel{~}{k}ca^{}(\lambda ,𝐤,t)a(\lambda ,𝐤,t)},`$ (2)
$`H_{\mathrm{mat}}`$ $`=`$ $`{\displaystyle \text{d}^3k\mathrm{}\stackrel{~}{\omega }_0b^{}(\lambda ,𝐤,t)b(\lambda ,𝐤,t)},`$ (3)
$`H_{\mathrm{bath}}`$ $`=`$ $`{\displaystyle \text{d}^3k_0^{\mathrm{}}𝑑\omega \mathrm{}\omega b_\omega ^{}(\lambda ,𝐤,t)b_\omega (\lambda ,𝐤,t)},`$ (4)
$`H_{\mathrm{int}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \text{d}^3k_0^{\mathrm{}}𝑑\omega \mathrm{}V(\omega )\left[b(\lambda ,𝐤,t)+b^{}(\lambda ,𝐤,t)\right]\left[b_\omega ^{}(\lambda ,𝐤,t)+b_\omega (\lambda ,𝐤,t)\right]}`$ (6)
$`+{\displaystyle \frac{i}{2}}{\displaystyle \text{d}^3k\mathrm{}\omega _c\sqrt{\frac{\stackrel{~}{\omega }_0}{\stackrel{~}{k}c}}\left[a(\lambda ,𝐤,t)+a^{}(\lambda ,𝐤,t)\right]\left[b^{}(\lambda ,𝐤,t)b(\lambda ,𝐤,t)\right]}.`$
We use the same notations as in . In particular, $`\stackrel{~}{k}c`$ stands for $`\sqrt{k^2c^2+\omega _c^2}`$, where the frequency $`\omega _c`$ equals $`\alpha /\sqrt{\rho \epsilon _0}`$, with $`\alpha `$ the coupling constant between field and matter, and $`\rho `$ the density. The resonance frequency $`\omega _0`$ of the polarization field is renormalized to $`\stackrel{~}{\omega }_0`$, which is the positive-frequency solution of

$$\stackrel{~}{\omega }_0^2=\omega _0^2+\stackrel{~}{\omega }_0_0^{\mathrm{}}\text{d}\omega V^2(\omega )/\omega .$$
(7)
The $`k`$-integrals in the Hamiltonian are understood to also denote a summation over the two transverse polarization directions labeled by $`\lambda `$. The creation and annihilation operators satisfy standard bosonic commutation relations. The Heisenberg equations of motion for the bath annihilation operators are:

$$\dot{b}_\omega (\lambda ,𝐤,t)=\frac{i}{2}V(\omega )\left[b(\lambda ,𝐤,t)+b^{}(\lambda ,𝐤,t)\right]i\omega b_\omega (\lambda ,𝐤,t),$$
(8)
and similarly for the creation operators. In the following we drop the $`(\lambda ,𝐤)`$-labels. We solve implicitly for the bath variables, as was done in  in a classical treatment of the model:

$$b_\omega (t)=\frac{i}{2}V(\omega )_0^t\text{d}t^{}\left[b(t^{})+b^{}(t^{})\right]e^{i\omega (tt^{})}+b_\omega (0)e^{i\omega t}.$$
(9)
The annihilation operators are defined in terms of the (transverse) physical fields:

$`a(t)`$ $`=`$ $`\sqrt{{\displaystyle \frac{\epsilon _0}{2\mathrm{}\stackrel{~}{k}c}}}\left[\stackrel{~}{k}cA(t)iE(t)\right],`$ (10)
$`b(t)`$ $`=`$ $`\sqrt{{\displaystyle \frac{\rho }{2\mathrm{}\stackrel{~}{\omega }_0}}}\left[\stackrel{~}{\omega }_0X(t)+{\displaystyle \frac{i}{\rho }}P(t)\right],`$ (11)
and similarly for the creation operators. Here $`A`$ and $`E`$ are the vector potential and the electric field, $`X`$ the polarization field and $`P`$ its canonical conjugate. Insertion of the solution (9) and its Hermitian conjugate in the equations of motion gives:

$`\dot{E}(t)`$ $`=`$ $`(\stackrel{~}{k}c)^2A(t)+(\omega _c^2/\alpha )P(t),`$ (12)
$`\dot{A}(t)`$ $`=`$ $`E(t),`$ (13)
$`\dot{X}(t)`$ $`=`$ $`(\omega _c/\alpha )^2\epsilon _0P(t)+(\omega _c^2/\alpha )\epsilon _0A(t),`$ (14)
$`\dot{P}(t)`$ $`=`$ $`\alpha ^2\stackrel{~}{\omega }_0^2/(\epsilon _0\omega _c^2)X(t)+\alpha ^2\stackrel{~}{\omega }_0/(2\epsilon _0\omega _c^2){\displaystyle _0^t}\text{d}t^{}F(tt^{})X(t^{})B(t).`$ (15)
In the last equation, the bath operator $`B(t)`$ is defined as

$$B(t)\sqrt{\frac{\mathrm{}\stackrel{~}{\omega }_0\rho }{2}}_0^{\mathrm{}}\text{d}\omega _1V(\omega _1)\left[b_{\omega _1}(0)e^{i\omega _1t}+b_{\omega _1}^{}(0)e^{i\omega _1t}\right],$$
(16)
whereas the the function $`F`$ in the convolution in (12) is:

$$F(t)2_0^{\mathrm{}}\text{d}\omega _1V^2(\omega _1)\mathrm{sin}(\omega _1t).$$
(17)
We get a system of algebraic equations by taking the Laplace transform, which we denote by a bar:

$$\left(\begin{array}{cccc}p& \stackrel{~}{k}^2c^2& 0& \omega _c^2/\alpha \\ 1& p& 0& 0\\ 0& \epsilon _0\omega _c^2/\alpha & p& \epsilon _0\omega _c^2/\alpha ^2\\ 0& 0& \alpha ^2\stackrel{~}{\omega }_0^2/\epsilon _0\omega _c^2\left[1\overline{F}(p)/(2\stackrel{~}{\omega }_0)\right]& p\end{array}\right)\left(\begin{array}{c}\overline{E}(p)\\ \overline{A}(p)\\ \overline{X}(p)\\ \overline{P}(p)\end{array}\right)=\left(\begin{array}{c}E(0)\\ A(0)\\ X(0)\\ P(0)\overline{B}(p)\end{array}\right).$$
(18)
Through the operator $`\overline{B}(p)`$ the bath remains part of the system of equations: this is as far as we can “integrate out” the bath variables.

Now we can determine the dielectric function $`\epsilon (\omega )`$, which is a classical quantity, by putting the determinant of the $`(4\times 4)`$ coefficient matrix to zero. The determinant gives the dispersion relation

$$\overline{D}(p)\overline{\epsilon }(p)p^2+k^2c^2=0,$$
(19)
with the “Laplace dielectric function”

$$\overline{\epsilon }(p)=1+\frac{\omega _c^2}{p^2+\stackrel{~}{\omega }_0^2\frac{1}{2}\stackrel{~}{\omega }_0\overline{F}(p)}.$$
(20)
The function $`\overline{F}(p)`$ is the Laplace transform of $`F(t)`$, which was defined in equation (17). From this we find the dielectric function

$$\epsilon (\omega )=\overline{\epsilon }(i\omega +\eta )=1\frac{\omega _c^2}{\omega ^2\stackrel{~}{\omega }_0^2+\frac{1}{2}\stackrel{~}{\omega }_0F(\omega )},$$
(21)
with infinitesimal positive $`\eta `$ and

$$F(\omega )\overline{F}(i\omega +\eta )=_0^{\mathrm{}}\text{d}\omega _1V^2(\omega _1)\left(\frac{1}{\omega _1\omega i\eta }+\frac{1}{\omega _1+\omega +i\eta }\right).$$
(22)
The difference between $`F(\omega )`$ and $`F(t)`$ is denoted by their arguments. The dielectric function satisfies the Kramers-Kronig-relations and has the property of a response function that $`\epsilon (\omega ^{})`$ equals $`\epsilon ^{}(\omega )`$. It can be shown that it has no poles in the upper half plane, provided that the integral in (7) exists. Previous authors  assumed that the analytical continuation of $`V^2(\omega )`$ to negative frequencies is anti-symmetrical in frequency. Then (21) reduces to the dielectric constant in , where it was shown to be identical to the more complicated expression in .

We combine (18) and (20) and write the Laplace fields in terms of the fields at time $`t=0`$, with coefficients that are functions of the Laplace dielectric function $`\overline{\epsilon }(p)`$ and susceptibility $`\overline{\chi }(p)=\overline{\epsilon }(p)1`$. For the electric field we find:

$`\overline{E}(p)`$ $`=`$ $`\overline{D}^1(p)\{\text{}pE(0)+[p^2\overline{\chi }(p)+k^2c^2]A(0)`$ (24)
$`+{\displaystyle \frac{\alpha }{\epsilon _0}}p[{\displaystyle \frac{p^2}{\omega _c^2}}\overline{\chi }(p)1]X(0)+{\displaystyle \frac{1}{\alpha }}p^2\overline{\chi }(p)[P(0)\overline{B}(p)]\}.`$
The other Laplace operators can be found in the same way and are listed in the Appendix. The inverse Laplace transform gives the fields at time $`t`$ in terms of the fields at time $`t=0`$:

$$E(t)=M_{EE}(t)E(0)+M_{EA}(t)A(0)+M_{EX}(t)X(0)+M_{EP}(t)P(0)+B_E(t),$$
(25)
where, for instance,

$$M_{EE}(t)=\frac{1}{2\pi i}_i\mathrm{}^i\mathrm{}\text{d}pe^{pt}\overline{D}^1(p)p.$$
(26)
The operator $`B_E(t)`$ in equation (25) is the contribution of the $`t=0`$ bath operators to the electric field. This term will be analyzed in more detail in section IV.

The equal time commutation relations of the field operators are

$$[A(\lambda ,𝐤,t),\epsilon _0E(\lambda ^{},𝐤^{},t)]=[X(\lambda ,𝐤,t),P(\lambda ^{},𝐤^{},t)]=i\mathrm{}\delta _{\lambda \lambda ^{}}\delta (𝐤𝐤^{}).$$
(27)
All other inequivalent combinations of operators commute. In particular $`A`$ and $`X`$ are independent canonical variables. Hence, we have the property $`[A,D]=[A,\epsilon _0E]`$, with the displacement field $`D`$ defined as $`\epsilon _0E\alpha X`$. With the help of (25) and (27), we can also calculate non-equal time commutators, for example:

$$[E(\lambda ,𝐤,t),E(\lambda ^{},𝐤^{},0)]=M_{EA}(t)[A(\lambda ,𝐤,0),E(\lambda ^{},𝐤^{},0)]=\frac{i\mathrm{}}{\epsilon _0}M_{EA}(t)\delta _{\lambda \lambda ^{}}\delta (𝐤𝐤^{}).$$
(28)
In principle we have solved the complete time evolution of the field operators. In section III we analyze in more detail their short-time behavior, whereas in section IV we consider the long-time limit.

## III Short-time limit: sum rules

For fixed $`k`$, the zeroes of the dispersion relation (19) are the poles of the integrand in (26). We assume that they are simple first-order poles and rewrite the integral (26) as an integral over frequencies $`\omega =ip`$. Then, using contour integration in the lower frequency half plane, we find the coefficients for the electric field:

$`M_{EE}(t)`$ $`=`$ $`{\displaystyle \underset{j}{}}\text{Re}\left[{\displaystyle \frac{v_{p,j}v_{g,j}}{c^2}}e^{i\mathrm{\Omega }_jt}\right],`$ (29)
$`M_{EA}(t)`$ $`=`$ $`kc{\displaystyle \underset{j}{}}\text{Im}\left[{\displaystyle \frac{v_{p,j}^2v_{g,j}}{c^3}}e^{i\mathrm{\Omega }_jt}\right],`$ (30)
$`M_{EX}(t)`$ $`=`$ $`{\displaystyle \frac{\alpha k^2c^2}{\omega _c^2\epsilon _0}}{\displaystyle \underset{j}{}}\text{Re}\left[{\displaystyle \frac{v_{p,j}v_{g,j}}{c^2}}\left(1{\displaystyle \frac{v_{p,j}^2}{c^2}}+{\displaystyle \frac{\omega _c^2}{k^2c^2}}\right)e^{i\mathrm{\Omega }_jt}\right],`$ (31)
$`M_{EP}(t)`$ $`=`$ $`{\displaystyle \frac{kc}{\alpha }}{\displaystyle \underset{j}{}}\text{Im}\left[{\displaystyle \frac{v_{g,j}}{c}}(1{\displaystyle \frac{v_{p,j}^2}{c^2}})e^{i\mathrm{\Omega }_jt}\right].`$ (32)
Some details of the calculation and a list of coefficients $`M_{mn}(t)`$ of other operators can be found in the Appendix. In these expressions, the frequencies $`\mathrm{\Omega }_j=\mathrm{\Omega }_j(k)`$ are the complex-frequency solutions of the dispersion relation $`\omega ^2\epsilon (\omega )k^2c^2=0`$. All $`\mathrm{\Omega }_j(k)`$ have a negative imaginary part. Since $`\epsilon (\omega ^{})=\epsilon ^{}(\omega )`$, it follows that $`\mathrm{\Omega }_j^{}(k)`$ is also a solution of the dispersion relation. We can choose $`\mathrm{\Omega }_j(k)`$ to be the solution with a positive real part. The summation over $`j`$ is a summation over all the polariton branches of the medium. For each branch, the complex phase velocity is defined as $`v_{p,j}(k)=\mathrm{\Omega }_j(k)/k`$ and the group velocity as $`v_{g,j}(k)=\text{d}\mathrm{\Omega }_j(k)/\text{d}k`$. For convenience, we leave out their explicit $`k`$-dependence in the following.

From equation (25) we can see that the “diagonal” coefficient $`M_{EE}(t)`$ in (29) should have the value $`1`$ at time $`t=0`$ and the “off-diagonal” coefficients $`M_{EA}(0),M_{EX}(0)`$, etc. should have the value $`0`$. The coefficients of the other field operators should also follow this rule. If these constraints are satisfied the non-equal time commutators like (28) get the right equal-time limits as well. The coefficients (29) can only have the right $`t=0`$ limits, if certain velocity sum rules are satisfied.

Velocity sum rules can be derived in a systematic way by evaluating the following two types of integrals:

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\text{d}\omega {\displaystyle \frac{(\omega +i\delta )^n}{\epsilon (\omega )\omega ^2k^2c^2}}`$ $`\text{for}n=1,0,1,`$ (33)
$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\text{d}\omega {\displaystyle \frac{(\omega +i\delta )^m\chi (\omega )}{\epsilon (\omega )\omega ^2k^2c^2}}`$ $`\text{for}m=1,0,1,2,3.`$ (34)
Here $`\epsilon (\omega )`$ is an arbitrary dielectric function that satisfies the Kramers-Kronig relations, so it is not necessarily of the specific form (21). The integrals can be evaluated using contour integration in the complex frequency plane. We can close the contours either in the upper or in the lower half plane. Equating the two answers gives a velocity sum rule. In this way one finds for all wavevectors $`k`$:

$`{\displaystyle \underset{j}{}}\text{Re}(v_{g,j}v_{p,j}/c^2)=1,`$ (35)
$`{\displaystyle \underset{j}{}}\text{Re}(v_{g,j}/v_{p,j})=1.`$ (36)
These sum rules can be found from (33) with $`n=1`$ and $`n=1`$, respectively. Both relations have been obtained before . The second was coined the Huttner-Barnett sum rule in , because of its importance in phenomenological quantum theories of dielectrics. A second group of sum rules has the form

$$\underset{j}{}\text{Im}(v_{g,j}v_{p,j}^{2q})=0q=1,0,1.$$
(37)
The rules with $`q=1`$ and $`q=1`$ follow from (34) with $`m=0`$ and $`m=2`$, respectively; the case with $`q=0`$ follows from (33) with $`n=0`$.

All of these sum rules are independent of any specific form of the dielectric function, as long as it satisfies the Kramers-Kronig relations. Other sum rules do depend on the behavior of $`\epsilon (\omega )`$ for high or low frequencies. For example, from (34) with $`m=1`$ we find:

$$\underset{j}{}\text{Re}(c^2v_{g,j}/v_{p,j}^3)=\epsilon (0).$$
(38)
This sum rule depends on the static limit of the dielectric function. For conductors the dielectric function is singular at $`\omega =0`$ , but for dielectric functions which can be found from the damped-polariton model, $`\epsilon (0)`$ is finite. Two other sum rules can be derived when for high frequencies $`\omega ^2\chi (\omega )`$ approaches a constant value that we name $`\omega _{\mathrm{lim}}^2`$. From (34) with $`m=3`$ we then find

$$\underset{j}{}\text{Re}(v_{g,j}v_{p,j}^3/c^4)=1+(\omega _{\mathrm{lim}}/kc)^2.$$
(39)
Moreover, if $`\omega ^2\chi (\omega )+\omega _{\mathrm{lim}}^2`$ falls off faster than $`\omega ^1`$, then the integral

$$_{\mathrm{}}^{\mathrm{}}\text{d}\omega \frac{\omega ^2[\omega ^2\chi (\omega )+\omega _{\mathrm{lim}}^2]}{\epsilon (\omega )\omega ^2k^2c^2}$$
(40)
produces the sum rule (37) with $`q=2`$.

Returning now to the time-dependent coefficients (29) (and the other ones in the Appendix), one finds by inspection that one needs all the above sum rules except (38) to prove that the coefficients have the right limits for $`t=0`$. In particular, from equation (21) it follows that the frequency $`\omega _{\mathrm{lim}}`$ as defined above exists in the damped-polariton model and equals $`\omega _c`$. Then with (36) and (39) we see that indeed one has $`M_{EX}(0)=0`$ in (29).

It is easy to prove the above sum rules in the following one-resonance model:

$$\epsilon (\omega )=1\frac{\omega _c^2}{\omega ^2\omega _0^2}.$$
(41)
This $`\epsilon (\omega )`$ is real and violates the Kramers-Kronig relations, but it can be considered as a limiting case of an acceptable dielectric function. The high-frequency limit of $`\omega ^2\chi (\omega )`$ indeed equals $`\omega _c^2`$. The two sum rules (35), (36) were shown to be valid for this model  and we want to check (39) as well. The dispersion relation is

$$\omega ^4(\omega _0^2+\omega _c^2+k^2c^2)\omega ^2+k^2c^2\omega _0^2=0,$$
(42)
which has two (real) solutions $`\mathrm{\Omega }_+^2`$ and $`\mathrm{\Omega }_{}^2`$ with sum $`(\omega _0^2+\omega _c^2+k^2c^2)`$ and product $`k^2c^2\omega _0^2`$. It follows that for all $`k`$

$$v_{p,+}^3v_{g,+}+v_{p,}^3v_{g,}=\frac{1}{4k^3}\frac{\text{d}}{\text{d}k}(\mathrm{\Omega }_+^4+\mathrm{\Omega }_{}^4)=\frac{1}{4k^3}\frac{\text{d}}{\text{d}k}\left[(\mathrm{\Omega }_+^2+\mathrm{\Omega }_{}^2)^22\mathrm{\Omega }_+^2\mathrm{\Omega }_{}^2\right]=\left(1+\frac{\omega _c^2}{k^2c^2}\right)c^4,$$
(43)
in agreement with (39). The other sum rules can also be checked for this simple model. The sum rules (37) obviously hold, because all group and phase velocities are real in this model. In models that respect the Kramers-Kronig-relations, these sum rules are nontrivial.

## IV Long-time limit

### A Field and medium operators

The coefficients $`M_{EE}(t)`$ etc. in (29) damp out exponentially in time. Every polariton branch has its own characteristic damping time $`\tau _j(k)=1/(\text{Im}\mathrm{\Omega }_j(k))`$. After a few times the maximum characteristic damping period, with the maximum taken over all branches, the exponentially damped coefficients can be neglected. We call this the long-time limit. The speed at which it is attained, depends on $`\epsilon (\omega )`$ and on $`k`$. For long times, only the bath operator $`B_E(t)`$ in (25) survives, because it has poles on the imaginary axis in the complex $`p`$-plane:

$$B_E(t)=\frac{1}{2\pi i\omega _c}\sqrt{\frac{\mathrm{}\stackrel{~}{\omega }_0}{2\epsilon _0}}_0^{\mathrm{}}\text{d}\omega _1V(\omega _1)_i\mathrm{}^i\mathrm{}\text{d}pe^{pt}\frac{p^2\overline{\chi }(p)}{\overline{\epsilon }(p)p^2+k^2c^2}\left[\frac{b_{\omega _1}(0)}{p+i\omega _1}+\frac{b_{\omega _1}^{}(0)}{pi\omega _1}\right].$$
(44)
Hence, in the long-time limit, all field operators are functions of the initial bath operators alone. For the electric field we find

$$E(t)E_l(t)=\frac{1}{\omega _c}\sqrt{\frac{\mathrm{}\stackrel{~}{\omega }_0}{2\epsilon _0}}_0^{\mathrm{}}\text{d}\omega _1V(\omega _1)\left[\frac{\omega _1^2\chi (\omega _1)b_{\omega _1}(0)e^{i\omega _1t}}{\epsilon (\omega _1)\omega _1^2k^2c^2}+\frac{\omega _1^2\chi ^{}(\omega _1)b_{\omega _1}^{}(0)e^{i\omega _1t}}{\epsilon ^{}(\omega _1)\omega _1^2k^2c^2}\right],$$
(45)
where the subscript $`l`$ denotes the long-time limit. The temporal (and spatial) Fourier components of the long-time solutions are:

$`E_l^+(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{\omega _c}}\sqrt{{\displaystyle \frac{\mathrm{}\stackrel{~}{\omega }_0}{2\epsilon _0}}}{\displaystyle \frac{V(\omega )\omega ^2\chi (\omega )b_\omega (0)}{\epsilon (\omega )\omega ^2k^2c^2}},`$ (46)
$`A_l^+(\omega )`$ $`=`$ $`{\displaystyle \frac{i}{\omega _c}}\sqrt{{\displaystyle \frac{\mathrm{}\stackrel{~}{\omega }_0}{2\epsilon _0}}}{\displaystyle \frac{V(\omega )\omega \chi (\omega )b_\omega (0)}{\epsilon (\omega )\omega ^2k^2c^2}},`$ (47)
$`X_l^+(\omega )`$ $`=`$ $`{\displaystyle \frac{1}{\alpha \omega _c}}\sqrt{{\displaystyle \frac{\mathrm{}\stackrel{~}{\omega }_0\epsilon _0}{2}}}{\displaystyle \frac{V(\omega )(\omega ^2k^2c^2)\chi (\omega )b_\omega (0)}{\epsilon (\omega )\omega ^2k^2c^2}},`$ (48)
$`P_l^+(\omega )`$ $`=`$ $`{\displaystyle \frac{i\alpha }{\omega _c^3}}\sqrt{{\displaystyle \frac{\mathrm{}\stackrel{~}{\omega }_0}{2\epsilon _0}}}{\displaystyle \frac{V(\omega )(\omega ^2k^2c^2\omega _c^2)\omega \chi (\omega )b_\omega (0)}{\epsilon (\omega )\omega ^2k^2c^2}},`$ (49)
where the superscript $`+`$ denotes the positive-frequency component of the operator. For future reference we also give the long-time limit of the electric field operator as a function of position and time:

$`𝐄_l(𝐫,t)=\sqrt{{\displaystyle \frac{\mathrm{}\stackrel{~}{\omega }_0}{2(2\pi )^3\epsilon _0\omega _c^2}}}`$ (51)
$`\times {\displaystyle }\text{d}𝐤{\displaystyle \underset{\lambda =1,2}{}}𝐞_\lambda (𝐤){\displaystyle _0^{\mathrm{}}}\text{d}\omega _1[{\displaystyle \frac{V(\omega _1)\omega _1^2\chi (\omega _1)b_{\omega _1}(\lambda ,𝐤,0)}{\epsilon (\omega _1)\omega _1^2k^2c^2}}e^{i(𝐤𝐫\omega _1t)}+\text{H.c.}].\text{ }`$
Similar expressions can be given for the other operators. Notice that these long-time solutions indeed are solutions of the equations of motion (12) and of the Maxwell equations. The canonical commutation relations (27) should be preserved in this long-time limit. Also, the non-equal time commutation relations like in equation (28) should be time-translation invariant. The commutation relations can be verified with the equality

$$\frac{\pi \stackrel{~}{\omega }_0}{2\omega _c^2}V^2(\omega )\left|\chi (\omega )\right|^2=\text{Im}\chi (\omega )=\text{Im}\epsilon (\omega )\epsilon _i(\omega ),$$
(52)
which follows from equations (21) and (22). Since $`\epsilon _i(\omega )`$ is anti-symmetric in $`\omega `$, all commutators can be shown to be proportional to integrals over the whole real frequency axis. Contour integration then leads to the required results.

The solutions found above can be related to those obtained by explicit diagonalization of the full Hamiltonian of the model. In  this diagonalization was carried out by using Fano’s technique. In that way the field and medium operators were written in terms of the diagonalizing annihilation operators (called $`C(𝐤,\omega )`$ in ) and the corresponding creation operators. If one replaces the bath annihilation operators $`b_\omega (𝐤,0)`$ in the long-time solutions (46) by the diagonalizing annihilation operators $`C(𝐤,\omega )`$, and if one makes similar replacements for the creation operators, the expressions for the field and medium operators in  are recovered.

The long-time solutions can be interpreted as follows: when the dielectric medium is prepared in a state that is not an eigenstate of the Hamiltonian and if the coupling $`V(\omega )`$ is nonzero for all frequencies, then the medium tends to an equilibrium that is determined by the state of the bath. The time it takes for this equilibrium to settle down is the time after which the long-time solutions can be used for the field operators. So one can always use the long-time solutions in the calculations, unless the medium has been specially prepared in a non-equilibrium state a short time before one does the experiment. The interpretation of the long-time solution will become clearer in section VI where we calculate spontaneous emission.

In summary, for times long after $`t=0`$, all field operators can be expressed solely in terms of the bath operators at time $`t=0`$. The time evolution is governed by the bath Hamiltonian alone. The field operators still satisfy Maxwell’s equations and the canonical commutation relations. Classical expressions for the Maxwell fields would have died exponentially to zero in this long time limit.

### B Relation with phenomenological theories

The long-time solutions of the field operators can be related to expressions in phenomenological theories, as we will show presently. In phenomenological quantum mechanical theories of homogeneous absorbing dielectrics , a noise current density operator $`𝐉`$ is added to the Maxwell equations in order to preserve the field commutation relations:

$`\times 𝐄^+(𝐫,\omega )`$ $`=`$ $`i\omega 𝐁^+(𝐫,\omega ),`$ (53)
$`\times 𝐁^+(𝐫,\omega )`$ $`=`$ $`i\omega \mu _0\stackrel{~}{𝐃}^+(𝐫,\omega )+\mu _0𝐉^+(𝐫,\omega ).`$ (54)
The displacement field $`\stackrel{~}{𝐃}^+`$ in the last equation is defined in terms of the electric field and the dielectric function as

$$\stackrel{~}{𝐃}^+(𝐫,\omega )=\epsilon _0\epsilon (\omega )𝐄^+(𝐫,\omega ).$$
(55)
We write $`\stackrel{~}{𝐃}`$ to stress the difference with the microscopic displacement field $`𝐃`$ in section II. After taking the spatial Fourier transform, and using $`𝐁^+=\times 𝐀^+`$ and $`𝐄^+=i\omega 𝐀^+`$, so that the first of the Maxwell equations is satisfied, one finds from the second equation:

$$\left[\omega ^2\epsilon (\omega )k^2c^2\right]A^+(\lambda ,𝐤,\omega )=\frac{1}{\epsilon _0}J^+(\lambda ,𝐤,\omega ).$$
(56)
The vector potential and all Maxwell fields can be calculated in terms of the noise current density $`J`$. The canonical commutation relations are preserved, if for the noise current one chooses :

$$[J^+(\lambda ,𝐤,\omega ),[J^+(\lambda ^{},𝐤^{},\omega ^{})]^{}]=\frac{\mathrm{}\omega ^2\epsilon _0\epsilon _i(\omega )}{\pi }\delta _{\lambda \lambda ^{}}\delta (𝐤𝐤^{})\delta (\omega \omega ^{}).$$
(57)
Instead of using the noise current operator, one defines basic bosonic operators

$$f(\lambda ,𝐤,\omega )=\sqrt{\frac{\pi }{\mathrm{}\omega ^2\epsilon _0\epsilon _i(\omega )}}J^+(\lambda ,𝐤,\omega ),$$
(58)
so that these operators satisfy simple commutation relations:

$$[f(\lambda ,𝐤,\omega ),f^{}(\lambda ^{},𝐤^{},\omega ^{})]=\delta _{\lambda \lambda ^{}}\delta (𝐤𝐤^{})\delta (\omega \omega ^{}).$$
(59)
Now we turn to the long-time solutions of the field operators that we determined in section IV A. The long-time solution of the vector potential in (46) obviously is a solution of the following inhomogeneous wave equation:

$$\left[\epsilon (\omega )\omega ^2k^2c^2\right]A_l^+(\lambda ,𝐤,\omega )=\frac{i}{\omega _c}\sqrt{\frac{\mathrm{}\stackrel{~}{\omega }_0}{2\epsilon _0}}V(\omega )\omega \chi (\omega )b_\omega (\lambda ,𝐤,0).$$
(60)
This kind of equation is well-known in Langevin theories : the coupling to a bath gives a damping term (here: a complex dielectric constant) in the equations of motion of the system. Besides damping, there is an extra term that is neglected classically. This term is the quantum noise operator, which features the bath operators at time $`t=0`$.

The long-time solution (60) can justify the phenomenological equation (56), if we identify

$$f_l(\lambda ,𝐤,\omega )=i\frac{V(\omega )\chi (\omega )}{|V(\omega )\chi (\omega )|}b_\omega (\lambda ,𝐤,0),$$
(61)
where we used equation (52). We see that up to a phase factor, the bath operators $`b_\omega (\lambda ,𝐤,0)`$ from the microscopic theory serve as basic bosonic operators $`f(\lambda ,𝐤,\omega )`$ in the phenomenological theories. We want to stress that the identification (61) is only valid in the long-time limit when the medium is in equilibrium with the bath.

In section II we saw that $`\epsilon _0E`$ is the canonical conjugate field of $`A`$ and that $`[A,D]`$ gives the canonical result as well. Since we can make the identification (61), the same relations hold in the phenomenological theory that was described in this section. But now let us calculate the commutator $`[A,\stackrel{~}{D}]`$ with $`\stackrel{~}{D}^+`$ defined as in equation (55) and $`\stackrel{~}{D}^{}`$ as its Hermitian conjugate. We can use the long-time solutions, because the commutation relations are preserved:

$$[A(\lambda ,𝐤,t),\stackrel{~}{D}(\lambda ^{},𝐤^{},t)]=\frac{2i\mathrm{}}{\pi }\delta _{\lambda \lambda ^{}}\delta (𝐤𝐤^{})_0^{\mathrm{}}\text{d}\omega \frac{\epsilon _r(\omega )\epsilon _i(\omega )\omega ^3}{|\epsilon (\omega )\omega ^2k^2c^2|^2},$$
(62)
with $`ϵ_r(\omega )`$ the real part of the dielectric constant. The symmetry of the integrand enables us to rewrite the right-hand side as an integral over all real frequencies. When using contour integration, one cannot replace $`\epsilon ^{}(\omega )`$ by $`\epsilon (\omega ^{})`$, but the analytical continuation to complex frequencies of $`\epsilon ^{}(\omega )=\epsilon (\omega )`$ must be used instead:

$$[A(\lambda ,𝐤,t),\stackrel{~}{D}(\lambda ^{},𝐤^{},t)]=i\mathrm{}\delta _{\lambda \lambda ^{}}\delta (𝐤𝐤^{})\underset{j}{}\text{Re}\left[\epsilon (\mathrm{\Omega }_j)v_{p,j}v_{g,j}/c^2\right],$$
(63)
where we assumed as before that all poles of the dispersion relation are first-order poles. Note that $`\epsilon (\mathrm{\Omega }_j)`$ depends on the behavior of the dielectric function in the upper half plane. Contrary to a statement in , the commutator does not give the canonical result, because in general there is no sum rule for the right-hand side of the equation. In other words, $`(D\epsilon _0E)`$ is canonically independent from $`E`$, but $`(\stackrel{~}{D}\epsilon _0E)`$ is not. The operator $`(D\stackrel{~}{D})`$ is proportional to the Langevin noise term in the wave equation for the electric field.

Now let us neglect absorption at all frequencies. Strictly speaking, the limit $`\epsilon _i(\omega )0`$ is unphysical because it violates the Kramers-Kronig relations, but the limit is sometimes taken for dielectrics that show negligible absorption at optical frequencies . When $`\epsilon (\omega )`$ becomes real, the solutions $`\mathrm{\Omega }_j`$ become real and in that limit one has $`\epsilon (\mathrm{\Omega }_j)\epsilon (\mathrm{\Omega }_j)=(c/v_{p,j})^2`$. Inserting this in (63) and using the Huttner-Barnett sum rule $`_j\text{Re}(v_{g,j}/v_{p,j})=1`$, we immediately find the canonical result for $`[A,\stackrel{~}{D}]`$. We compare this with the results in , where the dielectric function is assumed to be real. There a phenomenological Lagrangian was introduced and the fields $`A`$ and $`\stackrel{~}{D}`$ were correctly identified as a canonical pair. The Huttner-Barnett sum rule was invoked to show that their commutator indeed had the canonical form. It was concluded that it is misleading that also $`[A,\epsilon _0E]`$ has the canonical form. Here we have learnt that this misleading result is not surprising: in the limit of real dielectric constants and only then, both $`[A,\epsilon _0E]`$ and $`[A,\stackrel{~}{D}]`$ can have the canonical form in the same gauge, the reason being that $`\stackrel{~}{D}`$ approaches $`D`$ in that limit.

## V Model dielectric functions

Phenomenological theories as discussed in section IV B have expressions for $`\epsilon (\omega )`$ as input. In practice, this input will be the outcome of measurements of the dielectric function. By choosing the appropriate microscopic coupling constants and resonance frequencies in the damped-polariton model, one can hope to find a given dielectric function, thus providing a connection with phenomenological theories. It was argued in  that the well-known Lorentz oscillator form of the dielectric function could not be found from the damped-polariton theory in this way. We shall reconsider this issue below.

A dielectric function that follows from the damped-polariton Hamiltonian (2) will have a single resonance, because there is only one resonance frequency $`\omega _0`$ in the matter fields. Experimentally, one may find more resonances in the $`\epsilon (\omega )`$. This should not be used as an objection to the damped-polariton model, because in principle one could easily extend the theory with more material resonances. In this section, we consider two of these one-resonance models.

### A The Lorentz oscillator model

We want to find microscopic coupling constants in the damped-polariton theory so that the resulting $`\epsilon (\omega )`$ has the following Lorentz oscillator form:

$$\epsilon _{\mathrm{Lor}}(\omega )=1\frac{\omega _{c,\mathrm{Lor}}^2}{\omega ^2\omega _{\mathrm{res}}^2+2i\omega \kappa _0}.$$
(64)
Here $`\omega _{\mathrm{res}}`$ is the resonance frequency of the medium and $`\omega _{c,\mathrm{Lor}}`$ is a frequency that is related to the coupling strength between the electromagnetic and the matter field. Identifying $`\epsilon (\omega )`$ from equation (21) with $`\epsilon _{\mathrm{Lor}}(\omega )`$, we find apart from the trivial identification $`\omega _c=\omega _{c,\mathrm{Lor}}`$

$$_0^{\mathrm{}}\text{d}\omega _1V^2(\omega _1)\left[\frac{1}{\omega _1\omega i\eta }+\frac{1}{\omega _1+\omega +i\eta }\right]=4\left(\frac{\omega }{\stackrel{~}{\omega }_0}\right)i\kappa _0+\mathrm{\Delta },$$
(65)
where the frequency shift $`\mathrm{\Delta }`$ is defined such that $`\omega _{\mathrm{res}}^2=\stackrel{~}{\omega }_0^2\stackrel{~}{\omega }_0\mathrm{\Delta }/2`$. The coupling $`V^2(\omega _1)`$ is fixed by the identification of the imaginary parts and for all frequencies it equals $`V^2(\omega _1)=4\kappa _0\omega _1/(\pi \stackrel{~}{\omega }_0)`$. However, if we insert this coupling in the equation for the real parts, we find that the frequency shift $`\mathrm{\Delta }`$ is infinitely large. Also, the renormalized frequency $`\stackrel{~}{\omega }_0`$ in equation (7) blows up. We can solve this problem by introducing a frequency cut-off in the coupling, namely $`V^2(\omega _1)=4\kappa (\omega _1)\omega _1/(\pi \stackrel{~}{\omega }_0)`$ with

$$\kappa (\omega _1)=\frac{\kappa _0\mathrm{\Omega }^2}{\mathrm{\Omega }^2+\omega _1^2}.$$
(66)
With this choice one finds $`\stackrel{~}{\omega }_0=\sqrt{\omega _0^2+2\kappa _0\mathrm{\Omega }}`$, which clearly has a strong dependence on the cut-off frequency. The shift $`\mathrm{\Delta }`$ becomes both finite and frequency-dependent:

$$\mathrm{\Delta }(\omega )=\frac{4}{\pi \stackrel{~}{\omega }_0}𝒫_0^{\mathrm{}}\text{d}\omega _1\omega _1\kappa (\omega _1)\left(\frac{1}{\omega _1\omega }+\frac{1}{\omega _1+\omega }\right)=4\kappa (\omega )\frac{\mathrm{\Omega }}{\stackrel{~}{\omega }_0}.$$
(67)
The principal value integral can be evaluated by means of contour integration in the complex frequency plane. In this way we arrive at the following expression for the dielectric function:

$$\epsilon (\omega )=1\frac{\omega _c^2}{\omega ^2\omega _0^22\mathrm{\Omega }[\kappa _0\kappa (\omega )]+2i\omega \kappa (\omega )}.$$
(68)
We can choose $`\mathrm{\Omega }`$ arbitrarily high (but finite). Quite unlike $`\stackrel{~}{\omega }_0`$, the optical resonance frequency $`\omega _{\mathrm{res}}`$ approaches $`\omega _0`$ from above, the higher we choose the cut-off, since one has: $`\omega _{\mathrm{res}}^2\omega _0^2+2\kappa _0\omega _0^2/\mathrm{\Omega }`$. Note that the dielectric function (68) has the right high-frequency limit $`\omega ^2\chi (\omega )\omega _c^2`$ as required in section III.

It is well-known that there are two branches of solutions of the dispersion relation when the dielectric function is of the form (64): there is an upper and a lower polariton branch. The dielectric function (68) gives rise to another branch: it has a purely imaginary frequency with magnitude of the order of the cut-off frequency. This “cut-off branch” has negligible $`k`$-dependence. In fact, the leading $`k`$-dependent term for large $`\mathrm{\Omega }`$ is $`2i\omega _c^2\kappa _0k^2c^2/\mathrm{\Omega }^4`$. Clearly, the group velocity on this branch is practically zero, so that the contribution of the cut-off branch to the velocity sum rules of section III can be neglected.

We conclude that high cut-off frequencies can be chosen such that in the optical frequency regime the dielectric function cannot be discerned from a Lorentz dielectric function with resonance frequency $`\omega _{\mathrm{res}}=\omega _0`$ and damping constant $`\kappa _0`$. The solutions of the dispersion relation of the upper and the lower polariton branch together satisfy the sum rules of section III.

### B The point scattering model

In general the dielectric function $`\epsilon (\omega )`$ describes the propagation of a coherent light beam in a fixed direction in an isotropic medium. A complex $`\epsilon (\omega )`$ means that there is extinction, which can be caused either by scattering or absorption, or both. The dielectric function does not contain information about the extinction mechanism. A well-known dielectric medium showing polariton behavior is the dilute gas, which can be described as a collection of point dipoles that scatter light independently. If only one type of elastic scatterers is present, each having only one resonance, then the dielectric function is given by :

$$\epsilon _{\mathrm{sc}}(\omega )=1\frac{4\pi c^2\mathrm{\Gamma }_en}{\omega ^2\omega _{\mathrm{res}}^2+\frac{2}{3}i\mathrm{\Gamma }_e\omega ^3/c},$$
(69)
where $`n=N/V`$ is the density of the scatterers (not to be confused with the refractive index $`n(\omega )`$) and $`\mathrm{\Gamma }_e=e^2/(4\pi \epsilon _0m_ec^2)`$ is the classical electron radius. This dielectric function can also be found if one supposes that the medium consists of classical harmonically bound point charges whose motion is described by the Abraham-Lorentz equation. The dielectric function (69) has the property that the corresponding $`T`$-matrix $`t(\omega )`$ satisfies the optical theorem, with $`t(\omega )`$ defined as $`\epsilon (\omega )=1nt(\omega )(c/\omega )^2`$. However, (69) is not a proper response function, since it has a pole near the very large positive imaginary frequency $`3ic/(2\mathrm{\Gamma }_e)`$. This can be related to the need for the a-causal phenomenon called pre-acceleration to avoid so-called runaway solutions of the Abraham-Lorentz equation .

Although we know that in the damped-polariton theory only proper response functions can be found, we proceed like in the previous subsection and try to find coupling constants that in the optical regime give rise to the dielectric function (69). Equating with (21) we get $`\omega _c^2=4\pi c^2\mathrm{\Gamma }_en`$ and $`V^2(\omega _1)=4\mathrm{\Gamma }(\omega _1)\omega _1^3/(3\pi \stackrel{~}{\omega }_0c)`$, with

$$\mathrm{\Gamma }(\omega _1)=\frac{\mathrm{\Gamma }_e\mathrm{\Omega }^4}{\mathrm{\Omega }^4+\omega _1^4}.$$
(70)
Here we have inserted a convenient frequency cut-off from the start in order to keep finite the frequency $`\stackrel{~}{\omega }_0`$ and the shift $`\mathrm{\Delta }`$. Contour integration gives $`\stackrel{~}{\omega }_0^2=\omega _0^2+\sqrt{2}\mathrm{\Gamma }_e\mathrm{\Omega }^3/(3c)`$, and

$$\mathrm{\Delta }(\omega )=\frac{2\sqrt{2}}{3\stackrel{~}{\omega }_0c}\mathrm{\Gamma }(\omega )\mathrm{\Omega }(\mathrm{\Omega }^2+\omega ^2).$$
(71)
The dielectric function has the form

$$\epsilon (\omega )=1\frac{4\pi c^2\mathrm{\Gamma }_en}{\omega ^2\omega _0^2+[\sqrt{2}\mathrm{\Omega }^3(\mathrm{\Omega }^2\omega ^2)/(3\omega ^2c)][\mathrm{\Gamma }_e\mathrm{\Gamma }(\omega )]+\frac{2}{3}i\mathrm{\Gamma }(\omega )\omega ^3/c}.$$
(72)
In this case, the resonance frequency shifts to frequencies lower than $`\omega _0`$ and the shift is larger for larger cut-off frequencies. However, since the classical electron radius is so much smaller than an optical wavelength, it is very well possible to choose a cut-off frequency such that $`\omega _0\mathrm{\Omega }c/\mathrm{\Gamma }_e`$. Then for optical frequencies, the dielectric function (72) is of the form (69). Note that for high frequencies $`\omega ^2\chi (\omega )\omega _c^2`$ for the dielectric function (72), but not for (69).

Again, the frequency cut-off introduces a cut-off branch. In Fig. 1 we plot the real parts of the three solutions $`\mathrm{\Omega }_j(k)`$ of the dispersion relation. As a measure of the damping, we introduce $`\kappa `$ which is given by $`\mathrm{\Gamma }_e\omega _0^2/(3c)`$. For purpose of presentation, the numerical values of both $`\omega _c`$ and $`\kappa `$ were chosen artificially large for a dilute gas. The frequencies on the cut-off branch are of the same magnitude as the cut-off frequency $`\mathrm{\Omega }`$, much higher than the optical regime. The imaginary parts of the upper and lower polariton branches are plotted in Fig. 2. The imaginary part of the cut-off branch is large negative and practically constant for parameters as given in Fig. 1. Again, since the group velocity on the cut-off branch is practically zero, the upper and lower polariton branches together satisfy the sum rules of section III. In particular, Fig. 2 illustrates that the upper and lower polariton group velocities $`v_{g,u}`$ and $`v_{g,l}`$ satisfy the sum rule $`\text{Im }(v_{g,u}+v_{g,l})=0`$.

The cut-off, which was necessary to produce the dielectric function in the damped-polariton theory, neatly removes the pre-acceleration behavior associated with a pole in the upper halfplane and leads to a good response function. The form of the coupling $`V(\omega )`$ given above (70) has the following physical interpretation. By equating the damped-polariton dielectric function with (69), we assumed that the dilute gas can be described as a homogeneous dielectric. The light scattering by the gas molecules can be accounted for by an absorptive coupling to the free electromagnetic field, as long as only single scattering of light is relevant. Then scattered light is lost for propagation in the original direction. If the matter-bath coupling is dipole coupling, then for optical frequencies the product $`V^2(\omega _1)/\omega _1`$ should be proportional to the density of states of the electromagnetic field, which goes quadratically in frequency. This is indeed the case.

## VI Spontaneous emission

The spontaneous emission rate in principle is a time-dependent quantity. In this section we investigate the transient dynamics of the spontaneous emission rate of a guest atom in an absorbing medium, when the transition frequency of the guest atom is close to a material resonance of the medium. We show how our results relate to previous treatments of spontaneous emission in absorbing dielectrics, where Fermi’s Golden Rule was used to show that the time-independent (equilibrium) value for the spontaneous emission rate equals $`\mathrm{\Gamma }_0\text{Re}[n(\omega _A)]`$ . Recently, local field effects have been included in quantum electrodynamical formulations of the problem , but we shall not focus on them in this paper.

We model the guest atom as a two-level atom with ground state $`|g`$ and excited state $`|e`$ and Hamiltonian $`H_A=\mathrm{}\omega _A|ee|`$. The medium (with fields and bath included) is described by the damped-polariton model, with Hamiltonian $`H_M`$ given by (1). The total Hamiltonian is $`H=H_0+V`$, with $`H_0=H_M+H_A`$ and $`V=𝝁_A𝐄(𝐫_A)`$, the dipole interaction between the atom and the medium; $`𝝁_A`$ is the atomic dipole moment operator and $`𝐄(𝐫_A)`$ is the electric field operator at the position $`𝐫_A`$ of the atom.

Suppose that the damped-polariton system is prepared at time $`0`$ in a state described by a density matrix $`\rho _M(0)`$. We do not assume that $`\rho _M(0)`$ commutes with $`H_M`$, nor that it factorizes into a product of a density operator for the bath and a density operator for the undamped-polariton system (as is often assumed for convenience). At time $`t_0>0`$ we bring the guest atom in its excited state, and couple it to the damped-polariton system. Using perturbation theory, one can calculate  the time-dependent probability that the guest atom has emitted a photon at time $`t>t_0`$. We define the derivative of this quantity as the instantaneous spontaneous emission rate $`\mathrm{\Gamma }(t)`$. It is given as

$$\mathrm{\Gamma }(t)=\frac{2}{\mathrm{}^2}\text{Re}_{t_0}^t\text{d}t^{}e^{i\omega _A(tt^{})}\text{Tr}[\rho _M(0)𝝁𝐄(𝐫_A,t)𝝁𝐄(𝐫_A,t^{})],$$
(73)
where $`𝝁`$ is now the dipole transition matrix element of the guest atom.

If the guest atom is excited a long time after the initial preparation of the medium, all transient effects in the electric field have damped out. Hence, the field may be replaced by its long-time limit $`𝐄_l(𝐫_A,t)`$, which is given in (51). Since $`𝐄_l`$ depends only on the bath operators at $`t=0`$, we may write (73) in the form:

$$\mathrm{\Gamma }(t)=\frac{2}{\mathrm{}^2}\text{Re}_{t_0}^t\text{d}t^{}e^{i\omega _A(tt^{})}\text{Tr}_{\mathrm{bath}}[\rho _{\mathrm{red}}(0)𝝁𝐄_l(𝐫_A,t)𝝁𝐄_l(𝐫_A,t^{})].$$
(74)
Here $`\rho _{\mathrm{red}}`$ is the reduced density matrix obtained by tracing out the electromagnetic and material degrees of freedom: $`\rho _{\mathrm{red}}(0)=\text{Tr}_{\mathrm{em},\mathrm{mat}}\rho _M(0)`$. For the special case that the initial density matrix $`\rho _M(0)`$ factorizes, the reduced density matrix is the bath density matrix $`\rho _{\mathrm{bath}}(0)`$ at $`t=0`$. In general, the initial state of the electromagnetic and material degrees of freedom at $`t=0`$ does not play a role in the emission rate.

Spontaneous emission in its pure form arises if the reduced density matrix describes the ground state of the bath. Let us assume this is indeed the case. Upon inserting (51) in (74) we can perform the $`t^{}`$-integral, the integrals over the wavevector and the summations over the polarization directions. This leads to

$$\mathrm{\Gamma }(t)=\frac{\mu ^2}{3\pi ^2\mathrm{}\epsilon _0c^3}\text{Re}_0^{\mathrm{}}\text{d}\omega \omega ^3n(\omega )\frac{\mathrm{sin}[(\omega \omega _A)(tt_0)]}{\omega \omega _A},$$
(75)
with $`n(\omega )=\sqrt{ϵ(\omega )}`$ the complex refractive index.

For times $`(tt_0)`$ that are large enough, one may replace $`\mathrm{sin}[(\omega \omega _A)(tt_0)]/(\omega \omega _A)`$ by $`\pi \delta (\omega \omega _A)`$. However, the time scale at which this replacement is valid, depends on the resonance structure of the refractive index $`n(\omega )`$. Since we want to study just this time scale, we will not make the replacement. To evaluate the integral we multiply the integrand by a convergence factor $`\mathrm{\Omega }^4/(\mathrm{\Omega }^4+\omega ^4)`$, with $`\mathrm{\Omega }\omega _A`$. The specific choice of the cut-off frequency $`\mathrm{\Omega }`$ will only affect $`\mathrm{\Gamma }(t)`$ at time differences $`tt_0`$ much smaller than a single optical cycle. We need to use a high-frequency cut-off at this point, because the dipole approximation is incorrect for high frequencies.

For the dielectric function, we take the Lorentz oscillator form (68), and we choose the cut-off frequency in that model to be identical to the one inserted in (75). In Fig. 3, we give the real part of the refractive index which clearly changes rapidly near $`\omega =\omega _0`$. It is a familiar figure and it shows that the refractive index does not change much while increasing the cut-off frequency $`\mathrm{\Omega }`$ from $`10\omega _0`$ to infinity. The density of radiative modes around the material resonance is proportional to $`\omega ^2\text{Re}[n(\omega )]`$.

With this model for the dielectric function and the parameters as in Fig. 3, we calculated $`\mathrm{\Gamma }(t)`$ in the case that the transition frequency $`\omega _A`$ exactly equals $`\omega _0`$. Since the integrand in (75) is rapidly fluctuating, it is expedient to use complex countour deformation to evaluate the integral. We add an infinitesimal positive imaginary part to the denominator and split the sine into two complex exponentials. The contour of the integral with $`\mathrm{exp}[i(\omega \omega _A)(tt_0)]`$ in the integrand is deformed towards the positive imaginary axis. The contribution from the pole arising from the convergence factor can be neglected at time scales $`tt_0\omega _A^1`$. Likewise, the integration contour of the integral with $`\mathrm{exp}[i(\omega \omega _A)(tt_0)]`$ is deformed towards the negative imaginary axis. Again, the pole contribution from the convergence factor is negligible. Further contributions, which cannot be neglected, arise from the branch cuts of $`n(\omega )=\sqrt{\epsilon (\omega )}`$ and from the pole at $`\omega _A`$. The latter contribution is easily evaluated and yields the equilibrium value $`\mathrm{\Gamma }(\mathrm{})=\mathrm{\Gamma }_0\text{Re}n(\omega _A)`$. In contrast, the branch cuts yield time-dependent contributions to $`\mathrm{\Gamma }(t)`$. For large $`\mathrm{\Omega }`$ they are situated at $`\omega _1=i\kappa _0+\sqrt{\omega _0^2\kappa _0^2}`$ and $`\omega _2=i\kappa _0+\sqrt{\omega _0^2+\omega _c^2\kappa _0^2}`$. Around $`\omega _1`$ and $`\omega _2`$, we can approximate the dielectric function by

$`\epsilon (\omega _1i\delta e^{i\phi })`$ $``$ $`i\omega _c^2e^{i\phi }/(2\delta \sqrt{\omega _0^2\kappa _0^2}),`$ (76)
$`\epsilon (\omega _2i\delta e^{i\phi })`$ $``$ $`2i\delta e^{i\phi }\sqrt{\omega _0^2+\omega _c^2\kappa _0^2}/\omega _c^2.`$ (77)
The branch cut at $`\omega _1`$ gives the following contribution to the spontaneous emission rate:

$$\frac{\mathrm{\Gamma }_0}{\pi \omega _A^3}\text{Re}\left[e^{i\pi /4\kappa _0(tt_0)i(\sqrt{\omega _0^2\kappa _0^2}\omega _A)(tt_0)}J(t)\right],$$
(78)
where $`J(t)`$ is defined as:

$$J(t)=_0^{\mathrm{}}\text{d}\lambda \frac{e^{\lambda (tt_0)}[\sqrt{\omega _0^2\kappa _0^2}i(\lambda +\kappa _0)]^3}{\sqrt{\lambda }[\sqrt{\omega _0^2\kappa _0^2}\omega _Ai(\lambda +\kappa _0)]}\left[\frac{2\omega _c^2\sqrt{\omega _0^2\kappa _0^2}+i\lambda (\lambda ^2+4\omega _0^2+\omega _c^24\kappa _0^2)}{4(\omega _0^2\kappa _0^2)+\lambda ^2}\right]^{1/2}.$$
(79)
The branch cut around $`\omega _2`$ gives a similar contribution.

The integrals arising from the branch cuts and from the imaginary axis can easily be evaluated numerically, since their integrands are no longer rapidly fluctuating. The result is the solid line in Fig. 4. We see that the spontaneous emission rate builds up until it finally reaches the time-independent equilibrium value $`\mathrm{\Gamma }_0\text{Re}n(\omega _A)`$.

The dashed line in Fig. 4 is an analytical approximation for $`\mathrm{\Gamma }(t)`$, which captures the main features of the time dependence, at least qualitatively. It is derived by retaining only the contribution (78) in the time dependent part of $`\mathrm{\Gamma }(t)`$, as this is dominant for large $`t`$. Moreover, we approximate $`J(t)`$ by the first term in its asymptotic expansion for large $`tt_0`$. In this way we arrive at the following approximate expression for $`\mathrm{\Gamma }(t)`$:

$$\mathrm{\Gamma }(t)\mathrm{\Gamma }_0\text{Re}\left[n(\omega _A)\frac{\omega _c(\omega _0^2\kappa _0^2)^{5/4}}{\sqrt{2\pi }\omega _A^3(\sqrt{\omega _0^2\kappa _0^2}\omega _Ai\kappa _0)}\frac{e^{i\pi /4\kappa _0(tt_0)i[\sqrt{\omega _0^2\kappa _0^2}\omega _A](tt_0)}}{(tt_0)^{1/2}}\right].$$
(80)
As explained, this approximation contains only the contribution from the branch cut at $`\omega _1`$; the branch cut at $`\omega _2`$ gives a faster decaying term, which goes like $`e^{\kappa _0(tt_0)}/(tt_0)^{3/2}`$. The contributions from the integrals along the imaginary axis decay even faster.

It can be seen from equation (80) that the amplitude of the time-dependent part of $`\mathrm{\Gamma }(t)`$ falls off as $`e^{\kappa _0(tt_0)}/(tt_0)^{1/2}`$ and also, that the amplitude of the extra term is largest around resonance, when $`\omega _A\sqrt{\omega _0^2\kappa _0^2}`$. Away from resonance oscillations with frequency $`\sqrt{\omega _0^2\kappa _0^2}\omega _A`$ are present. Fig. 4 shows the on-resonance case when the time-dependent term shows no oscillations, but has relatively large amplitude.

The main result of the present discussion is the time-dependence of the spontaneous emission rate. The time-independent value is not reached instantaneously, but at a time scale that is governed by the resonance characteristics of the medium. In fact, the smaller the resonance width $`\kappa _0`$, the longer it takes to reach the time-independent value. Typically, it takes $`\omega _0/\kappa _0`$ optical cycles, as follows from the exponential $`e^{\kappa _0(tt_0)}`$ in the approximate expression (80). For narrow resonances with $`\omega _0/\kappa _0`$ large, the transient dynamics may take a substantial amount of time.

## VII Discussion and conclusions

We have solved the equations of motion for the field operators in the damped polariton model using Laplace transformations. The solutions of the field and medium operators are the sum of a transient and a permanent part. The latter are expressed solely in terms of the initial bath operators. Long after the initial time all field and medium operators are functions of the bath operators alone, provided the coupling to the bath is nonzero for all frequencies. The long-time solutions satisfy quantum Langevin equations in which the initial bath operators figure as the quantum noise source. The same continuum that produces the absorption also forms the noise source that keeps the commutation relations in order. This is conceptually simpler than expressing the quantum Langevin noise in terms of the creation and annihilation operators that diagonalize the total Hamiltonian of the damped-polariton model.

The effects of the initial state of the field and medium variables on the expectation values at a later time are noticeable only during a short period that is determined by the characteristic relaxation times of the damped polariton modes. Once these transient effects have died out the expectation values are determined by the reduced density matrix which follows from the full density matrix at the initial time by taking the trace over the degrees of freedom of field and matter (without bath). If the full density matrix at the initial time factorizes, the reduced density matrix equals the initial bath density matrix.

The method of long-time solutions can be used for other dissipative quantum systems as well. For models in which the Hamiltonian can be diagonalized completely, it is an alternative to the Fano diagonalization technique. The latter can be quite complicated, whereas our long-time solutions are found after the simple inversion of a $`4\times 4`$ matrix, as one sees from section II and IV. More generally, the long-time method may be useful for dissipative systems with a bilinear coupling to a harmonic oscillator bath whose dynamics can be integrated out.

We employed the method of long-time solutions to study transient effects in a medium described by a Lorentz oscillator dielectric function. This dielectric function (and that of the point-scattering model as well) can be derived from the damped-polariton model by taking a suitable bath coupling. Although a cut-off procedure turns out to be indispensable, the essential physics in the optical regime can be represented adequately in this way. Once the connection with the damped-polariton model has been established, spontaneous emission processes by a guest atom in a Lorentz oscillator dielectric can be investigated by means of the long-time method. Although transient effects due to the initial preparation of the dielectric have damped out after a few medium relaxation periods, transient behaviour of a different type shows up in the initial stages of the decay process. This transient behaviour, which is related to the preparation of the guest atom in its excited state, leads to a non-exponential decay - or in other words to a time-dependent spontaneous emission rate - if the atomic transition frequency is near a resonance of the dielectric. The non-exponential dynamics takes place at time scales that are inversely proportional to the width of the resonance. As we have shown, the characteristics of the time-dependent decay rate can be captured in an analytic asymptotic expression of which the qualitative features are corroborated by numerical methods.

## Acknowledgements

We would like to thank Ad Lagendijk, Rudolf Sprik and Willem Vos for stimulating discussions. This work is part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie, which is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek.

## Laplace operators and time-dependent coefficients

In section II the electric field operator $`\overline{E}(p)`$ was given in terms of the operators at $`t=0`$. Here we give the analogous expressions for the other Laplace operators. Furthermore, we show how to evaluate the time-dependent coefficients $`M_{mn}(t)`$ like in (29) for the electric field operator. Finally, we list the expressions for the coefficients of the other operators.

The expression for $`\overline{E}(p)`$ in (24) has the following analogous expressions for the other Laplace operators:

$`\overline{A}(p)`$ $`=`$ $`{\displaystyle \frac{1}{\overline{D}(p)}}\{\text{}E(0)+pA(0)`$ (82)
$`{\displaystyle \frac{\alpha }{\epsilon _0}}[{\displaystyle \frac{p^2}{\omega _c^2}}\overline{\chi }(p)1]X(0){\displaystyle \frac{1}{\alpha }}p\overline{\chi }(p)[P(0)\overline{B}(p)]\},`$
$`\overline{X}(p)`$ $`=`$ $`{\displaystyle \frac{1}{\overline{D}(p)}}\{{\displaystyle \frac{\epsilon _0}{\alpha }}p\overline{\chi }(p)E(0)+{\displaystyle \frac{\epsilon _0}{\alpha }}p^2\overline{\chi }(p)A(0)`$ (84)
$`+({\displaystyle \frac{k^2c^2}{\omega _c^2}}+{\displaystyle \frac{p^2}{\omega _c^2}}+1)p\overline{\chi }(p)X(0)+{\displaystyle \frac{\epsilon _0}{\alpha ^2}}(p^2+k^2c^2)\overline{\chi }(p)[P(0)\overline{B}(p)]\},`$
$`\overline{P}(p)`$ $`=`$ $`{\displaystyle \frac{1}{\overline{D}(p)}}\{\alpha [{\displaystyle \frac{p^2}{\omega _c^2}}\overline{\chi }(p)1]E(0)+\alpha p[{\displaystyle \frac{p^2}{\omega _c^2}}\overline{\chi }(p)1]A(0)`$ (87)
$`+{\displaystyle \frac{\alpha ^2}{\epsilon _0}}\left[{\displaystyle \frac{p^2}{\omega _c^2}}\overline{\chi }(p)1\right]\left({\displaystyle \frac{k^2c^2}{\omega _c^2}}+{\displaystyle \frac{p^2}{\omega _c^2}}+1\right)X(0)`$
$`+p({\displaystyle \frac{k^2c^2}{\omega _c^2}}+{\displaystyle \frac{p^2}{\omega _c^2}}+1)\overline{\chi }(p)[P(0)\overline{B}(p)]\}.`$
If we now apply the inverse Laplace transformation to these expressions, we find the full time dependence of the operators $`A`$, $`X`$ and $`P`$. The inverse Laplace transformation is a contour integration over the Bromwich contour that includes the whole imaginary $`p`$-axis. After transforming to frequency variables, the contour includes poles from $`\overline{D}^1(p)`$, which are in the lower halfplane, and moreover poles on the real frequency axis arising from $`\overline{B}(p)`$. The latter are important in the calculation of the long-time solutions of the operators in section IV. However, in the calculation of the coefficients $`M_{mn}(t)`$, which we will discuss here, they play no role. Let us consider as an example the coefficient $`M_{AE}(t)`$:

$`M_{AE}(t)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{i\mathrm{}+\eta }^{i\mathrm{}+\eta }}\text{d}pe^{pt}\overline{D}^1(p)={\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\text{d}\omega {\displaystyle \frac{e^{i\omega t}}{\epsilon (\omega )\omega ^2k^2c^2}}`$ (88)
$`=`$ $`{\displaystyle \frac{1}{4\pi kc}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\text{d}\omega \left({\displaystyle \frac{e^{i\omega t}}{n(\omega )\omega kc}}{\displaystyle \frac{e^{i\omega t}}{n(\omega )\omega +kc}}\right)={\displaystyle \frac{1}{kc}}{\displaystyle \underset{j}{}}\text{Im}(e^{i\mathrm{\Omega }_jt}{\displaystyle \frac{v_{g,j}}{c}}).`$ (89)
Note that $`M_{AE}(t)`$ is exponentially damped because all $`\mathrm{\Omega }_j`$ in the exponentials have negative imaginary parts. The other coefficients can be calculated in a similar way. The results are:

$`M_{AA}(t)`$ $`=`$ $`M_{EE}(t),`$ (90)
$`M_{AX}(t)`$ $`=`$ $`{\displaystyle \frac{\alpha k}{\omega _c^2\epsilon _0c^4}}{\displaystyle \underset{j}{}}\text{Im}\left[e^{i\mathrm{\Omega }_jt}v_{g,j}(1{\displaystyle \frac{v_{p,j}^2}{c^2}}+{\displaystyle \frac{\omega _c^2}{k^2c^2}})\right],`$ (91)
$`M_{AP}(t)`$ $`=`$ $`{\displaystyle \frac{1}{\alpha }}{\displaystyle \underset{j}{}}\text{Re}\left[e^{i\mathrm{\Omega }_jt}{\displaystyle \frac{v_{g,j}}{c}}\left({\displaystyle \frac{v_{p,j}}{c}}{\displaystyle \frac{c}{v_{p,j}}}\right)\right],`$ (92)
$`M_{XE}(t)`$ $`=`$ $`\epsilon _0M_{AP}(t),`$ (93)
$`M_{XA}(t)`$ $`=`$ $`\epsilon _0M_{EP}(t),`$ (94)
$`M_{XX}(t)`$ $`=`$ $`{\displaystyle \frac{k^2c^2}{\omega _c^2}}{\displaystyle \underset{j}{}}\text{Re}\left[e^{i\mathrm{\Omega }_jt}{\displaystyle \frac{v_{g,j}}{c}}\left({\displaystyle \frac{v_{p,j}}{c}}{\displaystyle \frac{c}{v_{p,j}}}\right)\left(1{\displaystyle \frac{v_{p,j}^2}{c^2}}+{\displaystyle \frac{\omega _c^2}{k^2c^2}}\right)\right],`$ (95)
$`M_{XP}(t)`$ $`=`$ $`{\displaystyle \frac{\epsilon _0kc}{\alpha ^2}}{\displaystyle \underset{j}{}}\text{Im}\left[e^{i\mathrm{\Omega }_jt}{\displaystyle \frac{v_{g,j}}{c}}\left({\displaystyle \frac{v_{p,j}}{c}}{\displaystyle \frac{c}{v_{p,j}}}\right)^2\right],`$ (96)
$`M_{PE}(t)`$ $`=`$ $`\epsilon _0M_{AX}(t),`$ (97)
$`M_{PA}(t)`$ $`=`$ $`\epsilon _0M_{EX}(t),`$ (98)
$`M_{PX}(t)`$ $`=`$ $`{\displaystyle \frac{\alpha ^2k^3c^3}{\omega _c^4\epsilon _0}}{\displaystyle \underset{j}{}}\text{Im}\left[e^{i\mathrm{\Omega }_jt}{\displaystyle \frac{v_{g,j}}{c}}\left(1{\displaystyle \frac{v_{p,j}^2}{c^2}}+{\displaystyle \frac{\omega _c^2}{k^2c^2}}\right)^2\right],`$ (99)
$`M_{PP}(t)`$ $`=`$ $`M_{XX}(t).`$ (100)
With the sum rules discussed in section III, one can see that the “diagonal” coefficients in this list equal $`1`$ at time $`t=0`$, whereas the other coefficients have the initial value $`0`$.


format: text

text:
# Measuring Time-Dependence of Dark Energy Density from Type Ia Supernova Data

## 1 Introduction

Most of the present energy content of our universe is unknown (Bahcall et al. (1999)). Distance-redshift relations derived from cosmological standard candles at redshifts between zero and a few are the most sensitive probes of the equation of state of the universe, which allows us to constrain the energy content of the universe.

Type Ia supernovae (SNe Ia) are our best candidates for cosmological standard candles. They can be calibrated to have small dispersion in their intrinsic luminosities (Phillips (1993); Riess, Press, & Kirshner (1995)). The data from two independent observational teams, the High-$`z`$ SN Search (Schmidt et al.) and the Supernova Cosmology Project (Perlmutter et al.), seem to suggest that our universe has a significant vacuum energy content (Garnavich et al. (1998); Riess et al. (1998); Perlmutter et al. (1999)).

While a cosmological constant term (vacuum energy) in the Einstein equations provides the simplest explanation of the current SN Ia data, other forms of energy (quintessence, dark energy, etc) have also been studied (e.g. White (1998); Garnavich et al. (1998); Steinhardt, Wang, & Zlatev (1999); Efstathiou (1999); Podariu & Ratra (2000); Waga & Frieman (2000)) and are consistent with current data as well. Since cosmology has matured into a phenomenological science at the turn of the new millennium, observational data will dominate aesthetics in the selection of cosmological models.

So far, most cosmologists have assumed time-independent equations of state, i.e., power-law dark energy density (see §2), in the context of constraining the energy content of the universe. In this paper, we allow the dark energy density of the universe to be a free function of redshift, i.e., arbitrary equation of state. We derive constraints on the dark energy density from current SN data, and assess future prospects of measuring the dark energy density using simulated data from realistic future SN Ia surveys. We expect a model-independent measurement of the dark energy density as a function of redshift to be very useful in our quest of the unknown energy contents of the universe.

## 2 Parametrization of the dark energy density

In a smooth Friedmann-Robertson-Walker (FRW) universe, the metric is given by $`ds^2=dt^2a^2(t)[dr^2/(1kr^2)+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)]`$, where $`a(t)`$ is the cosmic scale factor, and $`k`$ is the global curvature parameter. The cosmological redshift $`z`$ is given by $`1+z=1/a`$.

To make model-independent measurements of the equation of state, we replace the vacuum energy density $`\rho _\mathrm{\Lambda }`$ with $`\rho _X=\rho _X^0f(z)`$ in the total matter density of the universe:

$`\rho (z)`$ $`=`$ $`\rho _m^0(1+z)^3+\rho _k^0(1+z)^2+\rho _X^0f(z)`$ (1)
$`=`$ $`\rho _c^0\left[\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_k(1+z)^2+\mathrm{\Omega }_Xf(z)\right],`$
where the superscript “0” indicates present values, $`f(z=0)=1`$, and $`\mathrm{\Omega }_k=1\mathrm{\Omega }_m\mathrm{\Omega }_X=k/H_0^2`$. If the unknown energy is due to a cosmological constant $`\mathrm{\Lambda }`$, $`f(z)=1`$. Clearly, the function $`f(z)`$ is a very good probe of the nature of the unknown energy.

The comoving distance $`r`$ is given by (Weinberg (1972))

$$r(z)=cH_0^1\frac{S(\kappa \mathrm{\Gamma })}{\kappa },\kappa \left|\mathrm{\Omega }_k\right|^{1/2},$$
(2)
$$\mathrm{\Gamma }(z;\mathrm{\Omega }_m,\mathrm{\Omega }_X,f)=_0^z𝑑z^{}\frac{1}{E(z^{})},$$
(3)
$$E(z)\left[\mathrm{\Omega }_m(1+z^{})^3+\mathrm{\Omega }_Xf(z^{})+\mathrm{\Omega }_k(1+z^{})^2\right]^{1/2},$$
(4)
where

$`S(x)`$ $`=`$ $`\mathrm{sinh}(x),\mathrm{\Omega }_k>0`$ (5)
$`=`$ $`x,\mathrm{\Omega }_k=0`$
$`=`$ $`\mathrm{sin}(x),\mathrm{\Omega }_k<0.`$
The angular diameter distance is given by $`d_A(z)=r(z)/(1+z)`$, and the luminosity distance is given by $`d_L(z)=(1+z)^2d_A(z)`$.

Einstein’s equations in a FRW metric, together with the first law of thermodynamics give us

$$(1+z)\frac{d\rho }{dz}=3(\rho +p).$$
(6)
For unknown energy $`\rho _X(z)`$, we find

$`\rho _X(z)`$ $`=`$ $`\rho _X^0f(z)`$
$`p_X(z)`$ $`=`$ $`\rho _X^0\left[{\displaystyle \frac{1}{3}}(1+z)f^{}(z)f(z)\right].`$ (7)
Now we can write the equation of state of the unknown energy as

$$w(z)\frac{p_X(z)}{\rho _X(z)}=\frac{1}{3}(1+z)\frac{f^{}(z)}{f(z)}1.$$
(8)
A constant equation of state corresponds to $`f(z)(1+z)^\alpha `$, where $`\alpha `$ is a constant. The values $`\alpha =0`$, $`\alpha =3`$, and $`\alpha =4`$ correspond to a cosmological constant, matter, and radiation respectively.

To obtain accelerated expansion, we need $`\rho +3p<0`$. Since $`\rho _X+3p_X=\rho _X^0\left[(1+z)f^{}2f\right]`$, this implies $`\alpha <2`$ for $`f(z)(1+z)^\alpha `$.

The weak energy condition states that for all physically reasonable classical matter, the energy density of matter as measured by any observer is nonnegative (Wald (1984)). For a perfect fluid, the weak energy condition will be satisfied if and only if

$$\rho +p0.$$
(9)
This leads to (see Eq.)

$$f^{}(z)0.$$
(10)
The weak energy condition imposes strong constraints on the jointly estimated cosmological parameters $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, and the dark energy density $`f(z)`$.

Note that the comoving distance $`r(z)`$ depends on the equation of state of X through

$$\mathrm{\Gamma }(z)_0^z𝑑z^{}\left[\mathrm{\Omega }_m(1+z^{})^3+\mathrm{\Omega }_Xf(z^{})+\mathrm{\Omega }_k(1+z^{})^2\right]^{1/2}=H_0_0^r\frac{dr}{\sqrt{1kr^2}},$$
(11)
Hence

$$\mathrm{\Gamma }^{}(z)=\left[\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_Xf(z)+\mathrm{\Omega }_k(1+z)^2\right]^{1/2}=\frac{H_0}{\left[1+\mathrm{\Omega }_kH_0^2r^2\right]^{1/2}}\frac{dr}{dz}.$$
(12)
To measure $`f(z)`$ directly from data, we need to evaluate the derivative of the distance $`r(z)`$ with respect to redshift $`z`$.

To avoid taking derivatives of noisy data, we can parametrize $`f(z)`$ with its values at $`n`$ equally spaced redshifts $`z_i`$, and assume that $`f(z)`$ is given by linear interpolations at other values of $`z`$. We write

$`f(z)=\left({\displaystyle \frac{z_iz}{z_iz_{i1}}}\right)f_{i1}+\left({\displaystyle \frac{zz_{i1}}{z_iz_{i1}}}\right)f_i,z_{i1}<zz_i,`$
$`z_0=0,z_n=z_{max};f_0=1,f_n=f_{n1}`$ (13)
where $`f_i`$ ($`i=1,2,\mathrm{},n1)`$ are independent variables to be estimated from data.

## 3 Parameter estimation

The measured distance modulus for a SN Ia is

$$\mu _0^{(l)}=\mu _p^{(l)}+ϵ^{(l)}$$
(14)
where $`\mu _p^{(l)}`$ is the theoretical prediction

$$\mu _p^{(l)}=5\mathrm{log}\left(\frac{d_L(z_l)}{\text{Mpc}}\right)+25,$$
(15)
and $`ϵ^{(l)}`$ is the uncertainty in the measurement, including observational errors and intrinsic scatters in the SN Ia absolute magnitudes.

Denoting all the parameters to be fitted as s, we can estimate s using a $`\chi ^2`$ statistic, with (Riess et al. 1998)

$$\chi ^2(\text{s})=\underset{l}{}\frac{\left[\mu _p^{(l)}(z_l|\text{s})\mu _0^{(l)}\right]^2}{\sigma _{\mu _0,l}^2+\sigma _{mz,l}^2}\underset{l}{}\frac{\left[\mu _p^{(l)}(z_l|\text{s})\mu _0^{(l)}\right]^2}{\sigma _l^2},$$
(16)
where $`\sigma _{\mu _0}`$ is the estimated measurement error of the distance modulus, and $`\sigma _{mz}`$ is the dispersion in the distance modulus due to the dispersion in galaxy redshift, $`\sigma _z`$, due to peculiar velocities and uncertainty in the galaxy redshift (for the Perlmutter et al. data, the dispersion due to peculiar velocities is included in $`\sigma _{m_B^{eff}}`$, i.e., $`\sigma _{\mu _0}`$). Since

$$\sigma _{mz}=\frac{5}{\mathrm{ln}10}\left(\frac{1}{d_L}\frac{d_L}{z}\right)\sigma _z,$$
(17)
$`\sigma _{mz}`$ depends on the parameters s. The probability density function (PDF) for the parameters s is

$$p(\text{s})\mathrm{exp}\left(\frac{\chi ^2}{2}\right).$$
(18)
The normalized PDF is obtained by dividing the above expression by its sum over all possible values of the parameters s.

In order to impose the weak energy condition $`f^{}(z)0`$, we compute the PDFs on a $`N`$-dimensional grid for $`N`$ parameters. The PDF of a given parameter $`s_i`$ is obtained by integrating over all possible values of the other $`N1`$ parameters. To reduce the computation time, we can integrate over the Hubble constant $`H_0`$ analytically, and define a modified $`\chi ^2`$ statistic, with

$$\stackrel{~}{\chi }^2\chi _{}^2\frac{C_1}{C_2}\left(C_1+\frac{2}{5}\mathrm{ln}10\right),$$
(19)
where

$`\chi _{}^2`$ $``$ $`{\displaystyle \underset{l}{}}{\displaystyle \frac{1}{\sigma _l^2}}\left(\mu _p^{(l)}\mu _0^{(l)}\right)^2,`$
$`C_1`$ $``$ $`{\displaystyle \underset{l}{}}{\displaystyle \frac{1}{\sigma _l^2}}\left(\mu _p^{(l)}\mu _0^{(l)}\right),`$
$`C_2`$ $``$ $`{\displaystyle \underset{l}{}}{\displaystyle \frac{1}{\sigma _l^2}},`$ (20)
where

$$\mu _p^{}\mu _p(h=h^{})=42.3845\mathrm{log}h^{}+5\mathrm{log}\left[H_0r(1+z)\right].$$
(21)
We take $`h^{}=0.65`$. Our results are independent of the choice of $`h^{}`$.

After experimenting with a number of different techniques, we developed an adaptive iteration method of estimating $`f(z)`$ based on the requirement that $`f^{}(z)0`$ (i.e., the weak energy condition is satisfied). Starting with the initial guess of $`f(z)=f(z=0)=1`$ (a cosmological constant), we iteratively build up $`f(z)`$ as parametrized by Eq.(2) while minimizing $`\chi ^2`$.

## 4 Constraints of $`\rho _X(z)`$ from current SNe Ia data

Wang (2000b) has combined the data of the High-$`z`$ SN Search team (Schmidt et al. (1998); Garnavich et al. (1998); Riess et al. (1998)) and the Supernova Cosmology Project (Perlmutter et al. (1999)), yielding a total of 92 SNe Ia. Using the method described in the previous section, we estimate $`\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_X`$, and $`f(z)`$ simultaneously by minimizing the modified $`\chi ^2`$ (see Eq.(19)).

Fig.1 shows the dimensionless dark energy density $`f(z)`$ (as parametrized by Eq.(2)) measured from this set of 92 SNe Ia. It is consistent with $`f(z)=1`$ (a cosmological constant), with large uncertainty beyond $`z\begin{array}{c}>\hfill \\ \hfill \end{array}0.5`$. Note that none of the error bars extend beneath $`f(z)=1`$, because we have imposed the weak energy condition, i.e., $`f^{}(z)0`$, which implies that $`f(z)f(z=0)=1`$. The estimated values of $`\mathrm{\Omega }_m=.3`$ (0, .9); $`\mathrm{\Omega }_X=1.7`$ (.4, 2.2) are consistent with previous results (Garnavich et al. (1998); Riess et al. (1998); Perlmutter et al. (1999); Wang 2000b ). The errors are estimated from the ranges of parameters for which $`\chi ^2=\chi _{min}^2+1`$.

Wang (2000b) found that when fit to a model with a cosmological constant as the dark energy, flux-averaging changes the best fit model to this data set of 92 SNe Ia. Without flux-averaging, the best fit model is a closed universe with $`\mathrm{\Omega }_m=0.7\pm 0.4`$, and a vacuum energy density fraction of $`\mathrm{\Omega }_\mathrm{\Lambda }=1.2\pm 0.5`$, consistent with the estimated values in Fig.1. The flux-averaged data yield $`\mathrm{\Omega }_m=0.3\pm 0.6`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7\pm 0.7`$. This difference may have resulted from the data containing large redshift dependent uncertainties, which would have caused the results from the unbinned data to be biased. The effect of flux-averaging on the best fit model assuming an arbitrary dimensionless dark energy density $`f(z)`$ will be studied elsewhere.

Future cosmic microwave background (CMB) space missions MAP (Bennett et al. (1997)) and Planck (De Zotti et al. (1999)), together with the galaxy redshift surveys SDSS (Gunn (1999)) and 2df (Dalton et al. (2000)), will give us exquisitely accurate measurements of the geometry of the universe and the matter density in the universe (Eisenstein, Hu, & Tegmark (1999); Turner & Tyson (1999); Wang, Spergel, & Strauss (1999)). SN data can provide the unique probe on the nature of dark energy by allowing us to measure how the dark energy density varies with time.

Current cosmic microwave background (CMB) anisotropy measurements seem to indicate that we live in a flat universe (de Bernardis et al. (2000); Balbi et al. (2000)). Cluster abundances strongly suggest a low matter density universe (Bahcall, Lubin, & Dorman (1995); Carlberg et al. (1996); Bahcall & Fan (1998)). $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ is the best fit model to current observational data. We will use $`\mathrm{\Omega }_m=0.3`$ and $`\mathrm{\Omega }_X=0.7`$ for our simulated data in the rest of this paper.

## 5 Measuring $`\rho _X(z)`$ from future SNe Ia data

A large number of SNe Ia at $`z\begin{array}{c}>\hfill \\ \hfill \end{array}1`$ is critical in resolving the important systematic uncertainties of SNe Ia as cosmological standard candles, such as dust (Aguirre (1999)), gravitational lensing (Kantowski, Vaughan, & Branch (1995); Wambsganss et al. (1997); Holz (1998); Metcalf & Silk (1999); Wang (1999); Barber et al. (2000)), and luminosity evolution (Drell, Loredo, & Wasserman (2000); Riess et al. (1999); Wang 2000b ), and in making SNe Ia useful probes of the dark energy content of the universe. The most efficient method of obtaining a large number of SNe Ia at $`z>1`$ is conducting a supernova pencil beam survey on a dedicated large aperture telescope with a square degree field of view (Wang 2000a ).

To study how well SN data can probe the dark energy density, let us consider two hypothetical dimensionless dark energy densities $`f_q(z)`$ and $`f_k(z)`$, given by

$`f_q(z)`$ $`=`$ $`{\displaystyle \frac{e^{1.5z}}{(1+z)^{1.5}}},w_q(z)=1+0.5z`$
$`f_k(z)`$ $`=`$ $`\mathrm{exp}[0.9(1e^z)],w_k(z)=0.3(1+z)e^z1`$ (22)
We have chosen $`f_q(z)`$ and $`f_k(z)`$ to represent quintessence ($`dw_q/dz>0`$) and k-essence ($`dw_k/dz<0`$) models respectively (Caldwell, Dave, & Steinhardt (1998); Armendariz-Picon, Mukhanov, & Steinhardt (2000)). Note that $`f_q(z)`$ and $`f_k(z)`$ satisfy the weak energy condition $`f^{}(z)0`$; they give an accelerating universe for $`z\begin{array}{c}<\hfill \\ \hfill \end{array}1.33`$ and $`z\begin{array}{c}<\hfill \\ \hfill \end{array}2`$ respectively.

A feasible SN pencil beam survey (either from ground222The SNe Ia at $`z\begin{array}{c}>\hfill \\ \hfill \end{array}1.5`$ will likely require follow up spectroscopy from space. or from space), with a square degree field of view and for an effective observational period of one year, can yield almost 2000 SNe Ia out to $`z=2`$ (Wang 2000a ). Let us combine the data from the SN pencil beam survey with SN data at smaller redshifts, so that there are a minimum of 50 SNe Ia per 0.1 redshift interval at any redshift. This yields a total of 1966 SNe Ia for the quintessence and 1898 SNe Ia for the k-essence model, up to $`z=2`$ and for $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_X=0.7`$. We simulate the data by placing perfect standard candles at random redshifts, with the total number per 0.1 redshift interval given as described above. Then we add intrinsic and observational dispersions which are Gaussian with zero mean and a variance of 0.20 magnitudes. A systematic shift in $`\mu _0`$ of $`dm_{sys}z`$ magnitudes is also added to mimic possible systematic errors as one goes to larger redshifts.

Realistic SN Ia data should contain gravitational lensing noise. Wang (2000a,b) has shown that flux averaging is important in reducing the bias due to lensing. The effect of gravitational lensing and of flux averaging in the context of a general equation of state will be investigated elsewhere.

Fig.2 illustrates how well we can recover the dark energy densities $`f_q(z)`$ and $`f_k(z)`$ in the absence of lensing noise and systematic shifts, when we apply our adaptive iteration method to 100 random data sets with a realistic dispersion of 0.2 magnitudes. To study the dependence of our results on the parametrization of $`f(z)`$ \[see Eq.(2)\], we show results for (a) $`n=10`$, and (b) $`n=6`$. The thick and thin solid lines are the assumed true $`f_q(z)`$ and $`f_k(z)`$ respectively. We have assumed that we know $`\mathrm{\Omega }_m+\mathrm{\Omega }_X=1`$. The mean and 1-$`\sigma `$ errors of $`f(z)`$ and $`\mathrm{\Omega }_m`$ are estimated from averaging over the 100 random samples. The error of estimated $`\mathrm{\Omega }_m`$ reflects the resolution of the program. For a given data set, the recovered $`f(z)`$ and $`\mathrm{\Omega }_m`$ are expected to fall within the errors with 68% probability. Clearly, the quintessence and k-essence models can be differentiated marginally for $`z\begin{array}{c}<\hfill \\ \hfill \end{array}1.5`$. For $`z\begin{array}{c}>\hfill \\ \hfill \end{array}1.5`$, the errors increase significantly, while the estimates become more biased, making it impossible to differentiate between the two models. However, the correct trend in the time variation of the dark energy density can be clearly detected out to $`z=2`$.

Fig.2 shows that the parametrization of $`f(z)`$ with $`n=10`$ yields less biased estimates than $`n=6`$. This is as expected, since $`f(z)`$ is more accurately parametrized as one increases $`n`$. However, the errors in the estimates increase with $`n`$ as well. One must choose an optimal $`n`$ such that $`f(z)`$ is adequately parametrized, while the errors on the estimated amplitudes of $`f(z)`$ are not too large to be useful. We’ve experimented with $`n>10`$ parametrizations of $`f(z)`$, and found that $`n=10`$ is a good choice.

The biased estimates of the dark energy density $`f_q(z)`$ and $`f_k(z)`$ in Fig.2 is mainly due to the bias in the estimate of $`\mathrm{\Omega }_m`$. Fig.3 shows $`f_q(z)`$ and $`f_k(z)`$ estimated assuming that we know $`\mathrm{\Omega }_m=0.3`$, with the same line types as in Fig.2. Even with this ideal assumption, it is only possible to marginally differentiate between the two models.

Fig.4 shows the effect of adding a systematic shift of $`dm_{sys}z`$ to 100 random data sets with a realistic dispersion of 0.2 magnitudes. The line types are the same as in Fig.2. We have added a systematic shift in $`\mu _0`$ of (a) $`0.01z`$ magnitudes; (b) $`0.05z`$ magnitudes. Clearly, systematic shifts can significantly increase the bias in the estimate of $`\mathrm{\Omega }_m`$ and the estimates of $`f(z)`$ for $`z\begin{array}{c}>\hfill \\ \hfill \end{array}1.2`$, while having little effect on the estimates of $`f(z)`$ at $`z\begin{array}{c}<\hfill \\ \hfill \end{array}1.2`$. Systematic errors as a function of $`z`$ may arise from intrinsic properties of the supernovae varying with $`z`$ such as progenitor evolution or dust characteristics changing with metalicity of the Universe. Accuracy may also be limited by the observations themselves in the case of k-corrections or selection biases. The estimated errors from these sources in the current surveys are between 2% and 10% (Schmidt et al. (1998)). It is perhaps not realistic to expect to measure both $`\mathrm{\Omega }_m`$ and $`f(z)`$ to high redshift accurately from SN Ia data alone.

It is clear from Fig.4 that even in the presence of plausible systematic effects, we can expect to measure the time-variation in the dark energy density $`f(z)`$ with reasonable accuracy to a redshift of about 1.2. The bias and the errors in the estimated $`f(z)`$ increase substantially beyond $`z=1.2`$.

We only applied our adaptive iteration method to 100 random samples, and with a resolution of $`\mathrm{\Delta }\mathrm{\Omega }_m=0.02`$, because this method takes several hours per sample on a fast Sun work station in finding the best fit $`f(z)`$ and $`\mathrm{\Omega }_m`$. It is work planned for the future to improve this promising method for application to much larger number of random samples, as well as adding $`\mathrm{\Omega }_X`$ as an estimated parameter.

## 6 Implication for dark energy models

Recently, there has been a great deal of activity in exploring the possibilities of the existence of exotic dark energy (Peebles & Ratra (1988); Frieman et al. (1995); Caldwell, Dave, & Steinhardt (1998); Sahni & Wang (2000)) in the universe. While the present observational data are consistent with the dark energy being a cosmological constant, they do not rule out alternatives in the form of various scalar fields.

It is important that we measure the time dependence of the dark energy density. If the dark energy density is measured to be constant in time within reasonable uncertainties, a cosmological constant should be favored, and more theoretical efforts should be directed toward the derivation of a cosmological constant from first principles. At the very least, this places strong constraint on the classes of scalar-field models for the dark energy. On the other hand, if the time-dependence of the dark energy density is established by the observational data, we would come to the exciting discovery of new physics in the universe.

We have studied two dark energy models $`f_q(z)`$ and $`f_k(z)`$ (see Eq.(5)), representing two general classes of dark energy models, quintessence (Caldwell, Dave, & Steinhardt (1998)) and k-essence (Armendariz-Picon, Mukhanov, & Steinhardt (2000)). This allows us to examine how well different models can be differentiated by realistic data, as well as the robust determination of the time-dependence of the dark energy density.

Realistic future SN data, as described in the previous section, has the potential of determining the time-dependence of the dark energy density. This will clearly have a dramatic impact on models of the dark energy.

## 7 Conclusions

To access the prospects of measuring the time variation in the equation of state, we have developed a promising adaptive iteration method that is powerful in extracting the dark energy density $`f(z)=\rho _X(z)/\rho _X(z=0)`$ from realistic data. This method is based on the requirement that the weak energy condition (energy density of matter is nonnegative for any observer) is satisfied.

We have found that current type Ia supernova (SNe Ia) data are consistent with a cosmological constant, with large uncertainties at $`z\begin{array}{c}>\hfill \\ \hfill \end{array}0.5`$. We show that $`\mathrm{\Omega }_m`$ (assuming a flat universe) and the dimensionless dark energy density $`f(z)=\rho _X(z)/\rho _X(z=0)`$ can be measured reasonably well to about $`z=1.5`$ using type Ia supernova data from realistic future SN Ia pencil beam surveys, provided that the weak energy condition (energy density of matter is nonnegative for any observer) is imposed. For $`z\begin{array}{c}>\hfill \\ \hfill \end{array}1.5`$, the errors increase significantly, while the estimates become more biased, making it impossible to differentiate between different models (say, quintessence and k-essence). However, the correct trend in the time-dependence of the dark energy density can be clearly detected out to $`z=2`$, even in the presence of plausible systematic effects. This would allow us to determine whether the dark energy is vacuum energy, or some exotic form of energy with a time-dependent density.

The simulated data we used are for a SN pencil beam survey (either from ground or from space) with a square degree field of view and for an effective observational period of one year (Wang 2000a ), combined with SN data at smaller redshifts, so that there are a minimum of 50 SNe Ia per 0.1 redshift interval at any redshift. Although the dispersion of 0.20 magnitudes (intrinsic plus observational) assumed in our simulated data is appropriate for ground-based surveys, we expect our results to apply qualitatively to space based SN pencil beam surveys as well, because our assumed dispersion of 0.20 magnitudes is dominated by intrinsic dispersion (about 0.17 magnitudes).

At the completion of this lengthy numerical study of the feasibility of measuring the time-dependence of the dark energy density from realistic SN Ia data, we became aware of the recent paper by Maor, Brustein, & Steinhardt (2000). They claimed that distance-redshift relations derived from SNe Ia and similar classical measures are poor methods for resolving the time-dependence or measuring the amplitude of the equation of state $`w_X(z)=p_X(z)/\rho _X(z)`$, and consequently no useful information can be obtained about the future fate of the universe. Our work has confirmed the difficulty of measuring the properties of the dark energy from realistic SN Ia data. However, we have found their assessment to be overly pessimistic. Their work indicates that it is impossible to tell whether the equation of state $`w_X(z)`$ varies in time (see also, Barger & Marfatia (2001)), but knowing that the dark energy density $`\rho _X(z)`$ varies in time is sufficient to rule out vacuum energy as dark energy, thus giving support to exotic dark energy models. Our work has shown that one can indeed clearly detect the time-dependence of $`\rho _X(z)`$ using realistic future SN Ia data.

The main problem reported by Moar et al. was the “smearing” effect of the multi-integral relation between the luminosity distance $`d_L(z)`$ and the equation of state $`w_X(z)`$. Instead of $`d_L(z)`$ and $`w_X(z)`$, our analysis uses the time derivative of the comoving distance $`r^{}(z)`$ and the dimensionless dark energy density $`f(z)`$. Thus our results are less affected by the smearing effect. Our method will be refined and made more efficient, and should become quite useful in analyzing future SN data.

In view of our results, it is important that reasonable yet substantial observational efforts are devoted to future SN Ia surveys, for example, a dedicated 4m telescope for a SN pencil beam survey (Wang 2000a ), combined with surveys of nearby SNe Ia. The total cost of such surveys would be modest compared to the great scientific return, the determination of the systematic uncertainties of SNe Ia as cosmological standard candles, and the measurement of the time-dependence in the dark energy density of the universe to constrain fundamental physics.

Acknowledgements

It is a pleasure for us to thank Grant Mathews for helpful discussions, and the referee for useful suggestions. PMG acknowledges support from NASA LTSA grant NAG-9364.