# Signatures of an eccentric disc cavity: Dust and gas in IRS 48

Josh Calcino1, Daniel J. Price2, Christophe Pinte2,3, Nienke van der Marel 4,
Enrico Ragusa5, Giovanni Dipierro5, Nicolás Cuello6,7, Valentin Christiaens2
1School of Mathematics and Physics, The University of Queensland, QLD 4072,
Australia
2Monash Centre for Astrophysics (MoCA) and School of Physics and Astronomy,
Monash University, Clayton Vic 3800, Australia
3Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France
4Herzberg Astronomy & Astrophysics Programs, National Research Council of
Canada, 5071 West Saanich Road, Victoria BC V9E 2E7, Canada
5Department of Physics and Astronomy, University of Leicester, Leicester,
United Kingdom
6Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Santiago,
Chile
7Núcleo Milenio de Formación Planetaria (NPF), Chile
Contact e-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We test the hypothesis that the disc cavity in the ‘transition disc’ Oph IRS
48 is carved by an unseen binary companion. We use 3D dust-gas smoothed-
particle hydrodynamics simulations to demonstrate that marginally coupled dust
grains concentrate in the gas over-density that forms in in the cavity around
a low binary mass ratio binary. This produces high contrast ratio dust
asymmetries at the cavity edge similar to those observed in the disc around
IRS 48 and other transition discs. This structure was previously assumed to be
a vortex. However, we show that the observed velocity map of IRS 48 displays a
peculiar asymmetry that is not predicted by the vortex hypothesis. We show the
unusual kinematics are naturally explained by the non-Keplerian flow of gas in
an eccentric circumbinary cavity. We further show that perturbations observed
in the isovelocity curves of IRS 48 may be explained as the product of the
dynamical interaction between the companion and the disc. The presence of a
$\sim$0.4 M⊙ companion at a $\sim$10 au separation can qualitatively explain
these observations. High spatial resolution line and continuum imaging should
be able to confirm this hypothesis.

###### keywords:

protoplanetary discs — circumstellar matter — stars: individual: Oph IRS 48 —
methods: numerical — hydrodynamics

††pubyear: 2019††pagerange: Signatures of an eccentric disc cavity: Dust and
gas in IRS 48–Signatures of an eccentric disc cavity: Dust and gas in IRS 48

##  1 Introduction

Protoplanetary discs with cleared out central cavities, i.e.  _transition
discs_ , are promising sites for characterising planet formation. The cleared
central regions may indicate the presence of companions (Strom et al., 1989;
Marsh & Mahoney, 1992, 1993), but photoevaporation may also explain this
feature (Alexander et al., 2014; Turner et al., 2014) — the latter hypothesis
giving transition discs their name.

Resolved observations of transition discs with the Atacama Large
Millimetre/submillimetre Array (ALMA) have revealed large asymmetries in
(sub-)millimeter continuum observations, with Oph IRS 48 (van der Marel et
al., 2013, 2015), HD 142527 (Casassus et al., 2015), and AB Aur (Tang et al.,
2017) being a few notable examples. Dust asymmetries around the cavity are
commonly attributed to the presence of a gap-edge vortex (van der Marel et
al., 2013; Zhu et al., 2014; Fuente et al., 2017). Vortices may form due to
Rossby wave instability (RWI) at the edge of a radial pressure bump in the
disc. They are thought to act as dust traps (Barge & Sommeria, 1995). The
vortex is generally assumed to be induced by a planetary-mass companion
internal to the vortex (e.g. see van der Marel et al., 2013; Fuente et al.,
2017), although the presence of dead zones in the disc have been shown to
produce vortensity gradients that can trigger the formation of vortices via
the RWI (Regály et al., 2012; Ruge et al., 2016). In any case, very low values
of disc viscosity are required for the RWI to be triggered and produce
vortices (Regály et al., 2012; Ataiee et al., 2013). Gas and dust simulations
of vortices have shown that dust grains can strongly concentrate near the
centre of the vortex, creating the asymmetries observed in IRS 48 and other
discs (Birnstiel et al., 2013; Fuente et al., 2017).

The appearance and stability of vortices depends strongly on the disc
viscosity (Ataiee et al., 2013), with a higher viscosity leading to shorter
vortex lifetimes (Fu et al., 2014a; Zhu et al., 2014), or no vortex formation
at all. The effect of dust feedback on the gas can also suppress vortex
formation, and disrupt vortices that may form in the disc (Inaba & Barge,
2006; Fu et al., 2014b), although recent 3D simulations have shown that this
might not be the case (Lyra et al., 2018). It has also been shown that vortex
lifetimes are significantly reduced when simulations are initiated with lower-
mass planets which undergo runaway gas accretion to reach their final mass
(Hammer et al., 2017). Despite these challenges, vortices are still the most
commonly invoked mechanism for explaining high contrast ratio dust asymmetries
in transition discs.

Vortices, however, are not the only pathway to generate asymmetries in the gas
and dust distribution in a circumstellar disc. By dropping the requirement of
low viscosity values, Jupiter mass planets have been shown to form eccentric
cavities characterized by overdense features at the apocentre of the disc
cavity, as a consequence of the clustering of eccentric orbits (Ataiee et al.,
2013); however, the density contrast ratio of these eccentric features is much
lower than those found in observations.

Another distinct mechanism is the presence of a (sub-)stellar companion that
launches tidal streams of gas across the cavity, producing a pileup of
material at its edge (Ragusa et al., 2017). More massive companions were found
to produce stronger asymmetries at the cavity edge. The important distinction
between this mechanism, and a purely eccentric disc as explored by Ataiee et
al. (2013) is that the over-density is comoving with the gas around the
cavity. This feature has been observed in numerical simulations focusing on
black hole circumbinary disc simulations (Farris et al., 2014; D’Orazio et
al., 2016; Ragusa et al., 2016; Muñoz & Lai, 2016; Miranda et al., 2017), and
is seen to persist after several thousand orbits of the central binary.

The hypothesis of a binary-disc interaction was tested by Price et al. (2018b)
for the transition disc HD142527 — a disc which shows a strong azimuthal
asymmetry in dust continuum emission (Casassus et al., 2013). Price et al.
(2018b) demonstrated that many of the morphological features of HD 142527 can
be attributed to the observed companion (Biller et al., 2012; Lacour et al.,
2016), and most importantly, showed that vortices are not the only contender
in explaining dust asymmetries.

Thus far, numerical simulations have demonstrated that binary companions can
generate dust asymmetries with a contrast ratio of $\sim$10 (Ragusa et al.,
2017; Price et al., 2018b). However, transition discs have been observed with
contrast ratios in far excess of this value. The most famous example is that
of IRS 48, which has a contrast ratio greater than a factor of $\sim$100 (van
der Marel et al., 2013, 2015).

In this paper we demonstrate that large dust grains can concentrate in the gas
over-density that develops around high mass ratio companions on circular, co-
planar orbits. This enhances the dust contrast in these systems. We
investigate the kinematic signatures that this companion imprints on the
circumbinary disc and show how these kinematic signatures may be used to
distinguish between the vortex and circumbinary hypotheses in IRS 48.

##  2 Observations of IRS 48

###  2.1 Gas structure

Near infrared observations of IRS 48 by Brown et al. (2012) revealed a ring-
like structure of CO rovibrational emission approximately 30 au from the
central star. More recent observations of CO isotopologue rotational lines
taken with ALMA, tracing the colder gas, showed that CO is strongly depleted
within this ring-like structure by at least a factor of 100 (Bruderer et al.,
2014; van der Marel et al., 2016a), and the 13CO 6–5 emission extends out to
roughly 90 au (van der Marel et al., 2016a).

###  2.2 Large dust grain asymmetry

Large (sub-mm to cm) dust grains observed in IRS 48 show a high contrast ratio
dust asymmetry (van der Marel et al., 2013, 2015), that peaks in intensity at
approximately 60 au from the centre of the disk. The azimuthal extent of the
dust emission also appears to change with dust grain size, with larger grains
more azimuthally trapped than smaller grains (van der Marel et al., 2015, in
prep). Any hypothesis that attempts to explain the nature of the dust trap in
IRS 48 must reproduce this phenomenon.

###  2.3 Ring-like dust features

Ring-like features are observed across a range of wavelengths in the near
infrared within and co-spatial with the large dust grain cavity. Beginning
with the inner disk, there is evidence of dust emission from small grains and
polycyclic aromatic hydrocarbons in a ring-like feature within the CO cavity
(Geers et al., 2007; Schworer et al., 2017; Birchall et al., 2019).

Further from the central star an East-West asymmetry has also been observed in
18.6 $\mu$m thermal emission (Geers et al., 2007; Honda et al., 2018). It was
claimed by Honda et al. (2018) that this thermal emission is not due to
differences in the thermal structure between the East and West side of the
disk. However, it is possible that there are temperature differences less than
$\sim$ 45 K.

##  3 Methods

We performed 3D smoothed particle hydrodynamical (SPH) simulations using the
code phantom (Price et al., 2018a). We explore the evolution of the gas disc,
adding dust grains once the system has reached a quasi-steady state. The dust
is coupled to the gas via aerodynamic drag, but we neglect the back reaction
of the dust on the gas since the evolution of the gas disc is dominated by the
gravitational interaction between the central star and companion.

We model the dust particles using a two-fluid approach, where the gas and dust
are modelled as two separate SPH populations (Laibe & Price, 2012). This
method is efficient for simulating dust grains with high Stokes number which
are decoupled from the gas. This is an important difference between our dusty
circumbinary disc simulations compared with those of Ragusa et al. (2017), who
used a one-fluid approach assuming dust grains with low Stokes number (Laibe &
Price, 2014; Price & Laibe, 2015).

In our dust simulations we model the dynamics of two grain populations with
grain sizes $s_{\mathrm{grain}}$ of 100 $\mu$m and 1 mm, respectively, to
study the concentration of dust grains in the pressure maximum. We assume a
grain density $\rho_{\mathrm{grain}}=3$ g/cm3. Since the Stokes number is the
physically relevant parameter that governs the dust evolution, the simulation
results would remain identical if we were to scale the size and density of the
grains such that $s_{\mathrm{grain}}\rho_{\mathrm{grain}}$ remains constant
(see Eq. 3 of Dipierro et al., 2015). For this reason we refer to the 100 $\mu
m$ and 1 mm as ‘medium’ and ‘large’ grains, respectively.

We used a separate simulation for each dust species. Since we neglect the back
reaction of the dust onto the gas the final result would be the same if we
simulated both grain sizes together. Separate simulations are necessary simply
because simultaneous evolution of multiple types of dust particles is not yet
available in the code.

We model the primary and secondary stars as sink particles (Bate et al.,
1995), which interact with each other and the gas and dust disc
gravitationally. Particles can accrete onto both sink particles provided the
SPH particles around the sink are gravitationally bound and within a specified
accretion radius.

###  3.1 Initial Conditions

Our simulations contain a central protostar with a mass of 2 M⊙, consistent
with the mass estimated for IRS 48 (Brown et al., 2012). We place a companion
with a secondary-to-primary mass ratio $M_{\textrm{s}}/M_{\textrm{p}}=0.2$,
giving a companion mass of 0.4 M⊙. We use accretion radii of 2 au and 0.5 au
for the primary and companion stars respectively. These sink radii sizes mean
that we do not resolve the accretion discs around each individual sink, but
this is necessary to minimise computational time. The companion is initialised
with a semi-major axis of 10 au and is initially on a circular, co-planar
orbit. We discuss the consistency of a companion with this mass and semi-major
axis in section 5.1. The binary is free to evolve during the simulation, but
given our small initial gas mass, the change in the binary orbital elements
remains small.

We use $5\times 10^{6}$ SPH particles to model the gas disc, assuming a gas
mass of $5\times 10^{-4}$ M⊙, consistent with gas mass estimate of $5.5\times
10^{-4}M_{\odot}$ inside 90 au of IRS 48 (van der Marel et al., 2016a). We
found in hindsight that increasing the gas mass by a factor of 5 was needed to
reproduce the CO flux (see below). We discuss how this increase in gas mass
remains consistent with the mass estimate by van der Marel et al. (2016a) in
section 4.2. The gas is initially set up in Keplerian rotation in a co-planar
annulus with the inner and outer radius as $R_{\textrm{in}}=15$ au and
$R_{\textrm{out}}=40$ au, respectively. The radial extent of the gas in the
disc dramatically increases as the simulations progress. We set the Shakura-
Sunyaev alpha viscosity to $\alpha_{\textrm{SS}}\approx 1.5\times 10^{-3}$ by
employing a constant SPH artificial viscosity parameter $\alpha_{\rm
AV}=0.11$. The surface density of the gas is initialised as a power law with
$\Sigma\propto R^{-p}$, with $p=1$. The temperature profile of the disc is
locally isothermal with $T(R)\propto(R/R_{\textrm{in}})^{-2q}$ and $q=0.25$.
The aspect ratio of the disc is determined by the $q$ index where
$H/R=H_{\textrm{in}}/R_{\textrm{in}}(R/R_{\textrm{in}})^{\frac{1}{2}-q}$, and
$H_{\textrm{in}}=0.05$.

The gas disc is allowed to evolve and reach a quasi-steady state, taking on
the order of $\sim$500 binary orbits to do so. To study the effects of dust
grain migration and concentration in the dust trap, after 500 orbits we added
dust grains with an initial density distribution set from the gas density
distribution at that time assuming a dust to gas ratio of 1:100 for each grain
size. Since we neglect the back reaction of the dust onto the gas, the amount
of dust we add to the simulation is arbitrary.

We use a gas to dust particle ratio of 30:1 to prevent the dust particles from
getting trapped under the gas particle smoothing length (Laibe & Price, 2012).
The gas+dust simulations are run for a further $\sim$30 orbits of the binary.
This is a rather short period of time, but the addition of dust particles
decreases the timestep size of the simulations drastically. This limits the
duration for which we are able to track the evolution of the dust grains.

###  3.2 Radiative transfer and synthetic observations

We performed synthetic observations of our models using the Monte Carlo
radiative transfer code mcfost (Pinte et al., 2006; Pinte et al., 2009).

We employed a dust model based on the result of the SPH simulation by assuming
a dust population with a power-law grain size distribution given by
$dn/ds\propto s^{-m}$ between $s_{\min}$ to $s_{\max}$. We use the simulation
that contains the medium sized dust grains since very little continuum
emission comes from the much larger dust grains at the relevant wavelengths.
The dust opacity was computed assuming spherical and homogeneous dust grains.
There is evidence of grain growth in the dust trap of IRS 48 (van der Marel et
al., 2015), so for our radiative transfer calculations we assume
$s_{\min}=1\,\mu m$ and $s_{\max}=5\,mm$, with $m=3.0$. The total gas mass is
taken from our SPH simulations. The dust mass is obtained by assuming a gas to
dust ratio of 10, consistent with dust mass estimates in IRS 48 (van der Marel
et al., 2016b). We found that to reproduce the integrated 13CO emission, we
needed to increase the gas mass by a factor of $\sim$5, so that the effective
disc mass is $2.5\times 10^{-3}$M⊙, assuming a 13CO-to-H2 abundance ratio of
$2\times 10^{-6}$. When producing our integrated 13CO emission images, we use
this increased gas mass and use a gas to dust ratio of 10. The dust grain size
and density is scaled such that the Stokes number of the grains remains
unchanged when we scale the gas mass. Since the gas mass is small compared to
the mass of the binary, the feedback of the gas onto the binary is negligible.

We assume that the dust and gas are in thermal equilibrium and the dust
opacities are independent of temperature. We model the emission of the central
sources as black bodies, assuming an effective temperatures of
$T_{\textrm{p}}=10000$ K (Brown et al., 2012) and $T_{\textrm{s}}=3800$ K.

We used $10^{8}$ Monte Carlo photon packets to compute the temperature and
specific intensities. Images were then produced by ray-tracing the computed
source function.

We assumed an inclination of $i=50^{\circ}$, a position angle of
$\textrm{PA}=100^{\circ}$ (Bruderer et al., 2014), and a source distance of
134 pc (Gaia Collaboration et al., 2018) to create our simulated images of IRS
48. When comparing our simulated images with observations of IRS 48, we
convolved our images with a Gaussian beam to match the beam size of the
observations.

##  4 Results

Figure 1: Gas surface density of our circumbinary disc simulation from 500–506
binary orbits ($P_{B}$). The gas over-density orbits around the eccentric
cavity at roughly the Keplerian velocity, and is not the result of a pile up
of material at the apocentre of the eccentric disc, as in the case studied by
Ataiee et al. (2013). Accretion flows on the binary stars appear as ring-like
structures, and are generated every orbit of the binary companion. Outside of
the cavity, the influence of the binary companion is still evident through the
spiral density waves that are propagating through the disc.  Figure 2: Gas and
dust column densities in our SPH simulations after $\sim$522 binary orbits.
Dust was added after 500 binary orbits, once the eccentric cavity had reached
a quasi steady-state, with dust density initially proportional to the gas
surface density with a ratio of 1:100. The scale between the gas and the dust
is shifted by a factor of 100. Therefore, we can directly see that the dust
has concentrated if its surface density is elevated compared to that of the
gas, which is clearly the case for the large dust grains. Note that since we
are neglecting the back reaction of the dust onto the gas, the total amount of
dust we add into the simulation is arbitrary. The dust to gas ratio seen in
this Figure does not reflect the dust to gas ratio used when conducting the
radiative transfer modelling.

Figure 1 shows the gas surface density after $\sim$500 binary orbits
($P_{B}$). The over-dense feature seen in the gas and dust is comoving with
the fluid around the cavity, and is therefore different from the “traffic jam”
of material that forms around an eccentric disc (Ataiee et al., 2013). The
over-density also lacks vortical motion (Ragusa et al., 2017), and so is
distinct from a vortex. The precise mechanism that leads to the formation of
the cavity-edge over-density, its longevity and dependence on disc parameters,
are not fully understood. We discovered that the instability leading to the
cavity-edge over-density takes longer to initiate when the initial disc inner
radius is large, but the instability does still form once the inner disc
viscously evolves and reduces in radius. So although we used a close inner
radius of $R_{\textrm{in}}=15$ au, it would not have made a large difference
to the final outcome if we used a larger inner radius. We use a smaller inner
radius since we spend less time evolving the system to reach the stage of
developing the instability.

We show the outcome of adding dust into our simulations in Figure 2. Since the
gas surface density changes by about an order of magnitude around the
circumbinary disc, the Stokes number for each dust grain varies depending on
its location in the disc. The medium sized dust grains (middle panel) have a
Stokes number between 0.01—0.1, whereas the large dust grains have a Stokes
number 0.1–1, depending if the grains are located within ot outside of the gas
over-density. Despite only being present in the simulation for $\sim$22
$P_{B}$, there is rapid migration and concentration in the gas pressure
maximum for the large dust grains. Our main limitation preventing longer
evolution time was the drastic increase in computation expense as the large
dust grains concentrated. We were able to simulate the medium dust grains for
a much longer period and observed an increased level of radial drift and
concentration in the gas pressure maximum for these grains. We discuss this
limitation more in Section 5.2.

It is evident in Figure 2 that the dust grains have concentrated inside of the
gas over-density. The initial dust distribution was set to be proportional to
the surface density of the gas, with a gas to dust ratio of 100:1. The
colorbar is scaled by this factor in Figure 2 between the gas (left panel) and
dust (middle and right panels). Therefore, we can see that the dust has
concentrated if the surface density has increased, which is most obvious in
the large dust grains.

One revealing feature of Figure 2 is the difference in azimuthal extent
between the medium and large dust grains. The larger dust grains are more
radially and azimuthally concentrated. This feature is in agreement with the
more compact azimuthal distribution observed at longer wavelengths by van der
Marel et al. (2015). We show the dust grain surface density of our simulations
smoothed to a beam size of 0.1" in Figure 3, assuming the simulation is placed
at the at a distance of 134 pc. This Figure makes it more clear that we are
able to achieve a high dust contrast ratio approaching 100:1 in the medium
sized dust grains, and exceeding 100:1 in the large dust grains.

Figure 3: Medium and large dust grain surface density distributions smoothed
using a a beam size of 0.1" assuming the disc is observed at the distance of
IRS 48. Dust grain simulations are the same as in middle and right panels of
Figure 2. Medium sized grains concentrate inside the gas over-density, where a
contrast ratio approaching 100:1 is achieved. Large dust grains concentrate
more than the medium sized dust grains, exceeding a maximum contrast ratio of
100:1. The asymmetry in the medium dust grain surface density is more evident
than in Figure 2. Large dust grains display a more symmetric distribution.
Similar asymmetry is seen in multi-wavelength observations of IRS 48 (van der
Marel et al., 2015, in prep.).

Our simulated dust cavity is smaller ($\sim$20 au) than observed in IRS 48.
The gas cavity that develops in a circumbinary disc is sensitive to the mass
ratio, semi-major axis, eccentricity, and inclination of the companion
(Miranda & Lai, 2015; Thun et al., 2017). We did not explore a full range of
possible orbits for the companion, but exploring this parameter space is a
promising avenue to resolve this discrepancy. Given this, we scaled our code
units prior to conducting the radiative transfer calculations to better match
the observed cavity size in IRS 48. The code units are scaled such that the
length unit was increased by 30%, and the time unit scaled by an amount to
match the change in length scale, according to equation 33 in Price et al.
(2018a). This means that the semi-major axis of our binary companion is scaled
to 13 au.

###  4.1 Dust Continuum

Figure 4: Comparison of our synthetic 690 GHz continuum emission observations
(left panel) with the 690 GHz continuum emission presented in van der Marel et
al. (2015); van der Marel et al. (2016a) (right panel). All scales are in
normalised flux. The left panel shows the dust continuum emission when we
match the kinematics of our simulation with that of IRS 48 (see Figure 6). Our
synthetic observations have been convolved with a Gaussian beam to match the
observational data (white ellipse in the bottom left of each panel). Our
simulated cavity size is smaller than that of IRS 48, however better agreement
could be obtained by allowing the gas only simulation to advance for longer
before adding in dust. The cavity size also depends on the orbit of the
companion. Choosing a higher mass ratio for the companion (either by
increasing the companion mass, or decreasing the primary mass), or increasing
the eccentricity of the companion orbit, are two ways that we could increase
the cavity size. We discuss this point further in section 5.2.

Figure 4 compares our synthetic observations (left panel) to the 690 GHz
continuum emission presented in van der Marel et al. (2015); van der Marel et
al. (2016a) (right panel). The medium sized dust grains have been allowed to
evolve for $\sim 45P_{B}$ prior to generating the synthetic observations. We
recreate the asymmetric continuum emission observed in IRS 48 with modest
agreement in density contrast and spatial distribution of the dust grains.
Some notable differences between our simulated images, and the observed dust
emission in IRS 48, are the radial distance and extent of the emission. We
discuss ways in which we can increase the cavity size of our circumbinary disc
in section 5.2. The radial extent of the dust trap in IRS 48 is resolved at
the wavelength of the observations in Figure 4. A substantial amount of the
emission in these observations could be coming from smaller dust grains which
have a lower Stokes number than the grains in our simulations. From our
simulations containing the medium and large size dust grains, we can see that
grains with Stokes numbers closer to unity tend to be more radially and
azimuthally concentrated. Following the evolution of smaller dust grains in
our simulations may relieve this discrepancy since they are expected to have a
larger radial extent.

As the dust over-density orbits around the eccentric cavity, its observational
appearance changes, closely following the changes in the gas over-density
morphology seen in the panel of Figure 1. The left panel of Figure 4 shows our
best match to the observations of the right panel. In this simulated image,
the dust grains are travelling towards the apastron of their orbit, and is in
a similar position to the gas over-density in the bottom right panel of Figure
1.

Although our synthetic observations in Figure 4 are not able to completely
match the contrast ratio observed in IRS 48, we have plausible explanations
for how this may be achieved. Firstly, in our simulations we studied dust
grains that are marginally coupled to the gas motion. We observed rapid radial
and azimuthal concentration of the dust into the gas over-density for the 1mm
dust grains. If it were feasible to run these simulations for a longer period
of time, we see no reason why this trend would not continue, provided the dust
back-reaction on the gas remains negligible. We discuss this more in section
5.2.

It is also reasonable to expect that modelling dust growth and fragmentation
would affect the dust contrast ratio for a given grain size. Dust growth and
fragmentation is largely a function of the relative velocity of the dust
grains (Testi et al., 2014). Grains cease growing if their relative velocity
exceeds the fragmentation barrier. In an eccentric disc the relative velocity
of dust grains would likely be a function of position around their orbit of
the cavity. This may lead to either growth or fragmentation of dust within the
gas over-density of the circumbinary disc, depending on where it is on its
orbit around the cavity. How this affects the appearance of the azimuthal
asymmetry remains unknown.

###  4.2 Integrated 13CO Emission

Figure 5: Comparison of our 13CO (6-5) emission (left panel) with a comparison
to the 13CO (6-5) emission of IRS 48 (right panel) shown in van der Marel et
al. (2016a). The 13CO (6-5) emission in IRS 48 displays a strong drop in
brightness temperature co-spatial with the dust trap. This is likely due to
the 13CO emission being absorbed by the dust in the dust trap. We found that
in order to reproduce this effect, a large amount of dust is required.

Figure 5 (left panel) shows the integrated 13CO 6–5 emission from our
simulation, with a comparison to the integrated 13CO emission of IRS 48
presented in van der Marel et al. (2016a). One feature worth discussing is the
observed drop in 13CO brightness temperature co-located with the dust emission
asymmetry. This feature is likely caused by the high concentration of dust
grains at this location. The dust grains absorb the 13CO emission, reducing
the brightness temperature in this region. This effect has two consequences.

First, the strong absorption of 13CO emission implies a very large
concentration of dust at this location. In order for us to reproduce this
effect in our radiative transfer calculations, we required a low gas-to-dust
ratio of 10:1, consistent with the gas-to-dust ratio estimated by van der
Marel et al. (2016a). Note that this is a global gas-to-dust ratio, within the
dust trap it is almost certainly much higher. Such a vast amount of dust is
required in this location that in our radiative transfer modelling, the dust
over-density becomes optically thick. This means that the dust is actually
more concentrated in the dust trap of IRS 48 than is implied by the dust
emission observations. The level of 13CO emission absorption is also sensitive
to the choice of the dust grain population distribution.

Second, a depletion in 13CO emission will reduce the estimated amount of 13CO
in the disc. van der Marel et al. (2016a) estimated the disc mass of IRS 48 to
be $\sim 0.5$MJ, but this estimate would not account for the 13CO within the
dust trap that is not being detected due to the absorption of 13CO emission by
the dust grains. Therefore, the mass estimate by van der Marel et al. (2016a)
really measures the mass contained outside of the dust trap.

In order for us to be able to reproduce the dust absorption of the 13CO, and
match the brightness temperature of the 13CO, we needed to increase the gas
mass of our simulations by a factor of $\sim 5$. The majority of the disc mass
in our simulation is contained within the gas over-density (see left panel of
Figure 2), and therefore this increase in gas mass is not inconsistent with
the estimated disc mass outside of the dust trap by van der Marel et al.
(2016a).

We recreate the 13CO emission absorption feature in the synthetic image of our
simulation in the left panel of Figure 5. Since our dust distribution does not
perfectly match that observed in IRS 48, the 13CO absorption feature is also
not perfectly recreated. However we demonstrate that our model is capable of
reproducing this feature.

We also note that the location of the 13CO emission absorption feature (and
hence the dust over-density) is located outside of an inner 13CO ring. This is
agreement with the observation of an inner CO ring observed in IRS 48 (Brown
et al., 2012; van der Marel et al., 2016a).

###  4.3 Velocity Maps

Figure 6: Comparison of our synthetic velocity map (left panel) with the
velocity 13CO (6–5) map of IRS 48 shown in van der Marel et al. (2016a) (right
panel). We convolve our synthetic map with a Gaussian beam matching the
observational resolution (black ellipse in the bottom left of each panel). The
kinematic map of IRS 48 shows asymmetry between the blue and red wings. In our
synthetic map this is caused by the eccentric disc cavity and fast flowing gas
inside the cavity. If the eccentricity vector of the eccentric disc was
pointing along our line of sight, the kinematic map would show a more
symmetric blue and red wing. As the eccentricity vector is offset to the line
of sight in our synthetic map, an asymmetric blue and red wing is seen.

A peculiar asymmetry is observed in the kinematic map of 13CO (6–5) in IRS 48
(Figure 6, right panel). If this asymmetry is not an observational artefact,
then it must be explained in models of vortices and eccentric discs alike.

The left panel of Figure 6 shows our simulated kinematics. Our model naturally
produces asymmetric kinematics due to the eccentric orbit of the gas. We also
found that the position of the companion relative to the eccentric disc can
leave an imprint in the velocity maps. This signal is much smaller than that
of the eccentric disc itself, but may be observable in high spatial resolution
molecular emission line observations. The asymmetry is observed from most
azimuthal viewing angles, and becomes weaker when the eccentricity vector of
the disc is pointing towards or away from the line of sight. In our model of
IRS 48, the eccentricity vector of the disc is orientated nearly in the East-
West direction. High velocity accretion flows onto the stars inside the cavity
produces the central East-facing point in the red wing of the velocity maps in
Figure 6. How exactly the kinematic maps change due to the orientation and
inclination of the disc should be explored in future work.

With the sensitive nature of vortex formation in mind, it must be established
whether or not they could form in discs that produce the observed kinematics
of IRS 48. Recent work (Pérez et al., 2018; Huang et al., 2018) demonstrates
that simulated discs that contain vortices show kinematic maps that are mostly
symmetric, and the signature from the vortex is much smaller than what we
predict for a circumbinary disc. For the vortex interpretation of the dust
trap in IRS 48 to remain viable, it must be able to predict the asymmetry
observed in the kinematic maps as well.

###  4.4 Channel Maps

Figure 7: Comparison between our simulated channel maps convolved with two
different beam sizes (top and middle rows), with a spectral resolution of
0.489 km/s, and the observed channel maps of IRS 48 (bottom row) presented in
van der Marel et al. (2016b) for selected channels. Observed channel maps show
deviations from Keplerian rotation. In particular, the angles at which the
winglets in each channel propagate are not what is expected for a disk in
circular rotation. For example, the rest frame channel of the system should
have both winglets propagating at same angle as the position angle of the
system ($10^{\circ}$ in this case), but it can be seen that they are not. The
channel maps of our circumbinary disc simulation also show these features. The
left and right wing are more compact in the observations than they are in our
simulations, indicating our simulated total binary mass is higher than that in
IRS 48, or that the inclination of the system is not as high as currently
estimated (50∘).

Additional evidence for a circumbinary disc in IRS 48 is found in the
individual channel maps of the 13CO (6–5) emission line observations presented
in van der Marel et al. (2016a). A companion of stellar mass will cause
perturbations in the isovelocity contours of a disc. We also found that the
gas over-density propagates spiral density waves as it orbits around the
circumbinary disc, which can leave a distinct signature in the disc
kinematics. In the middle panel of Figure 7 we show channel maps from our
circumbinary disc simulation at a spectral resolution of 0.489 km/s, and
convolved with the same beam size as the observations of IRS 48 in the bottom
panel. Even with the companion located at $\sim$ 10 au, its presence is
discernible in the kinematics of the gas at the cavity edge and outer disc.
The top panel of Figure 7 shows the kinematics of our simulations as they
would appear with a smaller beam size of 0.05". The ‘kinks’ observed in the
iso-velocity curves are being generated by some combination of the companion
perturbing the disc, and the gas over-density propagating spiral density
waves.

The most convincing evidence of the disc around IRS 48 being eccentric is the
different angles at which the northern and southern winglets propagate, most
noticeable in the channel that corresponds to the rest frame velocity of IRS
48 (middle column of Figure 7). If the disc around IRS 48 were in circular
Keplerian rotation, both winglets would be pointing in an anti-parallel
direction, and the angle at which they point would be equal to the position
angle of the disc (estimated to be 10∘). This is not the case here. One
explanation for this is that the disc around IRS 48 is has an eccentric shape,
and hence a velocity field that deviates strongly away from the velocity field
of a disc with gas in a circular orbit.

Our simulation reproduces this feature due to the eccentricity that develops
in a circumbinary disc. We found that rotating the simulation changes the
angle at which the northern and southern winglets propagate. Placing the disc
such that the eccentricity vector of the disc is pointing roughly along the
East-West direction gives good agreement. Note that our simulated channel maps
are the same as those used to produce the velocity map shown in Figure 6, and
also the 13CO 6-5 integrated emission in Figure 5.

The blobby appearance in the southern wing of our simulation and in the
observations arises due to the 13CO emission absorption discussed in section
4.2.

The channel maps of our circumbinary disc simulation agree with the observed
channel maps of IRS 48. The isovelocity curves in IRS 48 appear to be
perturbed in a similar fashion to what is observed in our circumbinary disk
model. Observing these features at higher spatial resolution would provide
very strong evidence in favour of our circumbinary disc hypothesis of IRS 48.

Recent work by Huang et al. (2018) made predictions for the imprint a vortex
would leave on the channel maps and velocity map of IRS 48. The channel maps
presented in their paper (their Figure 4) do not produce the peculiar
structures seen in the channels maps of IRS 48. Furthermore, the velocity map
presented (third row of their Figure 6) does not display an asymmetry between
the blue and red-shifted sides of the disk. This does not exclude the
possibility that these features can be produced by vortex models, but it shows
that more work is required in understanding whether or not they are capable of
being produced. The circumbinary disk model is able to reproduce both of these
features.

One disagreement between our simulated channel maps and the observed channel
maps of IRS 48 is the rate at which the left and right wings close, most
noticeable in the left and right-most panels. This could be due to our sink
particle masses being higher than the binary mass in IRS 48. We could achieve
better agreement with the observed channel maps by reducing the mass of the
two stars. This would not affect the development of the eccentric disc since
this is sensitive to the mass ratio, and not the absolute mass, of the stellar
companions (Miranda & Lai, 2015).

##  5 Discussion

###  5.1 Consistency with Observations

Binary searches in IRS 48 have resulted in no detections in the $\sim$18-800
au regime down to $m_{K}$=10.1 (Ratzka et al., 2005) and no detections in the
$\sim$2.5-1400 au regime down to $m_{K}$=9.8 (Simon et al., 1995). These
limits were converted to mass estimates on the order of 100-150 MJ (Wright et
al., 2015), below our companion mass of 0.4 M⊙ at 13 au. However, IRS 48 is
strongly extincted, making the detection limits highly uncertain and a massive
companion may remain hidden (Wu & Sheehan, 2017; Schworer et al., 2017),
potentially explaining why a direct detection has not yet been made. However,
a recent companion search in IRS 48 has detected an asymmetry in the
brightness profile of an inner dusty disc around IRS 48 at $\sim$15 au, a
location consistent with the companion in our simulation (Birchall et al.,
2019).

Ring-like structures have been detected in the inner regions of the disc,
close to our hypothetical stellar companion (Geers et al., 2007; Brown et al.,
2012; Bruderer et al., 2014; Birchall et al., 2019). Accretion flows from the
outer disc fall inside the cavity, producing ring-like structures, which can
be seen in Figure 1. In this Figure we show how the gas over-density at the
cavity edge orbits at roughly the Keplerian velocity. The inner binary orbits
on a shorter timescale, and the position and structure of the accreation flows
inside the cavity vary on this timescale.

Furthermore, both the primary and companion stars will develop their own
accretion discs. We do not resolve these structures in our simulations since
they largely fall within the accretion radii of our sink particles,
particularly for the primary star which has a larger sink radius. Studying the
evolution and structure of the circumprimary and circumsecondary discs with
the circumbinary disc is difficult with SPH due to the high numerical
viscosity and low resolution present at low particle numbers.

###  5.2 Numerical Limitations

In order to get a better agreement with the observed cavity size in IRS 48, we
needed to increase the scale length of our simulations. Simulations of
circumbinary discs by Thun et al. (2017) indicate that a cavity is opened on a
dynamical timescale of 10–20 orbits, but also that the eventual cavity size is
set on the viscous timescale ($\sim$several thousand binary orbits). Given
that we are running 3 dimensional simulations, it is challenging to run our
simulations for a duration approaching the viscous timescale. We expect that
our cavity would continue to increase in size and may eventually be consistent
with the cavity size observed in IRS 48 (see Figure 9 of Thun et al., 2017),
and the increase in scale length would no longer be necessary. Other works
investigate the evolution of circumbinary discs over much longer timescales
(e.g. Miranda et al., 2017; Thun et al., 2017), however these tend to be done
with two dimensional simulations which are cheaper.

Adding dust grains also increases the computational time of our simulations,
since the timestep is now limited to be a fraction of the stopping time of the
simulated dust species. Furthermore, the stopping time drops as the dust
grains concentrate in the pressure maximum and around the cavity. The net
result is that eventually we are no longer able to follow the evolution of the
dust, particularly for the grains which migrate the fastest (e.g. the large
grains in our simulation). As mentioned previously, we were able to simulate
the medium dust grains for a longer duration, which was due to the fact that
they took longer to concentrate and for a small portion to become trapped
under the gas smoothing length.

Adjusting parameters such as the companion/primary mass ratio, the semi-major
axis and eccentricity of the companion could also provide us with better
agreement to the observed cavity size in IRS 48. However, performing parameter
space searches with three-dimensional dusty simulations is difficult. Even
without an extensive exploration of the parameter space, we can reproduce
features of the dust trap such as contrast ratio and azimuthal trapping
(Figures 2 and 4), the asymmetric velocity map (Figure 6), and substructure in
the channel maps (Figure 7). For this reason we believe that the simulation
limitations discussed in this section do not substantially change our
conclusion that IRS 48 is a circumbinary disc.

##  6 Conclusions

Our model of IRS 48 demonstrates that high contrast ratio asymmetries in dust
continuum emission can be created in a circumbinary disc. Our simulations only
include a stellar companion on a circular, co-planar orbit. Predictions of our
circumbinary disc hypothesis are as follows:

  1. 1.

We predict that there is a stellar companion in IRS 48 with an approximate
mass of $\sim$ 0.4$M_{\odot}$ and with a semi-major axis of $\sim$ 10 au.
These values may differ if the companion is on an eccentric and/or inclined
orbit, but for now we limit our analysis by assuming the companion to be on a
nearly circular, co-planar orbit.

  2. 2.

The interaction of this companion with the primary star causes the disc around
the binary system to become eccentric, and an over-density to form which is
co-moving with the gas around the cavity (Ragusa et al., 2017). Owing to the
eccentric shape of the disc, the position of the primary star should be offset
from the projected centre of the gas disc.

  3. 3.

Due to the non-Keplerian flow of the gas and dust in the eccentric disc, we
predict that the velocity maps of IRS 48 will be asymmetric (see Figure 6),
with the degree of asymmetry growing with higher spatial resolution
observations. We demonstrate that the isovelocity curves in the individual
channel maps will display deviations and kinks (see Figure 7). The deviations
should be more easily detected in higher spatial resolution observations.

It is not clear how the dust grain distribution and velocity fields of the
disc would differ with a companion on a more complicated orbit, although a
first analysis was presented in Price et al. (2018b). This large parameter
space of companion properties could be explored to yield a better match to
observations of IRS 48.

The mechanism for producing high contrast ratio dust horseshoes demonstrated
in this paper has consequences that extend beyond IRS 48. If it is shown that
our proposed stellar companion in IRS 48 indeed exists, high contrast ratio
dust horseshoes can no longer be thought of as sign-posts for planets. Rather,
they may be features of eccentric circumbinary discs around low mass ratio
($\gtrsim$0.2) binary stars.

## Acknowledgements

We thank the anonymous referee for their useful comments which improved the
quality of this manuscript. J.C. acknowledges an Australian Government
Research Training Program Scholarship. C.P., D.J.P. and V.C. acknowledge
funding from the Australian Research Council via FT170100040, FT130100034, and
DP180104235. NC acknowledges funding from FONDECYT grant 3170680 and from
CONICYT project Basal AFB-170002. E.R. and G.D. acknowledge funding from the
European Research Council under the European Union’s Horizon 2020 programme
(grant No 681601). We used SPLASH (Price, 2007).

## References

  * Alexander et al. (2014) Alexander R., Pascucci I., Andrews S., Armitage P., Cieza L., 2014, Protostars and Planets VI, pp 475–496
  * Ataiee et al. (2013) Ataiee S., Pinilla P., Zsom A., Dullemond C. P., Dominik C., Ghanbari J., 2013, A&A, 553, L3
  * Barge & Sommeria (1995) Barge P., Sommeria J., 1995, A&A, 295, L1
  * Bate et al. (1995) Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362
  * Biller et al. (2012) Biller B., et al., 2012, ApJ, 753, L38
  * Birchall et al. (2019) Birchall E. K., et al., 2019, MNRAS, p. 1010
  * Birnstiel et al. (2013) Birnstiel T., Dullemond C. P., Pinilla P., 2013, A&A, 550, L8
  * Brown et al. (2012) Brown J. M., Herczeg G. J., Pontoppidan K. M., van Dishoeck E. F., 2012, ApJ, 744, 116
  * Bruderer et al. (2014) Bruderer S., van der Marel N., van Dishoeck E. F., van Kempen T. A., 2014, A&A, 562, A26
  * Casassus et al. (2013) Casassus S., et al., 2013, Nature, 493, 191
  * Casassus et al. (2015) Casassus S., et al., 2015, ApJ, 812, 126
  * D’Orazio et al. (2016) D’Orazio D. J., Haiman Z., Duffell P., MacFadyen A., Farris B., 2016, MNRAS, 459, 2379
  * Dipierro et al. (2015) Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015, MNRAS, 453, L73
  * Farris et al. (2014) Farris B. D., Duffell P., MacFadyen A. I., Haiman Z., 2014, ApJ, 783, 134
  * Fu et al. (2014a) Fu W., Li H., Lubow S., Li S., 2014a, ApJ, 788, L41
  * Fu et al. (2014b) Fu W., Li H., Lubow S., Li S., Liang E., 2014b, ApJ, 795, L39
  * Fuente et al. (2017) Fuente A., et al., 2017, ApJ, 846, L3
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Geers et al. (2007) Geers V. C., Pontoppidan K. M., van Dishoeck E. F., Dullemond C. P., Augereau J. C., Merín B., Oliveira I., Pel J. W., 2007, A&A, 469, L35
  * Hammer et al. (2017) Hammer M., Kratter K. M., Lin M.-K., 2017, MNRAS, 466, 3533
  * Honda et al. (2018) Honda M., et al., 2018, PASJ, 70, 44
  * Huang et al. (2018) Huang P., Isella A., Li H., Li S., Ji J., 2018, ApJ, 867, 3
  * Inaba & Barge (2006) Inaba S., Barge P., 2006, ApJ, 649, 415
  * Lacour et al. (2016) Lacour S., et al., 2016, A&A, 590, A90
  * Laibe & Price (2012) Laibe G., Price D. J., 2012, MNRAS, 420, 2345
  * Laibe & Price (2014) Laibe G., Price D. J., 2014, MNRAS, 440, 2136
  * Lyra et al. (2018) Lyra W., Raettig N., Klahr H., 2018, Research Notes of the American Astronomical Society, 2, 195
  * Marsh & Mahoney (1992) Marsh K. A., Mahoney M. J., 1992, ApJ, 395, L115
  * Marsh & Mahoney (1993) Marsh K. A., Mahoney M. J., 1993, ApJ, 405, L71
  * Miranda & Lai (2015) Miranda R., Lai D., 2015, MNRAS, 452, 2396
  * Miranda et al. (2017) Miranda R., Muñoz D. J., Lai D., 2017, MNRAS, 466, 1170
  * Muñoz & Lai (2016) Muñoz D. J., Lai D., 2016, ApJ, 827, 43
  * Pérez et al. (2018) Pérez S., Casassus S., Benítez- Llambay P., 2018, MNRAS, 480, L12
  * Pinte et al. (2006) Pinte C., Ménard F., Duchêne G., Bastien P., 2006, A&A, 459, 797
  * Pinte et al. (2009) Pinte C., Harries T. J., Min M., Watson A. M., Dullemond C. P., Woitke P., Ménard F., Durán-Rojas M. C., 2009, A&A, 498, 967
  * Price (2007) Price D. J., 2007, Publ. Astron. Soc. Australia, 24, 159
  * Price & Laibe (2015) Price D. J., Laibe G., 2015, MNRAS, 451, 813
  * Price et al. (2018a) Price D. J., et al., 2018a, Publications of the Astronomical Society of Australia, 35, e031
  * Price et al. (2018b) Price D. J., et al., 2018b, MNRAS, 477, 1270
  * Ragusa et al. (2016) Ragusa E., Lodato G., Price D. J., 2016, MNRAS, 460, 1243
  * Ragusa et al. (2017) Ragusa E., Dipierro G., Lodato G., Laibe G., Price D. J., 2017, MNRAS, 464, 1449
  * Ratzka et al. (2005) Ratzka T., Köhler R., Leinert C., 2005, Astronomy and Astrophysics, 437, 611
  * Regály et al. (2012) Regály Z., Juhász A., Sándor Z., Dullemond C. P., 2012, MNRAS, 419, 1701
  * Ruge et al. (2016) Ruge J. P., Flock M., Wolf S., Dzyurkevich N., Fromang S., Henning T., Klahr H., Meheut H., 2016, A&A, 590, A17
  * Schworer et al. (2017) Schworer G., et al., 2017, The Astrophysical Journal, 842, 77
  * Simon et al. (1995) Simon M., et al., 1995, ApJ, 443, 625
  * Strom et al. (1989) Strom K. M., Strom S. E., Edwards S., Cabrit S., Skrutskie M. F., 1989, AJ, 97, 1451
  * Tang et al. (2017) Tang Y.-W., et al., 2017, ApJ, 840, 32
  * Testi et al. (2014) Testi L., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 339 (arXiv:1402.1354), doi:10.2458/azu_uapress_9780816531240-ch015
  * Thun et al. (2017) Thun D., Kley W., Picogna G., 2017, A&A, 604, A102
  * Turner et al. (2014) Turner N. J., Fromang S., Gammie C., Klahr H., Lesur G., Wardle M., Bai X.-N., 2014, Protostars and Planets VI, pp 411–432
  * Wright et al. (2015) Wright C. M., et al., 2015, MNRAS, 453, 414
  * Wu & Sheehan (2017) Wu Y.-L., Sheehan P. D., 2017, ApJ, 846, L26
  * Zhu et al. (2014) Zhu Z., Stone J. M., Rafikov R. R., Bai X.-n., 2014, ApJ, 785, 122
  * van der Marel et al. (2013) van der Marel N., et al., 2013, Science, 340, 1199
  * van der Marel et al. (2015) van der Marel N., Pinilla P., Tobin J., van Kempen T., Andrews S., Ricci L., Birnstiel T., 2015, ApJ, 810, L7
  * van der Marel et al. (2016a) van der Marel N., van Dishoeck E. F., Bruderer S., Andrews S. M., Pontoppidan K. M., Herczeg G. J., van Kempen T., Miotello A., 2016a, A&A, 585, A58
  * van der Marel et al. (2016b) van der Marel N., Cazzoletti P., Pinilla P., Garufi A., 2016b, ApJ, 832, 178


# MDF: Magnetic Particle Imaging Data Format

T. Knopp1,2, T. Viereck3, G. Bringout4, M. Ahlborg5, A. von Gladiss5, C.
Kaethner5, A. Neumann5, P. Vogel6, J. Rahmer7, M. Möddel1,2

1Section for Biomedical Imaging, University Medical Center Hamburg-Eppendorf,
Germany
2Institute for Biomedical Imaging, Hamburg University of Technology, Germany
3Institute of Electrical Measurement and Fundamental Electrical Engineering,
TU Braunschweig, Germany
4Physikalisch-Technische Bundesanstalt, Berlin, Germany
5Institute of Medical Engineering, University of Lübeck, Germany
6Department of Experimental Physics 5 (Biophysics), University of Würzburg,
Germany
7Philips GmbH Innovative Technologies, Research Laboratories, Hamburg, Germany

###### Abstract

Magnetic particle imaging (MPI) is a tomographic method to determine the
spatio-temporal distribution of magnetic nanoparticles. In this document, a
file format for the standardized storage of MPI and magnetic particle
spectroscopy (MPS) data is introduced. The aim of the Magnetic Particle
Imaging Data Format (MDF) is to provide a coherent way of exchanging MPI and
MPS data acquired with different devices worldwide. The focus of the MDF is on
sequence parameters, measurement data, calibration data, and reconstruction
data. The format is based on the hierarchical document format in version 5
(HDF5). This document discusses the MDF version 2.0.0, which is not backward
compatible with version 1.x.y.

Version 2.0.0

##  1 Introduction

The purpose of this document is to introduce a file format for exchanging
Magnetic Particle Imaging (MPI) and Magnetic Particle Spectroscopy (MPS) data.
The Magnetic Particle Imaging Data Format (MDF) is based on the hierarchical
document format (HDF) in version 5 [1]. HDF5 is able to store multiple
datasets within a single file providing a powerful and flexible data
container. To allow an easy exchange of MPI data, one has to specify a naming
scheme within HDF5 files which is the purpose of this document. In order to
create and access HDF5 data, an Open Source C library is available that
provides dynamic access from most programming languages. MATLAB supports HDF5
by its functions h5read and h5write. For Python, the h5py package exists, and
the Julia programming language provides access to HDF5 files via the HDF5
package. For languages based on the .NET framework, the HDF5DotNet library is
available.

The MDF is mainly focused on storing measurement data, calibration data, or
reconstruction data together with the corresponding sequence parameters and
metadata. Even though it is possible to combine measurement data and
reconstruction data into a single file, it is recommended to use a single file
for each of the following dataset types:

  1. 1.

Measurement data

  2. 2.

Calibration data

  3. 3.

Reconstruction data

###  1.1 Datatypes & Storage Order

MPI parameters are stored as regular HDF5 datasets. HDF5 attributes are not
used in the current specification of the MDF. For most datasets, a fixed
datatype is used, for example the drive-field amplitudes are stored as
H5T_NATIVE_DOUBLE values. Whenever a data type has a big- and little-endian
version, the little-endian data type should be used. For our convenience, we
refer to the HDF5 datatypes H5T_STRING, H5T_NATIVE_DOUBLE, and
H5T_NATIVE_INT64 as String, Float64 and Int64. Boolean data is stored as
H5T_NATIVE_INT8, which we refer to as Int8.

The datatype of the MPI measurement and calibration data offers more freedom
and is denoted by Number, which can be any of the following HDF5 data types:
H5T_NATIVE_FLOAT, H5T_NATIVE_DOUBLE, H5T_NATIVE_INT8, H5T_NATIVE_INT16,
H5T_NATIVE_INT32, H5T_NATIVE_INT64 or a complex number as defined below.

Since storing complex data in HDF5 is not standardized, we define a compound
datatype H5T_COMPOUND in HDF5 with fields r and i using one of the above
mentioned floating point types to represent the real and the imaginary part of
a complex number. This representation was chosen because it is also the
default behavior for complex numbers in Python using numpy and h5py.
Complex128 refers to the complex compound type with base type
H5T_NATIVE_DOUBLE.

For later identification of a data set, we store three Universally Unique
Identifiers (UUIDs) (RFC 4122) [2] in its canonical textual representation as
32 hexadecimal digits, displayed in five groups separated by hyphens
8-4-4-4-12 as for example ee94cb6d-febf-47d9-bec9-e3afa59bfaf8. For the
generation of the UUIDs, we recommend to use version 4 of the UUID
specification.

Whenever multidimensional data is stored dimensions are arranged in a way that
cache is utilized best for fast reconstruction or fast frame selection for
example. The leading dimension in the MDF specification is slowest to access
and the last dimension is the fastest to access (contiguous memory access).
That also means that dependent on the memory layout of your programming
language the order of the dimensions will be just as in the MDF specifications
(row major) or reversed (column major). Please take this into account when
reading or writing MDF files, which includes the usage of HDF5 viewers.

###  1.2 Units

With one exception, physical quantities are given in SI units. The magnetic
field strength is reported in T$\mu_{0}^{-1}=4\pi$Am${}^{-1}\mu_{0}^{-1}$.
This convention has been proposed in the first MPI publication [3] and
consistently been used in most MPI related publications. The aim of this
convention is to report the numbers on a Tesla scale, which most readers with
a background in magnetic resonance imaging are familiar with.

###  1.3 Optional Parameters

The MDF has 8 main groups in the root directory and 2 sub-groups. We
distinguish between optional and non-optional groups as well as optional, non-
optional, and conditional parameters.

Any optional parameter can be omitted, whereas any non-optional parameter in a
non-optional group is mandatory. Conditional parameters are linked to Boolean
parameters and have to be provided, if these parameters are true and can be
omitted if they are false. If a parameter is optional, non-optional, or
conditional is indicated by yes, no, or the corresponding Boolean parameter.

If a group is optional all of its parameters may be omitted, if this group is
not used. The groups /, /study, /experiment, /scanner, /acquisition contain
mostly metadata and are mandatory. The /tracer group is only mandatory if
magnetic material has been placed in the MPI system. The groups /measurement,
/calibration, and /reconstruction are all optional. In case of calibration
measurements, the /calibration group is mandatory. The reconstruction data is
stored in /reconstruction.

###  1.4 User defined parameters

It is possible to store user defined parameters in an MDF. For instance one
may want to store the temperature of the room in which your MPI device is
operated. In this case, one is free to add new parameters to any of the
existing groups. Moreover, if necessary, one can also introduce new groups. In
order to be able to distinguish these datasets and groups from the specified
ones, it is mandatory to use the prefix “_” for all parameters and groups. As
an example, one could add a new group /_room that includes the dataset
_temperature. Using the prefix “_” will ensure that the stored dataset is
compatible with future versions of the MDF.

###  1.5 Naming Convention

Several parameters within the MDF are linked in dimensionality. We use short
variable names to indicate these connections. The following table describes
the meaning of each variable name used in this specification.  Variable Number
of $A$ tracer materials/injections for multi-color MPI $N$ acquired frames,
same as a spatial position for calibration $O$ acquired frames w/o background
frames $J$ periods within one frame $Y$ partitions of each patch position $C$
receive channels $D$ drive-field channels $F$ frequencies describing the
drive-field waveform $V$ points sampled at receiver during one drive-field
period $W$ sampling points containing processed data ($W=V$ if no frequency
selection or bandwidth reduction has been applied) $K$ frequencies describing
the processed data ($K=V/2+1$ if no frequency selection or bandwidth reduction
has been applied) $Q$ frames in the reconstructed MPI data set $P$ voxels in
the reconstructed MPI data set $S$ channels in the reconstructed MPI data set

###  1.6 Contact

If you find mistakes in this document or the specified file format or if you
want to discuss extensions or improvements to this specification, please open
an issue on GitHub:
https://github.com/MagneticParticleImaging/MDF
As the file format is versionized it will be possible to extend it for future
needs of MPI. The current version discussed in this document is version 2.0.0.

###  1.7 arXiv

As of version 1.0.1, the most recent release of these specifications can also
be also found at:
https://arxiv.org/abs/1602.06072
If you use MDF, please cite us using the arXiv reference, which is also
available for download as MDF.bib from GitHub.

###  1.8 Code examples

If you want to get a basic impression of how to handle MDF files you can visit
the gitub repository of the MDF project:
https://github.com/MagneticParticleImaging/MDF
There you will find the directory example, which contains a code example
written in Julia, MATLAB, and Python. More details can be found in the README
of the repository.

###  1.9 Reference Implementation

A reference implementation for a high level MDF access is available at:
https://github.com/MagneticParticleImaging/MPIFiles.jl
MPIFiles.jl is a package written in the programming language Julia [4, 5, 6].
It can read MDF V1, MDF V2, and the dataformat of Bruker MPI systems using a
common interface. It also provides functions to convert MDF V1 or Bruker MPI
files into MDF V2.

##  2 Data (group: /)

##### Remarks:

Within the root group, the metadata about the file itself is stored. Within
several subgroups, the metadata about the experimental setting, the MPI
tracer, and the MPI scanner can be provided. The actual data is stored in
dedicated groups about measurement data and reconstruction data.

Parameter | Type | Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
version | String | 1 | 2.0.0 | no |  Version of the file format
uuid | String | 1 | 3170fdf8-f8e1-4cbf-ac73-41520b41f6ee | no |  Universally Unique Identifier (RFC 4122) of MDF file
time | String | 1 | yyyy-mm-ddThh:mm:ss.ms | no |  UTC creation time of MDF data set

###  2.1 Study Description (group: /study/, non-optional)

##### Remarks:

A study is supposed to group a series of experiments to support, refute, or
validate a hypothesis. The study group should contain name, number, uuid, and
description of the study.

Parameter | Type | Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
name | String | 1 |  | no |  Name of the study
number | Int64 | 1 |  | no |  Number of the study
uuid | String | 1 | 295258fe-b650-4e5f-96db-b83f11089a6c | no |  Universally Unique Identifier (RFC 4122) of study
description | String | 1 |  | no |  Short description of the study

###  2.2 Experiment Description (group: /experiment/, non-optional)

##### Remarks:

For each experiment within a study a name, number, uuid, and description have
to be provided. Additionally, the name of the subject imaged and the flag
isSimulation indicating if data has been obtained via simulations have to be
stored.

Parameter | Type | Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
name | String | 1 |  | no |  Experiment name
number | Int64 | 1 |  | no |  Experiment number within study
uuid | String | 1 | f96dbc48-1ebd-44c7-b04d-1b45da054693 | no |  Universally Unique Identifier (RFC 4122) of experiment
description | String | 1 |  | no |  Short description of the experiment
subject | String | 1 |  | no |  Name of the subject that was imaged
isSimulation | Int8 | 1 |  | no |  Flag indicating if the data in this file is simulated rather than measured

###  2.3 Tracer Parameters (group: /tracer/, optional)

##### Remarks:

The tracer parameter group contains information about the MPI tracers used
during the experiment. For each of the $A$ tracers name, batch, vendor,
volume, and molar concentration of solute per liter must be provided.
Additionally, the time point of injection can be stored.

This version of the MDF can handle two basic scenarios. In the first one,
static tracer phantoms are used. In this case, the phantom contains $A$
distinct tracers. For example, these might be particles of different core
sizes, mobile or immobilized particles. In this case, injectionTime is not
used. In the second case, $A$ boli (e.g. pulsed boli) are administrated during
the measurement, in which case the approximate administration volume, tracer
type and time point of injection can be provided. Note that the injection
clock recording the injection time should be synchronized with the clock,
which provides the starting time of the measurement.

In case of a background measurement with no applied tracers in the scanner,
the optional tracer group can be omitted.

Parameter | Type | Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
name | String | $A$ |  | no |  Name of tracer used in experiment
batch | String | $A$ |  | no |  Batch of tracer
vendor | String | $A$ |  | no |  Name of tracer supplier
volume | Float64 | $A$ | L | no |  Total volume of applied tracer
concentration | Float64 | $A$ | mol(solute)/L | no |  Molar concentration of solute per litre
solute | String | $A$ |  | no |  Solute, e.g. Fe
injectionTime | String | $A$ | yyyy-mm-ddThh:mm:ss.ms | yes |  UTC time at which tracer injection started

###  2.4 Scanner Parameters (group: /scanner/, non-optional)

##### Remarks:

The scanner parameter group contains information about the MPI scanner used,
such as name, manufacturer, boreSize, field topology, facility where the
scanner is installed, and the operator.

Parameter | Type | Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
name | String | 1 |  | no |  Scanner name
facility | String | 1 |  | no |  Facility where the MPI scanner is installed
operator | String | 1 |  | no |  User who operates the MPI scanner
manufacturer | String | 1 |  | no |  Scanner manufacturer
topology | String | 1 |  | no |  Scanner topology (e.g. FFP, FFL, MPS)
boreSize | Float64 | 1 | m | yes |  Diameter of the bore

###  2.5 Acquisition Parameters (group: /acquisition/, non-optional)

##### Remarks:

The acquisition parameter group can describe different imaging protocols and
trajectory settings. The corresponding data is organized into general
information, a subgroup containing information on the $D$ excitation channels,
and a subgroup containing information on the $C$ receive channels.

The term frame refers to the data collected during one acquisition period.
Usually all data within a frame will be combined to reconstruct a single MPI
image/tomogram. In the simplest scenario, this data is acquired during one
drivefield/cycle. If block-averaging is applied the amount of data captured
stays the same while the acquisition time increases by a factor of
numAverages. In this document we will refer to the product of one
drivefield/cycle and numAverages as one drive-field period.

In case that the static selection field is changed over time one will measure
several drive-field periods each with a different gradient field setting. In a
stepped multi-patch sequence the gradient field will be spatially shifted
and/or scaled in between drive-field periods and remain constant within. In
such a scenario a full frame consists of numPeriodsPerFrame ($J$) periods. The
gradient scaling or shift within each period are described by the fields
/acquisition/gradient and /acquisition/offsetField. The former has $J\times
Y\times 3\times 3$ entries whereas the later is of dimension $J\times Y\times
3$. While $J$ labels the drive-field periods, $Y$ can be used to describe
dynamic multi-patch sequences where the gradient and offset field change
during a drive-field period. Each period is discretized into $Y$ equidistantly
spaced time intervals. The values in /acquisition/gradient and
/acquisition/offsetField describe the gradient and offset field at the
beginning of each time interval. In case $Y=1$ the fields at the beginning of
the period $J$ are provided, which can be used to describe a static multi-
patch sequence. Higher numbers of $Y$ allow the description of an arbitrarily
fine grained gradient and offset field sequence.

It should be noted that the combined gradient and offset field can be
expressed locally around the scanner center
$\mathbf{r}_{\text{c}}\in\mathbb{R}^{3}$ via the Taylor expansion

$\mathbf{H}(\mathbf{r})=\mathbf{H}(\mathbf{r}_{\text{c}})+\mathbf{J}_{\mathbf{H}}(\mathbf{r}_{\text{c}})(\mathbf{r}-\mathbf{r}_{\text{c}})+\dots,$
(1)

where $\mathbf{J}_{\mathbf{H}}(\mathbf{r}_{\text{c}})$ is the Jacobian of
$\mathbf{H}$ at $\mathbf{r}_{\text{c}}$. /acquisition/offsetField stores the
field $\mathbf{H}(\mathbf{r}_{\text{c}})$ and /acquisition/gradient stores the
$3\times 3$ matrix $\mathbf{J}_{\mathbf{H}}$. Both parameters are for most
scanner topologies sufficient to uniquely describe the applied field sequence.
In case of non-linear fields and/or a higher degree of freedom in the number
of applied field generators, one may need additional custom fields in order to
uniquely define the sequence parameters.

Parameter | Type | Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
startTime | String | 1 | yyyy-mm-ddThh:mm:ss.ms | no |  UTC start time of MPI measurement
numAverages | Int64 | 1 | 1 | no |  Number of block averages per drive-field period.
numPeriodsPerFrame | Int64 | 1 | 1 | no |  Number of drive-field periods within a frame denoted by $J$
numFrames | Int64 | 1 | 1 | no |  Number of available measurement frames $N$
gradient | Float64 | $J\times Y\times 3\times 3$ | Tm${}^{-1}\mu_{0}^{-1}$ | yes |  Gradient strength of the selection field in $x$, $y$, and $z$ directions
offsetField | Float64 | $J\times Y\times 3$ | T$\mu_{0}^{-1}$ | yes |  Offset field applied

####  2.5.1 Drive Field Parameters (group: /acquisition/drivefield/, non-
optional)

##### Remarks:

The drive field subgroup describes the excitation details of the imaging
protocol. On the lowest level, each MPI scanner contains $D$ channels for
excitation. Since most drive-field parameters may change from period to
period, they have a leading dimension $J$.

These excitation signals are usually sinusoidal and can be described by $D$
amplitudes (drive field strengths), $D$ phases, a base frequency, and $D$
dividers. In a more general setting, the generated drive field of channel $d$
can be described by

$H_{d}(t)=\sum_{l=1}^{F}A_{l}\Lambda_{l}(2\pi f_{l}t+\varphi_{l})$

where $F$ is the number of frequencies on the channel, $A_{l}$ is the drive-
field strength, $\phi_{l}$ is the phase, $f_{l}$ is the frequency
(baseFrequency/dividerl), and $\Lambda_{l}$ is the waveform. The waveform is
specified by a dedicated parameter waveform. It can be set to sine, triangle,
or custom.

Parameter | Type | Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
numChannels | Int64 | 1 | 1 | no |  Number of drive field channels, denoted by $D$
strength | Float64 | $J\times D\times F$ | T$\mu_{0}^{-1}$ | no |  Applied drive field strength
phase | Float64 | $J\times D\times F$ | rad, $[-\pi,\pi)$ | no |  Applied drive field phase $\varphi$
baseFrequency | Float64 | 1 | Hz | no |  Base frequency to derive drive field frequencies
divider | Int64 | $D\times F$ | 1 | no |  Divider of the baseFrequency to determine the drive field frequencies
waveform | String | $D\times F$ | 1 | no |  Waveform type: sine, triangle or custom
cycle | Float64 | 1 | s | no |  Trajectory cycle is determined by lcm(divider)/baseFrequency. It will not change when averaging was applied. The duration for measuring the $V$ data points (i.e. the drive-field period) is given by the product of period and numAverages

####  2.5.2 Receiver (group: /acquisition/receiver/, non-optional)

##### Remarks:

The receiver subgroup describes the details on the MPI receiver. For a multi-
patch sequence, it is assumed, that the signal acquisition only takes place
during particle excitation. During each drive-field cycle, $C$ receive
channels record some quantity related to the magnetization dynamic. In most
cases these, will be continuous voltage signals induced into the $C$ receive
coils, which are proportional to the change of the particle magnetization.
This analog signal will usually be converted by some sort of analog to digital
converter (ADC) to a discrete series of integer numbers $r_{ci}$ for each
channel $c$. To map the these values to the MPI measurement signal
$u_{ci}^{ADC}$, one has to scale the numbers $r_{ci}$ and add an offset factor

$u_{ci}^{ADC}=a_{c}r_{ci}+b_{c}.$

Here, $a_{c}$ and $b_{c}$ are the scaling and offset factors corresponding to
channel $c$, which can be stored in dataConversionFactor. In case the
conversion was already performed the dataConversionFactor can be ommited.

The MPI measurement signal is acquired at $V$ equidistant time points. For
inductive measurement systems the signal is usually not measured directly at
the receive coils but amplified and filtered first, which may damp and distort
the signal. Therefore, a transfer function can be stored in the parameter
transferFunction, which relates the Fourier domain voltage induced at the
receive coil $\hat{u}_{k}^{\text{coil}}$ and the Fourier domain voltage
$\hat{u}_{k}^{\text{ADC}}$ measured at the ADC by

$\hat{u}_{k}^{\text{ADC}}=a_{k}\hat{u}_{k}^{\text{coil}},\quad k=1,\dots,K.$

Here, $a_{k}$ are the unitless parameters stored in transferFunction for each
receive channel individually.

For MPS systems one can additionally store a parameter that maps the induced
voltage to the mean magnetic moment of a magnetic nanoparticle located at the
center of the scanner. More precisely, in each receive coil a projection of
the mean magnetic moment onto the coil sensitivity is measured. The relation
of this projection and the voltage at the receive coil in frequency space
representation is given by

$\displaystyle\hat{u}_{k}^{\text{coil}}=2\pi\textrm{i}k\beta\hat{m}_{k}^{\text{proj}},\quad
k=1,\dots,K.$

where $\hat{m}_{k}^{\text{proj}}$ is the orthogonal projection of the magnetic
moment onto the coil sensitivity at the scanner center and $\beta$ is the
channel dependent conversion factor that is stored in the parameter
inductionFactor.

Parameter | Type | Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
numChannels | Int64 | 1 |  | no |  Number of receive channels $C$
bandwidth | Float64 | $1$ | Hz | no |  Bandwidth of the receiver unit
numSamplingPoints | Int64 | $1$ |  | no |  Number of sampling points during one period, denoted by $V$
unit | String | $1$ |  | no |  SI unit of the measured quantity, usually Voltage V
dataConversionFactor | Float64 | $C\times 2$ | unit | yes |  Dimension less scaling factor and offset $(a_{c},b_{c})$ to convert raw data into a physical quantity with corresponding unit of measurement unit
transferFunction | Complex128 | $C\times K$ |  | yes |  Transfer function of the receive channels in Fourier domain. unit is the field from the /measurement group
inductionFactor | Float64 | $C$ |  unit A-1m-2 | yes |  Induction factor mapping the projection of the magnetic moment to the voltage in the receive coil.

###  2.6 Measurement (group: /measurement/, optional)

##### Remarks:

MPI data is usually acquired by a series of measurements and optional
background measurements. Here, we refer to background measurements as MPI data
captured, when any signal generating material, e.g. a phantom or a delta
sample is removed from the scanner bore. Initially, all data is available in
time domain, where the data of a single frame consists of the signal recorded
for all periods in each receive channel, i.e. $J\times C\times V$ data points
per set with the temporal index being the fastest to access. If several
measurements are acquired (indicated by numFrames), the frame dimension is the
slowest to access. Along this dimension, the frames are ordered with respect
to the time at which they were acquired starting with the measurement acquired
first and stopping with the measurement acquired last. We refer to this data
as raw measurement data. In Fourier representation, each frame would be stored
by $J\times C\times K$ complex data points and $K=V/2+1$.

Often it is not convenient to store the raw data but to perform certain
processing steps and store the processed data. These steps may lead to a
reduction of the number of sampling points from $V$ to $W$ or a corresponding
reduction of frequency components $K$ depending on the final representation in
which the raw/processed data is stored. The most common processing steps are:

  1. 1.

Spectral leakage correction, which may be applied to ensure that each
individual frame is periodic.

  2. 2.

Background correction, where the background signal in subtracted.

  3. 3.

Fourier transformation bringing the data from time into the Fourier
representation an storing them in Fourier representation.

  4. 4.

Transfer function correction to obtain the magnetic moment or induced voltage
that has been measured.

  5. 5.

Frequency selection to reduce the number of frequency components, e.g.
bandwidth reduction or selection of high signal frequency components.

  6. 6.

Dimension permutation, which is usually applied to Fourier transformed data
exchanging the storing order of the data for fast access to the frames.

  7. 7.

Frame permutation to reorder the frames within the data set.

For each of the steps above there is a corresponding flag within this group
indicating if the corresponding processing step has been carried out.

During processing one might want to keep track which of the final $N$ frames
belong to background measurements and which do not. Therefore, the binary mask
isBackgroundFrame should be used. If frequency selection has been performed,
frequencySelection stores the $K$ frequency components (subset) selected from
the set of acquired frequency components. If performed, a frame permutation
can be described by a bijective mapping
$\sigma:\left\\{1,2,\dots,N\right\\}\rightarrow\left\\{1,2,\dots,N\right\\}$
of the set of frame indices to itself. If such a permutation is performed,
$\sigma$ is stored in the one-line notation as $\sigma(1)$, $\sigma(2)$,
$\dots$, $\sigma(N)$ in framePermutation.

Parameter | Type |  Dim | Optional |  Description
---|---|---|---|---
data | Number |  $N\times J\times C\times K$ or $J\times C\times K\times N$ or $N\times J\times C\times W$ or $J\times C\times W\times N$ | no |  Measured data at a specific processing stage
isBackgroundFrame | Int8 |  $N$ | no |  Mask indicating for each of the $N$ frames if it is a background measurement (true) or not
isSpectralLeakageCorrected | Int8 |  1 | no |  Flag, if spectral leakage correction has been applied
isBackgroundCorrected | Int8 |  1 | no |  Flag, if the background has been subtracted
isFourierTransformed | Int8 |  1 | no |  Flag, if the data is stored in frequency space
isTransferFunctionCorrected | Int8 |  1 | no |  Flag, if the data has been corrected by the transferFunction
isFrequencySelection | Int8 |  1 | no |  Flag, if only a subset of frequencies has been selected and stored, see frequencySelection
isFastFrameAxis | Int8 |  1 | no |  Flag, if the frame dimension $N$ has been moved to the last dimension
isFramePermution | Int8 |  1 | no |  Flag, if the order of frames has been changed, see framePermutation
frequencySelection | Int64 |  K | isFrequencySelection |  Indices of selected frequency components
framePermutation | Int64 |  $N$ | isFramePermutation |  Indices of original frame order

###  2.7 Calibration (group: /calibration/, optional)

##### Remarks:

The calibration group describes a calibration measurement (system matrix),
although it does not hold the data itself. Each of the calibration
measurements is taken with a calibration sample (delta sample) at a grid
centered position inside the FOV of the device. If available, background
measurements are taken with the delta sample outside of the FOV of the
scanner. Usually, not the raw measurements are stored but processed data,
where at least averaging, Fourier transformation, frame permutation and
transposition has been performed yielding a total of $N$ processed frames. Of
these frames $O$ frames correspond to the calibration scans at the $O$ spatial
positions, whereas the other frames correspond to background measurements
taken throughout the calibration process. All processing steps are documented
in the /measurement group.

If the measurements were taken on a regular grid of size $N_{x}\times
N_{y}\times N_{z}$, the permutation is usually done such that measurements are
ordered with respect to their $x$ position first, second with respect to their
$y$ position, and last with respect to their $z$ position. Background
measurements are collected at the end in /measurement/data, which in
combination with reordering of the measurements allows a fast access to the
system matrix. If a different storage order is used this can be documented
using the optional parameter order. For non-regular sampling points, there is
the possibility to explicitly store all $O$ positions. If the calibration
measurement is performed in a MPS system, the spatial positions are usually
emulated by applying offset fields, which can be stored in offsetFields.

Parameter | Type |  Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
method | String |  1 |  | no |  Method used to obtain calibration data. Can for instance be robot, hybrid, or simulation
size | Int64 |  $3$ |  | yes |  Number of voxels in each dimension, inner product is $O$
order | String |  1 |  | yes |  Ordering of the dimensions, default is xyz
positions | Float64 |  $O\times 3$ | m | yes |  Position of each of the grid points, stored as ($x$, $y$, $z$) triples
offsetFields | Float64 |  $O\times 3$ | T$\mu_{0}^{-1}$ | yes |  Applied offset field strength to emulate a spatial position ($x$, $y$, $z$)
deltaSampleSize | Float64 |  3 | m | yes |  Size of the delta sample used for calibration scan
fieldOfView | Float64 |  $3$ | m | yes |  Field of view of the system matrix
fieldOfViewCenter | Float64 |  $3$ | m | yes |  Center of the system matrix (relative to origin/center)
snr | Float64 |  $J\times C\times K$ |  | yes |  Signal-to-noise estimate for recorded frequency components

###  2.8 Reconstruction Results (group: /reconstruction/, optional)

Reconstruction results are stored using the parameter data inside this group.
data contains a $Q\times P\times S$ array, where $Q$ denotes the number of
reconstructed frames within the data set, $P$ denotes the number of voxels and
$S$ the number of multispectral channels. If no multispectral reconstruction
is performed, then one may set $S=1$. Depending on the reconstruction the grid
of the reconstruction data can be different from the system matrix grid.
Hence, grid parameters are mirrored in the /reconstruction group.

For analysis of the MPI tomograms, it is often required to know which parts of
the reconstructed tomogram have been covered by the trajectory of the field
free region. In MPI, one refers to the non-covered region as overscan region.
Therefore, the optional binary field isOverscanRegion stores for each voxel if
it is part of the overscan region. If no voxel lies within the overscan
region, isOverscanRegion may be omitted.

Parameter | Type | Dim | Unit/Format | Optional |  Description
---|---|---|---|---|---
data | Number | $Q\times P\times S$ |  | no |  Reconstructed data
size | Int64 | $3$ |  | yes |  Number of voxels in each dimension, inner product is $P$
order | String | 1 |  | yes |  Ordering of the dimensions, default is xyz
positions | Float64 | $P\times 3$ | m | yes |  Position of each of the grid points, stored as ($x$, $y$, $z$) tripels
fieldOfView | Float64 | $3$ | m | yes |  Field of view of reconstructed data
fieldOfViewCenter | Float64 | $3$ | m | yes |  Center of the reconstructed data (relative to scanner origin/center)
isOverscanRegion | Int8 | $P$ |  | yes |  Mask indicating for each of the $P$ voxels if it is part of the overscan region (true) or not

##  3 Changelog

###  3.1 v2.0.0

  * •

Version 2 of the MDF is a major update breaking backwards compatibility with
v1.x. The major update was necessary due to several shortcomings in the v1.x.

  * •

The naming of parameters was made more consistent. Furthermore, some
parameters moved from one group into another.

  * •

Defined a complex datatype using a HDF5 compound type.

  * •

In v1.x it was not possible to store background data. This functionality has
been added in v2.

  * •

We simplified the measurement group and made it much more expressive. In v1.x
it was not entirely clear, which processing steps have been applied to the
measurement data in the stored dataset. The measurement group now contains
several flags that precisely document the state of the stored data. Using this
it is now possible to also store calibration data in the measurement group.
The calibration group in turn only stores metadata about a calibration
experiment while the actual data is store in the measurement group.

  * •

Updated Affiliations in the MDF specification.

  * •

Improved the general descriptions of fields and groups.

  * •

In v1.x the MDF allowed many fields to have varying dimensions depending on
the context. As of version 2.0.0 only one field offers this freedom. This
change should make implementations handling MDF files less complex.

  * •

Number has been introduced as a generic type.

  * •

Added a table listing all variable names used in the descriptions of
parameters.

  * •

Added a section describing the possibility to add user defined parameters to
MDF files.

  * •

Added a description for optional and non-optional groups and conditional,
optional, and non-optional datasets.

  * •

Added a short section on the code examples on the Github repository.

  * •

Support for triangle wave forms has been added.

  * •

Support for multiple excitation frequencies on a drive-field channel has been
added.

  * •

Added the dimension $A$ to all fields of the tracer group to be able to
describe settings where multiple tracers are used or tracers are administered
multiple times.

  * •

Added the possibility to store the tracer concentration also for non iron
based tracer materials by adding the /tracer/solute field and redefining the
field tracer/concentration.

  * •

Improved documentation for the storage of ADC transfer functions. It is also
possible now to store the measurement data as integer data and use a
dataConversionFactor to describe the mapping to a physical representation
(e.g. Volt)

  * •

Added support for receive coil to ADC transfer functions.

  * •

Added support for mean magnetic moment to receive coil voltage transfer
functions.

  * •

Split the /study group into the /experiment group and the study group. This
allows to provide more fine grained information on study and experiment.

  * •

Added possibility to mark the overscan region.

  * •

Added new section changelog to the MDF documentation to record the development
of the MDF.

  * •

Updated README.md and MDF.bib in the github repository.

  * •

Updated code examples in the github repository.

  * •

Added a section on the MDF reference implementation MPIFiles.jl. Since sanity
checks will be covered by this package, the description on sanity checks has
been removed.

###  3.2 v1.0.5

  * •

Added the possibility to store different channels of reconstructed data.

  * •

Added support for receive channels with different characteristics (e.g.
bandwidth).

  * •

Made dataset /acquisition/receiver/frequencies optional.

  * •

Extended the description on the data types, which are used to store data.

  * •

Added references for Julia and HDF5 to the specifications.

###  3.3 v1.0.4

  * •

Clarify that HDF5 datasets are used to store MPI parameters.

###  3.4 v1.0.3

  * •

Updated Affiliations in the MDF specification.

  * •

Included data download into the Python and Matlab example code.

  * •

Changes in the Python and Matlab example code to be better comparable to the
Julia example code.

###  3.5 v1.0.2

  * •

Added reference to arXiv paper and bibtex file for reference.

###  3.6 v1.0.1

  * •

A sanity check within the Julia code shipped alongside the specifications.

  * •

An update to the specification documenting the availability of a sanity check.

  * •

Updated MDF files on https://www.tuhh.de/ibi/research/mpi-data-format.html.

  * •

Updated documentation to the Julia, Matlab and Python reconstruction scripts.

  * •

Improved Julia reconstruction script, automatically downloading the required
MDF files.

## References

  * [1] The HDF Group.  Hierarchical Data Format, version 5, 1997-2016.  http://www.hdfgroup.org/HDF5/.
  * [2] Paul J Leach, Michael Mealling, and Rich Salz.  A universally unique identifier (uuid) urn namespace.  2005\.
  * [3] Bernhard Gleich and Jürgen Weizenecker.  Tomographic imaging using the nonlinear response of magnetic particles.  Nature, 435(7046):1214–1217, 2005.
  * [4] Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and Alan Edelman.  Julia: A fast dynamic language for technical computing.  CoRR, abs/1209.5145, 2012.
  * [5] Jeff Bezanson, Jiahao Chen, Stefan Karpinski, Viral B. Shah, and Alan Edelman.  Array operators using multiple dispatch: a design methodology for array implementations in dynamic languages.  CoRR, abs/1407.3845, 2014.
  * [6] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B. Shah.  Julia: A fresh approach to numerical computing.  CoRR, abs/1411.1607, 2014.