format: text

text:
# Implementing initial state radiation for lepton induced processes in AMEGIC++

## 1 Introduction

Initial State Radiation (ISR) is the most important QED correction to the Born cross section. For instance, in the case of electron positron annihilations at energies beyond the $`Z`$–threshold it leads to the so–called radiative return causing a hugely increased cross section. There, via multiple emission of photons, the electron positron pair tends to reside on the $`Z`$–resonance (returning on it) picking up the propagator in its resonance region. In this way, the correction term becomes several times as large as the Born term at the initial center of mass energy. Evidently, the effect of ISR is of great significance for the interpretation of experimental results. However, since initial state photons tend to be collinear and soft they easily escape detection. In the detector, they tend to leave only imprints by a reduced overall energy and some netto transverse momentum and – possibly – by a boost of the visible final state along the beam axis. Due to this rather indirect method of measuring it, the precise simulation of ISR is crucial for most experimental analyses.

Several strategies exist for this task:

1. The structure function approach

The idea behind this approach is to simulate ISR by a probability density to find an electron with reduced momentum inside an incoming electron. Accordingly, the missing momentum is then attributed to photon radiation. That way the total cross section $`\sigma _{\mathrm{tot}}`$ including ISR is defined as the leading order cross section $`\sigma _\mathrm{B}`$ convoluted with two structure functions $`\mathrm{\Gamma }_{ee}`$:

$$\sigma _{\mathrm{tot}}=𝑑x_1𝑑x_2\mathrm{\Gamma }_{ee}(x_1,s)\mathrm{\Gamma }_{ee}(x_2,s)\sigma _\mathrm{B}(x_1x_2s).$$
(1)
Within the structure functions logarithmically enhanced contributions are exponentiated and thus re-summed to all orders with additional terms covering non–factorisable QED effects up to a finite order $`\alpha ^n`$. Additionally, terms describing weak and strong effects can be taken into account. However, since the structure function is integrated over all photon energies and transverse momenta, it is limited to situations in which only the overall effects of energy reduction and longitudinal boost via ISR are of interest. Arbitrary phase space cuts are difficult to implement. Specific photonic observables like for instance the number of photons above some energy threshold or their transverse momenta are beyond the reach of this approach.
2. The parton shower approach

In this approach, the electron radiation is followed step by step from the comparably low scales of the order of the electron mass connected to the incoming electrons to the high energy scale of the actual hard process. After setting up an electron distribution function at the initial scale, any individual photon emission is governed by the well–known splitting function

$$P_{ee}(z)=\frac{\alpha }{2\pi }\frac{1+z^2}{1z}$$
(2)
leading to a Dokshitzer–Gribov–Lipatov–Altarelli–Parisi (DGLAP) evolution equation for the structure function. Technically, the original structure function is being reproduced statistically at each scale via backward evolution and re-weighting, where the analytic integration is replaced by Monte Carlo integration. Apparently, photons are generated explicitly which allows for instance the application of arbitrary phase space cuts. However, one of the weaknesses of this approach is that it re-sums only the logarithmic enhanced contributions in the DGLAP equation. Even though the radiation of the hardest photon can be corrected with the first order matrix element, it is not clear how to systematically improve to higher perturbative order, i.e. how to make sure that matrix element expressions for the emission of additional photons are recovered for arbitrary processes and arbitrary numbers of photons.
3. The approach by Yennie, Frautschi and Suura (YFS)

Within this approach the explicit generation of photons can be corrected systematically to all orders in the coupling constant. The basic idea here is to introduce some arbitrary infrared cut–off on the photon energies and treat low energy real photons as “un-resolvable”. Their respective contribution will cancel the emerging virtual infrared divergences to all orders leading to finite factorisable terms which can be exponentiated. In this framework correction terms for harder photons can be easily introduced. Since this is the approach we are going to follow we postpone a more thorough discussion.

In this paper we want to give a description of the YFS method implemented in AMEGIC++ . The outline is as follows: In section 2, after briefly reviewing the YFS approach and a short introduction to the nomenclature, we will introduce our new approach to calculate matrix element corrections, where we make use of the ability of AMEGIC++ to calculate (almost) arbitrary tree level matrix elements automatically. In fact, it can be anticipated that our method might allow for the construction of a more generic simulation code for QED initial state radiation in arbitrary processes in the near future. In section 3 we will discuss the actual implementation of our approach and the interplay of ISR with the calculations of matrix elements and cross sections. We justify our method in section 4 by confronting our results for the radiative return to the $`Z`$-pole with experimental data from the LEP experiment. The generality of our program is illustrated by the application to the $`s`$–channel production of Higgs bosons at a possible muon collider. We will conclude with summarizing remarks.

## 2 Method

In this section we want to discuss our new method to include the effect of exact matrix elements in the generation process of initial state photons along the lines of the YFS approach. We will begin with a mini–review of this method suitable for Monte Carlo implementation.

### Mini-review of the YFS approach

The idea underlying the YFS approach is to separate the phase space for the emission of real photons into two regions via a cut–off $`ϵ`$ on the energy fractions, such that photons are coined infrared if their energies $`\omega _\gamma <ϵE_{\mathrm{beam}}`$. The contribution of these infrared photons is then used to cancel the virtual infrared divergences order by order in $`\alpha `$. The remainders of this procedure factorize and can be exponentiated into an universal factor, the YFS form factor:

$`F^{\mathrm{YFS}}(\epsilon )`$ $`=`$ $`\mathrm{exp}\left[2\alpha (B+\stackrel{~}{B})\right]`$ (3)
$`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{\alpha }{\pi }}\left({\displaystyle \frac{1}{2}}\mathrm{ln}{\displaystyle \frac{s}{m^2}}1+{\displaystyle \frac{\pi ^2}{3}}\right)+{\displaystyle \frac{2\alpha }{\pi }}\left(\mathrm{ln}{\displaystyle \frac{s}{m^2}}1\right)\mathrm{ln}\epsilon \right].`$
The virtual part is given by

$$\alpha B=\frac{d^Dk}{(2\pi )^D}\frac{1}{k^2}S(k)$$
(4)
whereas the real contribution reads

$$\alpha \stackrel{~}{B}=\underset{\omega <ϵE}{}\frac{d^{D1}k}{(2\pi )^{D1}2\omega }\stackrel{~}{S}(k).$$
(5)
The functions $`S(k)`$ and $`\stackrel{~}{S}(k)`$ denote the well known universal factorizing “radiation factors” or “eikonal factors” for virtual and real photons for a pair of two external charged lines with four–momenta $`p_1`$ and $`p_2`$

$`S(k)`$ $`=`$ $`{\displaystyle \frac{\alpha }{8\pi }}\left[{\displaystyle \frac{2p_1^\mu k^\mu }{k^22kp_1}}{\displaystyle \frac{2p_2^\mu k^\mu }{k^22kp_2}}\right]^2,`$
$`\stackrel{~}{S}(k)`$ $`=`$ $`{\displaystyle \frac{\alpha }{4\pi }}\left[{\displaystyle \frac{p_2^\mu }{kp_2}}{\displaystyle \frac{p_1^\mu }{kp_1}}\right]^2.`$ (6)
The emission of the visible real photons can then be corrected systematically to all orders in $`\alpha `$ to reproduce exact results as given by the corresponding matrix elements. To be more specific, let us consider as an example $`22`$ processes of the type $`l^+(p_1)l^{}(p_2)f(q_1)\overline{f}(q_2)`$, where the $`p_i`$ and $`q_j`$ label the four–momenta of the external particles. The total cross section including arbitrary numbers of real or virtual photons can be written as

$`\sigma `$ $`=`$ $`F^{\mathrm{YFS}}(\epsilon ){\displaystyle }{\displaystyle \frac{d^3\stackrel{}{q}_1}{q_1^0}}{\displaystyle \frac{d^3\stackrel{}{q}_2}{q_2^0}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\{{\displaystyle \frac{1}{n!}}\left[{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{d^3\stackrel{}{k}_i}{k_i^0}}\stackrel{~}{S}(p_1,p_2,k_i)\mathrm{\Theta }({\displaystyle \frac{2k_i^0}{\sqrt{s}}}\epsilon )\right]`$
$`\delta ^4\left(p_1+p_2q_1q_2{\displaystyle \underset{i=1}{\overset{n}{}}}k_i\right)`$
$`(\beta _0+{\displaystyle \underset{j=1}{\overset{n}{}}}{\displaystyle \frac{\beta _1(k_j)}{\stackrel{~}{S}(k_j)}}+{\displaystyle \underset{j,l=1,jl}{\overset{n}{}}}{\displaystyle \frac{\beta _2(k_j,k_l)}{\stackrel{~}{S}(k_j)\stackrel{~}{S}(k_l)}}+\mathrm{})\}.`$
The YFS form factor $`F^{\mathrm{YFS}}(\epsilon )`$ covers the contribution of factorizing soft real and virtual photons to all orders. The integral over the phase space of the final state particles consists of the integral over the two outgoing momenta $`q_1`$ and $`q_2`$ plus a sum over all possible numbers of photons $`k_i`$ with energy fractions above the resolution threshold $`ϵ`$. This constraint is reflected in the $`\mathrm{\Theta }`$-functions. The conservation of total four–momentum is enforced by the $`\delta `$-function in the second line. Finally, the last line includes the matrix element corrections. The $`\beta _n`$ denote infrared safe combinations of cross sections and eikonal factors with $`n`$ additional real photons as well as all finite contributions of any number of virtual photon loops. In practice, the order of virtual contributions is limited. In $`\beta _n^{(l)}`$ the superscript $`l`$ denotes the total number of virtual and real photons taken into account, i.e. the order of perturbation theory. From there we can read off the number of virtual photons, given by $`(ln)`$.

For a treatment exact up to first order $`\alpha `$, it is sufficient to determine $`\beta _0^{(0)}`$, $`\beta _0^{(1)}`$ and $`\beta _1^{(1)}`$, i.e. the cross sections for $`l^+(p_1)l^{}(p_2)f(q_1)\overline{f}(q_2)`$ up to one loop and for $`l^+(p_1)l^{}(p_2)f(q_1)\overline{f}(q_2)\gamma (k_i)`$ at the tree level. They are defined as

$`\beta _0^{(0)}`$ $`=`$ $`\rho _0^{(0)},\beta _0^{(1)}=\rho _0^{(1)},`$
$`\beta _1^{(1)}(k_i)`$ $`=`$ $`\rho _1^{(1)}(k_i)\stackrel{~}{S}(k_i)\beta _0^{(0)}.`$ (8)
The $`\rho _n^{(l)}`$ are the fully differential cross sections for the hard process $`l^+(p_1)l^{}(p_2)f(q_1)\overline{f}(q_2)`$ with $`n`$ real photons added and up to $`l`$th order in $`\alpha `$, and $`\stackrel{~}{S}`$ denotes the eikonal factor as defined in Eq.(2).

Analyzing the master equation, Eq.(2), we realize that

1. the number of photons that are explicitly produced is not constrained,
2. the resolved photons are allowed to have transverse momentum,
3. the cross section is independent of $`ϵ`$ (in fact, only the number of explicitly produced photons depends on the resolution scale), and
4. that in contrast to the parton shower approach any photon emission is corrected by the appropriate matrix element(s) to the chosen order in $`\alpha `$ including interference effects in multi photon final states.

### The MC algorithm

A Monte Carlo implementation for the YFS approach outlined above was described in  for the first time. The proposed procedure there was

1. to chose the c.m. energy with the form factor $`F^{\mathrm{YFS}}`$,
2. to determine the number of resolvable photons,
3. to supply each of the photons with a four–momentum $`k`$ according to the eikonal factor $`\stackrel{~}{S}(k)`$,
4. to reject all momenta outside their kinematical limits, and
5. to correct the configuration chosen with an additional weight stemming from the exact matrix element.

A more detailed description of the Monte Carlo steps can be found in Appendix A.

### Implementing matrix element corrections

The corrections on the photon distributions and the cross sections due to exact matrix elements can be read off from Eq.(2). Up to first order in $`\alpha `$ the corresponding weight reads:

$`w_{\mathrm{ME}}={\displaystyle \frac{1}{\sigma _B(s^{})}}\left\{\beta _0^{(1)}(p_i,q_j)+{\displaystyle \underset{n=1}{\overset{N_\gamma }{}}}{\displaystyle \frac{\beta _1^{(1)}(p_i,q_j,k_n)}{\stackrel{~}{S}(p_i,k_n)}}\right\},`$ (9)
where again the $`p_i`$ are the incoming momenta, the $`q_j`$ are the outgoing momenta and the $`k_n`$ are the momenta of the explicitly generated photons. In the following we will again restrict ourselves to the description of the process $`e^+e^{}f\overline{f}`$, i.e. a process with two incoming and two outgoing momenta only. However, the extension to more final state particle is straightforward.

An implementation of this weight might seem trivial, however, one serious question remains: How do we calculate matrix elements in situations, where the number of generated photons exceeds the number incorporated in the matrix element? Matrix elements are only unambiguously defined if all incorporated momenta comply with four–momentum conservation. This leads to a problem in cases when some of the generated photon momenta are omitted.

In general, it is solved by applying an unique description of how to evaluate $`\beta `$ functions. Two approaches have been developed in the past to deal with this situation, namely the projection method  and the extrapolation method . In the following we will briefly discuss both methods in order to point out their advantages and disadvantages. Our own solution will be presented afterwards.

The projection method defines for each $`\beta _n`$ function a mapping of the four–momenta from the full phase space with $`N_\gamma `$ photons to the restricted phase space with only $`n`$ photons. Therefore, $`N_\gamma n`$ momenta are completely discarded. The four–momenta of the remaining particles are transformed such that they again comply with four–momentum conservation. In general, this is achieved by a sequence of scalings and boosts applied in individual combinations on the incoming and outgoing particles. For instance, for the leading order term $`\beta _0`$ the projected momenta $`p_i^{},q_j^{}`$ are obtained by boosting the incoming momenta $`p_1`$ and $`p_2`$ into their common c.m. frame and scaling them to fulfill the condition $`(p_1^{}+p_2^{})^2=(q_1+q_2)^2=s^{}`$. The outgoing momenta $`q_1^{}`$ and $`q_2^{}`$ are thus defined by boosting the original $`q_j`$ into their common c.m. frame. A projection for higher order corrections is more involved and will be discussed at the end of this section in the framework of our new approach.

The extrapolation method, on the other hand, defines the matrix element in any phase space with additional photons, i.e. the momenta are kept untouched but the matrix element is rewritten in a suitable set of variables like e.g. angles and energies. Obviously, there is no universal way to construct this transformation and therefore this has to be done for each matrix element separately.

As an example we briefly review the corresponding algorithm in the program KoralZ, that is dealing with processes of the type $`e^+e^{}f\overline{f}`$ with $`fe^{}`$. The leading order term $`\beta _0`$ is written as

$$\beta _0(p_i,q_j;k_1,\mathrm{},k_{N_\gamma })=\frac{1}{2}\left[\frac{d\sigma _0}{d\mathrm{\Omega }}(s^{},\mathrm{cos}\vartheta _1)+\frac{d\sigma _0}{d\mathrm{\Omega }}(s^{},\mathrm{cos}\vartheta _2)\right],$$
(10)
where

$$\mathrm{cos}\vartheta _1=\frac{\stackrel{}{p}_1\stackrel{}{q}_1}{|\stackrel{}{p}_1||\stackrel{}{q}_1|}=\frac{\stackrel{}{p}_1\stackrel{}{q}_2}{|\stackrel{}{p}_1||\stackrel{}{q}_2|}$$
(11)
and

$$\mathrm{cos}\vartheta _2=\frac{\stackrel{}{p}_2\stackrel{}{q}_2}{|\stackrel{}{p}_2||\stackrel{}{q}_2|}=\frac{\stackrel{}{p}_2\stackrel{}{q}_2}{|\stackrel{}{p}_2||\stackrel{}{q}_2|}$$
(12)
are the angles between equal and opposite charged particles, respectively. In the presence of additional photons these angles differ, but $`\beta _0`$ is well calculable for any number of photons, since the variables $`s^{}=(q_1+q_2)^2`$, $`\mathrm{cos}\vartheta _1`$, and $`\mathrm{cos}\vartheta _2`$ are defined in all cases. For the calculation of $`\beta _1`$ the variable $`s=(p_1+p_2)^2`$ and the polar angle of the photon $`\mathrm{cos}\vartheta _\gamma `$ are used in addition to $`\mathrm{cos}\vartheta _1`$, $`\mathrm{cos}\vartheta _2`$ and $`s^{}`$.

However, any approach to define matrix element corrections has to comply with two conditions:

* First, the method has to guarantee that the original result is retained in case the neglected “spectator” photons have total zero four–momentum. In particular, this must hold true in the case of zero spectators.
* Second, the procedure should not give rise to new large logarithms at any order beyond the perturbative order up to which the weight is defined.

At this stage we would like to present a new method to define the $`\beta `$ functions. We aim at a general approach that takes full advantage of the automatic calculation of tree-level matrix elements in AMEGIC++. We do not want to interfere with the automatic calculation as such, which relies on the construction of four–momenta and their explicit conservation. Hence, the method of choice is rather a realization of the projection method.

Remember that the first order infrared finite correction term is

$$\beta _1(p_i,q_j,k_n)=\rho _1^{(1)}(p_i,q_j,k_n)\stackrel{~}{S}(p_i,k_n)\beta _0^{(0)}(p_i,q_j).$$

For this correction term the projection method is well defined as long as there is exactly one additional photon. Then we can apply the PM as introduced for the leading order term $`\beta _0`$ by redefining $`\beta _1`$ as

$$\beta _1(p_i,q_j,k_n)=\rho _1^{(1)}(p_i,q_j,k_n)\stackrel{~}{S}(p_i,k_n)\beta _0^{(0)}(p_i^{},q_j^{}).$$

For more than one photon the arguments of the differential cross section $`\rho _1^{(1)}`$ have to be projected as well. Naively, one might want to apply the same technique as for $`\beta _0`$ with the photon $`k_i`$ included in the outgoing momenta, but it turns out that this is not sufficient to describe events with additional hard photons.

Our solution to this problem is the following:

1. The second term, consisting of the eikonal factor multiplied with the Born term cross section is still treated as described above. Here, we observe that

* the scale of the evaluation of the eikonal is $`s`$, since it is evaluated with the original vectors $`p_1`$ and $`p_2`$. This means that the effect of all other photons is ignored which is yields the correct behavior in order to cancel the eikonal factors in the YFS exponent.
    * The scale of the Born cross section, i.e. the scale of the $`s`$–channel propagator, is $`s^{}`$.
2. In the first term we have to obtain a similar structure. This is achieved by the following steps:

* To calculate the photon part at the scale $`s`$, it is necessary to leave the incoming momenta as well as the photon momentum untouched, but the outgoing momenta have to be modified. The new momenta $`\widehat{q}_1`$ and $`\widehat{q}_2`$ are fixed by the constraint

    $`p_1+p_2`$ $`=`$ $`\widehat{q}_1+\widehat{q}_2+k_n.`$
    Thus, the differential cross section with one photon is given by

    $$\rho _1^{(1)}(p_1,p_2,\widehat{q}_1,\widehat{q}_2,k_n).$$
    (13)
    * In the previous step, the scale of the $`s`$–channel propagator has been changed to $`\widehat{s}=(\widehat{q}_1+\widehat{q}_2)^2`$ which can modify the cross section quite drastically. We cure this by a correction factor

    $$C=\frac{\beta _0(p_1^{},p_2^{},q_1^{},q_2^{})}{\beta _0(\widehat{p}_1,\widehat{p}_2,\widehat{q}_1,\widehat{q}_2)}.$$
    (14)
    Thus, the right peak structure is restored. The new incoming momenta are naturally defined via $`\widehat{s}=(\widehat{p}_1+\widehat{p}_2)^2`$.

Altogether, $`\beta _1`$ is defined as

$$\beta _1(p_i,q_j,k_n)=\frac{\beta _0(p_i^{},q_j^{})}{\beta _0(\widehat{q}_i,\widehat{q}_j)}\rho _1(p_i,\widehat{q}_j,k_n)\stackrel{~}{S}(p_i,k_n)\beta _0(p_i^{},q_j^{}).$$
(15)
## 3 Implementation

In this section we want to discuss some aspects of the interplay of AMEGIC++ with its ISR module. A more thorough description of AMEGIC++ is available in the manual , here we just list the basic steps to calculate cross sections:

1. As a first step, AMEGIC++ sets up the model in which the calculation is going to take place. This includes the definition of particles, the initialization of Feynman rules in terms of vertices and – possibly – the determination of the widths of unstable particles. This last step already relies on the construction and evaluation of corresponding Feynman diagrams and their integration.
2. The construction of Feynman diagrams is performed in three steps:

* The construction of “empty” topologies, i.e. without the specification of lines and vertices,
    * their mapping on “filled” Feynman diagrams with specified propagators and vertices,
    * and, finally, their translation into helicity amplitudes which can be stored in library files.
3. The integration of the diagrams is done with help of an adaptive multichannel method . The channels are constructed by inspection of the kinematical structure of corresponding Feynman diagrams.

Implementing ISR in the framework of YFS affects the last two steps. First, AMEGIC++ has to provide one additional matrix element to calculate the real photon matrix element corrections. The second step is therfore extended as follows:

* An “ISR–Photon” is defined, i.e. an additional particle type is introduced, with all properties of an ordinary photon but exclusively coupling to the charged particles in the initial state.
* One “ISR–Photon” is added to the list of outgoing particles. Correspondingly, “empty” topologies with one additional leg are constructed.
* The Feynman diagrams for this additional process are generated and translated into helicity amplitudes as pointed out above.

During the integration, step 3, the following changes take place:

* The reduced c.m. energy $`\sqrt{s^{}}`$ is determined, and a number of ISR photons is generated along the lines of the algorithm in Appendix A. With this energy a set of incoming momenta is defined in their c.m. frame.
* The momenta of the final state are generated and the leading order cross section is calculated with a multichannel method. So far the effect of ISR is only manifest by the reduced c.m. energy.
* The ISR weight is calculated using the generated final state momenta. Therefore, the matrix element prepared in step 2 with one additional photon as well as the Born matrix element are used according to our method introduced in section 2. The obtained weight is combined with the phase space weight.
* Finally, a boost of all incoming and outgoing momenta is performed in a way that all momenta, i.e. including the ISR photons, comply with four–momentum conservation. Note that this step does not modify the result of the integrated cross section but it is mandatory, when AMEGIC++ is used inside a Monte Carlo generator.

## 4 Results

To demonstrate the quality of our approach we will discuss two examples:

* Processes of the type $`e^+e^{}\gamma ^{},Z^{}f\overline{f}`$ at LEP2 energies including comparison with both experimental data and other generators (KoralZ).
* Processes of the type $`\mu ^+\mu ^{}\gamma ^{},Z^{},h^{}f\overline{f}`$ at higher energies including the effect of propagating Higgs bosons, both from the Standard Model and its minimal supersymmetric extension (MSSM).

### $`e^+e^{}f\overline{f}`$

The effect of ISR on processes of the type $`e^+e^{}f\overline{f}`$ at LEP2 energies is tremendous. Due to the “Radiative Return” huge numbers of events actually reside on the $`Z^0`$ resonance instead of the beam energy.

The effect is twofold: First, the visible energy in the detector tends to be reduced to energies of the order of the $`Z`$ mass accompanied by an increase of the total cross section. Secondly, the forward–backward asymmetry due to the axial coupling of the $`Z`$ boson tends to be shifted correspondingly.

We will highlight this effect in a series of figures. Fig. 2 exhibits the total cross section for $`e^+e^{}\mu ^+\mu ^{}(\gamma )`$ in comparison with data from the DELPHI collaboration . Fig. 2 depicts the combined forward–backward asymmetry in $`\mu `$ and $`\tau `$ pair production at LEP2 compared to data taken by OPAL . Finally, Fig. 3 shows the differential cross section with respect to the cosine of the angle between equally charged incoming and outgoing leptons compared to data from L3 . To guide the eye, we have included the uncorrected Born level prediction with a black line in the first two figures.

In all three cases different cuts on $`s^{}`$ have been applied, and in all cases the agreement of our generator with data is very good. Note that we did not take into account any FSR effects.

In Figs. 5 and 5 the distribution of the effective c.m. energy $`\sqrt{s^{}}`$ of the $`\mu ^+\mu ^{}`$ pair and the distribution of photon energies are shown in comparison with results from KoralZ (second order). We find remarkably good agreement over almost the whole energy range in the $`\sqrt{s^{}}`$ distribution. The energy distribution is dominated by the contribution of the hardest photon which is equally well described. The huge peak at 72.5 GeV corresponds to the radiation of one photon to reach the $`Z`$–resonance. Similarly, the edges at 48.9 GeV and 36.25 GeV correspond to a $`Z`$ return via emissions of two hard photons. In this region our result differs slightly from KoralZ, since the interference effects become more pronounced. However, since the first order matrix element is used for the correction of every explicitly generated photon (not only the hardest) some of the higher order matrix element effects are already included. Note that we have concentrated on ISR effects only and therefore switched off FSR and weak corrections in KoralZ. In case only the first order correction of KoralZ is used, the results would be nearly indistinguishable.

### Muon Collider

As another application we consider a possible muon collider . Its main purpose will be the precise measurement of the properties of the Higgs boson once it is found. The muons are used for two reasons: First, lepton colliders provide comparably pure initial states allowing for spectroscopical measurements (like for instance the exact determination of branching ratios of rare decays like $`h\gamma \gamma `$). Second, muons have a non-negligible coupling to the Higgs boson which in turn will be produced as an $`s`$–channel resonance.

In Fig. 7 and Fig. 7 the total cross section and the forward–backward asymmetry for $`\mu ^+\mu ^{}b\overline{b}`$ are plotted with and without ISR effects assuming a Standard Model Higgs boson with a mass of 115 GeV. The small width of roughly 4 MeV is an experimental challenge, e.g. effects of beam resolution have to be considered. The main effect of ISR is to increase the background due to radiative return events to the $`Z^0`$ peak. However, this background can be reduced via an appropriate $`s^{}`$ cut. Of course, a radiative return to the Higgs peak is also possible, but since the cross section is tiny the effect is small and only visible close to the peak.

Since our method is general, we can discuss further the extension to any scalar boson propagating in the $`s`$–channel. As an example we consider the case of Higgs bosons in the Minimal Supersymmetric extension of the Standard Model (MSSM) . In this model there are two more neutral Higgs bosons which tend to be considerably heavier than the light one. For the following discussion let us constrain ourselves to a supergravity inspired supersymmetry breaking scenario  with parameters to be found in the Appendix B. In this scenario we have a mass of $`m_{H_0}=395.0`$ GeV for the CP–even and and a mass of $`m_{A_0}=394.5`$ GeV for the CP–odd Higgs boson. The total cross section can be found in Fig. 9 and the influence on the forward–backward asymmetry is given in Fig. 9. At this point in the mSUGRA parameter space the individual peaks can not be resolved and the two Higgs bosons appear as one single, comparably broad, resonance only. The effect of ISR is a slight shift of the peak to larger c.m. energies and an increase of background due to radiative return to the lower lying $`Z`$ boson resonance. However, the latter effect can again be reduced by a suitable cut on the effective c.m. energy $`s^{}`$.

In striking contrast to the case of the lighter SM Higgs boson, in this case the radiative return to the peak structure from higher energies is important, see Figs. 11 and 11. Note that the contribution of a photon, a $`Z^0`$ boson, and all three Higgs bosons propagating in the $`s`$–channel as well as all interference terms have been taken into account. The effect of the interplay of the different terms becomes apparent in angular distributions, cf. Fig. 12, where different cuts on the c.m. energy have been applied. The Born term is dominated by the off–peak $`Z`$ boson contribution and exhibits therefore a slight angular dependence. When taking the radiative return to the heavy Higgs bosons peak into account, the distribution is only shifted due to the angular blind scalar coupling structure. The predominating contribution of the return to the $`Z`$ resonance restores the angular distribution already known from Fig. 3.

## 5 Conclusion and Outlook

In this paper we presented an implementation of the successful YFS algorithm into the matrix element generator AMEGIC++. To make full use of the possibilities offered by automatic matrix element generation with help of the helicity method and fully automatic integration of the matrix element we developed a new method to deal with the first order matrix element correction to multi photon emissions. We proved the quality of our method by two examples:

1. For lepton pair production at energies at and above the $`Z`$ boson resonance, we have confronted our simulator with experimental data from the LEP experiment and another ISR generator, KoralZ. In both cases we found encouraging agreement.
2. The application of our approach to the $`s`$–channel production of Higgs bosons at a muon collider has illustrated the generality of our approach. Let us note that extensions to other bosons in the $`s`$–channel with other couplings is straightforward and will be addressed in the near future.
3. Furthermore, the interface of AMEGIC++ with the parton shower implemented in APACIC++ provides a powerful tool for the simulation of hadronic final states including the effect of both ISR and the final state parton shower.

We feel obliged to comment briefly on possible improvements:

* Definitely, effects of Final State Radiation and its interference with ISR are missing and should be implemented. Additionally, some types of one–loop corrections like QCD corrections in the final state, electroweak corrections or box diagrams have not been implemented yet <sup>1</sup><sup>1</sup>1During the time of writing we became aware of , where an ISR calculation for a Higgs boson production at a muon collider was considered. This was based on the structure function approach and already includes FSR, IF–interference and QCD effects. Their findings are consistent with our Monte Carlo results.. Also, an extension to second order corrections should be done to achieve the precision needed for a proper interpretation of the anticipated precision data.
* Apart from $`s`$–channel processes, at $`e^+e^{}`$ colliders processes like $`ZZ`$ production play a crucial role and deserve a similar treatment. In the case of for instance $`WW`$ production an extension of our method is not a trivial task.

## Acknowledgements

We would like to thank K. Hamacher, U. Flagmeyer, U. Müller, A. Behrmann from the Delphi collaboration and P. Ward and D. Ward from the Opal collaboration for pleasant conversation and great help with the experimental data. F.K. gratefully acknowledges financial help of the DAAD. Furthermore, A.S. thanks BMBF for funding.

This work was supported in part by EU Fourth Framework Programme “Training and Mobility of Researchers”, Network “Quantum Chromodynamics and the Deep Structure of Elementary Particles”, contract FMRX-CT98-0194 (DG 12 - MIHT).

## Appendix A The MC algorithm

In this appendix we would like to discuss the Monte Carlo procedure for the generation of ISR photons, further details of the first three steps can be found in .

1. Choice of the reduced c.m. energy squared

First, the reduced c.m. energy squared is determined. This is done according to a “crude distribution” that consists of a fast approximation for the energy dependence of the cross section. This approximation already exhibits the correct peak structure as given by the $`s`$–channel propagators. Additionally, it includes the leading logarithmic photon factor. The corresponding probability density in terms of the energy fraction $`\nu =1s^{}/s`$ reads

$`f(\nu )={\displaystyle \frac{\gamma }{(1\nu )\nu ^{1\gamma }}}\mathrm{\Psi }(\nu )`$ (16)
with

$`\mathrm{\Psi }(\nu )=(1\nu )\sigma _{\mathrm{Born}}(s^{})\widehat{𝒥}(\nu )\text{and}\gamma ={\displaystyle \frac{2\alpha }{\pi }}\left(\mathrm{ln}{\displaystyle \frac{s}{m^2}}1\right)\mathrm{ln}\epsilon .`$ (17)
The energy fraction $`\nu `$ can be determined according to $`f(\nu )`$ by standard Monte Carlo sampling techniques. In the equation above $`\widehat{𝒥}(\nu )`$ is an approximate form of the Jacobian that will be discussed in step 3(b).
2. Choice of number of resolvable photons

Having determined the effective center of mass energy squared $`s^{}=s(1\nu )`$, we proceed by choosing the number of resolvable photons which in turn have to be generated explicitly. In the logarithmic approximation the photon number $`N_\gamma `$ follows a shifted Poisson distribution:

$$f(N_\gamma 1)=e^\mu \frac{\mu ^{N_\gamma 1}}{(N_\gamma 1)!}.$$
(18)
Obviously, its mean value depends on both $`s^{}`$ and the photon resolution parameter $`\epsilon `$. It is given by

$`\mu `$ $`=`$ $`{\displaystyle \frac{d^3\stackrel{}{k}}{k^0}\stackrel{~}{S}(p_1,p_2,k)\mathrm{\Theta }\left(\frac{2k^0}{\sqrt{s}}\epsilon \right)\mathrm{\Theta }\left(\nu \frac{2k^0}{\sqrt{s}}\right)}`$ (19)
$`=`$ $`{\displaystyle \frac{2\alpha }{\pi }}\mathrm{ln}\left({\displaystyle \frac{s}{m_e}}\right)\mathrm{ln}\left({\displaystyle \frac{\nu }{\epsilon }}\right)\mathrm{\Theta }\left(\nu \epsilon \right).`$ (20)
This can be identified in the master formula, Eq.(2), considering the additional constraint on the photon energies (the second $`\mathrm{\Theta }`$–function) not to exceed $`s^{}`$.
3. Construction of four–momenta

The construction of four–momenta for the photons is achieved in two steps:

1. Angles: The photon angles are now individually appointed to each photon with the help of the eikonal $`\stackrel{~}{S}`$. Clearly, $`\stackrel{~}{S}`$ exhibits no dependence on the azimuthal angle $`\phi `$, hence this angle is distributed uniformly. In contrast, $`\stackrel{~}{S}`$ depends strongly on the polar angle $`\vartheta `$, which in turn is chosen according to

    $$\overline{f}(\mathrm{cos}\vartheta )=\frac{1}{(1\beta \mathrm{cos}\vartheta )(1+\beta \mathrm{cos}\vartheta )}\text{with}\beta =\sqrt{1\frac{4m^2}{s}}.$$
    (21)
    Note that terms proportional to the lepton mass have been discarded at this stage. Therefore, a correction weight is introduced for each photon

    $$w_{\mathrm{angle}}(\vartheta )=\frac{f(\vartheta )}{\overline{f}(\vartheta )},$$
    (22)
    where the exact distribution $`f(\vartheta )`$ is given by

    $$f(\vartheta )=\overline{f}(\mathrm{cos}\vartheta )\frac{m^2}{s2m^2}\left(\frac{1}{(1\beta \mathrm{cos}\vartheta )^2}+\frac{1}{(1+\beta \mathrm{cos}\vartheta )^2}\right).$$
    (23)
    Thus, the weight to reproduce the correct mass–dependence reads

    $$w_{\mathrm{mass}}=\underset{i=1}{\overset{N_\gamma }{}}w_{\mathrm{angle}}(\vartheta _i).$$
    (24)
    2. Energies: The next step is the assignment of energies to each of the photons. In case that exactly one photon above the threshold has been generated, the energy of this photon is already completely constrained by $`\nu `$ and the chosen photon angles due to overall four–momentum conservation. In the case that more than one photon have been generated, the additional photon energies are chosen according to $`dk^0/k^0`$, i.e. the eikonal $`\stackrel{~}{S}`$. Then a scaling of all photons is performed in order to achieve four–momentum conservation. Thereby an Jacobian

    $$𝒥(K,\nu )=\frac{1}{2}\left(1+\frac{1}{\sqrt{1A\nu }}\right)$$
    (25)
    with the abbreviations

    $$A=\frac{K^2P^2}{(KP)^2},K=\underset{n=1}{\overset{N_\gamma }{}}k_n$$
    (26)
    is introduced. It had been taken into account already in step 1 as $`\widehat{𝒥}`$ assuming $`A=1`$, since the photon momenta have not been know at this time. A correction weight

    $$w_{\mathrm{jac}}=\frac{𝒥(K,\nu )}{\widehat{𝒥}(\nu )}$$
    (27)
    is introduced to cure this simplification. Furthermore, all photon energies have to be above the limit $`\epsilon \sqrt{s}/2`$ after rescaling, in order to avoid double counting. This leads to an additional weight

    $$w_\epsilon =\mathrm{\Theta }\left(k_n^0\frac{\sqrt{s}}{2}\epsilon \right).$$
    (28)
4. Matrix element correction

The next step is the calculation of the matrix element correction. This is done by introducing an additional weight according to

$$w_{\mathrm{ME}}=\frac{1}{\sigma _0(s^{})}\left\{\beta _0^{(1)}(p_i,q_j)+\underset{n=1}{\overset{N_\gamma }{}}\frac{\beta _1^{(1)}(p_i,q_j,k_n)}{\stackrel{~}{S}(p_i,k_n)}\right\}.$$
(29)
5. Re-weighting

The last step is the calculation of the final weight in order to cure the approximations made during the choice of the reduced c.m. energy (step 1) and the determination of the photon momenta (step 3) as well as to incorporate the matrix element corrections

$$w_{\mathrm{MC}}=w_{\mathrm{mass}}w_{\mathrm{jac}}w_\epsilon w_{\mathrm{ME}}.$$
(30)
This weight can now be used in a rejection method to produce unweighted events.

The ISR simulation in AMEGIC++ differs from this general description mainly in steps 1 and 4.

In step 1, a multichannel method is used to determine $`s^{}`$ to incorporate complicated peak structures. It can therefore be understood as the first part of the multichannel determination of the final state momenta. A library of fast $`f\overline{f}f^{}\overline{f}^{}`$ cross sections for arbitrary scalar and vector $`s`$–channel resonances enables the simulation of many processes in the SM and beyond. Due to the rejection mechanism these cross sections can be also used for multi-jet productions.

In step 4 the evaluation of the matrix element correction is modified as described in section 2, thereby taking full advantage of the features of AMEGIC++.

## Appendix B Program Parameters

The parameters used for the calculation of results presented in section 4 are given in Tab. 2 and Tab. 2.

To include higher order effects in the Yukawa coupling we have used a running b–quark mass according to the well known relation

$$m_\mu =m_{\mu _0}\left(\frac{\alpha _s(\mu )}{\alpha _s\mu _0}\right)^{\frac{\gamma _{m,0}}{\beta _0}}$$
(31)
with

$$\gamma _{m,0}=1,\beta _0=\frac{1}{12}(332n_f).$$
(32)


text:
# Strange attractor simulated on a quantum computer

11institutetext: Laboratoire de Physique Quantique, UMR 5626 du CNRS, Université Paul Sabatier, F-31062 Toulouse Cedex 4, France (Received:)

## Abstract

We show that dissipative classical dynamics converging to a strange attractor can be simulated on a quantum computer. Such quantum computations allow to investigate efficiently the small scale structure of strange attractors, yielding new information inaccessible to classical computers. This opens new possibilities for quantum simulations of various dissipative processes in nature.

Starting from the work of Lorenz lorenz , it has been realized that the dynamics of many various dissipative systems converges to so-called strange attractors ruelle . These objects are characterized by fractal dimensions and chaotic unstable dynamics of individual trajectories (see e.g. ott ; lichtenberg ). They appear in nature in very different contexts, including applications to turbulence and weather forecast lorenz ; ruelle , molecular dynamics hoover , synchronization pikovsky , chaotic chemical reactions belousov , multimode solid state lasers roy  and complex dynamics in ecological systems ecology ; plankton  and physiology glass . The efficient numerical simulation of such dissipative systems can therefore lead to many important practical applications.

Recently, it has been understood that quantum mechanics allows to perform computations in a fundamentally new way (see for review e.g. divi ; ekert ; steane ; preskill ; nielsen ). Indeed, quantum parallelism can enormously accelerate the computation and provide new information inaccessible to classical computers. Well-known examples are Shor’s factorization algorithm shor , which is exponentially faster than any known classical method, and Grover’s search algorithm grover , where the gain is polynomial. Even if important progress has been achieved during the last years, still it is essential to find new areas where quantum processors might give access to new information unreachable classically. Especially interesting are applications to dissipative systems with irreversible dynamics leading to a loss of information.

In this Letter we analyze how classical dissipative dynamics can be simulated on such quantum processors. To this aim, we study a simple deterministic model where dynamics converges to a strange attractor, and show that it can be efficiently simulated on a quantum computer. Even if the dynamics on the attractor is unstable, dissipative and irreversible, a realistic quantum computer divi ; ekert ; steane ; preskill ; zoller ; kane  can simulate it in a reversible way, and, already with 70 qubits, will provide access to new informations inaccessible for modern supercomputers.

To study how a quantum computer can simulate dissipative dynamics leading to a strange attractor, we choose the deterministic map given by:

$$\overline{y}=y/2+x\text{(mod}\text{2)},\overline{x}=x+\overline{y}\text{(mod}\text{1)}.$$
(1)
where $`0.5x<0.5`$, $`1y<1`$ and bars note the new values of variables. The map has one positive $`\lambda _+`$ and one negative $`\lambda _{}`$ Lyapunov exponents ($`\lambda _\pm =\mathrm{ln}[(5\pm \sqrt{17})/2]`$) so that the dynamics converges to a strange attractor with Hausdorff ($`D_H`$) and information ($`D_I`$) dimensions kaplan ; grassberger ; procaccia  $`D_HD_I=1+\lambda _+/|\lambda _{}|1.543`$.

To implement this map for a computer simulation it is necessary to discretize the phase space. We choose the natural discretization in the binary representation of coordinates $`(x,y)`$, so that the dynamics takes place on a regular square lattice with $`N\times 2N`$ points, with $`N=2^{n_q}`$ and $`n_q`$ integer. The division in (1) is realized by shift and truncation of the last binary digit. With this procedure, the dissipation generates a discretized irreversible map, displaying a discretized strange attractor (see Fig.1 top) which approaches the continuous one for large $`N`$. Such a map can be implemented efficiently on a quantum computer. For that, an initial image with $`N_d`$ points is coded in the wave function $`|\mathrm{\Psi }_0>=_{i,j}a_{ij}|x_i>|y_j>|0>|0>`$, where $`a_{i,j}=0`$ or $`1/\sqrt{N_d}`$. Here the two registers $`|x_i>`$ and $`|y_j>`$ with $`n_q`$ and $`n_q+1`$ qubits hold the values of the coordinates $`x`$ and $`y`$ ($`x_i=0.5+i/N,i=0,\mathrm{},N1`$ and $`y_j=1+j/N,j=0,\mathrm{},2N1`$). The third register with $`n_q1`$ qubits is used as workspace for modular additions, and the last one collects the truncated last digits generated by the divisions (“garbage”). We start with the simplified algorithm for which the garbage gets one digit at each map iteration so that $`t`$ iterations need $`t`$ qubits in the fourth register. The size of this register can be significantly reduced using a more refined algorithm we will describe later.

The algorithm starts with the initial state $`|\mathrm{\Psi }_0>`$; first it places the last qubit of $`|y>`$ in the garbage register, and uses $`n_q`$ swap gates to shift the $`|y>`$ qubits and obtain $`y/2`$. Then a modular addition is implemented in the way described in barenco  to add $`x`$ to $`y/2`$. After that, another modular addition adds the second register to the first. In total, this requires $`17n_q10`$ quantum operations using Toffoli, control-not and swap gates in contrast to $`O(2^{2n_q})`$ operations for the classical algorithm. It is interesting to note that the algorithm allows to restore the initial state: inverse map iterations are performed ($`x=\overline{x}\overline{y}`$, $`y/2=\overline{y}x`$) and $`y`$ is restored from $`y/2`$ using a qubit stored in the garbage register. This requires a similar number of operations as the forward iterations. In principle this can be done on a classical computer in $`O(2^{2n_q})`$ operations, but this requires additional exponentially large memory which stores about $`2^{2n_q}t`$ bits for $`N_dN^2`$, contrary to only $`t`$ qubits used by the quantum computer.

Fig.1 shows the dynamics generated by the discretized map (1) simulated on a quantum computer with exact and noisy unitary gates with imprecisions of amplitude $`ϵ`$. Due to the dissipative nature of the map (1), the initial image rapidly converges towards the strange attractor (already $`t=5`$ is enough). Even in the presence of relatively strong noise, the fractal structure of the attractor is well-preserved and the initial image can be reliably recovered after backwards iterations, despite the exponential instability of the classical dynamics (1). The precision of computation can be quantitatively characterized through the fidelity $`f`$ defined at a given moment of time as the projection of the quantum state in presence of gate imperfections on the exact state without imperfections. The global properties of the initial image can be recovered even at relatively low fidelity values.

Even if the quantum algorithm performs one map iteration only in $`O(n_q)`$ operations, it is important to take into account the measurement procedure that allows to extract efficiently the information coded in the wave function. Indeed the number of points in Fig.1 grows exponentially with $`n_q`$ and an exponential number of measurements is required to obtain the full density distribution. However, certain characteristics can be extracted in a polynomial number of measurements, providing new information inaccessible for classical computation. An example of such a quantity is the spectrum of phase space correlation functions, defined as $`C(t,k_{x,y})=_{x_0,y_0}\mathrm{exp}(2i\pi (x(t,x_0,y_0)+y(t,x_0,y_0)))\mathrm{exp}(2i\pi (k_xx_0+k_yy_0))`$, where the sum runs over the points $`(x_0,y_0)`$ of the initial distribution, and $`(x(t,x_0,y_0),y(t,x_0,y_0))`$ is the position of $`(x_0,y_0)`$ at time $`t`$. Such correlation functions have been studied for chaotic systems, where they determine various kinetic coefficients, for example the diffusion rate (see e.g. lichtenberg  p.328). Due to chaos, the function $`C(t,k_{x,y})`$ has significant values at exponentially high harmonics $`k_{x,y}\mathrm{exp}(|\lambda _{}|t)`$ which rapidly reaches harmonics of order $`N`$. In the theory of classical chaotic dynamics it is well-known that the information about such harmonics is very hard to access, since exponentially small scales should be explored, which can be done only with exponentially many trajectories lichtenberg ; peres . On the contrary, the quantum computation of $`C(t,k_{x,y})`$ can be done efficiently. For that, one makes $`t`$ iterations of (1), and creates the state $`a_{x,y}\mathrm{exp}(2i\pi (x+y))|x>|y>|0>|g(t)>`$. The preparation of this state is easily done by applying $`2n_q+1`$ one-qubit rotations to the first two registers. Then the garbage $`g(t)`$ is erased by iterating the map backwards $`t`$ times, that at the same time returns the coefficients $`a_{x,y}`$ to their original values. This creates the state $`a_{x_0,y_0}\mathrm{exp}(2i\pi (x(t,x_0,y_0)+y(t,x_0,y_0))|x_0>|y_0>|0>|0>`$, keeping phases unchanged. The whole procedure is sometimes called “phase kickback”. After that the application of a two-dimensional quantum Fourier transform ekert  yields in $`O(n_q^2)`$ operations the state $`C(t,k_{x,y})|k_x>|k_y>|0>|0>`$. A polynomial number of measurements yields the principal peaks, or enables to obtain a coarse-grained image of the spectral density $`|C(t,k_{x,y})|^2`$ in the Fourier space. Indeed, independently of $`n_q`$, one can measure the first $`n_f`$ and $`n_f+1`$ qubits of the $`x`$ and $`y`$ registers respectively, that gives integrated probability inside $`2^{2n_f+1}`$ cells. Fig.2 displays the spectral density for the case of Fig.1. It shows that new information about the coarse-grained spectral density can be obtained efficiently. Indeed, patterns are clearly present in Fig.2 and they vary irregularly with $`n_q`$ (compare Fig.2 middle left and bottom left). This confirms the nontrivial nature of information provided by the coarse-grained density $`|C(t,k_{x,y})|^2`$. Although the spectral density is more sensitive to noise than the distribution in Fig.1 still the patterns remain well-defined even in the presence of relatively strong errors (see Fig.2).

It is important to stress that even modern supercomputers are unable to find the properties of the spectral density for $`n_q20`$. Indeed, as is shown in Fig.3, a classical Monte Carlo algorithm requires an exponentially large number of trajectories $`M`$ ($`M=O(2^{2n_q})`$) to obtain the coarse-grained spectral density at fixed $`n_f`$ with fixed accuracy. In contrast, the quantum computation requires a number of measurements $`M`$ independent of $`n_q`$ (each measurement is done after $`t`$ map iterations and one Fourier transform which needs $`O(n_q^2)`$ quantum gates).

To study the effect of noisy gates in a more quantitative way, we show in Fig.4 the dependence of the fidelity $`f(t)`$ of the quantum computation of spectral amplitudes $`C(t,k_{x,y})`$ on the noise amplitude $`ϵ`$ and total number of gates $`n_g`$ applied ($`n_g=t(44n_q14)+(n_q+2)^22`$). The data show the global scaling law $`1f(t)ϵ^2n_g/14`$, valid for moderate $`ϵ`$. The physical origin of this scaling law is related to the fact that for randomly fluctuating unitary gates the loss of probability from the exact state is of the order of $`ϵ^2`$ for each gate operation. This law determines a time scale $`t_f1/(6ϵ^2n_q)`$ up to which a reliable quantum computation is possible ($`f0.5`$). Beyond $`t_f`$ the decoherence destroys the accuracy of the quantum computation and the results become strongly distorted, see e.g. Fig.2 bottom right. We note that a similar time scale appears in quantum computation of Shor’s algorithm on a realistic quantum computer zurek . The computation beyond the scale $`t_f`$ is possible but requires the application of quantum error-correcting codes, at the cost of additional qubits and gates (see e.g. steane ; preskill  and refs. therein). Without error correction, it is still possible to improve the fidelity of the final state by measuring the third and fourth registers. Indeed for exact computation they are at zero after the backwards iterations while with noisy gates the error probability grows as $`W_g3.5ϵ^2t`$ (Fig.4 inset). The measurements of these registers allow to select the correct states and increase the fidelity by a factor $`1/(1W_g)`$, e.g. for the case of Fig.1 bottom this procedure gives $`f=0.22`$.

The above algorithm is optimal for not very large times $`t`$. If one is interested in simulating the dynamics on the attractor for large $`t`$, then the size of the garbage register can be significantly reduced. Indeed, at any $`t`$, copying the results $`x(t),y(t)`$ in two additional registers and reversing the sequence of gates allows to erase the garbage and reproduce $`x_0,y_0`$. This procedure can be done recursively following the strategy of “reversible pebble game”, the description of which can be found in preskill . In the simplest version, $`n_t`$ qubits in the garbage register (plus the two additional registers for $`x_0,y_0`$) allow to perform map iterations up to $`t2^{n_t}`$. This gives only a polynomial increase in the number of elementary quantum operations, being proportional to $`n_g^{1.58}`$. The procedure becomes cost-effective for $`tn_q`$.

The algorithm described above can be generalized to other dissipative maps. For example, modular multiplications can be performed in $`O(n_q^2)`$ operations, as described in barenco . This allows to simulate efficiently the map (1) with $`x^2`$ term in first equation and also the Hénon attractor $`\overline{x}=y+1ax^2,\overline{y}=bx`$ lichtenberg . Such an algorithm can also be adapted to perform the finite-step integration of the Lorenz system $`\dot{x}=\sigma (xy),\dot{y}=xz+rxy,\dot{z}=xybz`$ lorenz . However the simulation of the dissipative dynamics of these models requires more qubits than for (1).

Thus on the example of the map (1), we have shown that a quantum computer can efficiently simulate dissipative irreversible dynamics. A quantum processor with 70 qubits will be able to provide new information about small-scale structures of strange attractors inaccessible to modern supercomputers.

We thank CalMiP in Toulouse and IDRIS in Orsay for access to their supercomputers which were used to simulate quantum computations. This work was supported in part by the EC RTN contract HPRN-CT-2000-0156 and also by the NSA and ARDA under ARO contract No. DAAD19-01-1-0553.


text:
# On the stable rank of algebras of operator fields over an n-cube

## 1. Introduction

Rieffel  introduced the notion of *stable rank* for $`C^{}`$-algebras as the noncommutative version of complex dimension of ordinary topological spaces. It turns out that the stable rank of a unital $`C^{}`$-algebra is the same as its Bass stable rank (see ). The purpose of this paper is to study stable rank for continuous field $`C^{}`$-algebras. Our main result is

###### Theorem 1.1.

Let $`𝒜`$ be a unital maximal full algebra of operator fields with base space $`[0,1]^k`$ and fibre algebras $`\{𝒜_t\}_{t[0,1]^k}`$. Then the stable rank of $`𝒜`$ satisfies the following inequality:

$$sr(𝒜)sup_{t[0,1]^k}sr(C([0,1]^k)𝒜_t).$$

One of the key ingredients of the proof of the above result is Nistor’s notion of *absolute connected stable rank* (see ). Recall that for a unital $`C^{}`$-algebra $`𝒜`$, the absoluted connected stable rank of $`𝒜`$ is numerically the same as the stable rank of the tensor product $`C[0,1]𝒜`$.

As an application of our theorem, we will compute the stable ranks of the universal $`C^{}`$-algebras of the (possibly higher rank) discrete Heisenberg groups. Recall that for a positive integer $`n`$, the discrete Heisenberg group of rank $`2n+1`$ is the group of all $`2n+1`$ by $`2n+1`$ upper triangular matrices with integer entries, with ones on the diagonal, and zero entries on all but the first row, last column and the diagonal. The rank $`2n+1`$ discrete Heisenberg group is naturally a lattice subgroup of the $`2n+1`$-dimensional Heisenberg Lie group.

The stable ranks of the universal $`C^{}`$-algebras of various type I Lie groups have been extensively studied (see , , ). In particular, we point out that for a simply connected nilpotent Lie group $`G`$, the stable rank of the universal $`C^{}`$-algebra $`C^{}(G)`$ of $`G`$ has been computed by Sudo and Takai (see ). Roughly speaking, they showed that the stable rank of $`C^{}(G)`$ is controlled by the ordinary topological dimension of the space of one-dimensional representations of $`G`$.

Recently, the stable ranks of a class of non-type I solvable Lie groups (which include the Mautner group) have also been computed (see ).

Our computations for the discrete Heisenberg groups constitute a class of interesting nontrivial examples of the stable rank of the universal $`C^{}`$-algebra of a non-type I discrete group.

In a later paper, we will apply the techniques developped in this paper to compute the stable ranks of arbitrary finitely generated, torsion-free two-step nilpotent groups.

A general reference for stable rank of $`C^{}`$-algebras is . General references for algebras of operator fields are ,  and .

In what follows, for a $`C^{}`$-algebra $`𝒜`$, the notation “$`sr(𝒜)`$” will always mean the stable rank of $`𝒜`$. If, in addition, $`𝒜`$ is unital and $`N`$ is a positive integer, then $`Lg_N(𝒜)`$ is the set of $`N`$-tuples $`(a_1,a_2,\mathrm{},a_N)`$ in $`𝒜^N`$ such that $`_{i=1}^N(a_i)^{}a_i`$ is an invertible element in $`𝒜`$.

## 2. Main results

###### Lemma 2.1.

Let $`𝒜`$ be a unital maximal full algebra of operator fields with base space $`[0,1]`$ and fibre algebras $`\{𝒜_t\}_t`$. Then the stable rank of $`𝒜`$ satisfies the inequality

$$sr(𝒜)sup_{t[0,1]}sr(C[0,1]𝒜_t).$$

###### Proof.

In this proof, we will always let “$``$” denote the continuity structure for the continuous field decomposition of $`𝒜`$ in the hypothesis.

Suppose that $`M=sup_{t[0,1]}sr(C[0,1]𝒜_t)`$ is a finite number. Let an $`M`$-tuple $`(a_1,a_2,\mathrm{},a_M)𝒜^M`$ and a positive real number $`ϵ>0`$ be given. Since $`[0,1]`$ is compact, let $`I_1,I_2,\mathrm{},I_N`$ be a finite set of open intervals which cover $`[0,1]`$ and for each $`i=1,2,\mathrm{},N`$, let $`(f_{i,1},f_{i,2},\mathrm{},f_{i,M})`$ be an $`M`$-tuple in $`𝒜^M`$ such that

1. $`_{j=1}^Mf_{i,j}(t)^{}f_{i,j}(t)`$ is an invertible element of $`𝒜_t`$ for $`tI_i`$. Here the operator field $`tf_{i,j}(t)`$ for $`t[0,1]`$, is the representation of $`f_{i,j}`$ as a continuous operator field in $`𝒜`$. And
2. $`f_{i,j}a_j<ϵ`$, for all $`i`$, $`j`$.

We can choose such intervals $`I_i`$ and such elements $`f_{i,j}`$, since our hypothesis implies that the stable rank of each fibre $`𝒜_t`$ is less than or equal to $`M`$, and by the existence and continuity of operator fields in a full algebra of operator fields (see the definition of full algebra of operator fields in ,  or ). Also, one needs to use the fact that any element of a unital $`C^{}`$-algebra that is sufficiently close to the unit is invertible. Finally, to make $`f_{i,j}`$ uniformly within $`ϵ`$ of $`a_j`$ (for all $`i`$, $`j`$), one needs the maximality of the algebra of operator fields $`𝒜`$ (see  Proposition 1 and  Theorem 1.1).

For simplicity, we will assume that there are only two intervals ($`I_1`$ and $`I_2`$) and only two $`M`$-tuples ( $`(f_{1,1},f_{1,2},\mathrm{},f_{1,M})`$ and $`(f_{2,1},f_{2,2},\mathrm{},f_{2,M})`$ ). We may additionally assume that neither interval is contained in the closure of the other, and that their intersection is a continuum. Our goal is to “connect” the two $`M`$-tuples over the intersection $`\overline{I_1I_2}`$ to get an element of $`Lg_M(𝒜)`$ which approximates $`(a_1,a_2,\mathrm{},a_M)`$ within $`ϵ`$. The proof for more than two intervals is an iteration of this procedure (after appropriate contraction, removal or addition of intervals…).

Our procedure for “connecting” the two $`M`$-tuples over $`\overline{I_1I_2}`$ will involve constructing a sequence of strictly increasing points $`\{t_n\}_{n=1}^{\mathrm{}}`$ in $`I_1I_2`$, and constructing sequences of operator fields $`\{\alpha _j^n\}_{n=1}^{\mathrm{}}`$ ($`j=1,2,\mathrm{},M`$) satisfying:

1. $`\alpha _j^n`$ is a continuous operator field in $`𝒜(n)`$, for all $`j`$, $`n`$. Here $`𝒜(n)`$ is the unital maximal full algebra of operator fields gotten by restricting $`𝒜`$ to the interval $`[t_n,t_{n+1}]`$. (In particular, the continuity structure for $`𝒜(n)`$ is gotten by taking the restriction to $`[t_n,t_{n+1}]`$ of all the fields in $``$).
2. $`\alpha _j^n`$ is within $`ϵ`$ of the restricion of $`a_j`$ to $`[t_n,t_{n+1}]`$, for all $`j`$, $`n`$.
3. $`\alpha _j^1(t_1)=f_{1,j}(t_1)`$ and $`\alpha _j^{n+1}(t_{n+1})=\alpha _j^n(t_{n+1})`$, for all $`j`$, $`n`$.
4. $`\alpha _j^{n+1}`$ is within $`ϵ/2^n`$ of the restriction of $`f_{2,j}`$ to $`[t_{n+1},t_{n+2}]`$, for all $`j`$, $`n`$. And
5. $`(\alpha _1^n,\alpha _2^n,\mathrm{},\alpha _M^n)`$ is in $`Lg_M(𝒜(n))`$; or equivalently, $`_{j=1}^M\alpha _j^n(t)^{}\alpha _j^n(t)`$ is an invertible element of $`𝒜_t`$ for all $`t[t_n,t_{n+1}]`$, for all $`n`$ (the proof of equivalence is a small spectral theory argument which uses the fact that $`[t_n,t_{n+1}]`$ is compact).

Henceforth, we will let “(\*)” denote conditions (1) - (5).

By  Lemma 2.4 and our hypothesis for $`M`$, we have that $`Lg_M(𝒜_t)\{(b_1,b_2,\mathrm{},b_M)𝒜_t^M:a_i(t)b_i<ϵ\text{ }i\}`$ is a nonempty connected open set for all $`t[0,1]`$. Hence, fixing $`t_1I_1I_2`$, there are continuous paths $`\gamma _j:[0,1]Lg_M(𝒜_{t_1})`$, $`j=1,2,\mathrm{},M`$, such that a) $`\gamma _j(0)=f_{1,j}(t_1)`$ and $`\gamma _j(1)=f_{2,j}(t_1)`$ for all $`j`$, and b) $`\gamma _j(s)`$ is within $`ϵ`$ of $`a_j(t_1)`$ for all $`s[0,1]`$ and for all $`j`$.

Note that since $`[0,1]`$ is compact, $`(\gamma _1,\gamma _2,\mathrm{},\gamma _M)`$ is in $`Lg_M(C[0,1]𝒜_{t_1})`$. Now consider the unital maximal full algebra of operator fields $`\stackrel{~}{𝒜}`$, which has base space $`[0,1]`$ and fibre algebras $`\{C[0,1]𝒜_t\}_{t[0,1]}`$. The continuity structure $`\stackrel{~}{}`$ for $`\stackrel{~}{𝒜}`$ consists of all operator fields of the form $`t_{i=1}^Nf_ic_i(t)`$, $`t[0,1]`$. Here the $`f_i`$s are in $`C[0,1]`$, the $`c_i`$s are continuous operator fields in $`𝒜`$ (with respect to the continuity structure $``$), and $`N`$ is a nonnegative integer. Now since $`\stackrel{~}{𝒜}`$ is a full algebra of operator fields, there are continuous operator fields $`t\gamma _j(.,t)`$ in $`\stackrel{~}{𝒜}`$ such that $`\gamma _j(.,t_1)=\gamma _j(.)`$, $`j=1,2,\mathrm{},M`$. (Here, the second variable of $`\gamma _j(.,.)`$ ranges over the base space of $`\stackrel{~}{𝒜}`$, and for $`t[0,1]`$, $`\gamma _j(.,t)`$ is an element of the fibre algebra $`C[0,1]𝒜_t`$ of $`\stackrel{~}{𝒜}`$).

One can show that the operator field $`t\gamma _j(1,t)`$ is a continuous operator field in $`𝒜`$ (with continuity structure $``$) for $`j=1,..,M`$. Also, one may view the field $`ta_j(t)`$ as an operator field in $`\stackrel{~}{𝒜}`$ in the natural way for $`j=1,..,M`$. From these and the fact that (as elements of $`C[0,1]𝒜_{t_1}`$) $`\gamma _j(.,t_1)`$ is within $`ϵ`$ of $`a_j(t_1)`$, we can find an open neighbourhood $`V`$ of $`t_1`$ with $`VI_1I_2`$ such that for all $`tV`$, a) as elements of $`C[0,1]𝒜_t`$, $`\gamma _j(.,t)`$ is within $`ϵ`$ of $`a_j(t)`$, $`j=1,2,\mathrm{},M`$, and b) $`\gamma _j(1,t)`$ is within $`ϵ/2`$ of $`f_{2,j}(t)`$. Furthermore, by continuity and since elements in a unital $`C^{}`$-algebra which are sufficiently close to the unit are invertible, we may assume that $`V`$ is sufficiently small so that for all $`tV`$, $`_{j=1}^M\gamma _j(.,t)^{}\gamma _j(.,t)`$ is an invertible element of $`C[0,1]𝒜_t`$.

One can show that for any continuous function $`g:[0,1][0,1]`$, the operator field $`t\gamma _j(g(t),t)`$ is a continuous field in $`𝒜`$, for $`j=1,2,\mathrm{},M`$. Hence, we can pick $`t_2V`$ such that $`t_1<t_2`$; and we can let $`\alpha _j^1(t)=_{df}\gamma _j((tt_1)/(t_2t_1))`$ for all $`t[t_1,t_2]`$, $`j=1,2,..,M`$. This will give us the first members of the sequences in (\*).

Now we can repeat almost exactly the same argument as before, replacing $`t_1`$ by the point $`t_2`$ and replacing $`f_{1,j}(t_1)`$ by $`\alpha _j^1(t_2)`$, $`j=1,2,\mathrm{},M`$. In this manner, we get a point $`t_3`$ and operator fields $`\alpha _j^2`$, $`j=1,2,\mathrm{},M`$, which will be the next members of the sequences in (\*). We need only note two minor modifications that are needed in the argument: a) firstly, one has to use the fact that for all $`t`$, the set

$$Lg_M(𝒜_t)\{(b_1,b_2,\mathrm{},b_M)𝒜_{t}^{}{}_{}{}^{M}:a_j(t)b_j<ϵ\text{ and }f_{2,j}(t)b_j<ϵ/2,\text{ for all }j\}$$

is a connected open set (see  Lemma 2.4), and b) when choosing the corresponding neighbourhood about $`t_2`$, one must make it sufficiently small so that the corresponding quantities which result will also be sufficiently small in order to fulfill condition (4) in (\*).

Repeating this process ad infinitum (making the appropriate modifications at each step), we get a sequence of points $`\{t_n\}_{n=1}^{\mathrm{}}`$ and sequences of continuous operator fields $`\{\alpha _j^n\}_{n=1}^{\mathrm{}}`$, $`j=1,2,\mathrm{},M`$, which fulfill the conditions in (\*). Now let $`\stackrel{~}{t}=lim_n\mathrm{}t_n`$. For $`j=1,2,\mathrm{},M`$, let $`\alpha _j`$ be the continuous operator field in $`𝒜`$ defined by

1. $`\alpha _j(t)=f_{1,j}(t)`$ for $`t[0,t_1]`$,
2. $`\alpha _j(t)=\alpha _j^n(t)`$ for $`t[t_n,t_{n+1}]`$, and
3. $`\alpha _j(t)=f_{2,j}(t)`$ for all $`t[\stackrel{~}{t},1]`$.

Then $`(\alpha _1,\alpha _2,\mathrm{},\alpha _M)Lg_M(𝒜)`$, and for all $`j=1,2,\mathrm{},M`$, $`\alpha _j`$ approximates $`a_j`$ within $`ϵ`$. ∎

###### Proof of Theorem 1.1.

We proceed by induction. The base case $`k=1`$ has already been dealt with in Lemma 2.1. We now do the induction step, supposing that $`k2`$. By  Theorem 4, let $`\pi :Prim(𝒜)[0,1]^k`$ be the continuous open surjection corresponding to the continuous field decomposition of $`𝒜`$ in the hypothesis (Here $`Prim(𝒜)`$ is the primitive ideal space of $`𝒜`$).

Since $`k2`$, the map $`p:[0,1]^k[0,1]^{k1}`$, given by projecting onto the first $`k1`$ coordinates, is a continuous open surjection. Hence, the composition $`p\pi :Prim(𝒜)[0,1]^{k1}`$ is a continuous open surjection. Hence by  Theorem 4, we can realize $`𝒜`$ as a unital maximal full algebra of operator fields with base space $`[0,1]^{k1}`$ and fibre algebras, say, $`\{_s\}_{s[0,1]^{k1}}`$. Hence by the induction hypothesis, $`sr(𝒜)sup_{s[0,1]}sr(C([0,1]^{k1})_s)`$.

But for each $`s`$, the fibre algebra $`_s`$ can be realized as a maximal full algebra of operator fields with base space $`p^1(s)=\{s\}\times [0,1]`$ and fibre algebras $`\{𝒜_r\}_{r\{s\}\times [0,1]}`$ (By  Theorem 4, $`_s`$ is isomorphic to $`𝒜/I`$ where $`I=(p\pi )^1(s)`$. Using this fact, one can construct the natural continuous open surjection of $`Prim(_s)`$ onto $`\{s\}\times [0,1]`$). Let us suppose that this continuous field decomposition of $`_s`$ is given by a continuity structure $`𝒢`$. Then $`C([0,1]^{k1})_s`$ can be realized as a unital maximal full algebra of operator fields with base space $`[0,1]`$ and fibre algebras $`\{C([0,1]^{k1})𝒜_r\}_{r\{s\}\times [0,1]}`$. Here the continuity structure consists of operator fields of the form $`r_{i=1}^Nf_ib_i(r)`$, for $`r\{s\}\times [0,1]`$. Here the $`f_i`$s are in $`C([0,1]^{k1})`$, the $`b_i`$s are continuous fields in $`_s`$ (with respect to the continuity structure $`𝒢`$) and $`N`$ is a nonnegative integer.

Hence by the induction hypothesis, for each $`s`$, we have that $`sr(C([0,1]^{k1})_s)sup_{r\{s\}\times [0,1]}sr(C([0,1]^k)𝒜_r)`$. It follows, then, that $`sr(𝒜)sup_{t[0,1]^k}sr(C([0,1]^k)𝒜_t)`$. ∎

We note that the statements of Lemma 2.1 and Theorem 1.1 would still be true if the unit interval $`[0,1]`$ was replaced by the circle $`S^1`$ (and if the $`k`$-cube $`[0,1]^k`$ was replaced by the $`k`$-torus $`𝕋^k`$). The proofs would be exactly the same. In other words, we have that

###### Corollary 2.2.

Let $`𝒜`$ be a unital maximal full algebra of operator fields with base space the $`k`$-torus $`𝕋^k`$ and fibre algebras $`\{𝒜_t\}_{t𝕋^k}`$. Then the stable rank of $`𝒜`$ satisfies

$$sr(𝒜)sup_{t𝕋^k}sr(C([0,1]^k)𝒜_t).$$

###### Theorem 2.3.

Let $`H_{2n+1}^{}`$ be the discrete Heisenberg group of rank $`2n+1`$. Let $`C^{}(H_{2n+1}^{})`$ be the univeral $`C^{}`$-algebra of $`H_{2n+1}^{}`$. Then $`sr(C^{}(H_{2n+1}^{}))=n+1`$.

###### Proof.

By  Proposition 1.7 and Theorem 4.3, and the fact that $`C(𝕋^{2n})`$ is a quotient of $`C^{}(H_{2n+1}^{})`$, the stable rank of $`C^{}(H_{2n+1}^{})`$ is greater than or equal to $`n+1`$. By  and  Theorem 3.4, $`C^{}(H_{2n+1}^{})`$ can be realized as a unital maximal full algebra of operator fields with base space the $`1`$-torus $`𝕋`$ and fibre algebras, say, $`\{𝒜_t\}_{t𝕋}`$. Hence, by the Corollary 2.2, the stable rank of $`C^{}(H_{2n+1}^{})`$ satisfies $`sr(C^{}(H_{2n+1}^{}))sup_{t𝕋}sr(C[0,1]𝒜_t)`$.

Now by  and  Theorem 3.4, each fibre algebra $`𝒜_t`$ can be realized as a unital maximal full algebra of operator fields with base space a torus with dimension less than or equal to $`2n`$ (the zero-dimensional torus being a point) and fibre algebras, say, $`\{B_s^t\}_{s𝕋^{l_t}}`$ where $`𝕋^{l_t}`$ is the base in the continuous field decomposition of $`𝒜_t`$ ($`l_t2n`$). Moreover, for each $`t𝕋`$ and for each $`s𝕋^{l_t}`$, $`_s^t`$ is isomorphic to either a full matrix algebra $`𝕄_n()`$ or $`𝕄_n()(^n𝔸_\theta )`$, where $`^n𝔸_\theta `$ is the $`n`$ times tensor product of a fixed irrational rotation algebra $`𝔸_\theta `$ with irrational rotation $`\theta `$. Now by  each irrational rotation algebra can be decomposed as an inductive limit of building blocks of the form $`M_m(C(𝕋))𝕄_n(C(𝕋))`$ (the integers $`m`$ and $`n`$ get arbitrarily large as we move along building blocks in the inductive limit). Hence, since $`C[0,1]𝒜_t`$ can be realized as a maximal full algebra of operator fields with base space $`𝕋^{l_t}`$ and fibre algebras $`\{C[0,1]_s^t\}_{s𝕋}`$, it follows, by Corollary 2.2,  Proposition 1.7 and Theorems 5.1 and 6.1, that $`sr(C[0,1]𝒜_t)n+1`$, for all $`t`$. Hence, $`sr(C^{}(H_{2n+1}^{})=n+1`$. ∎

## References

* J. Anderson and W. Paschke, ‘The rotation algebra’, *Houston J. Math.*, 15, (1989), 1 - 26.
* L. Baggett and J. Packer, ‘The primitive ideal space of two-step nilpotent group $`C^{}`$-algebras’, *J. Funct. Anal.*, 124, (1994), 389-426.
* K. R. Davidson, *$`C^{}`$-algebras by example*, Fields Institute Monographs, 6.
* G. A. Elliott and D. E. Evans, ‘The structure of the irrational rotation $`C^{}`$-algebra’, *Ann. of Math. (2)*, 138, no. 3, (1993), 477-501.
* J. M. G. Fell, ‘The structure of algebras of operator fields’, *Acta Math.*, 106, (1961), 233-280.
* N. E. Hassan, ‘Rangs Stables de certaines extensions’, *J. London Math. Soc.*, (2), 52, (1995), 605-624.
* R. Herman, L. N. Vaserstein, ‘The stable rank of $`C^{}`$-algebras’, *Invent. Math.*, 77, (1984), 553-555.
* R-Y Lee, ‘On $`C^{}`$-algebras of operator fields’, *Indiana U. Math. Journal*, 25, no. 4, (1976), 303-314.
* S. T. Lee and J. Packer, ‘Twisted group $`C^{}`$-algebras for two-step nilpotent and generalized discrete Heisenberg group’, *J. Operator Th.*, 34, no. 1, (1995), 91-124.
* V. Nistor, ‘Stable range for the tensor products of extensions of $`K`$ by $`C(X)`$’, *J. Operator Th.*, 16, (1986), 387-396.
* J. Packer and I. Raeburn, ‘On the structure of twisted group $`C^{}`$-algebras’, *Trans. AMS*, 334, no. 2, (1992), 685-718.
* M. A. Rieffel, ‘Dimension and stable rank in the $`K`$-theory of $`C^{}`$-algebras’, *Proc. London Math. Soc.*, (3), 46, (1983), 301-333.
* T. Sudo, ‘Dimension theory of group $`C^{}`$-algebras of connected Lie groups of type I’, *J. Math. Soc. Japan*, 52, no. 3, (2000), 583-590.
* T. Sudo, ’Structure of group $`C^{}`$-algebras of Lie semi-direct products $`\times |`$’, *J. Operator Theory*, 46, (2001), 25-38.
* T. Sudo and H. Takai, ‘Stable rank of the $`C^{}`$-algebras of nilpotent Lie groups’, *Internat. J. Math.*, 6, no. 3, (1995), 439-446.
* T. Sudo and H. Takai, ‘Stable rank of the $`C^{}`$-algebras of solvable Lie groups of type I’, *J. Operator Theory*, 38, no. 1, (1997), 67-86.
* J. Tomiyama, ‘Topological representation of $`C^{}`$-algebras’, *Tohoku Math. J.*, (2), 14, (1962), 187-204.