diff --git a/animate/interpolation.py b/animate/interpolation.py
index 54d2ddb..fb15eef 100644
--- a/animate/interpolation.py
+++ b/animate/interpolation.py
@@ -105,7 +105,7 @@ def project(source, target_space, **kwargs):
     projection operator when applied to :class:`firedrake.cofunction.Cofunction`\s.
 
     For details on the approach for achieving boundedness through mass lumping and
-    post-processing, see :cite:`FPP+:2009`.
+    post-processing, see :cite:`Farrell:2009`.
 
     :arg source: the function to be transferred
     :type source: :class:`firedrake.function.Function` or
@@ -296,7 +296,7 @@ def _validate_matching_spaces(Vs, Vt):
 def clement_interpolant(source, target_space=None, boundary=False):
     r"""
     Compute the Clement interpolant of a :math:`\mathbb P0` source field, i.e. take the
-    volume average over neighbouring cells at each vertex. See :cite:`Cle:75`.
+    volume average over neighbouring cells at each vertex. See :cite:`Clement:1975`.
 
     :arg source: the :math:`\mathbb P0` source field
     :kwarg target_space: the :math:`\mathbb P1` space to interpolate into
diff --git a/animate/metric.py b/animate/metric.py
index 9a7ac25..8069433 100644
--- a/animate/metric.py
+++ b/animate/metric.py
@@ -799,9 +799,9 @@ def compute_isotropic_dwr_metric(
         Compute an isotropic metric from some error indicator using an element-based
         formulation.
 
-        The formulation is based on that presented in :cite:`CPB:13`. Note that
+        The formulation is based on that presented in :cite:`Carpio:2013`. Note that
         normalisation is implicit in the metric construction and involves the
-        `convergence_rate` parameter, named :math:`alpha` in :cite:`CPB:13`.
+        `convergence_rate` parameter, named :math:`alpha` in :cite:`Carpio:2013`.
 
         Whilst an element-based formulation is used to derive the metric, the result is
         projected into :math:`\mathbb P1` space, by default.
@@ -834,9 +834,9 @@ def compute_anisotropic_dwr_metric(
         r"""
         Compute an anisotropic metric from some error indicator, given a Hessian field.
 
-        The formulation used is based on that presented in :cite:`CPB:13`. Note that
+        The formulation used is based on that presented in :cite:`Carpio:2013`. Note that
         normalisation is implicit in the metric construction and involves the
-        `convergence_rate` parameter, named :math:`alpha` in :cite:`CPB:13`.
+        `convergence_rate` parameter, named :math:`alpha` in :cite:`Carpio:2013`.
 
         If a Hessian is not provided then an isotropic formulation is used.
 
@@ -936,7 +936,7 @@ def compute_weighted_hessian_metric(
         Compute a vertex-wise anisotropic metric from a list of error indicators, given
         a list of corresponding Hessian fields.
 
-        The formulation used is based on that presented in :cite:`PPP+:06`. It is
+        The formulation used is based on that presented in :cite:`Power:2006`. It is
         assumed that the error indicators have been constructed in the appropriate way.
 
         :arg error_indicators: list of error indicators
@@ -987,7 +987,7 @@ def determine_metric_complexity(H_interior, H_boundary, target, p, **kwargs):
     Solve an algebraic problem to obtain coefficients for the interior and boundary
     metrics to obtain a given metric complexity.
 
-    See :cite:`LDA:10` for details. Note that we use a slightly different formulation
+    See :cite:Loseille:2010` for details. Note that we use a slightly different formulation
     here.
 
     :arg H_interior: Hessian component from domain interior
diff --git a/demos/bubble_shear.py b/demos/bubble_shear.py
index bf54c10..66c27db 100644
--- a/demos/bubble_shear.py
+++ b/demos/bubble_shear.py
@@ -511,3 +511,9 @@ def adapt_metric_advection(mesh, t_start, t_end, c):
 # to maintain numerical stability (look at CFL condition).
 #
 # This demo can also be accessed as a `Python script <bubble_shear.py>`__.
+#
+#
+# .. rubric:: References
+#
+# .. bibliography::
+#    :filter: docname in docnames
diff --git a/demos/combining_metrics.py b/demos/combining_metrics.py
index adcd7d8..384ee4a 100644
--- a/demos/combining_metrics.py
+++ b/demos/combining_metrics.py
@@ -23,9 +23,8 @@
 # intersection of two metrics produces a metric whose associated ellipsoid is
 # the largest ellipsoid that can fit within the ellipsoids of the two metrics:
 #
-# .. figure:: combining_ellipse_intersection.jpg
-#    :figwidth: 70%
-#    :align: center
+# .. raw:: html
+#   :file: combining_ellipse_intersection.svg
 #
 # Below we first set up two metrics: Metric 1 asks for a medium resolution of
 # :math:`hm=0.025` in the left, and a coarse resolution of :math:`hc=0.1` in
diff --git a/demos/demo_references.bib b/demos/demo_references.bib
deleted file mode 100644
index be1e854..0000000
--- a/demos/demo_references.bib
+++ /dev/null
@@ -1,24 +0,0 @@
-@article{FPP+:2009,
-  title={Conservative interpolation between unstructured meshes via supermesh construction},
-  author={Farrell, Patrick E and Piggott, Matthew D and Pain, Christopher C and Gorman, Gerard J and Wilson, Cian R},
-  journal={Computer methods in applied mechanics and engineering},
-  volume={198},
-  number={33-36},
-  pages={2632--2642},
-  year={2009},
-  publisher={Elsevier}
-}
-
-@article{Barral:2016,
-  title={Anisotropic mesh adaptation in Firedrake with PETSc DMPlex},
-  author={Barral, Nicolas and Knepley, Matthew G and Lange, Michael and Piggott, Matthew D and Gorman, Gerard J},
-  journal={arXiv preprint arXiv:1610.09874},
-  year={2016}
-}
-
-@phdthesis{Wilson:2010,
-  title={Modelling multiple-material flows on adaptive unstructured meshes},
-  author={Wilson, Cian},
-  year={2010},
-  school={Imperial College London}
-}
\ No newline at end of file
diff --git a/demos/ping_pong.py b/demos/ping_pong.py
index 08e4757..5391bf1 100644
--- a/demos/ping_pong.py
+++ b/demos/ping_pong.py
@@ -175,7 +175,7 @@
 # :math:`\mathbb P1` fields specifically, it is possible to achieve this using 'mass
 # lumping'. Recall the linear system above. Lumping amounts to replacing the mass matrix
 # :math:`\underline{\mathbf{M}_B}` with a diagonal matrix, whose entries correspond to
-# the sums over the corresponding mass matrix rows. (See :cite:`FPP+:2009` for details.)
+# the sums over the corresponding mass matrix rows. (See :cite:`Farrell:2009` for details.)
 #
 # The resulting bounded projection operator can be used in Animate by passing
 # ``bounded=True`` to the ``project`` function. ::
@@ -247,7 +247,8 @@
 #
 # This demo can also be accessed as a `Python script <ping_pong.py>`__.
 #
+#
 # .. rubric:: References
 #
-# .. bibliography:: demo_references.bib
+# .. bibliography::
 #    :filter: docname in docnames