R Romero

R López-Ibáñez

E A Dalchiele

J R Ramos-Barrado

F Martín

D Leinen

Compositional and physico-optical characterization of 05% Al-doped zinc oxide films prepared by chemical spray pyrolysis

Keywords: 81.15.Rs, 82.80.Pv, 68.53.Dv, 78.66.Fd thin films, chemical synthesis, doping, optical properties

come   L'archive ouverte pluridisciplinaire

Introduction

Al-doped ZnO is a multifunctional material. In the form of thin layers, it has been studied for different applications as transparent conductive oxide, infrared mirror, antireflection layer, gas sensor, etc [1]. Al-doped ZnO films are commonly prepared by sputtering or CVD [2][3][4][5]. However, the preparation by chemical spray pyrolysis (CSP) is an interesting alternative due to the potentially very low installation and maintenance cost of the technique. Chemical spraying can be done in open air without the need of special chambers and pumps for a controlled atmosphere or vacuum. Furthermore, CSP can be considered as environmentally friendly since mainly aqueous precursor solutions and air as driving gas are used. Sometimes, depending on the substrate and the material to be deposited, it might be necessary to add small quantities of ethanol (typically 50 ml per liter) and/or acetic acid (1 -10 ml per liter) to the aqueous precursor solution in order to stabilize the precursor (∼10 -2 M) in the aqueous solution by forming stable acylates [6] and to enhance the wettability of the solution to the substrate [7]. Further advantages of spray pyrolysis are: compact or porous films can be obtained depending on the spray conditions; doping can easily be performed by introduction of the doping element into the precursor solution; CSP in open atmosphere can easily be up-scaled to large area deposition [8]. A review on the possibilities of CSP was given by Patil [9].

For CSP, the selection criteria for a precursor in aqueous spray solution is determined by its chemical stability, at least over the period of deposition time, its ease of availability and its cost. This leads directly, in case of ZnO films, to the selection of zinc acetate and zinc chloride as precursors. Zinc alcoxides could also be used, but have to be stabilized in an aqueous spray solution and are more expensive. Previously, we reported about pure ZnO films deposited on silicon by spray pyrolysis using zinc acetate precursor. We showed that a maximum growth rate of 15 nm/min can be reached for substrate temperatures around 543 K [10] and we studied the electrical properties of those films in detail in the temperature range 223 -373 K [11]. We observed different preferred orientations in wurtzite type crystal structure using chloride or acetate zinc precursors and effects on crystal texture with increasing Al doping, as published recently [12]. In the present paper, we study the effect on film composition and physicooptical properties of the films when introducing increasing amounts of Al (0 -5 %) into the precursor solutions.

Experimental

In CSP, the precursor solution is sprayed onto the heated substrate. In our lab scale spray station, the substrate is fixed to a metal plate heated using an electrically resistive wiring, at the desired temperature of pyrolysis. During spraying, the substrate is moved underneath the spray nozzle, fixed in position, back and forth at a constant frequency, to obtain a deposit of homogeneous thickness. A stream of compressed air gas through the nozzle, free of oil, is used for the atomization of the solution to very fine droplets. Zinc acetate (Zn(CH 3 COO)

= = I V t I V t 53 . 4 2 ln π ρ (2)
where t is the film thickness of the sample, in our case determined from the optical transmission spectra (see below). The resistivity was highest of about 10 Ωm for the undoped ZnO films and lowest of about 10 -2 Ωm for 3% Al doped films.

Results and discussion

Surface analytical characterization

The chemical composition at the surface and in-depth was determined by XPS combined with 4 keV Ar + sputtering for ZnO:Al thin films obtained from both types of Zn precursors and for Al doping between 0 % and 5 %. As an example, Figure 1 shows the XPS spectra of an ZnO:Al film using zinc acetate precursor and 4 % aluminium doping, of its surface and after 1 min Ar + sputtering. A unique Zn 2p 3/2 peak and the principal contribution in the O 1s peak is observed at 1021.8 eV and 530.3 eV in binding energy, respectively, both typical for ZnO [START_REF] Moulder | Handbook of X-ray Photoelectron Spectroscopy[END_REF]. Furthermore, we observe in the O 1s peak another contribution around 531.7 eV which together with the C 1s peak at 288.7 eV clearly indicates the presence of high amounts of carboxylic species at the film surface [START_REF] Moulder | Handbook of X-ray Photoelectron Spectroscopy[END_REF]. In case of using Zn chloride precursor, spectra are similar but also chlorine is detected. The Cl 2p 3/2 peak is observed at 198.8 eV, a value near to the one documented for ZnCl 2 [START_REF] Moulder | Handbook of X-ray Photoelectron Spectroscopy[END_REF]. After sputtering the film surface, the Cl 2p signal persists showing that chlorine is homogeneously distributed with depth at a concentration of approximately 1 at% (spectra not shown here).

The amount of Al in the films was determined from the Al 2p peak, centred at 73.8 eV.

Table 1 shows as an example the XPS atomic concentrations determined for both type of samples, i.e. obtained from Zn acetate and Zn chloride precursor. As seen in Figure 1, the intensities of Zn and O increase after mild sputtering due to the elimination of the surface contamination in form of aliphatic carbon and carboxylic species. According to Table 1, Zn and O reach nearly the stoichiometric ratio of ZnO. However, the Al 2p peak is found in more or less the same intensity indicating that Al is present in excess at the film surface relative to the amount found inside the film. So, there seems to be some segregation of Al to the film surface during film growth. Figure 1 includes also the Al KLL Auger peak measured at a kinetic energy of 1388.4 eV. The Al 2p and Al KLL energies clearly show that aluminium is incorporated in the film in oxidised state [START_REF] Moulder | Handbook of X-ray Photoelectron Spectroscopy[END_REF].

The modified Auger parameter α' Al determined from the energy position of both peaks is a measure for charge polarisability [START_REF] Wagner | Faraday Discuss[END_REF] in the surroundings of Al atoms in the ZnO film. The value of α' Al is 1462.2 eV, slightly higher than in Al 2 O 3 and suggests that aluminium is in a substoichiometric environment relative to Al 2 O 3 [START_REF] Moulder | Handbook of X-ray Photoelectron Spectroscopy[END_REF]. This can be expected for the ZnO lattice which offers less O atoms for each cation compared to the Al 2 O 3 lattice.

As a quantitative summary, Figure 2 This shows that only mild sputtering is applicable for the analysis of the films.

As hydrogen cannot be detected with XPS, IR hemispherical reflectance measurements were performed for both types of samples with and without Al doping. Figure 3 represents the corresponding reflectance spectra between 2 and 4 µm in comparison to the fused silica substrate itself. For these measurements the samples were positioned on the sample port of an integrating gold sphere, covering its back with the gold cover plate. Thus we observe nearly total reflectance for that MIR region where fused silica is still transparent unless showing a typical sharp absorption at about 2.74 µm due to occlusion of water molecules remaining from silica fusion [START_REF] Socrates | Infrared and Raman characteristic group frequencies[END_REF]. Water incorporated in the bulk of films, prepared by spray pyrolysis from aqueous precursor solutions, would show a broad absorption due to the O-H stretch, from about 2.8 to 3.5 µm [START_REF] López | [END_REF]. The absence of such an absorption feature in Figure 3 shows that no water inclusions are present in the bulk of the films, giving evidence of an effective pyrolysis, curing and densification of the films at a substrate temperature of 623 K.

Optical characterization

Figure 4 shows the optical transmission spectra of A) undoped ZnO films and B) 3% Al doped ZnO films for increasing film thickness, i.e. increasing spray time, in all cases obtained from the Zn-acetate precursor (films obtained from Zn-chloride precursor have similar transmission spectra). As seen in Figure 4, ZnO and ZnO:Al films are highly transparent. The overall transmittance τ of the films was calculated between 400 nm and 900 nm as τ = T(λ)dλ / T glass (λ)dλ. For ZnO films, the transmittance τ is about 0.92, nearly independent of the film thickness, whereas for ZnO:Al films, the transmittance τ decreases with film thickness from 0.93 to 0.86. These values of transmittance are in the class of the best obtained for ZnO:Al films produced with spray pyrolysis [19][20][21][22] and comparable to other highly transparent ZnO films produced with more expensive physical or chemical techniques [23,24]. The decrease in transmittance at higher film thickness in ZnO:Al films may be attributed to the increased scattering of photons by crystal defects created by the Al doping [12,20] and a higher surface roughness in ZnO:Al films compared to ZnO films, as observed by Li et al. with atomic force microscopy [25]. A qualitative difference between Figure 4 A) and B) is that oscillations in the transmission spectra observed for undoped films, due to an interference effect at the interfaces (air / transparent film) and (transparent film / transparent substrate), are not observed in ZnO:Al films. It shows that Al doping leads to scattering of light and thus to the suppression of any optical interference. This may be attributed to a higher surface roughness in ZnO:Al films [25] and, in particular, the change in crystal growth from highly oriented in ZnO films to random with smaller average crystallite size in ZnO:Al films as observed previously with XRD [12].

In case of the undoped ZnO films, the film thickness has been calculated from transmission maxima and minima according to Swanepoel's method [26] given for transparent films on transparent substrates. In order to improve the precision of the method, in particular, when the films are very thin and a transmission maximum or minimum is difficult to determine, a curve has been fitted to each transmission spectrum with thickness d as basic fitting parameter, yielding thus thickness values of more confidence. Table 2 shows the ZnO film thickness values obtained. The index of refraction has also been calculated according to Swanepoel's method from the envelopes of transmission maxima and minima taking into account the whole set of transmission spectra of ZnO films on glass represented in Figure 4. For the transparent region, from 900 nm to 600 nm, the index of refraction is almost constant, increasing weakly from 1.9 to 2.0 whereas in the weak to medium absorbing region, 600 nm to 400 nm, the index of refraction increases more strongly up to about 2.15. These values are in agreement to data published on ZnO and Al-doped ZnO films studied with ellipsometry [27].

The effect of increasing film thickness in Al doped and undoped films can also be 2. We have obtained for this wavelength a mean value of (154000 ± 6000) cm -1 , which is in excellent agreement with absorption coefficient measurements derived from transmission spectra at room temperature of ZnO single crystal epitaxial films deposited by pulsed laser deposition and film thickness determination with surface profile measurements [29]. It agrees also well with the absorption coefficient observed for energies above the absorption edge in Ga doped ZnO standing free films [28]. Once, the mean absorption coefficient of the ZnO films was determined, we have used Lambert-Beer's law to estimate the film thickness of the 3 % Al doped ZnO:Al films. Those values are shown in the last two columns of Table 2.

The indicated error was calculated according to error propagation and depends mainly on the mean deviation <α>. However, for very low transmission values (thicker films) systematic errors may become important and exceed the indicated error values. In this context, we would like to point out that this method gives the possibility of film thickness estimation for thin films from optical transmission spectra where no interference phenomena are observed. Furthermore, a mean deposition rate of about 3.5 nm per minute can be deduced from Table 2 for both types of films. However, for very large spray times (>90 min) a reduction in deposition rate can be noticed. In order to confirm the correctness of ZnO:Al film thickness estimation, the cross section produced by fracture of the 3 % Al doped ZnO:Al film of 90 min spraying was examined with SEM. As seen in Figure 5, the measured thickness values indicated in the enlarged section of the cross section of the film agree well with the one given in Table 2 and thus show consistency for the aforementioned method of thickness estimation for ZnO:Al films within the given error ranges.

Figure 6 shows the absorption coefficient α versus photon energy hν of both series of films, calculated from each transmission spectrum shown in Figure 4 and the corresponding film thickness value given in Table 2. As can be seen, the absorption coefficient curves group together for the ZnO films (Figure 6 A)) and for the 3 % Al doped ZnO:Al films (Figure 6 B)), not showing any dependence on the film thickness.

This shows that the films are very homogeneous and reproducible. Only for the thinnest ZnO film (15 min spraying), a higher absorption coefficient was obtained, which might be due to the fact that the interface between film and glass substrate becomes more important for very thin films during their initial stages of film growth by the spray pyrolysis technique. Furthermore, the optical band gap E g was determined for each film, Finally, comparing films of similar thickness but different Al concentrations, Figure 8 shows the transmission spectra of ZnO films, obtained for 90 min spraying (Zn-acetate precursor) and increasing amounts of Al doping (0 % to 5 % AlCl 3 ). Figure 8 shows again how interference effects are suppressed in the transmission spectra when doping with Al. In addition, we observe a decrease of steepness in the ZnO band edge and its shift to higher values for increasing Al concentration. This can be seen more clearly in Figure 9 which represents the corresponding absorption coefficient curves calculated from the transmission spectra of Figure 8. The film thickness values have been determined as explained above, using the transmission of each film above the absorption edge and the absorption coefficient α = 154000 cm -1 . Figure 9 shows clearly how the absorption edge shifts to higher energies and broadens with increasing Al concentration. Both, broadening and blue shift can be explained with the Burstein-Moss effect [30,[START_REF] Moss | Optical Properties of Semiconductors[END_REF]: extra electrons introduced by Al doping into the semiconductor fill up the bottom of the conduction band making thus available for excitation only sites of higher energy E(k) which leads to an increase of the optical band gap and to a broadening of the absorption edge through a larger transition distribution of E(k) at room temperature.

Figure 10 shows the corresponding band gap values obtained for increasing Al concentration by extrapolation to zero absorption, fitting a straight line to the absorption edge in (αhν) 2 versus hν representation. In case of 0 % and 3 % Al doping the E g mean values aforementioned from Figure 5 are represented. In Figure 10, one can observe a nearly linear increase of the optical band gap with increasing Al concentration in the films. This blue shift of the optical band gap with increasing Al concentration has been reported previously for magnetron sputtered ZnO:Al films [START_REF] Sernelius | [END_REF] but not studied so far in films obtained with the spray pyrolysis technique and considering such a large range of Al doping (0 -5%). However, even for a fixed Al concentration, a large variation of band gap values can be found in the literatures, depending on the type of measurement and on the method of band gap evaluation. In some published work [33], the point of inflection at the band edge of the transmission spectra is directly used for optical band gap estimation. However, this can only be done for the comparison of films of the same film thickness, since the band edge in the transmission spectra is highly sensitive to changes in film thickness.

Table 3 shows a collection of bang gap values E g taken from the literature for ZnO and

Al doped ZnO films. The method of determining E g is also indicated for each value given in Table 3. As one can see, the spread of E g values is tremendous but also general tendencies depending on the method can be noticed. The lowest E g values have been obtained using transmission measurements, and if available together with reflection measurements, in order to calculate directly the absorption coefficient knowing the film thickness and thus determine the band gap value E g , as in this work, according to direct type transitions, (αhν) 2 = A(E -E g ), for (αhν) 2 = 0. Somewhat higher values can be found in work where ellipsometry data were used and complex dielectric function, refractive index and extinction coefficient were calculated and from the latter the absorption coefficient according to α = 4πk/λ, and then E g as above. A comparison of E g values using both methods was given by Li et al. [27] and, more or less, 0.1 eV higher values were obtained using ellipsometry compared to transmission data in ZnO and Al-doped ZnO films (see Table 3). Less difference in E g values were obtained for a ZnO film of 357 nm thickness studied by Senadim et al. [34] who compared also two methods: on the one hand, the absorption coefficient calculated from the transmission spectrum, as explained above, and on the other hand, measuring the reflectance spectrum and calculating, by the use of the Kramers-Kronig transformation, the dielectric function, refractive index and extinction coefficient and from the latter the absorption coefficient.

Comparatively even higher band gap values were obtained looking at the excitonic features in the absorption spectrum which are visible in single crystal ZnO films even at room temperature due to the high exciton binding energy of about 60 meV in ZnO [29],

compared to the thermal energy k B T. Three papers using excitonic features are mentioned in the upper part of Table 3, two of them using transmission spectra data and one of them ellipsometry data. Contrary to what has been mentioned above using the absorption edge for E g determination, in this case, the Eg value obtained from ellipsometry data is lower than the one based on transmission data.

Furthermore, in Al doped ZnO films, a shift to very high E g values with increasing Al concentration was reported by Sernelius et al. [START_REF] Sernelius | [END_REF], not comparable, at least according to our knowledge, to band gap shifts observed in other work for similar Al concentration.

Besides all aforementioned, the band gap in ZnO and Al doped ZnO films depends on the method of film preparation, the quality of crystalline film which can be obtained by the method, the substrate on which the film is deposited and the temperature of film growth, substrate, and/or post annealing. Various examples can be found in the recollection of data given in Table 3. Higher substrate temperatures in magnetron sputtered films lead to an increase in crystalline quality and somewhat higher E g values [25]. Non crystalline substrates as fused silica lead to lower E g values compared to crystalline substrates, as C or R plane sapphire [35]. Furthermore, the steepness of the absorption edge was found to be considerably less when using fused silica as substrate and was assigned to smaller grain sizes and more defects [35], which is consistent to our observations of absorption edge broadening with increasing Al concentration as seen in Figure 9.

Finally we would like to mention that we have obtained for our ZnO films on glass band gap values which are in the upper range of the corresponding ones tabulated in Table 3.

This indicates a good crystalline quality of our films. In addition, constant E g values for films of different thickness as seen in Figure 7 for ZnO and 3 % Al doped ZnO:Al films show good reproducibility of our films prepared by spray pyrolysis. However, the optical characterization has also shown that it is important to have films of a certain film thickness when evaluating the effects of Al doping on the optical properties of the films.

Conclusions

  represents the Al concentration in the films determined with XPS versus the Al concentration in the precursor solution. concentrations of various series of film preparation, with and without mild surface etching by 4 keV Ar +sputtering, are presented. The precursor for ZnO was Zn acetate or Zn chloride. In all cases, the precursor for Al doping was AlCl 3 . The principal information of Figure2is that we obtain for each series a linear relationship between the Al concentration in the film and the Al concentration in the precursor solution. This shows that Al doping in spray pyrolysis deposited zinc oxide films can easily be controlled. The as-deposited films reveal an excess of aluminium at the film surface whereas shortly sputtered films show a nearly one to one ratio. Further sputtering as used for depth profiling leads to preferential loss of Al atoms and thus to a constant Al concentration of about 1 % in the ion bombardment altered layer, independently of the initial Al concentration in the film.

5 %

 5 observed in the transmission spectra for energies beyond the ZnO:Al absorption edge, not exceeding the absorption edge of the glass substrate. According to Lambert-Beer's law of extinction the transmission shrinks for that energies with increasing film thickness coming to almost zero transmission for thicker films. That is equally observed for Al doped and undoped ZnO films although the ZnO:Al absorption edge shifts slightly with Al concentration as discussed below. In that energy region, the transmission around 340 nm, equivalent to 3.65 eV, (position of the minima of the first derivative of the transmission spectra of ZnO films; position that shifts slightly, as the absorption edge does, in ZnO:Al films with increasing Al concentration; figure not shown) is directly related to the amount of ZnO:Al deposited onto the glass substrate and can thus be used to estimate the film thickness, at least for low film thicknesses, of the Al doped films. To do so, we suppose that the absorption coefficient above the absorption edge of each film does not change considerably between Al doped and undoped ZnO films, similarly to what has been observed in Ga doped ZnO films up to Ga doping [28]. We have calculated the absorption coefficient of ZnO at 340 nm from the transmission values of the undoped films of different thickness according to α = [ln(T glass /T)]/d and indicated in Table

Figure 7 ,

 7 Figure 7, showing an almost constant value for each series of films independent of the

  Highly transparent Al doped zinc oxide films without water inclusions can be obtained by chemical spray pyrolysis spraying aqueous zinc solutions. The amount of Al doping in the films shows to be proportional to the one in the spray solution and thus can be controlled easily. 0 to 5 % Al doping shifts the ZnO band gap by about 50 meV to higher values and leads to absorption edge broadening. The absorption coefficient of films of different thickness and Al-doping was calculated from the transmission spectra. The thickness of the ZnO:Al films was estimated from the absorption of each film using the absorption coefficient of ZnO films. The thickness estimated for the ZnO:Al film of 3 % Al doping and 90 min spray time agrees well with its thickness observed by SEM. Comparison of the absorption coefficient of films of different spray time, i.e. film thickness, showed good reproducibility in the optical properties of the films. The authors are grateful to the Spanish Ministry of Education and Science (MEC) for funding (TEC2007-60996) and Consolider Ingenio (CSD2008-00023) and to the Junta de Andalucia through the research group FQM-192.

28 TFigure 1 :

 281 Figure Captions:

Figure 2 :

 2 Figure 2: Al concentration in the film determined with XPS versus Al concentration in the precursor solution, both as 100%•Al/(Al+Zn), of various series of film preparation with and without 4 keV Ar + surface etching. Straight lines are linear fits to the data points of each series.

Figure 3 :

 3 Figure 3: IR hemispherical reflectance spectra of pure and 3 % Al doped ZnO films from acetate or chloride precursor, deposited on fused silica. The reflectance spectrum of the substrate is also represented for comparison.

Figure 4 :

 4 Figure 4: Optical transmission spectra of A) undoped ZnO films on glass substrate and B) 3% Al doped ZnO films on glass substrate, for increasing spray time deposition (in minutes).

Figure 5 :

 5 Figure 5: SEM image of the cross section produced by fracture of the 3 % Al doped ZnO:Al film of 30 min spraying onto glass substrate (its transmission spectrum is shown in Figure 4 B)). For the enlarged area of the film cross section, three values (nm) of observed film thickness are indicated.

Figure 6 :

 6 Figure 6: Absorption coefficient curves calculated from the transmission spectra shown in Figure 4 and film thickness values given in Table 2. A) undoped ZnO films on glass and B) 3% Al doped ZnO films on glass.

Figure 7 :

 7 Figure 7: Optical band gap values E g versus film thickness of the undoped ZnO films and the 3% Al doped ZnO films.

Figure 8 :

 8 Figure 8: Optical transmission spectra of Al doped ZnO films on glass substrate for increasing Al concentration. Estimated film thickness: 334 nm (0 % Al), 268 nm (1 % Al), 252 nm (2 % Al), 285 nm (3 % Al), 253 nm (4 % Al) and 234 nm (5 % Al).

Figure 9 :Figure 8

 98 Figure 9: Absorption coefficient curves calculated from the transmission spectra of Figure 8 of 0 -5 % Al doped ZnO films.

Figure 10 :

 10 Figure 10: Optical band gap values E g versus Al concentration (x) in the Al doped ZnO films.

  2

  •2H 2 O) and zinc chloride (ZnCl 2 ) in concentrations of 10 -2 M, fixing the number of metallic atoms in pure water, were the precursors for ZnO. Al doping was performed by adding to the precursor solution aluminium chloride (AlCl 3 ) The pressure in the analysis chamber was about 10 -7 Pa. The sputter rate at 4 keV Ar + and beam rastering of 3×3 mm 2 is assumed to be in oxide material approximately 3 nm/min as determined in Ta 2 O 5 under identical sputter

	conditions. Binding energies were referenced to the C 1s peak at 284.8 eV [13]. Atomic
	concentrations were determined from C 1s, Cl 2p, O 1s, Zn 2p 3/2, Si 2p and Al 2p XPS
	peak areas using Shirley background subtraction [14] and sensitivity factors provided by
	the spectrometer manufacturer (Physical Electronics, Eden Prairie, MN).
	Optical transmission measurements were made in the UV-VIS range (200 nm -900 nm)
	up to 5 % calculated as	Zn	Al +	Al	100 •	%	. Substrates were glass and fused silica of size 2
	cm by 2 cm. The injection rate of the precursor solution using a syringe pump was 50
	cm 3 /h. The air pressure was 1.5 bar, resulting in an air flux through the nozzle of
	approximately 20 l/min. The nozzle to substrate distance was 20 cm, the substrate
	temperature 623 K (350ºC). Samples of various spraying times, i.e. different film
	thickness, were prepared.

The chemical composition was studied by XPS combined with 4 keV Ar + sputter depth profiling with a PHI 5700 equipment. A standard X-ray source, 15 kV, 300 W, Mg K α (1253.6 eV) was used.

with a Lambda 19 spectrometer from Perkin-Elmer. Spectral hemispherical reflectance measurements were made in the MIR range (2 µm -17µm) with an IFS66 from Bruker with external integrating gold sphere of 20 cm in diameter. Scanning electron microscopy (SEM) pictures were obtained with a JEOL JSM 5300 apparatus. The electrical resistivity ρ of the films was measured at room temperature on various sites of each sample with a home-made four point probe conductivity tester and calculated after the formula for thin films

[15]

Table 1 :

 1 XPS atomic concentrations of ZnO:Al films using Zn acetate or Zn chloride precursors. Zn/O ratio is given after surface etching.

	Sample	4 keV Ar +	Zn	O	Al	C	Cl Zn/O
		sputter time /
		min
	ZnO:Al Film from Fig. 1	0	11.1 32.7 2.2 54.0	-	-
	Zn acetate; 4 % Al	1	41.9 50.7 2.5	4.9	-	0.83
	ZnO:Al Film	0	11.9 34.7 2.3 50.5 0.6	-
	Zn chloride; 3.5 % Al	1	46.0 50.4 2.9	-	0.7 0.85

Table 2 :

 2 A) Absorption coefficient α calculated at 340 nm of ZnO films on glass for increasing spray time (film thickness d determined after Swanepoel[26]). The mean value α for ZnO films is given below. B) Film thicknesses d and estimated error ∆d of

	ZnO:Al (3%Al) films for increasing spray times, calculated after Lambert-Beer's law of
	extinction as given below using the mean value α . T are the measured transmission
	values of ZnO and ZnO:Al films at 340 nm and 335 nm, respectively; T glass is 0.87 and
	0.86, respectively.
		A) ZnO		B) ZnO:Al3%
	t /min	T	d / nm	α / nm -1	T	d / nm	∆d / nm
	15	0.3960	40	0.01900	0.4276	45	5
	30	0.2051	90	0.01606	0.2487	80	10
	45	0.1042	132	0.01607	0.0797	150	10
	60	0.0525	176	0.01595	0.0365	200	10
	75	0.0196	247	0.01536	0.0203	240	10
	90	0.0061	334	0.01485	0.0105	280	20
	105	0.0042	361	0.01476	0.0099	290	20
	120	0.0037	368	0.01481	0.0084	300	20
		α = (0.0154 ± 0.0006) nm -1	d = [ln(T glass /T)] ⁄ α

Table 3 :

 3 Band gap values from the literature for ZnO and Al-doped ZnO films.Indicated is also the type of substrate, the manner of film preparation, the measurement / method used for band gap determination and the Al concentration.

	film / preparation / author (year) ZnO single crystal film / PVD / Liang	method T α (excitonic feature) 0 Al conc.	E g (eV) 3.445 (4.2 K)
	and Hoffe (1968) [36]				3.40
	ZnO single crystal on sapphire / PLD	T α (excitonic feature) 0	3.455 (77 K)
	/ Muth et al. (1999) [29]				3.40
	ZnO single crystal on sapphire / PLD	ε α (excitonic feature) 0	3.372 (295 K)
	/ Washington (1998) [37]
	Al:ZnO on Si(100) / RF(ZnO) + DC(Al) magnetron sputtering (800ºC post annealing) / Li et al. (2008) [27]	T α, ε α		0 4 (wt%)	3.26 (T) 3.34 (ε) 3.32 (T)
					3.45 (ε)
	Al:ZnO on CaF 2 / RF(ZnO) +	T + R dielectric	0	3.40 (dα/dE)
	DC(Al) magnetron sputtering /	function (Fresnel
	Sernelius et al. (1988) [32]	equations)	α	2.14 (at%)	3.90 (dα/dE)
	ZnO on sapphire and fused silica /	T α		0	3.29 (0001)
	PLD at 600ºC / Srikant et al. (1997)				3.32 (1102)
	[35]				3.26 (f.s.)
					3.28 (f.s.) p.a. 800ºC
	ZnO on glass / PFCVAD / Senadim
	et al. (2006) [34]

% Al

Compositional and optical characterization of Al doped ZnO films prepared by CSP

Acknowledgements


Oliver Rinne
email: o.rinne@damtp.cam.ac.uk

An axisymmetric evolution code for the Einstein equations on hyperboloidal slices

Keywords: numbers: 04.20.Cv, 04.20.Ha, 04.25.D-, 04.30.-w

We present the first stable dynamical numerical evolutions of the Einstein equations in terms of a conformally rescaled metric on hyperboloidal hypersurfaces extending to future null infinity. Axisymmetry is imposed in order to reduce the computational cost. The formulation is based on an earlier axisymmetric evolution scheme, adapted to time slices of constant mean curvature. Ideas from a previous study by Moncrief and the author are applied in order to regularize the formally singular evolution equations at future null infinity. Longterm stable and convergent evolutions of Schwarzschild spacetime are obtained, including a gravitational perturbation. The Bondi news function is evaluated at future null infinity.

Introduction

A major problem in numerical relativity that is currently attracting considerable interest is the global treatment of entire asymptotically flat spacetimes. This is relevant both to the modelling of the gravitational radiation emitted by compact astrophysical sources, and to several questions in mathematical relativity such as the nonlinear stability of black holes and the quantitative asymptotic behaviour of perturbations of such spacetimes.

In the standard approach based on the Cauchy formulation, one truncates the spatial domain at a finite distance from the source, where boundary conditions must be supplied (see [START_REF] Sarbach | Absorbing boundary conditions for Einstein's field equations[END_REF] for a review). Apart from being compatible with the Einstein constraint equations and yielding a well-posed initial-boundary value problem, these boundary conditions must carefully take into account the propagation of gravitational radiation. For instance, crude choices of boundary conditions are known to lead to incorrect results on radiation tails in the evolution of perturbed black hole spacetimes [START_REF] Allen | Radiation tails and boundary conditions for black hole evolutions[END_REF][START_REF] Dafermos | A note on boundary value problems for black hole evolutions[END_REF]. Recently, improved absorbing boundary conditions that minimize spurious reflections of gravitational radiation have been derived and implemented [START_REF] Buchman | Towards absorbing outer boundaries in general relativity[END_REF][START_REF] Buchman | Improved outer boundary conditions for Einstein's field equations[END_REF][START_REF] Rinne | Implementation of higher-order absorbing boundary conditions for the Einstein equations[END_REF]. While constituting a significant improvement, this approach still relies on the assumption that the Einstein equations may be linearized about a given background spacetime so that incoming and outgoing radiation can be identified-an assumption that may not always be valid [START_REF] Deadman | Numerical relativity and asymptotic flatness[END_REF]. In the full nonlinear theory there is no satisfactory Confidential: not for distribution. Submitted to IOP Publishing for peer review 8 December 2009 quasi-local definition of a gravitational energy flux at a finite distance from the source [START_REF] Szabados | Quasi-local energy-momentum and angular momentum in general relativity Living[END_REF]. This fact also complicates the extraction of gravitational radiation in an invariant manner.

Gravitational radiation is only unambiguously defined at future null infinity I + [START_REF] Bondi | Gravitational waves in general relativity VII. Waves from axi-symmetric isolated systems[END_REF][START_REF] Sachs | Gravitational waves in general relativity VIII. Waves in asymptotically flat space-time[END_REF]. Thus it is very desirable to include I + in the computational domain. There are several ways of doing this. One approach is characteristic evolution, whereby spacetime is foliated by null hypersurfaces that can be compactified towards I + (see [START_REF] Winicour | Characteristic evolution and matching Living[END_REF] for a review). This works well in the far field but characteristic foliations are ill behaved in strong-field regions due to the formation of caustics in the null congruences generating the hypersurfaces. A compromise is Cauchy-characteristic matching, where a Cauchy evolution in the interior is matched to a characteristic evolution in an exterior region. In [START_REF] Reisswig | Unambiguous determination of gravitational waveforms from binary black hole mergers[END_REF] this has recently been applied to the a posteriori extraction of gravitational radiation from a Cauchy evolution (with finite outer boundary) of a binary black hole coalescence.

A more efficient solution is a foliation of spacetime into hyperboloidal hypersurfaces, which are everywhere spacelike but approach future null infinity rather than spacelike infinity. Such surfaces are not Cauchy surfaces-the foliation covers only the region to the future of some initial hyperboloidal hypersurface but not, for instance, the region around spacelike infinity. This is not a serious disadvantage for the applications we have in mind as we are mainly interested in computing the gravitational radiation at I + and studying the late-time behaviour of perturbations. Hyperboloidal foliations have the advantage of being as flexible as standard Cauchy foliations in the interior, so that a variety of gauge conditions can be used that have been known to be successful in strong-field situations. In [START_REF] Friedrich | Cauchy problems for the conformal vacuum field equations in general relativity[END_REF] Friedrich reformulated the Einstein equations for a conformally rescaled metric in such a way that they are manifestly regular at I + . This symmetric hyperbolic system is larger than the actual Einstein equations, e.g. it contains evolution equations for the Weyl curvature. For reviews of the theoretical development as well as numerical implementations, we refer the reader to [START_REF] Frauendiener | Conformal infinity Living[END_REF][START_REF] Husa | Problems and successes in the numerical approach to the conformal field equations[END_REF][START_REF] Husa | Numerical relativity with the conformal field equations[END_REF]. Perhaps most relevant to the present work is an axisymmetric reduction of the regular conformal field equations by Frauendiener and Hein [START_REF] Frauendiener | Numerical evolution of axisymmetric, isolated systems in general relativity[END_REF] that was applied to the numerical evolution of Minkowski spacetime and the boost-axisymmetric solutions of Bičák and Schmidt [START_REF] Bičák | Asymptotically flat radiative space-times with boost-rotation symmetry: The general structure[END_REF].

Despite these successes, most recent work in numerical relativity has used formulations of the Einstein equations (usually on truncated Cauchy hypersurfaces) rather than the regular conformal field equations. For this reason, it is worth investigating whether some formulation of the Einstein equations themselves can be adapted to hyperboloidal foliations reaching out to I + . As for the regular conformal field equations, one applies a conformal transformation to the metric in order to map I + to a finite coordinate location in a compactified coordinate system, where the conformal factor vanishes. When written in terms of the conformal metric, the Einstein equations contain inverse powers of the conformal factor, terms that are singular at I + . Numerical implementations must face the question of how to deal with these formally singular terms, which may be potential sources of instabilities.

One proposal for solving the Einstein equations on hyperboloidal slices was developed by Zenginoglu [START_REF] Zenginoglu | Hyperbolodial evolution with the Einstein equations[END_REF]. It combines a conformal transformation of the spacetime metric (where the conformal factor can essentially be freely specified) with a certain choice of gauge source functions for the Einstein equations in (generalized) harmonic gauge [START_REF] Friedrich | On the hyperbolicity of Einstein's and other gauge field equations[END_REF]. Preliminary numerical results on spherically symmetric evolutions of Schwarzschild spacetime were reported in section 2.5 of [START_REF] Zenginoglu | A conformal approach to numerical calculations of asymptotically flat spacetimes[END_REF]. However, it is not yet clear how the formally singular terms at I + are to be evaluated numerically in more general situations.

Here we follow a different approach due to Moncrief that is based on an ADMlike [START_REF] Arnowitt | The dynamics of general relativity Gravitation: an introduction to current research[END_REF] formulation of the Einstein equations with a constant-mean-curvature (CMC) slicing condition. In [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] we showed explicitly how the formally singular terms in the evolution equations at I + can be evaluated there in terms of conformally regular geometric data. The constraint equations that occur in this approach were solved numerically in [START_REF] Buchman | Black hole initial data on hyperboloidal slices[END_REF] and initial data on CMC slices describing various configurations of single and binary black holes were constructed. The present paper is concerned with the evolution problem. We assume spacetime to be axisymmetric in order to speed up the code so that different ways of treating I + numerically can be experimented with. We apply the ideas of [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] to an earlier axisymmetric evolution scheme on Cauchy slices [START_REF] Rinne | Constrained evolution in axisymmetry and the gravitational collapse of prolate Brill waves[END_REF] (see also [START_REF] Garfinkle | Numerical evolution of Brill waves[END_REF][START_REF] Choptuik | An axisymmetric gravitational collapse code[END_REF] for related schemes). The objective of the paper is to show that long-term stable evolutions are possible, including gravitational radiation.

The paper is organized as follows. In section 2 we explain our notation and gauge choices and obtain the constraint and evolution equations. Particular emphasis is placed on how the formally singular evolution equations can be regularized at I + . In section 3 we describe our numerical method for solving these equations. The numerical results are presented in section 4. As a first test problem we evolve Schwarzschild spacetime and demonstrate the convergence of the numerical solution. Next, a gravitational perturbation is included. A regular expression at I + for the Bondi news function [START_REF] Bondi | Gravitational waves in general relativity VII. Waves from axi-symmetric isolated systems[END_REF] describing the outgoing gravitational radiation is derived and evaluated numerically. We conclude and discuss possible directions for future work in section 5.

Formulation

Definitions and gauge choices

We assume that spacetime is axisymmetric with Killing vector ξ. For simplicity we also assume here that ξ is hypersurface orthogonal. Spherical polar coordinates t, r, θ, φ are chosen such that ξ = ∂/∂φ.

The spacetime metric (4) g αβ is written as (4)

g αβ = ψ -2(4) γ αβ , ( 1
)
where ψ is the conformal factor and (4) γ αβ the conformal metric. Greek indices α, β, . . . are spacetime indices (t, r, θ, φ). The conformal factor ψ is positive and approaches zero at future null infinity, which we put at a fixed coordinate location r = 1. The conformal metric (4) γ αβ is assumed to be regular there. We decompose the metric in ADM form as (4)

g αβ dx α dx β = -α 2 dt 2 + g ij (dx i + β i dt)(dx j + β j dt), (2)
where g ij is the spatial metric, α the lapse function and β i the shift vector. Latin indices i, j, . . . are spatial (r, θ, φ). Similarly, we write the conformal spacetime metric as (4)

γ αβ dx α dx β = -α 2 dt 2 + γ ij (dx i + β i dt)(dx j + β j dt), (3)
where γ ij = ψ 2 g ij is the conformal spatial metric and α = ψα the conformal lapse.

The spatial coordinates are chosen such that β φ = 0 and the conformal spatial metric takes the form γ ij dx i dx j = e 2η sin θ (dr 2 + r 2 dθ 2 ) + r 2 sin 2 θ dφ 2 .

(

) 4
This is sometimes called quasi-isotropic gauge. We note that this choice of gauge differs from the spatial harmonic gauge used in [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF]. ‡ The extrinsic curvature is defined by

∂ t g ij = -2αK ij + L β g ij , ( 5
)
where L denotes the Lie derivative. Preservation of the conditions g θθ = r 2 g rr and g rθ = 0 under the evolution equation ( 5) implies

β r ,r -β θ ,θ -r -1 β r = α(K r r -K θ θ ) ≡ α Ũ, (6)
β r ,θ + r 2 β θ ,r = 2αK r θ ≡ 2α Kr θ . ( 7
)
Similarly as in [START_REF] Rinne | Constrained evolution in axisymmetry and the gravitational collapse of prolate Brill waves[END_REF] we have defined the following components of the extrinsic curvature,

K r r -K θ θ ≡ ψ Ũ, K r θ = r -2 K θ r ≡ ψ Kr θ , (8)
and, motivated by regularity concerns on the axis of symmetry [START_REF] Rinne | Constrained evolution in axisymmetry and the gravitational collapse of prolate Brill waves[END_REF][START_REF] Garfinkle | Numerical evolution of Brill waves[END_REF][START_REF] Choptuik | An axisymmetric gravitational collapse code[END_REF], we also set

K θ θ -K φ φ ≡ ψ sin θ W . ( 9
)
It is easy to see that the quantities Ũ , Kr θ and W defined in ( 8) and ( 9) have the same conformal weight as the traceless momentum

π tr ij = √ g(K ij -1 3 g ij K l l ) (10)
used in [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF], where g denotes the determinant of g ij . Note also that because Kr θ , W and Ũ are either off-diagonal or differences of diagonal components of the extrinsic curvature, they can in fact be regarded as the three independent components of the traceless part of the extrinsic curvature in our case.

As in [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] we require the time slices to have constant mean curvature,

g ij K ij ≡ -K = const. ( 11
)
Note the slightly awkward sign in the definition of the constant K; with our sign convention for the extrinsic curvature [START_REF] Buchman | Improved outer boundary conditions for Einstein's field equations[END_REF] we require K to be positive so that the slices reach future null infinity. Preservation of [START_REF] Winicour | Characteristic evolution and matching Living[END_REF] under the Einstein evolution equation for the extrinsic curvature yields the following elliptic equation for the conformal lapse,

0 = α,rr + 2r -1 α,r + r -2 (α ,θθ + c s α,θ ) -1 2 sα(η ,rr + r -1 η ,r + r -2 η ,θθ + 2r -2 c s η ,θ -r -2 η) + 3 2 α P r (P r -2 Ãr ) + r -2 P θ (P θ -2 Ãθ ) -1 6 ψ -2 αe 2sη K 2 -5 2 αe 2sη 1 3 ( Ũ + 1 2 s W ) 2 + 1 4 (s W ) 2 + r -2 ( Kr θ ) 2 . ( 12
)
Here and in the following we use the abbreviations

s ≡ sin θ, c ≡ cos θ, P B ≡ ψ -1 ψ ,B , ÃB ≡ α-1 α,B , (13)
upper-case Latin indices A, B, . . . ranging over r and θ. ‡ It is, in a sense, its two-dimensional analogue: if we define a two-dimensional metric H AB by the part of the line element (4) orthogonal to the Killing vector, H AB dx A dx B = e 2η sin θ (dr 2 + r 2 dθ 2 ), then

V C ≡ H AB (Γ C AB -Γ ˚C AB ) = 0,

Constraint and evolution equations

As usual in ADM-like reductions, the Einstein equations split into constraint and evolution equations. The Hamiltonian constraint is

0 = ψ ,rr + 2r -1 ψ ,r + r -2 (ψ ,θθ + c s ψ ,θ ) -3 2 ψ -1 (ψ ,r ) 2 + r -2 (ψ ,θ ) 2 -1 2 sψ(η ,rr + r -1 η ,r + r -2 η ,θθ + 2r -2 c s η ,θ -r -2 η) -1 2 ψe 2sη 1 3 ( Ũ + 1 2 s W ) 2 + 1 4 (s W ) 2 + r -2 ( Kr θ ) 2 + 1 6 ψ -1 e 2sη K 2 . ( 14
)
The momentum constraints are

0 = C r ≡ 2 3 Ũ,r + 1 3 s W,r + r -2 Kr θ,θ + r -2 Kr θ ( c s + 2cη + 2sη ,θ -2P θ ) + Ũ(sη ,r -4 3 P r + 2r -1 ) + s W (r -1 -2 3 P r ), (15)
0 = C θ ≡ -1 3 Ũ,θ + 1 3 s W,θ + Kr θ,r + 2 Kr θ (sη ,r -P r + r -1 ) + Ũ(-cη -sη ,θ + 2 3 P θ ) + W ( 4 3 c -2 3 sP θ ). ( 16
)
We have the following evolution equations,

ψ ,t = β r ψ ,r + β θ ψ ,θ -( c s β θ + r -1 β r ) -1 3 α(K + ψ( Ũ + 2s W ), (17)
η ,t = β r η ,r + β θ η ,θ + c s β θ η + (s -1 β θ ) ,θ -α W , (18) W,t = β r W,r + β θ W,θ + (2 c s β θ + r -1 β r ) W +e -2sη r -2 -(s -1 α,θ ) ,θ + 2αψ -1 (s -1 ψ ,θ ) ,θ -αe -2sη η ,rr + 2r -1 η ,r + r -2 [η ,θθ -η + c(s -1 η) ,θ ] +( Ãr -2P r )η ,r -r -2 s -1 ( Ãθ -2P θ )(sη ,θ + cη) -2s -1 Kr θ β θ ,r + 2 3 αs -1 s W (s W + 1 2 Ũ -ψ -1 K) + 3r -2 ( Kr θ ) 2 (19) Kr θ ,t = β r Kr θ ,r + β θ Kr θ,θ + c s β θ Kr θ + e -2sη (-α ,rθ + 2αψ -1 ψ ,rθ ) +αe -2sη ( Ãr -2P r + r -1 )(sη ,θ + cη) + ( Ãθ -2P θ )(sη ,r + r -1 ) +cη ,r + 2 3 α Kr θ (s W -ψ -1 K + 2 Ũ) -r 2 β θ ,r Ũ, (20)
Ũ,t = β r Ũ,r + β θ Ũ,θ + Ũ( c s β θ + r -1 β r ) +e -2sη [-α ,rr + r -2 α,θθ + 2αψ -1 (ψ ,rr -r -2 ψ ,θθ )] +αe -2sη [( Ãr -2P r )(2sη ,r + r -1 ) + 2sr -1 η ,r -2r -2 ( Ãθ -2P θ )(sη ,θ + cη + c s )] + 1 3 α Ũ (2s W -2ψ -1 K + Ũ) + 4 Kr θ (β θ ,r -αr -2 Kr θ ). ( 21
)
We observe that all the equations are manifestly regular at the axis of symmetry θ = 0 when the parities of the fields with respect to θ = 0 are taken into account (table 1). These can be inferred from the general behaviour of regular axisymmetric tensor fields [START_REF] Rinne | A strongly hyperbolic and regular reduction of Einstein's equations for axisymmetric spacetimes[END_REF].

Partially constrained evolution and hyperbolicity

We determine the conformal factor ψ by solving the momentum constraint [START_REF] Frauendiener | Conformal infinity Living[END_REF] rather than the evolution equation [START_REF] Frauendiener | Numerical evolution of axisymmetric, isolated systems in general relativity[END_REF]. On the other hand, the extrinsic curvature variables W , Kr θ and Ũ are evolved using the evolution equations ( 19)-( 21) and the Table 1. Angular parities of the fundamental variables: even (+) or odd (-).

The parities about θ = 0 follow from regularity on axis, those about θ = π/2 from the reflection symmetry we impose.

Variable α 15) and ( 16) are not solved explicitly. Hence we adopt a partially constrained evolution scheme.

β r β θ ψ η W Kr θ Ũ θ = 0 + + -+ -- - + θ = π/2 + + -+ + + - + momentum constraints (
Of course the constraints are preserved by the evolution equations, which imply the following subsidiary system,

C t ,t β r C t ,r + β θ C t ,θ -1 2 αψ(C r ,r + r -2 C θ ,θ ), (22)
C r ,t

β r C r ,r + β θ C r ,θ + 4 3 e -2sη αψ -1 C t ,r , (23)
C θ ,r β r C θ ,r + β θ C θ ,θ + 4 3 e -2sη αψ -1 C t ,θ , (24) (25)
where the symbol signifies that only the principal parts have been displayed, and the remaining terms are homogeneous in the constraints. If the Hamiltonian constraint is enforced exactly, C t = 0, then the remaining evolution equations for C r and C θ are clearly hyperbolic (they are advection equations along the shift). If the Hamiltonian constraint is not enforced then the full constraint evolution system is not hyperbolic (the characteristic speeds

λ r = β r ± 2 3 e -sη α i, λ θ = β θ ± 2 3 r -1 e -sη α i (26)
are complex). This might be cured by adding multiples of the Hamiltonian constraint to the evolution equations; however, we have not found such a constraint addition that both yields a hyperbolic constraint evolution system and is compatible with regularity of the main evolution equations on the axis of symmetry. Hence a free evolution scheme that solves none of the constraint equations during the evolution does not appear to be feasible. However, our formulation includes three elliptic equations for the gauge variables α and β A anyway and solving one additional elliptic equation for the conformal factor is not overly expensive.

If we assume that α, β A and ψ are given by the solution of the elliptic equations, then the pair of evolution equations ( 18) and ( 19) forms a wave equation to principal parts,

(∂ t -β A ∂ A ) 2 η α2 e -2sη (η ,rr + r -2 η ,θθ ), ( 27
)
and the evolution equations ( 20) and ( 21) are advection equations along the shift to principal parts. Hence the evolution equations are hyperbolic. A rigorous proof of the well posedness of the present mixed elliptic-hyperbolic hyperboloidal initial value problem remains to be carried out. For a similar mixed elliptic-hyperbolic formulation on a Cauchy foliation, such a proof was given in in [START_REF] Andersson | Elliptic-hyperbolic systems and the Einstein equations[END_REF].

Regularity at future null infinity

The CMC slicing condition [START_REF] Reisswig | Unambiguous determination of gravitational waveforms from binary black hole mergers[END_REF], the constraint equations ( 14)-( 16) and the evolution equations for the extrinsic curvature ( 19)-( 21) are formally singular at I + , where the conformal factor ψ vanishes. Our elliptic solver (section 3.4) is well adapted to such degenerate elliptic equations. However, the right-hand sides of the evolution equations cannot be evaluated at I + in their present form. In [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] two methods for obtaining regular forms of the evolution equations at I + were presented: the first based on Taylor expansions, the second on an argument due to Penrose [START_REF] Penrose | Zero rest-mass fields including gravitation: Asymptotic behaviour[END_REF]. For the numerical results presented in section 4 we followed the second approach but we have also implemented the first approach, as discussed further below.

Penrose showed that provided certain smoothness assumptions hold (in particular, the conformal metric must be C 3 up to the boundary), the conformal Weyl tensor (i.e., corresponding to the conformal metric (4) γ αβ ) vanishes at I + . This implies [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] that the electric part of the physical Weyl tensor (i.e., corresponding to (4) g αβ ) also vanishes at I + ,

E αγ = n β n δ C αβγδ [ (4) g] = 0, ( 28
)
where n α is the unit normal to the t = const slices and = denotes equality at I + . It is straightforward to compute the Weyl tensor for the metric (4) γ αβ defined in ( 3) and ( 4). Any time derivatives are substituted using the evolution equations. (Here a computer algebra program is very useful; we used REDUCE [START_REF] Hearn | REDUCE User's Manual 3[END_REF].) Setting E αγ = 0 immediately gives regular expressions for the formally singular terms in ( 19)-( 21).

Thus we obtain the following manifestly regular evolution equations at

I + , W,t = β r W,r + β θ W,θ + (2 c s β θ + r -1 β r ) W -e -2sη r -2 (s -1 α,θ ) ,θ +αe -2sη η ,rr + 2r -1 η ,r + r -2 [η ,θθ -η + c(s -1 η) ,θ ] -A r η ,r + r -2 s -1 A θ (sη ,θ + cη) -2s -1 Kr θ β θ ,r + α[4s -1 r -2 ( Kr θ ) 2 -Ũ W ], ( 29
) Kr θ ,t = β r Kr θ ,r + β θ Kr θ,θ + c s β θ Kr θ -e -2sη α,rθ +αe -2sη ( Ãr -r -1 )(sη ,θ + cη) + Ãθ (sη ,r + r -1 ) -cη ,r +2α Kr θ (s W + Ũ ) -r 2 β θ ,r Ũ, ( 30
) Ũ,t = β r Ũ,r + β θ Ũ,θ + Ũ( c s β θ + r -1 β r ) + e -2sη (-α ,rr + r -2 α,θθ ) +αe -2sη [ Ãr (2sη ,r + r -1 ) -2r -2 ( Ãθ -c s )(sη ,θ + cη) -2sr -1 η ,r ] +α Ũ (2s W + Ũ ) + 4 Kr θ (β θ ,r -αr -2 Kr θ ). (31)
The other approach to regularity at I + discussed in [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] is based on an expansion of all the fields in (finite) Taylor series with respect to radius r about I + . Substituting these expansions in the singular constraint equations, one obtains conditions on the fields and their radial derivatives that can be used in order to evaluate the formally singular terms in the evolution equations. This approach is more general in that it does not require Penrose's smoothness assumption.

As necessary conditions for a regular evolution at I + , we obtained in [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] a set of conditions that had been discovered in [START_REF] Andersson | On the regularity of solutions to the Yamabe equation and the existence of smooth hyperboloidal initial data for Einstein's field equations[END_REF], namely that the shear of I + and the radial components π tr ri of the traceless momentum vanish. In the axisymmetric formulation we consider here, these conditions read W =e -sη η ,r , Ũ = 1 2 e -sη sη ,r ,

Kr θ = 0. ( 32
)
One might be tempted to impose these conditions as Dirichlet conditions on the extrinsic curvature at I + , thus avoiding to evaluate the singular evolution equations there. However, numerical evolutions with these Dirichlet boundary conditions were found to be unstable. This is perhaps not surprising because at I + all the fields are purely outgoing and hence imposing boundary conditions there is ill posed (although of course the conditions [START_REF] Andersson | On the regularity of solutions to the Yamabe equation and the existence of smooth hyperboloidal initial data for Einstein's field equations[END_REF] are satisfied for a given solution of the field equations). It is worthwhile to observe though, as we already did in [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF], that the evolution equations ( 19)-( 21) contain remarkable damping terms that cause potential violations of the regularity conditions [START_REF] Andersson | On the regularity of solutions to the Yamabe equation and the existence of smooth hyperboloidal initial data for Einstein's field equations[END_REF] to decay exponentially in the neighbourhood of I + :

u ,t = . . . -2 3 αψ -1 Ku ( 33
)
for each u ∈ { W , Kr θ , Ũ }. The term has the "right sign" because K > 0. Instead of imposing the regularity conditions [START_REF] Andersson | On the regularity of solutions to the Yamabe equation and the existence of smooth hyperboloidal initial data for Einstein's field equations[END_REF] directly, however, we can extract more information from the Taylor expansions that allows us to evaluate the singular terms in the evolution equations directly [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF]. The resulting modified evolution equations at I + are detailed for the current axisymmetric scheme in Appendix B. Their main difference to ( 29)-( 31) lies in the presence of terms involving first radial derivatives of the extrinsic curvature at I + (cf. equation (B.9)). We found this version of the evolution equations at I + to be numerically stable as well, and the results agree with those in section 4 within numerical error (see Appendix B for details).

It is worth noting here that the quasi-isotropic spatial gauge we use differs from the spatially harmonic gauge of [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF]. Nevertheless, we were able to apply the same procedure for evaluating the formally singular terms at I + . This provides some substance to the claim in [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] that the precise form of the spatial gauge conditions is inessential for the method to be applicable.

Numerical method

Having obtained a set of gauge conditions, constraints and evolution equations, we now describe how we solve them numerically. First we summarize the evolution scheme, i.e. which variables are evolved using which equations. Next we specify our choice for the computational domain and the discretization of the partial differential equations. We then focus on the two main aspects of the code, the time integration scheme and the elliptic solver. The code has been written in C++. Unless otherwise indicated all routines have been developed from scratch.

Evolution scheme

Our fundamental variables are α, β A , ψ, η, W , Kr θ and Ũ . The variables η, W , Kr θ and Ũ are evolved using the evolution equations ( 18)- [START_REF] Zenginoglu | A conformal approach to numerical calculations of asymptotically flat spacetimes[END_REF]. At each time step we solve the CMC slicing condition [START_REF] Reisswig | Unambiguous determination of gravitational waveforms from binary black hole mergers[END_REF] for α, the quasi-isotropic gauge conditions ( 6) and ( 7) for β A , and the Hamiltonian constraint ( 14) for ψ. Details on how the quasi-isotropic gauge conditions are treated are given in Appendix A.1. The momentum constraints are only solved initially (Appendix A.2).

Domain and discretization

For the applications considered in this paper, the spatial numerical domain is taken to be the quarter of an annulus, r min r 1 and 0 θ π/2. The inner boundary at r = r min will be taken to lie within a black hole event horizon (black hole excision), and the outer boundary at r = 1 corresponds to I + . The axis of symmetry is at θ = 0, and we impose an additional reflection symmetry about the plane θ = π/2.

The fields typically have the steepest gradients close to the black hole and hence it is advisable to use a non-uniform radial grid in order to provide more numerical resolution towards the inner boundary. We introduce a quadratic map

r(x) = Qx 2 + (1 -r min -Q)x + r min , ( 34
)
where 0 Q < 1 is a constant. This map satisfies r(0) = r min and r(1) = 1. We use a uniform grid in x, i.e. we place grid points at x i = i/N r , 0 i N r , where N r is the number of radial grid points. Note in particular that there are grid points right at the inner boundary and at I + . In the angular direction we use a uniform staggered grid:

θ j = π 2 (j -1 2 )/N θ , 1 i N θ .
Two layers of ghost points are added on either side, corresponding to j = -1, 0 and j = N θ+1 , N θ+2 . Values at these ghost points are set according to the angular parities of the various fields, which are listed in table 1. E.g. for a field u that is odd about θ = 0 we set u i0 = -u i,1 and u i,-1 = -u i2 , and for an even field we set u i0 = u i,1 and u i,-1 = u i2 . Here u ij denotes the value of u at the grid point with indices (i, j).

The spatial derivatives are discretized using fourth-order accurate finite differences. One-sided differences are used in the radial direction at the boundaries; everywhere else centred differences are used (see Appendix C for details).

Time integration

A fourth-order Runge-Kutta method is used in order to integrate the evolution equations forward in time. At each substep of this method, the elliptic equations are solved as described below in section 3.4. Since all the characteristics point towards the exterior of the domain both at the inner excision boundary and at I + , no boundary conditions are imposed there-the evolution equations are discretized using one-sided finite differences. At I + the evolution equations ( 19)- [START_REF] Zenginoglu | A conformal approach to numerical calculations of asymptotically flat spacetimes[END_REF] for the extrinsic curvature are replaced with their regularized versions (29)- [START_REF] Hearn | REDUCE User's Manual 3[END_REF]. Kreiss-Oliger dissipation [START_REF] Kreiss | Methods for the approximate solution of time dependent problems Global Atmospheric Research Programme[END_REF] is added in order to ensure stability (Appendix C). No modifications to the dissipation operators or special smoothing operations at I + have been required.

Elliptic solver and boundary conditions

The elliptic equations are solved by means of a multigrid method [START_REF] Brandt | Multilevel adaptive solutions to boundary value problems[END_REF][START_REF] Briggs | A Multigrid Tutorial[END_REF] using Full Approximation Storage (FAS). Because of the singular nature of the equations at I + (r = 1), we have found it necessary to use line relaxation in the radial direction, i.e. all points on a line θ = const are solved simultaneously. With our fourth-order finite-difference discretization this requires solving a penta-diagonal linear system of equations, for which a fast algorithm exist (we use the routines bandec and banbks in [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF]). For the Hamiltonian constraint ( 14) an outer Newton-Raphson iteration is employed for each such line solve; all remaining elliptic equations are linear. The lines θ = const are then traversed and updated in ascending order j = 1, 2, . . . , N θ (Gauss-Seidel relaxation in the θ-direction). We have found this method to be very robust and effective, leading to a reduction of the residual of about an order of magnitude per W-cycle. No problems with convergence at the singular boundary were encountered.

The elliptic equations require boundary conditions at the inner r = r min and outer r = 1 boundary (in addition to the parity boundary conditions at θ = 0, π/2), which we shall now specify.

At the outer boundary (I + ) we take the angular shift to vanish, β θ = 0. The conformal factor vanishes there by definition, ψ = 0. Preservation of this condition under the evolution equation ( 17) yields a boundary condition for the conformal lapse, α =e sη β r .

(

) 35
We remark that these boundary conditions for the lapse and shift ensure that the normal n α to the time slices is null at I + . The radial shift β r appearing in ( 35) cannot be chosen freely but is determined by the quasi-isotropic gauge conditions. After each Gauss-Seidel relaxation sweep we integrate [START_REF] Deadman | Numerical relativity and asymptotic flatness[END_REF] along the boundary to determine β r there. The integration constant can be fixed e.g. by requiring that the average of β r along the boundary agree with the value obtained from the exact Schwarzschild solution that we consider in section 4.

At the inner excision boundary, we also take β θ to vanish, and we freeze α to the value corresponding to the exact Schwarzschild solution. However, the conformal factor ψ cannot be chosen freely. This is because ψ obeys an evolution equation with outgoing characteristic speed at the excision boundary, equation [START_REF] Frauendiener | Numerical evolution of axisymmetric, isolated systems in general relativity[END_REF]. In order to determine the correct value of ψ there, we keep a copy of ψ that is evolved using [START_REF] Frauendiener | Numerical evolution of axisymmetric, isolated systems in general relativity[END_REF]. At the end of each timestep this copy of ψ is reset to agree with the solution of the Hamiltonian constraint. This procedure was found to be crucial for ensuring that the momentum constraints are preserved during the numerical evolution.

Numerical results

In this section, we apply our code to the evolution of a (perturbed) Schwarzschild black hole. First we focus on the unperturbed Schwarzschild solution and demonstrate that long-term stable and convergent evolutions can be obtained. Next we include a small gravitational perturbation. The outgoing gravitational radiation as represented by the Bondi news function is computed at I + .

Schwarzschild spacetime

First we describe how the exact solution is computed, which is needed both for the boundary conditions of some of the elliptic equations (section 3.4) and in order to compare with the numerical solution during the evolution. The Schwarzschild solution in CMC slicing was derived in [START_REF] Brill D R | K-surfaces in the Schwarzschild space-time and the construction of lattice cosmologies[END_REF][START_REF] Malec | Constant mean curvature slices in the extended Schwarzschild solution and the collapse of the lapse[END_REF]. Its line element is

ds 2 = -1 - 2M r dt 2 + 1 f 2 dr 2 - 2a f dt dr + r2 (dθ 2 +sin 2 θ dφ 2 ), ( 36
)
where

f (r) = 1 - 2M r + a 2 1/2 , a(r) = K r 3 - C r2 , ( 37
)
and M (mass), K (mean curvature) and C are constants.

As in [START_REF] Buchman | Black hole initial data on hyperboloidal slices[END_REF] we introduce a new radial coordinate r (the coordinate used in our evolution scheme) by requiring the spatial metric to be conformally flat in the new coordinates. Hence r obeys dr dr

= r rf (r) . ( 38
)
This ordinary differential equation is integrated numerically to high precision using Mathematica. The boundary condition is that r → 1 as r → ∞.

The metric variables defined in section 2.1 are then given by

ψ = r r , α = rf r , β r = - ra r , β θ = 0, η = 0. ( 39
)
The extrinsic curvature variables are found to be

Ũ = - 3C r2 r , W = 0, Kr θ = 0, ( 40
)
and the vector V A introduced in order to solve the momentum constraints initially

(Appendix A.2) is V r = C r 2 , V θ = 0. ( 41
)
For the numerical evolutions presented in this paper we choose M = 1, K = 1/2 and C = 2 (although we have obtained stable evolutions across a wide range of parameter values). The inner boundary is placed at r min = 0.05, just inside the horizon at r = 2 ⇔ r = 0.0635.

We run the simulation with two different numerical resolutions (N r , N θ ) = (64, 8) and (128, 16). The radial grid is nonuniform with Q = 3/4 in (34) in both cases.

The time step ∆t is limited by the Courant-Friedrichs-Lewy condition

∆t < min ∆r |v r ± | , ∆θ |v θ ± | , ( 42
)
where ∆r is the grid spacing and v r ± are the incoming and outgoing characteristic speeds of the evolution system ( 17)-( 21) in the r-direction (similarly in the θdirection), and the minimum is taken over all grid points. For our choice of parameters the time step is constrained by the outgoing radial characteristic speed at

I + , v r + = -β r + e -sη α = 2K 3 (43)
for the Schwarzschild solution. The time step is taken to be about 3/4 of the maximum allowed value computed in this way: for the lower resolution we use ∆t = 0.06.

During the evolution, we compute the difference of the numerical solution with respect to the exact solution for each of the evolved variables and take its discrete L 2 norm. This is divided by the norm of the exact solution for the corresponding variable whenever the latter does not vanish. Of these results a vector 2-norm (square root of the sum of the squares) is taken. The resulting quantity represents a measure of the overall relative error of the numerical solution as a function of time and is plotted in figure 1. The error grows only linearly with time. It decreases by a factor of about 15 on average as the resolution is doubled, close to the expected convergence factor of 2 4 = 16 for a fourth-order accurate scheme. The main contribution to the error comes from the region close to the black hole horizon where the fields have the steepest gradients. The numerical solution remains smooth all the way up to I + during the entire evolution.

Perturbed Schwarzschild spacetime

Next, we include a small gravitational perturbation in our evolutions. This is done by taking the variable η, which vanishes for the exact Schwarzschild solution, to be a Gaussian initially, The factor sin θ has been included for parity reasons (table 1). We choose A = 10 -4 , r c = 0.5 and σ = 0.05. We stress that these initial data are evolved with the full nonlinear Einstein equations. The linearized problem has been studied extensively; in particular, numerical evolutions of the Bardeen-Press equation [START_REF] Zenginoglu | Gravitational perturbations of Schwarzschild spacetime at null infinity and the hyperboloidal initial value problem[END_REF] and the Regge-Wheeler-Zerilli equations [START_REF] Zenginoglu | Asymptotics of black hole perturbations[END_REF] have been performed on hyperboloidal slices.

η = A sin θ exp - (r -r c ) 2 2σ 2 . ( 44
)
In order to demonstrate the consistency and convergence of our code, we monitor the momentum constraints, which are not solved explicitly during the evolution. For the unperturbed Schwarzschild solution, each term in the θ-momentum constraint [START_REF] Husa | Numerical relativity with the conformal field equations[END_REF] vanishes separately. This is not true for the r-momentum constraint [START_REF] Husa | Problems and successes in the numerical approach to the conformal field equations[END_REF]; it vanishes only because of a subtle cancellation of the terms 2 3 Ũ,r + Ũ(- 4 3 P r + 2r -1 ) = 0 [START_REF] Jr | Kinematics and dynamics of general relativity Sources of gravitational radiation[END_REF] for the Schwarzschild solution. Each separate term in this equation is large (in particular close to the inner boundary) and the violation of the r-momentum constraint in the perturbed evolutions is dominated by the finite-differencing error of these terms for the Schwarzschild background. Therefore we have divided the r-momentum constraint by the norm of a typical term in (45) (roughly 10 2 ) in order to give a better idea of the relative size of the constraint violation. The results are shown in figure 2. The violation of the momentum constraints decreases by just over an order of magnitude as the resolution is doubled, reasonably close to the expected convergence factor of 2 4 = 16.

Of prime interest is the gravitational radiation emitted by the perturbed black hole as measured at I + . It can be represented by the Bondi news function [START_REF] Bondi | Gravitational waves in general relativity VII. Waves from axi-symmetric isolated systems[END_REF]. This function was defined originally by studying the asymptotic behaviour of the metric in a special coordinate system. This coordinate system will in general not agree with the coordinates used in a numerical evolution. Fortunately, Stewart [START_REF] Stewart | Numerical relativity III. The Bondi mass revisited[END_REF] has developed a method for computing the news in an arbitrary coordinate system, working as we do here in terms of a conformally related metric. We begin by specifying a Newman-Penrose [START_REF] Newman | An approach to gravitational radiation by a method of spin coefficients[END_REF] null tetrad (l, n, m, m) that is adapted to I + in the sense that 15)-( 16) as functions of time for two different resolutions (Nr, N θ ) = (64, 8) (dashed lines) and (128, 16) (solid lines). The rmomentum constraint (left panel) has been normalized by the size of a typical term in [START_REF] Jr | Kinematics and dynamics of general relativity Sources of gravitational radiation[END_REF] as explained in the text.

m α ∂ α ψ = 0. ( 46

For definiteness we choose

m = 1 √ 2r e -sη ∂ ∂θ + i csc θ ∂ ∂φ . ( 47
)
Apart from the requirement (46) the tetrad is arbitrary. The resulting expressions will be invariant under spin and boost transformations. From the conformal spacetime metric (4) γ αβ specified in (3) and (4), we compute the conformal Ricci tensor (4) Rαβ . The time derivatives in this expression are substituted by using the evolution equations ( 17)- [START_REF] Zenginoglu | A conformal approach to numerical calculations of asymptotically flat spacetimes[END_REF]. The Bondi news function is now given by N = Φ02 ≡ mα mβ [START_REF] Buchman | Towards absorbing outer boundaries in general relativity[END_REF] Rαβ .

(48)

The resulting expression is not yet regular at I + but it can be written in a regular form by using the regularity conditions [START_REF] Andersson | On the regularity of solutions to the Yamabe equation and the existence of smooth hyperboloidal initial data for Einstein's field equations[END_REF], the boundary condition [START_REF] Briggs | A Multigrid Tutorial[END_REF] and expressions (B.2) and (B.3) for the radial derivatives of the conformal factor. Ultimately we obtain N = -e -2sη s(η ,rr + η ,θθ ) -s 2 (η ,r ) 2 + cη ,θ -s -1 η -e -sη α-1 β θ cη ,r -e -sη s W,r .(49)

The Bondi news function has spin weight -2 and hence we may expand it in spinweighted spherical harmonics [START_REF] Newman | Note on the Bondi-Metzner-Sachs group[END_REF] -2 Y m with m = 0 because of the axisymmetry,

N (t, θ)| I + = ∞ =2 N (t) -2 Y 0 (θ). ( 50
)
The = 2 contribution is plotted in figure 3. The higher-contributions are considerably smaller (by at least a factor 10 -4 ) and have not yet converged for the numerical resolutions used here.

After an initial peak when the outgoing part of the wave reaches I + , we expect the news to consist of a ringing phase due to the quasi-normal mode excitations of the black hole, followed by a polynomial decay ("tail") due to the backscatter of the radiation off the curved background spacetime. These features are apparent in figure 3. However, whereas the ringing phase has converged quite well, the tail is still highly resolution-dependent. Considerably higher resolutions and/or a higher order of accuracy of the discretization will be needed in order to study the tail behaviour. We expect the quasi-normal mode ringing to have the form

N ∝ e -κ t sin(ω t + ϕ ). (51)
The decay rate κ and frequency ω are estimated by performing a χ 2 fit of (51) to the numerical data for the news at I + in the time interval 60 t 120. We find κ 2 = 0.0893 ± 0.0009 and ω 2 = 0.3738 ± 0.0012. These values correspond to the high-resolution run and the errors have been estimated by comparing with the lowresolution run. This numerical result is consistent with the linearized-theory result from Leaver's continued-fraction method [START_REF] Leaver | An analytic representation for the quasi-normal modes of Kerr black holes[END_REF], κ 2 = 0.08896 and ω 2 = 0.37367.

Conclusions

The objective of this paper was to demonstrate that stable numerical evolutions of the Einstein equations on hyperboloidal slices extending to future null infinity I + can be obtained. While there has been extensive numerical work [START_REF] Frauendiener | Conformal infinity Living[END_REF][START_REF] Husa | Problems and successes in the numerical approach to the conformal field equations[END_REF][START_REF] Husa | Numerical relativity with the conformal field equations[END_REF] using the regular conformal field equations [START_REF] Friedrich | Cauchy problems for the conformal vacuum field equations in general relativity[END_REF] on such foliations, this is, as far as the author is aware, the first result using directly the Einstein equations in a dynamical situation. The advantage of this approach is that it is more similar to formulations of the Einstein equations on (truncated) Cauchy hypersurfaces that have been used successfully in numerical relativity for a long time.

In this paper, we assumed spacetime to be axisymmetric in order to be able to experiment with different ways of treating I + numerically in a relatively short computational time. The axisymmetric reduction of the Einstein equations we used (section 2) is similar to previous works [START_REF] Rinne | Constrained evolution in axisymmetry and the gravitational collapse of prolate Brill waves[END_REF][START_REF] Garfinkle | Numerical evolution of Brill waves[END_REF][START_REF] Choptuik | An axisymmetric gravitational collapse code[END_REF] but with a constant-mean-curvature slicing condition instead of maximal slicing.

The main new difficulty in adapting a formulation of the Einstein equations on Cauchy slices to hyperboloidal slices is the treatment of the formally singular terms in the evolution equations at I + . We applied two different methods [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] in order to obtain regular versions of these equations at I + (section 2.4). For the evolutions presented in this paper, we used Penrose's [START_REF] Penrose | Zero rest-mass fields including gravitation: Asymptotic behaviour[END_REF] result on the vanishing of the conformal Weyl tensor at I + . For this to hold, certain smoothness assumptions on the conformal metric at I + must be satisfied. An alternate method that does not require these assumptions is based on Taylor-expanding the fields near I + and evaluating the singular elliptic equations order by order; with this information, the apparently singular terms in the evolution equations can be evaluated (Appendix B). This alternate form of the evolution equations at I + was found to be stable as well and to yield comparable results to within numerical error. Future work should focus on this latter approach as it requires less restrictive smoothness assumptions and is most likely compatible with the so-called polylogarithmic behaviour at I + [START_REF] Andersson | On the regularity of solutions to the Yamabe equation and the existence of smooth hyperboloidal initial data for Einstein's field equations[END_REF].

Our numerical implementation (section 3) uses fairly standard techniques: fourthorder finite differences, the method of lines with a fourth-order Runge-Kutta scheme for the evolution equations, and a multigrid elliptic solver. The only non-standard modification to the multigrid solver is the use of a line relaxation in the radial direction in order to deal with the singular elliptic equations. With this method, no problems with the convergence of multigrid at the singular boundary were encountered. For the evolution equations, the usual Kreiss-Oliger dissipation [START_REF] Kreiss | Methods for the approximate solution of time dependent problems Global Atmospheric Research Programme[END_REF] sufficed to ensure stability; no special smoothing operations at I + were required.

We were able to evolve Schwarzschild spacetime for 1000M (M being the black hole mass) without any signs of instabilities and with nearly perfect fourth-order convergence (section 4). Next we included a small gravitational perturbation. In order to demonstrate the consistency and convergence of the scheme we monitored the momentum constraints during the evolution. We derived a regular expression for the Bondi news function at I + [START_REF] Bondi | Gravitational waves in general relativity VII. Waves from axi-symmetric isolated systems[END_REF][START_REF] Stewart | Numerical relativity III. The Bondi mass revisited[END_REF] representing the outgoing gravitational radiation. From this we computed the frequency and decay rate of the quasi-normal mode radiation and found good agreement with the tabulated values from linearized theory. The numerical resolutions we were able to reach were not yet high enough in order to study the subsequent polynomial decay (the "tail").

Currently the code is too slow in order to make runs with higher resolutions feasible (the run at the higher resolution in section 4.1 took about two weeks on a single processor). Not much effort has been spent on optimizing the computational efficiency of the code and there is certainly some room for improvement here. It seems likely, however, that a finite-difference method of higher than fourth order or a pseudospectral method will need to be used in order to reach the computational accuracy needed to study tails.

Although our multigrid solver has optimal computational complexity (linear in the number of unknowns) [START_REF] Brandt | Multilevel adaptive solutions to boundary value problems[END_REF][START_REF] Briggs | A Multigrid Tutorial[END_REF], the fact that our formulation contains elliptic equations is a major drawback on the computational speed. This is likely to be even more severe in the case without spacetime symmetries. From this point of view, a hyperbolic formulation such as the one by Zenginoglu [START_REF] Zenginoglu | Hyperbolodial evolution with the Einstein equations[END_REF] seems advantageous. It would be interesting to see whether the methods of [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF] can be applied to that formulation as well so that regular versions of the evolution equations at I + are obtained.

In this first study we imposed axisymmetry in order to be able to experiment with different numerical treatments of I + more quickly. The next step will be the extension to spacetimes without any symmetries. As far as the numerical stability of I + is concerned, we do not expect any new difficulties here because unlike spherical symmetry, axisymmetry already shares with the non-symmetric case what we believe are the key challenges, namely the existence of a dimension tangential to the conformal boundary and the presence of gravitational radiation. It remains to be seen which formulation of the Einstein equations will be more successful for this purpose. We expect the ideas and methods developed in this work to be applicable to a variety of formulations and gauge choices. = 2 contribution to the Bondi news function (49) at I + for an evolution using regularization at I + based on Taylor expansions (solid line) as compared with the approach based on Penrose's smoothness assumption (dashed line) that was used in section 4. The numerical resolution is (Nr, N θ ) = (128, 16). [START_REF] Briggs | A Multigrid Tutorial[END_REF], we obtain the following regular forms of the evolution equations at I + . W,t = e -2sη α c(s -1 η) ,θ + η ,θθ + η ,rr -s(η ,r ) 2 -η + α,θ ( c s η + η ,θ ) -(s -1 α,θ ) ,θ + e -sη α W,r + e -sη β θ [η ,r (cη -2 c s + sη ,θ ) -η ,rθ ], (B.9) Kr θ ,t = 0, (B.10) Ũ,t = - Figure B1 shows the Bondi news function at I + for a numerical evolution using the above Taylor-regularized evolution equations at I + . The result agrees well with the evolution based on the Penrose regularization technique that was discussed in section 2.4 and used for the evolution shown in figure 3. The discrepancy for times t 150 can be attributed to numerical error-figure 3 shows that the solution is still highly resolution-dependent at those late times. The corresponding curves in figures 1 and 2 virtually overlap for the two regularization techniques.

  where Γ C AB are the Christoffel symbols of H AB and Γ ˚C AB are the Christoffel symbols of the flat metric (H AB with η = 0). [The author thanks V. Moncrief for pointing this out.]

Figure 1 .

 1 Figure 1. Relative error of the numerical solution with respect to the exact solution for an evolution of Schwarzschild spacetime, computed as described in the text. Two different numerical resolutions are shown, (Nr, N θ ) = (64, 8) (dashed line) and (128, 16) (solid line).

Figure 2 .

 2 Figure 2. Momentum constraints (15)-(16) as functions of time for two different resolutions (Nr, N θ ) = (64, 8) (dashed lines) and (128, 16) (solid lines). The rmomentum constraint (left panel) has been normalized by the size of a typical term in[START_REF] Jr | Kinematics and dynamics of general relativity Sources of gravitational radiation[END_REF] as explained in the text.

Figure 3 .

 3 Figure 3. = 2 contribution to the Bondi news function (49), evaluated at I + , as a function of time. Numerical results for two different resolutions are shown, (Nr, N θ ) = (64, 8) (dashed line) and (128, 16) (solid line).

  Figure B1.= 2 contribution to the Bondi news function (49) at I + for an evolution using regularization at I + based on Taylor expansions (solid line) as compared with the approach based on Penrose's smoothness assumption (dashed line) that was used in section 4. The numerical resolution is (Nr, N θ ) = (128, 16).

1 2 s

 2 W,t . (B.11) Notice in particular how (B.10) and (B.11) are compatible with the regularity conditions (32).

Acknowledgment

I am grateful to Vincent Moncrief for numerous fruitful discussions. Helpful interactions with James Bardeen, Luisa Buchman, Harald Pfeiffer, Olivier Sarbach, John Stewart and Anıl Zenginoglu are much appreciated. I am supported by a Research Fellowship at King's College Cambridge.

Appendix A. Further details on the solution of the elliptic equations

Appendix A. [START_REF] Sarbach | Absorbing boundary conditions for Einstein's field equations[END_REF]

. Quasi-isotropic gauge conditions

The quasi-isotropic gauge conditions ( 6) and [START_REF] Deadman | Numerical relativity and asymptotic flatness[END_REF] are formally equivalent to the (inhomogeneous) Cauchy-Riemann equations. In order to solve them, one usually combines derivatives of these equations to obtain two decoupled Poisson equations, which can easily be solved. However, care must be taken so that the solution to the second-order Poisson equations is in fact also a solution to the original first-order Cauchy-Riemann equations (in general the two may differ by integration constants). We use the following approach.

The quasi-isotropic gauge conditions ( 6) and ( 7) are

Derivatives of these are combined to obtain a scalar Poisson equation,

This equation is solved for β θ using multigrid, with homogeneous Dirichlet boundary conditions. However, we do not solve the corresponding Poisson equation for β r . Instead, we observe that once β θ is known we can immediately integrate (A.1) and (A.2) to find β r . First, (A.2) is integrated along the outer boundary I + to determine β r there. The integration constant can be chosen arbitrarily, e.g. by requiring the average of β r along I + to agree with its value for the exact Schwarzschild solution. Next, (A.1) is integrated from each point on the outer boundary towards the interior. In this way we ensure that S r = 0 everywhere. We have also enforced S θ = 0 at the boundary. What remains to be shown is that S θ = 0 in the interior. But this follows immediately from (A.3): given that we already have S r = 0 everywhere, this equation reduces to S θ ,r = 0, and together with the initial condition S θ = 0 this implies S θ = 0 everywhere.

Apart from ensuring that the first-order conditions (A.1) and (A.2) are indeed satisfied, we also avoid solving another elliptic equation (for β r ) and thus save computational time.

Appendix A.2. Momentum constraints

For the intial data only, the momentum constraints are solved using a method similar to the one due to York [START_REF] Jr | Kinematics and dynamics of general relativity Sources of gravitational radiation[END_REF]. We introduce a vector V A and write the extrinsic curvature variables as

Note the similarity of (A.4) and (A.5) with ( 6) and [START_REF] Deadman | Numerical relativity and asymptotic flatness[END_REF]. The field Ŵ can be freely specified. The substitutions (A.4)-(A.6) completely decouple the momentum constraints ( 15) and ( 16) from the remaining elliptic equations:

We take V A to agree with the corresponding vector computed from the exact Schwarzschild solution at the radial boundaries, in particular we have V A = 0 at I + , and V θ also vanishes at the inner boundary. The angular parity of V A is identical with that of β A (table 1).

Appendix B. Alternate regularized evolution equations

In this appendix we derive alternate regular forms of the evolution equations at I + using the Taylor expansion method developed in [START_REF] Moncrief | Regularity of the Einstein equations at future null infinity[END_REF]. Within a hyperboloidal hypersurface t = const, all the fields are expanded in finite Taylor series about

These series are substituted in the singular elliptic equations and evaluated order by order. Again we have found a computer algebra programme [START_REF] Hearn | REDUCE User's Manual 3[END_REF] to be very helpful here. From the Hamiltonian constraint ( 14) we obtain

From the momentum constraints ( 15) and ( 16) we obtain, respectively,

Substituting the Taylor expansions (B.1) in the right-hand sides of the evolution equations ( 19)-( 21) and using the above results (B.2)-(B.8) for the radial derivatives of the fields, together with the regularity conditions [START_REF] Andersson | On the regularity of solutions to the Yamabe equation and the existence of smooth hyperboloidal initial data for Einstein's field equations[END_REF] and the boundary condition

Appendix C. Finite-difference operators

In order to compute derivatives with respect to r, we first compute finite-difference approximations to derivatives with respect to the map coordinate x and then apply the chain rule,

Here dx/dr and d 2 x/dr 2 are computed analytically from the given form of the mapping r(x), equation [START_REF] Brandt | Multilevel adaptive solutions to boundary value problems[END_REF]. Fourth-order accurate finite difference operators are used with respect to both x and θ. It suffices to display their one-dimensional forms here. The centred finite difference operators are

where h denotes the (uniform) grid spacing and the indices refer to the grid points where the quantities are evaluated. The forward differences are

and expressions for the backward differences follow by symmetry.

In the x-direction we use centred differences for 2 i N r -2, forward differences at i = 0, 1 and backward differences at i = N r , N r -1. In the θ-direction centred differences are used at all interior grid points 1 j N θ (note that ghost points are provided at j = -1, 0 and j = N θ + 1, N θ + 2).

The expression (s -1 u) ,θ is discretized by finite-differencing the function s -1 u with respect to θ. No such special differencing is used for (s -1 u ,θ ) ,θ , which we simply expand into s -1 u ,θθ -s -2 cu ,θ and then discretize in the standard way. Mixed derivatives u ,xθ are discretized by applying the finite-difference operators with respect to x and θ subsequently in the obvious way.

We use Kreiss-Oliger dissipation [START_REF] Kreiss | Methods for the approximate solution of time dependent problems Global Atmospheric Research Programme[END_REF], whereby the operator

is applied to each variable u both in the x and the θ direction and the results added to the discretized right-hand side of the evolution equation for u. Here D ± denote the first-order forward and backward finite-differencing operators, and a typical choice of the dissipation parameter is = 1/2. In the x-direction the above dissipation operator is added at grid points 2 i N r -2 and in the θ-direction at grid points 2 j N θ -1. Including j = 2, N θ -1 has been found to be essential in order to eliminate an angular instability.