text:
# Searching for continuous gravitational wave signals - The hierarchical Hough transform algorithm -

## 1 Introduction

It is expected that spinning neutron stars, whether isolated or in a binary system, emit continuous gravitational waves. For example, they emit if their shape presents some deviations from axisymmetry. The signals are though expected to be weak: for example from the measured spin-down rate of the period of the Crab or of the Vela pulsars one can place an upper limit to the waves amplitude of $`10^{24}`$. For the Crab pulsar this is a factor of a few higher than the $`99\%`$ confidence limit for a coherent $`10^7`$ s search with initial LIGO. In general, a coherent integration for extended periods of time - of the order of a few months - is necessary in order to accumulate signal to noise sufficient for detection. However, if the phase evolution of the signal is not known in advance, many coherent searches are necessary, each for every signal template resolvable in parameter space. As the time baseline of the coherent integration grows, the number of templates that one must search over grows dramatically and the computational requirements soon become prohibitive. It is then necessary to explore alternative search strategies . Hierarchical algorithms are the best way to explore large parameter space volumes with limited computational resources. In this talk I shall present the hierarchical algorithm that is being developed and implemented at the Albert Einstein Institute to search for continuous gravitational wave signals emitted by isolated sources. Somewhat similar strategies are also being pursued by the Virgo Rome group for the VIRGO experiment and by the UWM (University of Wisconsin, Milwaukee) group for the LIGO experiment. These approaches share many common features and the AEI group is actively collaborating with both the Virgo-Rome group and the UWM group.

## 2 The hierarchical Hough transform search strategy

The search strategy consists of three main steps:

* a first coherent search over data segments lasting about a day - it is foreseen that the total observation time will be of the order of a few (one to four) months.
* An incoherent search over the total observation time.
* Finally, again a coherent search with a time-baseline that is a sizeable fraction of the total observation time.

The latter search is limited to the regions in parameter space that have been selected in the previous stage. This is why hierarchical strategies are so convenient at investigating large parameter volumes: in these circumstances the minimum signal detectable with a given confidence is larger than that one could detect with a priori information on the source parameters. Thus, one does not loose in sensitivity by using in the first stages of the search less sensitive and more computationally affordable methods. With these one can identify which are the regions in parameter space that it is worth investigating with a higher sensitivity .

The procedure has been designed to be performed on a cluster of loosely coupled processors. Each node will run the same code and search a different region in parameter space. As a matter of fact, each node will search a different intrinsic emission frequency band. The signal parameters are: intrinsic emission frequency, position in the sky and a set of spin-down parameters. The maximum order of the spin-down parameters depends on the class of sources that one is investigating. The minimum spin down age of a class of sources embodies all the information necessary to define how many spin down parameters are necessary to describe the frequency evolution of the system.

Distributing the computational load with respect to the intrinsic frequency search parameter is convenient because it allows a natural distribution of the data set among the nodes. This is immediately understood if one reasons in the frequency domain. The maximum Doppler shift that a monochromatic signal will be subject to, with respect to the emission frequency $`f_0`$, is $`10^4f_0`$. This is the order of magnitude of the maximum frequency band of the data around $`f_0`$ that we need to access to search for a signal emitting at that frequency. The techniques that we have developed for our search make it possible to take advantage of this fact and distribute the computational load among the different processors according to the parameter intrinsic frequency and by distributing to the processors different data sets. This reduces some potential throughput problems since every node must access, at any time, only a small fraction of the total data set. For a 4 months search over a band of 500 Hz the total data that has to be accessible by any node is between 400MB and 4GB, depending on the total number of nodes (400MB refers to 100 nodes, 4GB to 10 nodes). HDs of few tens of GBs are completely standard.

## 3 The coherent stage

In this stage longer time-baseline Fourier transforms (DeFTs, hereafter) are produced from a set of shorter time-baseline Fourier transforms (SFTs, hereafter). At the same time a demodulation procedure takes place so that, if a signal having the appropriate parameters were present in the data, its demodulated spectra would present a sharp peak whose amplitude would be no more than $`5\%`$ lower than the total signal power. By appropriate signal parameters we mean here parameters that match those that define the phase evolution used in the demodulation procedure. In this case the frequency at which the peak appears is the intrinsic emission frequency of the signal: the procedure has removed the effects of the Doppler frequency modulation and of the intrinsic frequency evolution (spin-down).

The continuous wavy line in fig. (1) shows the position of a set of spectral peaks, as a function of time. These are the peaks of the signal in a set of SFTs. The circles in the same figure show the positions of the peaks in the demodulated spectra. Clearly the number of peaks has diminished by a factor $`T_{DeFT}/T_{SFT}`$, where $`T_{DeFT}`$ is the time-baseline of the DeFTs and $`T_{SFT}`$ the time-baseline of the SFTs. In this case $`T_{DeFT}/T_{SFT}=21`$.

There are a number of ways to perform this demodulation, which is a well known technique, for example in the field of radio astronomy. We have developed a method that combines only a few data points from each SFT in order to produce the demodulated one. The method has been described in  and . The performance of this method is comparable to that of other methods (such as stroboscopic resampling, ) but it has the advantage of dealing with a much smaller data set. In particular, if the maximum expected shift of the instantaneous frequency during $`T_{SFT}`$ amounts to less than $`\frac{1}{2T_{SFT}}`$, one can only consider about 16 data points per SFT and loose no more than $`5\%`$ in signal power.

The DeFT in the $`b`$-th frequency bin, $`\widehat{x}_b`$, is computed as follows:

$$\widehat{x}_b=\underset{\alpha =0}{\overset{M1}{}}\underset{k=0}{\overset{N1}{}}\stackrel{~}{x}_{\alpha \beta }\frac{1}{N}\underset{j=0}{\overset{N1}{}}e^{2\pi i(\mathrm{\Phi }_{\alpha jb}\frac{jk}{N})}.$$
(1)
With $`\stackrel{~}{x}_{ak}`$ we indicate the $`k`$-th frequency bin of the FFT of the $`\alpha `$-th data set. There are $`M`$ such SFTs, each with $`N`$ frequency bins. $`0<j<N1`$ is the time stamp of the short time-baseline data chunks. $`\mathrm{\Phi }_{\alpha jb}`$ is the expected phase at time $`t=N\alpha +j`$ for an intrinsic emission frequency $`\frac{b}{T_{DeFT}}`$. $`\mathrm{\Phi }`$ parametrically depends on the source location and on the spin-down parameter values of the template source for which one is demodulating. If the phase evolution can be described as linear in $`t`$ during the time duration $`T_{SFT}`$, so that a linear Taylor expansion of $`\mathrm{\Phi }`$ around the mid time point of every SFT time data chunk can be used, for large values of $`N`$ the summation over $`j`$ in eq. (1) can be performed and eq. (1) can be rewritten as

$$\widehat{x}_b=\underset{\alpha =0}{\overset{M1}{}}e^{iy_\alpha }\underset{k=0}{\overset{N1}{}}\stackrel{~}{x}_{\alpha \beta }P_{\alpha k}(b),$$
(2)
with

$`P_{\alpha k}(b)={\displaystyle \frac{\mathrm{sin}x^{}}{x^{}}}i{\displaystyle \frac{1\mathrm{cos}x^{}}{x^{}}}`$ (3)
$`x^{}={\displaystyle \underset{s}{}}f_ssB_{s\alpha }k`$ (4)
$`y_\alpha ={\displaystyle \underset{s}{}}f_sA_{s\alpha }.`$ (5)
In the previous expressions $`f_s`$ indicate the spin-down parameters of the different orders (labeled by the index $`s`$111Note that when $`s=0`$ this quantity denotes the intrinsic frequency, thus has to be computed for the value corresponding to the index $`b`$.), $`A_{s\alpha }`$ and $`B_{s\alpha }`$ are functions that depend on the phase evolution and so their value depends on $`\alpha `$ and on the spin-down parameters. The function $`P_{\alpha k}`$ is peaked around $`x^{}=0`$. Thus in the summation over $`k`$ in eq. (2) one only needs to consider a few values (NTERMS) of $`k`$ around $`k^{}`$ such that $`x^{}(k^{})=0`$. Eq. (2) can then be rewritten as

$$\widehat{x}_b=\underset{\alpha =0}{\overset{M1}{}}e^{iy_\alpha }\underset{k=k^{}\pm NTERMS}{}\stackrel{~}{x}_{\alpha \beta }P_{\alpha k}(b).$$
(6)
As said before, if $`NTERMS`$ is 8 the power loss due to this approximation is less than $`5\%`$. In principle then, to compute a DeFT for a single intrinsic frequency using $`M`$ SFTs, only $`64M`$ bytes of data are needed. If $`T_{DeFT}=1`$ day and $`T_{SFT}=1`$ hour, then this amounts to less than $`1.6`$ kbytes. The point is that in the frequency domain the useful data for a given demodulation are localized and thus one can work with those only and forget about the rest of the data set. This is more involved to implement if one works in the time domain where the information is distributed across the entire data set.

## 4 The incoherent stage

This stage of the search works on spectra and not on Fourier transforms, and this is why it is called “incoherent”. With it, we piece together the information from all the demodulated spectra and select the regions in parameter space that it is worth inspecting coherently with a longer time-baseline than $`T_{DeFT}`$. The way this is done is by building an histogram in parameter space: high clustering in any pixel in such space indicates that there is suspect consistency of the data set with what would be expected if a signal were present in the data having the parameters associated with that pixel.

The transformation that we use to go from the spectra to the histogram is the Hough transform. The basic idea of the Hough transform is very simple yet ingenious , . Suppose that one wanted to find a linear trajectory in an ensemble of points lying in the $`x,y`$ plane and wanted to estimate the parameters of the trajectory. Any such linear trajectory can be described by an equation of the kind $`y=ax+b`$. Uncertainties on $`x`$ and $`y`$ translate in gridding the parameter space $`a`$ and $`b`$ \- suppose there is $`N_a`$ and $`N_b`$ values of $`a`$ and $`b`$ respectively in the resulting grid. For every value of $`a`$ and for every data point $`x,y`$ it is possible to derive a value of $`b`$: the number count in the corresponding $`a,b`$ pixel in the parameter plane should then be enhanced by unity. This is how one would perform the Hough transform in this simple case. To understand what the outcome would be, suppose that the only data points were on a line. In the histogram each data point would contribute a count in the right pixel and a count in other random $`N_a`$ pixels, different from data point to data point. So, after having transformed all the data points, we would find that the average count in any pixel is $`\frac{N_p}{N_b}`$, whereas in the signal pixel it is $`N_p`$. Thus, we would observe a significant clustering in the pixel corresponding to the line parameters. If in the data set there were points that did not belong to the line, i.e. noise, these would raise the number count uniformly in the histogram, with the effect of diminishing the significance of the clustering level in the signal parameter pixel.

In our case the $`x,y`$ plane is a time-frequency plane. The points that appear in it correspond to the frequency bins in the spectra where a significant peak222We have developed techniques to define what a “significant peak” is. We do not want to go into this detail now, though. It suffices for the purposes of this presentation to know that a sensible way to define a significant peak exists. has been identified. If a signal is present in the data, due to the residual modulation, the signal peak will appear displaced from the intrinsic frequency bin with a pattern in time that depends on the mismatch between the trial signal parameters used in the demodulation and the actual signal parameters. Identifying suspect peak patterns that are consistent with what we would expect from a putative source is the goal of this stage of the procedure.

Suppose that the demodulation has been performed for a given sky position $`\widehat{N}`$ and spin-down parameters $`F_s`$ and that a source is located at a different position $`\widehat{n}`$ and that it has a different set of spin-down parameters $`f_s`$. The equation that governs the Hough transform in this general case can be cast in the following form:

$$\nu F_0=\stackrel{}{\xi }(t)(\widehat{n}\widehat{N}),$$
(7)
with:

* $`\nu `$ the measured frequency at time $`t`$ \- namely the frequency of the peak registered at time $`t`$
* $`f_0`$ the intrinsic frequency of the signal at a time $`t_0`$ (usually the beginning of the observation)
* $`F_0=f_0+_s\mathrm{\Delta }f_s[T_{\widehat{N}}(t)T_{\widehat{N}}(t_0)]^s`$
* $`\mathrm{\Delta }f_s=f_sF_s`$, the residual spin-down
* $`T_{\widehat{N}}(t)`$, the time at the solar system barycenter for the sky position defined by the versor $`\widehat{N}`$
* $`\stackrel{}{\xi }(t)=\{F_0+_sF_s[T_{\widehat{N}}(t)T_{\widehat{N}}(t_0)]^s\}\frac{\stackrel{}{v}(t)}{c}+\{_ssF_s[T_{\widehat{N}}(t)T_{\widehat{N}}(t_0)]^{s1}\}\frac{\stackrel{}{x}(t)\stackrel{}{x}(t_0)}{c}`$, a vector that depends on the demodulation parameters, on the position of the detector ($`\stackrel{}{x}`$) and on the relative source-detector velocity ($`\stackrel{}{v}`$).

For a given value of the intrinsic emission frequency and of the residual spin-down parameter $`\mathrm{\Delta }f_s`$, eq. (7) defines, for any given peak registered with a time stamp $`t`$, the locum of points were a source consistent with these parameters would be if it were to produce in the demodulated spectra labeled by the time stamp $`t`$ a peak at $`\nu `$. Such locum of points is derived by solving eq. (7) with respect to $`\widehat{n}`$. This results in a circle on the celestial sphere centered where the vector $`\stackrel{}{\xi }(t)`$ is pointing and with aperture angle $`\varphi `$ such that $`\mathrm{cos}\varphi =\frac{\nu F_0+\stackrel{}{\xi }(t)\widehat{N}}{|\stackrel{}{\xi }(t)|}`$. For every registered peak it does not suffice to enhance the number count in the pixels in the histogram that satisfy eq. (7). In fact, the registered frequencies all carry an uncertainty, due to the finite resolution of the spectrum. As a consequence the locum of points that is consistent with a given set of source parameters and a registered peak is not a circle but it lies inside an annulus, whose width depends on $`F_0`$, on the frequency resolution of the demodulated spectra and on the angle between $`\stackrel{}{\xi }`$ and $`\widehat{N}`$. The borders of such annulus satisfy an equation like eq. (7) with $`\nu \nu +\frac{1}{2T_{DeFT}}`$ and $`\nu \nu \frac{1}{2T_{DeFT}}`$.

Direct solution of eq. (7) involves, once all the variables are set, solving a linear equation in a trigonometric function of the ecliptic longitude, for every possible value of the other coordinate, the latitude. Thus, in this form, it is rather computationally costly. An alternative consists in reverting the argument and instead of transforming peaks into pixels, doing the contrary. In practice this means that given a pixel - defined by $`\widehat{n}`$ \- one computes $`\mathrm{cos}\varphi `$. This is easily done as a scalar product of the cartesian components of $`\widehat{n}`$ and $`\stackrel{}{\xi }`$. Then one uses eq. (7) to find what value of $`\nu `$ this corresponds to and then checks in the time-frequency diagram if a peak was registered at that frequency or not. If so, the number count can be enhanced in that pixel. The gain with respect to the direct solution for $`\widehat{n}`$ is that one does not deal with inverse trigonometric functions. The draw-back is that one is systematically inspecting a large number of pixels that are not part of any annulus, at any given time. A further gain can be obtained in these circumstances by a zooming approach. The concept is still that of the pixel-to-peak approach but one does not start by considering a pixel but rather a larger area. Thus it is possible to associate with this region not just a single frequency value but a range of frequencies. If peaks have been registered that lie in such frequency range then it means that there are annuli intersecting the region, otherwise it means that there are not. If there are no peaks associated with this area the procedure simply does not update the number count in that area and passes on to analyze another region. Note that the frequency interval associated with the area can be degenerate, namely reduce to a single frequency bin. In this case if a peak was registered at such frequency then the number count in all the pixels in the area can be enhanced by one and another region can be analyzed. If one is not in this degenerate case and if there are peaks associated with the region, then the procedure that we have just described is iterated on sub-patches inside the region. This process takes place up to a definite sub-patch size and from there on the pixel-by-pixel approach is used. Preliminary timing results of this method are encouraging when compared to the theoretical optimal performance of , especially for the case where first order spin down parameters are included in the search. Alternative schemes are also under study. These comprise look-up table approaches, peak-to-pixel approaches with a simplified master equation valid for small regions around $`\widehat{N}`$ and hybrid approaches where we try to take advantage of the positive features of the different methods while minimizing the drawbacks.

A different Hough histogram must be produced for every value of the intrinsic frequency and of the spin-down parameter residual. The “brute force” approach to this consists in scanning this parameter set and in producing for each of the parameter values a Hough histogram using all the peaks in the time-frequency diagram. But there exists an alternative way to do this. Let us define a partial Hough map (PHM) as being a Hough histogram generated by using peaks that come from the same spectrum, given a value of the intrinsic frequency and supposing no spin-down residual. The corresponding total Hough histogram is simply the sum of the PHMs coming from the entire set of spectra. Moreover, if every partial Hough map has been constructed by using the same intrinsic frequency, then the resulting total Hough map refers to the same spin down parameter value as the one used in the demodulation (no residual, as said before) and that intrinsic frequency. But note that the effect of a non-zero residual spin-down in the signal is the same as having a time varying intrinsic frequency. This suggests a strategy to re-use PHMs computed for different $`f_0`$s at different times in order to compute a Hough histogram for a non-zero spin down residual. If one has thus a collection of PHMs referring to different intrinsic frequencies and different times one can construct the Hough map for any $`f_0`$ and spin-down residual value simply by adding together at different times PHM corresponding to different $`f_0`$s. In other words, by choosing in a suitable manner the intrinsic frequencies of the different partial maps to be summed, one can investigate the parameter space around the coherent demodulation spin-down parameter value. This approach saves computations because it recognizes that the same sky locations can refer to different values of $`f_0`$ and spin-down residuals, and it avoids having to re-determine such sky locations more than once.

## 5 Unresolved issues

In principle there exist correlations in parameter space that allow considerable computational savings , . Given a certain coherent demodulation time-baseline, for every time label it is possible to identify an optimal gridding of the parameter space. It is though not yet clear how much this gridding changes with time. Understanding this is crucial because the incoherent stage pieces together information from demodulated spectra with different time stamps. For this to work, the spectra must have been demodulated for the same physical parameters. Else, the physical parameter space region to be considered for the incoherent stage is the intersection of the different regions and an appropriate book-keeping must be set up in order to ultimately cover the whole parameter space while avoiding redundand demodulations or excessive storage requirements for the set of demodulated spectra.

With time also the gridding of the parameter space for the incoherent stage changes. This means that the Hough transform pixel size in principle changes depending on what time stamp is attached to the peak that one is considering. In our simulations and operation count estimates we have always used the highest possible resolution, thus often over-resolving the parameter space and thus doing more computations than strictly necessary. The zooming approach adapts well to varying resolution grids, especially if the “PHM strategy” outlined in the previous paragraph is employed. This happens because each PHM refers to a single time-label, thus has its own resolution and the number count in each can be determined using that resolution, thus with the minimum number of operations. In order to ultimately sum together the PHMs, one can construct over-resolved maps from the lower resolution ones and obtain, in the end, the actual physical resolution over the total observing time in the different regions of the sky.

## 6 Development plans

The GEO detector is expected to have its first science data run in ’mid 2001. For that time the search code that performs the search for continuous signals must be in place. It is true that a latency exists between the time the first data are produced and the time the analysis is carried out. In fact, the total set must be acquired before the analysis is initiated. But we plan to use this time not for algorithm development but rather to understand the noise of the detector and refine our data preparation procedures.

We aim at having a hierarchical code working on a single node by the beginning of 2001. As stated in the introduction, the procedure has been designed to run on a set of loosely coupled processors. We are not very concerned, thus, with the generalization of the single-node code to many machines. This is mainly a computational load distribution problem and, as far as we can see, there will be no reason to do this dynamically. Provided that the data preparation is effective, the intrinsic frequency band that each node has to analyze can be decided a priori and the data distributed accordingly once and for all.

## Acknowledgments

The most complicated part of the algorithm is the efficient implementation of the Hough transform. We are developing this in collaboration with S. Frasca, F. Massaioli, C. Palomba and G. Amati at the University of Rome “La Sapienza”. Part of the work presented here is fruit of our joint efforts. We also acknowledge useful discussions with our GEO colleagues A. Vecchio, C. Cutler and B. Sathyaprakash and our colleagues at UWM who are involved in the development of another hierarchical search algorithm, B. Allen, B. Brady, T. Creighton and B. Owen. We would like to thank S. Berukoff who is currently coding the coherent search algorithm, following the LAL (LIGO algorithm library) specifications and helping in the general LAL-compliant set-up of the whole procedure.


format: text

text:
# Optimization of semiconductor quantum devices by evolutionary search

## Abstract

A novel simulation strategy is proposed to search for semiconductor quantum devices which are optimized with respect to required performances. Based on evolutionary programming, a tecnique implementing the paradigm of genetic algorithms to more complex data structures than strings of bits, the proposed algorithm is able to deal with quantum devices with preset non-trivial constraints (transition energies, geometrical requirements, etc.). Therefore, our approach allows for automatic design, thus avoiding costly by-hand optimizations. We demonstrate the advantages of the proposed algorithm by a relevant and non-trivial application, the optimization of a second-harmonic-generation device working under resonance conditions.

In addition to the new fundamental physics which can be investigated in lower-than-3 dimensions , the enormous interest in semiconductor nanostructures has been strongly driven by the possibility of engineering semiconductor-based quantum devices (SBQD’s) ; In particular, among the most successful applications one should mention unipolar (i.e., intersubband) lasers —the so-called quantum-cascade lasers— as well as a variety of semiconductor-based photodetectors .

The design and optimization of SBQD’s is often grounded on an expensive trial-and-error process: after a “promising” structure has been identified, an ad-hoc optimization is performed with respect to a few structural parameters. In this Letter we propose a novel simulation strategy based on evolutionary programming (EP) , to efficiently search in a large parameter space for optimal SBQD’s with respect to required performances: quantum-device optimization by evolutionary search (Q-DOES).

The optimization of SBQD’s is a complex multidimensional maximization problem of a fitness function containing the quality factors of the device; the hardness of the problem grows rapidly with the number of free parameters (number and width of the constituent layers, alloy and dopant concentrations, etc.); furthermore, it may be required that the device works under certain conditions (e.g., with a high optical efficiency in a well-defined energy range) leading to a non-trivial constrained problem. In this respect, genetic algorithms (GAs)  are often invoked for complex optimization tasks; in a typical GA, a pool of possible solutions to the problem —the population— is evolved following a genetic paradigm: solutions belonging to the population —coded as strings of bits, called chromosomes— are allowed to mate in pairs in order to cross their genomes (genetic crossover), while random mutations (bit flips) take place at a given rate to ensure genetic diversity. At the heart of the algorithm is selective pressure: the fittest chromosomes are given more chances to mate and pass their genomes to the next generation.

For complex data structures the bit-string representation is not usually the most natural one. In such cases, the genetic paradigm may be implemented more efficiently by designing context-dependent genetic operators acting on non-conventional chromosomes representing the “natural” data structure of the problem with improved efficency, a strategy often referred to as evolutionary programming . In EP, genetic operators may have a direct physical interpretation, and may be optimally designed to exploit the physical insight of the problem . Finally, in EP inclusion of constraints is usually simpler.

Let us consider a conventional multilayered heterostructure formed by a fixed number, $`N_l`$, of layers of a AxB1-x alloy; layer $`i`$ is characterized by the A concentration, $`x_i`$, and the number of AxB1-x monolayers, $`n_i`$. We map the structure onto a chromosome, each layer being in a one-to-one correspondence with a gene. Hence, we design each chromosome as a vector $`(g_1,\mathrm{},g_{N_l})`$, where $`g_i`$ contains all informations concerning a given layer; in the present case $`g_i=(x_i,n_i)`$ . Often, structural parameters have physical limitations ; such constraints define the solution domain. It is important (and easy) to design crossover and mutation operators which automatically produce offsprings within the solution domain.

We next discuss how to define the fitness function. For a SBQD the quality factors can be written as a functional of the free-carrier wavefunctions $`\psi _n`$ . From these ingredients, we can define a functional $`𝒱[\{\psi _n\}]`$ which measures the performance of the device to be optimized; its evaluation may require the calculation of energy levels, optical matrix elements, scattering rates, etc., depending on the specific device requirements. Non-trivial (non-linear) constraints, such as multiple inter-subband transition energies, may be dealt with by means of a penalty functional which is a measure over the parameter space, $`𝒫[\{\psi _n\}]0`$, going to zero when all constraints are satisfied. In our implementation we always design $`𝒱`$ to be minimum for the fittest structures; more specifically, we minimize

$$𝒱\left[\{\psi _n\}\right]\times \mathrm{exp}\left\{\pm 𝒫\left[\{\psi _n\}\right]/\delta p(t)\right\};$$
(1)
$`\delta p(t)`$ is a parameter which sets the strength of the constraints and is, in general, a function of the simulation time $`t`$; $`+()`$ is taken in the exponential factor if $`𝒱`$ is positive (negative) definite, at least in the region of interest.

The proposed Q-DOES algorithm proceeds as follows. We “evolve” a population of $`N_p`$ solutions, which is randomly initialized; at each simulation step we randomly choose $`\overline{N}_p<N_p`$ solutions which undergo, in pairs, different types of crossover operations (see Tab. 1) with probabilities $`p_c`$, designed so that if parents belong to the desired solution domain, so do the offsprings. Then, several mutation operators are applied (see Tab. 1), with probabilities $`p_m`$, to randomly chosen $`g_i`$’s genes of the offsprings. Each mutation operator has a different purpose: e.g., boundary and uniform mutations ensure ergodicity, while non-uniform mutations perform a local search around the solution coded by the chromosomes which reach the later stages of the simulation. When $`\overline{N}_p`$ offsprings have been generated by crossover and mutation, we convert the alloy concentrations and layer widths coded in their chromosomes —genotype— to the corresponding confinement potentials and space-dependent effective masses —phenotype— and we evaluate the fitness function of the offsprings; in the present implementation the electron wavefunctions $`\psi _n`$ and energy levels $`ϵ_n`$ are calculated within the conventional effective mass approximation (EMA) . Finally, the set of $`N_p+\overline{N}_p`$ chromosomes are ranked according to their fitness-function values, and the fittest $`N_p`$ chromosomes are passed to the next generation.

Let us consider a technologically relevant and highly non-trivial optimization task: the design of a SBQD for second harmonic generation (SHG) working under resonance condition. Schematically, this is an equally spaced three-level structure, where the lowest inter-subband transition energies, $`ϵ_2ϵ_1`$ = $`ϵ_3ϵ_2`$, are in resonance with a pumping radiation, $`\mathrm{}\omega `$, and light emission takes place at $`ϵ_3ϵ_1=2\mathrm{}\omega `$. At resonance, the emitted power is proportional to  $`\mu =|\mu _{12}\mu _{23}\mu _{31}|`$, where $`\mu _{ij}=\psi _i^{}(z)z\psi _j(z)𝑑z`$ are the intersubband dipole matrix elements corresponding to the individual optical transitions. Due to quantum confinement, these quantities can be greatly enhanced in a properly designed asymmetric heterostructure with respect to bulk systems; however, to maximize $`\mu `$ varying the confinement potential is a challenging task, not only because the resonance condition should be preserved in the process, but also because $`\mu `$ is unbounded, and may become very large if resonance conditions are relaxed. Therefore, this example is ideally suited to test the proposed simulation scheme applied to a “hard” problem. This SHG optimization task has been undertaken so far only for extremely simplified potential profiles  and by means of analytical methods .

In order to maximize the SHG emitted power, we set $`𝒱=\mu `$ and we evolve a large population of type-I ($`x0.4`$) AlxGa1-xAs heterostructures according to the Q-DOES algorithm described above. Resonance conditions with the CO2 laser radiation was enforced by setting $`𝒫=[\left(ϵ_2ϵ_1\mathrm{}\omega \right)^2+\left(ϵ_3ϵ_2\mathrm{}\omega \right)^2]^{1/2}`$. The optimized structure produced by the algorithm is shown in Fig. 1: it corresponds to an estimated SHG power $`\mu =4.10\text{nm}^3`$, while the intersubband transition energies fit $`\mathrm{}\omega `$ within numerical accuracy; this compares well with (and is actually $`5\%`$ better than) the largest estimated value 3.91 nm3  obtained so far, to our knowledge, for the same $`\mathrm{}\omega `$, for an idealized continuously graded potential . It is important to stress that the remarkably simple structure shown in Fig. 1 has been obtained with a random initial population, i.e., no guess has been made on the shape of the potential.

The device in Fig. 1 is the best performing one in a small set of optimized structures corresponding to different local minima the alghoritm converged to during a large set of runs, with corresponding fitness function values falling within $`5\%`$ from the best value. The rate of convergence of our algorithm in a typical run is demonstrated in Fig. 2, where we show the values of the fitness function, Eq. (1), and of the SHG intensity, $`\mu `$, of the full population. During the initial stages, the SHG intensity oscillates strongly, exploring a large part of the parameter space, until the search focuses on regions where the energy constraints are approximately satisfied. There, the exponential factor in (1) is $`1`$, and the SHG intensity coincides with the fitness function (in the figure the two vertical axes are shifted for clarity). Eventually, the distribution of the fitness function values gets narrower, as the population becomes degenerate, i.e., all chromosomes become similar; at this stage, we always find the constraints to be satisfied to a very high degree.

Although demonstrated with a simple, yet realistic problem, the algorithm is independent on the detailed structure of the genes, and more complex SBQD’s can be dealt with using basically the same scheme; natural candidates include quantum cascade lasers, mid- and far-infrared photodetectors, cavity resonators, etc. Besides, unlike recently proposed analytical optimization methods , a more accurate (microscopic) description of electronic states may be implemented using, e.g., tight-binding or pseudo-potential Hamiltonians. Obviously, dealing with more complex structures and/or Hamiltonians may increase the computational cost of the algorithm (which was not an issue in the present implementation), which is dominated by the evaluation of the fitness function; however, it should be noted that, in addition to the rapid convergence demonstrated above, the evaluation of the fitness function for different chromosomes is an intrinsically parallel task which may take advantage of a parallel architecture.


text:
# Casimir Effect, Hawking Radiation and Trace Anomaly

## 1 Introduction

The Casimir effect is one of the most interesting manifestations of nontrivial properties of the vacuum state in quantum field theory . Since its first prediction by Casimir in 1948 , this effect has been investigated for various cases of boundary geometries and various types of field. The Casimir effect can be viewed as a polarization of vacuum by boundary conditions. Another type of vacuum polarization arises in the case of external gravitational field. In the static situation, the disturbance of the quantum state induces vacuum energy and stress, but no particles are created. The creation of particles from the vacuum takes place due to the interaction with dynamical external constraints. For example the motion of a single reflecting boundary (mirror) can create particles , the creation of particles by time-dependent external gravitational field is another example of dynamical external constraints.
It has been shown  that particle creation by black hole in four dimension is as a consequence of the Casimir effect for spherical shell. It has been shown that the only existence of the horizon and of the barrier in the effective potential is sufficient to compel the black hole to emit black-body radiation with temperature that exactly coincides with the standard result for Hawking radiation. In , the results for the accelerated-mirror have been used to prove above statement. To see more about relation between moving mirrors and black holes refer to
Another relation between Casimir effect and Schwarzchild black hole thermodynamic is thermodynamic instability. Widom et al  showed the black hole capacity is negative, then an increase in its energy decreases its temperature. They showed that the electrodynamic Casimir effect can also produce thermodynamic instability.
In this paper the Casimir energy for massless scalar field in two dimensional Schwarzschild black hole for two parallel conducting plates with Dirichlet boundary conditions is calculated. The Casimir energy is obtained by imposing general requirements. Calculation of renormalized stress tensor for massless scalar field in two-dimensional Schwarzschild background has been done in  (for a pedagogical review see ). Knowing of Casimir energy in flat space and trace anomaly help us to calculate renormalized stress tensor. To see similar calculation in background of static domain wall refer to . In other situation, besides two previous quantities, we use Hawking radiation, which also has a contribution in stress tensor. Therefore the renormalized stress tensor is not unique and depends on the vacuum under consideration. Our paper is organised as follows: In section 2 general properties of stress tensor are discussed. Then in section 3, vacuum expectation value of stress tensor in two dimensions is obtained. Finally in section 4, we conclude and summarize the results.

## 2 General properties of stress tensor

In semiclassical framework for yielding a sensible theory of back reaction Wald  has developed an axiomatic approach. There one tries to obtain an expression for the renormalized $`T_{\mu \nu }`$ from the properties (axioms) which it must fulfill. The axioms for the renormalized energy momentum tensor are as follow:

1-For off-diagonal elements, standard result should be obtained.

2-In Minkowski space-time, standard result should be obtained.

3-Expectation values of energy momentum are conserved.

4-Causality holds .

5-Energy momentum tensor contains no local curvature tensor depending on derivatives of the metric higher than second order.

Two prescriptions that satisfy the first four axioms can differ at most by a conserved local curvature term. Wald, , showed any prescription for renormalized $`T_{\mu \nu }`$ which is consistent with axioms 1-4 must yield the given trace up to the addition of the trace of conserved local curvature. It must be noted that trace anomalies in stress-tensor, (that is, the nonvanishing $`T_\mu ^\mu `$ for a conformally invariant field after renormalization) are originated from some quantum behavior . In two dimensional space time one can show that a trace-free stress tensor can not be consistent with conservation and causality if particle creation occurs. A trace-free, conserved stress tensor in two dimensions must always remain zero if it is initially zero. One can show that the ”Davies-Fulling-Unruh” formula  for the stress tensor of scalar field which yields an anomalous trace , $`T_\mu ^\mu =\frac{R}{24\pi }`$, is the unique one which is consistent with the above axioms. In four dimensions, just as in two dimensions, a trace-free stress tensor which agree with the formal expression for the matrix elements between orthogonal states can not be compatible with both conservation laws and causality . It must be noted that, as Wald showed , with Hadamard regularization in massless case axiom (5) can not be satisfied unless we introduce a new fundamental length scale for nature. Regarding all these axioms, thus we are able to get an unambiguous prescription for calculation of stress tensor.

## 3 Vacuum expectation values of stress tensor

Our background shows a Schwarzschild black hole with following metric

$$d^2s=(12\frac{m}{r})dt^2+(12\frac{m}{r})^1dr^2+r^2(d\theta ^2+\mathrm{sin}^2d\phi ^2).$$
(1)
Now we reduce dimension of spacetime to two,

$$d^2s=(1\frac{2m}{r})dt^2+(1\frac{2m}{r})^1dr^2.$$
(2)
The metric (2) can be written in conformal form

$$d^2s=\mathrm{\Omega }(r)(dt^2+dr^2),$$
(3)
with

$$\mathrm{\Omega }(r)=1\frac{2m}{r},\frac{dr}{dr^{}}=\mathrm{\Omega }(r).$$
(4)
From now on, our main goal is the determination of a general form of conserved energy momentum tensor with regarding trace anomaly for the metric $`(2)`$. For the non zero Christoffel symbols of the metric $`(2)`$ we have in $`(t,r^{})`$ coordinates:

$$\mathrm{\Gamma }_{tt}^r^{}=\mathrm{\Gamma }_{tr^{}}^t=\mathrm{\Gamma }_{r^{}r^{}}^t=\mathrm{\Gamma }_{r^{}r^{}}^r^{}=\frac{1}{2}\frac{d\mathrm{\Omega }(r)}{dr}.$$
(5)
Then the conservation equation takes the form

$$_r^{}T_t^r^{}+\mathrm{\Gamma }_{tr^{}}^tT_t^r^{}\mathrm{\Gamma }_{tt}^r^{}T_r^{}^t=0,$$
(6)
$$_r^{}T_r^{}^r^{}+\mathrm{\Gamma }_{tr^{}}^tT_r^{}^r^{}\mathrm{\Gamma }_{tr^{}}^tT_t^t=0,$$
(7)
in which

$$T_r^{}^t=T_t^r^{},T_t^t=T_\alpha ^\alpha T_r^{}^r^{}$$
(8)
and $`T_\alpha ^\alpha `$ is anomalous trace in two dimension. Using equations $`(58)`$ it can be shown that

$$\frac{d}{dr}(\mathrm{\Omega }(r)T_t^r^{})=0$$
(9)
and

$$\frac{d}{dr}(\mathrm{\Omega }(r)T^{r_r^{}^{}})=\frac{1}{2}(\frac{d\mathrm{\Omega }(r)}{dr})T_\alpha ^\alpha .$$
(10)
Then equation $`(9)`$ leads to

$$T_t^r^{}=\alpha \mathrm{\Omega }^1(r),$$
(11)
where $`\alpha `$ is a constant of integration. The solution of Eq.(10) may be written in the following form

$$T_r^{}^r^{}(r)=(H(r)+\beta )\mathrm{\Omega }^1(r),$$
(12)
where

$$H(r)=1/2_l^rT_\alpha ^\alpha (r^{})\frac{d}{dr^{}}\mathrm{\Omega }(r^{})𝑑r^{},$$
(13)
with $`l`$ being an arbitrary scale of length and considering

$$T_\alpha ^\alpha =\frac{R}{24\pi }=\frac{m}{6\pi r^3}.$$
(14)
The function $`H(r)`$ produces the non-local contribution of the trace $`T_\alpha ^\alpha (x)`$ to the energy momentum tensor. Given a length scale $`l`$, one can show that,

$$lr_b=2m,$$
(15)
$$H(r)=\frac{m^2}{24\pi }(\frac{1}{16m^4}\frac{1}{r^4}).$$
(16)
Using the equations $`(8),(11)`$and $`(12)`$ it can be shown that in $`(t,r^{})`$ coordinates, we have the following most general form for the stress tensor in our interested background

$$T_\nu ^\mu (r)=\left(\begin{array}{cc}T_\alpha ^\alpha \mathrm{\Omega }(r)^1H(r)& 0\\ 0& \mathrm{\Omega }(r)^1H(r)\end{array}\right)+\mathrm{\Omega }^1\left(\begin{array}{cc}\beta & \alpha \\ \alpha & \beta \end{array}\right).$$
(17)
We will obtain two constants $`\alpha `$ and $`\beta `$ by imposing the second axiom of renormalization scheme. If we put $`m=0`$, then we reach the special case of Minkowski space-time. Here we introduce the state vector $`|C>`$ which is analog for the Boulware state . In the case of existence of boundary the Minkowski limt of $`|C>`$ is not the Minkowski state $`|M>`$. In this limit $`<C|T_{\mu \nu }|C>`$ is non zero and shows the effects of boundary conditions on vacuum of scalar filed. We have chosen the Dirichlet boundary condition $`\varphi (r_1)=\varphi (r_2)=0`$. The standard regularized casimir stress tensor in Minkowskian space-time is as follows

$$reg<0|T_{\mu \nu }|0>=\frac{\pi }{24a^2}\delta _{\mu \nu },$$
(18)
in which $`a`$ is the proper distance between the plates. We can obtain it using mode summation with the help, for example, Abel-Plana summation formula  without introducing the cut off. Comparing to $`(17)`$ we obtain

$$\beta =\frac{\pi }{24a^2}\frac{1}{384\pi m^2},\alpha =0.$$
(19)
Now we consider the Hartle-Hawking state $`|H>`$ . This state is not empty at infinity,even in absence of boundary conditions on quantum filed, corresponds to a thermal distribution of quanta at the Hawking temperature $`T=\frac{1}{8\pi m}`$. In fact,the state $`|H>`$ is related to a black hole in equilibrium with an infinite reservoir of black body radiation. In absence of boundary conditions stress tensor at infinity is equal to

$$<H|T_\nu ^\mu |H>=\frac{\pi T^2}{12}\left(\begin{array}{cc}2& 0\\ 0& 2\end{array}\right).$$
(20)
In the presence of the boundary conditions, (18) have to be add to the above relation. In this case at $`r\mathrm{}`$ we have

$$<H|T_\nu ^\mu |H>=(\frac{1}{384\pi m^2}\frac{\pi }{24a^2})\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),$$
(21)
therefore

$$\beta =\frac{\pi }{24a^2}\frac{1}{192\pi m^2},\alpha =0.$$
(22)
Difference between (22)and(19) is the existence of the bath of thermal radiation at temperature $`T`$.
In order to calculate the contribution from the Hawking evaporation process to the Casimir energy(total vacuum energy), for this special geometry,we introduce the last quantum state which is a convenient candidate for vacuum . This state is called as Unruh state $`|U>`$. In the limit $`r\mathrm{}`$, this state corresponds to the outgoing flux of a black body radiation in the black hole temperature $`T`$. The stress tensor in the limit $`r\mathrm{}`$ in the absence of the boundary conditions is in the form:

$$<U|T_\nu ^\mu |U>=\frac{\pi T^2}{12}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right).$$
(23)
In the presence of the boundary conditions the expression (18) should be added to the relation(23) which is a new stress tensor. Comparing the new stress with (17) one obtains:

$$\alpha =\frac{1}{768\pi m^2},\beta =\frac{\pi }{24a^2}\frac{1}{256\pi m^2}.$$
(24)
Finally stress tensor(17) in Boulware ,Hartle-Hawking and Unruh states is

$$<B|T_\nu ^\mu (r)|B>=\left(\begin{array}{cc}T_\alpha ^\alpha \mathrm{\Omega }^1(r)H(r)& 0\\ 0& \mathrm{\Omega }^1H(r)\end{array}\right)+\mathrm{\Omega }^1(\frac{\pi }{24a^2}+\frac{1}{384\pi m^2})\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(25)
$$<H|T_\nu ^\mu (r)|H>=\left(\begin{array}{cc}T_\alpha ^\alpha \mathrm{\Omega }^1(r)H(r)& 0\\ 0& \mathrm{\Omega }^1H(r)\end{array}\right)+\mathrm{\Omega }^1(\frac{\pi }{24a^2}+\frac{1}{192\pi m^2})\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(26)
$$<U|T_\nu ^\mu (r)|U>=\left(\begin{array}{cc}T_\alpha ^\alpha \mathrm{\Omega }^1(r)H(r)& 0\\ 0& \mathrm{\Omega }^1H(r)\end{array}\right)+\mathrm{\Omega }^1\left(\begin{array}{cc}\frac{\pi }{24a^2}+\frac{1}{256\pi m^2}& \frac{1}{768\pi m^2}\\ \frac{1}{768\pi m^2}& \frac{\pi }{24a^2}\frac{1}{256\pi m^2}\end{array}\right)$$
(27)
. Here the presence of

$$\frac{\pi }{24a^2}\mathrm{\Omega }^1\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(28)
is due to the boundary conditions. Therefore (25),(26) and (27)are separable as follows:

$$<B|T_\nu ^\mu |B>=<B|T_\nu ^{(g)\mu }|B>+<B|T_\nu ^{(b)\mu }|B>.$$
(29)
$$<H|T_\nu ^\mu |H>=<H|T_\nu ^{(g)\mu }|H>+<H|T_\nu ^{(b)\mu }|H>+<H|T_\nu ^{(t)\mu }|H>.$$
(30)
$$<U|T_\nu ^\mu |U>=<U|T_\nu ^{(g)\mu }|U>+<U|T_\nu ^{(b)\mu }|U>+<U|T_\nu ^{(r)\mu }|U>.$$
(31)
where $`<T_\nu ^{(g)\mu }>`$ stands for gravitational, $`<T_\nu ^{(b)\mu }>`$ stands for boundary part, $`<T_\nu ^{(t)\mu }>`$ corresponds to the bath of thermal radiation at temperature $`T`$, and $`<T_\nu ^{(r)\mu }>`$ is Hawking radiation contribution. It should be noted that trace anomaly has contribution in the first term only which comes from the background not boundary effect. However it has contribution to the total Casimir energy momentum tensor. In the regions $`r<r_1`$and $`r>r_2`$the boundary part is zero and only gravitational polarization part is present.
The vacuum boundary part pressures acting on plates are

$$P_b^{(1,2)}=P_b(r=r_{1,2})=<T_1^{(b)1}(r=r_{1,2})>=\mathrm{\Omega }^1(r_{1,2})\frac{\pi }{24a^2},$$
(32)
and have attractive nature. The effective pressure created by other parts in (25-27) are the same from the both sides on the plates, and hence leads to the zero effective force.

## 4 Conclusion

We have found the renormalized energy-momentum tensor for massless scalar field on background of 1+1 dimensional Schwarzschild black hole for two parallel plates with Dirichlet boundary conditions, by making use of general properties of stress tensor only. We propose that if we know the stress tensor for a given boundary in Minkowski space-time, the Casimir effect in gravitational background can be calculated. We have found direct relation between trace anomaly and total Casimir energy. In addition, by considering the Hawking radiation for observer far from black hole, who is the same as Minkowski observer, this radiation contributes to the Casimir effect.
In this paper we have derived three renormalized energy-momentum tensors for our case of study. This is due to selecting three vacuum states for calculation. If we consider Boulware vacuum, stress tensor will have two parts: boundary part and gravitational part. But using Hartle-Hawking and Unruh vacuums will result in another term added to stress tensor, which respectively corresponds to a bath of thermal radiation and Hawking radiation.
It seems that similar results can be obtain for four-dimensional Schwarzschild black hole .

Acknowledgement
I would like to thank Prof A.Saharian for useful discussion, Prof D.V.Vassilevich for the reference  introduce me and Mr.A.Rezakhani for reading manuscript. .


text:
# The Discovery of Argon in Comet C/1995 O1 (Hale-Bopp)

## 1 Introduction

As well preserved, ancient relics of the chemistry present at the formation stage of the outer solar system, comets provide one of the most valuable tools available for understanding the formation processes and conditions extant when the planets formed, $``$4.5 Gyr ago (Mumma et al., 1993).

Among the most long-standing observational goals for understanding cometary origins has been the search for cometary noble gases. Owing to their combination of high volatility and disaffinity to chemical reactions, noble gases provide a key diagnostic to the thermal history of cometary ices. More specifically, the members of the He, Ne, Ar, Kr sequence display successively higher sublimation temperatures; as such, they provide a series of thermometers that can be exploited to constrain the thermal history, and therefore the sites of cometary origins. Although the interpretation of noble gas abundances in cometary comae is complicated by the details of their trapping and release efficiencies in cometary ice (Owen et al., 1991), their detection has nonetheless been highly desired. Unfortunately, however, although significant upper limits revealing He depletions of $``$104 in comet Austin (Stern et al., 1992) and Ne depletions of $``$25 in comet Hale-Bopp (Krasnopolsky et al., 1997) have been obtained, no detection of any cometary noble gas emmanating from within a comet has previously been obtained.

He, Ne, Ar, and Kr each have resonance transitions in the far and extreme ultraviolet. Among the noble gases, argon offers a particularly good combination of comparatively high cosmogonic abundance, moderate sublimation temperature, and good UV resonance fluorescence efficiency; together these properties suggested some time ago (Stern et al., 1992; Mumma et al., 1993) that argon may be the easiest noble gas to detect in comets. We took advantage of the apparition of the unusually active and bright comet C/1995 O1 (Hale-Bopp) in order to conduct a new and more sensitive search for argon in comets.

## 2 The 950–1100 Å Spectrum of Comet Hale-Bopp

With the objectives of (i) obtaining a general survey spectrum of Hale-Bopp in the EUV near its perihelion, and (ii) making a more sensitive search for argon than had previously been attempted for any comet, we launched the EUVS sounding rocket telescope/spectrometer payload on a suborbital mission timed to coincide with the comet’s perihelion. This timing also corresponded closely to the epoch of Hale-Bopp’s peak in activity.

The 184 kg EUVS payload (Slater et al., 1995) consists of a 40-cm diameter grazing incidence telescope, a long slit Rowland circle spectrograph, and its accompanying Ranicon 2-D microchannel plate detector, power system, and telemetry electronics. The telescope is a diamond-turned f/15 Wolter type II grazing incidence design, with a 30-cm aperture (Cash et al., 1989). The primary mirror is Ni coated; the secondary is SiC coated. For the Hale-Bopp flight, EUVS was configured to study the bandpass from 820–1100 Å; the characteristic effective area of the instrument in this bandpass is 0.5 cm2.

EUVS was launched to observe Hale-Bopp on a NASA Black Brant IX sounding rocket from White Sands, New Mexico, at 03:25 UT on 30 Mar 1997. At this time Hale-Bopp was less than 48 hours from perihelion, with a geocentric distance of $`\mathrm{\Delta }`$=1.34 AU, a heliocentric distance of $`R`$=0.915 AU, and a heliocentric radial velocity of $`\dot{R}`$=–1.20 km s-1. The launch occurred in darkness, with the Sun and the comet 116 deg and 76.1 deg from the zenith, respectively. The launch vehicle performed nominally, lifting the instrument to a peak altitude of 314 km. During the flight EUVS remained in the Earth’s shadow at all times; the line of sight to the comet remained above local horizontal. These two factors effectively eliminated telluric dayglow emissions from the spectra that EUVS obtained. The data we describe below were obtained during the 195 sec period when the payload was above 200 km, where telluric absorption of the EUV is also negligible (Paxton and Anderson, 1992).

EUVS successfully obtained spectra to survey Hale-Bopp in the EUV bandpass using both its 130$`\times `$38.6 arcsec (12.5 Å) medium resolution EUVS (MRES) and 130x6.4 arcsec high resolution (2 Å) EUVS (HRES) slits; the centers of these two slits are located 130 arcsec apart. Based on count rate measurements received in real time as telemetry, we elected to maximize SNR at the expense of spectral resolution and place the center of brightness of Hale-Bopp’s coma in the MRES; the adjacent HRES slit was therefore located some 1.27$`\times `$105 km away, along the tail axis of the comet. The EUVS Hale-Bopp dataset contains both strong H Ly $`\beta `$, and the first ever spectroscopic detections of other neutral H and O features in the EUV, cometary O+, and cometary argon (Stern et al., 1998). In this Letter, we concentrate on the 950–1100 Å portion of the EUVS data (see Figure 1) in order to focus on the evidence for, and implications of, Ar in Hale-Bopp; a future publication will discuss the other emissions in the EUVS Hale-Bopp dataset.

The EUVS wavelength scale was established by a least-squares fit to a series of Pt-lamp lines imaged onto the detector before and after the flight; these pre- and post-flight wavelength calibrations concurred to within 3 Å. During data reduction, accurate co-registration of the wavelength scales for the MRES and HRES data was achieved using the centroid of the prominent 1025.7 Å Lyman $`\beta `$ emission line in Hale-Bopp. The effective area of EUVS was calibrated in flight as a function of wavelength using the UV-bright B0.5 IV star, $`ϵ`$ Per. A post-flight laboratory effective area calibration was also obtained, using an O+/Ar resonance source at five wavelengths across the EUVS bandpass. The two calibrations are in good ($`\pm `$35%) agreement.

Using the effective area of the instrument and the solid angle of the slit(s), we converted the Hale-Bopp raw counts spectra to the brightness spectra shown in Figure 1. The brightness error bars shown in Figure 1 were computed at each plotted wavelength from the standard deviation of the mean of the uncertainties due to counting statistics (including the term due to the nominal background subtraction), but neglect the systematic uncertainty in effective area; errors due to ephemeris and pointing geometry uncertainties ($`<<`$1 deg) are completely negligible.

The MRES data in Figure 1 reveals what we believe are the 1048 and 1066 Å ArI features blended together just redward of the far stronger, 1026 Å HI/OI blend. Additional emissions are seen between 975–990 Å and 1080–1097 Å. Although these emissions are clearly statistically significant, their identification is not presently secure; therefore, the analysis of these emissions and the marginally optically thick H Ly $`\beta `$ signal EUVS detected will be discussed in a second publication. Regarding the latter, here we simply note that the Lyman $`\beta `$ feature is blended with a contribution from an OI \[3d3D1,2,0$`>`$2s22p4 3P2,1,0\] triplet at 1025.7, 1027.4, and 1028.7 Å, which, owing to radiative pumping by Ly $`\beta `$, is difficult to precisely model, but which to first order contributes $``$15 R; the telluric background Ly$`\beta `$ brightness is estimated to contribute some 20 R above the background. Sky background spectra obtained while EUVS was maneuvering toward Hale-Bopp contain only a weak instrumental background and H I Ly $`\beta `$ emission generated by a combination of geocoronal and interplanetary medium hydrogen. The fact that no other features are present in EUVS sky background spectra, combined with the observing altitude and geometry discussed above, virtually eliminates the possibility of the features seen at Hale-Bopp as being due to telluric contamination.

Returning now to the main subject of this paper, additional evidence for the Ar identification comes in three forms. First, the width of the MRES feature extending from 1045–1070 Å is consistent with a pair of blended lines 19 Å apart (such as the two Ar lines). Second, the sense of the asymmetry seen in the 1045–1070 Å emission is consistent with the fact that the 1048 Å line’s resonance fluorescence efficiency in sunlight (i.e., its g-factor) is 2.4 times higher than the 1066 Å line’s. Third, as shown in the righthand panel of Figure 1, there is also both clear evidence for a 1048 Å Ar feature (present some 3.2 $`\sigma `$ above the local background) and an hint of the 1066 Å feature in the higher-resolution HRES data, thereby all but eliminating pathological detection cases (e.g., due to flat field effects or count statistics variations) of the Ar features in the main slit.

Since the MRES data has considerably better count levels than the HRES data, we used our line detection software to retrieve individual, background-subtracted gaussians (with widths established from calibration line sources) for the two suspected MRES dataset argon features. This software, which obtains best fit wavelengths and integrated brightnesses gaussians for features detected above a specified background, retrieved brightnesses and statistical error bars of 24$`\pm `$8 R centered at 1049 Å, and 12$`\pm 5`$ R at 1066 Å (1 R=106 ph cm-2 s-1 emitted into 4$`\pi `$ sr). It is worth pointing out that the 2.0$`\pm `$0.6 ratio of these two retrieved brightnesses corresponds well (within our statistics) to the 2.4:1 value predicted by their g-factors, and line formation theory for optically thin emission.111The adopted background level we are using is based on the entire 825–1100 Å MRES spectrum. If we instead adopt an upward sloping background that better corresponds to a slope set by the three lowest pixel brightnesses in the 950–1100 Å MRES spectrum, then one retrieves 22 R at 1049 Å, and 10 R at 1067 Å. Although we do not advocate such a choice of background, we do note that it does more closely reflect the 2.4:1 brightness ratio expected on the basis of the g-factors of the two Ar I features.

## 3 Argon and Its Implications

In order to interpret the Ar emissions, we first converted the signal brightness $`B`$ in Rayleighs into a slit-average column density $`\overline{N}`$ along the line of sight for each Ar line.222On the date of our observations, the projected MRES slit size at the distance of the comet was $``$37,600$`\times `$126,600 km. The Ar photoionization scale length at 0.915 AU ($``$2.7$`\times `$106 km) is large compared to the EUVS slit. The conversion to column densities in the slit was performed according to the optically thin approximation: $`\overline{N_i}=(10^6B_i)/g_i`$, where $`B_i`$ is the brightness of line $`i`$, and $`g_i`$ is the photon fluorescence rate (i.e., ph s-1) of that emission at the comet. We adopted resonance fluorescence efficiencies at 1 AU of 5.3$`\times `$10-8 s-1 and 2.2$`\times `$10-8 s-1 for the 1048 and 1066 Å lines (after adjustment for the variation in solar FUV using the daily F10.7 solar flux index at the time of our flight)333Owing to the fact that we observed Hale-Bopp essentially at perihelion, when its heliocentric radial velocity was consequently small, no Swings effect correction was made.; these g-values (derived from (Meier, 1991; Parker et al., 1998) were then adjusted to reflect Hale-Bopp’s 0.915 AU heliocentric distance.

To derive an Ar production rate ($`Q_{Ar}`$), we constructed a simple model of the coma’s Ar distribution, assuming a spherically symmetric, steady-state radial outflow diverging from a point source. To calculate QAr, we assumed an argon outflow velocity in equilibrium with the H2O outflow, i.e., $`v_i`$=1.25 km s-1 at 0.915 AU. Based on the MRES slit brightnesses (and counting+background brightness errors) of the two argon lines, an error-weighted average Ar production rate of 1.1$`\pm `$0.3$`\times `$1029 s-1, was derived.

To interpret the physical meaning of the EUVS Ar production rate obtained above, we ratio the derived Ar production rates to the O production rate of the comet. H2O is the dominant O-bearing molecule in cometary ices, so we derived the O production rate $`Q_O`$ on the established 1.1$`\times `$1031 s-1 perihelion H2O production rate (Biver et al., 1997; Colom et al., 1997). To account for the oxygen in other O-bearing species (CO, CO2, CHONs, SiOx) when deriving the Ar production ratio \[Ar/O\], we adopt $`Q_{\mathrm{oxygen}}=1.5Q_{\mathrm{H}_2\mathrm{O}}`$ (based on Giotto NMS data adjusted for the higher CO abundance of Hale-Bopp; W. Huebner, pers. comm.). Based on this and the $`Q_{Ar}`$ estimate given above, we thus derive an estimate of \[Ar/O\]=0.0058$`\pm `$0.0017 in Hale-Bopp’s coma. In what follows we assume that this ratio of production rates is indicative of the coma’s Ar/O abundance ratio. The most recent cosmogonic (i.e., solar) \[Ar/O\] abundance ratio gives \[Ar/O\]=0.00372 (Grevesse and Sauval, 1998). Figure 2 presents the derived correspondence between the error-weighted average brightness $`B`$ of the two Ar lines, and the quantity $`Q_{Ar}`$/$`Q_O`$, ratioed to its cosmogonic value.

Although we recognize the difficulty of connecting coma to nuclear abundances with precision because one does not at present know exactly how gases are stored in the nucleus444As a result, inferring absolute abundance ratios with high confidence will likely remain problematic until cryogenic comet nucleus samples can be returned to Earth., we interpret our results at their face value: Hale-Bopp’s coma appears to be enriched in Ar relative to cosmogonic proportions by 1.8$`\pm `$0.5:1.555CO is almost as volatile as argon. It is intriguing that Hale-Bopp’s coma is enriched in Ar (this paper) and also has a notably high CO abundance (Biver et al., 1997). However, we must note that the error estimate given here derives purely from the EUVS MRES counting statistics. Our effective area calibration uncertainty (35%) could reduce this ratio to $``$1.2, but it is equally likely to increase the ratio to 2.4; so too, there is a dispersion in the various Hale-Bopp $`Q_{\mathrm{H}_2\mathrm{O}}`$ measurements at perihelion of 20%, which could reduce this ratio to $``$1.4, but it is equally likely to increase the ratio to 2.2, and a quoted 16% uncertainty in the solar Ar/O ratio (Grevesse and Sauval, 1998). Adding all of these errors in quadrature, we find Hale-Bopp’s coma appears to have an Ar/O ratio relative to cosmogonic proportions that is 1.8$`\pm `$0.96:1.

Because noble gases by their nature do not participate in either the cometary ice chemistry or coma chemistry, noble gases are parent species in the coma, and their abundance in the neutral coma should reflect their bulk abundance within the cometary nucleus. We assume this is the case in what follows.666Though near-surface thermal segregation effects, which are difficult to accurately predict at present, could in principle affect the accuracy of this correspondence.

The lack of depletion of Hale-Bopp’s Ar relative to cosmogonic proportions is interesting in several respects. The first of these concerns the thermal history of Hale-Bopp itself. The fact that Ar is not depleted in Hale-Bopp indicates that the deep interior of Hale-Bopp’s nucleus cannot have ever reached equilibrium temperatures of perhaps 35–40 K, else the Ar would have been lost to sublimation (Bar-Nun and Owen, 1998). Neon, which is even more volatile than argon, sublimates vigorously at temperatures of 16–20 K. The EUVE satellite returned spectra of Hale-Bopp in 1996 showing Ne is $`>`$25 times depleted relative to cosmogonic proportions (Krasnopolsky et al., 1997); this is strong evidence that Hale-Bopp’s deep interior has been warmed above 20 K.

Together, the EUVE satellite and EUVS rocket experiment results trap Hale-Bopp’s bulk internal temperature as having exceeded the Ne sublimation loss range of 16–20 K, but not the Ar sublimation range of approximately 35–40 K. These temperature constraints are consistent with others Hale-Bopp retrieved from H2O ortho-para (Crovisier et al., 1997) and D/H ratios (Blake et al., 1999). We interpret our result as indicating that either the solar nebula was far colder and far richer in Ar than models typically predict (Lunine et al., 2000), or that Hale-Bopp was formed in the cold, Kuiper Belt region (i.e., well beyond the 20–30 AU Uranus-Neptune zone) and was then subsequently ejected to the Oort Cloud (its more recent dynamical home) without ever spending much time in the (warmer) Jupiter-Saturn zone, or both. Although either implication is contrary to “conventional wisdom,” we note that the transport of comets from the Kuiper Belt region (usually after Neptune-induced inward evolution) to the Oort Cloud has been detected in recent dynamical simulations of Oort Cloud formation.(Dones et al., 2000).

This first detection of a native cometary noble gas, argon, whets our appetite for more measurements, in order to assess how typical Hale-Bopp is in its Ar abundance. We also look forward to detection of other noble gases in comets, particularly Ne and Kr, and to the eventual determination of whether the ices within individual comets contain varying proportions of Ar and other chemically inert tracers as a guide to understanding the heterogeneity vs. homogeneity of material within cometary nuclei, and therefore better still understanding their origins.

We thank NASA mission manager Anel Flores, NASA Sounding Rocket branch, NASA Headquarters discipline scientist Henry Brinton, and the entire NASA/contractor team that made these time-critical observations a success. At SwRI, we thank our support staff, including Tom Booker, John McDonald, Clarence McGinnis, David Orth, Bill Tomlinson, and Roy Welch. At the University of Colorado, we thank Gary Kushner and Scott McDonald for assistance with our instrument calibration. Leslie Young provided a careful and insightful reading of this manuscript; David Grinspoon and Walter Huebner provided illuminating discussions. An anonymous referee made a series of helpful comments as well. This work was supported by NAG5-5006.


text:
# Real-Time-RG Analysis of the Dynamics of the Spin-Boson Model

\[

## Abstract

Using a real-time renormalization group method we determine the complete dynamics of the spin-boson model with ohmic dissipation for coupling strengths $`\alpha 0.10.2`$. We calculate the relaxation and dephasing time, the static susceptibility and correlation functions. Our results are consistent with quantum Monte Carlo simulations and the Shiba relation. We present for the first time reliable results for finite cutoff and finite bias in a regime where perturbation theory in $`\alpha `$ or in tunneling breaks down. Furthermore, an unambigious comparism to results from the Kondo model is achieved.

\]

Introduction. The spin-boson model is one of the most fundamental quantum dissipative systems . It plays an important role in describing defect tunneling in solids , quantum tunneling between flux states in a SQUID , electron tunneling between quantum dots , or electron transfer in chemical and biological reactions . The model provides a nontrivial description of dissipation in a quantum system and has attracted much interest due to its simplicity and applicability. It consists of a two-state (spin) system coupled linearly to an infinite-dimensional harmonic oscillator (boson) bath with the Hamiltonian $`H=H_0+H_B+V`$, where $`H_0=\frac{\mathrm{\Delta }}{2}\sigma _x+\frac{ϵ}{2}\sigma _z`$, $`H_B=_q\omega _qa_q^{}a_q`$, and $`V=\frac{\sigma _z}{2}_qg_q(a_q^{}+a_q)`$. Here, $`\sigma _x`$ and $`\sigma _z`$ are the usual Pauli matrices. $`\mathrm{\Delta }`$ is the tunnel matrix element and $`ϵ`$ the energy difference between the two states. $`a_q^{}`$ ($`a_q`$) creates (annihilates) a boson with energy $`\mathrm{}\omega _q`$. The quantity $`g_q`$ represents the coupling strength of the two-state system to the coordinate of the $`q^{th}`$ oscillator. The coupling to the environment is completely defined by the spectral density $`J(\omega )=\pi _qg_q^2\delta (\omega \omega _q)`$, which is usually parametrized by $`J(\omega )=2\pi \alpha \omega ^{n+1}e^{\omega /D}`$. The case $`n=0`$ corresponds to the Ohmic bath, which we want to consider in the following. Therefore the coupling to the bath is now characterized by the coupling strength $`\alpha `$ and the high-energy cutoff $`D`$.

Except for some special parameter values an exact solution of the spin-boson model is not known . Most of the studies are based on perturbative approaches in $`\alpha `$  or $`\mathrm{\Delta }`$ . The latter is known as the noninteracting blip approximation (NIBA) and gives very reliable results at zero bias for the diagonal matrix elements of the reduced density matrix $`p(t)`$ of the two-level system. Recent real-time quantum Monte Carlo (QMC) simulations  provide also reliable information on the nondiagonal elements of $`p(t)`$, but they were only applied to the symmetric case and the correct long-time behaviour has not been checked. Flow equation methods based on infinitesimal unitary transformations  have reproduced the Shiba-relation up to $`\alpha 0.0250.05`$ (with an error of $`310\%`$), but they only addressed spectral properties of the system. Other approaches try to use a mapping of the spin-boson model on the anisotropic Kondo model , and solve the latter exactly using numerical renormalization group (NRG) , Bethe ansatz , or conformal field theory (CFT) . However, NRG and Bethe ansatz provide only spectral properties or dynamics at very short time scales , and CFT has solved so far only the unbiased case $`ϵ=0`$ for the diagonal elements of $`p(t)`$. Furthermore, and most importantly, the mapping on the Kondo model can not be proven rigorously, and the relation of the parameters is not precisely known . It is known that the mapping is incorrect for finite cutoff D, but it is at least established that the scaling behaviour agrees with that of the spin-boson model .

In this paper, we will present for the first time a solution of the complete dynamics of the spin-boson model for $`\alpha 0.10.2`$. We study the diagonal and nondiagonal part of the reduced density matrix $`p(t)`$ starting from an arbitrary nonequilibrium state $`p_0`$. From the asymptotic behaviour we determine the relaxation and dephasing time. Furthermore, we calculate the spin susceptibility and correlation functions. Especially, and in contrast to many other methods, we solve directly the spin-boson model and present results at finite cutoff $`D`$ and finite bias $`ϵ`$. Therefore, our results provide for the first time the possibility for a quantitative and unambigious comparism to results obtained from the mapping of the spin-boson model on the Kondo model. We find that for $`ϵ=0`$ the relaxation parameters agree with those of CFT and in the scaling limit the susceptibility agrees rather well with Bethe-ansatz results, but for the latter deviations occur at finite bias. To demonstrate the reliability of our results, we show the consistency with chromostochastic quantum dynamics (CSQD) , and check the Shiba-relation as well as the scaling behaviour.

To obtain our results, we will use a recently developed real-time renormalization group (RTRG) method . This method has been successfully applied to equilibrium problems , and to the study of nonequilibrium stationary states . Here we will apply it for the first time to the dynamics of the reduced density matrix, and generalize it to the evaluation of correlation functions as well. Since this technique is rather straightforward and easily generalized to other models, this may open a new possibility for the study of various kinds of dissipative quantum systems, like e.g. many-level systems, magnetic nanoparticles interacting with phonons, coupled quantum dots or other kinds of dissipative environments.

Kinetic equation and RTRG approach. The RTRG approach is based on a kinetic equation for the density matrix, see Ref.  for details. We only mention the main steps:

1. The time evolution of the reduced density matrix is written in Liouville space as $`p(t)=\mathrm{Tr}_B\mathrm{exp}(iLt)p_0\rho _B^{eq}`$, where $`\mathrm{Tr}_B`$ denotes the trace over the bath degrees of freedom, and $`L=[H,]=L_0+L_B+L_V`$ is the Liouville operator. The initial density matrix is assumed to decouple into an arbitrary nonequilibrium distribution $`p_0`$ for the two-level system and an equilibrium distribution $`\rho _B^{eq}`$ for the oscillator bath.

2. The propagator $`\mathrm{exp}(iLt)`$ is expanded in the interaction part $`L_V`$, and the trace $`\mathrm{Tr}_B`$ over the bath degrees of freedom is performed by application of Wick’s theorem. In this way one obtains a series of terms where vertices $`G^p`$ of the local system (originating from the $`\sigma _z/2`$ factor in $`V`$) are connected by pair contractions $`\gamma ^{pp^{}}(t)`$ of the bath. Here, $`p=\pm `$ indicates wether the interaction takes place on the forward or the backward propagator. With $`j=_qg_q(a_q^{}+a_q)`$ and $`\gamma (t)=Tr_B\rho _B^{eq}j(t)j`$, we obtain $`\gamma ^{pp^{}}(t)=R(t)+ip^{}S(t)`$ with $`\gamma (t)=R(t)+iS(t)`$ and

$`R(t)`$ $`=`$ $`2\alpha \mathrm{Re}\{(\pi T)^2{\displaystyle \frac{1}{\mathrm{sinh}^2(\pi T(ti/D))}}\},`$ (1)
$`S(t)`$ $`=`$ $`2\alpha \mathrm{Im}\{{\displaystyle \frac{1}{(ti/D)^2}}\},`$ (2)
where we restricted ourselves to the physically relevant situation $`DT`$.

3. From the diagrammatic language one can derive a formally exact kinetic equation for $`p(t)`$

$$\dot{p}(t)+iL_0p(t)=_0^t𝑑t^{}\mathrm{\Sigma }(tt^{})p(t^{}),$$
(3)
where $`\mathrm{\Sigma }(tt^{})`$ is a superoperator acting on $`p(t^{})`$ and is defined by the sum over all irreducible diagrams. In Laplace space we get the formal solution $`p(z)=\mathrm{\Pi }(z)p_0`$ with $`\mathrm{\Pi }(z)=i/(zL_0i\mathrm{\Sigma }(z))`$.

4. The kernel $`\mathrm{\Sigma }(z)`$ is calculated by a renormalization group procedure. Short time scales of $`\gamma (t)`$ are integrated out first by introducing a short-time cutoff $`t_c`$ into the correlation function $`\gamma (t)\gamma (t,t_c)`$. In each renormalization group step, the time scales between $`t_c`$ and $`t_c+dt_c`$ are integrated out, starting from $`t_c=0`$ and ending at $`t_c=\mathrm{}`$. As a consequence, one generates RG equations for $`\mathrm{\Sigma }(z)`$, $`L_0`$, $`G^p`$, and the two boundary vertex operators $`A^p`$ and $`B^p`$ (defined as the rightmost and leftmost vertex of the kernel $`\mathrm{\Sigma }(z)`$). Within the scheme of a perturbative RG analysis, generation of multiple vertex-operators is neglected here, as in Ref. . For the present model, we choose the cutoff dependence either as $`\gamma (t,t_c)=\gamma (t)\mathrm{\Theta }(tt_c)`$ (choice I) or

$$\gamma (t,t_c)=\frac{d}{dt}\left(\stackrel{~}{R}(t)\mathrm{\Theta }(tt_c)\right)+iS(t)\mathrm{\Theta }(tt_c),$$
(4)
with $`R(t)=(d/dt)\stackrel{~}{R}(t)`$ (choice II). It turns out that choice I is better for $`t_c\mathrm{\Delta }`$, whereas choice II is needed for small $`t_c`$ in order to avoid unphysical linear dependencies on the cutoff $`D`$. At $`t_c^0`$ we use a smooth crossover from choice I to choice II. Solving the RG equations numerically gives the kernel $`\mathrm{\Sigma }(z)`$ in Laplace space, and the reduced density matrix can be deduced. The static susceptibility follows from $`\chi _0=(d/dϵ)\mathrm{Tr}_0\sigma _zp_{st}`$, where $`\mathrm{Tr}_0`$ is the trace over the local system, and $`p_{st}`$ denotes the stationary solution $`p_{st}=lim_t\mathrm{}p(t)`$.

5. To calculate correlation functions of the form $`\chi (t)=\frac{1}{2}\mathrm{Tr}\{[\sigma _z(t),\sigma _z]\rho ^{eq}\}`$ we have to generalize the RG procedure. With $`C=G^++G^{}`$ we get in Laplace space

$$\chi (z)=\mathrm{Tr}_0\{\sigma _z\mathrm{\Pi }(z)(C+\mathrm{\Sigma }_C(z))p_{st}\},$$
(5)
where $`\mathrm{\Sigma }_C(z)`$ is analogously defined to $`\mathrm{\Sigma }(z)`$ but contains the “vertex” $`C`$ at any time point. Within the framework of Ref. , the RG equation for $`\mathrm{\Sigma }_C(z)`$ can be easily derived and reads

$`{\displaystyle \frac{d}{dt_c}}\mathrm{\Sigma }_C(z)`$ $`=`$ $`{\displaystyle \underset{pp^{}}{}}{\displaystyle _0^{\mathrm{}}}𝑑t{\displaystyle _0^t}𝑑t^{}{\displaystyle \frac{d}{dt_c}}\gamma ^{pp^{}}(t,t_c)`$ (7)
$`e^{iz(tt^{})}A^pe^{iL_0(tt^{})}Ce^{iL_0t^{}}\overline{B}^p^{},`$
where $`\overline{B}^p=B^p|_{z=0}`$.

Results. Figs. 1 and 2 show the time evolution of $`p(t)`$ for the unbiased and biased case. Initially, the two-level system is prepared in the spin-up state ($`u`$). For $`ϵ=0`$, we achieve a very good agreement with CSQD. Only the real parts of the nondiagonal elements, which correspond to $`\sigma _x`$, exhibit a deviation of approximately $`5\%`$. However, the CSQD can not give an accurate error for $`\sigma _x`$ . The diagonal elements oscillate in time with the frequency $`\mathrm{\Delta }_r\mathrm{\Delta }(\mathrm{\Delta }/D)^{\alpha /(1\alpha )}`$. This agrees with the renormalized tunnel matrix element, which is the characteristic energy scale of the spin boson model . In Fig. 1 we obtain $`\mathrm{\Delta }_r=0.626\mathrm{\Delta }`$ (CFT gives $`\mathrm{\Delta }_r=0.625\mathrm{\Delta }`$). In contrast, the real part of the nondiagonal elements is purely decaying for $`ϵ=0`$ (only for $`ϵ0`$ it also exhibits oscillations). In the scaling limit, i.e. for $`D,\mathrm{\Delta }\mathrm{}`$ such that $`\mathrm{\Delta }_r=\mathrm{const}`$, the diagonal elements are universal, i.e. they only depend on $`\mathrm{\Delta }_rt`$. The nondiagonal elements however have an extra factor of $`\mathrm{\Delta }_r/\mathrm{\Delta }`$ (see inset in Fig. 1) which is consistent with . The long-time behaviour is given by an exponential decay. For $`ϵ=0`$ the decay constants of the diagonal (nondiagonal) elements correspond to the relaxation (dephasing) time $`\tau ^{rel}`$ ($`\tau ^{dep}`$), see insets in Fig. 2. For $`\tau ^{rel}`$ we find good agreement with CFT.For $`ϵ0`$ the diagonal (nondiagonal) elements also contain an incoherent (coherent) decay with $`\tau ^{dep}`$ ($`\tau ^{rel}`$). We note that our zero-temperature results at finite bias and the dephasing time are the first ones presented in the literature.

The static susceptibility for the unbiased and biased case is shown in Figs. 3 and 4. Clearly, the susceptibility depends strongly on the cutoff $`D`$ (see right inset in Fig. 3). However, in the scaling limit the susceptibility is universal. We normalize our results using the zero-temperature value $`\chi _0(T=0)`$. For $`ϵ=0`$, the susceptibility is a monotonous function, whereas a local maximum is obtained for the biased case. A quantitative comparism with recent Bethe ansatz results for the Kondo model shows a very good agreement for the unbiased case (see left inset in Fig. 3). However, for $`ϵ0`$, there are deviations which demonstrates that the commonly used relation of the parameters of the Kondo model to the ones of the spin-boson model has to be taken with care. In the high-temperature limit we obtain the correct $`(1/2T)`$ law, which is independent of $`\alpha `$, $`ϵ`$ and $`D`$ .

The result for the imaginary part of the dynamic susceptibility $`\chi ^{\prime \prime }(\omega )=_{\mathrm{}}^{\mathrm{}}𝑑\omega e^{i\omega t}\frac{1}{2}\mathrm{Tr}\{[\sigma _z(t),\sigma _z]\rho ^{eq}\}`$ is shown in Fig. 5, with $`S(\omega )=\chi ^{\prime \prime }(\omega )/\omega `$. The Shiba-relation provides an exact equation for the low frequency limit of $`\chi ^{\prime \prime }`$ at $`T=0`$

$$\underset{\omega 0}{lim}S(\omega )=2\pi \alpha \chi _0^2.$$
(8)
We tested Eq. (8) for different $`\alpha `$, $`ϵ`$ and $`D`$, see Tab. I. For $`\alpha 0.1`$, we achieved a very good agreement (error smaller than $`5\%`$). Like the static susceptibility, the correlation function depends strongly on the cutoff $`D`$, but in the scaling limit, the normalized quantity $`S(\omega )/S(0)`$ depends only on $`\omega /\omega _{max}`$ (see inset in Fig. 5), where $`\omega _{max}\mathrm{\Delta }_r`$ denotes the frequency where $`S(\omega )`$ is maximum. We note that our results for the correlation function are the first ones presented for the spin-boson model for $`\alpha =0.1`$. NRG results  are very accurate for low frequency but fail for $`\omega \mathrm{\Delta }`$ , flow equation methods have already an error of $`25\%`$ concerning the Shiba-relation , and CSQD does not provide a check of the Shiba-relation and only gives data for the unbiased case $`ϵ=0`$ .

In summary, we investigated the spin-boson model for ohmic dissipation using real-time renormalization group. For coupling parameters $`\alpha 0.10.2`$ we achieved a full solution of the static and dynamical properties. We calculated the real-time evolution of all matrix elements of the reduced density matrix together with the static susceptibility and correlation functions. In contrast to previous works we are not restricted to zero bias or the scaling limit. The restriction in $`\alpha `$ is due to the fact that we neglected double vertex diagrams. Our results show that the RTRG method is a very flexible tool to treat various kinds of dissipative quantum systems.

Acknowledgments. We acknowledge useful discussions with K. Schönhammer and U. Weiss. We also thank J. T. Stockburger and T. A. Costi for their numerical data. This work was supported by the ”Deutsche Forschungsgemeinschaft” as part of ”SFB 345” (M.K.) and “SFB 195” (H.S.).


text:
# The Evolution of Ly𝛼 Absorbing Galaxies

## References

Bland-Hawthorn, J., Freeman, K. C., & Quinn, P. J. 1997, ApJ, 490,143

Chen, H.-W., Lanzetta, K. M., & Fernandez-Soto, A. 2000, ApJ, 533, 120

Davé, R., Hernquist, L., Katz, N., & Weinberg, D. 1999, ApJ, 511, 521

Haardt, F. & Madau, P. 1996, ApJ, 461, 20

Linder, S. M. 1998, ApJ, 495, 637

Linder, S. M. 2000, ApJ, 529, 644

Weymann, R. J. et al. 1998, ApJ, 506, 1