Lattice QCD per principianti, 12/N

Quisquilie e pinzellàcchere

Avevo promesso almeno un bonus track, ossia come all'atto pratico si portano avanti le simulazioni. L'affare è altamente non banale e spero di riuscire a rendere almeno l'idea. Ci sono altre questioni pratiche che ho lasciato in sospeso. Per esempio come si trattano i gradi di libertà fermionici; li si deve prendere con le molle, con estrema cautela, perché causano casini: non si può discretizzare in maniera ingenua l'equazione di Dirac. Varie ed eventuali. Vamos.

Hybrid Monte Carlo

Ricapitolo l'algoritmo di Metropolis per chi se lo fosse perso:

1. Prendi una qualsiasi variabile di campo del tuo sistema; nel nostro caso, una variabile di link, quindi una matrice $\text{SU}(3)$;
2. Modificala un po';
3. Calcola il cambiamento di azione $\Delta S$;
4. Accetta la nuova configurazione con probabilità $\text{min}\{1,\exp(-\Delta S)\}$, e in caso di mancata accettazione tieni la vecchia configurazione;
5. Torna al punto 1 un numero sufficiente di volte.

Ora, il problema è che la nostra azione di QCD è: $$ S = S_g - \sum_f\log\text{det}\,M_f $$ in cui $S_g$ è l'azione di pura gauge e $M_f$ è la matrice di accoppiamento fermionica per un quark di sapore $f$. Questo significa che se applicassimo in maniera pedissequa l'algoritmo di Metropolis, ogni volta che cambiamo un link dovremmo calcolare come si modifica il determinante fermionico, che ovviamente dipende dai campi di gauge. Questo comporta un numero di operazioni $O(N^3)$, dove $N$ è il numero di link del nostro reticolo. Poiché uno sweep Monte Carlo richiede $O(N)$ operazioni, ci troviamo alle prese con un algoritmo che richiede $O(N^4)$ operazioni per un singolo sweep. Decisamente no buono.

La prima parte della soluzione l'abbiamo adombrata nell'ottava puntata di questa serie, quando abbiamo scritto $$ \det M =  \int [{\cal D}\phi^\dagger][{\cal D}\phi] \text{e}^{-\phi^\dagger M^{-1} \phi} $$ in cui $\phi$ è un campo bosonico (complesso) un po' speciale, perché oltre a $3$ indici di colore possiede anche $4$ indici aggiuntivi di spin che si contraggono con gli indici di Dirac presenti in $M^{-1}$. Insomma è quel che tecnicamente si chiama uno pseudofermione. Per semplicità di notazione d'ora in poi considererò un solo indice di flavour.

Voi direte: "Bel guadagno! Prima all'esponenziale avevamo un'onesta matrice estremamente sparsa, e l'abbiamo scambiata con la sua inversa, che al contrario introduce un'interazione estremamente non-locale tra gli pseudofermioni!" Sebbene questo sia vero, il vantaggio diventa chiaro quando introduciamo i momenti coniugati alle variabili di link, $\pi_\mu(x)$, e un termine cinetico aggiuntivo nel peso di Boltzmann. Possiamo quindi definire una "vera hamiltoniana" $$ H = \frac{1}{2}\pi^2 + S_g -\phi^\dagger M^{-1} \phi $$ e si può dimostrare che il valore d'aspettazione nell'ensemble microcanonico definito da questa hamiltoniana coincide con il valore d'aspettazione canonico definito prima. Possiamo quindi utilizzare le equazioni del moto di Hamilton e scrivere un tranquillo programma di Dinamica Molecolare. I nostri "vecchi" valori d'aspettazione coincidono ora con la media temporale su un tempo fittizio $\tau$ d'evoluzione del sistema. Notate che non occorre introdurre i momenti coniugati di $\phi$ in quanto $\phi$ già entra nell'hamiltoniana con un termine "cinetico" (gaussiano). L'algoritmo diventa dunque:

1. Scegli i campi $\pi_\mu(x)$ con distribuzione gaussiana;
2. Scegli un campo $\eta = M^{-1}\phi$ (quindi $\phi = M\eta$) con distribuzione gaussiana;
3. Integra le equazioni del moto di hamilton, numericamente, per un numero di passi pari a $\tau/\delta\tau$, in cui $\tau \simeq 1$ e $\delta\tau$ passo d'integrazione; per ogni passo d'integrazione occorre calcolare $M^{-1}\phi$, che è molto molto più semplice rispetto al calcolo dell'intero determinante;
4. Calcola i tuoi valori d'aspettazione e torna al punto 1 tante volte quanto necessario.

In realtà se usassimo questo algoritmo, che è pura Dinamica Molecolare, introdurremmo un errore dovuto all'integrazione numerica, perché per quanto riduciamo $\delta\tau$ non riusciremo mai a conservare esattamente l'hamiltoniana del sistema. Quindi tra il passo 3 e il passo 4 inseriamo un nuovo step: 

• Calcola $\delta H$ e accetta la nuova configurazione con probabilità uguale a $\text{min}\{1,\exp(-\delta H)\}$; se la configurazione non è accettata tieni la vecchia.

Questo algoritmo è noto con il nome di Hybrid Monte Carlo.

Fermion doubling

In una delle puntate precedenti abbiamo visto che per motivi di hermiticità la migliore definizione di derivata discreta, per variabili fermioniche, è $$ \partial_\mu\psi(x) = \frac{\psi(x+\mu)-\psi(x-\mu)}{2a} $$ Se ora prendiamo un fermione massless, libero, cioè poniamo tutti i link $U_\mu(x) = 1$, e calcoliamo l'inverso del propagatore fermionico nello spazio dei momenti, tramite una semplice trasformata di Fourier otteniamo: $$ S^{-1}(p) = \frac{i}{a}\gamma_\mu\sin(ap_\mu) $$ (ricordate la somma sugli indici ripetuti). Se mandiamo $a \to 0$ facendo il limite otteniamo l'usuale risultato nel continuo: $$ S^{-1}(p) = i\gamma_\mu p_\mu $$ (ricordate che siamo nell'euclideo e non c'è bisogno di alzare e abbassare gli indici). Purtroppo c'è un problema: poiché siamo su reticolo (consideriamolo pure infinito in tutte e quattro le direzioni, tanto è una discussione di principio), il momento $p_\mu$ ha dei limiti precisi: $ -\pi \le ap_\mu \lt \pi $. Questa è chiamata zona (o cella) di Brillouin. Vedete subito che c'è un problema, perché il propagatore ha un polo non solo per $ap_\mu = (0,0,0,0)$ ma anche in tutti gli angoli della zona di Brillouin, ossia in tutti i punti in cui sostituite uno o più $\pi$ al posto di uno o più $0$ nel quadrivettore $ap_\mu$. Quindi, riassumo: avete un fermione massless nel continuo, e lo riconoscete dal fatto che il suo propagatore ha un polo a momento nullo. Lo mettete su reticolo e ottenete $16$ fermioni, uno per ogni angolo della zona di Brillouin (in $d$ dimensioni ottenete $2^d$ fermioni). Cercate di tornare al continuo e questi $16$ fermioni vi rimangono sul groppone. È un bel problema, soprattutto se considerate che se ora aggiungete un termine di massa, questi $16$ fermioni avranno tutti la stessa massa, e di certo non vogliamo $16$ quark up, $16$ quark down eccetera...

Potreste pensare che tutto ciò discende dal modo in cui abbiamo scelto di discretizzare la derivata, e che se fossimo stati più furbi avremmo saputo evitare questo problema. In realtà le cose sono un po' più complicate di così. Esiste un teorema no-go, dovuto a Nielsen e Ninomiya, che dice questo: non è possibile scrivere un'azione fermionica discreta, su reticolo, che allo stesso tempo:

1. sia locale;
2. sia chiralmente invariante;
3. non abbia il problema del fermion doubling.

È una situazione obiettivamente imbarazzante. Se rinunciamo alla località rinunciamo anche alla possibilità concreta di fare simulazioni. D'altro canto la simmetria chirale è quella che garantisce che la rinormalizzazione della massa dei quark sia solo moltiplicativa: in altre parole, se ponete $m_0 = 0$, ossia massa nuda nulla, nella lagrangiana, questo vi garantisce che la massa rimane nulla anche sotto rinormalizzazione. D'altronde non possiamo tenerci il fermion doubling. Per fortuna esistono varie soluzioni al problema, ognuna delle quali comporta complicazioni tecniche. 

La prima soluzione è quella di rinunciare alla simmetria chirale esplicita su reticolo: occorre inserire nell'azione fermionica un termine irrilevante, di dimensione $5$, che quindi non distrugge la rinormalizzabilità della teoria ma rompe esplicitamente la simmetria chirale, che viene recuperata solo nel limite del continuo. Questi sono i così detti fermioni di Wilson. Rinunciare alla chiralità significa che sotto rinormalizzazione la massa dei quark acquista un termine che diverge linearmente come $a^{-1}$, e occorre fare fine-tuning della massa bare per cancellarlo. In pratica la massa rinormalizzata dei quark non è nulla quando è nulla la massa nuda, ma per un certo valore critico della massa nuda. Per trovare questo valore critico si sfrutta la PCAC, ossia che la massa al quadrato del pione, per piccole masse dei quark, è proporzionale alla massa dei quark. Quindi si estrapola nella massa quadra del pione come nel grafico sotto.

Un'altra possibile soluzione è quella di accettare almeno parzialmente il doubling, ma di spalmare i gradi di libertà fermionici su un ipercubo elementare, in modo che su ogni sito compaia un solo indice di spin: in questo modo si riesce a ridurre il fattore di doubling da $16$ a $4$. Questo approccio va sotto il nome di fermioni staggered, ed è molto amato dagli americani perché comporta tempi di calcolo minori, ma ha parecchi punti deboli. Sebbene sia vero che la simmetria chirale sia parzialmente conservata, e che la rinormalizzazione della massa dei quark continui a essere moltiplicativa, come nel continuo, abbiamo pur sempre $4$ fermioni diversi con la stessa massa (ognuno dei $4$ tipi viene chiamato, in inglese, taste - confronta con flavour) e quel che otteniamo è un mescolamento spin-taste che rende complicato capire quali siano poi gli stati fisici della teoria.

Esistono soluzioni più arzigogolate, che comportano per esempio l'introduzione di una fittizia quinta dimensione (domain wall fermions), che sono forse le più pulite da un punto di vista teorico ma comportano tempi di calcolo enormemente dilatati.