diff --git a/GUIDE.md b/GUIDE.md
index bc447b5..945a42e 100644
--- a/GUIDE.md
+++ b/GUIDE.md
@@ -984,7 +984,8 @@ The function is included in the `dpd` program,
 and the expected value of _P_ (labelled EOS below)
 is printed for comparison at the end.
 We give results obtained by both
-the Lowe thermostat (L) and the Shardlow algorithm (S).
+the [Lowe thermostat](https://doi.org/10.1209/epl/i1999-00365-x) (L) 
+and the [Shardlow algorithm](https://doi.org/10.1137/S1064827501392879) (S).
 We take the default values of _a_ ρ/T=75, and of other parameters not mentioned below.
 
  _T_   | ρ | _P_ (EOS) | _P_ (L)   | _P_ (S)
diff --git a/dpd_module.f90 b/dpd_module.f90
index 6918852..6b4285e 100644
--- a/dpd_module.f90
+++ b/dpd_module.f90
@@ -213,6 +213,7 @@ SUBROUTINE lowe ( box, temperature, gamma_step )
     REAL, INTENT(in) :: gamma_step  ! Pair selection probability = Gamma * timestep
 
     ! Apply pairwise Lowe-Andersen thermostat to velocities stored in array v
+    ! CP Lowe, Europhys Lett, 47, 145 (1999).
 
     ! It is assumed that positions in the array r are in units where box = 1
     ! and that the array ij contains a list of all np pairs within range
@@ -262,6 +263,7 @@ SUBROUTINE shardlow ( box, temperature, gamma_step )
     REAL, INTENT(in) :: gamma_step  ! Gamma * timestep
 
     ! Implements the Shardlow integration algorithm for velocities stored in array v
+    ! T Shardlow, SIAM J Sci Comput, 24, 1267 (2003).
 
     ! It is assumed that positions in the array r are in units where box = 1
     ! and that the array ij contains a list of all np pairs within range
diff --git a/python_examples/GUIDE.md b/python_examples/GUIDE.md
index e8c1de2..a3dec6a 100644
--- a/python_examples/GUIDE.md
+++ b/python_examples/GUIDE.md
@@ -890,7 +890,8 @@ The function is included in the `dpd` program,
 and the expected value of _P_ (labelled EOS below)
 is printed for comparison at the end.
 We give results obtained by both
-the Lowe thermostat (L) and the Shardlow algorithm (S).
+the [Lowe thermostat](https://doi.org/10.1209/epl/i1999-00365-x) (L) 
+and the [Shardlow algorithm](https://doi.org/10.1137/S1064827501392879) (S).
 Results below were obtained with runs of 10 blocks, each 10000 steps
 as for the Fortran examples (program default is 10×1000).
 We take the default values of _a_ ρ/T=75, and of other parameters not mentioned below.
diff --git a/python_examples/dpd_module.py b/python_examples/dpd_module.py
index 748c035..7eb7765 100644
--- a/python_examples/dpd_module.py
+++ b/python_examples/dpd_module.py
@@ -131,6 +131,8 @@ def lowe ( box, temperature, gamma_step, v, pairs ):
     (but necessary, since the operations are supposed to be sequential).
     """
 
+    # CP Lowe, Europhys Lett, 47, 145 (1999).
+    
     # It is assumed that the pairs list contains all pairs within range
 
     import numpy as np
@@ -157,6 +159,8 @@ def shardlow ( box, temperature, gamma_step, v, pairs ):
     (but necessary, since the operations are supposed to be sequential).
     """
 
+    # T Shardlow, SIAM J Sci Comput, 24, 1267 (2003).
+    
     # It is assumed that the pairs list contains all pairs within range
 
     import numpy as np