diff --git a/.codespell-ignores b/.codespell-ignores
index 5ba18fe76fc..16679f80d6c 100644
--- a/.codespell-ignores
+++ b/.codespell-ignores
@@ -1,3 +1,4 @@
 chec
 usera
 nd
+studi
diff --git a/docs/bibliography.rst b/docs/bibliography.rst
index f6f7bb8679a..2e42522b1ff 100644
--- a/docs/bibliography.rst
+++ b/docs/bibliography.rst
@@ -2,38 +2,5 @@
 References
 **********
 
-.. [chaudhuri93] B. B. Chaudhuri und P. Kundu. "Optimum circular fit
-    to weighted data in multi-dimensional space". Pattern
-    Recognition Letters 14.1 (1993), S. 1–6
 
-.. [chalmecalvet2013] R. Chalme-Calvet, M. de Naurois, J.-P. Tavernet
-    "Muon efficiency of the H.E.S.S" telescope. AtmoHEAD Workshop, 2013.
-    https://arxiv.org/pdf/1403.4550.pdf
-
-.. [denaurois2009] M de Naurois, L. Roland
-    "A high performance likelihood reconstruction of gamma-rays for imaging
-    atmospheric Cherenkov telescopes". Astroparticle Physics 32.5, (2009)
-    https://arxiv.org/abs/0907.2610
-
-.. [temme2016] Temme T.F., "On the hunt for photons: analysis of Crab Nebula
-    data obtainedby the first G-APD Cherenkov telescope" 2016
-    DOI: 10.17877/DE290R-17773
-
-.. [factcleaning] https://github.com/fact-project/fact-tools
-
-.. [ctatopleveldatamodel] K. Kosack et al, "Top-level Data Model", CTAO Computing
-    Department Internal Documentation, CTA-SPE-OSO-000000-0001, v1A.
-
-.. [gasparetto] Gasparetto T., "Development of a computing farm for Cloud
-    computing on GPU Development and optimisation of data-analysis
-    methodologies for the Cherenkov Telescope Array",
-    Ph.D. thesis, Academic Year 2018/2019.
-
-.. [lts_regression] P. J. Rousseeuw, K. van Driessen, "Computing LTS Regression for Large Data Sets",
-    Data Mining and Knowledge Discovery, 12, 29-45, DOI: 10.1007/s10618-005-0024-4
-
-.. [hillas_uncertainties] W. Wittek, "Image Parameters", MAGIC-TDAS internal note 02-03, (2002).
-
-.. [cta_sb_ob] Oya et al, Scheduling Block Data Model Specification, CTA-SPE-COM-000000-0003, Issue 1, Rev. c
-
-.. [cta_r1event] CTA R1/Event Data Model Specification, version 1 revision C (2022).
+.. bibliography::
diff --git a/docs/conf.py b/docs/conf.py
index 61e104d8cd7..af107516ae3 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -53,6 +53,7 @@
     "sphinx.ext.mathjax",
     "sphinx_automodapi.automodapi",
     "sphinx_automodapi.smart_resolver",
+    "sphinxcontrib.bibtex",
     "nbsphinx",
     "matplotlib.sphinxext.plot_directive",
     "numpydoc",
@@ -399,3 +400,10 @@ def setup(app):
     "iminuit": ("https://scikit-hep.org/iminuit/", None),
     "traitlets": ("https://traitlets.readthedocs.io/en/stable/", None),
 }
+
+
+# workaround for sphinx-gallery-conf not being cacheable and warning resulting in docs failure
+suppress_warnings = ["config.cache"]
+
+bibtex_bibfiles = ["references.bib"]
+bibtex_encoding = "utf8"
diff --git a/docs/references.bib b/docs/references.bib
new file mode 100644
index 00000000000..aa40c827f29
--- /dev/null
+++ b/docs/references.bib
@@ -0,0 +1,138 @@
+@article{chaudhuri93,
+  title={Optimum circular fit to weighted data in multi-dimensional space},
+  author={Chaudhuri, B. B. and Kundu, P.},
+  journal={Pattern Recognition Letters},
+  volume={14},
+  number={1},
+  pages={1--6},
+  year={1993},
+  doi={10.1016/0167-8655(93)90126-X},
+}
+
+@inproceedings{chalmecalvet2013,
+  title={Muon efficiency of the H.E.S.S. telescope},
+  author={Chalme-Calvet, R. and de Naurois, M. and Tavernet, J.-P.},
+  booktitle={Proceedings, Atmospheric Monitoring for High-Energy Astroparticle Detectors (AtmoHEAD)},
+  location={Saclay, France},
+  year={2013},
+  eprint={1403.4550},
+  archivePrefix={arXiv},
+  primaryClass={astro-ph.IM},
+}
+
+@article{denaurois2009,
+  title={A high performance likelihood reconstruction of gamma rays for imaging atmospheric Cherenkov telescopes},
+  author={De Naurois, M. and Rolland, L.},
+  journal={Astroparticle Physics},
+  volume={32},
+  number={5},
+  pages={231--252},
+  year={2009},
+  doi={10.1016/j.astropartphys.2009.09.001},
+}
+
+@phdthesis{phd-temme,
+  title={On the hunt for photons: analysis of Crab Nebula data obtained by the first G-APD Cherenkov telescope},
+  doi={10.17877/DE290R-17773},
+  author={Temme, Thomas Fabian},
+  school={TU Dortmund},
+  year={2016},
+}
+
+
+@misc{fact-tools,
+  title={FACT-Tools},
+  author={Brügge, K. and Nöthe, M. and Buß, J. and others},
+  url={https://github.com/fact-project/fact-tools},
+  year=2019,
+  version={v1.1.3},
+  doi={10.5281/zenodo.3529353},
+}%
+}
+
+
+@manual{ctao-top-level-data-model,
+  title={CTA Top-level Data Model Specification},
+  date={2022-04-30},
+  version={1b},
+  number={CTA-SPE-OSO-000000-0001},
+  organization={CTAO},
+  url={https://redmine.cta-observatory.org/dmsf/files/16325},
+}
+
+@phdthesis{phd-gasparetto,
+  author={Gasparetto, T.},
+  title={Development of a computing farm for Cloud computing on GPU - Development and optimisation of data-analysis methodologies for the Cherenkov Telescope Array},
+  year={2020},
+  school={Università degli Studi di Trieste},
+  url={https://arts.units.it/handle/11368/2963769},
+}
+
+@article{lts-regression,
+  title={Computing LTS regression for large data sets},
+  author={Rousseeuw, P. J. and Van Driessen, K.},
+  journal={Data mining and knowledge discovery},
+  volume={12},
+  pages={29--45},
+  year={2006},
+  doi={10.1007/s10618-005-0024-4},
+}
+
+%
+@techreport{hillas-uncertainties,
+  title={Image Parameters},
+  author={M. Wittek},
+  organization={MAGIC},
+  number={TDAS-02-03},
+  year={2002},
+  note={MAGIC internal note},
+}
+
+@manual{cta-sb-ob-data-model,
+  title={Scheduling Block Data Model Specification},
+  date={2022-12-20},
+  version={1e},
+  number={CTA-SPE-COM-000000-0003},
+  organization={CTAO},
+  url={https://redmine.cta-observatory.org/dmsf/files/15587},
+}
+
+@manual{ctao-r1-event-data-model,
+  title={CTAO R1/Event Data Model Specification},
+  date={2023-08-10},
+  version={2a},
+  number={CTA-SPE-COM-000000-0002},
+  organization={CTAO},
+  url={https://redmine.cta-observatory.org/dmsf/files/8627},
+}
+
+@article{muon-review,
+  title={Using muon rings for the calibration of the Cherenkov telescope array: a systematic review of the method and its potential accuracy},
+  author={Gaug, M. and Fegan, S. and Mitchell, A. M. W. and Maccarone, M. and Mineo, T. and Okumura, A.},
+  journal={The Astrophysical Journal Supplement Series},
+  volume={243},
+  number={1},
+  pages={11},
+  year={2019},
+  doi={10.3847/1538-4365/ab2123},
+}
+
+@article{hofmann-1999-comparison,
+  title={Comparison of techniques to reconstruct VHE gamma-ray showers from multiple stereoscopic Cherenkov images},
+  author={Hofmann, W. and Jung, I. and Konopelko, A. and Krawczynski, H. and Lampeitl, H. and Pühlhofer, G.},
+  journal={Astroparticle Physics},
+  volume={12},
+  number={3},
+  pages={135--143},
+  year={1999},
+  doi={10.1016/S0927-6505(99)00084-5},
+}
+
+
+@article{bright-star-catalog,
+  title={The Bright star catalogue},
+  author={Hoffleit, D. and Jaschek, C.},
+  journal={New Haven},
+  year={1991},
+  url={https://heasarc.gsfc.nasa.gov/W3Browse/star-catalog/bsc5p.html},
+}
diff --git a/environment.yml b/environment.yml
index f098cd53c1e..d632eb51d6e 100644
--- a/environment.yml
+++ b/environment.yml
@@ -39,6 +39,7 @@ dependencies:
   - sphinx-automodapi
   - pydata-sphinx-theme
   - sphinx-design
+  - sphinxcontrib-bibtex
   - tomli
   - towncrier
   - tqdm
diff --git a/examples/core/table_writer_reader.py b/examples/core/table_writer_reader.py
index f2f7c2f2e79..7a4ce4d1bce 100644
--- a/examples/core/table_writer_reader.py
+++ b/examples/core/table_writer_reader.py
@@ -230,6 +230,11 @@ def create_stream(n_event):
 
 table.attrs
 
+######################################################################
+# close the file
+#
+h5.close()
+
 
 ######################################################################
 # Reading one-row-at-a-time:
diff --git a/pyproject.toml b/pyproject.toml
index c0b8e248540..3aa09e3cf69 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -69,6 +69,7 @@ docs = [
     "numpydoc",
     "sphinx-design",
     "sphinx-gallery",
+    "sphinxcontrib-bibtex",
     "jupyter",
     "notebook",
     "graphviz",
diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py
index 1a76d141a1d..067cd8b6c6e 100644
--- a/src/ctapipe/containers.py
+++ b/src/ctapipe/containers.py
@@ -142,7 +142,7 @@ class CoordinateFrameType(enum.Enum):
 
 
 class EventType(enum.Enum):
-    """Enum of EventTypes as defined in the CTA Data Model [cta_r1event]_"""
+    """Enum of EventTypes as defined in :cite:p:`ctao-r1-event-data-model`"""
 
     # calibrations are 0-15
     FLATFIELD = 0
@@ -1424,13 +1424,14 @@ class ArrayEventContainer(Container):
 
 
 class SchedulingBlockContainer(Container):
-    """Stores information about the scheduling block. This is a simplified
-    version of the SB model, only storing what is necessary for analysis. From
-    [cta_sb_ob]_
+    """Stores information about the scheduling block.
 
+    This is a simplified version of the SB model, only storing what is necessary for analysis.
+    From :cite:p:`cta-sb-ob-data-model`.
     """
 
     default_prefix = ""
+
     sb_id = Field(UNKNOWN_ID, "Scheduling block ID", type=np.uint64)
     sb_type = Field(
         SchedulingBlockType.UNKNOWN,
diff --git a/src/ctapipe/image/cleaning.py b/src/ctapipe/image/cleaning.py
index c5845090254..6dc8ec54c62 100644
--- a/src/ctapipe/image/cleaning.py
+++ b/src/ctapipe/image/cleaning.py
@@ -333,8 +333,7 @@ def fact_image_cleaning(
 
     References
     ----------
-    See  [temme2016]_ and for implementation [factcleaning]_
-
+    See :cite:p:`phd-temme` and for implementation :cite:p:`fact-tools`.
     """
 
     # Step 1
diff --git a/src/ctapipe/image/muon/fitting.py b/src/ctapipe/image/muon/fitting.py
index 24dc21ebc88..2fce59f326f 100644
--- a/src/ctapipe/image/muon/fitting.py
+++ b/src/ctapipe/image/muon/fitting.py
@@ -12,7 +12,7 @@ def kundu_chaudhuri_circle_fit(x, y, weights):
     Fast and reliable analytical circle fitting method previously used
     in the H.E.S.S. experiment for muon identification
 
-    Implementation based on [chaudhuri93]_
+    Implementation based on :cite:p:`chaudhuri93`
 
     Parameters
     ----------
diff --git a/src/ctapipe/image/muon/intensity_fitter.py b/src/ctapipe/image/muon/intensity_fitter.py
index e169aeb42a4..b579a72a8c6 100644
--- a/src/ctapipe/image/muon/intensity_fitter.py
+++ b/src/ctapipe/image/muon/intensity_fitter.py
@@ -239,11 +239,8 @@ def image_prediction_no_units(
     min_lambda_m=300e-9,
     max_lambda_m=600e-9,
 ):
-    """Function for producing the expected image for a given set of trial
-    muon parameters without using astropy units but expecting the input to
-    be in the correct ones.
-
-    See [chalmecalvet2013]_
+    """
+    Unit-less version of `image_prediction`.
     """
 
     # First produce angular position of each pixel w.r.t muon center
@@ -431,6 +428,18 @@ def create_initial_guess(center_x, center_y, radius, telescope_description):
 
 
 class MuonIntensityFitter(TelescopeComponent):
+    """
+    Fit muon ring images with a theoretical model to estimate optical efficiency.
+
+    Function for producing the expected image for a given set of trial
+    muon parameters without using astropy units but expecting the input to
+    be in the correct ones.
+
+    The image prediction function is currently modeled after :cite:p:`chalmecalvet2013`.
+
+    For more information, also see :cite:p:`muon-review`.
+    """
+
     spe_width = FloatTelescopeParameter(
         help="Width of a single photo electron distribution", default_value=0.5
     ).tag(config=True)
diff --git a/src/ctapipe/image/pixel_likelihood.py b/src/ctapipe/image/pixel_likelihood.py
index 7ac0f3ecc44..3bc71e0e70c 100644
--- a/src/ctapipe/image/pixel_likelihood.py
+++ b/src/ctapipe/image/pixel_likelihood.py
@@ -1,7 +1,7 @@
 """
 Class for calculation of likelihood of a pixel expectation, given the pixel amplitude,
 the level of noise in the pixel and the photoelectron resolution.
-This calculation is taken from [denaurois2009]_.
+This calculation is taken from :cite:p:`denaurois2009`.
 
 The likelihood is essentially a poissonian convolved with a gaussian, at low signal
 a full possonian approach must be adopted, which requires the sum of contributions
@@ -52,7 +52,7 @@ class PixelLikelihoodError(RuntimeError):
 def neg_log_likelihood_approx(image, prediction, spe_width, pedestal):
     """Calculate negative log likelihood for telescope.
 
-    Gaussian approximation from [denaurois2009]_, p. 22 (equation between (24) and (25)).
+    Gaussian approximation from :cite:p:`denaurois2009`, p. 22 (equation between (24) and (25)).
 
     Simplification:
 
@@ -105,7 +105,7 @@ def neg_log_likelihood_numeric(
 ):
     """
     Calculate likelihood of prediction given the measured signal,
-    full numerical integration from [denaurois2009]_.
+    full numerical integration from :cite:p:`denaurois2009`.
 
     Parameters
     ----------
diff --git a/src/ctapipe/io/metadata.py b/src/ctapipe/io/metadata.py
index 655d9b45a7a..b4fd670144f 100644
--- a/src/ctapipe/io/metadata.py
+++ b/src/ctapipe/io/metadata.py
@@ -1,11 +1,11 @@
 """
-Management of CTA Reference Metadata, as defined in the CTA Top-Level Data Model
-document [ctatopleveldatamodel]_ , version 1A. This information is required to be
-attached to the header of any
-files generated.
+Management of CTA Reference Metadata.
+
+Definitions from :cite:`ctao-top-level-data-model`.
+This information is required to be attached to the header of any files generated.
 
 The class Reference collects all required reference metadata, and can be turned into a
-flat dictionary.  The user should try to fill out all fields, or use a helper to fill
+flat dictionary. The user should try to fill out all fields, or use a helper to fill
 them (as in `Activity.from_provenance()`)
 
 .. code-block:: python
diff --git a/src/ctapipe/reco/hillas_intersection.py b/src/ctapipe/reco/hillas_intersection.py
index 9304b94c259..3c8a1809edd 100644
--- a/src/ctapipe/reco/hillas_intersection.py
+++ b/src/ctapipe/reco/hillas_intersection.py
@@ -71,7 +71,7 @@ def _far_outside_fov(fov_lat, fov_lon):
 class HillasIntersection(HillasGeometryReconstructor):
     """
     This class is a simple re-implementation of Hillas parameter based event
-    reconstruction. e.g. https://arxiv.org/abs/astro-ph/0607333
+    reconstruction. See algorithm I of :cite:p:`hofmann-1999-comparison`.
 
     In this case the Hillas parameters are all constructed in the shared
     angular (Nominal) system. Direction reconstruction is performed by
@@ -80,14 +80,12 @@ class HillasIntersection(HillasGeometryReconstructor):
     reconstruction is performed by performing the same procedure in the
     tilted ground system.
 
-    The height of maximum is reconstructed by the projection os the image
+    The height of maximum is reconstructed by the projection of the image
     centroid onto the shower axis, taking the weighted average of all images.
 
     Uncertainties on the positions are provided by taking the spread of the
     crossing points, however this means that no uncertainty can be provided
     for multiplicity 2 events.
-
-    Note: only input from CameraFrame is currently supported
     """
 
     atmosphere_profile_name = traits.CaselessStrEnum(
diff --git a/src/ctapipe/reco/hillas_reconstructor.py b/src/ctapipe/reco/hillas_reconstructor.py
index 0f82e2bc51d..3906122047d 100644
--- a/src/ctapipe/reco/hillas_reconstructor.py
+++ b/src/ctapipe/reco/hillas_reconstructor.py
@@ -212,7 +212,7 @@ def __call__(self, event):
 
     def initialize_arrays(self, event, hillas_dict):
         """
-        Creates flat arrays of needed quantities from the event structure
+        Creates flat arrays of needed quantities from the event structure.
 
         Parameters
         ----------
@@ -227,8 +227,8 @@ def initialize_arrays(self, event, hillas_dict):
 
         Notes
         -----
-        The part of the algorithm taking into account divergent pointing
-        mode and the correction to the psi angle is explained in [gasparetto]_
+        The part of the algorithm taking into account divergent pointing mode and
+        the correction to the psi angle is explained in :cite:p:`phd-gasparetto`,
         section 7.1.4.
         """
         # get one telescope id to check what frame to use
@@ -378,10 +378,9 @@ def estimate_core_position(array_pointing, psi, positions):
 
         Notes
         -----
-        The part of the algorithm taking into account divergent pointing
-        mode and the usage of a corrected psi angle is explained in [gasparetto]_
+        The part of the algorithm taking into account divergent pointing mode and
+        the usage of a corrected psi angle is explained in :cite:p:`phd-gasparetto`
         section 7.1.4.
-
         """
 
         # Since psi has been recalculated in the fake CameraFrame
diff --git a/src/ctapipe/utils/astro.py b/src/ctapipe/utils/astro.py
index 4181f69f201..43ed10d1c64 100644
--- a/src/ctapipe/utils/astro.py
+++ b/src/ctapipe/utils/astro.py
@@ -12,13 +12,14 @@
 
 def get_bright_stars(pointing=None, radius=None, magnitude_cut=None):
     """
-    Returns an astropy table containing star positions above a given magnitude within
-    a given radius around a position in the sky, using the Yale bright star catalog
-    which needs to be present in the ctapipe-extra package. The included Yale bright
-    star catalog contains all 9096 stars, excluding the Nova objects present in the
-    original catalog from  Hoffleit & Jaschek (1991),
-    https://adsabs.harvard.edu/abs/1991bsc..book.....H, and is complete down to
-    magnitude ~6.5, while the faintest included star has mag=7.96.
+    Get an astropy table of bright stars.
+
+    This function is using the Yale bright star catalog, available through ctapipe
+    data downloads.
+
+    The included Yale bright star catalog contains all 9096 stars, excluding the
+    Nova objects present in the original catalog and is complete down to magnitude
+    ~6.5, while the faintest included star has mag=7.96. :cite:p:`bright-star-catalog`
 
     Parameters
     ----------