2 Problemas

Talvez o maior problema enfrentado pelas simulações numéricas seja o recurso da suavização. Vimos que tal recurso nos é inevitável em simulações de galáxias, já que ainda não podemos realizar simulações com $10^{11}$ partículas. Por causa disso, as partículas em uma simulação não podem ser entendidas como a representação de estrelas individuais. Pelo contrário, devem ser consideradas como realizações estatísticas da distribuição de massa em um dado sistema. Enquanto pode parecer que a suavização tem apenas o efeito de diminuir a resolução espacial, e tem implicações apenas em pequena escala, alguns autores argumentam que a suavização pode provocar distúrbios em larga escala (ver, e.g., Gerber 1996). Isto porque a suavização modifica a energia potencial gravitacional de ligação do sistema, modificando a configuração de equilíbrio de sistemas auto-gravitantes.

Já foi mostrado que a suavização influencia a estabilidade local em discos, e Sellwood (1981) mostrou que o parâmetro de Toomre (1963, 1964) pode ser tão modificado por uma suavização imprópria que acaba suprimindo a instabilidade de barra.

Forças suavizadas implicam em que um sistema de partículas em equilíbrio não mais satisfaz o teorema do virial convencional, $2T+W=0$, onde $T$ é a energia cinética total e $W$ é a energia potencial total do sistema (ver, e.g., BT87). Entretanto, o virial de Clausius (ver Chandrasekhar 1957; Sellwood 1987) pode ser utilizado no lugar de $W$, recuperando o teorema do virial para qualquer lei de força. O virial de Clausius pode ser expresso por $\sum m\vec{x}\cdot\ddot{\vec{x}}$.

Merritt (1996) discute uma forma de se encontrar para um dado código acolisional a melhor escolha para o parâmetro de suavização. Se, por um lado, um valor excessivo para $\epsilon$ produz uma força muito distinta da Newtoniana, e acaba por alisar estruturas reais de pequena escala, por outro lado, um valor pequeno não satisfaz a necessidade de diminuir a granulosidade do sistema. Merritt chega a uma relação para um valor ótimo de $\epsilon$ em função de $N$, para uma esfera de Plummer:

$$\epsilon_{ot}\approx 1.1 \times N^{-0.28}.$$ (80)

E, para uma esfera de Hernquist (1990):

$$\epsilon_{ot}\approx 1.5 \times N^{-0.44}.$$ (81)

Athanassoula et al. (2000) discutem a otimização de $\epsilon$ para sistemas com distribuições de densidade distintas e chegam a resultados similares. Em particular, para sistemas triaxiais, o valor das razões axiais não interfere significativamente no valor de $\epsilon_{ot}$. No entanto, sistemas mais concentrados necessitam uma menor suavização. Isso significa que, na simulação de galáxias com diferentes componentes (e.g., bojo e disco), as forças na componente mais centralmente concentrada são representadas de forma mais pobre. O erro é maior quando a fração de massa na componente mais densa é maior.

Por outro lado, há outras formas de se estabelecer a suavização em uma simulação. Uma forma importante é a introduzida por Hernquist & Katz (1989). A suavização é feita através da seguinte função [considerando $\vec{F}=-m^2\vec{r}f(r)$ e $u=r/\epsilon$]:

$$f(r) = \left\{ \begin{array}{lll} 1/\epsilon{^3}[(4/3) - (6/... ...q u\leq 2$,} \\ 1/r^3, & \mbox{$u\geq 2$.} \end{array}\right.$$ (82)

Note que, nesta forma de suavização, a força passa a ser Newtoniana a partir de $r\geq 2\epsilon$. Ou seja, ao contrário da suavização padrão, esta é uma aproximação apenas para pequenas distâncias. Athanassoula et al. (2000) mostraram que esta técnica não produz resultados significativamente diferentes dos alcançados com a técnica padrão, apesar de serem mais precisos. Isso mostra que, embora existam, os problemas da suavização não devem ser exagerados.

Vale notar que alguns autores têm investido em códigos com parâmetros de suavização que variam no tempo e no espaço, embora haja controvérsia nos efeitos que esse esquema pode trazer. Gerber (1996) afirma que variações bruscas em $\epsilon$ podem mudar as condições de equilíbrio de um sistema.

Métodos distintos também sofrem problemas distintos. No caso do código hierárquico, um problema notável é a violação da terceira lei de Newton. Imagine uma partícula isolada A distante de uma partícula B, que faz parte de um agrupamento de partículas. A força de A em B é corretamente calculada, mas A somente sofre o efeito do agrupamento de partículas como um todo. Mesmo assim, a força total sobre A vinda do agrupamento é adequadamente calculada e, portanto, pode-se esperar que a evolução do sistema seja corretamente estabelecida. A precisão nos códigos hierárquicos aumenta quando se diminui o ângulo de abertura $\theta$, mas uma maneira mais eficiente de se melhorar os resultados é aumentando o número de partículas. No caso dos métodos com rede, especialmente em simulações tri-dimensionais, a geometria da rede pode introduzir problemas. Uma rede esférica, por exemplo, tem suas células com formato de cunha próximo aos pólos, e isso pode produzir sub-estruturas espúrias (ver, e.g., Klaric 1995). Outro problema é o fato de que partículas que escapam da rede deixam de entrar na contabilização das forças. Essa é uma poderosa vantagem dos códigos de árvore sobre os métodos com rede. Os primeiros são independentes da geometria, não perdem partículas, e não perdem tempo fazendo cálculos em células vazias.

Uma forma de verificar a exatidão dos resultados de simulações é verificar se a energia total do sistema e o seu centro de massa se conservam durante a execução do código. Diferenças inferiores a 1% são desejáveis. Uma forma mais robusta é verificar a reversibilidade do experimento. O experimento deve ser reversível, exceto por erros de arredondamento. Isto, no entanto, não é válido no esquema ``leapfrog'', que é intrinsicamente reversível.

Vejamos mais algumas considerações interessantes. Hernquist & Barnes (1990) fizeram uma avaliação a respeito da taxa de relaxação em diferentes métodos, tratando da evolução do mesmo sistema. Concluíram que os métodos comumente utilizados têm valores parecidos para estas taxas, o que significa que não há métodos mais ou menos colisionais que outros. Aarseth, Hénon & Wielen (1974) fizeram uma avaliação da evolução de uma esfera de Plummer auto-gravitante, também através de distintas técnicas. Verifica-se que o acordo entre diferentes métodos é muito bom, o que é encorajador.

Evidentemente, os sistemas estelares não são compostos apenas por estrelas. No caso de galáxias, por exemplo, o meio interestelar pode influenciar a evolução dinâmica do sistema como um todo. Isto é particularmente importante para galáxias de tipos morfológicos tardios, ricas em gás, e ainda mais para aglomerados de galáxias. Além disso, outros fenômenos ocorrem em um sistema estelar, além da gravitação. Em particular, a formação e a evolução de estrelas são detalhes importantes na evolução de um sistema, já que há, na primeira, consumo do meio interestelar e o surgimento de novas partículas, e, na segunda, injeção de matéria e energia para o meio interestelar. Em particular, uma técnica conhecida como SPH (``smoothed particle hydrodynamics'' - hidrodinâmica de partículas suavizadas) é utilizada para modelar o meio interestelar. A técnica também utiliza partículas, que, no entanto, são tratadas como um fluido que obedece às equações de um gás compressível. As propriedades termodinâmicas e hidrodinâmicas do fluido são atualizadas utilizando as equações de conservação da hidrodinâmica. Inclui-se viscosidade, aquecimento e esfriamento do gás, de forma artificial. A formação estelar é em geral mimetizada da seguinte forma. A cada passo de tempo, a taxa de formação estelar é calculada proporcionalmente à densidade local do gás. A quantidade de gás que ultrapassa um certo limite na densidade é aos poucos transformada em estrelas. A transformação vai retirando partículas que se comportam segundo as leis do SPH e adicionando partículas à parte estritamente gravitacional do código. Para mimetizar a injeção de energia por supernovas em regiões de formação estelar, pode-se adicionar energia cinética ao gás que se encontra na região onde há a formação de estrelas. Códigos mais completos são, evidentemente, mais complexos, e, muitas vezes, mais caros em termos de tempo de cálculo, ou possuem certas fraquezas ao desconsiderar outros aspectos importantes. Por outro lado, a mimetização de certos efeitos deve ser tradada com cuidado, para evitar efeitos espúrios. Certamente que devemos aperfeiçoar os códigos a fim de torná-los sempre mais realistas. Friaça e Terlevich (1998), por exemplo, utilizam um código que inclui evolução química e hidrodinâmica (ainda que seja apenas uni-dimensional) para estudar a formação e evolução de galáxias elípticas e sua relação com AGN's. Enfim, é evidente que cada problema pode ser considerado por códigos específicos, posto que há ainda muito o que desvendar. Assim, apesar de incompletos, os códigos $N$-corpos são uma poderosa ferramenta no estudo de sistemas estelares.

É interessante notar que, apesar de os tratamentos analíticos serem exatos, estes tratam dos sistemas estelares como se fossem sistemas contínuos, o que certamente não são. Por outro lado, como vimos, as simulações impõem uma granulosidade excessiva aos sistemas, muito maior do que nos sistemas reais. Esta diferença essencial produz resultados distintos, como é o caso da dispersão de velocidades necessária para estabilizar um disco estelar. O tratamento analítico prevê valores menores do que os necessários em simulações (Miller 1972). Esse resultado nos indica que nos aproximamos da realidade ao considerarmos ambos os tratamentos (analítico e numérico) e percebendo que a realidade deve parecer-se com o caminho do meio!