Multi-Messenger Analysis

Overview

Multi-messenger astronomy combines observations from different “messengers” (electromagnetic radiation, gravitational waves, neutrinos, cosmic rays) to provide a more complete picture of astrophysical transients. redback provides dedicated infrastructure for multi-messenger analysis through the MultiMessengerTransient class.

The key advantage of multi-messenger analysis is that different messengers can probe different aspects of the same physical system, and shared parameters can be jointly constrained by all available data. For example, in a binary neutron star merger:

Gravitational waves constrain masses, spins, and distance
Kilonova emission (optical/IR) constrains ejecta properties and viewing angle
Gamma-ray burst afterglow (X-ray/radio) constrains jet properties and viewing angle
Viewing angle and distance are shared across all messengers

The `MultiMessengerTransient` class

Basic Usage

The MultiMessengerTransient class provides a high-level interface for combining data from multiple messengers:

import redback
from redback.multimessenger import MultiMessengerTransient, MultiMessengerLikelihood

# Create transient objects for each messenger
optical_transient = redback.get_data.get_kilonova_data_from_open_transient_catalog_data(
    transient='AT2017gfo'
)

# Assuming we have X-ray and radio data as well
xray_transient = redback.transient.Transient(...)
radio_transient = redback.transient.Transient(...)

# Create multi-messenger object
mm_transient = MultiMessengerTransient(
    optical_transient=optical_transient,
    xray_transient=xray_transient,
    radio_transient=radio_transient,
    name='AT2017gfo'
)

Supported Messengers

The class supports multiple types of data:

Electromagnetic transients (optical, X-ray, radio, UV, infrared) - Provided as redback.transient.Transient objects - Can use any redback transient class (Kilonova, Afterglow, Supernova, etc.)
Gravitational waves - Provided as pre-constructed bilby.gw.likelihood objects - Use standard bilby.gw workflow to create the likelihood
Neutrinos - Provided as custom bilby.Likelihood objects
Custom likelihoods - Any custom likelihood inheriting from bilby.Likelihood

Joint Analysis

The fit_joint() method performs parameter estimation using all available data:

import bilby

# Define models for each messenger
models = {
    'optical': 'two_component_kilonova_model',
    'xray': 'tophat',  # afterglow model
    'radio': 'tophat'
}

# Define priors (including shared parameters)
priors = bilby.core.prior.PriorDict()
priors['viewing_angle'] = bilby.core.prior.Uniform(0, 1.57)
priors['redshift'] = 0.01  # Fixed

# Kilonova-specific parameters
priors['mej_1'] = bilby.core.prior.Uniform(0.01, 0.1)
priors['vej_1'] = bilby.core.prior.Uniform(0.1, 0.3)
# ... more priors ...

# Afterglow-specific parameters
priors['loge0'] = bilby.core.prior.Uniform(50, 54)
priors['logn0'] = bilby.core.prior.Uniform(-3, 2)
# ... more priors ...

# Specify model kwargs
model_kwargs = {
    'optical': {'output_format': 'magnitude'},
    'xray': {'output_format': 'flux_density', 'frequency': xray_freq},
    'radio': {'output_format': 'flux_density', 'frequency': radio_freq}
}

# Run joint analysis
result = mm_transient.fit_joint(
    models=models,
    priors=priors,
    shared_params=['viewing_angle', 'redshift'],
    model_kwargs=model_kwargs,
    nlive=2000,
    sampler='dynesty',
    outdir='./results_joint'
)

Joint results are returned as MultiMessengerResult objects. They keep standard result behaviour such as posterior access, evidence values, and plot_corner(), including the same corner-title formatting used by RedbackResult. They do not automatically reconstruct a single transient for plot_lightcurve(), plot_spectrum(), or related helpers; use the relevant transient object and redback.analysis plotting functions when plotting one component of the joint fit.

Upper Limits

Photometric, X-ray, radio, UV, and infrared transient objects passed directly to MultiMessengerTransient use redback likelihood construction. If a transient includes finite upper-limit values via its detections array, MultiMessengerTransient uses GaussianLikelihoodWithUpperLimits automatically for that messenger. Upper limits with missing numerical limit values are excluded with a warning, matching the single-transient fitting behaviour.

Shared Parameters

The shared_params argument validates parameters that are shared across messengers. A shared parameter must appear in at least two likelihoods using the same sampled parameter name; otherwise fit_joint() raises an error rather than silently running an unshared fit. If two models use different native names for the same physical quantity, provide parameter_mappings to expose a common sampled name.

Common shared parameters include:

viewing_angle: Observer’s viewing angle (affects both EM and GW signals)
luminosity_distance: Distance to source
redshift: Cosmological redshift
time_of_merger: Reference time for all emissions

For example, an afterglow model may call the viewing angle thv, while an optical model or GW likelihood may call it viewing_angle:

result = mm_transient.fit_joint(
    models={'optical': optical_model, 'xray': 'tophat'},
    priors=priors,
    shared_params=['viewing_angle'],
    parameter_mappings={'xray': {'viewing_angle': 'thv'}},
    model_kwargs=model_kwargs
)

Individual Fits for Comparison

For comparison with joint analysis, you can fit each messenger independently:

individual_models = {
    'optical': 'two_component_kilonova_model',
    'xray': 'tophat',
    'radio': 'tophat'
}

# Define separate priors for each messenger
optical_priors = bilby.core.prior.PriorDict()
# ... optical-specific priors ...

xray_priors = bilby.core.prior.PriorDict()
# ... X-ray-specific priors ...

individual_priors = {
    'optical': optical_priors,
    'xray': xray_priors,
    'radio': radio_priors
}

# Fit each independently
individual_results = mm_transient.fit_individual(
    models=individual_models,
    priors=individual_priors,
    model_kwargs=model_kwargs,
    nlive=2000
)

# Access individual results
optical_result = individual_results['optical']
xray_result = individual_results['xray']

If an individual fit should use the same sampled parameter name as the joint fit while the native model uses a different name, pass the same parameter_mappings dictionary used by fit_joint().

Joint Spectrum and Photometry Analysis

A common use case in transient astronomy is jointly fitting spectroscopic and photometric data from the same object. This combines:

Photometry: Multi-band lightcurves constraining time evolution and integrated properties
Spectroscopy: Detailed spectral energy distribution constraining temperature, composition, and line features

By fitting both jointly, physical parameters can be better constrained as the data types provide complementary information.

Basic Approach

Since Spectrum and Transient are different classes in redback, the recommended approach for joint spectrum + photometry fitting is to use custom likelihoods:

import redback
from redback.multimessenger import MultiMessengerTransient, MultiMessengerLikelihood
from redback.transient_models import supernova_models, spectral_models

# Load or create photometry data
photometry = redback.transient.Transient(
    time=phot_times,
    magnitude=mags,
    magnitude_err=mag_errs,
    bands=bands,
    data_mode='magnitude'
)

# Load or create spectrum at a specific epoch
spectrum = redback.transient.Spectrum(
    angstroms=wavelengths,
    flux_density=flux,
    flux_density_err=flux_err,
    time='3 days'
)

# Build likelihoods
phot_likelihood = redback.likelihoods.GaussianLikelihood(
    x=photometry.time,
    y=photometry.magnitude,
    sigma=photometry.magnitude_err,
    function=supernova_models.arnett,
    kwargs={'output_format': 'magnitude', 'bands': photometry.bands}
)

spec_likelihood = redback.likelihoods.GaussianLikelihood(
    x=spectrum.angstroms,
    y=spectrum.flux_density,
    sigma=spectrum.flux_density_err,
    function=spectral_models.blackbody_spectrum_at_z,
    kwargs={}
)

# Create multi-messenger object with custom likelihoods
mm_transient = MultiMessengerTransient(
    custom_likelihoods={
        'photometry': phot_likelihood,
        'spectrum': spec_likelihood
    },
    name='joint_spec_phot'
)

Setting Up Priors

Define priors for parameters used by both photometry and spectrum models:

import bilby

priors = bilby.core.prior.PriorDict()

# Shared parameters (constrained by both data types)
priors['redshift'] = 0.01  # Fixed if known
priors['f_nickel'] = bilby.core.prior.Uniform(0.01, 0.5)
priors['mej'] = bilby.core.prior.Uniform(0.1, 3.0)        # ejecta mass
priors['vej'] = bilby.core.prior.Uniform(5000, 20000)     # ejecta velocity
priors['kappa'] = bilby.core.prior.Uniform(0.5, 10.0)
priors['kappa_gamma'] = bilby.core.prior.Uniform(0.01, 1.0)
priors['temperature_floor'] = bilby.core.prior.Uniform(1000, 5000)

# Spectrum-specific parameters (epoch-dependent)
priors['temperature'] = bilby.core.prior.Uniform(3000, 10000)  # at t_spec
priors['r_phot'] = bilby.core.prior.LogUniform(1e13, 1e15)     # photosphere size

Running Joint Analysis

Use bilby.run_sampler directly with a joint likelihood:

import bilby

# Combine likelihoods
joint_likelihood = MultiMessengerLikelihood(
    phot_likelihood,
    spec_likelihood
)

# Run joint fit
result = bilby.run_sampler(
    likelihood=joint_likelihood,
    priors=priors,
    sampler='dynesty',
    nlive=2000,
    outdir='./results_joint_spec_phot',
    label='joint_spectrum_photometry'
)

Multiple Spectra at Different Epochs

If you have spectra at multiple epochs, create separate likelihoods with epoch-dependent parameters:

# Spectra at 1, 3, and 7 days
spectrum_1d = redback.transient.Spectrum(wave, flux_1d, err_1d, time='1 day')
spectrum_3d = redback.transient.Spectrum(wave, flux_3d, err_3d, time='3 days')
spectrum_7d = redback.transient.Spectrum(wave, flux_7d, err_7d, time='7 days')

# Create likelihoods with epoch-specific parameters
# Use lambda functions or wrappers to map parameters correctly
spec_1d_likelihood = redback.likelihoods.GaussianLikelihood(
    x=spectrum_1d.angstroms,
    y=spectrum_1d.flux_density,
    sigma=spectrum_1d.flux_density_err,
    function=lambda wave, temp_1d, r_phot_1d, **kw:
        spectral_models.blackbody_spectrum_at_z(
            wave, temp=temp_1d, rph=r_phot_1d, **kw
        )
)
# Similar for 3d and 7d...

# Set up priors for all epochs
priors = bilby.core.prior.PriorDict()
priors['mej'] = ...      # Shared across all
priors['temp_1d'] = ...  # Epoch-specific
priors['temp_3d'] = ...  # Epoch-specific
priors['temp_7d'] = ...  # Epoch-specific
priors['r_phot_1d'] = ...
priors['r_phot_3d'] = ...
priors['r_phot_7d'] = ...

# Combine all likelihoods
joint_likelihood = MultiMessengerLikelihood(
    phot_likelihood,
    spec_1d_likelihood,
    spec_3d_likelihood,
    spec_7d_likelihood
)

Best Practices

Time synchronization
- Ensure spectrum epochs align with photometry time grid
- Use consistent time reference (e.g., explosion time, discovery)
- Account for any time delays between different observations
Parameter consistency
- Share physical parameters (mass, velocity, composition)
- Allow epoch-dependent parameters to vary (temperature, radius)
- For temperature evolution, consider parametric models: \(T(t) \propto t^{-\alpha}\)
Model selection
- Start with simple models (blackbody spectrum + analytic lightcurve)
- Add complexity as justified by data quality
- For detailed analysis, use consistent radiative transfer for both photometry and spectra
Wavelength coverage
- Verify spectrum wavelength range overlaps with photometric bands
- Consider filter transmission functions when comparing
- Account for extinction/reddening consistently
Systematic uncertainties
- Include flux calibration uncertainties
- Consider host galaxy contamination
- Account for atmospheric/instrumental effects

Example Workflow

See examples/joint_spectrum_photometry_example.py for a complete worked example demonstrating:

Simulating multi-band photometry and spectrum
Building appropriate likelihoods
Setting up shared and epoch-specific priors
Running individual and joint fits for comparison
Handling multiple spectra at different epochs

Comparison with Individual Fits

Joint fitting typically provides:

Tighter constraints on shared parameters due to complementary information
Breaking degeneracies (e.g., ejecta mass vs. velocity)
Consistency checks - if joint fit has much lower evidence than individual fits, models may be inconsistent
Better predictions for unobserved wavelengths/epochs

Compare results:

# Run photometry-only fit
phot_result = bilby.run_sampler(phot_likelihood, priors=phot_priors, ...)

# Run spectrum-only fit
spec_result = bilby.run_sampler(spec_likelihood, priors=spec_priors, ...)

# Run joint fit
joint_result = bilby.run_sampler(joint_likelihood, priors=joint_priors, ...)

# Compare evidence
print(f"Photometry ln(Z): {phot_result.log_evidence:.2f}")
print(f"Spectrum ln(Z): {spec_result.log_evidence:.2f}")
print(f"Joint ln(Z): {joint_result.log_evidence:.2f}")

# Compare posteriors
import corner
fig = corner.corner(joint_result.posterior, color='blue', labels=param_labels)
corner.corner(phot_result.posterior, fig=fig, color='red')

Joint Galaxy and Transient Spectrum Analysis

When observing transients embedded in bright host galaxies, the observed spectrum contains contributions from both the transient and the underlying galaxy. Properly accounting for this contamination is critical for accurate parameter inference.

Why This Matters

Host galaxy contamination can significantly bias transient parameters:

Continuum shape: Galaxy stellar continuum affects transient temperature/SED fits
Line features: Galaxy emission/absorption lines can mimic or hide transient features
Flux normalization: Incorrect flux leads to wrong luminosity, radius estimates
Common scenarios: Supernovae (especially SNe Ia), TDEs, transients in galaxy centers

Combined Model Approach

The most straightforward approach is to model the observed spectrum as the sum of galaxy and transient components:

def combined_galaxy_transient_model(wavelength, redshift,
                                   galaxy_temperature, galaxy_luminosity,
                                   transient_temperature, transient_r_phot,
                                   **kwargs):
    '''Combined model: observed spectrum = galaxy + transient'''

    # Galaxy component (stellar continuum)
    galaxy_flux = galaxy_spectrum_model(
        wavelength, redshift, galaxy_temperature, galaxy_luminosity
    )

    # Transient component (e.g., SN photosphere)
    transient_flux = transient_spectrum_model(
        wavelength, redshift, transient_temperature, transient_r_phot
    )

    # Observed spectrum is the sum
    return galaxy_flux + transient_flux

# Create likelihood with combined model
combined_likelihood = redback.likelihoods.GaussianLikelihood(
    x=wavelengths,
    y=observed_flux,
    sigma=flux_errors,
    function=combined_galaxy_transient_model,
    kwargs={}
)

Setting Up Priors

Define priors for both galaxy and transient parameters:

priors = bilby.core.prior.PriorDict()

# Shared parameter: redshift (typically well-known from galaxy)
priors['redshift'] = bilby.core.prior.Gaussian(0.05, 0.001)

# Galaxy parameters
priors['galaxy_temperature'] = bilby.core.prior.Uniform(4000, 7000)  # K
priors['galaxy_luminosity'] = bilby.core.prior.Uniform(0.5, 10.0)

# Transient parameters
priors['transient_temperature'] = bilby.core.prior.Uniform(5000, 15000)  # K
priors['transient_r_phot'] = bilby.core.prior.LogUniform(1e14, 1e16)  # cm

The shared redshift provides a strong constraint linking both components.

Running the Analysis

result = bilby.run_sampler(
    likelihood=combined_likelihood,
    priors=priors,
    sampler='dynesty',
    nlive=1000,
    outdir='./results_galaxy_transient',
    label='joint_galaxy_transient'
)

# Extract best-fit components
best_params = result.posterior.iloc[result.posterior['log_likelihood'].idxmax()]

galaxy_component = galaxy_spectrum_model(wavelengths, **best_params)
transient_component = transient_spectrum_model(wavelengths, **best_params)

Comparison with Transient-Only Fit

Ignoring the galaxy leads to biased parameters. Compare fits with/without galaxy:

# Biased fit (transient only)
transient_only_likelihood = redback.likelihoods.GaussianLikelihood(
    x=wavelengths, y=observed_flux, sigma=flux_errors,
    function=transient_spectrum_model
)
result_biased = bilby.run_sampler(transient_only_likelihood, transient_priors, ...)

# Correct fit (galaxy + transient)
result_correct = bilby.run_sampler(combined_likelihood, full_priors, ...)

# Evidence comparison
print(f"Transient-only ln(Z): {result_biased.log_evidence:.2f}")
print(f"Joint model ln(Z): {result_correct.log_evidence:.2f}")
print(f"Bayes factor: {result_correct.log_evidence - result_biased.log_evidence:.2f}")

The joint model should show much higher evidence if galaxy contamination is significant.

Adding Galaxy Emission Lines

For more realistic galaxy modeling, include emission lines:

def galaxy_with_lines_model(wavelength, redshift,
                            galaxy_temperature, galaxy_luminosity,
                            h_alpha_flux, h_beta_flux, oiii_flux,
                            line_width, **kwargs):
    '''Galaxy model with continuum + emission lines'''

    # Stellar continuum
    continuum = galaxy_continuum_model(...)

    # Common emission lines
    h_alpha = gaussian_line(wavelength, 6563 * (1+redshift), h_alpha_flux, line_width)
    h_beta = gaussian_line(wavelength, 4861 * (1+redshift), h_beta_flux, line_width)
    oiii = gaussian_line(wavelength, 5007 * (1+redshift), oiii_flux, line_width)

    return continuum + h_alpha + h_beta + oiii

def gaussian_line(wavelength, center, amplitude, width):
    '''Gaussian emission line profile'''
    return amplitude * np.exp(-0.5 * ((wavelength - center) / width)**2)

Add priors for emission line parameters:

priors['h_alpha_flux'] = bilby.core.prior.LogUniform(1e-17, 1e-15)
priors['h_beta_flux'] = bilby.core.prior.LogUniform(1e-18, 1e-16)
priors['oiii_flux'] = bilby.core.prior.LogUniform(1e-18, 1e-16)
priors['line_width'] = bilby.core.prior.Uniform(1, 5)  # Angstroms

Best Practices

Pre-explosion data
- If available, use pre-explosion galaxy spectrum to constrain galaxy parameters
- Can fix galaxy parameters or use as informative priors
- Reduces degeneracies between galaxy and transient
Model selection
- Start simple: blackbody or power-law continuum for both components
- Add complexity as justified: emission lines, absorption features
- Use stellar population synthesis for galaxy (e.g., FSPS, BC03)
- Use radiative transfer for transient (e.g., TARDIS, CMFGEN)
When to use joint fitting
- Transient is bright relative to galaxy (SNR > 10)
- No pre-explosion spectrum available for template subtraction
- Transient and galaxy spectra overlap significantly
- Need physical model for both components
Alternative approaches
- Template subtraction: If pre-explosion spectrum exists, subtract scaled template
- Spectral decomposition: Use tools like STARLIGHT or pPXF first
- Spatial separation: Extract spatially separated spectra if PSF resolution allows
Systematic checks
- Compare to pre-explosion imaging/spectroscopy
- Verify galaxy parameters are consistent across multiple epochs
- Check transient parameters match photometric evolution
- Account for extinction (Galactic + host)
Validation
- Plot decomposition showing galaxy, transient, and total components
- Check residuals for systematic structure
- Compare posterior distributions with/without galaxy model
- Use simulations to verify parameter recovery

Example Workflow

See examples/joint_galaxy_transient_spectrum_example.py for a complete demonstration including:

Simulating composite spectrum (galaxy + transient)
Setting up combined model
Comparing fits with/without galaxy model
Visualizing spectral decomposition
Showing parameter biases when galaxy is ignored
Advanced topics: emission lines, systematic uncertainties

Common Applications

Supernova spectroscopy: Type Ia SNe in bright host galaxies require careful galaxy subtraction. Joint fitting allows proper separation and uncertainty propagation.
Tidal Disruption Events (TDEs): TDEs occur in galaxy centers, often with strong AGN or stellar contamination. Joint modeling separates transient from persistent emission.
Kilonovae: While typically in fainter hosts, nearby events may require galaxy modeling, especially in red bands where older stellar populations contribute.
AGN outbursts: Separating variable AGN component from host galaxy requires joint modeling of both persistent and transient features.

Advanced: Including Gravitational Wave Data

For GW+EM analysis, construct a gravitational wave likelihood using bilby.gw and pass it to the MultiMessengerTransient:

import bilby.gw

# Set up GW analysis (following bilby.gw workflow)
duration = 32
sampling_frequency = 2048

waveform_generator = bilby.gw.WaveformGenerator(
    duration=duration,
    sampling_frequency=sampling_frequency,
    frequency_domain_source_model=bilby.gw.source.lal_binary_neutron_star,
    waveform_arguments={'waveform_approximant': 'IMRPhenomPv2_NRTidal'}
)

# Set up interferometers
interferometers = bilby.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
interferometers.set_strain_data_from_power_spectral_densities(
    sampling_frequency=sampling_frequency,
    duration=duration
)

# Create GW likelihood
gw_likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
    interferometers=interferometers,
    waveform_generator=waveform_generator
)

# Create multi-messenger object with GW data
mm_transient = MultiMessengerTransient(
    optical_transient=optical_transient,
    xray_transient=xray_transient,
    gw_likelihood=gw_likelihood,
    name='GW170817'
)

# Define priors including GW parameters
priors = bilby.core.prior.PriorDict()

# GW parameters
priors['chirp_mass'] = bilby.core.prior.Gaussian(1.2, 0.1)
priors['mass_ratio'] = bilby.core.prior.Uniform(0.5, 1.0)
priors['luminosity_distance'] = bilby.gw.prior.UniformSourceFrame(10, 250)
priors['theta_jn'] = bilby.core.prior.Uniform(0, 3.14159)  # GW viewing angle
# ... other GW parameters ...

# EM parameters
priors['mej'] = bilby.core.prior.Uniform(0.01, 0.1)
# ... other EM parameters ...

# Run joint GW+EM analysis
result = mm_transient.fit_joint(
    models={'optical': kilonova_model, 'xray': afterglow_model},
    priors=priors,
    shared_params=['luminosity_distance', 'theta_jn'],  # Link GW and EM
    nlive=2000
)

Note that theta_jn in gravitational wave analysis often corresponds to the viewing angle in electromagnetic models, though the exact relationship depends on the jet geometry and assumptions.

Utility Functions

create_joint_prior

The create_joint_prior() utility helps construct prior dictionaries for joint analysis:

from redback.multimessenger import create_joint_prior

# Define individual priors
optical_priors = bilby.core.prior.PriorDict({
    'viewing_angle': bilby.core.prior.Uniform(0, 1.57),
    'mej': bilby.core.prior.Uniform(0.01, 0.1)
})

xray_priors = bilby.core.prior.PriorDict({
    'viewing_angle': bilby.core.prior.Uniform(0, 1.57),
    'logn0': bilby.core.prior.Uniform(-3, 2)
})

# Create joint prior with shared viewing_angle
joint_prior = create_joint_prior(
    individual_priors={'optical': optical_priors, 'xray': xray_priors},
    shared_params=['viewing_angle']
)

# Result:
# joint_prior['viewing_angle']  -> shared prior
# joint_prior['mej']            -> optical-specific, matching the model parameter name
# joint_prior['logn0']          -> X-ray-specific, matching the model parameter name

Non-shared parameters keep their likelihood parameter names. If the same non-shared name appears in multiple messenger priors, create_joint_prior() raises a ValueError; use distinct model parameter names or a small model wrapper in that case. If a requested shared parameter is not present in any individual prior and is not supplied through shared_param_priors, create_joint_prior() also raises a ValueError.

Adding/Removing Messengers

You can dynamically add or remove messengers:

# Add a new messenger
mm_transient.add_messenger('uv', transient=uv_transient)

# Or add with a pre-constructed likelihood
mm_transient.add_messenger('neutrino', likelihood=neutrino_likelihood)

# Remove a messenger
mm_transient.remove_messenger('radio')

Examples

Complete worked examples are available in examples/:

multimessenger_example.py
- Simulates optical kilonova + X-ray afterglow + radio afterglow
- Demonstrates individual vs. joint fitting
- Shows how shared parameters improve constraints
joint_spectrum_photometry_example.py
- Joint photometry + spectroscopy analysis
- Simulates multi-band lightcurves and spectrum at specific epoch
- Shows how to use custom likelihoods for different data types
- Demonstrates handling multiple spectra at different epochs
joint_galaxy_transient_spectrum_example.py
- Joint galaxy + transient spectrum analysis
- Shows how to separate host galaxy contamination from transient
- Demonstrates combined model approach (galaxy + transient)
- Compares biased (transient-only) vs. unbiased (joint) fits
- Includes emission line modeling
joint_grb_gw_example.py
- Joint GW + GRB afterglow analysis
- Shows integration with bilby.gw
- Demonstrates parameter linking between GW and EM

Best Practices

Start with individual fits

Fit each messenger independently first to:
- Verify your models are appropriate
- Check for any data quality issues
- Understand individual parameter constraints
- Use as comparison for joint analysis
Choose shared parameters carefully

Only share parameters that are physically the same across messengers:
- Viewing angle (if jet is aligned with binary orbital axis)
- Distance/redshift (always shared for same source)
- Event time (with appropriate time delays)
Use informative priors

For shared parameters, use priors that reflect physical constraints:
- Distance: use bilby.gw.prior.UniformSourceFrame for GW sources
- Viewing angle: consider constraints from non-detection in certain bands
Check model compatibility

Ensure your models are compatible in their parameterization:
- Same definition of viewing angle
- Consistent reference time
- Same distance definition (luminosity vs. comoving)
Validate with simulations

Test your analysis pipeline on simulated data:
- Inject known parameters
- Verify recovery in joint analysis
- Check that shared parameters are correctly constrained
Compare evidence

Use Bayes factors to compare:
- Joint model vs. independent models
- Different choices of shared parameters
- Different physical models

Common Pitfalls

Inconsistent parameter definitions

Different models may use different conventions (e.g., viewing angle from jet axis vs. orbital axis). Create wrapper functions to translate between conventions.
Time reference mismatch

Ensure all data use consistent time references (e.g., time of merger, trigger time).
Correlated systematic errors

Joint analysis assumes independent data. Be cautious if systematic errors are correlated across messengers.
Over-constraining with incompatible models

If your models are inconsistent or incomplete, joint analysis may give misleading results.

References

For more information on multi-messenger analysis:

Abbott et al. 2017 (GW170817): ApJL 848, L12
Coughlin et al. 2018 (Multi-messenger Bayesian PE): MNRAS 480, 3871
See also examples/joint_grb_gw_example.py for GW+EM implementation details
See also examples/analyse_spectrums.ipynb for basic spectrum fitting workflow

Multi-Messenger Analysis

Overview

The MultiMessengerTransient class

Basic Usage

Supported Messengers

Joint Analysis

Upper Limits

Shared Parameters

Individual Fits for Comparison

Joint Spectrum and Photometry Analysis

Basic Approach

Setting Up Priors

Running Joint Analysis

Multiple Spectra at Different Epochs

Best Practices

Example Workflow

Comparison with Individual Fits

Joint Galaxy and Transient Spectrum Analysis

Why This Matters

Combined Model Approach

Setting Up Priors

Running the Analysis

Comparison with Transient-Only Fit

Adding Galaxy Emission Lines

Best Practices

Example Workflow

Common Applications

Advanced: Including Gravitational Wave Data

Utility Functions

create_joint_prior

Adding/Removing Messengers

Examples

Best Practices

Common Pitfalls

References

The `MultiMessengerTransient` class