  We present a Bayesian framework for joint and coherent analyses of
multimessenger binary neutron star signals. The method, implemented in our
bajes infrastructure, incorporates a joint likelihood for multiple datasets,
support for various semi-analytical kilonova models and numerical-relativity
(NR) informed relations for the mass ejecta, as well as a technique to include
and marginalize over modeling uncertainties. As a first application, we analyze
the gravitational-wave GW170817 and the kilonova AT2017gfo data. These results
are then combined with the most recent X-ray pulsars analyses of PSR J0030+0451
and PSR J0740+6620 to obtain EOS constraints.Various constraints on the
mass-radius diagram and neutron star properties are then obtained by resampling
over a set of ten million parametrized EOS built under minimal assumptions. We
find that a joint and coherent approach improves the inference of the extrinsic
parameters (distance) and, among the instrinc parameters, the mass ratio. The
inclusion of NR informed relations strongly improves over the case of using an
agnostic prior on the intrinsic parameters. Comparing Bayes factors, we find
that the two observations are better explained by the common source hypothesis
only by assuming NR-informed relations. These relations break some of the
degeneracies in the employed kN models. The EOS inference folding-in PSR
J0952-0607 minimum-maximum mass, PSR J0030+0451 and PSR J0740+6620 data
constrains, among other quantities, the neutron star radius to
$R_{1.4}={12.30}^{+0.81}_{-0.56}$ km ($R_{1.4}={13.20}^{+0.91}_{-0.90}$ km) and
the maximum mass to $M_{max}={2.28}^{+0.25}_{-0.17}~{\rm M_\odot}$
($M_{max}={2.32}^{+0.30}_{-0.19}~{\rm M_\odot}$) where the ST+PDT (PDT-U)
analysis of Vinciguerra et a (2023) for PSR J0030+0451 is employed. Hence, the
systematics on PSR J0030+0451 data reduction currently dominate the mass-radius
diagram constraints.
