diff --git a/doc/source/reference/kwant.physics.rst b/doc/source/reference/kwant.physics.rst
index 90128888cb7aa734a978a45e3f6ed0d099d48ff5..38de2fee310dc4cb0a1646205970ad22e1a8a6dc 100644
--- a/doc/source/reference/kwant.physics.rst
+++ b/doc/source/reference/kwant.physics.rst
@@ -10,4 +10,4 @@ Leads
 
    Bands
    modes
-   self_energy
+   selfenergy
diff --git a/doc/source/whatsnew/0.3.rst b/doc/source/whatsnew/0.3.rst
index d9c5add361365526fb92f0ebaea6521cb25c1ddd..5ee9c7ae9771a556cb13e5d074498b6d5fdd7093 100644
--- a/doc/source/whatsnew/0.3.rst
+++ b/doc/source/whatsnew/0.3.rst
@@ -28,6 +28,8 @@ Some renames
 ------------
 * site groups are now called site families.  This affects all the names that
   used to contain "group" or "groups".
+* ``self_energy`` has been renamed to ``selfenergy`` in all cases, most notably
+  in `kwant.physics.selfenergy`.
 * ``wave_func`` has been renamed to `~kwant.solvers.default.wave_function`,
 * ``MonatomicLattice`` has been renamed to `~kwant.lattice.Monatomic`,
 * ``PolyatomicLattice`` has been renamed to `~kwant.lattice.Polyatomic`.
diff --git a/examples/square.py b/examples/square.py
index 14e0bd272b49972eb3e27556e2a819b479ea53c8..dbfd7e821a09cc40b544bba229c1d42187112317 100644
--- a/examples/square.py
+++ b/examples/square.py
@@ -6,7 +6,7 @@ from __future__ import division
 import numpy as np
 from matplotlib import pyplot
 import kwant
-from kwant.physics.selfenergy import square_self_energy
+from kwant.physics.selfenergy import square_selfenergy
 
 __all__ = ['System' ]
 
@@ -16,8 +16,8 @@ class Lead(object):
         self.t = t
         self.potential = potential
 
-    def self_energy(self, fermi_energy):
-        return square_self_energy(self.width, self.t, self.potential,
+    def selfenergy(self, fermi_energy):
+        return square_selfenergy(self.width, self.t, self.potential,
                                   fermi_energy)
 
 class System(kwant.system.FiniteSystem):
diff --git a/examples/tests/test_square.py b/examples/tests/test_square.py
index e49d67fffcd208e99c58b5c468c349f4d956a9e6..9abfbbb1aad05377e33f9d0bbeb491369faf7fa3 100644
--- a/examples/tests/test_square.py
+++ b/examples/tests/test_square.py
@@ -26,12 +26,12 @@ def test_hamiltonian():
             assert_almost_equal(m, m_herm)
             assert_almost_equal(m_herm, sys.hamiltonian(j, i))
 
-def test_self_energy():
+def test_selfenergy():
     sys = square.System((2, 4), 1)
     for lead in xrange(len(sys.lead_interfaces)):
         n_orb = sum(
             sys.num_orbitals(site) for site in sys.lead_interfaces[lead])
-        se = sys.self_energy(lead, 0)
+        se = sys.selfenergy(lead, 0)
         assert_equal(len(se.shape), 2)
         assert_equal(se.shape[0], se.shape[1])
         assert_equal(se.shape[0], n_orb)
diff --git a/kwant/builder.py b/kwant/builder.py
index 53b3bd62f3c7d15a11ec3394c95c95a07c149870..38f0511895119263c5f817d694b7956382655d6f 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -388,10 +388,10 @@ class Lead(object):
         Notes
         -----
         The finalized lead must at least have a single method
-        ``self_energy(energy)`` but it can be a full
+        ``selfenergy(energy)`` but it can be a full
         `kwant.system.InfiniteSystem` as well.
 
-        The method ``self_energy`` of the finalized lead must return a square
+        The method ``selfenergy`` of the finalized lead must return a square
         matrix of appropriate size.
 
         The order of interface sites is assumed to be preserved during
@@ -444,21 +444,21 @@ class SelfEnergy(Lead):
 
     Parameters
     ----------
-    self_energy_func : function
+    selfenergy_func : function
         Function which returns the self energy matrix for the interface sites
         given the energy and optionally a list of extra arguments.
     interface : sequence of `Site` instances
     """
-    def __init__(self, self_energy_func, interface):
-        self._self_energy_func = self_energy_func
+    def __init__(self, selfenergy_func, interface):
+        self._selfenergy_func = selfenergy_func
         self.interface = tuple(interface)
 
     def finalized(self):
         """Trivial finalization: the object is returned itself."""
         return self
 
-    def self_energy(self, energy, args=()):
-        return self._self_energy_func(energy, args)
+    def selfenergy(self, energy, args=()):
+        return self._selfenergy_func(energy, args)
 
 
 class ModesLead(Lead):
@@ -482,7 +482,7 @@ class ModesLead(Lead):
     def modes(self, energy, args=()):
         return self._modes_func(energy, args)
 
-    def self_energy(self, energy, args=()):
+    def selfenergy(self, energy, args=()):
         modes = self.modes(energy, args)
         return physics.selfenergy(modes=modes)
         
diff --git a/kwant/physics/leads.py b/kwant/physics/leads.py
index 370e4ecac3db7962dd6a0ff470b038b469fbee57..a9a630e999fe6e4789fdda2c16246b5d3073a96b 100644
--- a/kwant/physics/leads.py
+++ b/kwant/physics/leads.py
@@ -17,7 +17,7 @@ from .. import linalg as kla
 
 dot = np.dot
 
-__all__ = ['self_energy', 'modes', 'Modes']
+__all__ = ['selfenergy', 'modes', 'Modes']
 
 Linsys = namedtuple('Linsys', ['eigenproblem', 'v', 'extract', 'project'])
 
@@ -503,7 +503,7 @@ def modes(h_onslice, h_hop, tol=1e6):
     return Modes(vecs, vecslmbdainv, nmodes, v)
 
 
-def self_energy(lead_modes, tol=1e6):
+def selfenergy(lead_modes, tol=1e6):
     """
     Compute the self-energy generated by a lead.
 
@@ -539,7 +539,7 @@ def self_energy(lead_modes, tol=1e6):
         return la.solve(vecslmbdainv.T, vecs.T).T
 
 
-def square_self_energy(width, hopping, fermi_energy):
+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 96054f8674903aafd87fd2042671e8114ef06eb5..e182fa8c2a646dfcbfd5cb3fbb8407a8cf4b82f4 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.self_energy(leads.modes(h, t))
+modes_se = lambda h, t: leads.selfenergy(leads.modes(h, t))
 
 def h_slice(t, w, e):
     h = (4 * t - e) * np.identity(w)
@@ -25,7 +25,7 @@ def test_analytic_numeric():
     t = 0.78                    # hopping element
     e = 1.3                     # Fermi energy
 
-    assert_almost_equal(leads.square_self_energy(w, t, e),
+    assert_almost_equal(leads.square_selfenergy(w, t, e),
                         modes_se(h_slice(t, w, e), -t * np.identity(w)))
 
 
@@ -48,8 +48,8 @@ def test_regular_fully_degenerate():
     h_onslice[w:, w:] = h_onslice_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
-    g[:w, :w] = leads.square_self_energy(w, t, e)
-    g[w:, w:] = leads.square_self_energy(w, t, e)
+    g[:w, :w] = leads.square_selfenergy(w, t, e)
+    g[w:, w:] = leads.square_selfenergy(w, t, e)
 
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
 
@@ -78,8 +78,8 @@ def test_regular_degenerate_with_crossing():
     h_onslice[w:, w:] = -h_onslice_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
-    g[:w, :w] = leads.square_self_energy(w, t, e)
-    g[w:, w:] = -np.conj(leads.square_self_energy(w, t, e))
+    g[:w, :w] = leads.square_selfenergy(w, t, e)
+    g[w:, w:] = -np.conj(leads.square_selfenergy(w, t, e))
 
     assert_almost_equal(g, modes_se(h_onslice, hop))
 
@@ -104,7 +104,7 @@ def test_singular():
     h_onslice[:w, w:] = h_hop_s
     h_onslice[w:, :w] = h_hop_s
     h_onslice[w:, w:] = h_onslice_s
-    g = leads.square_self_energy(w, t, e)
+    g = leads.square_selfenergy(w, t, e)
 
     print np.round(g, 5) / np.round(modes_se(h_onslice, h_hop), 5)
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
@@ -132,7 +132,7 @@ def test_singular_but_square():
     h_onslice[w:, w:] = h_onslice_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
-    g[:w, :w] = leads.square_self_energy(w, t, e)
+    g[:w, :w] = leads.square_selfenergy(w, t, e)
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
 
 def test_singular_fully_degenerate():
@@ -163,8 +163,8 @@ def test_singular_fully_degenerate():
     h_onslice[3*w:4*w, 3*w:4*w] = h_onslice_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
-    g[:w, :w] = leads.square_self_energy(w, t, e)
-    g[w:, w:] = leads.square_self_energy(w, t, e)
+    g[:w, :w] = leads.square_selfenergy(w, t, e)
+    g[w:, w:] = leads.square_selfenergy(w, t, e)
 
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
 
@@ -197,8 +197,8 @@ def test_singular_degenerate_with_crossing():
     h_onslice[3*w:4*w, 3*w:4*w] = -h_onslice_s
 
     g = np.zeros((2*w, 2*w), dtype=complex)
-    g[:w, :w] = leads.square_self_energy(w, t, e)
-    g[w:, w:] = -np.conj(leads.square_self_energy(w, t, e))
+    g[:w, :w] = leads.square_selfenergy(w, t, e)
+    g[w:, w:] = -np.conj(leads.square_selfenergy(w, t, e))
 
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
 
diff --git a/kwant/physics/tests/test_noise.py b/kwant/physics/tests/test_noise.py
index bf5f1d58f5a67ee9d58e1d25fce1b8a68c954d4f..4a9c25bf24c439801f74a8740ab315cfa4f71956 100644
--- a/kwant/physics/tests/test_noise.py
+++ b/kwant/physics/tests/test_noise.py
@@ -43,9 +43,9 @@ class LeadWithOnlySelfEnergy(object):
     def __init__(self, lead):
         self.lead = lead
 
-    def self_energy(self, energy, args=()):
+    def selfenergy(self, energy, args=()):
         assert args == ()
-        return self.lead.self_energy(energy)
+        return self.lead.selfenergy(energy)
 
 
 def test_twoterminal():
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index 55a095606007f426c98dfe0965e4c8dae836289f..7ec7b61ca66de46eb8bcf521d600577771df624a 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -224,7 +224,7 @@ class SparseSolver(object):
                     # See comment about zero-shaped sparse matrices at the top.
                     rhs.append(np.zeros((lhs.shape[1], 0)))
             else:
-                sigma = lead.self_energy(energy, args)
+                sigma = lead.selfenergy(energy, args)
                 lead_info.append(sigma)
                 indices = np.r_[tuple(slice(offsets[i], offsets[i + 1])
                                       for i in interface)]
diff --git a/kwant/solvers/tests/_test_sparse.py b/kwant/solvers/tests/_test_sparse.py
index 38460be5b3897d4320fcc5908a0c9eac79a971ee..16c6ac9007d66d2a06302d41bf9ab4a79d0c306a 100644
--- a/kwant/solvers/tests/_test_sparse.py
+++ b/kwant/solvers/tests/_test_sparse.py
@@ -21,9 +21,9 @@ class LeadWithOnlySelfEnergy(object):
     def __init__(self, lead):
         self.lead = lead
 
-    def self_energy(self, energy, args=()):
+    def selfenergy(self, energy, args=()):
         assert args == ()
-        return self.lead.self_energy(energy)
+        return self.lead.selfenergy(energy)
 
 
 # Test output sanity: that an error is raised if no output is requested,
@@ -273,7 +273,7 @@ def test_tricky_singular_hopping(solve):
 
 # Test equivalence between self-energy and scattering matrix representations.
 # Also check that transmission works.
-def test_self_energy(solve):
+def test_selfenergy(solve):
     np.random.seed(4)
     system = kwant.Builder()
     left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
@@ -316,7 +316,7 @@ def test_self_energy(solve):
     assert_almost_equal(t_should_be, sol.transmission(1, 0))
 
 
-def test_self_energy_reflection(solve):
+def test_selfenergy_reflection(solve):
     np.random.seed(4)
     system = kwant.Builder()
     left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
diff --git a/kwant/solvers/tests/test_mumps.py b/kwant/solvers/tests/test_mumps.py
index 01f7d6cc58ca762f412ef1ff849c7615ea9b9bfd..2b1b8b15d1467285a00d8aa443bc3fcf4fd2b827 100644
--- a/kwant/solvers/tests/test_mumps.py
+++ b/kwant/solvers/tests/test_mumps.py
@@ -81,19 +81,19 @@ def test_tricky_singular_hopping():
 
 
 @skipif(_no_mumps)
-def test_self_energy():
+def test_selfenergy():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_self_energy(solve)
+        _test_sparse.test_selfenergy(solve)
 
 
 @skipif(_no_mumps)
-def test_self_energy_reflection():
+def test_selfenergy_reflection():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_self_energy_reflection(solve)
+        _test_sparse.test_selfenergy_reflection(solve)
 
 
 @skipif(_no_mumps)
diff --git a/kwant/solvers/tests/test_sparse.py b/kwant/solvers/tests/test_sparse.py
index 9f3edd694d9d8101f0810250906d4fb6ad86daa9..46b58a16587af765bd5bf84ecb836e81c36b2847 100644
--- a/kwant/solvers/tests/test_sparse.py
+++ b/kwant/solvers/tests/test_sparse.py
@@ -39,12 +39,12 @@ def test_tricky_singular_hopping():
     _test_sparse.test_tricky_singular_hopping(solve)
 
 
-def test_self_energy():
-    _test_sparse.test_self_energy(solve)
+def test_selfenergy():
+    _test_sparse.test_selfenergy(solve)
 
 
-def test_self_energy_reflection():
-    _test_sparse.test_self_energy_reflection(solve)
+def test_selfenergy_reflection():
+    _test_sparse.test_selfenergy_reflection(solve)
 
 
 def test_very_singular_leads():
diff --git a/kwant/system.py b/kwant/system.py
index 5e557d2b5c789e1745af0ab613fc02dd34dfd23a..cfc1e413752e087f39aadf6b014a417c8bed972c 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -69,7 +69,7 @@ class FiniteSystem(System):
     ------------------
     leads : sequence of lead objects
         Each lead has to provide a method 
-        ``self_energy(energy, args)`` or ``modes(energy, args)``.
+        ``selfenergy(energy, args)`` or ``modes(energy, args)``.
     lead_interfaces : sequence of sequences of integers
         Each sub-sequence contains the indices of the system sites to which the
         lead is connected.
@@ -78,7 +78,7 @@ class FiniteSystem(System):
     -----
     The length of `leads` must be equal to the length of `lead_interfaces`.
 
-    For lead ``n``, the method leads[n].self_energy must return a square matrix
+    For lead ``n``, the method leads[n].selfenergy must return a square matrix
     whose size is ``sum(self.num_orbitals(neighbor)`` for neighbor in
     self.lead_interfaces[n])``. The output format for ``leads[n].modes`` is more
     complicated, and it should match the output of `~kwant.physics.modes`.
@@ -161,10 +161,10 @@ class InfiniteSystem(System):
         ham.flat[::ham.shape[0] + 1] -= energy
         return physics.modes(ham, self.inter_slice_hopping(args=args))
 
-    def self_energy(self, energy, args=()):
+    def selfenergy(self, energy, args=()):
         """Return self-energy of a lead.
 
         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.self_energy(self.modes(energy, args=()))
+        return physics.selfenergy(self.modes(energy, args=()))