Welcome to concord’s documentation!

concord implements some commonly-used approaches for determining system parameters for neutron stars from thermonuclear bursts, and fully accounting for astrophysical uncertainties. The code includes the capability for quantitative comparison of observed burst properties with the predictions of numerical models, and the tools presented here are intended to make comprehensive model-observation comparisons straightforward.

This code is under active development, but the v1.0.0 release is associated with a companion paper accepted by Astrophysical Journal Supplements (see Galloway et al. 2022, also available at arXiv:2210.03598). A preprint of the paper is available in the doc subdirectory of the repository

Getting started

Clone the repository and install using pip:

python3 -m pip install

You can then import the repository and use the functions. Here’s a very simple example, to find the peak luminosity of a burst from 4U 0513+40 measured by RXTE, as part of the MINBAR sample. The first part calculates the isotropic luminosity, neglecting the uncertainty in both the peak flux and the distance:

import concord as cd
import astropy.units as u
F_pk, e_F_pk = 21.72, 0.6 # 1E-9 erg/cm^2/s bolometric; MINBAR #3443
d = (10.32, 0.24, 0.20) # asymmetric errors from Watkins et al. 2015
l_iso = cd.luminosity( F_pk, dist=d[0], isotropic=True )
print (l_iso)
# 2.767771097997098e+38 erg / s

The second part takes into account both the uncertainties in the peak flux and distance (including the asymmetric errors), and also includes the model-predicted effect of the high system inclination (>80 degrees):

l_asym = cd.luminosity( (F_pk, e_F_pk), dist=d, burst=True, imin=80, imax=90, fulldist=True)
lc = l_asym['lum'].pdf_percentiles([50, 50 - cd.CONF / 2, 50 + cd.CONF / 2])
l_unit = 1e38*u.erg/u.s
print ('''\nIsotropic luminosity is {:.2f}e38 erg/s
  Taking into account anisotropy, ({:.2f}-{:.2f}+{:.2f})e38 erg/s'''.format(l_iso/l_unit, lc[0]/l_unit, (lc[0]-lc[1])/l_unit, (lc[2]-lc[0])/l_unit))
# Isotropic luminosity is 2.77e38 erg/s
#   Taking into account anisotropy, (4.85-0.43+0.51)e38 erg/s

With the fulldist=True option, the function returns a dictionary with the Monte-Carlo generated distribution of the result (key lum) and all the intermediate quantities.

Check the Inferring burster properties notebook, which includes the examples used in the paper, for additional demonstrations of usage.

Obtaining burst data

You’ll also need some burst measurements or data; at their simplest, most of the functions operate on measurable quantities (with uncertainties), like the peak burst flux, burst fluence, persistent flux etc.

You can access data from over 7000 bursts from 85 sources via MINBAR, by also installing the Python repository from https://github.com/outs1der/minbar.

You can also use your own time-resolved spectroscopic data to generate an concord.burstclass.ObservedBurst object, and then run the functions on that. You might need to define routines to read in your own data, depending on the format.

You can download some reference bursts from http://burst.sci.monash.edu/reference, and put the data (including the .dat files and the table2.tex file) in the concord/data directory.

You can use the example concord.burstclass.KeplerBurst class to read in some model bursts from http://burst.sci.monash.edu/kepler.

Alternatively you can use that class as an example to write your own appropriate for your model results.

How the functions work

There are three principal components; the functions in utils.py, the anisotropy treatment in diskmodel.py and the observed and model burst classes in burstclass.py.

The functions in utils.py provide the basic functionality as described in Functions. In order to treat the uncertainties in the observed quantities, we adopt a Monte-Carlo (MC) approach via the astropy Distribution package. Input values (including measured quantities) can be provided in four different ways:

  • scalar

  • value with symmetric error, represented as a tuple (value, error)

  • value with asymmetric error, as a 3-element tuple (value, lo_error, hi_error)

  • an arbitrary array.

Quantities with error estimates will be converted to a Distribution object of user-defined size, and with a symmetric or asymmetric normal distribution. In this way we can represent a wide range of probability distribution functions (PDFs) for the input parameters. These Distribution objects can then be used in calculations as for scalars, hence providing uncertainty propagation at the cost of additional computation.

Input values can be provided with units, or if units are absent, will be assumed to have the standard units for MINBAR quantities; for flux, \(10^{-9}\ \rm{erg\,cm^{-2}\,s^{-1}}\), burst fluence \(10^{-6}\ \rm{erg\,cm^{-2}}\), recurrence time hr, and so on.

All functions can operate assuming isotropic distributions for the burst or persistent emission (isotropic=True), but by default will include the possible effects of anisotropic emission via the diskmodel.py routine. This routine incorporates the modelling of He & Keek (2016, ApJ 819, #47), via ASCII tables provided with the code.

There are three options for the treatment of emission anisotropy:

  • isotropic=True - assumes both the burst and persistent emission is isotropic (\(\xi_b=\xi_p=1\))

  • isotropic=False - incorporates anisotropy, by drawing a set of inclination values and calculating the H-fraction for each value. Inclinations are uniform on the sphere (up to a maximum value of \(i=72^\circ\), based on the absence of dips) and the value returned is the mean and 1-sigma error, or the complete distributions (with the fulldist option)

  • inclination=<i> - calculate for a specific value of inclination

Where observed or model burst lightcurves are available, the classes in burstclass.py provide a way to represent those observations and perform various standard analyses, including observation-model comparisons. The ObservedBurst class can be instantiated from an ASCII file giving the burst flux as a function of time, or in a number of other ways.

An example KeplerBurst model burst class is provided, which offers a number of ways of reading in Kepler burst runs, and which can be adapted to outputs of different codes.

The reference bursts provided by Galloway et al. (2017, PASA 34, e109) can be read in directly, provided the data is included in the data directory.

Functions

concord.utils.L_Edd(F_pk, dist=<Quantity 8. kpc>, nsamp=1000, isotropic=False, inclination=None, imin=0.0, imax=75.0, dip=False, conf=68.0, fulldist=False)

This routine estimates the Eddington luminosity of the source, based on the measured peak flux of a radius-expansion burst and the distance (or some distribution thereof. The default isotropic=False option also takes into account the likely effect of anisotropic burst emission, based on the models provided by concord

This routine has been supplanted by concord.utils.luminosity(), which it now calls; the only thing not incorporated into the new routine is the “simple” estimate of the isotropic luminosity error,

L_Edd_iso_err = L_Edd_iso * np.sqrt((F_pk_err / F_pk) ** 2 + (2. * dist_err / dist) ** 2)

Burst emission is assumed (burst=True); apart from the F_pk value, parameters are as for concord.utils.luminosity()

Parameters:

F_pk – peak burst flux, MINBAR units assumed if not present

Returns:

dictionary including calculation results and assumed distributions (as for luminosity)

concord.utils.Q_nuc(Xbar, quadratic=True, old_relation=False, coeff=False)

This function implements the approximation to the nuclear energy generation rate Q_nuc given the mean hydrogen fraction in the fuel layer, Xbar, determined by Goodwin et al. 2019 (ApJ 870, #64). The value returned is in units of MeV/nucleon:

>>> q = Q_nuc(0.73)
5.35517578

There are three versions of the approximation, which can be selected by using the quadratic and old_relation flags; by default the most precise (and recent) version is used:

>>> q = Q_nuc(0.73, quadratic=False)
5.758715
>>> q = Q_nuc(0.73, old_relation=True)
4.52

If you want the coefficients, rather than the actual Q_nuc value, you can use the coeff flag (in which case the first argument is ignored):

>>> q_0, q_1 = Q_nuc(0.0,quadratic=False,coeff=True)
[1.3455, 6.0455]
Parameters:

  • Xbar – average hydrogen mass fraction of the burst column

  • quadratic – set to True to use the (more accurate) quadratic approximation

  • old_relation – set to True to use the earlier approximation

  • coeff – set to True to return the coefficients used, in which case Xbar is ignored

Returns:

Q_nuc value, in MeV/nucleon, or the coefficients if coeff=True

concord.utils.X_0(Xbar, zcno, tdel, opz=1.259, debug=False, old_relation=False)

Routine (extracted from concord.utils.hfrac()) to determine the accreted H fraction X_0 given the average H-fraction at ignition, Xbar, the CNO metallicity zcno, and the burst recurrence time, tdel.

The time to exhaust the accreted hydrogen is calculated according to the formula derived by Lampe et al. (2016, ApJ 819, #46); if the old_relation=True flag is set, the alternative prefactor quoted by Galloway et al. (2004 ApJ 601, 466) is used instead

Parameters:

  • Xbar – average hydrogen mass fraction of the burst column

  • zcno – CNO mass fraction in the burst fuel

  • tdel – burst recurrence time

  • opz – NS surface redshift, 1+z

  • old_relation – use the older expression for the time to exhaust the accreted hydrogen

  • debug – display debugging messages

Returns:

inferred hydrogen mass fraction X_0 of the accreted fuel

concord.utils.alpha(_tdel, _fluen, _fper, c_bol=None, nsamp=1000, conf=68.0, fulldist=False)

Routine to calculate alpha, the observed ratio of persistent to burst energy (integrated over the recurrence time, tdel) from the input measurables

Usage:

>>> alpha = cd.alpha(2.681, 0.381, 3.72, 1.45)
>>> print (alpha)
136.64233701

>>> alpha = cd.alpha((2.681, 0.007), (0.381, 0.003), (3.72, 0.18), (1.45, 0.09))
>>> print (alpha)
[136.23913536   9.88340505  11.21502717]
Parameters:

  • tdel – burst recurrence time

  • fluen – burst fluence

  • fper – persistent flux

  • c_bol – bolometric correction on persistent flux

Returns:

alpha-values, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if fulldist=True) a dictionary with the distributions of the alpha value and the input or adopted distributions of the intermediate values

concord.utils.asym_norm(m, sigm=None, sigp=None, nsamp=1000, statistics='max', positive=False, model=1)

Draw samples from an asymmetric error distribution characterised by a mean and upper and lower 68% confidence intervals (sigp and sigm, respectively). Used primarily to generate distributions from quantities with asymmetric errors, by concord.utils.value_to_dist().

Follows the treatment of Barlow (2003)

With the positive flag set, it will continue drawing samples until all are > 0. Note that the resulting distribution may not quite have the required shape

Parameters:

  • m – mean (central) value for the distribution

  • sigm – lower 68th-percentile uncertainty

  • sigp – upper 68th-percentile uncertainty

  • nsamp – number of samples required

  • statistics – convention for input stats for asymmetric distributions

  • positive – ensure all the samples are positive (may affect the distribution)

  • model – implementation of the distribution function; only one option is currently implemented

Returns:

distribution with nsamp values having the desired shape

concord.utils.calc_mr(g, opz)

‘ Calculates neutron star mass and radius given a surface gravity and redshift. This routine is required when (for example) comparing a model computed at a fixed gravity, with observations requiring a particular value of redshift (to stretch the model lightcurve to match the observed one). The combination of surface gravity and redshift implies unique mass and radius, and the constraints on 1+z can be translated to constraints on M, R

Parameters:

  • g – surface gravity, in units equivalent to cm/s**2

  • opz – surface redshift 1+z

Returns:

neutron star mass in g, radius in cm

concord.utils.check_M_R_opz(M, R, opz)

Utility routine to check the consistency of the mass, radius and redshift passed to various routines (e.g. mdot)

Parameters:

  • M – neutron star mass, or None

  • R – neutron star radius, or None

  • opz – 1+z where z is the surface gravitational redshift, or None

Returns:

boolean, True if the values are consistent, False otherwise

concord.utils.create_logger()

Create a logger instance where messages are sent. See https://docs.python.org/3/library/logging.html

Returns:

logger instance

concord.utils.decode_LaTeX(string, delim='pm')

This function converts a LaTeX numerical value (with error) to one or more floating point values. It is expected that the number is formatted as $3.1\pm1.2$, and the default delimiter assumes this scenario

Will not work on asymmetric errors, expressed (for example) as $1.4_{-0.1}^{+0.2}$

Example usage:

>>> cd.decode_LaTeX('$1.4\pm0.3')
(1.4, 0.3)
Parameters:

  • string – string with LaTeX numerical value and error

  • delim – separator of number and uncertainty; defaults to \pm

Returns:

tuple giving the value and error(s)

concord.utils.dist(_F_pk, nsamp=None, dip=False, empirical=False, X=0.0, M=<Quantity 1.4 solMass>, opz=1.259, T_e=0.0, isotropic=False, inclination=None, imin=0.0, imax=75.0, model='he16_a', conf=68.0, fulldist=False, plot=False)

This routine estimates the distance to the source, based on the measured peak flux of a radius-expansion burst. Two main options are available; empirical=False (default), in which case the Eddington luminosity is calculated as equation 7 from Galloway et al. (2008, ApJS 179, 360); or empirical=True, in which case the estimate of Kuulkers et al. 2003 (A&A 399, 663) is used. The default isotropic=False option also takes into account the likely effect of anisotropic burst emission, based on the models provided by concord

Example usage:

>>> cd.dist((30., 3.), isotropic=True, empirical=True)
<Quantity [10.27267949,  0.50629732,  0.59722829] kpc>
>>> cd.dist((30., 2., 5.), isotropic=False, empirical=True)
<Quantity [10.75035916,  1.39721634,  1.461003  ] kpc>
>>> cd.dist((30., 3.), isotropic=False, empirical=True, fulldist=True, nsamp=100)
{'dist': <QuantityDistribution [11.12852694, 11.4435149 ,  9.79809775,  9.82146105,  9.34573615,
        10.4330567 , 10.99421133, 10.66816313,  8.83667458, 12.83249538,
        14.15062276, 13.26883293,  9.80702496,  8.72915705, 12.04546987,
            .
            .
            .
Parameters:

  • _F_pk – peak burst flux

  • nsamp – number of samples required for the distributions

  • dip – set to True if the source is a “dipper”

  • empirical – flag to use the empirical Eddington luminosity

  • X – hydrogen mass fraction in the atmosphere

  • M – neutron-star mass

  • opz – neutron-star surface redshift

  • T_e – temperature factor used in the theoretical L_Edd expression

  • isotropic – set to True to assume isotropic emission

  • inclination – inclination value (or distribution)

  • imin – minimum of allowed inclination range

  • imax – maximum of allowed inclination range

  • conf – confidence interval for output limits

  • fulldist – set to True to output full distributions for all parameters

  • plot – set to True to generate a simple plot

Returns:

inferred distance values, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if fulldist=True) a dictionary with the distributions of the distance and the input or adopted distributions of the other parameters

concord.utils.g(M, R, Newt=False, units='cm/s^2')

This function calculates the surface gravity given a mass M and radius R, using the Newtonian expression if the flag Newt is set to True

The result is returned in units of cm/s^2 by default

Parameters:

  • M – neutron-star mass

  • R – neutron-star radius

  • Newt – default False returns the GR value, set to True to calculate the Newtonian value

  • units – the required units

Returns:

the surface gravity, in the required units

concord.utils.hfrac(_tdel, _alpha=None, fper=None, fluen=None, c_bol=None, zcno=0.02, opz=1.259, old_relation=False, isotropic=False, inclination=None, imin=0.0, imax=75.0, model='he16_a', conf=68.0, fulldist=False, nsamp=None, debug=False)

Estimates the H-fraction at burst ignition, based on the burst properties In the absence of the alpha-value(s), you need to supply the persistent flux and burst fluence, along with the recurrence time, so that alpha can be calculated

There is also a mode for calculating a single value of the parameters, based on a single value of alpha, tdel and the inclination; this mode is used by the function itself, within a loop over the inclination values.

Example usage:

>>> import astropy.units as u
>>> cd.hfrac(2.5*u.hr, 140., inclination=30.)
** WARNING ** assuming inclination in degrees
(<Quantity 0.19564785>, <Quantity 0.26656582>, 30.0)

>>> cd.hfrac((2.681, 0.007), fluen=(0.381, 0.003),
        fper=(3.72, 0.18), c_bol=(1.45, 0.09),nsamp=100,fulldist=True)
{'xbar': NdarrayDistribution([ 0.15246606,  0.20964612,  0.14169592,  0.11998812,  0.14293726,
    0.13757633, -0.0850957 ,  0.27077167,  0.01620781,  0.13673547,
   -0.00821256,  0.10193499, -0.08151998, -0.0370918 ,  0.26400253,
    0.1813988 , -0.04371571,  0.14432454,  0.28422351, -0.04962202,
            .
            .
            .
 'X_0': NdarrayDistribution([ 0.22852126,  0.28567535,  0.21799341,  0.19609461,  0.21895133,
            .
            .
            .
 'i': <QuantityDistribution [45.58518673, 15.61321579, 36.30632308, 48.18062101, 45.46515182,
            .
            .
            .
Parameters:

  • tdel – burst recurrence time

  • alpha – burst alpha, which (if not supplied) can be calculated instead from the fluence, fper, and c_bol

  • fper – persistent flux

  • fluen – burst fluence

  • c_bol – bolometric correction on persistent flux

  • opz – neutron-star surface redshift

  • zcno – CNO mass fraction in the burst fuel

  • old_relation – flag to use the old (incorrect) relation for Q_nuc

  • isotropic – set to True to assume isotropic emission

  • inclination – inclination value (or distribution)

  • imin – minimum of allowed inclination range

  • imax – maximum of allowed inclination range

  • conf – confidence interval for output limits

  • model – model string of He & Keek (2016), defaults to “A”

  • conf – confidence % level for limits, ignored if fulldist=True

  • fulldist – set to True to return the distributions of each parameter

  • debug – display debugging messages

Returns:

a tuple of values of Xbar, X_0 and inclination; or a dictionary including distributions of each of the values, as well as the distributions adopted for the observables

concord.utils.homogenize_params(theta, nsamp=None)

This method is used by the MC routines, to ensure consistency of the parameter set. We want to make sure that each parameter (if a distribution) has the same length as any other parameter (unless it’s a scalar). Each parameter is labeled and provided with a value (and possibly also error), and units. Inclination is also provided with the isotropy flag, and the angle limits.

If none of the input parameters are already distributions, the provided number of samples nsamp will be used as the dimensions of the output distributions (this may be redundant, since we now also do the inclination distributions here)

Usage:

(par1, par2, … nsamp ) = homogenize_params( dictionary_of_input_params_and_units, nsamp )

Example usage:

>>> F_pk, _nsamp = homogenize_params( {'fpeak': ((33., 2.), cd.MINBAR_FLUX_UNIT)}, 100)

>>> F_pk, _inclination, _nsamp = cd.homogenize_params( {'fpeak': ((33., 2.), cd.MINBAR_FLUX_UNIT),
                               'incl': (None, u.deg, isotropic, 0.0, 75.)}, 100)
Parameters:

  • theta – dictionary with parameters and units

  • nsamp – desired size of distribution objects, if not already set by one or more of the parameters

Returns:

tuple with all the parameters in the same order they’re passed, plus the actual number of samples per parameter

concord.utils.intvl_to_errors(perc)

This function converts a confidence interval defined by a 3-element array to the more useful version with errors .

Parameters:

perc – the array of percentiles (central, lower limit, upper limit)

Returns:

3-element array with (central value, lower uncertainty, upper uncertainty)

concord.utils.iso_dist(nsamp=1000, imin=0.0, imax=75.0)

Routine to generate an isotropic distribution of inclinations (angle from the system rotation axis to the line of sight) from imin up to some maximum value, defaulting to 75 degrees. This value corresponds to the likely maximum possible for a non-dipping source.

Parameters:

  • nsamp – number of samples required

  • imin – lower limit of distribution in degrees (default is zero)

  • imax – upper limit of distribution in degrees (default is IMAX_NDIP)

Returns:

distribution of inclinations between imin, imax

concord.utils.len_dist(d)

Utility routine to replace the len function for distributions, and also return sensible values for scalars rather than throwing an error

Parameters:

d – object (scalar, array or distribution)

Returns:

length of the object (0 if empty/None)

concord.utils.lum_to_flux(_lum, dist=None, c_bol=None, nsamp=None, burst=True, dip=False, isotropic=False, inclination=None, imin=0.0, imax=75.0, model='he16_a', conf=68.0, fulldist=False, plot=False)

This routine converts luminosity to flux, based on the provided distance Basically the inverse of the concord.utils.luminosity() function

Parameters:

  • lum – luminosity to convert to flux

  • dist – distance to the source

  • c_bol – bolometric correction to apply

  • nsamp – number of samples

  • burst – set to True for burst emission, to select the model anisotropy for bursts

  • dip – set to True for a dipping source (deprecated)

  • isotropic – set to True if isotropic value required

  • inclination – system inclination or distribution thereof

  • imin – minimum inclination for generated distribution

  • imax – maximum inclination for generated distribution

  • model – model string of He & Keek (2016), defaults to “A”

  • conf – confidence % level for limits, ignored if fulldist=True

  • fulldist – set to True to return the distributions of each parameter

  • plot – plots the resulting distributions

Returns:

inferred flux, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if fulldist=True) a dictionary with the distributions of the flux and the input or adopted distributions of the intermediate values

concord.utils.luminosity(_F_X, dist=None, c_bol=None, nsamp=None, burst=True, dip=False, isotropic=False, inclination=None, imin=0.0, imax=75.0, model='he16_a', conf=68.0, fulldist=False, plot=False)

Calculate the inferred luminosity given a measured X-ray flux and distance The flux, error and distance can be single values, or arrays (in which case it’s assumed that they’re distributions)

This is a more sophisticated version of the L_Edd routine (forgot I had that)

Example usage:

Calculate the isotropic luminosity corresponding to a flux of 3e-9 erg/cm^2/s at 7.3 kpc

>>> import concord as cd
cd.luminosity(3e-9,7.3,isotropic=True)

Calculate the range of luminosities corresponding to a persistent flux of 3e-9 erg/cm^2/s at 7.3 kpc, and assuming isotropic inclination distribution (i < 75 deg)

>>> cd.luminosity(3e-9,7.3,burst=False)

Calculate the range of luminosities corresponding to a burst flux of 3e-8, with uncertainty 1e-9, and an inclination of 45-60 degrees, plot (and return) the resulting distribution

>>> cd.luminosity((3e-8,1e-9),7.3,imin=45,imax=60,plot=True,fulldist=True)
Parameters:

  • F_X – X-ray flux, MINBAR units assumed if not present

  • dist – distance to the source

  • c_bol – bolometric correction to apply

  • nsamp – number of samples

  • burst – set to True for burst emission, to select the model anisotropy for bursts

  • dip – set to True for a dipping source (deprecated)

  • isotropic – set to True if isotropic value required

  • inclination – system inclination or distribution thereof

  • imin – minimum inclination for generated distribution

  • imax – maximum inclination for generated distribution

  • model – model string of He & Keek (2016), defaults to “A”

  • conf – confidence % level for limits, ignored if fulldist=True

  • fulldist – set to True to return the distributions of each parameter

  • plot – plots the resulting distributions

Returns:

luminosity, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if fulldist=True) a dictionary with the distributions of the luminosity and the input or adopted distributions of the intermediate values

concord.utils.mdot(_F_per, _dist, c_bol=None, M=None, R=None, opz=None, isotropic=False, inclination=None, imin=0.0, imax=75.0, dip=False, model='he16_a', nsamp=None, conf=68.0, fulldist=False)

Routine to estimate the mdot given a (bolometric) persistent flux, distance, and inclination. This calculation was adapted initially from equation 2 of Galloway et al. (2008, ApJS 179, 360), and uses an approximation to \(Q_\rm{grav} = c^2 z/(1+z) \approx GM_{\rm NS}/R_{\rm NS}\) which is good to about 10% for a typical neutron star

Usage:

>>> import concord as cd
# calculate the mdot corresponding to a flux of 1e-9 at 10kpc, for a 10km
# 1.4 M_sun neutron star (this is the prefactor for equation 2 in Galloway
# et al. 2008)
>>> import astropy.constants as c
>>> cd.mdot(1., 10., M=1.4*c.M_sun, R=10.*u.km, isotropic=True)
WARNING:homogenize_params:no bolometric correction applied
<Quantity 6691.30392224 g / (cm2 s)>
Parameters:

  • F_per – persistent flux, MINBAR units assumed if not present

  • dist – distance to the source

  • c_bol – bolometric correction to apply

  • M – neutron-star mass

  • R – neutron-star radius

  • opz – 1+z where z is the surface gravitational redshift, or None

  • isotropic – set to True if isotropic value required

  • inclination – system inclination or distribution thereof

  • imin – minimum inclination for generated distribution

  • imax – maximum inclination for generated distribution

  • dip – set to True for a dipping source (deprecated)

  • model – model string of He & Keek (2016), defaults to “A”

  • nsamp – number of samples to generate

  • conf – confidence % level for limits, ignored if fulldist=True

  • fulldist – set to True to return the distributions of each parameter

Returns:

mdot-values, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if fulldist=True) a dictionary with the distributions of the mdot value and the input or adopted distributions of the intermediate values

concord.utils.redshift(M, R)

This function calculates the gravitational redshift 1+z

Parameters:

  • M – neutron-star mass

  • R – neutron-star radius

Returns:

1+z

concord.utils.solve_radius(M, R_Newt, eta=1e-06)

This routine determines the GR radius given the NS mass and Newtonian radius, assuming the GR and Newtonian masses are identical

Solving is tricky so we just use an iterative approach

Parameters:

  • M – neutron-star mass

  • R_Newt – Newtonian neutron-star radius

  • eta – tolerance for convergence

Returns:

neutron-star radius R in GR

concord.utils.tdel_dist(nburst, exp, nsamp=1000)

Function to generate a synthetic distribution of recurrence times, given nburst observed events over a total exposure of exp. The resulting distribution may approximate the PDF of the underlying burst rate, but does assume that the bursts are independent, which is generally not the case

Parameters:

  • nburst – number of events detected

  • exp – total exposure time

  • nsamp – number of samples to generate

Returns:

distribution giving the PDF of the average recurrence time

concord.utils.value_to_dist(_num, nsamp=1000, unit=None, statistics='max', verbose=False)

This method converts a measurement to a distribution, to allow flexibility in how values are implemented in the various routines. Primarily used by concord.utils.homogenize_distributions().

Intended to be very flexible in terms of the way the units are provided, but unfortunately np.shape for a tuple with elements having units, does not work! So have to pass the units on the entire object

Also have an issue with combining distributions with and without units; “…the NdarrayDistribution will not combine with Quantity objects containing units” see https://docs.astropy.org/en/stable/uncertainty for more details

Parameters:

  • num – scalar/array to convert to distribution

  • statistics – convention for input stats for asymmetric distributions

Returns:

astropy distribution object

Example usage:

>>> y = cd.value_to_dist(3.) # scalars are kept as single values
>>> z = cd.value_to_dist((3.,0.1),nsamp=10) # generate 10 samples from a normal distribution around 3.0 with st. dev 0.1
>>> a = cd.value_to_dist((3.,0.5,0.1)) # generate default number of samples from an asymmetric Gaussian
>>> b = cd.value_to_dist(a.distribution) # convert an array to a distribution

with units:

>>> z = cd.value_to_dist((3.,0.1)*u.hr,nsamp=10) # generate 10 samples from a normal distribution around 3.0 with st. dev 0.1
concord.utils.yign(_E_b, dist=None, nsamp=None, R=<Quantity 11.2 km>, opz=1.259, Xbar=0.7, quadratic=False, old_relation=False, isotropic=False, inclination=None, imin=0.0, imax=75.0, dip=False, model='he16_a', conf=68.0, fulldist=False)

Calculate the burst column from the burst fluence, adapted initially from equation 4 of Galloway et al. (2008, ApJS 179, 360)

Parameters:

  • _E_b – burst fluence

  • _dist – distance ot bursting source

  • R – neutron star radius

  • opz – 1+z, surface redshift

  • Xbar – mean H-fraction at ignition

  • quadratic – flag to use the quadratic expression for Q_nuc

  • old_relation – flag to use the old (incorrect) relation for Q_nuc

  • isotropic – assume isotropic flux (or not)

  • inclination – inclination value (or array; in degrees)

  • imin – minimum inclination (degrees)

  • imax – maximum inclination (degrees)

  • dip – whether or not the source is a dipper

  • model – model for anisotropy calculation

  • nsamp – number of samples to generate

  • conf – confidence interval for uncertainties; or

  • fulldist – output the full distribution of calculated values

Returns:

inferred ignition column values, either a scalar (if all the inputs are also scalars); a central value and confidence range; or (if fulldist=True) a dictionary with the distributions of the ignition column and the input or adopted distributions of the other parameters

Classes

class concord.burstclass.KeplerBurst(filename=None, run_id=None, path=None, source='gs1826', grid_table=None, batch=None, run=None, **kwargs)

Example simulated burst class. Apart from the lightcurve (which is defined with time, lumin, and lumin_err columns), the additional (minimal) attributes requred are:

Filename:

source file name

Tdel, tdel_err:

recurrence time and error (hr)

Lacc:

accretion luminosity (in units of Mdot_Edd)

R_NS:

model-assumed radius of the neutron star (GR)

G:

surface gravity assumed for the run

An example call is as follows:

>>> c_loZ=KeplerBurst(filename='mean1.data',path='kepler', lAcc=0.1164,Z=0.005,H=0.7, tdel=4.06/opz,tdel_err=0.17/opz, g = 1.858e+14*u.cm/u.s**2, R_NS=11.2*u.km)

If you provide a batch and run value, like so:

>>> c = KeplerBurst(batch=2,run=9)

or (for the older results) a run_id value, like so:

>>> c = KeplerBurst(run_id='a028',path='../../burst/reference')

the appropriate model lightcurve will be read in, and the model parameters will be populated from the table

info()

Display information about the KeplerBurst

class concord.burstclass.Lightcurve(*args, **kwargs)

The fundamental burst object, which has attributes like time (array), flux/luminosity etc. and methods including plot. Minimal attributes are:

Time, dt:

array of times and time bin duration

Timepixr:

set to 0.0 (default) if time is the start of the bin, 0.5 for midpoint

Flux, flux_err:

flux and error (no units assumed)

Lumin, lumin_err:

luminosity and error

Normally only one of flux or luminosity would be supplied

Example: simulated burst (giving time and luminosity)

>>> Lightcurve(self, time=d['time']*u.s, lumin=d['luminosity']*u.erg/u.s, lumin_err=d['u_luminosity']*u.erg/u.s)

Example: observed burst giving flux

>>> Lightcurve(self, time=d['col1']*u.s, dt=d['col2']*u.s, flux=d['col3']*1e-9*u.erg/u.cm**2/u.s, flux_err=d['col4']*1e-9*u.erg/u.cm**2/u.s)
fluence(plot=False, warnings=True)

Calculate the fluence for a lightcurve, following the approach from get_burst_data, as used for MINBAR

This code relies on the presence of the dt_nogap and good attributes, which should be calculated when the Lightcurve object is defined

Parameters:

  • plot – set to True to display a diagnostic plot

  • warnings – set to False to suppress display of warnings

Returns:

the burst fluence and uncertainty

observe(param=[<Quantity 6.1 kpc>, <Quantity 60. deg>, 1.26, <Quantity -10. s>], obs=None, disc_model='he16_a', c_bol=1.0)

Convert a luminosity profile to a simulated observation, with plausible errors. See Keek & Heger (2011, ApJ 743, #189) for the background calculations

This method might make more sense as part of the :py:class:`KeplerBurst`_ class (or a parent SimulatedBurst class), but for now we include this method as part of the lightcurve class, so that future classes of model bursts can access it

Parameters:

  • param – tuple with the distance, inclination, redshift and time offset

  • obs – template observed burst to use for the time bins. If not supplied, uniform 0.25 s bins will be used, covering the burst

  • disc_model – choice of disc model for the anisotropy

  • c_bol – bolometric correction to adjust the burst fluxes

Returns:

concord.burstclass.ObservedBurst object

plot(yerror=True, obs_color='b', model_color='g', **kwargs)

Plot the lightcurve, accommodating both flux and luminosities

Parameters:

  • yerror – display flux/luminosity errors (default True)

  • obs_color – set the colour for an observed burst

  • model_color – set the colour for a simulated burst

Returns:

the plot

print()

Print the basic parameters of the lightcurve. Can be called in turn by the info commands for concord.burstclass.KeplerBurst, concord.burstclass.ObservedBurst to display the properties of the parent class… once those are written

Returns:

write(filename='lightcurve.csv', addhdr=None)

Write the lightcurve to a file

Parameters:

  • filename – filename to save the data to

  • addhdr – additional text to include as a header (default None)

class concord.burstclass.ObservedBurst(time, dt, flux, flux_err, **kwargs)

Observed burst class. Apart from the lightcurve (which is defined with time, flux, and flux_err columns), the additional attributes set are:

Filename:

source file name

Tdel, tdel_err:

recurrence time and error (hr)

Comments:

ASCII file header text

Table_file:

LaTeX file of table 2 from Galloway et al. (2017)

Table:

contents of table 2

Row:

entry in the table corresponding to this burst

Fper, fper_err:

persistent flux level (1e-9 erg/cm^2/s)

C_bol:

bolometric correction

Mdot:

accretion rate (total, not per unit area)

Fluen_table:

burst fluence

F_pk:

burst peak flux

Alpha:

alpha value

The units for selected parameters are all stored as astropy.units objects

compare(mburst, param=[<Quantity 6.1 kpc>, <Quantity 60. deg>, 1.26, <Quantity -10. s>], breakdown=False, plot=False, subplot=True, weights={'fluxwt': 1.0, 'tdelwt': 2500.0}, disc_model='he16_a', debug=False)

This is the key method for burst observation-model lightcurve comparisons, and generates a likelihood that can be used for MCMC (for example). It can be used to plot the observations with the models rescaled by the appropriate parameters.

weights give the relative weight to the recurrence time and persistent flux for the likelihood. Since you have many more points in the lightcurve, you may want to weight these greater than one so that the MCMC code will try to match those preferentially

Parameters:

  • mburst – model burst for comparison, for example an concord.burstclass.KeplerBurst object

  • param – tuple with the distance, inclination, redshift and time offset

  • breakdown – display the breakdown of the likelihood between the recurrence time, persistent flux and lightcurve

  • plot – set to True to show a plot of the comparison

  • subplot – include a subplot with the recurrence time comparison (default True)

  • weights – dict giving the relative weight to the recurrence time (key tdelwt) and persistent flux (fluxwt) for the likelihood.

  • disc_model – choice of disc model for the anisotropy

  • debug – display debugging information

classmethod from_file(filename, path='.', **kwargs)

Method to read in a file and populate the relevant arrays to create an concord.burstclass.ObservedBurst object This routine will read in any ascii lightcurve file which matches the format of the “reference” bursts, i.e. with columns:

'Time [s]' 'dt [s]' 'flux [10^-9 erg/cm^2/s]' 'flux error [10^-9 erg/cm^2/s]' 'blackbody temperature kT [keV]' 'kT error [keV]' 'blackbody normalisation K_bb [(km/d_10kpc)^2]' 'K_bb error [(km/d_10kpc)^2]' chi-sq
  -1.750  0.500   1.63    0.054   1.865  0.108   13.891   4.793  0.917
  -1.250  0.500   2.88    1.005   1.862  0.246   22.220   4.443  1.034
  -0.750  0.500   4.38    1.107   1.902  0.093   30.247   2.943  1.089
  -0.250  0.500   6.57    0.463   1.909  0.080   46.936   6.969  0.849
    .
    .
    .

Calling approach: >>> obs = ObservedBurst.from_file(filename) >>> obs = ObservedBurst.from_file(‘gs1826-24_3.530h.dat’,’example_data’)

info()

Display information about the ObservedBurst

Returns:

classmethod minbar(id, remote=False)

Method to generate an observed burst from MINBAR, and populate the relevant arrays. Requires the minbar repo to be present (obviously)

Calling approach:

>>> obs = ObservedBurst.minbar(2203)
Parameters:

  • id – the burst ID in MINBAR DR1

  • remote – if true, force download from the remote repo

Returns:

classmethod ref(source, dt)

Method to read in a reference burst and populate the relevant arrays to create an concord.burstclass.ObservedBurst

Calling approach:

>>> obs = ObservedBurst.ref('GS 1826-24', 3.5)
concord.burstclass.fper(mburst, param, c_bol=1.0)

Calculates the persistent flux, based on the supplied mdot, redshift etc. Earlier we passed the individual parameters, but easiest to just provide the model burst

concord.burstclass.lhoodClass(params, obs, model, **kwargs)

Calculate the likelihood related to one or more model-observation comparisons The corresponding call to emcee will (necessarily) look something like this:

>>> sampler = emcee.EnsembleSampler(nwalkers, ndim, lhoodClass, args=[obs, models, weights])

Use the kwargs construction to pass the additional parameters weights, disc_model to the compare method if required

Parameters:

  • params – tuple with the distance, inclination, redshift and as many time offsets as there are bursts to compare

  • obs – one or more :py:class:`concord.burstclass.ObservedBurst`s

  • model – one or more model bursts, one for each of the observed bursts; for example a concord.burstclass.KeplerBurst

Returns:

likelihood value

concord.burstclass.modelFunc(p, obs, model, disc_model)

This function performs the stretching and rescaling of the (model predicted) burst lightcurve based on the input parameter set The calculations are based on those in Appendix B of Keek & Heger (2011, ApJ 743, #189)

Parameters:

  • p – parameter array, (for now) a tuple with the appropriate units; ( distance, inclination, redshift, time offset )

  • obs – observed burst to use for the time bins

  • model – model burst to be stretched and rescaled

  • disc_model – choice of disc model for the anisotropy

concord.burstclass.plot_comparison(obs, models, param=None, sampler=None, ibest=None)

This routine will plot each observation vs. each model, for a given set of parameters. If you supply an emcee sampler object, it will find the highest-likelihood sample and use that

concord.burstclass.plot_contours(sampler, parameters=['$d$', '$i$', '$1+z$'], ignore=10, plot_size=6)

Simple routine to plot contours of the walkers, ignoring some initial fraction of the steps (the “burn-in” phase)

Documentation is here https://samreay.github.io/ChainConsumer/index.html

Indices and tables