diff --git a/doc/src/bond_gaussian.rst b/doc/src/bond_gaussian.rst
index affc39bc43..f23d3eabdc 100644
--- a/doc/src/bond_gaussian.rst
+++ b/doc/src/bond_gaussian.rst
@@ -27,6 +27,9 @@ The *gaussian* bond style uses the potential:
 
    E = -k_B T ln\left(\sum_{i=1}^{n} \frac{A_i}{w_i \sqrt{\pi/2}} exp\left( \frac{-(r-r_{i})^2}{w_i^2})\right) \right)
 
+This analytical form is a suitable potential for obtaining
+mesoscale effective force fields which can reproduce target atomistic distributions :ref:`(Milano) <Milano0> 
+
 The following coefficients must be defined for each bond type via the
 :doc:`bond_coeff <bond_coeff>` command as in the example above, or in
 the data file or restart files read by the :doc:`read_data <read_data>`
@@ -60,3 +63,8 @@ Default
 
 none
 
+----------
+
+.. _Milano0:
+
+**(Milano)** G. Milano, S. Goudeau, F. Müller-Plathe, J. Polym. Sci. B Polym. Phys. 43, 871 (2005).