Resonancias de J1903+0327

Parte I — de qué se trata

El acoplamiento

Si la materia oscura es un campo ultraliviano —masas de 10⁻²² eV, unas 30 órdenes por debajo de un electrón— su longitud de onda de de Broglie es del tamaño de una galaxia y el campo se comporta como una onda clásica que oscila con frecuencia igual a su masa, φ ∝ cos(m t). No es un gas de partículas: es un campo coherente que baña al sistema solar y a cualquier púlsar de la galaxia.

Acá consideramos un campo vectorial (espín 1) acoplado al número bariónico menos leptónico, B−L. Ese acoplamiento le da a cada cuerpo una carga proporcional a su número de neutrones, y por lo tanto una fuerza extra que oscila en el tiempo:

f⃗(t) = − 𝓕 sin(m t̃ + Υ₀) n̂ 𝓕 = g · (δ_c / m_N) · √(2 ρ_DM) · ρ La fuerza apunta en la dirección n̂ del campo y oscila con frecuencia m. Su amplitud lleva el acoplamiento g —lo que queremos acotar—, la asimetría de carga entre las dos estrellas (δ_c), la densidad local de materia oscura, y la amplitud local del campo ρ, que es estocástica.

Como g aparece lineal en la fuerza, todo el problema es: ¿cuán grande puede ser g sin que ya hubiéramos visto el efecto?

El efecto sobre un púlsar binario

Un púlsar de milisegundo es un reloj: emite un pulso por rotación con una estabilidad comparable a la de un reloj atómico. Si está en un sistema binario, el tiempo de llegada de cada pulso depende de dónde está en su órbita. Medimos los tiempos de llegada (TOAs) durante años, les restamos el modelo orbital que mejor ajusta, y lo que queda son los residuos: microsegundos.

La fuerza extra perturba la órbita. La forma estándar de escribirlo son las ecuaciones osculantes (de Gauss), que dan cómo evoluciona cada elemento orbital bajo una perturbación. Por ejemplo, la excentricidad:

ė = ( √(1−e²) / (a ω_b) ) [ f_r sinθ + f_θ (cosθ + cos E) ] Con f_r y f_θ las componentes radial y tangencial de la fuerza, θ la anomalía verdadera y E la excéntrica. Hay una ecuación así para cada uno de los cinco elementos que el timing mide.

Cada elemento perturbado mueve los residuos a través de la derivada que el ajuste de timing ya conoce:

r(t) = ∑_k ( ∂r / ∂p_k ) δp_k(t) con p_k = los parámetros orbitales P_b, x, e, ω, T₀. Eso define el template: la forma que tendría el residuo si el campo existiera. Es lineal en g.

La inferencia

Con el template en mano, el resto es estadística estándar. Los residuos son gaussianos con una matriz de covarianza C que incluye el ruido blanco medido, el ruido rojo del púlsar y la marginalización sobre el modelo de timing:

ln L(d | θ) = −½ (r − h(θ))ᵀ C⁻¹ (r − h(θ)) − ½ ln det(2πC) p(θ | d) ∝ L(d | θ) π(θ) ← Bayes El vector θ tiene 7 componentes: el acoplamiento g, la amplitud ρ y la fase γ del campo, dos ángulos para su dirección, y dos parámetros de ruido rojo. Los priors π son los físicos: Rayleigh en ρ, uniforme en la fase, isótropo en la dirección.

De los siete parámetros sólo nos interesa uno. Marginalizamos los otros seis, lo que deja la evidencia marginal Z(g), y de ahí sale la posterior de g y su percentil 95: la cota.

Z(g) = ∫ L(d | g, ψ) π(ψ) dψ ψ = los otros seis p(log g | d) ∝ 10^{log g} · Z(g) → cota al 95 %

Hacer esa integral de seis dimensiones es lo caro. El método de siempre es una cadena de Markov (MCMC): se pasea por el espacio de los siete parámetros y el histograma de dónde estuvo aproxima la posterior. Es lo que veníamos usando, una cadena por cada masa del barrido.

Las resonancias

Acá está lo interesante. La fuerza oscila con frecuencia m; la órbita, con ω_b. Cuando m se acerca a un múltiplo de ω_b, el efecto deja de promediarse a cero órbita a órbita y se acumula: los elementos orbitales derivan linealmente en el tiempo. Sobre T_obs = 10.7 años eso amplifica la señal en un factor ω_bT_obs ≈ 260.

δa|_sec = (2 𝓕 sinϑ · t̃ / ω_b) [ J'_N(Ne) cosφ cosΥ₀ − (√(1−e²)/e) J_N(Ne) sinφ sinΥ₀ ] La respuesta secular en el semieje, justo en m = Nω_b. El factor t̃ es el crecimiento lineal. Los J_N son funciones de Bessel: en una órbita excéntrica resuenan TODOS los armónicos, no sólo el fundamental. Hay una expresión así para cada elemento orbital.

J1903+0327 es el binario más excéntrico del wideband de NANOGrav 15 yr, con e = 0.4367. Por eso lo elegimos: tiene seis resonancias utilizables en vez de una.

Izquierda: el residuo inducido a lo largo de cuatro órbitas, mucho mayor en resonancia. Derecha: la amplitud por órbita crece lineal en resonancia y es plana fuera. — **La señal.** A la izquierda, cómo se ve el residuo inducido en cuatro órbitas. A la derecha, su amplitud órbita por órbita: en resonancia crece lineal desde la época de referencia; fuera, se queda chata. Es el mismo acoplamiento en los tres casos — sólo cambia la masa del campo.

El pico es angosto. Su ancho relativo es 1/(n ω_b T_obs), o sea 0.0039 para el fundamental: cuatro partes en mil. Por eso el barrido en masa se mide en anchos de pico, no en eV.

Parte II — el problema

Un antipico donde debería haber un pico

Al hacer el barrido en masa aparecía algo raro. Acercándose a la resonancia la cota mejora, como corresponde: hay más señal. Pero justo en el centro exacto m = nω_b la cota empeora, y bastante. Un antipico, apuntando para el lado equivocado.

Y en ese mismo punto la cadena de Markov se porta mal. La forma de medirlo es el ESS, effective sample size: una cadena de 40 000 pasos no da 40 000 muestras independientes, porque cada paso está correlacionado con el anterior. El ESS es el número de muestras efectivamente independientes. Con ESS de unas pocas decenas, un percentil 95 es puro ruido.

La pregunta del proyecto

¿El antipico es físico, o es que la cadena no llega?

Con una sola herramienta no se puede contestar: si el método falla justo donde queremos mirar, no hay con qué contrastar. Hacía falta calcular la misma posterior por un camino que no muestree.

Parte III — el desarrollo

Primero: ¿el template está bien?

Antes de dudar del muestreo hay que descartar que el problema sea la señal. Las expresiones resonantes de arriba son formas cerradas, obtenidas desarrollando en series de Fourier–Bessel. Se validan contra lo único que no puede estar de acuerdo por casualidad: las ecuaciones osculantes exactas, evaluadas sin ninguna expansión sobre la órbita kepleriana.

Las tasas en serie de Bessel reproducen el lado derecho de Gauss crudo. ~10⁻¹⁵
Integradas durante 15 años, las formas cerradas coinciden con la integración numérica de las osculantes, en los ocho canales orbitales, en resonancia y fuera. al piso de la cuadratura: ≤6×10⁻⁷ (10⁻⁹–10⁻⁸ en las resonancias de J1903)
El template que usa el pipeline reproduce esas formas cerradas. 1.6×10⁻¹²

El template no es el problema.

El método sin muestreo: tres reducciones exactas

La idea es hacer la integral de seis dimensiones de frente en vez de muestrearla. Suena imposible, pero el problema tiene tres propiedades que lo permiten. Ninguna de las tres es una aproximación.

1. Toda la información de una masa cabe en 27 números

La log-verosimilitud es exactamente cuadrática en la señal, y la señal es bilineal en los ángulos: se descompone en seis templates fijos que no dependen de ellos. Entonces, a ruido fijo:

h = g ρ ∑_k c_k(ϑ,γ,φ) U_k ln L = ln L₀ + (gρ) c·d − ½ (gρ)² cᵀGc con d_k = U_kᵀC⁻¹r (6 números) y G_kl = U_kᵀC⁻¹U_l (21). Calculados una vez, dan la verosimilitud en cualquier punto del espacio de parámetros, sin volver a tocar los datos.

2. La amplitud y su fase son un gaussiano, y salen en forma cerrada

El prior de la amplitud es una Rayleigh y el de la fase es uniforme. Eso quiere decir que el par (ρcosγ, ρsinγ) es un gaussiano bidimensional. Y la señal es lineal en ese par. Una integral gaussiana con un exponente lineal más cuadrático tiene forma cerrada: esas dos dimensiones se van exactamente.

No es una comodidad, es lo que hace que el cálculo funcione: intentando hacer esa integral por cuadratura numérica, la cota no converge — en resonancia el integrando es un pico agudísimo justo en esa dirección.

3. Queda una sola integral, sobre la esfera

Lo que sobrevive es un promedio sobre la dirección del campo, y eso se hace con una regla producto que es exacta para polinomios esféricos. Converge geométricamente, y coincide con un Monte Carlo uniforme independiente.

El resultado: la misma posterior, sin ningún sorteo. Y sin la limitación de la cadena, se puede evaluar también el centro exacto del pico.

Cómo se verificó

Esto es una secuencia: si el primer eslabón falla, los demás no significan nada.

El template contra las ecuaciones osculantes exactas, en los tres pasos de arriba. al piso de la cuadratura numérica; 1.6×10⁻¹² el del pipeline
La bilinealidad, que es la hipótesis de la que cuelga todo: 60 direcciones al azar por masa. peor error relativo 2.4×10⁻¹⁴
El núcleo (d, G) contra la verosimilitud del pipeline, en 25 puntos al azar por masa. peor |Δln L| = 1.4×10⁻¹⁰ — el eslabón crítico
La integral gaussiana, dirección por dirección, contra cuadratura numérica bruta con 120 000 nodos en la fase. 7×10⁻¹⁵
La malla esférica, por refinamiento con extrapolación del límite, y contra tres Monte Carlo de 10⁶ direcciones. producción a <0.01 dex del límite

Parte IV — resultados

Las dos curvas

Sobre las 180 masas en común, los dos métodos coinciden con una mediana de −0.062 dex y un RMS de 0.250 dex. Y el desacuerdo se achica al alargar la cadena, que es la firma de que lo que queda es ruido de muestreo:

iteraciones por cadena	ESS mediano	RMS contra la cuadratura
10 000	42.6	0.362 dex
40 000	141.6	0.250 dex

La mejor cota del barrido

log₁₀ g < −22.300

al 95 % de credibilidad, en 4.9679×10⁻²² eV. El barrido cubre 186 masas entre 5.0294×10⁻²³ eV y 3.0176×10⁻²⁰ eV.

El centro del pico

Seis paneles, uno por armónico, con la cota contra el apartamiento del centro medido en anchos de pico. Perfil simétrico con mínimos a ±3.2 anchos y un pico hacia arriba en el centro. — **Zoom en cada armónico.** El apartamiento va en anchos de pico. La estrella naranja es el centro exacto, que sólo la cuadratura evalúa. El perfil es simétrico, con la mejor cota a ±3.16 anchos y un hundimiento claro justo en el centro.

armónico	m [eV]	mejor del pico	centro exacto	penalidad [dex]
n = 1	5.0294×10⁻²²	−22.300	−20.670	1.63
n = 2	1.0059×10⁻²¹	−21.913	−20.001	1.91
n = 3	1.5088×10⁻²¹	−21.582	−19.620	1.96
n = 4	2.0117×10⁻²¹	−21.274	−19.218	2.06
n = 5	2.5147×10⁻²¹	−20.979	−18.796	2.18
n = 6	3.0176×10⁻²¹	−20.693	−18.649	2.04

Por qué: una fase que el púlsar no ve

El campo tiene una orientación (dos ángulos) y una fase. Con la cuadratura se puede hacer algo que con la cadena no: fijar todos los ángulos y calcular la cota para esa configuración exacta, marginalizando sólo la amplitud — que también sale en forma cerrada. Y el resultado es claro: el problema no está en hacia dónde apunta el campo, sino en su fase.

Para cada orientación, la señal depende de la fase a través de dos combinaciones: una que el púlsar registra bien y otra — la fase ciega — que registra mal. Es como empujar una hamaca: en fase, cada empujón suma; a 90°, lo que empujás en media órbita lo devolvés en la otra. Lejos del centro las dos fases se ven parecido; justo en resonancia la ciega se vuelve casi invisible.

La cota al 95 % con todos los ángulos fijos, contra la distancia a la fase ciega en escala logarítmica. En el centro exacto la curva está alrededor de −23 lejos de la fase ciega y sube hasta −20.1 dentro de unas décimas de grado; la marginalizada, −20.67, queda arriba de casi toda la curva. A 3.2 anchos la curva es plana en −22.5. — **La cota con los ángulos fijos, contra la fase.** En el centro exacto (naranja) la cota típica es −23.06, pero a menos de unas décimas de grado de la fase ciega sube hasta −20.12. A 3.2 anchos del centro (gris) la fase casi no importa. Las líneas punteadas son las cotas marginalizadas.

La marginalizada la fija la fase ciega

0.25 %

del espacio de ángulos da, en el centro exacto, una cota más débil que la marginalizada (−20.67). En todo el resto la cota a ángulos fijos es más fuerte. A 3.2 anchos esa fracción es 3.4 %.

Cómo una franja tan angosta mueve la cota dos décadas: un embudo

La cuadratura también da la posterior exacta en el plano (acople, fase), sin sortear nada. En el centro tiene forma de embudo: a acople grande sólo sobrevive una franja de fase cada vez más angosta. Pero se angosta como 1/g, que es justo lo que compensa el prior plano en g: la posterior de log g queda con una meseta de varias décadas, y el percentil 95 se va para arriba.

Arriba, la posterior de log g: en el centro una meseta ancha y un histograma de las cadenas ruidoso; a 3.2 anchos una campana que el histograma sigue bien. Abajo, la posterior exacta en el plano log g contra distancia a la fase ciega: en el centro una banda diagonal que baja hasta centésimas de grado a acople grande; a 3.2 anchos, una mancha compacta. Las muestras MCMC, en azul, se superponen. — **La posterior exacta en el plano (acople, fase)**, con las muestras de seis cadenas MCMC encima (azul). Arriba, la marginal en log g.

Y ese embudo explica la cadena. Es la geometría clásica en la que un muestreador no mezcla: los pasos que sirven en la boca son enormes para el cuello. Las cadenas lo siguen — en el centro, sus muestras de mayor acople están a una mediana de 0.08° de la fase ciega, contra 17° a 3.2 anchos — pero no lo recorren bien. Con 100 000 pasos, en el centro el ESS mediano es 70 y seis semillas dan cotas separadas por 1.27 dex; a 3.2 anchos, ESS 376 y 0.40 dex.

La prueba: sacando la franja, el antipico desaparece

Si la franja ciega es la causa, excluirla del prior tiene que borrar el antipico, y excluir la misma franja fuera del centro casi no tiene que cambiar nada.

El antipico, cota en el centro menos cota a 3.2 anchos, contra el ancho de la franja de fase excluida, para los seis armónicos. Sin excluir nada vale entre 1.6 y 2.2 dex; al excluir unos pocos grados cae a cero en todos. — **El antipico contra el ancho de la franja excluida** alrededor de la fase ciega, para los seis armónicos (naranja, n = 1). Arriba, qué fracción del *prior* de la fase se saca. Cada punto compara dos cotas calculadas con los mismos ajustes.

Con todas las fases el antipico vale 1.7 a 2.3 dex. Sacando |δ| < 5° — 5.6 % del prior de la fase — queda en 0.2 dex, y sacando |δ| < 10° en 0.06 dex, en los seis armónicos por igual. El control a 3.2 anchos casi no se mueve.

De dónde sale la fase ciega

Izquierda: el acople al que la fase fuerte y la ciega darían señal/ruido 1, y la cota, contra el apartamiento del centro. La fuerte es plana; la ciega y la cota suben juntas hacia el centro. Derecha: la brecha entre la fase ciega y la fuerte crece hacia el centro tanto con el modelo completo como con la señal sola, sin ningún ajuste. — **Izquierda:** la cota sigue a la fase ciega; la fuerte ni se entera de la resonancia. **Derecha:** la brecha entre las dos fases, con el modelo completo y con la señal sola, sin ningún ajuste.

La ceguera ya está en la señal: sin ningún ajuste de timing ni modelo de ruido, acercarse al centro vuelve a la fase ciega 1.5 dex más difícil de ver. Lo poco que queda de ella es lento, y el modelo de ruido de baja frecuencia (ruido rojo y variaciones de DM) se lo come: eso agrega otros 0.9 dex en el centro. Los parámetros orbitales del ajuste, incluida la deriva ω̇, no participan: sacarlos del modelo no cambia nada.

Conclusión

El antipico es físico, no un artefacto del muestreo. Aparece igual en un método que no muestrea y que evalúa la misma verosimilitud que la cadena. Y tiene una causa concreta y comprobable: en resonancia exacta, para cada orientación del campo hay una fase casi invisible, y esa franja angosta de fase es la que fija la cota. Sacándola, el antipico desaparece en los seis armónicos.
La cadena, ahí, no servía para decidirlo. La posterior tiene forma de embudo, y en un embudo un muestreador no mezcla: seis semillas dan seis respuestas separadas por más de un dex. No es que diera mal: no daba nada reproducible.
Ahora se puede explorar esa región. El método sin muestreo evalúa el centro exacto de cada pico — donde la cadena no entra — y coincide con ella donde sí mezcla.

Y de yapa: donde una cadena tarda decenas de minutos por masa, la cuadratura tarda segundos. El barrido completo de 186 masas corre en minutos.

Lo que no se afirma

Esta es la cota de un púlsar. La envolvente conjunta de los 50 sistemas del wideband es bastante más fuerte fuera de resonancia.
La derivada secular ω̇ está fiteada en el modelo y se deja en la matriz de diseño, sin congelar. Sacarla (y también P_b) no cambia la sensibilidad de ninguna de las dos fases, medido en el centro y a 3.16 anchos del fundamental.
La prueba de la franja cambia el prior de la fase: es un diagnóstico de qué fija la cota, no una cota alternativa.
El máximo de la log-razón de verosimilitud sobre todo el barrido es 5.52, que para 186 masas y una marginalización angular no es una detección. No se reporta ninguna.
No se afirma reproducibilidad bit a bit entre plataformas. La cuadratura es determinista — no hay ningún sorteo; el MCMC no.

Procedencia

Ningún número de esta página está tipeado a mano: todos salen de resultados/numeros.json, que produce scripts/comparar_metodos.py, y scripts/verificar_provenance.py chequea que sigan coincidiendo. El reporte en PDF usa el mismo registro por macros de LaTeX.

qué	lo produce	a partir de
curva por cuadratura	`scripts/barrido_cuadratura.py`	`src/cuadratura.py`, `data/`
curva por MCMC	`scripts/barrido_mcmc.py`	`src/pipeline_vectorial.py`, `data/`
las cadenas del centro, con muestras	`scripts/cadenas_ciegas.py`	6 semillas × 2 masas
las direcciones ciegas	`scripts/nucleos_ciegas.py`, `scripts/direcciones_ciegas.py`	`src/cuadratura.py`, `data/`
la grilla de masas	`src/grilla_masas.py`	ω_b y T_obs medidos sobre el `.tim`
las figuras	`comparar_metodos.py`, `figuras_explicativas.py`, `direcciones_ciegas.py`	los `.npz`
las validaciones	`tests/`	—

Los registros legibles por máquina están en structure/: claims.yaml mapea cada afirmación a su evidencia, figures.yaml cada figura a su script, scripts.yaml el camino crítico, y decisiones.md cada punto donde había más de un camino razonable, con el motivo de la elección. bash scripts/reproduce.sh corre todo, de los datos a las figuras.