.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/topology/plot_diabolical_points.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_topology_plot_diabolical_points.py: Diabolical Points: Conical Intersections in a Two-Parameter Family ====================================================================== Berry and Wilkinson (1984) asked when two energy levels of a real symmetric Hamiltonian -- the quantum billiard in a triangle, in their case -- can be degenerate. The answer goes back to von Neumann and Wigner: a real symmetric :math:`2\times2` block has *two* conditions for degeneracy (equal diagonal elements, vanishing off-diagonal one), so without a symmetry one parameter is not enough. Levels **avoid crossing** along any one-parameter path, and meet only at isolated points of a two-parameter plane. There the two sheets of the spectrum touch like the tips of two cones -- a double cone, the "diabolo" that gave **diabolical points** their name. A diabolical point is also a source of geometry: the real eigenvector of either level changes **sign** after one loop around it, a Berry phase of :math:`\pi` (the same mechanism as around graphene's Dirac points, but here with no lattice symmetry pinning the degeneracy in place). This example builds the tight-binding version of a *scalene triangle*: a triangular flake whose hoppings along two bond directions, :math:`t_1` and :math:`t_2`, play the role of the triangle's two angles, with a small potential on one site that breaks every mirror symmetry. It finds a diabolical point, checks that it is a cone, that one-parameter paths miss it, that the eigenvector flips sign around it (:func:`tbkit.exceptional.encircle`), and that it *moves* -- rather than disappears -- when the potential changes. .. GENERATED FROM PYTHON SOURCE LINES 30-39 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt from scipy.optimize import minimize import tbkit.lattices as lattices from tbkit.system import System from tbkit.plot import Plot import tbkit.exceptional as ex .. GENERATED FROM PYTHON SOURCE LINES 40-47 A triangular flake with three bond directions ----------------------------------------------------- Ten sites of the triangular lattice cut into an equilateral triangle. The Hamiltonian is split by bond direction (0, 60 and 120 degrees), so that the hoppings along two of them, :math:`t_1, t_2`, can be tuned; the third stays at :math:`t = -1`. A potential :math:`V` on one site next to a corner removes the remaining mirror symmetries. .. GENERATED FROM PYTHON SOURCE LINES 47-76 .. code-block:: Python lat = lattices.triangular() lat.get_lattice(n1=4, n2=4) frac_i = lat.coor['x'] - lat.coor['y'] / np.sqrt(3) frac_j = 2 * lat.coor['y'] / np.sqrt(3) lat.remove_sites([int(i) for i in np.nonzero(np.rint(frac_i + frac_j) > 3)[0]]) assert lat.sites == 10 blocks = [] for ang in (0., 60., 120.): sys = System(lat) sys.set_hopping([{'n': 1, 'ang': ang, 't': -1.}]) sys.get_ham() blocks.append(sys.ham.toarray().real) defect = int(np.argmin(np.hypot(lat.coor['x'] - 1., lat.coor['y']))) assert np.allclose(sum(blocks), sum(blocks).T) def flake(p, v=0.3): '''H(t1, t2): real symmetric, no symmetry left for v != 0.''' ham = blocks[0] + p[0] * blocks[1] + p[1] * blocks[2] ham[defect, defect] += v return ham sys_plot = System(lat) sys_plot.set_hopping([{'n': 1, 't': -1.}]) fig_lat = Plot(sys_plot).lattice(plt_hop=True, ms=14, figsize=(4, 3.6)) .. image-sg:: /api/gallery/topology/images/sphx_glr_plot_diabolical_points_001.png :alt: plot diabolical points :srcset: /api/gallery/topology/images/sphx_glr_plot_diabolical_points_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 77-85 Finding the diabolical point of levels 5 and 6 ------------------------------------------------------- Levels are counted from 0, in increasing energy. The gap :math:`E_6 - E_5` is minimized over the plane (its square is smooth: a quadratic form at the cone). It reaches zero -- to machine precision -- at an isolated point well away from the three lines :math:`t_1 = 1`, :math:`t_2 = 1`, :math:`t_1 = t_2` on which the unperturbed triangle would have a mirror symmetry. .. GENERATED FROM PYTHON SOURCE LINES 85-105 .. code-block:: Python N = 5 def gap(p, v=0.3): en = np.linalg.eigvalsh(flake(p, v)) return en[N + 1] - en[N] def find_dp(v, guess): res = minimize(lambda p: gap(p, v) ** 2, guess, method='Nelder-Mead', options={'xatol': 1e-13, 'fatol': 1e-32, 'maxiter': 6000}) return res.x dp = find_dp(0.3, (1.15, 0.45)) print('diabolical point at (t1, t2) = ({:.6f}, {:.6f}), gap = {:.1e}'.format(*dp, gap(dp))) assert gap(dp) < 1e-10 assert min(abs(dp[0] - 1.), abs(dp[1] - 1.), abs(dp[0] - dp[1])) > 0.1 .. rst-class:: sphx-glr-script-out .. code-block:: none diabolical point at (t1, t2) = (1.162474, 0.462589), gap = 0.0e+00 .. GENERATED FROM PYTHON SOURCE LINES 106-112 A cone: the gap grows linearly in every direction ---------------------------------------------------------- Near a diabolical point the two levels are :math:`\bar E \pm |\mathbf{q}|` for some linear map :math:`\mathbf{q}` of the displacement, so the gap is proportional to the distance -- with a slope that depends on the direction (an elliptic cone). .. GENERATED FROM PYTHON SOURCE LINES 112-121 .. code-block:: Python for theta in np.linspace(0., np.pi, 4, endpoint=False): u = np.array([np.cos(theta), np.sin(theta)]) slopes = [gap(dp + r * u) / r for r in (1e-3, 1e-4)] print('direction {:5.1f} deg: gap/r = {:.5f} (r = 1e-3), {:.5f} (r = 1e-4)' .format(np.degrees(theta), *slopes)) assert slopes[0] > 0.01 assert np.isclose(slopes[0], slopes[1], rtol=1e-2) .. rst-class:: sphx-glr-script-out .. code-block:: none direction 0.0 deg: gap/r = 0.76804 (r = 1e-3), 0.76750 (r = 1e-4) direction 45.0 deg: gap/r = 1.04507 (r = 1e-3), 1.04569 (r = 1e-4) direction 90.0 deg: gap/r = 2.22458 (r = 1e-3), 2.22464 (r = 1e-4) direction 135.0 deg: gap/r = 2.10816 (r = 1e-3), 2.10817 (r = 1e-4) .. GENERATED FROM PYTHON SOURCE LINES 122-127 One parameter is not enough: avoided crossings ------------------------------------------------------ Scanning :math:`t_1` at fixed :math:`t_2` gives an avoided crossing unless the line happens to go through the diabolical point exactly. The smallest gap along the line grows with the distance by which it misses. .. GENERATED FROM PYTHON SOURCE LINES 127-136 .. code-block:: Python t1_scan = np.linspace(dp[0] - 0.25, dp[0] + 0.25, 1001) min_gaps = {} for offset in (0., 0.01, 0.03): gaps = np.array([gap((t1, dp[1] + offset)) for t1 in t1_scan]) min_gaps[offset] = gaps.min() print('t2 = t2* + {:.2f}: smallest gap along the line = {:.4f}'.format(offset, gaps.min())) assert min_gaps[0.] < 2e-3 < min_gaps[0.01] < min_gaps[0.03] .. rst-class:: sphx-glr-script-out .. code-block:: none t2 = t2* + 0.00: smallest gap along the line = 0.0000 t2 = t2* + 0.01: smallest gap along the line = 0.0041 t2 = t2* + 0.03: smallest gap along the line = 0.0112 .. GENERATED FROM PYTHON SOURCE LINES 137-145 Berry phase pi: the eigenvector changes sign ----------------------------------------------------- Carry the eigenvector of level 5 around a circle by parallel transport (:func:`tbkit.exceptional.encircle`). Around the diabolical point it comes back as :math:`-|u\rangle`: the phase is :math:`\pi`, for either level. Around a circle that misses it, it comes back unchanged. For a real symmetric matrix this is literally a sign change of the real eigenvector, which cannot be removed by any choice of phases. .. GENERATED FROM PYTHON SOURCE LINES 145-158 .. code-block:: Python around = ex.encircle(flake, ex.circle(dp, 0.05), band=N, nk=200) upper = ex.encircle(flake, ex.circle(dp, 0.05), band=N + 1, nk=200) missing = ex.encircle(flake, ex.circle(dp + (0.2, 0.), 0.05), band=N, nk=200) print('phase around the DP: {:.6f} (level 5), {:.6f} (level 6); around nothing: {:.2e}' .format(around.phase, upper.phase, missing.phase)) assert around.final_band == N # the eigenvalue comes back: no swap at a DP assert np.isclose(abs(around.phase), np.pi, atol=1e-8) assert np.isclose(abs(upper.phase), np.pi, atol=1e-8) assert np.isclose(missing.phase, 0., atol=1e-8) # the vorticity of the pair is 0 (unlike an exceptional point's 1/2) assert abs(ex.vorticity(flake, ex.circle(dp, 0.05), (N, N + 1))) < 1e-10 .. rst-class:: sphx-glr-script-out .. code-block:: none phase around the DP: -3.141593 (level 5), -3.141593 (level 6); around nothing: -2.04e-32 .. GENERATED FROM PYTHON SOURCE LINES 159-165 Not pinned by symmetry: the point moves with the potential -------------------------------------------------------------------- A degeneracy forced by symmetry sits on a symmetric line or point. This one is *accidental*: change :math:`V` and it moves continuously across the plane, and it cannot disappear by a small perturbation -- it can only move, until it meets another diabolical point. .. GENERATED FROM PYTHON SOURCE LINES 165-175 .. code-block:: Python track = {} guess = dp for v in (0.3, 0.35, 0.4, 0.45, 0.5): track[v] = find_dp(v, guess) guess = track[v] print('V = {:.2f}: DP at ({:.4f}, {:.4f}), gap {:.1e}'.format(v, *track[v], gap(track[v], v))) assert gap(track[v], v) < 1e-10 assert np.linalg.norm(track[0.5] - track[0.3]) > 0.01 .. rst-class:: sphx-glr-script-out .. code-block:: none V = 0.30: DP at (1.1625, 0.4626), gap 0.0e+00 V = 0.35: DP at (1.1905, 0.4722), gap 1.1e-16 V = 0.40: DP at (1.2184, 0.4823), gap 1.1e-16 V = 0.45: DP at (1.2462, 0.4925), gap 1.1e-16 V = 0.50: DP at (1.2735, 0.5029), gap 0.0e+00 .. GENERATED FROM PYTHON SOURCE LINES 176-178 The double cone, and the flip of the eigenvector ------------------------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 178-211 .. code-block:: Python fig = plt.figure(figsize=(12, 4.6)) ax1 = fig.add_subplot(1, 3, 1, projection='3d') r = np.linspace(0., 0.08, 25) phi = np.linspace(0., 2 * np.pi, 49) R, PHI = np.meshgrid(r, phi) T1, T2 = dp[0] + R * np.cos(PHI), dp[1] + R * np.sin(PHI) levels = np.array([[np.linalg.eigvalsh(flake((a, b)))[N:N + 2] for a, b in zip(ra, rb)] for ra, rb in zip(T1, T2)]) for n, color in ((0, 'tab:blue'), (1, 'tab:red')): ax1.plot_surface(T1, T2, levels[..., n], color=color, alpha=0.7, linewidth=0) ax1.set_xlabel('$t_1$') ax1.set_ylabel('$t_2$') ax1.set_zlabel('$E$') ax1.set_title('A diabolo: levels 5 and 6') ax2 = fig.add_subplot(1, 3, 2) for offset, style in ((0., 'k-'), (0.01, 'b--'), (0.03, 'g:')): ens = np.array([np.linalg.eigvalsh(flake((t1, dp[1] + offset)))[N:N + 2] for t1 in t1_scan]) ax2.plot(t1_scan, ens, style, lw=1.2) ax2.set_xlabel('$t_1$') ax2.set_ylabel('$E$') ax2.set_title('Crossing only through the DP') ax3 = fig.add_subplot(1, 3, 3) overlap = [np.real(np.vdot(around.vector[0], v) / np.vdot(around.vector[0], around.vector[0])) for v in around.vector] ax3.plot(around.s, overlap, 'k') ax3.axhline(0., color='gray', lw=0.5) ax3.set_xlabel('loop parameter $s$') ax3.set_ylabel(r'$\langle u(0)|u(s)\rangle$') ax3.set_title('Eigenvector returns as $-u$') fig.tight_layout() .. image-sg:: /api/gallery/topology/images/sphx_glr_plot_diabolical_points_002.png :alt: A diabolo: levels 5 and 6, Crossing only through the DP, Eigenvector returns as $-u$ :srcset: /api/gallery/topology/images/sphx_glr_plot_diabolical_points_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.682 seconds) .. _sphx_glr_download_api_gallery_topology_plot_diabolical_points.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_diabolical_points.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_diabolical_points.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_diabolical_points.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_