La scoperta degli asteroidi near-Earth (NEA) al telescopio è un lavoro fondamentale per la mitigazione del rischio impatto degli asteroidi con la Terra. Tuttavia la semplice scoperta non basta a mitigare il rischio: è necessario conoscere anche l’evoluzione orbitale passata e futura del NEA per cercare di capire la sua origine e il suo futuro in relazione con la Terra. In questo articolo del blog vediamo come si possa calcolare la traiettoria di un NEA che si muove nel Sistema Solare, una volta note le condizioni iniziali di posizione e velocità. Come si vedrà il compito è tutt’altro che semplice e come esempio vedremo l’asteroide Apophis, che il 13 aprile 2029 farà un flyby molto stretto con la Terra.

Introduzione

Dopo la scoperta di un NEA la cosa più importante che si può fare nell’immediato è continuare a seguirlo per il maggior periodo di tempo possibile per determinarne con precisione l’orbita eliocentrica: solo conoscendo i 6 elementi orbitali ossia semiasse maggiore a, eccentricità e, inclinazione i, argomento del perielio omega, longitudine del nodo ascendente Omega e anomalia media M all’epoca T, è possibile propagare in avanti nel tempo la posizione dell’asteroide e calcolare la probabilità d’impatto dell’asteroide con la Terra.

Tuttavia c’è un problema, che spesso non viene messo in evidenza in campo divulgativo: le orbite kepleriane degli asteroidi, ad esempio quelle elencate nello Small-Body Database del JPL o nel database di NEODyS, sono valide solo per una certa epoca ossia per un dato istante di tempo. L’orbita di un corpo minore del Sistema Solare sarebbe perfettamente kepleriana ossia sempre uguale nel tempo solo se fosse l’unico oggetto a girare attorno al Sole, ma in realtà questo non accade: ci sono gli 8 pianeti – da Mercurio a Nettuno – più di un milione di asteroidi e alcune migliaia di comete e questo solo fra i corpi noti. Il risultato è che mentre un corpo minore orbita attorno al Sole perché è il corpo di gran lunga più massiccio di tutto il sistema, risente anche della forza gravitazionale esercitata dagli altri corpi (principalmente i pianeti) e questo ne altera la traiettoria nel tempo che non sarà più una ellissi kepleriana, ma una curva molto più complessa. Quindi gli elementi orbitali si riferiscono sempre alle orbite osculatrici kepleriane ossia a quelle traiettorie che verrebbero percorse se sul corpo, a partire dall’epoca dell’orbita, si esercitasse solo la forza di gravità del Sole, senza perturbatori. Ne deriva che anche le tre leggi empiriche di Kepler nel Sistema Solare sono solo delle buone approssimazioni della realtà e che valgono rigorosamente solo nel caso di due corpi.

In generale la forza di gravità esercitata sull’i-esimo corpo di massa m e con vettore di posizione r_i dagli N-1 che compongono il sistema sarà data dalla somma delle singole forze gravitazionali esercitate da tutti gli altri corpi. Usando la seconda legge della dinamica di Newton, F = ma, e l’espressione per la forza di gravità F = Gm₁m₂/r², si trova l’accelerazione che agisce sull’i-esimo corpo (i =1,…N):

Nell’Eq. precedente G = 6.672 10^-11 (N m²/kg) è la costante della gravitazione universale. Per trovare la funzione vettoriale r_i(t) che descrive la posizione dell’i-esimo corpo in funzione del tempo e la funzione vettoriale v_i(t) che ne descrive la velocità, bisogna quindi risolvere N equazioni differenziali tutte accoppiate fra di loro. A differenza delle più note equazioni algebriche, dove si deve determinare il valore di un’incognita (solitamente indicata con la lettera x), nel caso delle equazioni differenziali “l’incognita” è un’intera funzione come è, ad esempio, f(x) = x³+x+1, quindi si tratta di un problema che può essere enormemente più complesso e non sempre risolvibile analiticamente. Qui non è possibile dire di più, ma per chi volesse approfondire le nozioni di base del calcolo differenziale e integrale partendo da zero e senza liceo scientifico alle spalle, consiglio la lettura del libricino di Gustavo Bessiere: “Il calcolo differenziale e integrale reso facile e attraente” che, nonostante gli anni, resta sempre un ottimo punto di partenza per capire i concetti di base della matematica superiore.

Sfortunatamente il sistema di equazioni differenziali scritto sopra non è risolvibile analiticamente, quindi non esiste una funzione matematica in grado di descrivere il moto di un asteroide nel Sistema Solare, così come esiste nel caso i corpi siano solo 2, come ad esempio in un sistema Sole-pianeta. Per questo motivo, se si vuole prevedere la posizione di un asteroide nel tempo, è necessario ricorrere alle tecniche del calcolo numerico. Come si può procedere? Conoscendo posizione e velocità iniziali degli N-1 corpi e dell’asteroide al tempo T₀ si può calcolare numericamente posizione e velocità per il tempo T₁=T₀+dt, dove dt è un intervallo di tempo sufficientemente piccolo in modo tale da poter considerare le forze esercitate sull’asteroide come costanti. In questo modo si possono ottenere sia la nuova posizione r₁ = r₀ + v₀ dt, sia la nuova velocità v₁ = v₀ + A₀ dt e ripetere il ciclo per un dt successivo fino a ottenere la posizione e la velocità dell’asteroide per l’epoca finale desiderata. Naturalmente per ogni ciclo la posizione dei pianeti cambia perché si spostano e si sposta pure l’asteroide, quindi la forza gravitazionale che si esercita su quest’ultimo cambia a ogni step del calcolo. Arrivati all’epoca finale desiderata con i valori finali dei vettori r e v (per un totale di 6 numeri), si possono anche ottenere i 6 nuovi elementi orbitali a, e, i, omega, Omega e M che caratterizzano l’orbita osculatrice per l’epoca che – ovviamente – saranno diversi da quelli iniziali a seconda delle perturbazioni più o meno intense subite dall’asteroide fra l’epoca iniziale e quella finale. Fare questi calcoli manualmente è impensabile, per via dell’elevato numero di operazioni necessarie, ma è un compito perfetto per un computer in grado di eseguire un elevato numero di operazioni al secondo.

Fra i software liberamente disponibili in grado di effettuare un’integrazione numerica delle orbite dei corpi minori uno dei più user-friendly è senz’altro MERCURY6.2 di John E. Chambers. Si tratta di un integratore a N-corpi con interfaccia testuale che utilizza metodi come Bulirsh-Stoer, Everhart e altri. MERCURY6 è progettato per calcolare l’evoluzione orbitale di oggetti che si muovono nel campo gravitazionale di un corpo centrale massiccio, quindi si possono integrare le orbite dei corpi minori del Sistema Solare, oppure affrontare numericamente il problema dei 3 corpi oppure descrivere l’evoluzione di sistemi planetari extrasolari e così via. Il software include gli effetti delle forze gravitazionali di tipo newtoniano che si esercitano tra corpi che si presume siano masse puntiformi, anche se il corpo centrale può avere forma non sferica e, nel caso delle comete, può includere anche forze non gravitazionali. Il pacchetto MERCURY è composto da diversi driver e subroutine, insieme a un numero di file di input che vanno modificati per adattarli al sistema di cui si vuole simulare l’evoluzione, ma la configurazione è abbastanza semplice e il pacchetto corrente, scritto in Fortran, è configurato per utenti Linux e Mac OS X. Vediamo come installarlo su un sistema Linux e un esempio di utilizzo.

Installazione di Mercury6

Per prima cosa bisogna accertarsi di avere il compilatore gfortran installato nel proprio sistema. A questo scopo basta aprire un terminale e digitare which gfortran seguito da invio. Se non compare nulla va installato il compilatore digitando, sempre da terminale, sudo apt install gfortran. Il download dei pacchetti del compilatore Fortran e l’installazione nel sistema è questione di poche decine di secondi. A questo punto si può aprire il file d’archivio del pacchetto MERCURY scaricato da GitHub e mettere i file sorgente in una cartella del disco fisso. I file più importanti sono: mercury6_2.for (il programma per l’integrazione numerica delle orbite), element6.for (estrae in file separati gli elementi orbitali dei corpi di cui si è integrata l’orbita in funzione del tempo) e close6.for (per ciascun corpo dell’integrazione estrae le date e le distanza dei flyby con gli altri corpi). Sono questi i file sorgenti che vanno compilati da console usando il Makefile di MERCURY e digitando make buil. In questo modo vengono creati gli eseguibili mercury6, element6 e close6 che sono quelli che manderemo in esecuzione da console. Prima però dobbiamo sistemare i file di testo per l’input.

File di configurazione per l’input

I file di input sono importanti perché permettono di passare all’integratore tutte le informazioni necessarie per fare evolvere numericamente il proprio sistema di corpi e raccogliere i risultati. Vediamo i file più importanti, per una descrizione dettagliata si può consultare il file README.md incluso nel pacchetto di Mercury6.

big.in: questo file contiene i dati iniziali per tutti i corpi maggiori del sistema TRANNE che per il corpo centrale (il Sole se si sta facendo un’integrazione del Sistema Solare). Un corpo maggiore è definito come uno che perturba e interagisce con tutti gli altri oggetti durante l’integrazione e tipicamente si tratta dei pianeti. Tutte le righe che nel file iniziano con il carattere ) sono considerate commenti e verranno ignorate da mercury6, tuttavia, NON bisogna eliminare la prima riga di commento che inizia con )O+_06. I dati iniziali per l’integrazione numerica delle orbite possono essere di tre tipi:

1) Cartesiano: vanno specificate le tre componenti del vettore posizione iniziale e le tre componenti del vettore velocità iniziale rispetto al corpo centrale del sistema. La posizione va espressa in au (1 unità astronomica è pari a 149.597.870 km) mentre la velocità è espressa in au/d (1 unità astronomica al giorno è pari a 1731,457 km/s).

2) Asteroidale: in questo caso basta specificare gli elementi orbitali kepleriani osculatori: a = semiasse maggiore in au; e = eccentricità; i = inclinazione dell’orbita in gradi; omega = argomento del pericentro (perielio se il corpo centrale è il Sole) in gradi; Omega = longitudine nodo ascendente in gradi; M = anomalia media in gradi. L’anomalia media M, da cui si può ricavare l’anomalia vera v, fornisce la posizione del corpo sull’orbita.

3) Cometario: si tratta sempre di elementi kepleriani osculatori solo che al posto di a si fornisce q la distanza del pericentro (perielio se il corpo centrale è il Sole) e al posto dell’anomalia media si fornisce T, il tempo del passaggio al pericentro (perielio se il corpo centrale è il Sole).

Qualsiasi tipo di input si utilizzi, Cartesiano, Asteroidale o Cometario, va specificata anche l’epoca ossia il giorno giuliano per cui valgono gli elementi orbitali forniti. Nel file big.in fornito di default con Mercury6 si trovano i dati cartesiani per gli otto pianeti del Sistema Solare più il pianeta nano Plutone.

small.in: questo file contiene i dati iniziali dei corpi minori del sistema oggetto dell’integrazione. Un corpo minore è definito come uno che perturba e interagisce solo con i corpi maggiori durante l’integrazione. Quindi, i corpi minori si ignorano completamente: non si perturbano a vicenda e non possono scontrarsi tra loro. Se a questi oggetti non si assegna una massa, allora si comporteranno come particelle di prova senza perturbare i corpi maggiori. Se si vuole simulare il moto di asteroidi e comete nel Sistema Solare, considerate le masse trascurabili rispetto ai pianeti, è questa la scelta più efficiente da fare.

param.in: questo file contiene i parametri che controllano come viene eseguita l’integrazione. La prima opzione riguarda la scelta dell’algoritmo di integrazione, di default è selezionato il metodo Bulirsch-Stoer, che va bene un po’ per tutte le situazioni. La seconda opzione è il tempo di inizio dell’integrazione (in giorni). Questo tempo non deve necessariamente coincidere con l’epoca dei corpi maggiori, semplicemente l’integratore inizierà a produrre l’output a partire dalla data scelta dall’utente. Se si stanno integrando oggetti nel Sistema Solare, si può misurare il tempo in giorni giuliani (ad es. 1 gennaio 2000 = 2451544,5). Il terzo parametro è il tempo in cui terminerà l’integrazione (sempre in giorni). Se il tempo finale è inferiore a quello iniziale l’integrazione avviene indietro nel tempo, viceversa se il tempo finale è superiore a quello iniziale si va avanti nel tempo. Il parametro successivo è l’intervallo di output (in giorni), che determina la frequenza con cui il programma memorizzerà gli elementi orbitali per i corpi. I parametri successivi di questo file possono essere lasciati come sono. Infine, nei file elemen.in e close.in vanno elencati i nomi dei corpi minori di cui vogliamo gli elementi orbitali in funzione del tempo e le date dei flyby con i corpi maggiori. Questi file vengono letti dagli eseguibili element6 e close6, se non ci sono i nomi dei corpi minori desiderati non si produrrà alcun output. Per ripartire con una nuova integrazione occorre cancellare i vecchi file e si può fare rapidamente digitando make rm-gen nella console e premendo invio.

Un test su Apophis

Ora che abbiamo visto quali sono i principali file di configurazione di Mercury6 vediamo un esempio di funzionamento utilizzando i file di configurazione che sono dati di default. Come abbiamo già detto il file big.in contiene i pianeti del Sistema Solare più il pianeta nano Plutone. Nel file small.in ci sono due asteroidi near-Earth: Apollo e Apophis. Considerato che Apophis farà un flyby stretto con la Terra il 13 aprile 2029 verifichiamo questa previsione usando Mercury6. Gli elementi orbitali osculatori di Apophis li possiamo prendere dal JPL Small Body Database, oppure lasciare quelli di default. In param.in come data iniziale (start time) prendiamo il 25 febbraio 2023 (JD₀ = 2460000,5) e come data finale (stop time) il 17 ottobre 2047 (JD = 2469000,5). Come output interval scegliamo 10 giorni, così che gli elementi orbitali saranno salvati ogni dieci giorni a partire dal JD₀. Nei file elemen.in e close.in va messo il nome di Apophis nella stessa forma in cui si trova in small.in così da poter avere, terminata l’integrazione numerica del sistema, anche l’elenco degli elementi orbitali e dei flyby.

Ora, dal terminale aperto nella cartella dove abbiamo i file di configurazione e gli eseguibili di Mercury6 creati in precedenza, possiamo digitare ./mercury6 e premendo invio partirà l’integrazione della durata di qualche secondo. Con il procedere dell’integrazione verranno mostrate le variazioni dell’energia totale E e del momento angolare totale L. Visto che queste quantità sono conservare la loro variazione deve essere molto minore di zero ed è un modo per avere il check visivo che il calcolo numerico sta andando bene. Terminata l’integrazione numerica, digitando ./elements6 verrà creato il file contenente gli elementi orbitali di Apophis in funzione del tempo, infine digitando ./close6 verrà creato il file con l’elenco dei flyby. Alla fine dei calcoli la console dovrebbe avere l’aspetto di quella mostrata nella Fig. 2. Nella cartella di Mercury6 troveremo il file di testo APOPHIS.clo che mostrerà l’elenco dei flyby con i corpi maggiori (vedi Fig. 3) e il file di testo APOPHIS.aei che riporta l’elenco degli elementi orbitali in funzione del tempo. Il file con i dati numerici sugli elementi orbitali in funzione del tempo è meglio interpretabile se si usano per creare dei grafici: il risultato sarà simile a quello di Fig. 4.

Cosa succede se invece di andare avanti nel tempo, per vedere i futuri flyby di Apophis, si va un po’ indietro nel tempo fino a circa 54.000 anni fa? Quello che si trova è l’evoluzione orbitale di un tipico NEA soggetto alle perturbazioni gravitazionali dei pianeti Venere e Terra: a causa delle perturbazioni planetarie ci sono delle piccole oscillazioni del semiasse maggiore, delle oscillazioni abbastanza contenute sull’eccentricità e l’inclinazione dell’orbita (perturbazioni periodiche) e diverse rotazioni complete dell’argomento del perielio e della longitudine del nodo ascendente (perturbazioni secolari), che passano ripetutamente da 0° a 360° ricominciando il ciclo (vedi Fig. 5). Il valore dell’argomento del perielio, con il trascorrere del tempo, passa da 0° a 360° quindi la rotazione della linea degli apsidi (la retta che congiunge afelio e perielio), avviene in senso antiorario o diretto se vista dal polo nord dell’eclittica. La longitudine del nodo ascendente invece passa da 360° a 0°, quindi la linea dei nodi ruota in senso orario o retrogrado. La rotazione completa del perielio avviene in circa 10.000 anni, mentre quella del nodo ascendente è sugli 8000 anni. I minimi dell’oscillazione dell’eccentricità tendono a verificarsi per i valori di omega pari a 0° oppure a 180°. Il semiasse maggiore (e quindi il periodo orbitale), cambia poco perché è direttamente collegato all’energia totale dell’orbita di Apophis (somma di quella cinetica e di quella potenziale gravitazionale), quantità che tende a conservarsi. Le piccole oscillazioni del semiasse maggiore sono dovute al piccolo scambio di energia con i pianeti. Chiaramente più ci si allontana dall’epoca presente, più le incertezze degli elementi orbitali si amplificano determinando una non predicibilità della posizione dell’asteroide e della sua evoluzione orbitale, per questo motivo la Fig. 5 va intesa solo come un esempio della possibile evoluzione orbitale passata di Apophis.

Sappiamo che le orbite dei NEA sono caotiche ossia che con piccoli cambiamenti delle condizioni iniziali si ottengono delle differenze notevoli dello stato finale. Apophis non fa eccezione e si può vedere cambiando leggermente il valore del semiasse maggiore dell’epoca attuale da 0,92271626 au a 0,92278626 au e lasciando invariato tutto il resto. Questo equivale ad abbassarlo di appena 10.000 km ossia lo 0,007%. Il risultato è mostrato nella Fig. 6: come si vede l’andamento nel tempo degli elementi orbitali di Apophis ora è completamente diverso, anche se qualitativamente il comportamento è simile a prima con le perturbazioni periodiche e secolari. Questo è un esempio di orbita caotica, che non vuol dire che non sia calcolabile, ma che lo stato finale dipende in modo critico da quello iniziale.