diff --git a/oqupy/bath.py b/oqupy/bath.py
index 4266297..eabdc34 100644
--- a/oqupy/bath.py
+++ b/oqupy/bath.py
@@ -71,7 +71,7 @@ def __init__(
             self._coupling_operator = tmp_coupling_operator
             self._unitary = np.identity(self._dimension)
         else:
-            w, v = np.linalg.eig(tmp_coupling_operator)
+            w, v = np.linalg.eigh(tmp_coupling_operator)
             self._coupling_operator = np.diag(w)
             self._unitary = v
             assert np.allclose(tmp_coupling_operator, \
diff --git a/tests/coverage/bath_test.py b/tests/coverage/bath_test.py
index b4d6fc4..cfd22b1 100644
--- a/tests/coverage/bath_test.py
+++ b/tests/coverage/bath_test.py
@@ -19,7 +19,6 @@
 
 from oqupy.bath import Bath
 from oqupy.bath_correlations import PowerLawSD
-from oqupy import operators
 
 def test_bath():
     wc = 4.0
@@ -79,6 +78,26 @@ def test_bath():
     assert np.array_equal(bath_C.west_degeneracy_map,[0,0,0,0])
     assert np.array_equal(bath_C.north_degeneracy_map,[0,0,0,0])
 
+# Checks that the diagonalisation of the bath also works for degenerate
+# coupling operators (see issue #165)
+def test_bath_degenerate_coupling_op():
+    degenerat_coupling_op = np.array(
+        [[0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
+        [1.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
+        [1.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
+        [0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
+        [1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j],
+        [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
+        [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
+        [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j]])
+    correlations = PowerLawSD(
+        alpha=0.02,
+        zeta=1,
+        cutoff=0.5,
+        cutoff_type="exponential",
+        temperature=0.5)
+    Bath(degenerat_coupling_op, correlations) # this failed before PR #166
+
 def test_bath_bad_input():
     wc = 4.0
     alpha = 0.3