Struttura a grande scala e struttura
delle galassie
Giuseppe Murante
Gli elementi del problema

Teoria
– Un po’ di Cosmologia
– Modello gerarchico: cenni
– Tecniche N-Body: perche’, come, quando

Applicazioni: esperimenti N-Body
– Conseguenze del modello gerarchico sulla
struttura dei dischi galattici
– Risultati preliminari dello studio
dell’evoluzione di dischi galattici esponenziali
in ambiente cosmologico
Struttura a grande scala e struttura
delle galassie: teoria
Giuseppe Murante
23 aprile 2002
Cenni di Cosmologia
La dinamica e' nell'evoluzione del fattore di scala a(t)
Descrizione generale dell'Universo: eq. di Friedmann-
Lemaitre-Robertson-Walker, soluzioni non stazionarie
dell'eq. di campo di Einstein
Parametri:
3H 2 (t )
8G
H (t )  a (t ) / a (t ); H 0  100h km s 1 Mpc 1
 c (t ) 
m (t )   ( t ) /  c (t )
densita' critica
Cost. di Hubble
Parametro di densita'
  m     k
>0: espansione accelerata
k0: geometria non piatta
Densita‘ totale
L'evoluzione del modello (e
delle strutture) dipende anche
dal tipo di materia oscura.
La materia oscura

Se si misura la massa luminosa, si arriva
ad un centesimo della densità critica.
 …ma misure della massa necessaria per
tenere legate gravitazionalmente le
strutture più grandi che vediamo
forniscono valori più alti…
 Anche la stabilità del disco delle galassie,
e la forma delle curve di rotazione,
richiedono più massa...
 Inoltre: misure di lensing gravitazionale,
luminosita` X degli ammassi di galassie..
Materia non barionica

Questa massa mancante, la materia oscura,
è difficile che sia “materia normale”
(barionica) sia per motivi astrofisici (si
dovrebbe poterne dedurre l’esistenza da
osservazioni) che cosmologici (limiti dalla
teoria della nuclesintesi)
 Molte ipotesi su cosa sia e quanta ne
esista; nessuna osservazione sin’ora.
 Materia oscura fredda: priva di
dispersione termica di velocita’ (se fosse
calda –neutrini massivi- problemi a
formare le stutture ad alto redshift).
La formazione delle strutture

I cataloghi di
 Esempio: LCRS:
galassie e di
clusters di
galassie mostrano
strutture: vuoti,
filamenti e muri di
galassie…
 Come si originano
queste strutture? E
le galassie?
Modello di formazione


Teoria dominante (non unica): le piccole
disomogeneità iniziali, viste da COBE, collassano
gravitazionalmente generando le strutture osservate.
Le fluttuazioni (gaussiane) di densita’ sono descritte
dal loro spettro di potenza P(k).
 2 (r) 


1
2 2
ikr 2
2
n
P
(
k
)
e
k
dk
;
P
(
k
)

P
(
k
)
T
(
k
);
P
(
k
)

Ak
in
in

Il tipo preciso di strutture (numero, posizione,
dimensioni…) dipende dal tipo di materia oscura,
dalla quantità di materia totale, ed altri parametri
cosmologici, oltre a quelli globali.
I piu’ importanti parametri sono la funzione di
trasferimento T(k) e l’indice primordiale n, ma i
dettagli del collasso sono dati anche dal tipo di
materia oscura.
Evoluzione delle fluttuazioni

L’evoluzione gravitazionale della
materia non collisionale fredda e’
descritta dalle equazioni di VlasovPoisson:
2 ( x, t )  4a 2 (t ) ( x, t )
Poisson
f ( x, p, t )
p
f ( x, p, t )

 f ( x, p, t )  m ( x, t ) 
0
3
t
ma (t )
p
Vlasov
Funzione di distribuzione:
f(x,p,t) descrive la dinamica degli elementi di massa nello
spazio delle fasi. Inizialmente determinata da P(k).
Equazioni di Eulero-Poisson
b (t )
: densita' media dell'Universo al tempo t
: contrasto di densita'
 ( x, t )  [  ( x, t )  b (t )] / b (t )
( x, t ) : potenziale gravitazionale (peculiare)
f(x,p,t) : funzione di distribuzione degli elementi di massa
Se si richiede che f(x,p,t) sia funzione a singolo valore della
posizione (no multistream) si ottengono le eq. di Eulero-Poisson:
 ( x, t )
   v( x, t ) ( x, t )  0
t
 v ( x, t )
a (t )
 ( x, t )
2
v ( x, t )  (v ( x, t )   )v ( x, t )  
t
a (t )
t
 2  ( x, t )  4a 2 (t ) ( x, t )
Continuita'
Eulero
Poisson
Nota: entra il contrasto di densita' non la densita' globale: le equazioni
descrivono l'evoluzione delle strutture non quella dell'Universo
Studio della formazione

Osservativo: nuovi cataloghi, nuove
misure CMBR, osservazioni in cielo
profondo nell’infrarosso e nei raggi X,
lensing gravitazionale, effetto SZ…
 Teorico: calcolo analitico delle
caratteristiche dei “modelli cosmologici”
in regime lineare e semilineare; loro studio
con simulazioni numeriche al calcolatore.
Struttura degli aloni di materia oscura

Collasso sferico di una perturbazione
omogenea: possibili calcoli analitici a
partire dalle equazioni di Eulero-Poisson.
Leggi di scala, raggio e densita’ viriali.
 Struttura degli aloni
– Richiesto calcolo numerico
– Profili universali di densita’
– Momento angolare, triassialita’, velocita’
Numero di aloni: funzione di massa

Press&Schecter 1974: a partire dal collasso sferico
uniforme, calcolano il numero di aloni di massa M dato uno
spettro di fluttuazioni di densita’ P(k) (distribuite in modo
Gaussiano):
dn
2 b   c  d

e


dM
 M   ( M )  dM
 c2
2
4R3 b
; c 
;M
 (M )
3
Dove le fluttuazioni di densita’ sono date da:
c
 2 ( M )  4  WT2 (kR) P(k )k 2dk
• Formula ri-ricavata con approcci statistici tipo random-walk
(excursion sets: top-hat nello spazio di Fourier invece che collasso
sferico). Fornisce anche la distribuzione dei progenitori degli oggetti
appartenenti ad una data classe di massa.
•Buon accordo con N-Body; accordo ancora migliore con formule
modificate (es.: Sheth-Tormen 2000). Ma: l’identificazione oggetti nelle
simulazioni e’ un problema!
Numero di aloni: funzione di massa

W e’ un filtro nello spazio di Fourier che seleziona la scala R, in
questo caso un top-hat sferico:
 sin x

WT ( x )  3
 cos x  / x 2
 x


c e’ il contrasto di densita’ lineare di una fluttuazione di densita’
sferica collassata, c = 1.686(z)c, c=0.055(0.018) per cosmologie
critiche (senza costante cosmologica)
Nel caso di spettro di potenza a legge di potenza del tipo:
 k  k n M *  M 0 (1  z ) 6 /(3n ) 
2
Si ricava per la funzione di massa l’espressione:
2  b  n  M 

n( M , )dM 
1  
 M  3  M * ( ) 
3 n
 3
3 n
6
e
1 M

 
2  M * ( ) 
Modello gerarchico


Si dice “modello gerarchico” il modello in cui gli aloni si
formano prima a scale piccole e poi, per aggregazione
gravitazionale, a scale via via crescenti.
In contesto cosmologico, sotto l’ipotesi che la materia
oscura sia fredda, la formazione di strutture e’
gerarchica a qualunque scala alla quale lo spettro
di potenza abbia indice spettrale locale n-3 (es.:
Peebles 1980, “Large scale structure of the Universe)

Non e’ possibile confutare il modello gerarchico
senza confutare anche il modello di formazione
della struttura a grande scala dell’Universo, e di
conseguenza l’intero modello standard della
cosmologia!
Tecniche N-Body
Codici diretti
 Codici PM
 Codici misti
 Tree-codes

Calcolo parallelo
A causa dell’alta risoluzione richiesta, e’ necessario fare uso di
supercalcolo (grandi computers di solito paralleli)
Simulazione dell'interazione
gravitazionale
•Metodo
piu' semplice: un insieme di punti massivi
interagisce gravitazionalmente: forza di Newton:
N
Fi  
j 1
•Tutte
Gmi m j (r i  r j )
ri  r j
3
le particelle interagiscono con tutte le altre (codice
diretto)
•Discretizzazione del tempo: prima si calcolano tutte le
forze, poi si spostano tutte le particelle di un time-step.
•Condizioni al contorno naturalmente vuote.
Codice diretto: problemi
•Occorre
un controllo accurato del time-step per evitare
errori di integrazione.
•Occorre un algoritmo di integrazione raffinato per
ottenere accuratezza e controllo dell'errore numerico.
•Tanto piu' le particelle sono vicine, tanto piu' piccolo deve
essere il time-step; se e' troppo grande le particelle
possono schizzare via su traiettorie iperboliche:
collisionalita'.
•Il tempo di calcolo necessario e' proporzionale al
quadrato del numero di particelle:
TN
2
Softening
•Per
aggirare il problema della collisionalita' viene
introdotto un parametro di softening che modifica la forma
della forza, ad esempio:
N
Fi  
j 1
Gmi m j (r i  r j )
( ri  r j   )
•A scale
3
dell'ordine del softening, la forza non e' piu'
newtoniana, e neanche conservativa (!)
•Tests: i risultati della simulazione sono uguali a quelli di
un calcolo senza softening su scale
l  K , K  1.5
•
Vari tipi di softening (K varia!)
Regolarizzazione; time-step
variabile
•Non
sempre si puo' trascurare quello che
avviene a piccola scala: es.: globular clusters.
•Introdotta la regolarizzazione: trattazione
semianalitica degli incontri a 2 o 3 corpi.
•Introdotte tecniche di integrazione alle differenze
finite (Ahmad-Cohen 1973), time-step variabili
per ciascuna particella e metodi predictorcorrector per poterli utilizzare.
•Schema "standard" per un codice diretto o PP
(Particle-Particle)
Schema codice PP
t  min i (ti  ti )
Calcolo nuove posizioni
valutazione nuovo t
i
Forza
estrapolata(predictor)
o calcolata (corrector)
Calcolo coordinate (tutte le
particelle o solo i vicini)
Calcolo forza
(predictor)
sulla particella i
t i
No
piccolo?
Si
Aggiorno posizione
e tempo particella i
Ricalcolo vicini
ed intervallo temporale
Ricalcolo forza
(corrector)
N-Body in Cosmologia
•Dark matter --> fluido non collisionale freddo.
•Le "particelle" N-Body rappresentano un elemento
di fluido :
ne serve il maggior numero possibile
•Le particelle devono campionare il fluido:
assegnazione della massa
•Necessario il softening (il fluido non ha incontri a due corpi)
o una tecnica equivalente
•Coordinate comoventi: si segue l'evoluzione delle
perturbazioni di densita' (Vlasov-Poisson).
•I coefficienti delle equazioni di evoluzione gravitazionale
dipendono ora dal tempo attraverso il fattore di scala a(t)
•a(t) contiene il modello (e tutta la Relativita' Generale).
•Integrazione diretta del campo, nello spazio delle fasi,
numericamente troppo pesante!
Codici Particle-Mesh
•Utilizzo
di tecniche su griglia ed FFT (eq. di Poisson risolta
nello spazio di Fourier, integrazione nello spazio fisico). NO
collisionalita' "diretta".
•Condizioni al contorno naturalmente periodiche.
•Il potenziale e le forze vengono calcolate su griglia ed
interpolate alle particelle; il campo di densita' viene
ricalcolato a partire dalle particelle.
•Integratore: leap-frog secondo ordine con opportuni
coefficienti che dipendono dal tempo.
•I coefficienti dipendenti dal tempo contengono tutta la
cosmologia (modello e dinamica).
•Il tempo di calcolo e' dominato dal tempo impiegato per le
FFT:
T  N log 2 N
Assegnazione massa ed interpolazione
•Schemi di assegnazione della massa:
•NGP (Nearest Grid Point): la massa assegnata
alla cella
in cui si trova la particella. Campo discontinuo.
•CIC (Cloud In Cell): massa assegnata alle celle entro cui
una "nuvola" centrata sulla particella, di dimensione R, si
trova, in modo proporzionale alla distanza del centro della
cella (di lato L) dal centro della nuvola (la particella).
Campo continuo, derivata prima discontinua.
Normalmente R=L, 8 celle.
•TCS (Triangular Shaped Cloud): l'assegnazione non e'
lineare ma pesata in modo che sia continuo il campo e la
sua derivata prima. R=1.5L (27 celle).
•Forza interpolata con schema inverso rispetto
all'assegnazione, per evitare self-forces.
Funzioni di Green
•Per
risolvere l'eq. di Poisson nello spazio di Fourier
occorre calcolare la funzione di Green dell'operatore
Laplaciano in tale spazio
•"Simple-man green function":
C
G 2
k
e' riferita all'operatore su R3; bisogna
considerare che la derivazione avviene in maniera
discreta e tener conto dell'assegnazione delle quantita':
la piu' semplice funzione di Green risulta:
•Questa
G

1
L2 sin 2 2k x / L  sin 2 2k y / L  sin 2 2k z / L
(Hockney&Eastwood 1980, Efsthatiou, Davis, Frenk, White 1985)
Schema codice Particle-Mesh
Campionamento a particelle del campo
Calcolo del potenziale gravitazionale
(FFT)
Calcolo delle forze su griglia
Interpolazione forze alle particelle
Spostamento delle particelle
Il campo viene ricalcolato su griglia
Problemi codici PM
•La
risoluzione L della simulazione e' quella della griglia; i
risultati sono affidabili a scale ~3L
•"Collisionalita' di second'ordine": dalla discretezza della
griglia (dipende dallo schema di assegnazione della
massa)
•Errori di integrazione di ordine O(L^2)
•Non c'e' perfetta conservazione dell'energia (eq. di
Layzier-Irvine)
•Codici di questo tipo NON SONO ADATTI, senza
notevoli modifiche, a problemi che non abbiano
condizioni al contorno periodiche (es.: studio sistema
solare; simulazione di UN cluster di galassie)
•Problema delle repliche: no oggetti non-lineari su scale
maggiori di 1/10 del box-size
Codici derivati
•P3M: PM a grande scala, PP a piccola scala.
-VANTAGGI: risoluzione dell'ordine di un PP,
ma
molto piu' veloce.
-SVANTAGGI: softening  collisionalita'; difficolta' di
implementazione; interfaccia PM/PP; risoluzione in
massa << risoluzione in forza
•APM: PM adattativo, itera la procedura di costruzione
della griglia nelle sone sovradense
-VANTAGGI: risoluzione in forza dell'ordine di un PP
con la velocita' di un PM
-SVANTAGGI: difficolta' di implementazione;
tipicamente, ancora, risoluzione in massa <<
risoluzione in forza
•AP3M: misto delle tecniche di cui sopra
Codici piu' moderni
•Treecode:
~ sviluppo in multipoli, dati ordinati ad albero.
Consente tecniche multimassa.
-VANTAGGI: risoluzione dell'ordine di un PP, ma
molto piu' veloce; simulazione oggetti isolati.
-SVANTAGGI: softening  collisionalita'; difficolta' di
implementazione; no periodicita‘ (ma: Ewald
summation) ; minore controllo degli errori di
integrazione
•"ART": PM adattativo con tecniche ad albero, senza
ripetizione griglia, multimassa
-VANTAGGI: risoluzione dell'ordine di un PP; puo'
simulare un oggetto ad alta risoluzione in ambiente
periodico a bassa risoluzione
-SVANTAGGI: difficolta' di implementazione; meno
veloce di un [A]P3M; no SPH (ma: ora, tentativi di
implementazione, Kratsov 2002)
Treecode
Lo spazio viene suddiviso in celle, a loro volta
suddivise ancora, sinche' ogni cella finale contiene
una sola particella
•Ogni
particella interagisce con il baricentro delle celle di
primo livello, a meno che la loro dimensione sia vista
sotto un angolo maggiore del "parametro di tolleranza". In
questo caso interagisce coi baricentri delle celle di
secondo livello, e cosi' via.
•Colle particelle vicine c'e' iterazione diretta (si risolvono
le celle sino al livello delle particelle).
•Implementazione parallela pubblica: GADGET (Springel
et al. 2001)
Adaptive Refinement Tree
•Codice PM al livello "zero"
•Se in una cella si supera una
soglia di densita', questa
viene divisa in 8, indi:
l'eq. di Poisson viene trasformata in eq. di diffusione:

  
  2  

2
una soluzione iniziale "rilassa" alla soluzione di equilibrio
- soluzione dell'eq. di Poisson - per
 
•La
soluzione iniziale proviene dal livello di refinement
precedente - al primo livello, proviene dalla soluzione su
griglia.
Struttura a grande scala e struttura
delle galassie:
esperimenti N-Body
Giuseppe Murante
Struttura degli aloni: profili di densita’

Navarro, Frenk, White 1996,1997: simulazioni
numeriche: gli aloni di materia oscura mostrano,
in regime fortemente non-lineare, profili radiali di
densita’ indipendenti dal modello cosmologico e
dalla massa.
 Tali profili hanno forma:
 ( r )   crit
c
( r / rs )(1  r / rs ) 2
• Quello che varia con la massa e’ la concentrazione c:
200
c3
c  r200 / rs   c 
3 [ln( 1  c)  c /(1  c)]
R
: raggio viriale
Struttura degli aloni: profili di densita’
Gli aloni di DM mostrano quindi
cuspidi in densita’ a tutte le scale di
massa.
 Discussione in corso sulla forma
esatta della cuspide: Moore et al
1999, Ghigna et al 2000:

1
 (r ) 
3/ 2
3/ 2
(r / rs ) [1  (r / rs ) ]
Struttura degli aloni: altre caratteristiche

Momento angolare, caratterizzato da:

E
1/ 2
J
GM 5 / 2
 distribuito in modo lognormale, e
< >  0.05
 Triassialita’: gli aloni sono prolati
 Profili (…universali?) di velocita’

Galassie in ambiente cosmologico



Trattazione N-Body “classica” per esperimenti di
stabilita’ dei dischi: si usano aloni “artificiali”
(profili isotermi, King modificati…) in cui
immergere il disco per studiare stabilita’, barre..
Trattazione “cosmologica”: si usano profili NFW.
Effetti gravitazionali importanti:
– Frizione dinamica
– Risonanze
– Violent relaxation

Studi su:
– effetti di incontri con “satelliti” di DM e/o stellari
– streams di DM e loro eventuale individuazione
– trasferimento di momento angolare dai barioni alla DM

Tutto questo non concerne il problema della
formazione delle galassie ma la loro evoluzione.
Fragilita’ dei dischi galattici



Velasquez, White 1999: satellite stellare accresce
sul disco. Orbite prograde: riscaldamento del
disco; retrograde: tilt coerente
Navarro 2001: disk/bulge/halo + diverse centinaia
di satelliti DM con masse, densita’ ed orbite
derivate da simulazioni cosmologiche; nessuna
evidenza di riscaldamento in eccesso rispetto a
quello numerico del caso isolato
Mayer, Governato, Colpi, Quinn, Moore… 19982001: evoluzione di galassie nane satelliti nel
gruppo locale
Streams di materia oscura





Helmi, White, Springel 2002: studio dello spazio delle fasi di
un alone DM di massa ~ quello della via lattea.
Altissima risoluzione numerica (67 milioni di particelle), ma:
no stelle-gas; a queste scale dovrebbero influire sulla
dinamica.
La simulazione era originariamente quella di un ammasso
di galassie, riscalato grazie all’autosimilarita’ dello spettro
di potenza.
Risultato: non si dovrebbero vedere gli streams, la
distribuzione nello spazio delle fasi e’ uniforme, tranne che
per le particelle piu’ energetiche che hanno “memoria” della
storia di merging dell’alone
Questa e’ la maggior risoluzione oggi ottenibile; non e’
ancora tecnicamente possibile una simulazione del genere
con stelle e gas.
Trasferimento di momento angolare





Il problema: teoria e simulazioni numeriche
prevedono distribuzioni con forti cuspidi nei centri
dei dischi galattici; le osservazioni trovano cores
(“cusps-cores problem”).
Proposto meccanismo di trasferimento di momento
angolare dal disco all’alone tramite risonanza della
DM con la barra (Weiberg, Katz 2001)
..ma: altri autori affermano che l’efficienza e’ bassa
(Valenzuela, Klypin 2002; Athanassoula 2002; CMM)
Pattern speed: Valenzuela, Klypin: quasi nessun
rallentamento; CMM: piccolo rallentamento,
Athanassolua: sensibile rallentamento.
..durata barra in ambiente cosmologico? (El-Zant,
Shloshman 2002; CMM)
Un esempio: disco stellare immerso
in alone cosmologico
Novita’ rispetto ai casi precedenti: l’alone evolve
autoconsistentemente nella simulazione cosmologica
 Tecnica di immersione:
– Simulazione cosmologica, 25 Mpc di lato, in cui si identifica un
alone di DM di dimensioni galattiche
– Disco esponenziale, generato in equilibrio col potenziale
gravitazionale dell’alone prescelto ed immerso nella simulazione
– Il momento angolare del disco e dell’alone viene allineato.
– Parametri da studiare: MDM vs Mstars, redshift di immersione e
conseguente diverso stadio evolutivo dell’alone. Attualmente:
redshift immersione z=2
– Problema della risposta dell’alone (non-fisica) all’immersione
– Studio casi “isolati” per determinarne l’impatto


Simulazioni multimassa ad alta risoluzione: solo la zona di interesse
e’ campionata finemente. Codice: GADGET; calcolatori utilizzati:
Beowulf(OATo), t3e(cineca), sp3(cineca)
Curir, Mazzei,Murante 2002
Il disco galattico nell’alone cosmologico
MDM: Mstars= 1:1
MDM: Mstars= 3:1
Disco galattico a z=1.75. Si nota la maggiore
stabilita’ quando la DM abbia massa
percentualmente maggiore rispetto al disco.
Contorni di isodensita’ nel disco
(Caso 1:1)

Si nota la presenza di una barra stellare; la sua pattern speed
puo’ essere stimata.
Momento angolare e profilo di densita’ dell’alone
Trasferimento di momento angolare dal disco all’alone ed
irripidimento del profilo di densita’ dell’alone. Il primo e’
evidente negli andamenti opposti con il redshift mostrati in
figura; il secondo e’ appena apprezzabile. (Caso 1:1)
Conclusioni




Modello di formazione delle strutture su grande
scala “robusto”, inserito in ambito cosmologico,
successi osservativi.
Questo modello inizia ad essere applicato anche
allo studio della struttura delle galassie.
Ancora problemi di potenza di calcolo…
Diversa questione: galaxy formation,
autoconsistente, a partire dalla cosmologia
– Tecniche semianalitiche
– Tecniche numeriche con gas, SF, feedback
– …questo problema sara’ trattato nella
presentazione di A. Diaferio.
Profilo di densita’ di Navarro-Frenk-White
Profilo di Moore et al
Orbita prograda:
si nota
l’ispessimento
del disco
Orbita
retrograda: il
disco mostra
un tilt
(l’ispessimento
apparente qui e’
dovuto proprio
al tilt)
Riscaldamento del disco,
misurato con le
dispersioni di velocita’
in R,z. Non si apprezza
differenza tra modelli
con e senza
sottostruttura.
Svuotamento alone
Risonanza
orbitale


Weinberg, Katz 2001: forte trasferimento di momento
angolare dalla barra all’alone tramite risonanza.
E’ stato usato un codice simplettico ad alta precisione.
Funzione di massa integrale, CDM
(Brian, Norman 1998)

Valenzuela, Klypin 2002: formazione di barra in un disco
esponenziale immerso in un alone NFW. Poco trasferimento
di momento angolare dal disco all’alone.


El-Zant, Shloshman 2002: traiettorie particelle nella barra sono
caotiche
Il volume nello spazio delle fasi di queste traiettorie non coincide
con quello della barra: la barra si dissolve in un tempo t<1/H
In figura gli esponenti di Liapunov
per traiettorie inizialmente situate
nella barra.
Utilizzato un potenziale analitico
costante, tranne che per la barra
che, ruotando, ne genera uno
dipendente dal tempo
Traiettorie integrate direttamente:
NON e’ un N-body.

Helmi, White, Springel 2002.
In basso, la distribuzione in slices dello
spazio delle fasi delle particelle di DM
in un cubo di 2 kpc centrato a 8.5 kpc
dal centro dell’alone.
Le particelle evidenziate sono quelle ad
alta velocita’, provenienti da satelliti di
DM che hanno originato mergers
relativamente recenti con l’alone
principale
Condizioni iniziali per l’immersione del disco. I diversi colori mostrano
le differenti shell contenenti particelle via via piu’ massive. La zona
blu e’ ad alta risoluzione; quelle circostanti verdi e rosse
contengono due shells ciascuna.
Scarica

G. Murante