Spectralm ethodsi n generalr elativity- towardt he simulationo f 3D-gravitationalc ollapseo f neutrons tars S. Bonazzola J. Frieben* E. Gourgoulhon J.A. Marck Abstract ca! investigation.I t is further stimulated by the prospects of gravitational wave astronomyw hich will turn into an Severala pplicationso f spectralm ethodst o problemsr e- observational science toward the end of this decade thanks lated to the relativistic astrophysicso f compacto bjects to gravitational wave observatoriesl ike LIGO, VIRGO and are presented.B asedo n a proper definitiono f the analyti- GEO600t hat are nowu nderc onstruction[4 , 5, 6]. cal propertieso f regular tensorialf unctionsw e have devel- We use the (3+1)-formalismo f generalr elativity [7] oped a spectral method in a generals pherical-likec oordi- which consistsi n foliating spacetimei nto a sequenceo f nate system. The applicationsi nclude the investigationo f space-likeh ypersurfacews hich representc urved3 -spacea t sphericallys ymmetric neutron star collapsea s well as the a fixed coordinatet ime t. The fabric of spacetimei s then solutiono f the coupled2 D-Einstein-Maxwell equationsf or determinedb y the 3-metrich ij and four additionaql uan- magnetized,r apidly rotating neutron stars. In both cases tities, the lapsef unctionN and the shift vectorN i which the resulting codes are efficient and give results typically fix the propagationo f the space-likeh ypersurfacesin time severalo rdersO f magnitudem ore accuratet han equivalent and the changeo f the spatial coordinates ystemb etween codesb asedo n finite differences chemesW. e further report adjacent hypersurfaces.T his Hamilton type approacht o the current status of a 3D-code aiming at the simulation generalr elativity resultsi n a temporal first order evolution of non-axisymmetric neutron star collapsew here we have schemef or the dynamicalv ariablesw hich is completedb y chosen a tensor based numerical scheme. somec onstrainte quationsw hich ensuret he consistencyo f gravitationala nd matter fields. FurthermoreN and N i Key words: spectral methods, numerical relativity, neu- have to be determinedb y the choiceo f appropriateg auge tron stars, black holes,g ravitational collapse. conditionsw hich typically lead to elliptic equationst hat AMS subject classifications: 65M70, 65N35, 83-08. have to be solveda t each time step. For stationary con- figurations all time derivatives vanish and one obtains a systemo f couplede lliptic equationsf or the gravitational I Introduction fields. The efficients olutiono f elliptic equationsi s hence of central interest for us. Compacto bjectsi n astrophysicss ucha s neutrons tarsa nd Letu sc onsidae rc ovariaPnto issoenq uatioNnl ili- $ blackh olesa re subjectedt o the strongf ield regimeo f grav- in a conformallyf lat axisymmetrics pacew here the line itation and have hence to be treated within the framework element reads of generalr elativity. The growingi nterest in the numerical solutiono f the Einstein equationsf or astrophysicallyr ele- (1) dl2 - A4(r,O)( dr2+r2 d O2+2rs in2d0(cid:127)b 2). vant systemsh as given rise to a new branch of computa- tionalp hysic-s numericarle lativity[ 1,2 , 3]. Thisd evelop- The formere quationc an be rewrittent o yield a Poisson- like equationfo r N wherew e havei solatedt he flat space ment is due to the increasinglyp owerfulc omputationarle - LaplacianA ! and contributedth e curvaturete rmst o the sourcesw hich make thesep roblemsa ccessibleto a numeri- source. Here a denotes In A. 'D(cid:127)partement d'Astrophysique Relativiste et de Cosraologie (UPR 176 du CNRS), Observatoired e Paris, Sectiond e Meudon, (2)A IN= $ with$= A45-2(OraO1 rN+(cid:127)OoaOoN). F-92195 Meudon Cede, France, e-mail: [rieben@obspm.[r This equation has to be solvedb y iteration. The solution ICOSAHOM'95: Proceedingso f the Third International Con- ofA IN---(cid:127) ate achit eratiohna sh encteo b ea ccomplished ferenceo n Spectral and High Order Methods. ((cid:127))1996 Houston sufficientlfya sti n ordert o keept he total computationc ost Journal of Mathematics, University of Houston. at a reasonable level. 4 ICOSAHOM 95 After outliningt he basicf eatureso f our spectralm ethod x sini+J8c os(cid:127)8c osi4s'i n(cid:127)4'. [8, 9], we will proceedin a first stept o the investigation of black hole formation due to sphericallys ymmetricn eu- Having specifiedt he tensor componentsfi (cid:127)...iN with re- tron star collapsew hich has provedt he high aptitude of spect to the Cartesian frame we derive the components spectralm ethodsin this field [10, 11, 12]. The second related to the local orthonormalf rame of sphericalc oordi- part is devotedt o the study of axisymmetrics tationary nates by a non-singular coordinatet ransformation. rotating bodiesw hich has been applied to model rapidly In order to infer the analytical properties of a scalar rotatingn eutrons tars[ 13, 14]. This workh asb eene x- functioni t is usefult o rearranget he sum in (4). We first tended recently to include strong magneticf ields for the collecta ll the terms referringt o cosm 4' and sinm 4' respec- firstt ime into neutrons tar models[1 5].S peciael mphasis tively. We write in all casesh as been put on the extensiveu se of external M and intrinsict ests[ 16, 17, 18] of the self-consistencayn d (5) f(r,8 ,4'=) (cid:127) (a,,(r8, )c osrn4b',+(cid:127)( r,8 )s inm4') the attained accuracyo f the numerical results. The result- rn(cid:127)-O ing neutron star modelsp rovideu s with the requiredi nitial value modelsf or the investigationo f 3D-gravitational col- wherea m(r, 8) and bin(r,8 ) behavei denticallyi n the fur- lapse of neutron stars which will reveal the whole range ther procedure. We opt for cos1 8 and sin 18 as basisf unc- of gravitational wave emission associatedw ith this phe- tions in 8 which allowst he application of FFT-techniques nomenon. We give an overview about the inset of spectral for this transformation. An immediate conclusion from methodsi n this project which is currently in work. Here (4) and the casei -t-j eveni s that the coefficientas2 ,(cid:127)(r, 8) a new method for the efficient inversion of a generalized and b2,(cid:127) (r, 8) havet o be expandedi n terms of cos1 8w hile 3D-vector Poissone quation is a first major result. from the odd caset hat the expansiono f a2,(cid:127)+(cid:127) (r, 8) and b2,(cid:127)+i (r, 8) hast o be doneo n the set sin1 8. We therefore specify 2 Spectral methods in general rel- L ativity ((3) = coslO, l(cid:127). O 2.1 Coordinates and regularity conditions L (?) a2m+l8(r)= , (cid:127) &(cid:127),(cid:127),(cid:127)+(cid:127)(r) sin18. The space-likeh ypersurfacess temming from the former l=0 choiceo f the (3+1)-formalisma re conceivedto describe Note that due to the well definedp arity of a2,(cid:127) and some asymptotically flat space, containing a compact, thesec oefficient-s a priori only definedf or 0 < 8 < 7r- can mostly starlike object. The natural choicei s thus a sphe- be continueda nalytically to periodic functionso f 8 on the ricallikec oordinates ystem( r, 0, 4'). intervalI -w, The pseudo-singularitiews hich appear in this case can In the samem annera s beforew e find from (4) that the be overcomeb y a proper definition of regularity conditions polynomials& l(r) are symmetricw ith respectt o the inver- of tensorial quantities. A consequenat pplicationo f parity sionr -* -r for I evena nda ntisymmetrfiocr I odd] As rules derived from these conditionsa llows further to opti- basisf unction set in r we decidef or Chebyshevp olynomi- mize code efficiencya nd precision. als due to their superiorp ropertiesi n finite approximation We considert he related Cartesian type coordinate sys- schemeso f non-periodic functions and the availability of tem (x, y, z) = (r sin8 cos4 ',r sin8 sin4 ',r cos8 ). We define fast Chebyshevtr ansforms.S implificationso f the expan- a tensoriaql uantityT i(cid:127) ...iNi n sphericacl oordinate(sr, 8, 4') sion schemei n the presenceo f additional symmetriesa re to be regular,i f its componentsfi (cid:127)...iN with respectt o preeisedi n the followingp aragraphs. Cartesianc oordinate(sx , y, z) are regulari n the senseth at Any regularf unctiona dmits an expansiono f this kind, they can be expandedi nto a polynomials um of the type but regularity furthermorei mplies additional constraints N on the differentc oefficientsH. aving set up regulari nitial (3) data (cid:127)eeordingt o (3), regularityo f the involvedq uantities i,j,k=O is maintainedd uring a calculationb y the applicationo f whichc anb e written in termso f (r, 8, 4') as regularo perators- we herei gnoret he influenceo f numer- N ical effectsd ue to aliasingo r roundoffe rrors. (4) Let us eousider the eovariant derivative of a vector in sphericacl oordinatesF. or the scalarp otentialU (r, 8, 4')- i,j,k----O SpectralM ethodsI n GeneralR elativity 5 r sin (cid:127) cos whichc anb em odifiedre, placinsgi n20b y (1-cos20). L 21 M (s) (16) /+(,., EE E !(cid:127)0 k----O m----O we have x (1-cos20t )c os:(cid:127)Oc o2st -k(cid:127)bs ink (9) (Ur, Us,U (cid:127))'- (sin0 cos(cid:127)b , cosO cos(cid:127)b ,- sin(cid:127)b ). L 21 M (17) The covariandt erivativeU sl(cid:127) whichr eads !--0 k----O m--O x (1 - cos0:) (cid:127) cos:+* (cid:127) (cid:127) co2s(cid:127) - k (cid:127)b sink 1 1 (10) U01=(cid:127) rsinOO (cid:127) Us U(cid:127) rtan(cid:127) According to Sec. 2.1 we concludef or the decomposition can be rewritten to yield of a supersymmetricf unction: 1 1. f is (cid:127)r-periodic in the angular variable (cid:127)b, (11) UO-l-O-r- s in(cid:127) (OU(cid:127) o- cosUO(cid:127)) ß 2. a2,(cid:127) (r, O)h ast o be expandedin to a sumo f cosl Ow here A numericale valuationa ccordingto (8) and (11) reveals I is even in the symmetric case and odd in the anti- a perfectly regular behavior. symmetric one, 2.2 Supersymmetric case 3. 5t( r) is an evenp olynomiali n r for f symmetrica nd an odd polynomialf or f antisymmetric. Additional spatial symmetrieso f physicals ystemsi nvolve continuouss ymmetry operationsl ike rotation about a dis- Physicalp roblemsw hich imply the use of supersymmetric tinct axis with an associatedK illing vector field or discrete functions allow to restrict the computational domain to transformationsl ike inversiona t the equatorial plane z = 0, [0,(cid:127)r] in (cid:127)b and to [0,(cid:127)r /2] in 0 which leadst o an overall leadingt o distinct parity propertieso f the different tensor reductiono f the effectiveg rid size by a factor four. components. Fort hec omponenotfs a vectofri elda ssociatewdi tht he We defines upersymmet[r1y9 ]b y the followingb ehavior case of even supersymmetryw e conclude: U(cid:127) can be ex- of a scalar function: f is invariant with respect to inver- panded into a sum of cos2 10 , Uo into a sum of sin 210 siona t the z-axis, hencef (-x,-y,z)= f(x,y, z), whilef and U(cid:127) into a sum of sin(2/+ 1)0 whichm eanst hat Uo is [anti-]s ymmetriwc ith respectto reflectionat the equa- undergoesa changeo f sign by reflectiona t the equatorial torial plane z = 0, hencef ((cid:127), y, -z) = :El(z, y, z). Conse- planew hileU r and U(cid:127) remainu nchangedT. he expansion quently f can be expandedi nto a sum of the type of the radial part is done in terms of odd powerso f r for all components. L 2l M l-'O k--O m----O 2.3 Axisymmetric case L 21 M Axisymmetry restricts the set of scalar functions under (13f)_ (x,y,z-)' E E E ck(cid:127)(cid:127)x2(cid:127)-(cid:127)Y(cid:127)z2(cid:127)+(cid:127) oddc ase considerationto functionsa ccordingt o Sec. 2.1 where the 1=0 k(cid:127)O m=O only remainingte rm in (5) is ao(r, 0). Thus f can be writ- which definess ubsetso f the generals calarf unctionsi ntro- ten as L ducedi n Sec.2 .1. A write-up in termso f (r, 0, (cid:127)b) gives (18) f(r,0 )= E ai,o(rc)o s/0 L 21 M l----0 (14) = E E E where(cid:127)l ,o(r) is an evenf unctioni n r for I evena nd an odd l----O k----O m----O function for I odd. x sin2(cid:127)Oc os2'ncOo s2(cid:127)-k(cid:127)bs ink(cid:127)b, For the componentUsr and Uoo f a vectorf ield compat- L 21 M ible .witht he assumptiono f axisymmetryw e concludeU: r (15) /_(,-o, ,,)= E E hast o be expandedin termso f cosl O wheret he radial part 1----0 k----O m----O is evenf or I odd and odd for I even. Uo hast o be expanded x sin2(cid:127) cos2m+(cid:127)oco s:(cid:127)-k$ in terms of sinl O wherew e alsoh avea parity changei n r. 6 ICOSAHOM 95 The propertieso f the lackingc omponenUt (cid:127) whichw e Concerningt he hydrodynamicalp art we employ a set needt o handleN (cid:127) in the modelso f rotatingn eutrons tars of particularv ariablesw hich lead to equationsr essembling can be derivedf rom the Killing equationl inked to axisym- very closelyt heir Newtonian counterpartsi ncludinga gen- merry. A shorte xaminationre vealst hat U(cid:127) has to be eral relativistic generalizationo f the classicalE uler equa- expandedin terms of sinl O with a parity changei n r - it tion. We note the privilegedr ole of the Eulerian or local behavesi denticallya s Us. rest observerO 0 in this formulation. The hydrodynamical equationsr ead 2.4 Spherically symmetric case 1 (21) Otqœ-(cid:127)'(cid:127) O((cid:127) rS(e+p(cid:127)) )-V-O , For sakeo f completenesws e add the sphericallys ymmetric case. As a further restriction of the axisymmetricc asew e (22O)U , q-V(cid:127)O- (cid:127)U I (N keepf rom (18) only a0,0(r). f(r) thereforere ads œ+p(cid:127)O (cid:127)p+ uotp K (19) f(r) = (cid:127) asrks k Fs q-4 (cid:127)rrp, 1 (2a) O,D+ ;(cid:127) O(cid:127)(rSDW= )0 , which is an ordinary even polynomial in r while U(cid:127) is rep- resentedb y an odd polynomial in r. (24) Otsr + V(cid:127)O(cid:127)sr = 0 where we have introduced the energy density (cid:127) and the 3 Spherically symmetric neutron fluid velocity U measured by O0, the coordinate baryon star collapse density D and the entropy per baryon ss while V r de- notes the fluid coordinatev elocity. The Lorentz factor F 3.1 Basic equations isd efineads F = (1- US)-«.I n additioonn eh asto s olve The investigationo f neutron star equilibrium configura- (25) Om(r,t) = sœ (r,t), tions in sphericals ymmetry was the first problem in the fully general relativistic regime being solved by us by (26O) rff9A((rsm, t)(=r 't) Us) r(cid:127) q-4 7 rr [ p q -(e q P-) ] meanso f the (3+1)-formalism of generalr elativity and a spectraml ethod[ 10,8 , 9]. A favorablec hoiceo f the linee l- whereA (r, t) is relatedt o m(r, t) through emendt ss = 9,(cid:127) dx(cid:127) dx(cid:127) in thec aseo f sphericsayl mmetry is given by RGPS (Radial Gauge-PolarS licing) coordi- (27) A(r- ,(t1)2 ra)(- r«' t) nates[ 20]a nd reads (20) dss = -NSdts + ASdsr + r s( dOSs+i nS0bdse).. The neutron star matter is modeled as a perfect fluid, We stresst he particular nature of this problemw here the adoptinga realisticd ensem atter equationo f state. field variableA is not really a dynamicalq uantity. In fact it is uniquelyd etermineda t any momenta s well as the lapse 3.2 Numerical method functionN by the matter fieldsw hichh avet o be evolvedb y meanso f the hydrodynamicael quations.S incet he solution The initial value model for the dynamical calculations outside the star is known in advance to coincide with the is providedb y solvingt he Tolman-Oppenheimer-Volkoff static Schwarzschilds olutiono f a point masso f the same equations describinga spherically symmetric static star size accordingt o the Birkhoff theorem,w e benefit from a where each model is determinedb y the central value of doubles implification.F irst the wholet ime evolutioni s yet the pseudo-enthalpyH c. A first solutioni s obtainedb y in- determinedb y propagatingt he hydrodynamicavl ariables tegrationo f this systemo f ordinaryd ifferentiale quations and further the calculation can be restricted to the stellar whilea n overalnl umericaalc curacoy f the ordero f 10- (cid:127)4, interior. The interior solutionf or the gravitationalf ield adapted to the subsequenut se of a spectralm ethod, is has then to be matched to the analytical exterior one. We achievedb y iteration of the approximates olution. further note that due to the static character of the exterior The computationald omain is identical with the stellar solutionn o gravitationalw aves- whicho therwisew ouldb e interiord uringt he wholec alculation,t hankst o a comoving an observableo f most importance- are emitted. grid whoseo uter boundaryc oincidesw ith the star surface. SpectralM ethods In General Relativity 7 i (cid:127) I 0 2 4 6 0.2 0.4 0.6 0.8 1 (cid:127) / R.(0) Figure 1: Relative variation of the stellarr adiusw ith elaps- Figure2 : Profileso f the lapsef unctionN for t rangingf rom ingt imef ora stablee quilibriummo dewl ithn B<n (cid:127) rit and 0 to 7.296 ms. The locationo f the Schwarzschildra dius Rs [AM/Mmax[ = -5 x 10- 5. The hydrodynamictaiml escale is indicatedb y a verticalb ar. Directiono f increasingti me is indicated on the lower left. is downward. This maintains a constant spatial resolutiond uring the col- the beginning of the collapseu ntil the moment where the lapsea nd a fine samplingo f the steepeningg radientsn ear whole configurationi s highly relativistic and practically the star surface thanks to the accumulation of the Gaufi- frozenw e observedth e differentv ariablesto reproduceth e Lobatto points in this region. It furthermore minimizes analyticavla luesw ithina n erroro f bettert han 10- 5 [10], the advectivet erms, hencei mproving the numericala ccu- while the errors related to previouss tudiesb asedo n finite racy. All quantities are expanded in terms of Chebyshev differencmee thodasr eo f theo rdero f 10- 5 [23]. polynomialisn r mappingth e interval[ 0,R ,(t)] onto[ 0,1 ] Neutron star modelsn ear the maximumm assc onfig- taking into accountt heir analytical propertiesa ccordingt o uration - Mmax- 1.924Mo and R -- 10.678 km for the Sec.2 .4. The 2n do rders emi-implictiitm ei ntegratioenn - employedE OS - are interestingw ith respectt o their sta- surest he stability of the code which allowsu s to perform bility againstr adial perturbations.C onfigurationsw ith a simulationso f a duration of many dynamical timescales. centralb aryond ensitya pproachingth e critical one exhibit This ability is very important when studyingt he effects an oscillatoryb ehaviorw hichi s dominatedb y the funda- of perturbationso n equilibrium configurations.T his task mentalm odeo f oscillationo f increasingp eriodl ength. It is also favoredb y the fact that we integratet he original representsa uniform growinga nd shrinkingo f the entire systemo f equationsw ithout any artificial viscosityt o sta- star which is modulated by less important harmonicso f bilize the code while in addition the intrinsic viscosityo f higher order. Fig. 1 showst his temporal variation for a spectral methods is negligible. The ingoing characteristic stable configuration.W e stresst hat these oscillationsa re of the hydrodynamicals ystem at r--R, gives rise to one entirely driven by roundoff and discretization errors of a boundary condition which is chosent o fix the baryon den- totalo rdero f 10- (cid:127)ø - no externaflo rceh asb eena pplied sity at the star boundary. It is imposed by means of a to trigger this variability. Increasingt he central baryon (cid:127)-Lanczos procedure[2 1]o n the systema s a wholew hich densityb eyondt he critical value one enterst he branch of is the well posed mathematical procedure. unstablec onfigurations.T his time the fundamental mode starts a contraction of the star which results in an unlim- ited collapse.T hough the actual coordinatec hoicei s not 3.3 Results capablet o properlyc apturet he formationo f an appar- Among the varioust ests we have imposedo n our codeo ne ent horizon which would clearly reveal the formation of a describesth e collapseo f a homogeneoudsu st spherew hose black hole the evolutiono f metric potentialsa nd matter solutionw asg ivenb y Oppenheimearn dS nyder[ 22].F rom variablesg ivesa distinct indication of this event. Take for $ ICOSAHOM 95 densitiesi n the individual subdomainsw hich appearsd ur- ing the collapse.F our to five zonesa re typically needed to covert he densec entral core forming the new born neu- tron star and the outer layers of much lower density with the desiredr esolution. In contrastw ith the original code boundary conditionsa re imposedi n this caseb y meanso f a modified(cid:127)- -Lanczos schemea pplied in coefficienst pace. 4 Axisymmetric rotating relativi- stic bodies 4.1 Basic equations In order to study rapidly rotating neutron stars we have made the assumptionso f spacetimeb eing axisymmetric, t [msl stationarya nd asymptoticall.yf lat. We have further sup- poseds pacetimet o be circular, thus the absenceo f merid- Figure 3: Relative variation of the total baryon number ional currentsi n the sourceso f the gravitational field. In during the collapse. this cases pacetimec an be describedb y MSQI (Maximal Slicing-QuasIsi otropicc) oordinate[s1 3]w hichh avea line elemendt s2 = g(cid:127) dx(cid:127) dx(cid:127) of the form: instanceF ig. 2 which showst he time developmento f the lapse function N, measuringt he elapsedp roper time of (28) ds(cid:127) = -N(cid:127)dt (cid:127) + A4B- (cid:127) (drY+r 2dO(cid:127)) the local Eulerian observer Oo. Inside the Schwarzschild +A4B2r(cid:127) sin20(d 4-(cid:127) N(cid:127)dt) (cid:127). radiusi t tendst owardz erof or increasingc oordinatet ime t. This behavior- calledt he 'lapseo f the lapse'- is relatedt o Spacetimeis hencef ully determinedb y the four metricp o- the singularitya voidancep ropertyo f the chosenc oordinate tentialsN , N (cid:127), A and B. MSQI-coordinateasr e global gaugew hich keepst he spatial hypersurfacesfr om propa- coordinatesa nd lead to elliptic operatorsw hich admit a gatingi nto a formings ingularity.T he dynamicalt imescale consistentt reatment of boundary conditions and an effi- of the collapsei s the time elapsedf rom the beginningo f cient solutionb y spectralm ethods. The matter fields are the collapseu ntil the momentw heret he evolutiona ppears chosento model a perfect fluid where first a polytropic, to be frozen to a distant observer and has the value t = 7.4 hencea nalytical,e quationo f state wasu sed. The assump- ms for the consideredc onfiguration. Fig. 3 illustrates the tion of a perfectf luid reducest he equationso f motion to relativee rror committedo n the total baryonn umberw hich an algebraice quationf or the heat function H. From the is a conservedg lobalq uantity. In the early phasei t is con- Einstein equationso ne derivesf our elliptic equationsf or servedw ith a relativea ccuracoy f 2x10- s andw ith 4x10- 6 the variablesv = In N, N (cid:127) and during the violent stageso f the collapse. The deviation increaseusp to 5 x 10- 5 in the final phasew heres harp (29) G(r,O)-- N(r,O)A2(r,O)B(r,O), gradients form near the horizon. The proper working of the codei n the perturbativer egimeh as been recentlyc on- (30) ½(r,0 ) = v(r, + 2(r, - firreed by direct comparisonw ith a linear adiabatic code whose final form reads [12]. The calculationhsa ves howna veryg ooda greement of the frequencieso f the fundamental modes for the two A 4 (31)A sy - - (cid:127)-(cid:127)[4(cid:127)r(E-FSii)-F2(k(cid:127)-Fk(cid:127))]-O(cid:127)O((cid:127)-F2a+(cid:127)), oppositea pproachesT. he originalv ersiono f our codeh as beene xtendedl ater on to includet he processeosf neutrino (32(cid:127)3)(cid:127)- -167rB7N J-(cid:127) ---r sint(cid:127)ON(cid:127)O(6e(cid:127)-F3(cid:127)-(cid:127)), productiona nd transport during neutron star collapsein r sin order to compute the observablen eutrino emissionf or a distanto bserver[1 1]. A multi-domaine xtensiono f this _ NA (cid:127) (33) code is in use in order to simulate type II supernovaein sphericals ymmetry.I t is particularlyw ell suitedt o prop- A 4 erly captureth e highc ontrasut p to _(cid:127) 106o f the mean (34) Spectral .(cid:127)letl(cid:127)ods In ((cid:127)eneral Relativity 9 where we have introduced (cid:127)=ln A, (cid:127)=ln B and polynomials. The inner zone is mapped onto half the def- initioni nterval[ 0,1 ] whicha llowst o take into accountt he (35) (cid:127)(r, O)= r sinO((cid:127)(r,O), parity propertieso f regular functionsw ith respectt o the origin accordingt o Sec. 2.3. In the compactifiedz one the (36) J(cid:127)r(cid:127)(r,O) = rsin0N (cid:127)(r,O) expansionh as been done in the usual manner on the whole as well as the abridged notation definitionin terval[ -1, 1] of Chebyshepvo lynomials. The effective scheme which is based on a relaxation 1 methodw orksa s follows.W e considerr igidly rotating neu- (37) = + ao. tron starsw ith a polytropice quationo f state. A particular configurationis henced eterminedb y fixing the value H(cid:127) Further employedq uantitiesa re the total energyd ensity of the heat function at the centre of the star and its angu- E, thes trestse nsoSr ij, them omentudme nsityJ i andk (cid:127), lar velocity (cid:127). We start from very crude initial conditions ks whicha re relatedt o the extrinsicc urvaturete nsorK (cid:127)j. where all the metric quantities are set to their flat space A2, Asa nda s denotsec alaLr aplacianins t woa ndt hree values( c,, (cid:127), (cid:127), ½ and N (cid:127) = 0; G = 1) and the matter dimensionsa nd a vector Laplacian in three dimensionsr e- distribution is determinedb y a first approximate guess. spectively. While the GRV2 identity related to (34) holds for an Existencea nd uniquenesso f the solutiono f thesee lliptic exact solutiono f the Einstein equationsw e have to enforce equationasr e ensuredfo r physicallyre levanct ases[2 4,2 5, this consistencyr elation at each iteration step in order 20]. to avoid a logarithmicd ivergenceo f the approximateo ne. Note that (33) and (34) can be continueda nalytically This can be accomplishebdy modifying( 34) accordingt o to yieldg enuine2 D-Poissone quationsin the entire( r, 0)- plane. One infers immediately that (cid:127) exhibits a logarith- (38) A2 mic divergencef or r (cid:127) o(cid:127) unless the total integral over the sourceo f (34) vanishesid enticallyw hereasw e require (39c)r= (cid:127) 8 (cid:127)rBSA(cid:127)'(cid:127),-, --o(cid:127) '=(cid:127) (cid:127)-A -4 (cid:127).[3(k(cid:127)-.F((k:%(cid:127))]-, )2. (cid:127)[ r=(cid:127) = 0. This 2D-virial theoremo f generalr elativity {GHV2) [16]i s in somes ensere latedt o the classicaNl ew- tonian virial theorem and furnishesa consistencyc ondition At each iteration step A is choseni n such a way that the for any solutiono f the Einstein equationsw hich is compat- total sourcei ntegral vanishes. The final solution has to ible with our basic assumptions. It has to be taken into satisfy( 34) exactlyw hichi s equivalentto A= 1. The de- accountd uring the calculation and providesa strong con- viation of A from unity during the iteration measurest he sistencyc hecko f the numericals olution. violation of self-consistencoy f the approximates olution and can be usedt o monitor convergence. The sourceso f (31)-(34) exhibit somet erms involving 4.2 Numerical method simpleo peratorsl ike r, sinS, 0(cid:127), etc. which are accurately The mathematicapl roblemi nvolves( 31)-(34) and an al- computedi n coefficienst paceb eforee valuatingt he entire gebraicf irst integral equationf or the matter fields. Our expressionisn configurations pace.E xpansiono f the total numerical solutions are exact in the senset hat the gov- sourcesin terms of the angulare igenfunctionso f the differ- erninge quationsa re derivedf rom the full theory of gen- entL aplacian(sP (cid:127)ø( cos0 ), P(cid:127) (sin0 ) and( cos1 8,s in(cid:127) 8) for eral relativity without any analytical approximationw hile As,a s, andA (cid:127) respectivelelya)d tso a systemof O DEsin the numerical code solvest hese equationsi n all spacee x- the radial variable. The uniqueg lobal solutioni s obtained tendingt he numericali ntegrationt o spatial infinity which by appropriatel inear combinationso f the corresponding allowst o imposet he exact boundary conditionso f asymp- particular and homogeneouss olutionsi n order to match totical flatness on the gravitational fields as well as the the piecewises olutionsa t the grid interfacea nd to satisfy- properc alculationo f the sourcete rmso f (31)-(34) which the boundary conditionsa t r-o(cid:127). The new values of the fill all space. gravitationalf ield variablesa re then usedt o update the This is accomplishedb y the use of two grids, where the matter distributionb y meanso f the first integral equation first one coverst he stellar interior usingt he radial variable and the iteration can go on. r in thei nterva[l 0,R ], whilet heo uters paceis c ompactified thankst o the variablet ransformu = r -(cid:127) and in this way 4.3 Tests mapping[R ,o (cid:127)] ontot hef initei nterva[lR - (cid:127), 0]. Whilei n the O-variable a Fourier expansiona ccordingt o Sec. 2.1 We have subjectedo ur codet o two different kinds of tests is used,t he radial part is expandedi n terms of Chebyshev which ensuret he reliability of the numerical results.