diff --git a/.github/workflows/flake8.yml b/.github/workflows/flake8.yml
index 3b067cd2..efb81565 100644
--- a/.github/workflows/flake8.yml
+++ b/.github/workflows/flake8.yml
@@ -21,6 +21,6 @@ jobs:
       - name: flake8 Lint
         uses: py-actions/flake8@v1
         with:
-          ignore: "E501,W605"
+          ignore: "E501"
           exclude: "__init__.py, input/__init__.py"
           path: "pyerrors"
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 8ef8055e..6d95d672 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -30,6 +30,6 @@ pytest --cov=pyerrors --cov-report html
 ```
 The linter `flake8` is executed with the command
 ```
-flake8 --ignore=E501,W605 --exclude=__init__.py pyerrors
+flake8 --ignore=E501 --exclude=__init__.py pyerrors
 ```
 Please make sure that all tests are passed for a new pull requests.
diff --git a/pyerrors/obs.py b/pyerrors/obs.py
index 69af5cea..5cc33b69 100644
--- a/pyerrors/obs.py
+++ b/pyerrors/obs.py
@@ -1333,7 +1333,7 @@ def correlate(obs_a, obs_b):
 
 
 def covariance(obs, visualize=False, correlation=False, **kwargs):
-    """Calculates the covariance matrix of a set of observables.
+    r'''Calculates the covariance matrix of a set of observables.
 
     The gamma method has to be applied first to all observables.
 
@@ -1349,11 +1349,9 @@ def covariance(obs, visualize=False, correlation=False, **kwargs):
     Notes
     -----
     The covariance is estimated by calculating the correlation matrix assuming no autocorrelation and then rescaling the correlation matrix by the full errors including the previous gamma method estimate for the autocorrelation of the observables. For observables defined on a single ensemble this is equivalent to assuming that the integrated autocorrelation time of an off-diagonal element is equal to the geometric mean of the integrated autocorrelation times of the corresponding diagonal elements.
-    $$
-    \tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}
-    $$
+    $$\tau_{\mathrm{int}, ij}=\sqrt{\tau_{\mathrm{int}, i}\times \tau_{\mathrm{int}, j}}$$
     This construction ensures that the estimated covariance matrix is positive semi-definite (up to numerical rounding errors).
-    """
+    '''
 
     length = len(obs)
     cov = np.zeros((length, length))