diff --git a/doc/source/whatsnew/0.3.rst b/doc/source/whatsnew/0.3.rst
index 3ce75445831a44d4bf303b90694969495a5224f4..eaaf78ce93fbb80452ceff0c7dcc00bd85a1c85e 100644
--- a/doc/source/whatsnew/0.3.rst
+++ b/doc/source/whatsnew/0.3.rst
@@ -133,6 +133,3 @@ responsibility.
 The new class `~kwant.builder.ModesLead` allows to attach leads that have a
 custom way of calculating their modes (e.g. ideal leads) directly to a
 `~kwant.builder.Builder`.
-
-`~kwant.physics.selfenergy` now uses the output of `~kwant.physics.modes` as
-input instead of the slice Hamiltonian and inter-slice hopping.
diff --git a/kwant/builder.py b/kwant/builder.py
index 8dbc29a9f74d6e8f40d8771e7da1992d0860e6fe..6869c547e7af41f28346339cf1622dc403e735d3 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -481,7 +481,7 @@ class ModesLead(Lead):
 
     def selfenergy(self, energy, args=()):
         modes = self.modes(energy, args)
-        return physics.selfenergy(modes=modes)
+        return physics.selfenergy_from_modes(modes=modes)
 
 
 
diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
index a9a630e999fe6e4789fdda2c16246b5d3073a96b..dec942ef95b9c13dfda92ef11f317c50d1d522f4 100644
--- a/kwant/physics/leads.py
+++ b/kwant/physics/leads.py
@@ -503,20 +503,15 @@ def modes(h_onslice, h_hop, tol=1e6):
     return Modes(vecs, vecslmbdainv, nmodes, v)
 
 
-def selfenergy(lead_modes, tol=1e6):
+def selfenergy_from_modes(lead_modes):
     """
-    Compute the self-energy generated by a lead.
-
-    The lead is described either by the unit-cell Hamiltonian h_onslice and the
-    hopping matrix h_hop, or by its modes.
+    Compute the self-energy generated by lead modes.
 
     Parameters
     ----------
     lead_modes : Modes(vecs, vecslmbdainv, nmodes, v) a named tuple
         The modes in the lead, as calculated by
         `kwant.physics.modes`.
-    tol : double
-        Numerical tolerance used in several places.
 
     Returns
     -------
@@ -539,6 +534,35 @@ def selfenergy(lead_modes, tol=1e6):
         return la.solve(vecslmbdainv.T, vecs.T).T
 
 
+def selfenergy(h_onslice, h_hop, tol=1e6):
+    """
+    Compute the self-energy generated by the lead.
+
+    Parameters
+    ----------
+    h_onslice : numpy array, real or complex, shape (N,N) The unit cell
+        Hamiltonian of the lead slice.
+    h_hop : numpy array, real or complex, shape (N,M)
+        The hopping matrix from a lead slice to the one on which self-energy
+        has to be calculated (and any other hopping in the same direction).
+    tol : float
+        Tolerance for numerical error.
+
+    Returns
+    -------
+    Sigma : numpy array, real or complex, shape (M,M)
+        The computed self-energy. Note that even if `h_onslice` and `h_hop`
+        are both real, `Sigma` will typically be complex. (More precisely, if
+        there is a propagating mode, `Sigma` will definitely be complex.)
+
+    Notes
+    -----
+    For simplicity this function internally calculates the modes first.
+    This may cause a small slowdown, and can be improved if necessary.
+    """
+    return selfenergy_from_modes(modes(h_onslice, h_hop, tol))
+
+
 def square_selfenergy(width, hopping, fermi_energy):
     """
     Calculate analytically the self energy for a square lattice.
diff --git a/kwant/physics/tests/test_leads.py b/kwant/physics/tests/test_leads.py
index e182fa8c2a646dfcbfd5cb3fbb8407a8cf4b82f4..0ec651be4153635409cd93c9bcc8b22fd26ebf93 100644
--- a/kwant/physics/tests/test_leads.py
+++ b/kwant/physics/tests/test_leads.py
@@ -11,7 +11,7 @@ import numpy as np
 from numpy.testing import assert_almost_equal
 from kwant.physics import leads
 
-modes_se = lambda h, t: leads.selfenergy(leads.modes(h, t))
+modes_se = leads.selfenergy
 
 def h_slice(t, w, e):
     h = (4 * t - e) * np.identity(w)
diff --git a/kwant/system.py b/kwant/system.py
index fc5eb368179391585c184ad3958843b6201f4384..d07d54ed6fff84d80484fc3f04002fed16ccf708 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -148,10 +148,10 @@ class InfiniteSystem(System):
                                           sparse=sparse, args=args)
 
     def modes(self, energy, args=()):
-        """Return self-energy of a lead.
+        """Return mode decomposition of the lead
 
-        The returned matrix has the shape (n, n), where n is
-        ``sum(self.num_orbitals(i) for i in range(self.slice_size))``.
+        See documentation of `~kwant.physics.Modes` for the return
+        format details.
         """
         ham = self.slice_hamiltonian(args=args)
         shape = ham.shape
@@ -167,4 +167,10 @@ class InfiniteSystem(System):
         The returned matrix has the shape (n, n), where n is
         ``sum(self.num_orbitals(i) for i in range(self.slice_size))``.
         """
-        return physics.selfenergy(self.modes(energy, args=()))
+        ham = self.slice_hamiltonian(args=args)
+        shape = ham.shape
+        assert len(shape) == 2
+        assert shape[0] == shape[1]
+        # Subtract energy from the diagonal.
+        ham.flat[::ham.shape[0] + 1] -= energy
+        return physics.selfenergy(ham, self.inter_slice_hopping(args=args))