# Testing photoevaporation and MHD disk wind models through future high-
angular resolution radio observations: the case of TW Hydrae

Luca Ricci  Department of Physics and Astronomy, California State University
Northridge, 18111 Nordhoff Street, Northridge, CA 91330, USA  Sarah K. Harter
Department of Physics and Astronomy, California State University Northridge,
18111 Nordhoff Street, Northridge, CA 91330, USA  Barbara Ercolano  University
Observatory, Faculty of Physics, Ludwig-Maximilians-Universität München,
Scheinerstr. 1, 81679 Munich, Germany  Michael Weber  University Observatory,
Faculty of Physics, Ludwig-Maximilians-Universität München, Scheinerstr. 1,
81679 Munich, Germany

###### Abstract

We present theoretical predictions for the free-free emission at cm
wavelengths obtained from photoevaporation and MHD wind disk models adjusted
to the case of the TW Hydrae young stellar object. For this system, disk
photoevaporation with heating due to the high-energy photons from the star has
been proposed as a possible mechanism to open the gap observed in the dust
emission with ALMA. We show that the photoevaporation disk model predicts a
radial profile for the free-free emission that is made of two main spatial
components, one originated from the bound disk atmosphere at 0.5-1 au from the
star, and another more extended component from the photoevaporative wind at
larger disk radii. We also show that the stellar X-ray luminosity has a
significant impact on both these components. The predicted radio emission from
the MHD wind model has a smoother radial distribution which extends to closer
distances to the star than the photoevaporation case. We also show that a
future radio telescope such as the Next Generation Very Large Array (ngVLA)
would have enough sensitivity and angular resolution to spatially resolve the
main structures predicted by these models.

protoplanetary disks — circumstellar matter — planets and satellites:
formation

##  1 Introduction

Planets form from the gas and dust particles contained in disks orbiting young
stars. The overall time evolution and final dispersal of these planet-forming
disks play a major role in determining the properties of the planets born out
of these systems. Hence, understanding the physical processes which regulate
the late evolution and final dispersal of young disks is key for any theory of
planet formation (e.g., Alexander et al., 2014; Ercolano & Pascucci, 2017).

The X-ray photoevaporation model (Ercolano et al., 2008; Owen et al., 2010,
2011, 2012; Picogna et al., 2019) has been successful in explaining the
observed two-timescale and inside-out dispersal of most disks (Luhman et al.,
2010; Koepferl et al., 2013; Ercolano & Pascucci, 2017). In this model, soft
X-rays (0.1 keV $<E<$ 1 keV) from the central star can penetrate the disk
atmosphere and efficiently heat up the disk. In the regions of the disk in
which the sound speed of the heated gas exceeds the local escape velocity due
to stellar gravity, gas flows away from the disk, and a thermal wind is
centrifugally launched (e.g., Owen et al., 2010). When the mass loss rate due
to the wind and the mass accretion rate become comparable, a gap forms in the
disk, followed by an inner hole which quickly grows in stellocentric radius,
with the outer disk regions being dispersed from the inside-out.

Although these models can reproduce the spectral profile of emission lines
which are thought to originate in the disk wind (e.g., Ercolano & Owen, 2010,
2016; Weber et al., 2020), no conclusive evidence for a gap opened by
photoevaporation has been obtained yet. Interestingly, Ercolano et al. (2017)
have shown that models of X-ray photoevaporation can explain the gap observed
at a stellocentric radius of 1 au in the dust emission of the TW Hydrae disk
(Andrews et al., 2016) with the Atacama Large Millimeter/submillimeter Array
(ALMA), together with other properties of that system (e.g., mass accretion
rate).

As shown by Owen et al. (2013), hydrodynamic models of disk photoevaporation
coupled to radiative transfer post-processing techniques predict free-free
emission produced from a hot bound X-ray heated atmosphere and/or from the
photoevaporation wind with specific radial substructures (see also Pascucci et
al., 2014). In this work, we produce similar predictions for the free-free
radio emission from the X-ray photoevaporation models presented in Ercolano et
al. (2017) adjusted to the case of the TW Hydrae (TW Hya) disk.

A separate mechanism that has been proposed to contribute to the disk
dispersal, as well as to play a role in the redistribution of angular momentum
in the disk, are magneto-hydrodynamic (MHD) winds (Blandford & Payne, 1982;
Anderson et al., 2003; Bai et al., 2016). Contrary to the purely hydrodynamic
thermal winds driven by photoevaporation, the efficiency of MHD winds strongly
depends on the net magnetic flux retained in mature disks. Although the
evolution of MHD winds depend on a number of assumptions, some qualitative
differences with the photoevaporation winds have been identified in the
literature. For example, gas in MHD winds can be accelerated also at very
small stellocentric radii where photoevaporation winds are inefficient because
of the strong gravitational field of the star. Hence, it is interesting to
investigate whether these differences can be traced and tested also through
the radio emission of protoplanetary disks, which can be imaged at high
angular resolution using radio interferometry.

In this paper we follow the procedure outlined in Weber et al. (2020) to
calculate the radio emission expected from a MHD wind adapted to the case of
the TW Hya system, and compare these predictions to the expected radio
emission from models of X-ray driven photoevaporation.

We also show that a future radio interferometer like the Next Generation Very
Large Array (ngVLA, Murphy, 2018) would be capable of detecting and spatially
resolving several of the predicted substructures in the free-free emission at
cm wavelengths from the photoevaporation and MHD winds models. The comparison
with the theoretical predictions would provide key constraints to some
important parameters for models of disk dispersal, together with conclusive
evidence for the ability of either photoevaporation or MHD winds to produce
gaps and other substructures in the morphology of young disks.

Section 2 describes the main properties of the X-ray photoevaporation and MHD
wind models adapted to the case of TW Hya, and the methods used to derive the
model synthetic images at a wavelength of 1 cm. Section 3 shows the results of
these calculations, as well as predictions of future ngVLA observations to
detect and spatially resolve the expected disk substructures in the free-free
radio emission. Sections 4 and 5 outline the discussion and conclusions of
this work, respectively.

##  2 Methods: from disk simulations to ngVLA observations

In this Section we present the methods used to estimate the radio free-free
emission from ionized gas in the two contests of X-ray photoevaporation and
MHD disk models applied to the TW Hya disk. We then describe how we used the
results of these calculations to derive predictions for future radio
observations with the ngVLA. We start by summarizing the main properties of
the TW Hya system which are relevant for this work.

###  2.1 The case of TW Hydrae

TW Hya is the closest gas-rich young disk to Earth (distance $\approx$ 60.1
pc, Gaia Collaboration et al., 2018). With a stellar mass of $\approx
0.8-0.9~{}M_{\odot}$ (Huang et al., 2018), it is a member of the TW Hya
association, which has a relatively advanced age of $\sim 10$ Myr (e.g.,
Weinberger et al., 2013). Together with its nearly face-on orientation in the
sky (disk inclination $\approx$ 7 deg, Qi et al., 2004; Andrews et al., 2016),
these properties make this system a unique laboratory to investigate at the
best possible spatial resolution the imprints on the disk structure left by
the physical mechanisms which are responsible for the disk evolution and
dispersal.

Several hi-res optical observations (e.g., Debes et al., 2017) and
interferometric observations at mid-IR (e.g., Ratzka et al., 2007) and sub-
mm/mm wavelengths (e.g., Andrews et al., 2016) have resolved a variety of
different substructures at different stellocentric radii. In particular,
significant depletion of dust within the first couple of au from the star has
been inferred from the analysis presented in Calvet et al. (2002); Eisner et
al. (2006); Hughes et al. (2007); Menu et al. (2014), as well as from the
detection of a gap at 1 au from the star via direct imaging of the dust
continuum emission with ALMA at 0.87 mm (Andrews et al., 2016).

The low values inferred for the mass accretion rate of the TW Hya star ($\sim
4\times 10^{-10}-2\times 10^{-9}~{}M_{\odot}~{}\rm{yr}^{-1}$, Muzerolle et
al., 2000; Brickhouse et al., 2012; Manara et al., 2014) and its relatively
high X-ray luminosity ($\sim 2\times 10^{30}$ erg s-1) extrapolated in the 0.1
keV $<E<$ 10 keV spectral range (Robrade & Schmitt, 2006) make the inner
regions of this disk susceptible to significant mass loss by X-ray driven
photoevaporation (Owen et al., 2010). In fact, evidence of ongoing
photoevaporation in the TW Hya disk has been presented from the analysis of
the [NeII] 12.8 $\micron$ (Herczeg et al., 2007; Pascucci & Sterzik, 2009) and
[OI]6300 emission lines (Pascucci et al., 2011; Ercolano & Owen, 2016).

Ercolano et al. (2017) built on these results to show that also the dust
depletion observed at mid-IR to sub-mm wavelengths is consistent with the
expectations of X-ray photoevaporation disk models, making this the first
candidate object where a photoevaporative gap may have been imaged around the
time at which it is being created.

###  2.2 X-ray photoevaporation disk model

We used the methods described by Picogna et al. (2019) to obtain radiation-
hydrodynamical simulations of an X-ray photoevaporated disk around a young
star with properties similar to TW Hya as described by Ercolano et al. (2017).
We refer the reader to the cited articles for details about the numerical
methods and the input parameters. For the X-ray luminosity $L_{X}$ in the
models, we considered two values of $10^{29}$ and $10^{30}$ erg s-1,
respectively. Although the latter value is more in line with the observational
results for TW Hya (Section 2.1) we decided to run simulations also with a
lower value for an initial exploration of the dependence on $L_{X}$ of the
results from our models.

We then post-processed the grids using an improved version of the mocassin
radiative transfer code (Ercolano et al., 2003, 2005, 2008) to obtain images
at a wavelength of 1 cm. In Appendix A we describe the method used in this
work to reduce the Monte Carlo noise at long wavelengths.

###  2.3 MHD wind model

In order to estimate the spatial distribution of the free-free emission from a
MHD wind, we used the simple magnetocentrifugally driven wind model presented
in Weber et al. (2020). This is based on an analytical description with the
density and velocity structure obtained by Milliner et al. (2019) from the
Blandford & Payne (1982) axisymmetric self-similar solutions for a
magnetocentrifugally driven MHD wind for thin disks. More details can be found
in Milliner et al. (2019) and Weber et al. (2020). In particular, the model
considered here is model MHD-1 in Weber et al. (2020), with a mass-loss rate
of $10^{-8.6}~{}M_{\odot}/$yr and an X-ray luminosity for the star of
$L_{X}=2\times 10^{30}$ erg s-1.

Like for the X-ray photoevaporation disk model presented in the last section,
a synthetic image at a wavelength of 1 cm was obtained using the mocassin
code.

###  2.4 Simulations of the ngVLA observations

The model images obtained in Sections 2.2 and 2.3 were converted into
predictions for future observations with the ngVLA at $\lambda=$ 1 cm using
the CASA software package (McMullin et al., 2007).

Given the expected angular resolution ($\approx$ milliarcsec) and sensitivity
at cm wavelengths, a future ngVLA has the potential to spatially resolve
structures in the ionized gas emission in disks in nearby star forming
regions, as expected by models of disk photoevaporation and MHD disk winds.

To simulate the results of future interferometric observations with the ngVLA,
we adopted the same procedure as in Ricci et al. (2018) and Harter et al.
(2020), which used the SIMOBSERVE task to generate the visibility dataset in
the $(u,v)$ Fourier space, and the SIMNOISE task to add the noise by
corrupting the visibilities.

For the ngVLA simulations we considered the original ngVLA Rev B array
configuration with antennas distributed across the US Southwest and Mexico.
This configuration includes 214 antennas of 18 meter diameter, with baselines
up to 1000 km (Selina et al., 2018).

For the imaging of the interferometric visibilities we employed the CLEAN
algorithm with Briggs weighting, and adjusted the robust parameter to give a
reasonable synthesized beam and noise performance. In particular, the ngVLA
images were computed with a Briggs weighting scheme with robust parameter
$R=-2$ (uniform weighting). We also employed a multiscale clean approach to
better recover compact emission at both high brightness and larger and more
diffuse structures in the model. The disk was centered at the location of TW
Hydrae, i.e. RA(J2000)$=$11:01:51.90, Dec(J2000)$=$-34:42:17.0, and the
assumed distance is 60.1 pc (Gaia Collaboration et al., 2018).

Figure 1: Neutral hydrogen density map for the photoevaporation disk models
presented in Section 2.2 as a function of the radial distance from the star
and vertical coordinate from the disk midplane. Left and right panels are for
models with $L_{X}=10^{29}$ erg s-1 and $10^{30}$ erg s-1, respectively.
Vectors indicate the direction of the gas flow and their length represents the
speed of the flow with the reference speed of 10 km s-1 shown above the top
right corner of the right panel (relation between vector length and speed is
linear). The superimposed contour lines define the 90% emission region for the
free-free emission at 1 cm, after a 2D gaussian smoothing with $\sigma=2$ grid
cells to make the contour lines less noisy and confusing. Note that in both
panels the scaling of the x and y axes are different, and the significant
stretch of the y-axis was chosen to better highlight the wind structure in the
vertical direction.  Figure 2: Same as in Figure 1 but for the MHD wind model
presented in Section 2.3. The inner region within 1.5 au $\times$ 1.5 au in
the radial and vertical coordinates is shown in the overlaid panel on the top
left corner. The superimposed contour lines define the 90$\%$ emission region
for the free-free emission at 1 cm, after a 2D gaussian smoothing with
$\sigma=2$ grid cells to make the contour lines less noisy and confusing.

##  3 Results

Figures 1 and 2 present the density and velocity maps for the neutral hydrogen
as derived from the simulations presented in Sections 2.2 and 2.3 for the
photoevaporation and MHD wind models, respectively. The contour lines show the
regions in the $(r,z)$ plane where the bulk of the emission at 1 cm is
originated from.

The comparison between the X-ray photoevaporation and MHD wind models
highlights significant differences between the spatial distribution of the
emission. These differences are evident in Figure 3, which shows the radial
profiles for the surface brightness of the ionized gas emission at 1 cm for
three models calculated in this work.

The photoevaporation disk models predict emission with two main spatial
components. The first is localized at stellocentric distances of $\approx
0.5-0.9$ au, and corresponds to emission originated from the hot inner disk
atmosphere, which is gravitationally bound to the star-disk system. The second
component is spatially much broader, as it extends from about 1 to $4-5$ au
from the star. The radial shape of this component depends on the X-ray
luminosity of the star, as a higher luminosity ($L_{X}=10^{30}$ erg s-1)
produces significantly stronger emission, with a pronounced peak at a
stellocentric distance of about 3 au (Fig. 3, middle panel), than a model with
lower X-ray luminosity ($L_{X}=10^{29}$ erg s-1, Fig. 3, top panel). Contrary
to the first component which comes from the bound disk atmosphere, this
spatially broader component originates from ionized gas in the unbound wind.

Compared to the disk photoevaporation model, the predicted radial profile for
the MHD wind model is much smoother, and the emission extends to the very
inner edge of the grid (Fig. 3, bottom panel). This is a manifestation of the
fact that whereas a photoevaporative thermal wind cannot be launched very
close to the star, because of its intense gravity, a MHD wind can instead
develop also from those regions. We note here that the emission from the MHD
wind model is originated only from the wind as a disk is not present in the
analytical model presented in Section 2.3. The addition of a disk would not
modify the results shown here as the wind would strongly screen the disk from
ionising radiation.

The 2D maps are presented in Figure 4, which shows both the synthetic model
maps (left column) and the results of the simulated ngVLA observations (right
column). We note here that the rather pronounced multiple narrow rings which
are visible on the X-ray photoevaporation models, especially with
$L_{X}=10^{29}$ erg s-1, are an effect of the noise in the model (see Appendix
A), and therefore these small-scale substructures in the emission are not
physical. The total flux densities at 1 cm are 0.05, 0.15 and 1.8 mJy for the
photoevaporation model with $L_{X}=10^{29}$ and $10^{30}$ erg s-1, and the MHD
wind model, respectively.

The maps on the right column show how the ngVLA, with its sensitivity and
angular resolution of a few milliarcsec at 1 cm, would be able to resolve the
main structures expected in the continuum emission at 1 cm for the
photoevaporation and MHD wind models presented here for TW Hya.

##  4 Discussion

The results presented in the previous sections indicate that future radio
observations at high angular resolution with an instrument like the ngVLA can
spatially resolve sub-structures in the free-free gas emission expected by
models of disk photoevaporation and MHD winds.

The models of this study were adjusted to the specific case of TW Hyadrae, for
which the process of disk photoevaporation has been proposed to explain the
formation of the gap at 1 au from the star, together with other properties of
the TW Hya system (Ercolano & Pascucci, 2017). Figure 3 shows that disk
photoevaporation models predict a bright ring at 1 cm very close to the
location of the dust gap at sub-mm/mm wavelengths. Future radio observations
can test this prediction. Also, multi-wavelength observations at cm
wavelengths can be used to constrain the physical mechanism responsible for
the emission by measuring the spectral index (e.g., Rodmann et al., 2006).

Also the MHD wind model presented here and adjusted to the case of TW Hya
predicts significant emission at 1 cm from the ionized gas in the wind. The
hypotheses of photoevaporation-driven and MHD winds can be tested by
constraining the spatial profile of the emission, as radiation from an MHD
wind can extend to radii in the disk much closer to the star (Figure 3). This
is in general agreement with the interpretation of the different components
identified in the line profiles of typical wind diagnostics towards young
stars (e.g., Weber et al., 2020).

The results of our photoevaporation disk model with different values for the
X-ray luminosity of the star show different radial profiles for the surface
brightness at 1 cm, both in terms of the peak intensity and the radial shape
of the emission. For all the models calculated in this work, the bulk of the
emission is produced within about 5 au from the central star, corresponding to
about 0.08 arcsec at the distance of TW Hya. However, the intensity of the
first peak associated to emission from hot gas bound to the disk, as well as
the spatial extent of the emission from the photoevaporative wind further from
the star depends on the X-ray luminosity of the star itself. Given the
estimate for the X-ray luminosity for the TW Hya star, we expect the radio
emission to show two prominent peaks, at stellocentric distances of $\approx
0.5-1$ au and $3-3.5$ au, respectively. We note that the total flux from the
model with $L_{X}=10^{30}$ erg s-1 is more in line with the fluxes measured at
cm wavelengths, even though this system lacks strong constraints on the flux
from ionized gas at 1 cm (Pascucci et al., 2012). The expected flux from our
MHD wind model is a factor of $\sim 10\times$ higher than in the
photoevaporation case. It is worth noticing that our assumed mass loss rate is
comparable to some of the mass accretion rate estimates for the TW Hya young
star in the literature (e.g., Brickhouse et al., 2012). Although this is line
with some previous predictions from models of disk winds which account for
non-ideal MHD (e.g., Bai & Stone, 2013), other similar models with different
assumptions predict mass loss rates lower than mass accretion rates by factors
of $\sim 2-10$ (e.g., Gressel et al., 2020).

A broader investigation of the model parameters, not necessarily related to
the specific case of TW Hya, would be necessary to better characterize the
dependence of the radio-emission on the stellar X-ray luminosity, mass loss
rate, and other parameters of the disk-star system. Such investigation would
be important also for the interpretation of the radio emission already
observed towards several young low-mass stars in nearby star forming regions
(e.g., Rodmann et al., 2006; Macías et al., 2016; Ubach et al., 2017), as well
as to predict possible trends between some of the key physical quantities for
the star and disk, which would guide future surveys of young stellar objects
(YSOs) at radio wavelengths.

Although this work focuses on the case of TW Hya, which is a nearly face-on
disk, high-angular resolution observations of disks in nearly edge-on
orientations have the potential to resolve the vertical extent of the free-
free emission from the wind. This would provide another way to test the
predictions of X-ray photoevaporation and MHD wind models. Whereas computing
results at different inclination angles is beyond the scope of this study, we
plan on performing this investigation in a future work.

Figure 3:  Radial profiles of the surface brightness of the free-free emission
at a wavelength of 1 cm for the three models calculated in this work. From top
to bottom: X-ray photoevaporation model with $L_{X}=10^{29}$ erg s-1, X-ray
photoevaporation model with $L_{X}=10^{30}$ erg s-1, MHD wind models.

The disk models presented here neglect the contribution of the dust emission
at the wavelength of 1 cm considered in this study. This is justified by an
analysis of the dust continuum emission measured at 0.87 mm for the TW Hya
disk and presented in Andrews et al. (2016), which shows a peak brightness
temperature of 30 K. This can be considered as an upper limit for the
brightness temperature at 1 cm under the assumption that the dust emission is
optically thick even at cm wavelengths. The peak brightness would be only a
factor of 1.5 higher than the rms noise in the maps presented in Fig. 4 and in
general a factor $>10$ lower than the peak brightness temperature for the
predicted free-free emission for the disk models presented in this work.
Observations at wavelengths longer than 1 cm would further decrease the amount
of contamination from dust thermal emission, but would have the critical
disadvantage of providing poorer spatial resolution than at 1 cm. Since one of
our main goals is to quantify the potential of future hi-res observations to
test the model predictions for the spatial distribution of the disk free-free
emission, we leave to a future work the investigation of the free-free
emission at longer cm wavelengths.

Figure 4: Synthetic model images (left column) and ngVLA simulated
observations (right column) for the ionized gas continuum emission at 1 cm.
The top and middle panels show the maps from the photoevaporation disk models
with X-ray stellar luminosities of $10^{29}$ and $10^{30}$ erg s-1,
respectively. The bottom panels show the maps from the MHD disk models. The
unit for the surface brightness in each map of the models is Jy/pixel, where
the sides of the square pixels have sizes of 0.8 and 0.2 mas for the
photoevaporation and MHD models, respectively. In the ngVLA maps the
synthesized beam has sizes of $5.2~{}\rm{mas}\times 2.2~{}\rm{mas}$ with a
position angle of 5.3 degrees. The rms noise on each ngVLA map is 27 nJy/beam.
Note that the intensity color scales on the maps are different.

##  5 Conclusions

We have presented predictions for the free-free emission at 1 cm from models
of X-ray photoevaporation as well as of magnetocentrifugally driven MHD winds
with parameters adjusted to the case of the TW Hya YSO.

We showed that the X-ray photoevaporation models that have been proposed to
explain the opening of the gap observed by ALMA in the dust thermal emission
predict bright free-free emission in the region close to the gap in the TW Hya
disk. This emission is characterized by two main spatial components: one at
stellocentric distances of $\approx 0.5-1$ au, which is originated from hot
ionized gas in the inner disk atmosphere close to the star; one broader
component that extends to $\approx 4-5$ au from the star, which is associated
to unbound gas in the photoevaporative wind. The intensity and radial slope of
these components depend on the X-ray luminosity of the star.

The magnetocentrifugally driven MHD wind model predicts a spatially smoother
free-free emission from ionized gas, which extends to the closest regions to
the star in our model. Contrary to the case of disk photoevaporation where the
closest regions to the star are bound to the stellar gravity, the MHD wind in
our model is launched also from the innermost regions of the disk.

The predicted integrated fluxes from our models at 1 cm are in line with those
observed for TW Hya, as well as other YSOs in nearby star forming regions
using current radio facilities such as the Karl G. Jansky Very Large Array
(VLA) and the Australia Telescope Compact Array (ATCA). However, to spatially
resolve the predicted free-free emission and potentially distinguish between
the models presented here, observations providing angular resolutions below 10
mas at cm wavelengths are required. We showed that a future ngVLA, with its
current reference design, would have enough sensitivity and angular resolution
to detect and spatially resolve the main structures predicted by these models.

The results outlined in this work refer to the specific case of TW Hya. A
future investigation of the same models covering a broader region of the model
parameter space would be necessary to explore the dependence of these
predictions on the model parameters, and guide future observational surveys of
protoplanetary disks at radio wavelengths. These have the potential to shed
light on physical mechanisms which have been proposed to play a key role on
the evolution and dispersal of protoplanetary disks, and therefore on the
formation of planets.

We thank the anonymous referee for his/her comments which helped to clarify
the manuscript. This work was supported in part by the ngVLA Community Studies
program, coordinated by the National Radio Astronomy Observatory, which is a
facility of the National Science Foundation operated under cooperative
agreement by Associated Universities, Inc. BE acknowledges support from the
DFG Research Unit ”Transition Disks” (FOR 2634/1, ER 685/8-1). MW acknowledges
support from the DFG Research Unit ”Transition Disks” (FOR 2634/1, ER
685/11-1). This research was partially supported by the Excellence Cluster
ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy - EXC-2094 -
390783311. Part of the simulations have been carried out on the computing
facilities of the Computational Center for Particle and Astrophysics (C2PAP).

## References

  * Alexander et al. (2014) Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 475
  * Anderson et al. (2003) Anderson, J. M., Li, Z.-Y., Krasnopolsky, R., & Blandford, R. D. 2003, ApJ, 590, L107, doi: 10.1086/376824
  * Andrews et al. (2016) Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40, doi: 10.3847/2041-8205/820/2/L40
  * Bai & Stone (2013) Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76, doi: 10.1088/0004-637X/769/1/76
  * Bai et al. (2016) Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152, doi: 10.3847/0004-637X/818/2/152
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883, doi: 10.1093/mnras/199.4.883
  * Brickhouse et al. (2012) Brickhouse, N. S., Cranmer, S. R., Dupree, A. K., et al. 2012, ApJ, 760, L21, doi: 10.1088/2041-8205/760/2/L21
  * Calvet et al. (2002) Calvet, N., D’Alessio, P., Hartmann, L., et al. 2002, ApJ, 568, 1008, doi: 10.1086/339061
  * Debes et al. (2017) Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205, doi: 10.3847/1538-4357/835/2/205
  * Eisner et al. (2006) Eisner, J. A., Chiang, E. I., & Hillenbrand, L. A. 2006, ApJ, 637, L133, doi: 10.1086/500689
  * Ercolano et al. (2005) Ercolano, B., Barlow, M. J., & Storey, P. J. 2005, MNRAS, 362, 1038, doi: 10.1111/j.1365-2966.2005.09381.x
  * Ercolano et al. (2003) Ercolano, B., Barlow, M. J., Storey, P. J., & Liu, X. W. 2003, MNRAS, 340, 1136, doi: 10.1046/j.1365-8711.2003.06371.x
  * Ercolano et al. (2008) Ercolano, B., Drake, J. J., Raymond, J. C., & Clarke, C. C. 2008, ApJ, 688, 398, doi: 10.1086/590490
  * Ercolano & Owen (2010) Ercolano, B., & Owen, J. E. 2010, MNRAS, 406, 1553, doi: 10.1111/j.1365-2966.2010.16798.x
  * Ercolano & Owen (2016) —. 2016, MNRAS, 460, 3472, doi: 10.1093/mnras/stw1179
  * Ercolano & Pascucci (2017) Ercolano, B., & Pascucci, I. 2017, Royal Society Open Science, 4, 170114, doi: 10.1098/rsos.170114
  * Ercolano et al. (2017) Ercolano, B., Rosotti, G. P., Picogna, G., & Testi, L. 2017, MNRAS, 464, L95, doi: 10.1093/mnrasl/slw188
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gressel et al. (2020) Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126, doi: 10.3847/1538-4357/ab91b7
  * Harter et al. (2020) Harter, S. K., Ricci, L., Zhang, S., & Zhu, Z. 2020, ApJ, 905, 24, doi: 10.3847/1538-4357/abcafc
  * Herczeg et al. (2007) Herczeg, G. J., Najita, J. R., Hillenbrand, L. A., & Pascucci, I. 2007, ApJ, 670, 509, doi: 10.1086/521545
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42, doi: 10.3847/2041-8213/aaf740
  * Hughes et al. (2007) Hughes, A. M., Wilner, D. J., Calvet, N., et al. 2007, ApJ, 664, 536, doi: 10.1086/518885
  * Koepferl et al. (2013) Koepferl, C. M., Ercolano, B., Dale, J., et al. 2013, MNRAS, 428, 3327, doi: 10.1093/mnras/sts276
  * Luhman et al. (2010) Luhman, K. L., Allen, P. R., Espaillat, C., Hartmann, L., & Calvet, N. 2010, ApJS, 186, 111, doi: 10.1088/0067-0049/186/1/111
  * Macías et al. (2016) Macías, E., Anglada, G., Osorio, M., et al. 2016, ApJ, 829, 1, doi: 10.3847/0004-637X/829/1/1
  * Manara et al. (2014) Manara, C. F., Testi, L., Natta, A., et al. 2014, A&A, 568, A18, doi: 10.1051/0004-6361/201323318
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, Astronomical Society of the Pacific Conference Series, Vol. 376, CASA Architecture and Applications, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Menu et al. (2014) Menu, J., van Boekel, R., Henning, T., et al. 2014, A&A, 564, A93, doi: 10.1051/0004-6361/201322961
  * Milliner et al. (2019) Milliner, K., Matthews, J. H., Long, K. S., & Hartmann, L. 2019, MNRAS, 483, 1663, doi: 10.1093/mnras/sty3197
  * Murphy (2018) Murphy, E. 2018, Science with a Next Generation Very Large Array, Vol. 517
  * Muzerolle et al. (2000) Muzerolle, J., Calvet, N., Briceño, C., Hartmann, L., & Hillenbrand, L. 2000, ApJ, 535, L47, doi: 10.1086/312691
  * Owen et al. (2012) Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880, doi: 10.1111/j.1365-2966.2011.20337.x
  * Owen et al. (2011) Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13, doi: 10.1111/j.1365-2966.2010.17818.x
  * Owen et al. (2010) Owen, J. E., Ercolano, B., Clarke, C. J., & Alexand er, R. D. 2010, MNRAS, 401, 1415, doi: 10.1111/j.1365-2966.2009.15771.x
  * Owen et al. (2013) Owen, J. E., Scaife, A. M. M., & Ercolano, B. 2013, MNRAS, 434, 3378, doi: 10.1093/mnras/stt1254
  * Pascucci et al. (2012) Pascucci, I., Gorti, U., & Hollenbach, D. 2012, ApJ, 751, L42, doi: 10.1088/2041-8205/751/2/L42
  * Pascucci et al. (2014) Pascucci, I., Ricci, L., Gorti, U., et al. 2014, ApJ, 795, 1, doi: 10.1088/0004-637X/795/1/1
  * Pascucci & Sterzik (2009) Pascucci, I., & Sterzik, M. 2009, ApJ, 702, 724, doi: 10.1088/0004-637X/702/1/724
  * Pascucci et al. (2011) Pascucci, I., Sterzik, M., Alexander, R. D., et al. 2011, ApJ, 736, 13, doi: 10.1088/0004-637X/736/1/13
  * Picogna et al. (2019) Picogna, G., Ercolano, B., Owen, J. E., & Weber, M. L. 2019, MNRAS, 487, 691, doi: 10.1093/mnras/stz1166
  * Qi et al. (2004) Qi, C., Ho, P. T. P., Wilner, D. J., et al. 2004, ApJ, 616, L11, doi: 10.1086/421063
  * Ratzka et al. (2007) Ratzka, T., Leinert, C., Henning, T., et al. 2007, A&A, 471, 173, doi: 10.1051/0004-6361:20077357
  * Ricci et al. (2018) Ricci, L., Liu, S.-F., Isella, A., & Li, H. 2018, ApJ, 853, 110, doi: 10.3847/1538-4357/aaa546
  * Robrade & Schmitt (2006) Robrade, J., & Schmitt, J. H. M. M. 2006, A&A, 449, 737, doi: 10.1051/0004-6361:20054247
  * Rodmann et al. (2006) Rodmann, J., Henning, T., Chandler, C. J., Mundy, L. G., & Wilner, D. J. 2006, A&A, 446, 211, doi: 10.1051/0004-6361:20054038
  * Selina et al. (2018) Selina, R. J., Murphy, E. J., McKinnon, M., et al. 2018, in Astronomical Society of the Pacific Conference Series, Vol. 517, Science with a Next Generation Very Large Array, ed. E. Murphy, 15
  * Ubach et al. (2017) Ubach, C., Maddison, S. T., Wright, C. M., et al. 2017, MNRAS, 466, 4083, doi: 10.1093/mnras/stx012
  * Weber et al. (2020) Weber, M. L., Ercolano, B., Picogna, G., Hartmann, L., & Rodenkirch, P. J. 2020, MNRAS, 496, 223, doi: 10.1093/mnras/staa1549
  * Weinberger et al. (2013) Weinberger, A. J., Anglada-Escudé, G., & Boss, A. P. 2013, ApJ, 762, 118, doi: 10.1088/0004-637X/762/2/118

##  Appendix A Improved photon packet statistics in low emission spectral
regions

The standard version of the mocassin code expresses the radiation field in
term of photon packets, $\varepsilon(\nu)$, of constant energy, such that

$\varepsilon(\nu)=\varepsilon_{0}=n\,h\,\nu,$ (A1)

where $n$ is the number of photons contained in the packet and $\nu$ is the
frequency of the photons (Lucy 1999). When a photon of energy
$\varepsilon(\nu_{a})=\varepsilon_{0}$ is absorbed, it is immediately re-
emitted from the same location with frequency $\nu_{e}$, which is drawn
stochastically from the probability density function, based on the local
medium emissivities. The total energy of the packet remains the same, meaning
that in practice the number of photons contained in the re-emitted packet is
changed. Keeping the photon packet energy constant is a simple way of ensuring
conservation of energy throughout the simulation domain and thus accelerate
convergence (Lucy 1999).

Figure 5: Spectral energy distributions (SEDs) for the models shown in the
paper integrated over all lines of sight. The left panel shows that the
standard and flat PDF models agree on the global SEDs. The right panel shows
that the flat PDF method is able to achieve a higher signal to (Monte Carlo)
noise in the cm region, which is in the low energy tail of the SED.

The probability density function, $p(\nu)$, of an energy distribution $f(\nu)$
is given by

$p(\nu)=\frac{f(\nu)}{\int f(\nu)d\nu}.$ (A2)

The disadvantage of stochastically casting packets of constant energy from
$p(\nu)$ is that packets of frequencies where $f(\nu)$ is small will have a
very low probability of being cast. This means that the emerging spectral
energy distribution will have low signal to (Monte Carlo) noise in those low
$f(\nu)$ regions.

The cm-wavelength region in protoplanetary discs, which is dominated by free-
free emission of the gas, is one example where the standard mocassin approach
will fail (see Figure 5, blue and black lines). An alternative approach,
employed in this work, is to apply a ”flat” probability density function, i.e.
packets are assigned a random frequency when they are emitted, thus equally
sampling all spectral regions. This approach effecively means that packets are
being drawn from the following probability density function

$p(\nu)_{\rm{flat}}=\frac{1}{\int d\nu}.$ (A3)

In order to ensure energy conservation, the energy $\Delta E$ of the emitted
packet must then be weighted by the local probability density function
$p(\nu)$, such that

$\Delta\,E=\Delta\,E_{0}\frac{p(\nu)}{p(\nu)_{\rm{flat}}}=\Delta\,E_{0}\frac{f(\nu)}{\int
f(\nu)d\nu}\int d\nu.$ (A4)

In the case of the primary radiation field,
$\Delta\,E_{0}=\frac{L_{\rm{source}}}{N_{p}}$, where $L_{\rm{source}}$ is the
luminosity of the irradiating source and $N_{p}$ is the total number of energy
packets to be employed in the simulation. In the case of the secondary
radiation field $\Delta\,E_{0}$ is the energy of the absorbed energy packet.
The above ensures that the integral of all energy packets emerging from a
given location still reflects the local emission spectrum, thus enabling
conservation of energy at all locations.

Figure 5 compares the total emerging spectral energy distribution (SED,
integrated for all viewing angles) from our simulations using the standard
mocassin approach (black and blue lines) with the SED obtained using flat PDFs
(orange and green lines). The top panel shows the entire simulated frequency
range and demonstrates that the standard method and the flat PDF methods yield
the same results in terms of the global emerging SED. The bottom panel shows a
zoom of the same spectrum in the cm region, showing that the signal to (Monte
Carlo) noise is much lower in the SED obtained using the flat PDF method. The
SEDs were obtained using the same total number of energy packets, $N_{p}$.


We study behavioral metrics in an abstract coalgebraic setting. Given a coalgebra $\alpha\colon X \to FX$ in $\Set$, where the functor $F$ specifies the branching type, we define a framework for deriving pseudometrics on $X$ which measure the behavioral distance of states.

A first crucial step is the lifting of the functor $F$ on $\Set$ to a functor $\overline{F}$ in the category $\PMet$ of pseudometric spaces. We present two different approaches which can be viewed as generalizations of the Kantorovich and Wasserstein pseudometrics for probability measures. We show that the pseudometrics provided by the two approaches coincide on several natural examples, but in general they differ.

Then a final coalgebra for $F$ in $\Set$ can be endowed with a behavioral distance resulting as the smallest solution of a fixed-point equation, yielding the final $\overline{F}$-coalgebra in $\PMet$. The same technique, applied to an arbitrary coalgebra $\alpha\colon X \to FX$ in $\Set$,  provides the behavioral distance on $X$. Under some constraints we can prove that two states are at distance $0$ if and only if they are behaviorally equivalent.

§ INTRODUCTION

Increasingly, modelling formalisms are equipped with quantitative
information, such as probability, time or weight. Such quantitative
information should be taken into account when reasoning about
behavioral equivalence of system states, such as bisimilarity.  In
this setting the appropriate notion is not necessarily equivalence,
but a behavioral metric that measures the distance of the behavior of
two states. In a quantitative setting, it is often unreasonable to
assume that two states have exactly the same behavior, but it makes
sense to express that their behavior differs by some (small) value

The above considerations led to the study of behavioral metrics which
aims at quantifying the the distance between the behavior of
states. Since different states can have exactly the same behavior it
is quite natural to consider pseudometrics, which allow
different elements to be at zero distance.

Earlier contributions defined behavioral metrics in the setting of probabilistic systems <cit.> and of metric transition systems
Our aim is to generalize these ideas and to study behavioral metrics
in a general coalgebraic setting. The theory of coalgebra <cit.> is nowadays a well-established tool for defining and
reasoning about various state based transition systems such as
deterministic, nondeterministic, weighted or probabilistic automata.
Hence, it is the appropriate setting to ask and answer general
questions about behavioral metrics.

* How can we define behavioral metrics for transition
systems with different branching types? We provide a coalgebraic
framework in the category of pseudometric spaces $\PMet$ that allows
to define and reason about such metrics.

* Are the behavioral metrics canonical in some way? We
provide a natural way to define metrics by lifting functors from
$\Set$ to the category of pseudometric spaces.  In fact, we study
two liftings: the Kantorovich and the Wasserstein lifting and
observe that they coincide in many cases. This provides us with a
notion of canonicity and justification for the choice of metrics.

* Does the measurement of distances affect behavioral
equivalence? If we start by considering coalgebras in $\PMet$
(as, e.g., in <cit.>), it is not entirely clear
a priori whether the richer categorical structure influences the
notion of behavioral equivalence. In our setting we start with
coalgebras in $\Set$ and put distance measurements “on top”,
showing that, under some mild constraints, the original notion of
behavioral equivalence is not compromised, in the sense that two
states are behaviorally equivalent iff their distance is $0$.

* Are there generic algorithms to compute metrics?
Coalgebra is a valuable tool to define generic methods that can be
instantiated to concrete cases in order to obtain prototype
algorithms. In our case we give a (high-level) procedure for
computing behavioral distances on a given coalgebra, based
on determining the smallest solution of a fixed-point equation.

A central contribution of this paper is the lifting of a
functor $F$ from $\Set$ to $\PMet$. Given a pseudometric space
$(X,d)$, the goal is to define a suitable pseudometric on $FX$.  Such
liftings of metrics have been extensively studied in transportation
theory <cit.>, e.g. for the case of the
(discrete) probability distribution functor, which comes with a nice
analogy: assume several cities (with fixed distances between them) and
two probability distributions $s,t$ on cities, representing supply and
demand (in units of mass).  The distance between $s,t$ can be measured
in two ways: the first is to set up an optimal transportation plan
with minimal costs (in the following also called coupling) to
transport goods from cities with excess supply to cities with excess
demand. The cost of transport is determined by the product of mass and
distance. In this way we obtain the Wasserstein distance.  A different
view is to imagine a logistics firm that is commissioned to handle the
transport.  It sets prices for each city and buys and sells for this
price at every location.  However, it has to ensure that the price
function is nonexpansive, i.e., the difference of prices between two
cities is smaller than the distance of the cities, otherwise it will
not be worthwhile to outsource this task. This firm will attempt to
maximize its profit, which can be considered as the Kantorovich
distance of $s,t$. The Kantorovich-Rubinstein duality informs us that
these two views lead to the exactly same result, a very good argument
for the canonicity of this notion of distance.

It is our observation that these two notions of distance lifting can
analogously be defined for arbitrary functors $F$, leading to a rich
general theory.  The lifting has an evaluation function as
parameter. As concrete examples, besides the probability distribution
functor, we study the (finite) powerset functor (resulting in the
Hausdorff metric) and the coproduct and product bifunctors. In the
case of the product bifunctor we consider different evaluation
functions, each leading to a well-known product metric. The
Kantorovich-Rubinstein duality holds for these functors, but it does
not hold in general (we provide a counterexample).

After discussing functor liftings, we define coalgebraic behavioral
pseudometrics and answer the questions above. Specifically we
show how to compute distances on the final coalgebra as well as on
arbitrary coalgebras via fixed-point iteration and we prove that the
pseudometric obtained on the final coalgebra is indeed a metric.
In app:fibrations we discuss a fibrational perspective on our work and we compare with <cit.>.  All proofs for our results are in sec:proofs.

§ PRELIMINARIES, NOTATION & EVALUATION FUNCTIONS

We assume that the reader is familiar with the basic notions of category theory, especially with the definitions of functor, product, coproduct and weak pullbacks.

For a function $f\colon X \to Y$ and sets $A \subseteq X$, $B \subseteq Y$ we write $f[A] := \set{f(a) \mid a \in A}$ for the image of $A$ and $f^{-1}[B]= \set{a \in A \mid f(x) \in B}$ for the preimage of $B$. If $Y \subseteq [0,\infty]$ and $f,g\colon X \to Y$ are functions we write $f \leq g$ when $\forall x \in X:f(x)\leq g(x)$.

Given a natural number $n \in \N$ and a family $(X_i)_{i = 1}^n$ of sets $X_i$ we denote the projections of the (cartesian) product of the $X_i$ by $\pi_i^n\colon \prod_{i=1}^n X_i \to X_i,$ or just by $\pi_i$ if $n$ is clear from the context. For a source $(f_i\colon X \to X_i)_{i = 1}^n$ we denote the unique mediating arrow to the product by $\langle f_1,\dots,f_n\rangle \colon X \to \prod_{i=1}^{n}X_i$. Similarly, given a family of arrows $(f_i\colon X_i \to Y_i)_{i = 1}^n$, we write $f_1 \times \dots \times f_n = \langle f_1 \circ \pi_1,\dots,f_n \circ \pi_n \rangle\colon \prod_{i=1}^n X_i \to \prod_{i=1}^n Y_i$.

We quickly recap the basic ideas of coalgebras. Let $F$ be an endofunctor on the category $\Set$ of sets and functions. An $F$-coalgebra is just a function $\alpha \colon X \to FX$. Given another $F$-coalgebra $\beta\colon Y \to FY$ a coalgebra
homomorphism from $\alpha$ to $\beta$ is a function $f \colon A \to B$
such that $\beta \circ f = Ff \circ \alpha$. We call an $F$-coalgebra
$\kappa\colon \Omega \to F\Omega$ final if for any other
coalgebra $\alpha \colon X \to FX$ there is a unique coalgebra
homomorphism $\finalO \colon X \to \Omega$. The final coalgebra need
not exist but if it does it is unique up to isomorphism. It can be
considered as the universe of all possible behaviors. If we have an
endofunctor $F$ such that a final coalgebra $\kappa\colon \Omega \to
F\Omega$ exists then for any coalgebra $\alpha\colon X \to FX$ two
states $x_1, x_2 \in X$ are said to be behaviorally equivalent
if and only if $\final{x_1} = \final{x_2}$.

We now introduce some preliminaries about (pseudo)metric spaces. Our
(pseudo)metrics assume values in a closed interval $[0,\top]$, where
$\top \in (0,\infty]$ is a fixed maximal element (for our examples we
will use $\top = 1$ or $\top=\infty$). In this way the set of
(pseudo)metrics over a fixed set with pointwise order is a complete lattice (since $\reals$ is) and
the resulting category of pseudometric spaces is complete and

Given a set $X$, a pseudometric on $X$ is a function $d\colon X\times X\to \reals$ such that for all $x,y,z\in X$, the following axioms hold: $d(x,x) = 0$ (reflexivity), $d(x,y) = d(y,x)$ (symmetry), $d(x,z) \le d(x,y)+d(y,z)$ (triangle inequality). If additionally $d(x,y)=0$ implies $x=y$, $d$ is called a metric. A (pseudo)metric space is a pair $(X,d)$ where $X$ is a set and $d$ is a (pseudo)metric on $X$.

By $d_e\colon \reals^2\to \reals$ we denote the ordinary Euclidean distance on $\reals$, i.e., $d_e(x,y) = |x-y|$ for $x,y\in\reals\setminus\{\infty\}$, and – where appropriate – $d_e(x,\infty) = \infty$ if $x\neq \infty$ and $d_e(\infty,\infty) = 0$. Addition is defined in the usual way, in particular $x + \infty = \infty$ for $x\in\prealinf$.

Hereafter, we only consider those functions between two pseudometric spaces that do not increase the distances.

Let $(X,d_X)$, $(Y,d_Y)$ be pseudometric spaces. A function $f\colon X\to Y$ is called nonexpansive if $d_Y \circ (f \times f) \le d_X$. In this case we write[The nonexpansive functions correspond exactly to the Lipschitz-continuous functions with Lipschitz constant $L \leq 1$, which is the reason for writing $f \colon (X, d_X) \nonexpansiveTo (Y, d_Y)$.] $f \colon (X, d_X) \nonexpansiveTo (Y, d_Y)$. If equality holds, $f$ is called an isometry.

For our purposes it will turn out to be useful to consider the following alternative characterization of the triangle inequality using the concept of nonexpansive functions.

A symmetric function $d\colon X^2\to \reals$ with $d(x,x) = 0$ for all $x \in X$ satisfies the triangle inequality iff for all $x \in X$ the function $d(x,\_)\colon X \to \reals$ is nonexpansive.

As stated before, our definition of a pseudometric gives rise to a suitably rich category.

For a fixed $\top \in (0,\infty]$ we denote by $\PMet$ the category of all pseudometric spaces and nonexpansive functions.

This category is complete and cocomplete (see prop:PMetcomplete) and, in particular, it has products and coproducts as we will see in Examples <ref> and <ref>. We now introduce two motivating examples borrowed from <cit.> and <cit.>.

[Probabilistic Transition Systems and Behavioral Distance]

We regard probabilistic transition systems as coalgebras of the form $\alpha\colon X\to \Distributions (X+\mathbf{1})$, where $\Distributions $ is the probability distribution functor (with finite support) which maps a set $X$ to the set $\Distributions X = \{P\colon X\to [0,1]\mid \sum_{x\in X} P(x) = 1, P \text{ has finite support}\}$ and a function $f \colon X \to Y$ to the function $\Distributions f \colon \Distributions X \to \Distributions Y, P \mapsto \lambda y.\sum_{x \in   f^{-1}[\set{y}]}P(x)$. Here $\alpha(x)(y)$, for $x,y\in X$, denotes the probability of a transition from a state $x$ to $y$ and $\alpha(x)(\checkmark)$ stands for the probability of terminating from $x$ (we use $\checkmark$ for the single element of the set $\mathbf{1}$).

In <cit.> a metric for the continuous version of these systems is introduced, by considering a discount factor $c \in (0,1)$. In the discrete case we obtain the behavioral distance $d\colon X^2\to [0,1]$, defined as the least solution of the equation $d(x,y) = \overline{d}(\alpha(x),\alpha(y))$, where $x,y\in X$ and $\overline{d}\colon (\Distributions (X+\mathbf{1}))^2\to [0,1]$ is defined in two steps: First, $\hat{d}\colon (X+\mathbf{1})^2\to [0,1]$ is defined as $\hat{d}(x,y) = c\cdot d(x,y)$ if $x,y\in X$, $\hat{d}(\checkmark,\checkmark) = 0$ and $1$ otherwise. Then, for all $P_1,P_2\in \Distributions (X+\mathbf{1})$, $\overline{d}(P_1,P_2)$ is defined as the supremum of all values $\sum_{x \in X+\mathbf{1}} f(x)\cdot \big(P_1(x) - P_2(x)\big)$, with $f \colon (X+\mathbf{1},\hat{d}) \nonexpansiveTo ([0,1],d_e)$ being an arbitrary nonexpansive function.  As we will further discuss in exa:probability-distribution-functor, $\overline{d}$ is the Kantorovich pseudometric given by the space $(X+\mathbf{1},\hat{d})$.

We consider a concrete example from <cit.>, illustrated on the left of fig:example-ts. The behavioral distance of $u$ and $z$ is $d(u,z) = 1$ and hence $d(x,y) = c\cdot \epsilon$.

[Metric Transition Systems and Propositional Distances]

We give another example based on the notions of
<cit.>. A finite set $\Sigma =
\{r_1,\dots,r_n\}$ of propositions is given and each proposition
$r\in\Sigma$ is associated with a pseudometric space $(M_r,d_r)$. A
valuation $u$ is a function with domain $\Sigma$ that assigns to
each $r\in\Sigma$ an element of $M_r$. We denote the set of all
valuations by $\mathcal{U}[\Sigma]$. A metric transition system is a
tuple $M = (S,\tau,\Sigma,[\cdot])$ with a set $S$ of states, a
transition relation $\tau\subseteq S\times S$, a finite set $\Sigma$
of propositions and a valuation $[s]$ for each state $s\in S$. We
write $\tau(s)$ for $\{s'\in S\mid (s,s')\in\tau\}$ and require that
$\tau(s)$ is finite.

In <cit.> the propositional distance
between two valuations is given by $\LiftedFunctor{\mathit{pd}}(u,v)
= \linebreak \max_{r\in\Sigma} d_r(u(r),v(r))$ for
$u,v\in\mathcal{U}[\Sigma]$.  The (undirected) branching distance
$d\colon S\times S\to\preal$ is defined as the smallest fixed-point
of the following equation, where $s,t\in S$:

\begin{equation}\label{eq:Hausd}
d(s,t) = \max\{\LiftedFunctor{\mathit{pd}}([s],[t]),
\max_{s'\in\tau(s)}\min_{t'\in\tau(t)}
d(s',t'),\max_{t'\in\tau(t)}\min_{s'\in\tau(s)} d(s',t') \}
\end{equation}

Note that, apart from the first argument, this coincides with the Hausdorff distance.

We consider an example which appears similarly
in <cit.> (see
Figure <ref>, right) with a single proposition $r\in
\Sigma$, where $M_r = [0,1]$ is equipped with the Euclidean distance
$d_e$. According to (<ref>), $d(x_1,y_1)$ equals the
Hausdorff distance of the reals associated with the sets of
successors, which is $0.3$ (since this is the maximal distance of
any successor to the closest successor in the other set of
successors, here: the distance from $y_3$ to $x_3$).

In order to model such transition systems as coalgebras we define
the following $n$-ary auxiliary functor: $G(X_1,\dots,X_n) =
\{u\colon \Sigma\to X_1+\dots+X_n\mid u(r_i) \in X_i \}$.  Then
coalgebras are of the form $\alpha\colon S\to
G(M_{r_1},\dots,M_{r_n})\times \PowersetFinite(S)$, where
$\PowersetFinite$ is the finite powerset functor and
$\alpha(s) = ([s],\tau(s))$. As we will see later in
ex:metric-ts-2, the right-hand side of
(<ref>) can be seen as lifting a metric $d$ on $X$ to a
metric on $G(M_{r_1},\dots,M_{r_n})\times
\PowersetFinite(X)$.

[->,auto,node distance=1cm,semithick,>=latex]
every state=[fill=black!10,text=black,inner sep=0.5pt, minimum size=17pt]
every node=[text=blue]

[state] (x) $x$;
(center) [below=0.3cm of x] ;
[state] (u) [left=1cm of center]$u$;
[state,accepting] (z) [right=1cm of center] $z$;
[state] (y) [below=0.3 cm of center] $y$;

(x) edge node[above left] $\frac{1}{2}-\epsilon$ (u);
(x) edge node[above right] $\frac{1}{2}+\epsilon$ (z);
(y) edge node[below left] $\frac{1}{2}$ (u);
(y) edge node[below right] $\frac{1}{2}$ (z);
(u) edge[loop left] node[left] $1$ (u);

[->,auto,node distance=1cm,semithick,>=latex]
every state=[fill=black!10,text=black,inner sep=0.5pt, minimum size=17pt]
every node=[text=blue]

[state] (x1) $x_1$;
(p1) [right=0mm of x1] $0$;
[state] (x2) [below left=of x1] $x_2$;
(p2) [right=0mm of x2] $0.4$;
[state] (x3) [below right=of x1] $x_3$;
(p3) [left=0mm of x3] $0.7$;

[state] (y1) [right=3.5cm of x1] $y_1$;
(q1) [right=0mm of y1] $0$;
[state] (y2) [below left=of y1] $y_2$;
(q2) [right=0mm of y2] $0.5$;
[state] (y3) [below right=of y1] $y_3$;
(q3) [left=0mm of y3] $1$;

(x1) edge (x2);
(x1) edge (x3);
(x2) edge[loop left] (x2);
(x3) edge[loop right] (x3);

(y1) edge (y2);
(y1) edge (y3);
(y2) edge[loop left] (y2);
(y3) edge[loop right] (y3);

A probabilistic transition system (left) and a metric transition system (right).

Generalizing from the examples, we now establish a general framework for deriving such behavioral distances. In both cases, the crucial step is to find, for a functor $F$, a way to lift a pseudometric on $X$ to a pseudometric on $FX$. Based on this, one can set up a fixed-point equation and define behavioral distance as its smallest solution. Hence, in the next sections we describe how to lift an endofunctor $F$ on $\Set$ to an endofunctor on $\PMet$.

Let $U\colon \PMet \to \Set$ be the forgetful functor which maps every pseudometric space to its underlying set. A functor $\LiftedFunctor{F}\colon\PMet\to\PMet$ is called a lifting of a functor $F\colon \Set \to \Set$ if it satisfies $U\circ \LiftedFunctor{F} = F\circ U$.

It is not difficult to prove that such a lifting is always
monotone on pseudometrics over a common set, i.e. for any two
pseudometrics $d_1 \leq d_2$ on the same set $X$, we also have
$\LiftedMetric{F}{d_1} \leq \LiftedMetric{F}{d_2}$ where
$\LiftedMetric{F}{d_i}$ are the pseudometrics on $FX$ obtained by
applying $\LiftedFunctor{F}$ to $(X, d_i)$ (see prop:monotone).
Similarly to predicate lifting of coalgebraic modal logic <cit.>, liftings on $\PMet$ can be conveniently defined via an evaluation function.

Let $F$ be an endofunctor on $\Set$. An evaluation function for $F$ is a function $\ev_F\colon F\reals \to \reals$. Given such a function, we define the evaluation functor to be the endofunctor $\EvaluationFunctor{F}$ on $\Set/\reals$, the slice category[The slice category $\Set/\reals$ has as objects all functions $g\colon X\to\reals$ where $X$ is an arbitrary set. Given $g$ as before and $h\colon Y \to \reals$, an arrow from $g$ to $h$ is a function $f\colon X \to Y$ satisfying $h \circ f = g$.] over $\reals$, via $\EvaluationFunctor{F}(g) = \ev_F\circ Fg$ for all $g \in \Set/\reals$. On arrows $\EvaluationFunctor{F}$ is defined as $F$.

§ LIFTING FUNCTORS TO PSEUDOMETRIC SPACES À LA KANTOROVICH

Let us now consider an endofunctor $F$ on $\Set$ with an evaluation function $\ev_F$. Given a pseudometric space $(X,d)$, our first approach will be to take the smallest possible pseudometric $d^F$ on $FX$ such that, for all nonexpansive functions $f\colon (X,d) \nonexpansiveTo (\reals,d_e)$, also $\EvaluationFunctor{F}f\colon (FX,d^F)\nonexpansiveTo (\reals,d_e)$ is nonexpansive again, i.e. we want to ensure that for all $t_1,t_2 \in FX$ we have $d_e(\EvaluationFunctor{F}f(t_1), \EvaluationFunctor{F}f(t_2))\leq d^F(t_1,t_2)$. This idea immediately leads us to the next definition.

[Kantorovich Pseudometric & Kantorovich Lifting]definitiondef:kantorovich
Let $F\colon \Set \to \Set$ be a functor with an evaluation function $\ev_F$. For every pseudometric space $(X,d)$ the Kantorovich pseudometric on $FX$ is the function $\Kantorovich{F}{d}\colon FX\times FX\to \reals$, where for all $t_1,t_2 \in FX$:

\begin{align*}
\Kantorovich{F}{d}(t_1,t_2) := \sup \{d_e(\EvaluationFunctor{F}f(t_1),\EvaluationFunctor{F}f(t_2)) \mid f\colon (X,d) \nonexpansiveTo (\reals,d_e)\}\,.
\end{align*}

The Kantorovich lifting of the functor $F$ is the functor $\LiftedFunctor{F}\colon \PMet\to\PMet$ defined as $\LiftedFunctor{F}(X,d) = (FX,\Kantorovich{F}{d})$ and $\LiftedFunctor{F}f = Ff$.

It is easy to show that $\Kantorovich{F}{d}$ is indeed a
pseudometric. Since $\LiftedFunctor{F}$ inherits the preservation of
identities and composition of morphisms from $F$ we can prove that
nonexpansive functions are mapped to nonexpansive functions and
isometries to isometries.

The Kantorovich lifting $\LiftedFunctor{F}$ of a functor $F$ preserves isometries.

We chose the name Kantorovich because our definition is reminiscent of the Kantorovich pseudometric in probability theory. If we take the proper combination of functor and evaluation function, we can recover that pseudometric (in the discrete case) as the first instance for our framework.

[Probability Distribution Functor]exampleexa:probability-distribution-functor
We take $\top = 1$ and the probability distribution functor $\Distributions $ from ex:probabilistic-1. As evaluation function we take the expected value of the identity on $[0,1]$, i.e. for any $P \in \Distributions [0,1]$ we have $\ev_\Distributions (P) = \ExpectedValue{P}{\id_{[0,1]}} = \sum_{x\in [0,1]} x\cdot P(x)$ yielding $\EvaluationFunctor{\Distributions }g (P) = \ExpectedValue{P}{g} = \sum_{x\in [0,1]} g(x)\cdot P(x)$ for all $g\colon X \to [0,1]$. For every pseudometric space $(X,d)$ we obtain the discrete Kantorovich pseudometric $\Kantorovich{\Distributions }{d}\colon \Distributions X \times \Distributions X \to [0,1]$, defined as, for all $P_1, P_2 \in \Distributions X$, $\Kantorovich{\Distributions }{d}(P_1,P_2) = \sup\{ \sum_{x \in X} f(x)\cdot \big(P_1(x) - P_2(x)\big) \mid f \colon (X,d) \nonexpansiveTo ([0,1],d_e)\}$.

In general Kantorovich liftings do not preserve metrics, as shown by the following example.

Let $F\colon \Set \to \Set$ be given as $FX = X \times X$ on sets and $Ff = f \times f$ on functions and take $\top = \infty$, $\ev_F\colon F[0,\infty] \to [0,\infty], \ev_F(r_1, r_2) = r_1+r_2$. For a metric space $(X,d)$ with $|X|\geq 2$ let $t_1 = (x_1, x_2) \in FX$ with $x_1 \not = x_2$ and define $t_2 := (x_2, x_1)$. Clearly $t_1 \not = t_2$ but for every nonexpansive function $f\colon (X,d) \nonexpansiveTo (\reals, d_e)$ we have $\EvaluationFunctor{F}f(t_1) = f(x_1) + f(x_2) = f(x_2) + f(x_1) = \EvaluationFunctor{F}f(t_2)$ and thus $\Kantorovich{F}{d}(t_1,t_2) = 0$.

§ WASSERSTEIN PSEUDOMETRIC AND KANTOROVICH-RUBINSTEIN DUALITY

We have seen that our first lifting approach bears close resemblance to the original Kantorovich pseudometric on probability measures. In that context there exists another pseudometric, the Wasserstein pseudometric, which under certain conditions coincides with the Kantorovich pseudometric. We will define a generalized version of the Wasserstein pseudometric and compare it with our generalized Kantorovich pseudometric. To do that we first need to define how we can couple elements of $FX$.

Let $F \colon \Set \to \Set$ be a functor and $n \in \N$. Given a set $X$ and $t_i \in FX$ for $1 \leq i \leq n$ we call an element $t \in F(X^n)$ such that $F\pi_i(t) = t_i$ a coupling of the $t_i$ (with respect to $F$). We write $\Couplings{F}(t_1, t_2, \dots, t_n)$ for the set of all these couplings.

If $F$ preserves weak pullbacks, we can define new couplings based on given ones.

[Gluing Lemma]lemmalem:coupling
Let $F \colon \Set \to \Set$ be a weak pullback preserving functor, $X$ a set, $t_1, t_2, t_3 \in FX$, $t_{12} \in \Couplings{F}(t_1,t_2)$, and $t_{23} \in \Couplings{F}(t_2,t_3)$ be couplings. Then there is a coupling $t_{123} \in \Couplings{F}(t_1, t_2, t_3)$ such that $F(\langle \pi_1^3, \pi_2^3\rangle)(t_{123}) = t_{12}$ and $F(\langle\pi_2^3, \pi_3^3\rangle)(t_{123}) = t_{23}$.

This lemma already hints at the fact that our new lifting will only work for weak pullback preserving functors, which is a standard requirement in coalgebra. In addition to that we have to impose three extra conditions on the evaluation functions.

[Well-Behaved Evaluation Function]definitiondef:well-behaved
Let $\ev_F$ be an evaluation function for a functor $F\colon \Set \to \Set$. We call $\ev_F$ well-behaved if it satisfies the following conditions:

* $\EvaluationFunctor{F}$ is monotone, i.e., for $f,g\colon X\to\reals$ with
$f\le g$, we have $\EvaluationFunctor{F}f\le\EvaluationFunctor{F}g$.

* For each $t\in F(\reals^2)$ it holds that $d_e(\ev_F(t_1), \ev_F(t_2)) \leq \EvaluationFunctor{F}d_e(t)$ for $t_i : = F\pi_i(t)$.

* $\ev_F^{-1}[\set{0}] = Fi[F\{0\}]$ where $i \colon \set{0} \hookrightarrow\reals$ is the inclusion map.

While the first condition of this definition is quite natural, the other two need to be explained. Condition 2 is needed to ensure that $\EvaluationFunctor{F}\id_\reals = \ev_F\colon F\reals\to\reals$ is nonexpansive once $d_e$ is lifted to $F\reals$ (cf. the intuition behind the Kantorovich lifting, where we ensure that $\EvaluationFunctor{F}f$ is nonexpansive whenever $f$ is nonexpansive). Furthermore
Condition 3 intuitively says that exactly the elements of $F\{0\}$ are mapped to $0$ via $\ev_F$. Before we define the Wasserstein pseudometric and the corresponding lifting, we take a look at an example of a functor together with a well-behaved evaluation function.

[Finite Powerset Functor]exampleex:evaluation-max
Let $\top = \infty$. We take the finite powerset functor $\PowersetFinite$ with evaluation function $\max\colon \PowersetFinite(\prealinf) \to \prealinf$ with $\max \emptyset = 0$. This evaluation function is well-behaved whereas $\min\colon \PowersetFinite(\prealinf) \to \prealinf$ is not well-behaved.

Let $F\colon \Set \to \Set$ be a weak-pullback preserving functor with well-behaved evaluation function $\ev_F$. For every pseudometric space $(X,d)$ the Wasserstein pseudometric on $FX$ is the function $\Wasserstein{F}{d} \colon FX \times FX \to \reals$ given by, for all $t_1,t_2 \in FX$,

\begin{align*}
\Wasserstein{F}{d}(t_1, t_2) := \inf \{\EvaluationFunctor{F}d(t) \mid t \in \Couplings{F}(t_1,t_2)\}\,.
\end{align*}

We define the Wasserstein lifting of $F$ to be the functor $\LiftedFunctor{F}\colon \PMet\to\PMet$, $\LiftedFunctor{F}(X,d) = (FX,\Wasserstein{F}{d})$, $\LiftedFunctor{F}f = Ff$.

This time it is not straightforward to prove that $\Wasserstein{F}{d}$ is a pseudometric, so we explicitly provide the following result. Its proof relies on all properties of well-behavedness of $\ev_F$ and uses lem:coupling which explains why we need a weak pullback preserving functor.

The Wasserstein pseudometric is a well-defined pseudometric on $FX$.

In contrast to that, it is not hard to show functoriality of $\LiftedFunctor{F}$ and, as in the Kantorovich case, the lifted functor preserves isometries.

The Wasserstein lifting $\LiftedFunctor{F}$ of a functor $F$ preserves isometries.

In contrast to our previous approach, metrics are preserved in certain situations.

[Preservation of Metrics]propositionprop:preservation
Let $(X,d)$ be a metric space and $F$ be a functor. If the infimum in def:wasserstein is a minimum for all $t_1,t_2 \in FX$ where $\Wasserstein{F}{d}(t_1,t_2) = 0$ then $\Wasserstein{F}{d}$ is a metric, thus also $\LiftedFunctor{F}(X,d) = (FX,\Wasserstein{F}{d})$ is a metric space.

Please note that a similar restriction for the Kantorovich lifting (i.e. requiring that the supremum in def:kantorovich is a maximum) does not yield preservation of metrics: In exa:kantorovich-lifting-does-not-preserve-metrics the supremum is always a maximum but we do not get a metric.

Let us now compare both lifting approaches. Whenever it is defined, the Wasserstein pseudometric is an upper bound for the Kantorovich pseudometric.

Let  $F$  be a  weak pullback preserving functor with well-behaved evaluation function. Then for all pseudometric spaces $(X,d)$ it holds that $\Kantorovich{F}{d} \leq \Wasserstein{F}{d}$.

In general this inequality may be strict in general, as the following example shows.

The functor of
preserves weak pullbacks and the evaluation function is
well-behaved. We continue the example and take $t_1=(x_1,x_2)$,
$t_2=(x_2,x_1)$. The unique coupling $t \in \Gamma_F(t_1,t_2)$ is
$t=(x_1, x_2, x_2, x_1)$. Using that $d$ is a metric we conclude that
$\Wasserstein{F}{d}(t_1, t_2) = \EvaluationFunctor{F}d(t) = d(x_1,x_2)
+ d(x_2,x_1) = 2d(x_1,x_2) > 0 = \Kantorovich{F}{d}(t_1,t_2)$.

When the inequality can be replaced by an equality we will
in the following say that the Kantorovich-Rubinstein duality holds.
In this case we obtain a canonical notion of distance on $FX$, given a
pseudometric space $(X,d)$. To calculate the distance of $t_1,t_2\in
FX$ it is then enough to find a nonexpansive function $f\colon
(X,d)\nonexpansiveTo (\reals,d_e)$ and a coupling $t \in
\Gamma_F(t_1,t_2)$ such that
$d_e(\EvaluationFunctor{F}f(t_1),\EvaluationFunctor{F}f(t_2)) =
\EvaluationFunctor{F}d_e(t)$. Then, due to
prop:wasserstein-vs-kantorovich, this value equals
$\Kantorovich{F}{d}(t_1,t_2) = \Wasserstein{F}{d}(t_1,t_2)$.  We will
now take a look at some examples where the duality holds.

[Identity Functor]

Take $F = \Id$ with the identity evaluation map $\ev_\Id = \id_\reals$. For any $t_1,t_2\in X$, $t := (t_1, t_2)$ is the unique coupling of $t_1, t_2$. Hence, $\Wasserstein{F}{d}(t_1, t_2) = d(t_1,t_2)$. With the function $d(t_1,\_) \colon (X,d) \nonexpansiveTo (\reals,d_e)$ we obtain duality because we have $d(t_1,t_2) = d_e(d(t_1, t_1), d(t_1,t_2)) \leq \Kantorovich{F}{d}(t_1,t_2) \leq \Wasserstein{F}{d}(t_1,t_2) = d(t_1,t_2)$ and thus equality. Similarly, if we define $\ev_\Id(r) = c\cdot r$ for $r\in \reals$, $0<c\leq1$, the Kantorovich and Wasserstein liftings coincide and we obtain the discounted distance $\Kantorovich{F}{d}(t_1,t_2) = \Wasserstein{F}{d}(t_1,t_2) = c\cdot d(t_1,t_2)$.

[Probability Distribution Functor]exampleexa:probability-distribution-functor2
The functor $\Distributions $ of exa:probability-distribution-functor preserves weak pullbacks <cit.> and the evaluation function $\ev_\Distributions $ is well-behaved. We recover the usual Wasserstein pseudometric
$\Wasserstein{\Distributions }{d}(P_1,P_2)= \inf \{\sum_{x_1,x_2 \in X} d(x_1,x_2) \cdot P(x_1, x_2) \mid P \in \Couplings{\Distributions }(P_1, P_2)\}$
and the Kantorovich-Rubinstein duality <cit.> from transportation theory for the discrete case.

[Finite Powerset Functor & Hausdorff Pseudometric]exampleexa:hausdorff
Let $\top = \infty$, $F = \PowersetFinite$ with evaluation map $\ev_{\PowersetFinite}\colon \PowersetFinite(\prealinf) \to \prealinf$, $\ev_{\PowersetFinite}(R) = \max R$ with $\max \emptyset = 0$ (as in ex:evaluation-max). In this setting we obtain duality and both pseudometrics are equal to the Hausdorff pseudometric $d_H$ on $\PowersetFinite(X)$ which is defined as, for all $X_1,X_2\in \PowersetFinite(X)$,

\begin{align*}
d_H(X_1,X_2) = \max\left\{\max_{x_1\in X_1} \min_{x_2\in X_2} d(x_1,x_2), \max_{x_2\in X_2} \min_{x_1\in X_1} d(x_1,x_2) \right\}\,.
\end{align*}

Note that the distance is $\infty$, if either $X_1$ or $X_2$ is empty.

It is also illustrative to consider the countable powerset
functor.  Using the supremum as evaluation function, one obtains again
the Hausdorff pseudometric (with supremum/infimum replacing
maximum/minimum). However,  in this case the Hausdorff
distance of different countable sets might be $0$, even if we lift
a metric. This shows that in general the Wasserstein lifting does not preserve metrics but we need an extra condition, e.g. the one in prop:preservation.

§ LIFTING MULTIFUNCTORS

Our two approaches can easily be generalized[The details are spelled out in app:multifunctors, here we provide just the basic ideas.] to lift a multifunctor $F\colon\Set^n\to \Set$ (for $n \in \N$) in a similar sense as given by def:lifting to a multifunctor $\LiftedFunctor{F}\colon \PMet^n \to \Set$. The only difference is that we start with $n$ pseudometric spaces instead of one. Now we need an evaluation function $\ev_F\colon F(\reals, \dots, \reals) \to \reals$ which we call well-behaved if it satisfies conditions similar to def:well-behaved and which gives rise to an evaluation multifunctor $\EvaluationFunctor{F} \colon(\Set/\reals)^n \to \Set/\reals$. Given $t_1,t_2 \in F(X_1,\dots,X_n)$ we write again $\Couplings{F}(t_1,t_2) \subseteq F(X_1^2,\dots, X_n^2)$ for the set of couplings which is defined analogously to def:coupling. For pseudometrics $d_i\colon X_i^2\to \reals$, we can then define the Kantorovich/Wasserstein pseudometric $\KantorovichMulti{F}{d}{1,\dots,n}, \WassersteinMulti{F}{d}{1,\dots,n}\colon F(X_1,\dots,X_n)\times F(X_1,\dots,X_n)\to \reals$, as $\KantorovichMulti{F}{d}{1,\dots,n}(t_1,t_2) := \sup \{d_e(\EvaluationFunctor{F}(f_1,\dots,f_n)(t_1),\EvaluationFunctor{F}(f_1,\dots,f_n)(t_2)) \mid {f_i\colon (X_i,d_i) \nonexpansiveTo (\reals,d_e)}\}$ and $\WassersteinMulti{F}{d}{1,\dots,n}(t_1, t_2) := \inf \{\EvaluationFunctor{F}(d_1,\dots,d_n)(t) \mid t \in \Couplings{F}(t_1,t_2)\}$. This setting grants us access to new examples such as the product and the coproduct bifunctors.

[Product Bifunctor]exampleexa:product-bifunctor
For the product bifunctor $F\colon \Set^2 \to \Set$ where $F(X_1,X_2) = X_1\times X_2$ and $F(f_1,f_2) = f_1 \times f_2$ we consider the evaluation function $\max \colon\reals^2\to\reals$ and for fixed parameters $c_1,c_2 \in (0,1]$ and $p \in \mathbb{N}$ the function $\rho\colon\reals^2\to\reals$, $\rho(x_1,x_2) = (c_1 x_1^p + c_2 x_2^p)^{1/p}$. These functions are well-behaved, the Kantorovich-Rubinstein duality holds and the supremum [infimum] of the Kantorovich [Wasserstein] pseudometrics is always a maximum [minimum]. For the first function we obtain the $\infty$-product pseudometric $d_\infty((x_1,x_2),(y_1,y_2)) = \max(d_1(x_1,y_1), d_2(x_2,y_2))$ and for the other function the weighted $p$-product pseudometric $d_p((x_1,x_2),(y_1,y_2)) = (c_1d_1^p(x_1,y_1)+c_2d_2^p(x_2,y_2))^{1/p}$.

Note that the pseudometric space $(X_1\times X_2,d_\infty)$ is the usual binary (category theoretic) product of $(X_1,d_1)$ and $(X_2, d_2)$. Similarly, we can also obtain the binary coproduct.

[Coproduct Bifunctor]exampleexa:coproduct-bifunctor
For the coproduct bifunctor $F\colon \Set^2 \to \Set$, where $F(X_1,X_2) = X_1 + X_2 = X_1\times\{1\} \cup X_2\times \{2\}$ and $F(f_1,f_2) = f_1+f_2$ we take the evaluation function $\ev_F\colon\reals+\reals\to\reals$, $\ev_F(x,i) = x$. This function is well-behaved, the Kantorovich-Rubinstein duality holds and the supremum of the Kantorovich pseudometric is always a maximum whereas the infimum of the Wasserstein pseudometric is a minimum if and only if any coupling of the two elements exists. We obtain the coproduct pseudometric $d_+$ where $d_+((x_1,i_1), (x_2, i_2))$ is equal to $d_{i}(x_1,x_2)$ if $i_1=i_2=i$ and equal to $\top$ otherwise.

§ FINAL COALGEBRA & COALGEBRAIC BEHAVIORAL PSEUDOMETRICS

In this section we assume an arbitrary lifting $\LiftedFunctor{F}\colon \PMet \to \PMet$ of an endofunctor $F$ on $\Set$. For any pseudometric space $(X,d)$ we write
$\LiftedMetric{F}{d}$ for the pseudometric obtained by applying
$\LiftedFunctor{F}$ to $(X,d)$. Such a lifting can be obtained as
described earlier, but also by taking a lifted multifunctor and
fixing all parameters apart from one, or by the composition of such
functors. The following result ensures that if $\kappa \colon \Omega \to F
\Omega$ is a final $F$-coalgebra, then there is also a final
$\LiftedFunctor{F}$-coalgebra which is constructed by simply enriching
$\Omega$ with a pseudometric $d_\Omega$.

Let $\LiftedFunctor{F}\colon\PMet\to\PMet$ be a lifting of a functor $F\colon \Set\to \Set$ which has a final coalgebra $\kappa\colon\Omega\to F\Omega$. For every ordinal $i$ we construct a pseudometric $d_i\colon\Omega\times\Omega\to\reals$ as follows: $d_0 := 0$ is the zero pseudometric, $d_{i+1} := \LiftedMetric{F}{d_i}\circ(\kappa\times\kappa)$ for all ordinals $i$ and $d_j = \sup_{i<j} d_i$ for all limit ordinals $j$.
This sequence converges for some ordinal $\theta$, i.e $d_\theta = \LiftedMetric{F}{d_\theta}\circ(\kappa\times\kappa)$. Moreover $\kappa\colon (\Omega,d_\theta) \nonexpansiveTo (F\Omega,\LiftedMetric{F}{d_\theta})$ is the final $\LiftedFunctor{F}$-coalgebra.

We noted that for any set $X$, the set of pseudometrics over
$X$, with pointwise order, is a complete lattice. Moreover the lifting
$\LiftedFunctor{F}$ induces a monotone function $\_^F$ which maps any
pseudometric $d$ on $X$ to $d^F$ on $FX$. If, additionally, such function is
$\omega$-continuous, i.e., it preserves the supremum of
$\omega$-chains, the construction in thm:final-coalgebra
will converge in at most $\omega$ steps, i.e., $d_\theta=d_\omega$.
We show in prop:w-continuity that the liftings induced by the
finite powerset functor and the probability distribution functor with
finite support are $\omega$-continuous. The arguments used for
convergence here suggests a connection with the work
in <cit.>, which provides fixed-point results for metric
functors which are not locally contractive.

Beyond equivalences of states, in $\PMet$ we can measure the distance
of behaviors in the final coalgebra. More precisely, the
behavioral distance of two states $x,y\in X$ of some coalgebra
$\alpha\colon X\to FX$ is defined via the pseudometric
$\mathit{bd}(x,y) = d_\theta(\final{x},\final{y})$. Such distances can
be computed analogously to $d_\theta$ above, replacing $\kappa\colon
\Omega \to F \Omega$ by $\alpha$. This way we do not need to explore the entire final coalgebra (which might be too large) but can restrict to the interesting part.

Let the chain of the $d_i$ converge in $\theta$ steps and
$\overline{F}$ preserve isometries. Let furthermore $\alpha\colon X\to
FX$ be an arbitrary coalgebra. For all ordinals $i$ we define a pseudometric $e_i\colon X\times X\to \reals$ as follows: $e_0$ is the zero pseudometric, $e_{i+1} = e_i^F\circ (\alpha\times\alpha)$ for all ordinals $i$ and $e_j = \sup_{i<j} e_i$ for all limit ordinals $j$. Then we reach a fixed point after $\zeta \leq \theta$ steps, i.e. $e_\zeta = e_\zeta^F\circ (\alpha\times\alpha)$, such that $\mathit{bd} = e_\zeta$.

Since $d_\theta$ is a pseudometric, we have that if
$\final{x}=\final{y}$ then $bd(x,y)=0$. The other direction does
not hold in general: for this $d_\theta$ has to be a proper metric.
thm:d-omega-is-metric at the end of this section
provides sufficient conditions guaranteeing this property.

To this aim, we proceed by recalling the final coalgebra construction
via the final chain which was first presented in the dual setting
(free/initial algebra).

Let $\mathbf{C}$ be a category with terminal object $\one$ and limits of ordinal-indexed cochains. For any functor $F\colon \mathbf{C} \to \mathbf{C}$ the final chain consists of objects $W_i$ for all ordinals $i$ and connection morphisms $p_{i,j} \colon W_j \to W_i$ for all ordinals $i \le j$. The objects are defined as $W_0 := \one$, $W_{i+1} := FW_i$ for all ordinals $i$, and $W_j := \lim_{i<j}W_i$ for all limit ordinals $j$. The morphisms are determined by $p_{0,i} :=\ !\colon W_i \to \one$, $p_{i,i} = \id_{W_i}$ for all ordinals $i$, $p_{i+1,j+1} := Fp_{i,j}$ for all ordinals $i < j$ and if $j$ is a limit ordinal the $p_{i,j}$ are the morphisms of the limit cone. They satisfy $p_{i,k} = p_{i,j} \circ p_{j,k}$ for all ordinals $i \leq j \leq k$. We say that the chain converges in $\lambda$ steps if $p_{\lambda,\lambda+1}\colon W_{\lambda+1} \to W_\lambda$ is an iso.

This construction does not necessarily converge, but if it does, we get a final coalgebra.

Let $\mathbf{C}$ be a category with terminal object $\one$ and limits of ordinal-indexed cochains. If the final chain of a functor $F\colon \mathbf{C} \to \mathbf{C}$ converges in $\lambda$ steps then $p_{\lambda, \lambda+1}^{-1}\colon W_\lambda \to FW_\lambda$ is the final coalgebra.

We now show under which circumstances $d_\theta$ is a metric and how our construction relates to the construction of the final chain.

Let $\LiftedFunctor{F}\colon\PMet\to\PMet$ be a lifting of a functor $F\colon \Set\to \Set$ which has a final coalgebra $\kappa\colon\Omega\to F\Omega$. Assume that $\LiftedFunctor{F}$ preserves isometries and metrics, that the final chain for $F$ converges and the chain of the $d_i$ converges in $\theta$ steps. Then $d_\theta$ is a metric, i.e. for $x,y\in \Omega$ we have $d_\theta(x,y) = 0\iff x=y$.

We will now get back to the examples studied at the beginning of the paper (ex:probabilistic-1 and ex:metric-ts-1) and discuss in which sense they are instances of our framework.

[Probabilistic Transition System, revisited]

To model the behavioral distance from
ex:probabilistic-1 in our framework, we set $\top=1$
and proceed to lift the following three functors: we first consider
the identity functor $\Id$ with evaluation map $\ev_\Id\colon
[0,1]\to [0,1]$, $\ev_\Id(z) = c\cdot z$ in order to integrate the
discount (ex:kant-rubinst). Then, we take the
coproduct with the singleton metric space
(exa:coproduct-bifunctor).  The combination of the two
functors yields the discrete version of the refusal functor of
<cit.>, namely $\LiftedFunctor{R}(X,d) = (X+\mathbf{1},
\hat{d})$ where $\hat{d}$ is taken from
ex:probabilistic-1.  Finally, we lift the probability
distribution functor $\Distributions $ to obtain
$\LiftedFunctor{\Distributions }$
(exa:probability-distribution-functor). All functors
satisfy the Kantorovich-Rubinstein duality and preserve metrics.

It is readily seen that $\LiftedFunctor{\Distributions}(\LiftedFunctor{R}(X,d))=(\Distributions(X+1), \overline{d})$, where $\overline{d}$ is defined as in ex:probabilistic-1).  Then, the least solution of $d(x,y) = \overline{d}(\alpha(x),\alpha(y))$ can be computed as in thm:comp-dist.

[Metric Transition Systems, revisited]

To obtain propositional distances in metric transition systems we set $\top = \infty$. We also define, for the auxiliary functor $G$, an evaluation function $\ev_G\colon G(\prealinf,\dots,\prealinf)\to\prealinf$ with $\ev_G(u) = \max_{r\in\Sigma} u(r)$. Let $\LiftedFunctor{G}$ be the corresponding lifted functor. It can be shown, similarly to exa:product-bifunctor, that the Kantorovich-Rubinstein duality holds and metrics are preserved. We instantiate the given pseudometric spaces $(M_{r_i},d_{r_i})$ as parameters and obtain the functor $\LiftedFunctor{F}(X,d) = \LiftedFunctor{G}((M_{r_1},d_{r_1}),\dots,(M_{r_n},d_{r_n})) \times \LiftedFunctor{\PowersetFinite}(X,d)$ (for the lifting of the powerset functor see exa:hausdorff).  Then, via thm:comp-dist, we obtain exactly the least solution of (<ref>) in ex:metric-ts-1.

§ RELATED AND FUTURE WORK

The ideas for our framework are heavily influenced by work on quantitative variants of (bisimulation) equivalence of probabilistic systems. In that context at first Giacalone et al. <cit.> observed that probabilistic bisimulation <cit.> is too strong and therefore introduced a metric based on the notion of $\epsilon$-bisimulations.

Using a logical characterization of bisimulation for labelled Markov processes (LMP) <cit.>, Desharnais et al. defined a family of metrics between these LMPs <cit.> via functional expressions: if evaluated on a state of an LMP, such a functional expression measures the extent to which a formula is satisfied in that state.  A different, coalgebraic approach, which inspired ours, is used by van Breugel et al. <cit.>. As presented in more detail in the examples above, they define a pseudometric on probabilistic systems via the Kantorovich pseudometric for probability measures. Moreover, they show in <cit.> that this metric is related to the logical pseudometric by Desharnais et al.

Our framework provides a toolbox to determine behavioral distances for different types of transition systems modeled as coalgebras.

Moreover, the liftings introduced in this paper pave the way to extend several coalgebraic methods to reason about quantitative properties of systems.

For instance the bisimulation proof principle, which allows to check
behavioral equivalence, assumes a specific meaning in $\PMet$: every
coalgebra $\alpha \colon (X,d) \to \LiftedFunctor{F}(X,d)$
coinductively proves that the behavioral distance $bd$ of the
underlying $F$-coalgebra on $\Set$ is smaller or equal than $d$.
Indeed, since $\finalO$ is nonexpansive, $d \geq
d_\theta(\finalO,\finalO)=bd$. This principle, which has already
been stated in different formulations (see
e.g. <cit.>), can now be
enhanced via up-to techniques by exploiting the liftings
introduced in this paper and the coalgebraic understanding of such
enhancements given in <cit.>.

Since up-to techniques can exponentially improve algorithms for
equivalence-checking, we hope that they could also optimize some of
the algorithms for computing (or approximating) behavioral
distances <cit.>.
At this point, it is worth recalling that the Kantorovich-Rubinstein
duality has been exploited in <cit.> for defining one of these
algorithms: the characterization given by the Wasserstein metric
allows to reduce to linear programming.

Another line of research potentially stemming from our work concerns the
so-called abstract GSOS <cit.> which
provides abstract coalgebraic conditions ensuring compositionality of
behavioral equivalence (with respect to some operators). By taking our lifting to
$\PMet$, abstract GSOS guarantees the nonexpansiveness of behavioral
distance, a property that has captured the interest of several
researchers <cit.>.  The main
technical challenge would be to lift to $\PMet$ not only functors, but
also distributive laws. Lifting of distributive laws would also be
needed for defining linear behavioral distances, exploiting the
coalgebraic account of trace semantics based on Kleisli
categories <cit.>.

We finally observe that the chains of Theorems <ref> and <ref> can be understood in terms of fibrations along the lines of <cit.>. A detailed comparison with <cit.> can be found in app:fibrations.

The authors are grateful to Franck van Breugel, Neil Ghani and  Daniela Petrişan for several precious suggestions and inspiring discussions.
The second author acknowledges the support by project  .