Nipoti, Carlo

Italian

Almae Matris Studiorum Campus

Calcolo con Laboratorio II - a.a. 2007/2008 115/5/20083. Modelli numerici di galassie[Riferimenti bibliografici: Galactic Dynamics (Binney & Tremaine), NumericalRecipes (Press et al.)]Per studiare l’evoluzione temporale di un sistema a N-corpi e` necessario costruirele condizioni iniziali, cioe` assegnare posizioni, velocita` e masse delle particelle altempo t = 0. Consideriamo, per esempio, di voler rappresentare una galassia inequilibrio di massa totaleM con N particelle, ciascuna con massaM/N . Un modellodi galassia in equilibrio e` definito completamente dalla funzione di distribuzionef(x,v) nello spazio delle fasi. La f e` la massa per unita` di volume d3xd3v nellospazio delle fasi. Quindi la massa totale del sistema e`M =∫fd3xd3v (3.2)e la quantita` f(x,v)/M e` la probabilita` di trovare una stella nel volume infinitesimodi spazio delle fasi d3xd3v centrato in (x,v). Un modello a N-corpi e` unarealizzazione numerica della funzione di disribuzione. Quindi, per ottenere ilmodello a N -corpi dobbiamo generare per ogni particella 6 coordinate (x,v) condistribuzione di probabilita` f(x,v)/M . Consideriamo il caso di un modello sferico,che puo` rappresentare una galassia ellittica poco schiacciata (E0-E1) o un ammassoglobulare. Per semplicita` ci limitiamo al caso di funzione di distribuzione isotropaf(x,v) = f(E), doveE =12v2 + φ(x) (3.3)e` l’energia per unita` di massa della particella e φ e` il potenziale gravitazionale. Ilmodello ha distribuzione di densita`ρ(r) =∫fd3v (3.4)e profilo di massaM(r) = 4pi∫ r0r′2ρ(r′)dr′. (3.5)Nel caso di un sistema autogravitante con distribuzione di densita` ρ il potenzialegravitazionale φ e` dato dall’equazione di Poisson∇2φ = 4piGρ. (3.6)Generiamo prima le 3N coordinate spaziali in modo da riprodurre il profilo didensita` del sistema. Sia u un generico numero reale nell’intervallo 0 < u < 12 L.M. Atsrofisica e Cosmologia - Bolognagenerato con un generatore di numeri pseudocasuali (u e` detto uniform deviate).Se la distribuzione cumulativa di massa M(r) e` nota e facilmente invertibile (ilche accade in molti casi tipicamente considerati), l’assegnazione della coordinataradiale si puo` eseguire con il metodo della trasformazione (transformation method;si veda Numerical recipes, capitolo 7). La coordinata radiale ri di ogni particella siottiene invertendo la relazione u = M(ri)/M(rmax), dove rmax e` il raggio massimodella distribuzione. Le coordinate angolari ϕi e ϑi si ottengono semplicemente comeϕi = 2piu e ϑi = arccos[2u − 1]. Per l’assegnazione della coordinata ϑ si richiedeche la distribuzione sia uniforme in cosϑ perche´ l’elemento di volume in coordinatesferiche e` proporzionale a sinϑdϑ. Tipicamente un codice a N-corpi richiede leposizioni e le velocita` delle particelle in coordinate cartesiane, quindi trasformiamo(ri, ϕi, ϑi) in (xi, yi, zi).Nota la funzione di distribuzione (analitica o numerica), le componenti dellavelocita` delle particelle si possono ottenere con il metodo del rifiuto (rejectionmethod) di Von Neumann [si veda Numerical Recipes, capitolo 7]. Il metodo delrifiuto e` un metodo standard per generare numeri secondo una data funzione didistribuzione p(x) per la quale non sia calcolabile in modo semplice l’inversa P−1(x)della funzione di distribuzione cumulativaP (x) =∫ x0p(x′)dx′. (3.7)In generale, la funzione di distribuzione cumulativa per modelli di galassie none` facilmente invertibile, quindi il metodo del rifiuto e` molto utilizzato. Notoil potenziale gravitazionale del sistema φ(r), e` immediato calcolare il potenzialegravitazionale φi ≡ φ(ri) di una particella con raggio ri. Si generino quindi terne dinumeri casuali (tra -1 e 1) u1,i = 2u− 1, u2,i = 2u − 1, u3,i = 2u− 1, finche` non siottiene una terna tale che u21,i+u22,i+u23,i < 1. Tre candidate componenti di velocita`sono quindi vx,i =√−2φiu1,i, vy,i =√−2φiu2,i, vz,i =√−2φiu3,i (non prendiamoin considerazione terne con u21,i + u22,i + u23,i ≥ 1 perche` in tal caso la particella nonsarebbe gravitazionalmente legata al sistema, avrebbe cioe` energia positiva o nullaEi = v2i /2 + φi ≥ 0). Per decidere se accettare o meno le candidate componentidella velocita` generiamo un altro numero pseudocasuale u e lo confrontiamo conf(Ei)/fmax,i, dove fmax,i = f(φi): se u < f(Ei)/fmax,i accetto la velocita`, seu > f(Ei)/fmax,i rifiuto e riparto generando una nuova terna con u21,i+u22,i+u23,i < 1.Spesso si vuole generare un modello con distribuzione di densita` data, ma dicui non e` nota la funzione di distribuizione (un tipico esempio e` quello delle galassiea disco). In tal caso e` possibile utilizzare il metodo detto di approssimazionemaxwelliana locale (local Maxwellian approximation), basato sull’ipotesi che ladistribuzione di velocita` in un punto sia maxwelliana. Modelli generati con ilmetodo di approssimazione maxwelliana non sono modelli di equilibrio, ma di “quasiequilibrio”, nel senso che ci si aspetta che in pochi tempi dinamici raggiunganouna configurazione di equilibrio con distribuzione simile a quella assegnata. Nelmetodo dell’approssimazione maxwelliana, le velocita` di una particella di assegnatecoordinate spaziali xi vengono assegnate con una distribuzione Maxwelliana convarianza data dalla dispersione di velocita` nel punto xi. Le componenti del tensoredispersione di velocita` nel punto xi si calcolando risovendo le equazioni di Jeans. Nelcaso di un modello sferico isotropo, la componente radiale del tensore dispersione diCalcolo con Laboratorio II - a.a. 2007/2008 3velocita` si ottiene dall’intergrazione dell’equazione di Jeansdρ(r)σ2r(r)dr= −ρ(r)GM(r)r2. (3.8)Dato che ρσ2r = 0 per r =∞, la dispersione di velocita` al raggio r e`σ2r (r) =1ρ(r)∫ ∞rρ(x)GM(x)x2dx. (3.9)Poiche´ il modello e` isotropo, ciascuna componente cartesiana della velocita` hadispersione di velocita` σ(x) = σr(x). Quindi le componenti vx,i, vy,i, vz,i dellaparticella nel punto xi hanno una distribuzione gaussiana con media nulla evarianza σi ≡ σ(xi). Per generare le componenti cartesiane della velocita` dell’i-esima particella, si generano tre numeri pseudocasuali w1,i, w2,i, w3,i nell’intervallo[−∞,∞] con distribuzione gaussiana normalizzata (cioe` con media nulla e varianzaunitaria)p(w)dw =1√2pie−w2/2dw. (3.10)Si continuano a generare terne fino ad ottenerene una che soddisfi la condizionew ≡√w21,i + w22,i + w23,i <√−2φi/σi (particella legata). Quindi le componenti dellavelocita` si ottengono come vx,i = σiw1,i, vy,i = σiw2,i, vz,i = σiw3,i. Per generare unavariabile pseudocasuale con distribuzione gaussiana normalizzata (equazione 3.10)si puo` usare il metodo di Box-Muller (si veda Numerical recipes, capitolo 7).Nel metodo di Box-Muller si generano coppie di numeri psedocasuali u1 e u2con distribuzione uniforme nell’intervallo [−1, 1] fino ad ottenere una coppia conR = u21+ u22< 1. Da questi si possono ottenere coppie di numeri ciascuno condistribuzione gaussiana normalizzata w1 e w2 con le seguenti trasformazioniw1 = u1√−2 lnRR, (3.11)w2 = u2√−2 lnRR. (3.12)4 L.M. Atsrofisica e Cosmologia - Bologna3.1 Laboratorio: esercizio CUno dei modelli piu` semplici di sistemi stellari autogravitanti e` la sfera di Plummerisotropa, con funzione di distribuzionef(E) ={A(−E) 72 E < 00 E ≥ 0, (3.13)distribuzione di densita`ρ(r) =3M4pib31(1 + r2b2)52, (3.14)potenziale gravitazionaleφ(r) = −GMb1(1 + r2b2)12(3.15)e profilo di massaM(r) =M(r2r2 + b2) 32, (3.16)dove b e` il raggio di scala e M e` la massa totale.1. Si generi un modello a N-corpi (con N particelle di uguale massa) di una sferadi Plummer isotropa utilizzando la funzione di distribuzione con il metodo delrifiuto di von Neumann.2. Si generi un modello a N-corpi (con N particelle di uguale massa) di una sferadi Plummer isotropa utilizzando il metodo dell’approssimazione maxwellianalocale. [Suggerimento: si calcoli la dispersione di velocita` per n valori del raggio0 < rk < rmax (k = 1, ..., n) integrando numericamente l’equazione di Jeans,memorizzando quindi i valori in un array a n componenti. La dispersionedi velocita` alla posizione della i-esima particella ri viene quindi calcolata perinterpolazione.]3. Assumendo M = 1010M⊙ e b = 1 kpc, si tracci il grafico della posizione delleparticelle nei piani (x, y), (x, z), (y, z) e, nei due casi (metodo del rifiuto eapprossimazione maxwelliana) (x, vx), (y, vy) e (z, vz).4. Per il modello a N-corpi generato si calcolino:(a) la posizione del centro di massaxcm =1MN∑i=1mixi; (3.17)(b) i raggi r30, r50 e r70, tali che M(r30) = 0.3M , M(r50) = 0.5M eM(r70) = 0.7M ;(c) la velocita` del centro di massavcm =1MN∑i=1mivi; (3.18)Calcolo con Laboratorio II - a.a. 2007/2008 5(d) l’energia cineticaK =12N∑i=1mi[(vx,i − vx,cm)2 + (vy,i − vy,cm)2 + (vz,i − vz,cm)2]; (3.19)(e) la dispersione di velocita` viriale σV , definita daK =12Mσ2V ; (3.20)(f) l’energia potenzialeW = −N∑i=1∑j>iGmimj|xi − xj| ; (3.21)(g) il rapporto viriale |2K/W |;(h) per le proiezioni lungo gli assi x, y e z, il raggio efficace Re (taleche la massa proiettata entro Re sia 0.5M) e la dispersione di velocita`σp(R ≤ Re) proiettata entro un’apertura circolare di raggio Re: peresempio, proiettando lungo l’asse z,σ2p(R ≤ Re) =2MRi≤Re∑i=1,...,Nmi(vz,i − vz,cm)2, (3.22)dove Ri = x2i + y2i .

Appunti di lezione 3 + Esercizio C

Abstract

Similar works

Full text

Available Versions

Almae Matris Studiorum Campus

Almae Matris Studiorum Campus