README.md

@@ -2,7 +2,7 @@
 
 Solver for Usadel equations in 2D with SOC gauge fields.
 
-Uses a [discretization with local gauge invariance](doc/discretization.md).
+Uses a [discretization with local gauge invariance](doc/discretization.rst).
 
 ## Requirements
 
@@ -16,9 +16,19 @@ Assuming prerequisites are installed:
 
 `python3 -m pip install .`
 
-See also https://mesonbuild.com/Quick-guide.html#compiling-a-meson-project
+## Documentation
+
+Built during installation. Can be found in the installed package
+directory, or opened with:
+
+`python3 -c 'import usadelndsoc; usadelndsoc.docs()'`
+
+## Bibliography
+
+P. Virtanen, "Numerical saddle point finding for diffusive superconductor actions",
+https://arxiv.org/abs/2406.19894 (2024).
 
 ## Acknowledgments
 
 This work was supported by European Union's HORIZON-RIA programme (Grant
-Agreement No.~101135240 JOGATE).
+Agreement No. 101135240 JOGATE).

doc/api.rst

-===
-API
-===
+===========
+Library API
+===========
 
 .. tableofcontents::
 

doc/conf.py

-project = 'usadelndsoc'
-copyright = '2024, Pauli Virtanen'
-author = 'Pauli Virtanen'
-release = '0.1'
+import os, sys
 
-extensions = ['sphinx.ext.autodoc', 'sphinx.ext.napoleon']
-exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+project = "usadelndsoc"
+copyright = "2024, Pauli Virtanen"
+author = "Pauli Virtanen"
+release = "0.1"
 
-html_theme = 'nature'
+extensions = ["sphinx.ext.autodoc", "sphinx.ext.napoleon"]
+exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"]
+
+html_theme = "nature"
 
 autodoc_member_order = "groupwise"
 
-autodoc_default_options = {
-    'special-members': '__init__, __call__'
-}
+autodoc_default_options = {"special-members": "__init__, __call__"}
+
+
+#
+# Setup module loader for autodoc when building docs from meson
+#
+loader_modules = [("usadelndsoc", os.environ.get("SPHINX_MESON_SOURCE_DIR"))]
+loader_ext = {"usadelndsoc._core": os.environ.get("SPHINX_MESON_EXT_PATH")}
+
+if loader_modules[0][1] is not None:
+    import importlib.abc, importlib.machinery
+
+    class MetaFinder(importlib.abc.MetaPathFinder):
+        @classmethod
+        def find_spec(self, fullname, path, target=None):
+            for k, v in loader_modules:
+                if fullname == k or fullname.startswith(k + "."):
+                    if v:
+                        source_dir = v
+                        break
+            else:
+                return None
+
+            if loader_ext.get(fullname):
+                loader = importlib.machinery.ExtensionFileLoader
+                origin = loader_ext[fullname]
+                is_package = False
+            else:
+                loader = importlib.machinery.SourceFileLoader
+                origin = os.path.join(source_dir, *fullname.split(".")[1:])
+                is_package = os.path.isdir(origin)
+                origin += "/__init__.py" if is_package else ".py"
+
+            if not os.path.exists(origin):
+                return None
+
+            print(
+                f"IMPORT: load {fullname} from {origin} with {loader}", file=sys.stderr
+            )
+
+            loader = loader(fullname, origin)
+            spec = importlib.machinery.ModuleSpec(
+                fullname, loader, is_package=is_package, origin=origin
+            )
+            if is_package:
+                spec.submodule_search_locations = [os.path.dirname(origin)]
+                spec.has_location = True
+            return spec
+
+    sys.meta_path.insert(0, MetaFinder)
+
+import usadelndsoc, usadelndsoc.core

doc/examples.rst

+Examples
+========
+
+Usage example scripts from under `examples/` folder of the project:
+
+Supercurrent diode effect in Rashba superconductor
+--------------------------------------------------
+
+.. literalinclude:: ../examples/cpr_sns.py
+
+Anomalous current in bulk
+-------------------------
+
+.. literalinclude:: ../examples/bulk_j_an.py
+

doc/index.rst

@@ -6,13 +6,29 @@
 usadelndsoc
 ===========
 
+Solver for Usadel equations in 2D with SOC gauge fields.
+
 .. toctree::
    :maxdepth: 1
    :caption: Contents:
 
    discretization
+   examples
    api
 
+
+Bibliography
+------------
+
+P. Virtanen, "Numerical saddle point finding for diffusive superconductor actions",
+https://arxiv.org/abs/2406.19894 (2024).
+
+Acknowledgments
+---------------
+
+This work was supported by European Union's HORIZON-RIA programme (Grant
+Agreement No. 101135240 JOGATE).
+
 Indices and tables
 ==================
 

doc/meson.build

@@ -11,6 +11,7 @@ endif
 doc_files = files(
   'index.rst',
   'discretization.rst',
+  'examples.rst',
   'api.rst',
   'conf.py',
 )
@@ -24,10 +25,12 @@ doc_html = custom_target('doc',
     '@CURRENT_SOURCE_DIR@',     # source dir
     '@OUTPUT@',                 # output dir
   ],
+  env: {'SPHINX_MESON_EXT_PATH': usadelndsoc_pymod.full_path(),
+        'SPHINX_MESON_SOURCE_DIR': meson.project_source_root() / 'src' / 'usadelndsoc'},
   depend_files: [doc_files, usadelndsoc_sources],
   depends: [usadelndsoc_pymod],
   output: 'html',
   build_by_default: true,
   install: true,
-  install_dir: get_option('datadir') / 'doc' / 'usadelndsoc'
+  install_dir: py3.get_install_dir(subdir : 'usadelndsoc' / 'doc', pure: false)
 )

examples/bulk_j_an.py

@@ -95,7 +95,6 @@ def j_anom(T, h, n=5, alpha_soc=0.01, perturbative=False):
         gm = (ww + 1j * h) / np.sqrt((ww + 1j * h) ** 2 + 1)
         fp = 1 / np.sqrt((ww - 1j * h) ** 2 + 1)
         fm = 1 / np.sqrt((ww + 1j * h) ** 2 + 1)
-        # XXX: check factors of 2
         Jy_an = (gp - gm) * (1 + gp * gm + fp * fm) * (-1j * alpha_soc**3) * eta
         Jantot += pi * aa * Jy_an
 

examples/cpr_sns.py

 # -*- coding:utf-8; eval: (blacken-mode) -*-
 """
-Calculate anomalous current in bulk case
+Calculate SNS junction supercurrent diode effect
 """
 import joblib
 import numpy as np
@@ -156,7 +156,7 @@ def do(W_xi=12, multin=True):
         return (abs(Jp) - abs(Jm)) / (abs(Jp) + abs(Jm))
 
     ax = axs[1]
-    ax.plot(L / xi, 100 * eff(Jxs[0]))
+    ax.plot(L / xi, 100 * eff(Jxs[0]), "*-")
     if multin:
         ax.plot(L / xi, 100 * eff(Jxs[1]), "k:", alpha=0.25)
         ax.plot(L / xi, 100 * eff(Jxs[2]), "k:")

examples/edelst.py

+# -*- coding:utf-8; eval: (blacken-mode) -*-
+"""
+Calculate spin density in finite-size 2D Rashba superconductor case, generated by the inverse spin-galvanic effect.
+"""
+import os
+import joblib
+import numpy as np
+from numpy import pi
+import matplotlib.pyplot as plt
+from scipy.linalg import expm
+import logging
+from numpy import exp, pi
+
+import usadelndsoc
+import usadelndsoc.bcs
+from usadelndsoc.matsubara import get_matsubara_sum
+from usadelndsoc.solver import *
+from usadelndsoc.util import vectorize_parallel
+from usadelndsoc.plotting import plot_Delta_J_2d
+
+from rashba2 import Gamma_to_alpha
+
+
+mem = joblib.Memory("cache")
+
+
+def get_solver(soc_alpha, eta, Lx=10, Ly=10, nx=5, ny=5):
+    sol = Solver(nx=nx, ny=ny)
+    sol.mask[...] = MASK_NONE
+    sol.mask[0, :] = MASK_TERMINAL
+    sol.mask[-1, :] = MASK_TERMINAL
+
+    sol.Lx = Lx
+    sol.Ly = Ly
+
+    sol.Omega[...] = 0
+
+    sol.D = 1.0
+    sol.eta = eta
+
+    nx, ny = sol.shape
+    dx = sol.Lx / nx
+    dy = sol.Ly / ny
+
+    Ax = soc_alpha * np.kron(S_0, S_y)
+    Ay = -soc_alpha * np.kron(S_0, S_x)
+
+    sol.Ux[...] = expm(1j * dx * Ax)
+    sol.Uy[...] = expm(1j * dy * Ay)
+
+    return sol
+
+
+@vectorize_parallel(returns_object=True, noarray=True)
+def solve(T, phi, L_S=5, L=5, W=None, n=5, Gamma_DP=10, Gamma_ST=0.4, **solver_kw):
+    usadelndsoc.logger.setLevel(-20)
+
+    if W is None:
+        W = L
+
+    T_c0 = np.exp(np.euler_gamma) / np.pi
+    L = 2 * L_S + L
+
+    nx = n
+    ny = int(round(n * W / L))
+    nS = int(round(nx * L_S / L))
+
+    alpha_soc, eta = Gamma_to_alpha(Gamma_DP, Gamma_ST)
+
+    Delta0 = usadelndsoc.bcs.BCS_Delta(T, T_c0)
+
+    sol = get_solver(soc_alpha=alpha_soc, eta=eta, Lx=L, Ly=W, nx=nx, ny=ny)
+
+    T_c = np.zeros(sol.shape, dtype=float)
+    T_c[:nS] = T_c0
+    T_c[-nS:] = T_c0
+
+    sol.Delta[nS:-nS] = 0
+    sol.Delta[:nS] = np.eye(2) * Delta0 * np.exp(1j * phi / 2)
+    sol.Delta[-nS:] = np.eye(2) * Delta0 * np.exp(-1j * phi / 2)
+
+    solver_kw.setdefault("workers", max(1, os.cpu_count() // 2))
+
+    sol.self_consistency(T, T_c, **solver_kw)
+    # sol.matsubara_sums(T, 0, workers=max(1, os.cpu_count() // 2), **solver_kw)
+
+    sol.reset()
+    return sol
+
+
+class g_an:
+    """
+    Analytical solution for Q and spin density in infinite 2D strip of width L,
+    generated by inverse spin-galvanic effect due to charge supercurrent.
+    """
+
+    def __init__(self, Delta, alpha, A_x, L=2, eta=0.1):
+        self.Delta = Delta
+        self.alpha = alpha
+        self.L = L
+        self.eta = eta
+        self.A_x = A_x
+
+    def _coefs(self, omega):
+        Delta = self.Delta
+        alpha = self.alpha
+        L = self.L
+        eta = self.eta
+        Ax = self.A_x
+
+        E = np.sqrt(omega**2 + Delta**2)
+        kE = np.sqrt(2 * E)
+        kp = np.sqrt(
+            kE**2 - 2 * alpha**2 + 2j * alpha * np.sqrt(7 * alpha**2 + 4 * kE**2)
+        )
+        km = np.sqrt(
+            kE**2 - 2 * alpha**2 - 2j * alpha * np.sqrt(7 * alpha**2 + 4 * kE**2)
+        )
+
+        ap = -4 * alpha * kp / (kp**2 - 4 * alpha**2 - kE**2)
+        am = -4 * alpha * km / (km**2 - 4 * alpha**2 - kE**2)
+
+        Qby3 = -4j * eta * Ax * alpha**3 * Delta**2 / (2 * alpha**2 + E) / E**2
+        Qby1 = -Qby3 * omega / Delta
+
+        # [Qz, Qy]
+        cu = np.array(
+            [
+                [kp - 2 * alpha * ap, km - 2 * alpha * am],
+                [kp * ap + 2 * alpha, km * am + 2 * alpha],
+            ]
+        )
+        cd = np.array(
+            [
+                [-kp + 2 * alpha * ap, -km + 2 * alpha * am],
+                [kp * ap + 2 * alpha, km * am + 2 * alpha],
+            ]
+        )
+
+        f0 = Delta / E
+
+        JHyz2 = -16 * Ax * alpha**2 * f0
+        JHyzrhs = -eta / 4 * JHyz2
+
+        QJyz1 = -1j * JHyzrhs * omega / E
+        QJyz3 = +1j * JHyzrhs * Delta / E
+
+        v1 = np.array([-2j * alpha * 1j * Qby1 + QJyz1, 0 * E])
+        v3 = np.array([-2j * alpha * 1j * Qby3 + QJyz3, 0 * E])
+
+        v2 = np.array([0 * E, 0 * E])
+
+        v = np.array([v1, v2, v3])
+
+        cu = np.moveaxis(cu, [0, 1], [-2, -1])
+        cd = np.moveaxis(cd, [0, 1], [-2, -1])
+        v = np.moveaxis(v, 1, -1)
+
+        bu = np.linalg.solve(cu, v[..., None])[..., 0]
+        bd = np.linalg.solve(cd, v[..., None])[..., 0]
+
+        bu = np.moveaxis(bu, -1, 1)
+        bd = np.moveaxis(bd, -1, 1)
+        k = np.array([kp, km])
+        a = np.array([ap, am])
+
+        g03 = omega / E
+        g01 = Delta / E
+
+        Qby = np.array([Qby1, 0 * Qby1, Qby3])
+
+        return g03, g01, k, a, bu, bd, Qby
+
+    def _Q(self, x, omega):
+        x, omega = np.broadcast_arrays(x, omega)
+
+        g03, g01, k, a, bu, bd, Qby = self._coefs(omega)
+
+        vd = bd * np.exp(-k * (x + self.L / 2))
+        vu = bu * np.exp(k * (x - self.L / 2))
+
+        dQz = vd.sum(axis=1) + vu.sum(axis=1)
+        dQy = (vd * (-a)).sum(axis=1) + (vu * a).sum(axis=1)
+
+        return dQy, dQz
+
+    def _Dg(self, x, omega, bulk=False):
+        x, omega = np.broadcast_arrays(x, omega)
+
+        g03, g01, k, a, bu, bd, Qby = self._coefs(omega)
+
+        sl = (Ellipsis, None, None)
+
+        t1y = np.kron(S_x, S_y)
+        t2y = np.kron(S_y, S_y)
+        t3y = np.kron(S_z, S_y)
+
+        t1z = np.kron(S_x, S_z)
+        t2z = np.kron(S_y, S_z)
+        t3z = np.kron(S_z, S_z)
+
+        dQy, dQz = self._Q(x, omega)
+        if bulk:
+            dQy *= 0
+            dQz *= 0
+
+        Q = (
+            +(Qby[0] + dQy[0])[sl] * t1y
+            + (Qby[1] + dQy[1])[sl] * t2y
+            + (Qby[2] + dQy[2])[sl] * t3y
+            + (0 + dQz[0])[sl] * t1z
+            + (0 + dQz[1])[sl] * t2z
+            + (0 + dQz[2])[sl] * t3z
+        )
+
+        return Q
+
+    def g0(self, omega):
+        g03, g01, k, a, bu, bd, Qby = self._coefs(omega)
+
+        sl = (Ellipsis, None, None)
+        t1 = np.kron(S_x, S_0)
+        t3 = np.kron(S_z, S_0)
+        return t3 * g03[sl] + t1 * g01[sl]
+
+    def g(self, x, omega, bulk=False):
+        x, omega = np.broadcast_arrays(x, omega)
+        return self.g0(omega) + self._Dg(x, omega, bulk=bulk)
+
+    def S(self, x, T, bulk=False):
+        x = np.asarray(x)
+
+        E_0 = 50 * abs(self.Delta)
+        omega, a = get_matsubara_sum(T, E_0)
+
+        x = np.asarray(x)[..., None]
+
+        g = a[..., :, None, None] * self.g(x, omega, bulk=bulk)
+        s = 0.5j * np.pi * g.sum(axis=-3)
+
+        sx = tr(s @ np.kron(S_z, S_x)).real
+        sy = tr(s @ np.kron(S_z, S_y)).real
+        sz = tr(s @ np.kron(S_z, S_z)).real
+
+        return 0.5 * np.moveaxis(np.array([sx, sy, sz]), 0, -1)
+
+
+def plot_fig6(save=False):
+    """
+    Plot paper Fig. 6 for the inverse spin-galvanic effect.
+    """
+    import usadelndsoc.bcs
+    from matplotlib import colors
+
+    Gamma_ST = 0.4
+    alpha, eta = Gamma_to_alpha(Gamma_DP=10, Gamma_ST=Gamma_ST)
+
+    xi = 1
+
+    T_c0 = np.exp(np.euler_gamma) / np.pi
+    Delta01 = usadelndsoc.bcs.BCS_Delta(T=0.1, Tc=T_c0)
+    Delta03 = usadelndsoc.bcs.BCS_Delta(T=0.3, Tc=T_c0)
+    Delta04 = usadelndsoc.bcs.BCS_Delta(T=0.4, Tc=T_c0)
+
+    sol01 = solve(T=0.1, phi=0.05, n=90, L_S=2, L=0, W=5, workers=5, mem=mem).item()
+    sol03 = solve(T=0.3, phi=0.05, n=90, L_S=2, L=0, W=5, workers=5, mem=mem).item()
+    sol04 = solve(T=0.4, phi=0.05, n=90, L_S=2, L=0, W=5, workers=5, mem=mem).item()
+
+    sol = sol01
+
+    A_x = 0.05 / sol.Lx
+    print(A_x)
+
+    g01 = g_an(Delta=Delta01, alpha=alpha, A_x=A_x, L=sol.Ly, eta=eta)
+    g03 = g_an(Delta=Delta03, alpha=alpha, A_x=A_x, L=sol.Ly, eta=eta)
+    g04 = g_an(Delta=Delta04, alpha=alpha, A_x=A_x, L=sol.Ly, eta=eta)
+
+    vs = [(0.1, sol01, g01), (0.3, sol03, g03), (0.4, sol04, g04)]
+
+    y = np.linspace(-sol01.Ly / 2, sol01.Ly / 2, 303)
+
+    plt.clf()
+
+    plt.gcf().set_size_inches(1.2 * 86 / 25.4, 86 / 25.4 * 0.6)
+    plt.gcf().set_layout_engine("constrained")
+
+    ax1 = plt.subplot(121)
+
+    S0 = xi * A_x * Gamma_ST
+
+    for j, (T, sol, g) in enumerate(vs):
+        plt.plot(y / xi, sol.result.S_tot_s(0, y).real[:, 1] / S0, "-", label=str(T))
+        plt.plot(y / xi, g.S(y, T=T).real[:, 1] / S0, "k:", label="__nolabel__")
+
+    plt.ylabel(r"$S_y / (\nu \Gamma_{st} A_x \xi_0)$")
+    plt.xlabel(r"$y/\xi_0$")
+
+    ax2 = plt.subplot(122)
+
+    offs = [-0.05, 0, 0.05]
+
+    for j, (T, sol, g) in enumerate(vs):
+        plt.plot(
+            y / xi,
+            sol.result.S_tot_s(0, y).real[:, 2] / S0 + offs[j],
+            "-",
+            label=str(T),
+        )
+        plt.plot(
+            y / xi, g.S(y, T=T).real[:, 2] / S0 + offs[j], "k:", label="__nolabel__"
+        )
+
+    plt.ylabel(r"$S_z / (\nu \Gamma_{st} A_x \xi_0)$")
+    plt.xlabel(r"$y/\xi_0$")
+    plt.legend(
+        ncols=3,
+        handlelength=0.7,
+        handletextpad=0.4,
+        columnspacing=0.4,
+        loc="lower left",
+        title=r"$T / \Delta_0$",
+    )
+
+    if save:
+        plt.savefig("../figs/edelst.pdf")
+
+
+def main():
+    plot_fig6(save=True)
+
+
+if __name__ == "__main__":
+    main()

examples/rashba2.py

+# -*- coding:utf-8; eval: (blacken-mode) -*-
+"""
+Calculate spin density and anomalous current in finite-size 2D Rashba superconductor case under exchange field.
+"""
+import os
+import joblib
+import numpy as np
+from numpy import pi
+import matplotlib.pyplot as plt
+from scipy.linalg import expm
+import logging
+from numpy import exp, pi
+
+import usadelndsoc
+from usadelndsoc.matsubara import get_matsubara_sum
+from usadelndsoc.solver import *
+from usadelndsoc.util import vectorize_parallel
+from usadelndsoc.plotting import plot_Delta_J_2d
+
+
+mem = joblib.Memory("cache")
+
+
+def get_solver(soc_alpha, eta, Lx=10, Ly=10, h=0.5, D=1.0, nx=5, ny=5):
+    sol = Solver(nx=nx, ny=ny)
+    sol.mask[...] = MASK_NONE
+    sol.Lx = Lx
+    sol.Ly = Ly
+
+    sol.Omega[...] = 0
+    sol.Delta[:, :] = np.eye(2)
+
+    sol.D = D
+    sol.eta = eta
+
+    nx, ny = sol.shape
+    dx = sol.Lx / nx
+    dy = sol.Ly / ny
+
+    Ax = soc_alpha * np.kron(S_0, S_y)
+    Ay = -soc_alpha * np.kron(S_0, S_x)
+
+    sol.Ux[...] = expm(1j * dx * Ax)
+    sol.Uy[...] = expm(1j * dy * Ay)
+
+    sol.Omega[...] += h * np.kron(S_z, S_x)
+
+    return sol
+
+
+def Gamma_to_alpha(Gamma_DP, Gamma_ST):
+    r"""
+    Transform from (Gamma_DP, Gamma_ST) to (alpha_soc, eta).
+
+    .. math::
+
+       \Gamma_{DP} &= 2 \alpha^2 p_F^2 \tau \\
+       \Gamma_{ST} &= \Gamma_{DP} \frac{\alpha\tau}{\xi_0}
+
+    and :math:`\xi_0^2 = D/\Delta_0`.
+
+    The length unit is :math:`L_0 = \sqrt{D/E_0}` where
+    :math:`E_0` is the energy unit.
+    """
+
+    alpha_soc = (Gamma_DP / 4) ** 0.5
+    eta = Gamma_ST / (4 * alpha_soc**3)
+    return alpha_soc, eta
+
+
+@vectorize_parallel(returns_object=True, noarray=True)
+def j_anom(T, h, L=5, W=None, n=5, Gamma_DP=10, Gamma_ST=0.4, **solver_kw):
+    usadelndsoc.logger.setLevel(1)
+
+    if W is None:
+        W = L
+
+    T_c0 = 0.556
+
+    nx = n
+    ny = int(round(n * W / L))
+
+    alpha_soc, eta = Gamma_to_alpha(Gamma_DP, Gamma_ST)
+
+    sol = get_solver(soc_alpha=alpha_soc, eta=eta, h=h, D=1.0, Lx=L, Ly=W, nx=nx, ny=ny)
+    sol.self_consistency(T, T_c0, workers=max(1, os.cpu_count() // 2), **solver_kw)
+
+    sol.reset()
+    return sol
+
+
+class g_an:
+    """
+    Analytical solution for Q and spin density in infinite 2D strip of length L,
+    generated by exchange field.
+    """
+
+    def __init__(self, Delta, alpha, hx, hy, hz, L=2, eta=0.1):
+        self.Delta = Delta
+        self.h = (hx, hy, hz)
+        self.alpha = alpha
+        self.L = L
+        self.eta = eta
+
+    def _coefs(self, omega):
+        Delta = self.Delta
+        hx, hy, hz = self.h
+        alpha = self.alpha
+        L = self.L
+
+        E = np.sqrt(omega**2 + Delta**2)
+        kE = np.sqrt(2 * E)
+        kp = np.sqrt(
+            kE**2 - 2 * alpha**2 + 2j * alpha * np.sqrt(7 * alpha**2 + 4 * kE**2)
+        )
+        km = np.sqrt(
+            kE**2 - 2 * alpha**2 - 2j * alpha * np.sqrt(7 * alpha**2 + 4 * kE**2)
+        )
+
+        ap = -4 * alpha * kp / (kp**2 - 4 * alpha**2 - kE**2)
+        am = -4 * alpha * km / (km**2 - 4 * alpha**2 - kE**2)
+
+        Qbz = -hz / (2 * E + 8 * alpha**2)
+        Qbx = -hx / (2 * E + 4 * alpha**2)
+        Qby = -hy / (2 * E + 4 * alpha**2)
+
+        cr = np.array(
+            [
+                [kp - 2 * alpha * ap, km - 2 * alpha * am],
+                [kp * ap + 2 * alpha, km * am + 2 * alpha],
+            ]
+        )
+        cl = np.array(
+            [
+                [-kp + 2 * alpha * ap, -km + 2 * alpha * am],
+                [kp * ap + 2 * alpha, km * am + 2 * alpha],
+            ]
+        )
+        v = np.array([2 * alpha * Qbx, -2 * alpha * Qbz])
+
+        cl = np.moveaxis(cl, [0, 1], [-2, -1])
+        cr = np.moveaxis(cr, [0, 1], [-2, -1])
+        v = np.moveaxis(v, 0, -1)
+
+        bl = np.linalg.solve(cl, v)
+        br = np.linalg.solve(cr, v)
+
+        bl = np.moveaxis(bl, -1, 0)
+        br = np.moveaxis(br, -1, 0)
+        k = np.array([kp, km])
+        a = np.array([ap, am])
+
+        g03 = omega / E
+        g01 = Delta / E
+
+        return g03, g01, k, a, bl, br, Qbx, Qby, Qbz
+
+    def _Q(self, x, omega):
+        x, omega = np.broadcast_arrays(x, omega)
+
+        g03, g01, k, a, bl, br, Qbx, Qby, Qbz = self._coefs(omega)
+
+        vl = bl * np.exp(-k * (x + self.L / 2))
+        vr = br * np.exp(k * (x - self.L / 2))
+
+        dQz = vl.sum(axis=0) + vr.sum(axis=0)
+        dQx = (vl * (-a)).sum(axis=0) + (vr * a).sum(axis=0)
+
+        return dQx, dQz
+
+    def _dQ(self, x, omega):
+        x, omega = np.broadcast_arrays(x, omega)
+
+        g03, g01, k, a, bl, br, Qbx, Qby, Qbz = self._coefs(omega)
+
+        vl = bl * np.exp(-k * (x + self.L / 2)) * (-k)
+        vr = br * np.exp(k * (x - self.L / 2)) * k
+
+        dQz = vl.sum(axis=0) + vr.sum(axis=0)
+        dQx = (vl * (-a)).sum(axis=0) + (vr * a).sum(axis=0)
+        return dQx, dQz
+
+    def g0(self, omega):
+        g03, g01, k, a, bl, br, Qbx, Qby, Qbz = self._coefs(omega)
+
+        sl = (Ellipsis, None, None)
+        t3 = np.kron(S_z, S_0)
+        t2 = np.kron(S_y, S_0)
+        return t3 * g03[sl] + t2 * g01[sl]
+
+    def _Dg(self, x, omega, der=False):
+        x, omega = np.broadcast_arrays(x, omega)
+
+        g03, g01, k, a, bl, br, Qbx, Qby, Qbz = self._coefs(omega)
+
+        sl = (Ellipsis, None, None)
+
+        t3 = np.kron(S_z, S_0)
+        g0 = self.g0(omega)
+
+        if not der:
+            dQx, dQz = self._Q(x, omega)
+            f = 1
+        else:
+            dQx, dQz = self._dQ(x, omega)
+            f = 0
+
+        Q = (
+            (Qbx * f + dQx)[sl] * np.kron(S_0, S_x)
+            + (Qby * f + 0)[sl] * np.kron(S_0, S_y)
+            + (Qbz * f + dQz)[sl] * np.kron(S_0, S_z)
+        )
+
+        return 1j * g0 @ (t3 @ g0 - g0 @ t3) @ Q
+
+    def g(self, x, omega, bulk=False):
+        x, omega = np.broadcast_arrays(x, omega)
+        if bulk:
+            return self.g0(omega)
+        return self.g0(omega) + self._Dg(x, omega)
+
+    def dg(self, x, omega):
+        x, omega = np.broadcast_arrays(x, omega)
+        return self._Dg(x, omega, der=True)
+
+    def S(self, x, T, der=False, bulk=False):
+        x = np.asarray(x)
+
+        E_0 = 10 * max(np.linalg.norm(self.h), abs(self.Delta))
+        omega, a = get_matsubara_sum(T, E_0)
+
+        x = np.asarray(x)[..., None]
+
+        if not der:
+            g = a[..., :, None, None] * self.g(x, omega, bulk=bulk)
+        else:
+            g = a[..., :, None, None] * self.dg(x, omega, bulk=bulk)
+        s = -1j * np.pi * g.sum(axis=-3)
+
+        sx = tr(s @ np.kron(S_z, S_x)).real
+        sy = tr(s @ np.kron(S_z, S_y)).real
+        sz = tr(s @ np.kron(S_z, S_z)).real
+
+        return 0.5 * np.moveaxis(np.array([sx, sy, sz]), 0, -1)
+
+    def dS(self, x, T):
+        return self.S(x, T, der=True)
+
+    def jy_an_t(self, x, T):
+        S = self.S(x, T)
+        dS = self.dS(x, T)
+        jy = self.dS[..., 2] - 2 * self.alpha * S[..., 0]
+        return jy * 2 * self.alpha**2 * self.eta
+
+
+def solve(n, **kw):
+    """
+    Solve for spin-galvanic effect
+    """
+    h = 0.1
+    T = 0.2
+
+    sol = j_anom(T, h, n=n, mem=mem, **kw).item()
+    nx, ny = sol.shape
+
+    res = sol.result
+
+    dx = sol.x[1] - sol.x[0]
+    dy = sol.y[1] - sol.y[0]
+
+    x = np.linspace(sol.x[0] - dx / 2, sol.x[-1] + dx / 2, 8 * nx)
+    y = np.linspace(sol.y[0] - dy / 2, sol.y[-1] + dy / 2, 8 * ny)
+
+    I = res.J_tot_c(x[:, None], y[None, :], method="bilinear").real
+    S = res.S_tot_s(x[:, None], y[None, :]).real
+
+    plt.clf()
+
+    plt.gcf().set_size_inches(8, 6)
+
+    plt.subplot(211)
+    plot_Delta_J_2d(sol.x, sol.y, res.Delta.A[:, :, 0, 0], x, y, I[..., 0], I[..., 1])
+    plt.axis("equal")
+
+    plt.subplot(212)
+    plt.pcolormesh(x, y, S[..., 0].T, rasterized=True)
+    plt.axis("equal")
+
+    plt.savefig("rashba2.png")
+
+    return sol
+
+
+def plot_osc(
+    sol, save=False, T=0.2, Gamma_DP=10, Gamma_ST=1, h=0.1, vert=True, Slog=True
+):
+    import usadelndsoc.bcs
+    from matplotlib import colors
+
+    x = np.linspace(-sol.Lx / 2, sol.Lx / 2, 307)[1:-1]
+    y = np.linspace(-sol.Ly / 2, sol.Ly / 2, 303)
+
+    res = sol.result
+    j = res.J_tot_c(x[:, None], y[None, :]).real
+    s = res.S_tot_s(x[:, None], y[None, :]).real
+    ja = np.hypot(j[..., 0].real, j[..., 1].real)
+    s0 = res.S_tot_s(x[:, None], [0]).real
+    s00 = res.S_tot_s([0], [0]).real[0, 0]
+
+    xi = 1
+
+    jamin = 10 ** (np.floor(np.log10(np.quantile(ja.ravel() / ja.max(), 0.05))))
+
+    T_c0 = 0.556
+    Delta = usadelndsoc.bcs.BCS_Delta(T, T_c0, h)
+    alpha, eta = Gamma_to_alpha(Gamma_DP, Gamma_ST)
+    g = g_an(Delta=Delta, alpha=alpha, hx=-h, hy=0, hz=0, L=sol.Lx, eta=eta)
+
+    plt.clf()
+
+    plt.gcf().set_size_inches(1.2 * 86 / 25.4, 86 / 25.4 * 0.6)
+    plt.gcf().set_layout_engine("constrained")
+
+    if vert:
+        ax1 = plt.subplot(211)
+        plt.tick_params(axis="both", labelbottom=False)
+    else:
+        ax1 = plt.subplot(121)
+
+    plt.streamplot(
+        x / xi,
+        y / xi,
+        j[..., 0].T,
+        j[..., 1].T,
+        density=1.25,
+        linewidth=0.75,
+        arrowsize=0.6,
+        cmap="plasma",
+        color=ja.T / ja.max(),
+        norm=colors.LogNorm(vmin=jamin, vmax=1),
+    )
+
+    L = sol.Lx / xi
+    W = sol.Ly / xi
+    plt.plot([-L / 2, L / 2], [W / 2, W / 2], c="b", ls=":", lw=1)
+    plt.plot([-L / 2, L / 2], [-W / 2, -W / 2], c="b", ls=":", lw=1)
+    plt.plot([-L / 2, -L / 2], [-W / 2, W / 2], c="b", ls=":", lw=1)
+    plt.plot([L / 2, L / 2], [-W / 2, W / 2], c="b", ls=":", lw=1)
+
+    plt.ylabel(r"$y/\xi_0$")
+    plt.axis("equal")
+    cb = plt.colorbar()
+    plt.gca().set_title(r"$j/j_{\rm max}$")
+
+    if not vert:
+        plt.xlabel(r"$x/\xi_0$")
+
+    if vert:
+        ax2 = plt.subplot(212, sharex=ax1)
+    else:
+        ax2 = plt.subplot(122)
+
+    plt.axvline(-L / 2, ls=":", c="b")
+    plt.axvline(+L / 2, ls=":", c="b")
+
+    s2 = g.S(x, T=T)
+    s20 = g.S(0, T=T)[..., 0]
+    s22 = g.S(0, T=T)[..., 2]
+
+    if Slog:
+        plt.semilogy(x / xi, abs(s0[:, 0, 0] / s00 - 1), "-")
+        plt.semilogy(x / xi, abs(s2[:, 0] / s20 - 1), "k:")
+        plt.ylabel(r"$|S_x(x,0)-S_{x,\rm bulk}|/S_{x,\rm bulk}$")
+        plt.ylim(1e-5, 1)
+    else:
+        plt.plot(x / xi, (s0[:, 0, 0] / s00 - 1), "-", label="$S_x$")
+        plt.plot(x / xi, (s2[:, 0] / s20 - 1), "k:", label="__nolabel__")
+        plt.plot(x / xi, (s0[:, 0, 2] / s00), "-", label="$S_z$")
+        plt.plot(x / xi, (s2[:, 2] / s00), "k:", label="__nolabel__")
+        plt.ylabel(r"$[S_i(x,0)-S_{i,\rm bulk}]/S_{x,\rm bulk}$")
+        plt.legend()
+
+    plt.xlabel(r"$x/\xi_0$")
+
+    plt.xlim(1.05 * x[0] / xi, 1.05 * x[-1] / xi)
+
+    kE = np.sqrt(2j * pi * T)
+    kp = np.sqrt(kE**2 - 2 * alpha**2 + 2j * alpha * np.sqrt(7 * alpha**2 + 4 * kE**2))
+    km = np.sqrt(kE**2 - 2 * alpha**2 - 2j * alpha * np.sqrt(7 * alpha**2 + 4 * kE**2))
+    print(2 * pi / (xi * kp), 2 * pi / (xi * km))
+
+    if save:
+        plt.savefig("osc.pdf")
+
+
+def plot_fig2(save=False):
+    """
+    Plot paper Fig. 3 for the spin-galvanic effect in a square
+    """
+    sol = solve(n=60, tol=1e-8)
+    plot_osc(sol, vert=False, Slog=False)
+    if save:
+        plt.savefig("../figs/oscsq.pdf")
+    return sol
+
+
+def plot_fig5(save=False):
+    """
+    Plot paper Fig. 5 for the spin-galvanic effect in long strip
+    """
+    sol = solve(n=220, L=10, W=1, tol=1e-8)
+    plot_osc(sol)
+    if save:
+        plt.savefig("../figs/osc.pdf")
+    return sol
+
+
+def plot_fig4(sol=None, scl=9, save=False):
+    import pyvista as pv
+    import matplotlib.pyplot as plt
+    from PIL import Image
+
+    if sol is None:
+        sol = solve(n=60, tol=1e-8)
+
+    x = np.linspace(sol.x[0], sol.x[-1], 10)
+    y = np.linspace(sol.y[0], sol.y[-1], 10)
+    z = np.array([0])
+
+    X, Y, Z = np.broadcast_arrays(x[:, None], y[None, :], z)
+
+    xi = 1
+
+    S = sol.result.S_tot_s(x[:, None], y[None, :])
+    U = scl * S[:, :, 0].real
+    V = scl * S[:, :, 1].real
+    W = scl * S[:, :, 2].real
+
+    R = np.column_stack((U.T.ravel(), V.T.ravel(), W.T.ravel()))
+
+    grid = pv.RectilinearGrid(x / xi, y / xi, z / xi)
+    grid["magnitude"] = np.linalg.norm(R, axis=1)
+    grid["vectors"] = R
+    grid.set_active_vectors("vectors")
+
+    bgrid = pv.RectilinearGrid(x / xi, y / xi, z / xi)
+
+    arrows = grid.glyph(orient="vectors", factor=1.6)
+
+    pv.set_plot_theme("document")
+
+    figscale = 6 if save else 1
+
+    p = pv.Plotter(
+        notebook=0, window_size=(500 * figscale, 300 * figscale), off_screen=save
+    )
+    p.add_mesh(arrows)
+    p.add_mesh(bgrid, opacity=0.25, color="gray", show_edges=False)
+
+    bounds = [
+        -sol.Lx / 2 / xi,
+        sol.Lx / 2 / xi,
+        -sol.Ly / 2 / xi,
+        sol.Ly / 2 / xi,
+        -0.7,
+        0.7,
+    ]
+
+    p.view_isometric()
+    p.remove_scalar_bar()
+    p.show_grid(
+        mesh=arrows,
+        show_xlabels=True,
+        show_ylabels=True,
+        show_zlabels=False,
+        show_zaxis=False,
+        grid=True,
+        xtitle="x",
+        ytitle="y",
+        ztitle="",
+        font_size=9 * figscale,
+        bounds=bounds,
+        n_xlabels=3,
+        n_ylabels=3,
+        padding=0.0,
+        fmt="%d",
+        location="outer",
+    )
+
+    if save:
+        img = np.array(p.show(return_img=True))
+        alpha = img.view("u1,u1,u1")[:, :, 0] == np.array((255, 255, 255), "u1,u1,u1")
+        alpha = np.uint8(255) - np.uint8(255) * alpha
+        img = np.concatenate((img, alpha[:, :, None]), axis=-1)
+        img = Image.fromarray(img)
+        img = img.crop(img.getbbox())
+        if save == True:
+            filename = "../figs/sfield.png"
+        else:
+            filename = save
+        img.save(filename)
+    else:
+        p.show()
+
+
+def main():
+    plot_fig2(save=True)
+    plot_fig5(save=True)
+    plot_fig4(save=True)
+
+
+if __name__ == "__main__":
+    main()

meson.build

@@ -78,5 +78,5 @@ if cpp_extra_args != []
   add_global_arguments(cpp_extra_args, language: 'cpp')
 endif
 
-subdir('usadelndsoc')
+subdir('src')
 subdir('doc')

pyproject.toml

@@ -2,8 +2,11 @@
 build-backend = 'mesonpy'
 requires = [
   'meson-python',
-  'oldest-supported-numpy',
+  'numpy < 2',
   'pybind11',
+  'Sphinx',
+  'numpydoc',
+  'scipy',
 ]
 
 [project]
@@ -24,7 +27,7 @@ classifiers = [
 requires-python = '>=3.7'
 dependencies = [
   'scipy >= 1.0.0',
-  'numpy >= 1.19.0',
+  'numpy >= 1.19.0, < 2',
 ]
 
 [project.optional-dependencies]
@@ -38,7 +41,7 @@ filterwarnings = [
   'ignore:the imp module is deprecated:DeprecationWarning',
   'ignore::DeprecationWarning:usadel1',
   'ignore::DeprecationWarning:numdifftools',
-  'ignore::usadel1.nonlinearsolver.BroydenWarning'
+  'ignore:.*Broyden reset',
 ]
 addopts = ['--doctest-modules', '--log-level=ERROR']
 testpaths = ['usadelndsoc', 'tests']

usadelndsoc/__init__.py → src/usadelndsoc/__init__.py

@@ -61,5 +61,14 @@ class with_log_level(contextlib.ContextDecorator):
         return False
 
 
+def docs():
+    """Open documentation in web browser"""
+    import webbrowser
+    import importlib.resources
+
+    p = importlib.resources.files("usadelndsoc") / "doc" / "html" / "index.html"
+    webbrowser.open_new("file://" + str(p.absolute()))
+
+
 _init_logging()
 del _init_logging, logging, contextlib

usadelndsoc/meson.build → src/usadelndsoc/meson.build

@@ -9,9 +9,9 @@ usadelndsoc_sources = files(
       'util.py',
 )
 
-usadelndsoc_pymod = py3.extension_module(
+usadelndsoc_pymod  = py3.extension_module(
   '_core',
-  ['../src/core.cpp'],
+  ['../core.cpp'],
   dependencies : deps + [numpy_dep, py3_dep, pybind11_dep],
   install : true,
   subdir : 'usadelndsoc',
@@ -19,6 +19,6 @@ usadelndsoc_pymod = py3.extension_module(
 
 py3.install_sources(
   usadelndsoc_sources,
-  pure : false,
   subdir : 'usadelndsoc',
+  pure : false,
 )

usadelndsoc/solver.py → src/usadelndsoc/solver.py

@@ -255,9 +255,9 @@ class Solver:
         Set system spatial size in cells.
         """
         self._core.set_shape(nx, ny)
-        self._Phi = np.zeros((nx, ny, 2, 2, 2), dtype=np.complex_)
-        self._J_tot = np.zeros((nx, ny, 4, 4, 4), dtype=np.complex_)
-        self._S_tot = np.zeros((nx, ny, 4, 4), dtype=np.complex_)
+        self._Phi = np.zeros((nx, ny, 2, 2, 2), dtype=np.complex128)
+        self._J_tot = np.zeros((nx, ny, 4, 4, 4), dtype=np.complex128)
+        self._S_tot = np.zeros((nx, ny, 4, 4), dtype=np.complex128)
         self._ac_tot = np.array(0)
 
     shape = _core_property("shape")
@@ -643,10 +643,10 @@ class Solver:
         return True
 
     def _eval_rhs(self, Phi):
-        return self._core.grad(Phi).ravel().view(np.float_)
+        return self._core.grad(Phi).ravel().view(np.float64)
 
     def _eval_jac_mul(self, Phi, dPhi):
-        return self._core.hess_mul(Phi, dPhi).ravel().view(np.float_)
+        return self._core.hess_mul(Phi, dPhi).ravel().view(np.float64)
 
     def _eval_jac(self, Phi):
         if not self._core.hess_computed:
@@ -662,7 +662,7 @@ class Solver:
         d2 = np.r_[H.data.real, -H.data.imag, -H.data.imag, -H.data.real]
         return coo_matrix(
             (d2, (i2, j2)),
-            dtype=np.float_,
+            dtype=np.float64,
             shape=(H.shape[0] * 2, H.shape[1] * 2),
             copy=False,
         ).tocsc()
@@ -1440,7 +1440,7 @@ def _self_consistency(
         constraint_x = np.array([])
 
     A_shape = (z.size, z.size)
-    A_dtype = np.float_
+    A_dtype = np.float64
 
     outer_v = []
     count = 0
@@ -2194,7 +2194,7 @@ def _interp_flow_1d(y):
     y = np.asarray(y)
     n = y.shape[0]
     shape = y.shape[1:]
-    dtype = np.result_type(y.dtype, np.float_)
+    dtype = np.result_type(y.dtype, np.float64)
 
     # KKT
     # variables: [x[0], ..., x[2*n], eta[0], ... eta[n-1]]