kelvin@cluxnei:~/blog$ cat simulando-o-universo.md

2026-09-12 PT ~25 min de leitura read in English

Simulando o universo: a física por trás de 700 corpos numa aba do navegador

Este texto é uma exploração, não uma história. Vamos abrir o motor de uma simulação gravitacional e olhar cada peça: as equações e para que servem, os nomes reais das técnicas, os números medidos de desempenho, e a parte que quase ninguém escreve, que são os lugares onde a simulação mente de propósito, com o custo da mentira calculado.

A física, do começo

Antes de olhar o código, vale montar o vocabulário inteiro. São poucas fórmulas e cada uma responde a uma pergunta específica.

A única força que existe aqui

Tudo começa com a lei da gravitação universal de Newton. Dois corpos de massa m₁ e m₂, separados por uma distância r, se atraem com força:

Dois fatos fazem dessa fórmula simples um problema difícil. O primeiro é que a força cai com o quadrado da distância, mas não chega a zero nunca: todo corpo puxa todo corpo, para sempre. O segundo é que aceleração é força dividida por massa, então a aceleração que um corpo i sofre é a soma das contribuições de todos os outros:

Com dois corpos isso tem solução exata, que são as elipses de Kepler. Com três ou mais, não tem. Poincaré mostrou no fim do século XIX que não existe solução fechada geral para o problema de N corpos. É por isso que se simula: não porque a matemática seja desconhecida, mas porque a única forma de saber onde os corpos estarão daqui a um ano é avançar o tempo em passos pequenos e somar.

Densidade, raio e as quantidades derivadas

Um corpo não tem raio arbitrário. Dados massa e densidade, o raio é consequência de assumir uma esfera:

Amarrar o raio à densidade é a diferença entre um modelo e uma animação: a partir daqui, colisões, marés e densidades passam a significar alguma coisa.

A temperatura no centro de uma estrela sai de um argumento de equilíbrio: a pressão do gás tem que segurar o peso das camadas de cima. Disso vem uma estimativa conhecida como temperatura virial:

Aplicada ao Sol, ela dá 1,39 × 10⁷ K contra os 1,57 × 10⁷ K reais, ou seja 12% de erro com uma linha de conta. É boa o suficiente para decidir o que queima e o que não queima.

As distâncias que decidem o comportamento

Três raios aparecem o tempo todo, e cada um marca a fronteira de um regime diferente.

A esfera de Hill é a região onde a gravidade de um corpo ganha da estrela que ele orbita. É o território dele:

O raio de maré é onde a diferença de gravidade entre o lado próximo e o lado distante de um corpo supera a gravidade que o mantém inteiro. Passou disso, ele se despedaça:

E o raio de Schwarzschild é o horizonte de eventos de um buraco negro, a distância em que a velocidade de escape iguala a da luz:

Guarde esse último. A comparação entre r_s e as outras escalas da simulação vai decidir, mais adiante, o que vale a pena simular e o que é decoração.

Unidades: por que G = 4π²

No Sistema Internacional, a constante gravitacional é 6,674 × 10⁻¹¹ m³·kg⁻¹·s⁻². É um número horrível para trabalhar: as massas ficam na ordem de 10³⁰, as distâncias em 10¹¹, e um float64 gasta boa parte da precisão só carregando expoente.

A saída clássica em astronomia é trocar as unidades. Se você medir massa em massas solares, distância em unidades astronômicas e tempo em anos, algo bonito acontece.

Considere uma órbita circular. A aceleração centrípeta tem que ser exatamente a gravidade:

Agora imponha a definição da própria unidade astronômica: um corpo a 1 UA de uma estrela de 1 M☉ completa a volta em 1 ano. A velocidade dessa órbita é o perímetro dividido pelo período, ou seja 2π. Substituindo r = 1 e M = 1:

E isso entrega de graça o melhor teste que essa simulação tem. Está escrito como doutrina no cabeçalho de units.js:

Coloque um corpo a 1 UA com velocidade 2π UA/ano e ele tem que fechar a órbita em exatamente um ano.

É um teste falsificável que não depende de nenhuma outra parte do código estar certa. Se ele falhar, a base está errada e não adianta olhar mais nada.

A singularidade e o softening de Plummer

A lei de Newton tem um defeito fatal para quem simula: quando r → 0, a força vai para o infinito. Dois corpos que passam muito perto recebem uma aceleração absurda em um único passo de tempo e saem disparados, e a simulação explode. Não é um bug de implementação. É a fórmula sendo honesta sobre uma situação que o passo de tempo não consegue resolver.

A solução ingênua é limitar a aceleração num teto arbitrário. Funciona no sentido de que nada voa para o infinito, mas é fisicamente indefensável: toda vez que o limitador dispara, você está injetando ou removendo energia do sistema sem que nada físico tenha acontecido.

A solução usada aqui é o softening de Plummer, que adiciona um comprimento ε ao denominador:

Não há descontinuidade, não há energia criada do nada, e ainda economiza uma normalização de vetor que poderia gerar NaN quando a distância é exatamente zero.

O critério para escolher ε é simples de enunciar: ele precisa ser muito menor que as separações orbitais que interessam, ou distorceria as próprias órbitas que queremos medir. Um quinto de raio solar contra separações de ordem 1 UA passa folgado. O efeito medido no fechamento de uma órbita de teste é desprezível: erro de fase 1,030 × 10⁻³ rad com softening ligado contra 1,011 × 10⁻³ rad desligado.

O integrador, e como saber se ele está certo

Integrar significa: dado onde os corpos estão e a aceleração que sentem, descobrir onde estarão um instante depois. O método mais simples é o de Euler, que avança a posição usando a velocidade atual e a velocidade usando a aceleração atual. É de primeira ordem, e na prática as órbitas precessam visivelmente e a energia deriva.

O método usado aqui é o Velocity-Verlet, que custa quase o mesmo, já que basta guardar a aceleração do passo anterior, e é de segunda ordem:

Além da ordem, ele tem uma propriedade que importa muito em gravitação: é simplético. A energia não deriva sistematicamente para cima ou para baixo, ela oscila em torno do valor certo. Uma órbita fechada continua fechada depois de milhares de voltas.

A parte que realmente prova alguma coisa

"Parece funcionar" não é verificação. A forma correta de testar um integrador é medir sua ordem de convergência: se ele é de segunda ordem, dividir o passo por k tem que dividir o erro por k².

redução do passo:  35×
redução do erro:   1222×
esperado (35²):    1223

O integrador tem exatamente a ordem de convergência que a teoria prevê. Essa é uma afirmação que podia ter saído falsa, e é isso que a torna útil.

O passo de tempo é do cenário, não do programa

Um detalhe de engenharia com consequência física direta: não existe um passo de tempo global. Cada cenário publica o seu, porque uma binária apertada precisa de um passo muito menor que um disco de 20 UA, e reusar o do disco quebraria a binária visivelmente.

O critério vem da órbita mais interna, que é sempre a mais rápida. No disco padrão ela leva cerca de 0,35 ano, e o Velocity-Verlet quer bem mais de 20 passos por órbita:

FIXED_DT = 0,0035 ano   cerca de 100 passos na órbita interna

É também por isso que os cenários de buraco negro rodam muito mais devagar no seu navegador: perto de um corpo de 10⁶ M☉ o passo necessário é minúsculo, e cada segundo de relógio cobre muito menos tempo simulado.

Barnes-Hut: trocar exatidão por tempo

Calcular a força entre todos os pares de n corpos custa O(n²). Com 800 corpos são 319.600 pares por passo, e a 100 passos por órbita isso vira muita conta.

O algoritmo de Barnes-Hut resolve isso com uma observação simples: um aglomerado distante de corpos puxa aproximadamente igual a uma única massa pontual no centro de massa dele. Divide-se o espaço num octree, em que cada cubo se subdivide em oito até sobrar um corpo por célula, e ao calcular a força sobre um corpo percorre-se a árvore decidindo célula a célula pelo critério de abertura:

Isso derruba o custo para O(n log n). O preço é erro, e o erro é medido, não estimado:

Aquela última linha é outro teste falsificável: com θ tendendo a zero, a árvore deixa de aproximar qualquer coisa e tem que reproduzir a soma direta. Reproduz, com quinze casas.

A honestidade desconfortável

Barnes-Hut é a otimização mais citada em simulação de N corpos. Nesta escala, porém, o ganho medido foi:

800 corpos:    ~1,1×   o cruzamento fica bem aqui
~2000 corpos:  ~2,2×   a partir daqui compensa de verdade

Com a contagem padrão, a árvore quase não paga o próprio custo de construção. Ela está lá porque o número de corpos é ajustável e o usuário pode pedir 4000, mas seria desonesto vendê-la como a razão de a simulação ser rápida.

Por que toda estrela é somada exata

Aqui está uma decisão que vale mais que a otimização inteira. A árvore de Barnes-Hut é aplicada apenas aos corpos leves. Toda estrela é somada exatamente, corpo a corpo, fora da árvore.

O raciocínio: as estrelas carregam quase toda a massa do sistema. A órbita de cada corpo leve é, essencialmente, a órbita dele em torno desses poucos termos. Aproximar as estrelas por um multipolo colocaria o erro justamente onde ele domina o resultado, o que seria economizar no lugar errado. E como existem de um punhado a algumas centenas de estrelas, o custo é O(S·n) e exato, contra O(n log n) e aproximado para todo o resto.

As duas passagens tocam pares disjuntos: a árvore nunca vê uma estrela, e a soma direta cobre todo par que envolve uma. Nada é contado duas vezes, nada é esquecido.

Um buraco negro também entra nessa lista, e o comentário no código faz questão de dizer que isso não é uma afirmação de que ele brilha. É apenas o conjunto da soma direta, e um buraco é o objeto mais pesado de qualquer sistema que o contenha.

O momento que a árvore quebra, e a correção

Há um efeito colateral sutil. As forças calculadas pela árvore não são simétricas entre pares: o corpo A sente a célula que contém B, mas B sente uma célula diferente que contém A. A terceira lei de Newton deixa de valer exatamente, e sobra uma força líquida espúria sobre o sistema inteiro.

A correção é elegante. Some essa força líquida e subtraia-a como uma aceleração uniforme de todos os corpos que passaram pela árvore. Isso restaura a condição abaixo sem alterar nenhuma aceleração relativa, e portanto sem tocar em nenhuma órbita:

força líquida bruta da árvore:   6 × 10⁻³
depois da correção:              3,9 × 10⁻¹⁶

medido com 800 corpos, após 531 fusões

Composição: fusão em cadeia alfa

Cada corpo carrega uma composição química. Não um rótulo, mas frações de massa de oito espécies da cadeia alfa, que é a sequência de fusão que governa a vida de uma estrela:

H → He → C, O → Ne, Mg → Si → Fe

cada etapa precisa de uma temperatura central maior que a
anterior. O ferro é o fim da linha, porque fundir ferro
consome energia em vez de liberar.

A temperatura central de cada corpo é recalculada a cada passo pela fórmula virial da introdução, e ela decide qual etapa está acesa. Isso é um equilíbrio, não uma integração: o tempo térmico de uma estrela e o tempo dinâmico diferem por ordens de grandeza. O tempo de Kelvin-Helmholtz do Sol é de dezenas de milhões de anos, o dinâmico é de meia hora. Na escala do passo, a estrutura térmica é quase estática, e calcular o equilíbrio sai mais barato e mais estável que evoluí-lo.

O bug que só um teste pegaria

A cadeia como especificada inicialmente tinha estados absorventes: nitrogênio e magnésio não tinham saída. Um núcleo a 5 × 10⁹ K, quente o bastante para chegar ao ferro, terminava com 69% de magnésio e parava ali.

A correção foi adicionar as duas transições que faltavam: N → O, já que o nitrogênio é catalisador do ciclo CNO mas é destruído por captura alfa quando o hélio acende, e Mg → Si. O comentário que ficou no código resume o invariante:

O ferro é o único estado absorvente desta cadeia, e esse é justamente o ponto.

Outro caso do mesmo tipo: o limiar de ignição, se fosse só por temperatura, acenderia anãs marrons. A estimativa virial erra por um fator de cerca de 4 perto do limite, porque em objetos muito pequenos quem segura a estrela é a pressão de degenerescência dos elétrons, que não custa calor. Hoje a ignição é condicionada à massa, e o resultado é o correto:

0,079 M☉  →  anã marrom, T_c = 1,22 × 10⁷ K, não queima nada
0,080 M☉  →  estrela, queima

As licenças poéticas

Esta é a seção que eu mais queria escrever, porque é a que costuma ficar de fora. Toda simulação mente em algum lugar. A diferença entre um modelo e uma animação é se as mentiras estão marcadas e se o custo delas foi calculado.

1. O raio de captura, e a maior mentira do projeto

Em escala real, um planeta tem raio de ordem 10⁻⁵ UA contra separações orbitais de ordem 1 UA. A seção de choque geométrica é tão absurdamente pequena que um disco real leva de 10⁷ a 10⁸ anos para formar planetas. Numa aba de navegador, nada jamais aconteceria.

Então o raio de captura é multiplicado por 800. Isso comprime a formação planetária em alguns milhares de anos simulados. O que salva isso de ser um truque é o que vem junto:

• A dinâmica não é tocada. Órbitas, ressonâncias e espalhamento continuam exatos. Só o tamanho do alvo é exagerado.
• O custo foi medido: com fator 800, o raio de captura de um embrião a 1 UA é 3,3 vezes a própria esfera de Hill dele. Ou seja, corpos do disco interno se capturam bem fora do território onde deveriam mandar.
• A alternativa honesta está escrita ao lado: com fator 250 o raio de captura fica dentro da esfera de Hill em todo o disco interno, ao custo de 11 massas terrestres em 150 anos em vez de uma super-Terra em 27.

E a fronteira entre licença e física está demarcada em uma linha do código. O foco gravitacional, que aumenta a seção de choque por um fator (1 + v_esc²/v_rel²) e importa quando a velocidade relativa é baixa, é física real, não licença.

Há também uma mentira que foi conscientemente não cometida: a estrela central absorve a 10 raios solares, não a 800. O comentário explica por quê, já que 800 vezes daria 3,7 UA e ela comeria o disco inteiro.

2. O tamanho desenhado

Um raio solar é 0,00465 UA. Num disco de 20 UA, um planeta em escala real ocupa alguns milionésimos da cena: a tela seria preta. Os corpos são desenhados inflados, numa curva compressiva de expoente 0,55, para que um Júpiter não engula uma Terra pelo fator que engoliria de verdade.

O valor dessa escala foi ajustado por medição, não por gosto. Com o valor que parecia razoável de início, a faixa dinâmica inteira colapsava:

Ruptura por maré: o teste contra a natureza

Quando uma estrela passa perto demais de um buraco negro supermassivo, a diferença de gravidade entre o lado de perto e o lado de longe supera a gravidade própria dela e ela se despedaça. Metade dos destroços fica ligada e volta, a outra metade escapa. É um evento real, observado dezenas de vezes, e tem uma assinatura característica.

A pergunta de projeto era se dava para simular isso sem hidrodinâmica. A resposta veio de comparar o raio de maré com o comprimento de suavização, aquela régua de honestidade da seção anterior:

Num buraco supermassivo, a região onde a estrela se despedaça é centenas de vezes maior que a escala em que a simulação já deixou de resolver. Ou seja: não é preciso fingir. Ao cruzar o raio de maré, a estrela é substituída por 600 corpos que carregam a massa e a composição dela, e a física que já existia faz o resto. As partículas da frente orbitam mais rápido que as de trás, e o esticamento aparece sozinho.

A espaguetificação deixa de ser desenhada e passa a ser consequência.

Não funcionou de primeira, e o motivo importa

Na primeira tentativa, os fragmentos foram distribuídos ao longo de um raio solar. Isso os deixou a 1,5 × 10⁻⁵ UA uns dos outros, 65 vezes mais perto que o comprimento de suavização. E abaixo do softening, o kernel de Plummer não é gravidade: é uma mola linear.

O resultado foi imediato e inequívoco. A distribuição de energia colapsou em quatro picos em menos de 100 passos, e o expoente ajustado deu −0,41, com r² = 0,04. Não era uma imprecisão. Era o modelo descrevendo outra coisa.

A correção foi alargar o fluxo até o espaçamento médio ser alguns comprimentos de suavização. E foi classificada, no próprio código, como licença poética da mesma espécie do fator 800, porque é exatamente isso que ela é.

O resultado

A assinatura observacional de um evento de ruptura por maré é a taxa de retorno dos destroços seguindo uma lei de potência:

O expoente sai cerca de 2,5% mais íngreme que o teórico, de forma consistente, e a causa foi diagnosticada em vez de varrida para baixo do tapete: faltam retornos tardios no fim da janela, o que inclina o ajuste. O ajuste puramente kepleriano sobre as energias iniciais dá −1,663.

Onde não funciona

Com um buraco negro estelar de 10 M☉, a conservação continua exata e a divisão ainda é meio a meio, mas o resultado físico não presta. Nessa massa, o raio de absorção da própria estrela é cinco vezes maior que o raio de maré, então a "ruptura" dispara a nove raios de maré, longe de onde deveria. O cruzamento útil fica em torno de 10⁴ M☉.

O que foi deliberadamente recusado

A decisão mais madura do projeto inteiro é uma coisa que não foi implementada.

Perto de um buraco negro, a gravidade newtoniana está errada. Existe um atalho conhecido para capturar os efeitos relativísticos sem resolver Einstein, o potencial pseudo-newtoniano de Paczyński-Wiita:

Seria fácil de adicionar e renderia uma boa linha no README. Foi recusado, e o motivo está gravado no código com o pedido explícito de que ninguém o adicione por engano depois. A conta:

buraco de 10 M☉:
  r_s          = 1,97 × 10⁻⁸ UA
  ISCO (3·r_s) = 6 × 10⁻⁷ UA
  softening    = 1 × 10⁻³ UA

a ISCO fica 4 a 5 ordens de grandeza ABAIXO da escala
em que a força já é deliberadamente errada

Seria decoração pintada sobre um lugar que a simulação não resolve, ao custo de um caso especial no laço mais quente do programa e da conservação exata de momento. Para efeito de comparação: num buraco de 10⁶ M☉ a razão r_s / softening é 19,7, e aí sim o horizonte é resolvível. Por isso buracos negros supermassivos entraram e a relatividade não.

Pelo mesmo critério ficou de fora a hidrodinâmica de verdade (SPH). Precisaria de pressão, choques e uma equação de estado, e mudaria o caráter do projeto inteiro. Os fragmentos da ruptura são massas gravitacionais, e o texto do próprio cenário diz isso ao usuário.

O resumo

Se der para levar uma coisa só deste texto, que seja o critério que aparece em todas as seções: escolher a escala em que você está disposto a errar, e então recusar tudo o que acontece abaixo dela.

O comprimento de suavização define essa escala. Ele justifica o softening, justifica a recusa da relatividade, condena a primeira versão do fluxo de destroços e decide quais buracos negros valem a pena. Um único número, aplicado com consistência, faz mais pela honestidade de um modelo do que qualquer quantidade de física adicionada por cima.

A simulação está em cluxnei.dev/universe-simulation/ e roda inteira no seu navegador, sem servidor e sem instalação. O código está no GitHub, incluindo todos os comentários citados aqui, que são a melhor parte.

kelvin@cluxnei:~/blog$