diff --git a/tests/test_solver.py b/tests/test_solver.py
index 05688129faa7cd2ed3b3de9445589f5425281772..52d00dc61feef28aa0a73844a80d30a988321e85 100644
--- a/tests/test_solver.py
+++ b/tests/test_solver.py
@@ -61,7 +61,7 @@ def test_example_sns_1d_J_usadel1():
sol.Ly = 1
T = 0.2
- Delta, I0, _, success = sol.self_consistency(T=T, T_c0=0, tol=1e-10)
+ Delta, I0, _, _, success = sol.self_consistency(T=T, T_c0=0, tol=1e-10)
tau3 = np.diag([1, 1, -1, -1])
I0 = tr(I0[:-1, 0, 0] @ tau3)
@@ -101,7 +101,7 @@ def test_example_sss_1d_J_usadel1():
T_c0 = np.exp(np.euler_gamma) / np.pi
T = 0.2
- Delta, I0, _, success = sol.self_consistency(T=T, T_c0=T_c0, workers=1)
+ Delta, I0, _, _, success = sol.self_consistency(T=T, T_c0=T_c0, workers=1)
tau3 = np.diag([1, 1, -1, -1])
I0 = tr(I0[:-1, 0, 0] @ tau3)
@@ -361,7 +361,7 @@ def test_selfcons_h():
solver.Omega[...] += h * np.kron(S_z, S_y)
T_c0 = np.exp(np.euler_gamma) / np.pi
- Delta, I0, _, success = solver.self_consistency(T=T, T_c0=T_c0, workers=1)
+ Delta, I0, _, _, success = solver.self_consistency(T=T, T_c0=T_c0, workers=1)
Delta0 = BCS_Delta(T, T_c0, h=h)
Deltam = np.broadcast_to((Delta0 * S_0)[None, None, :, :], Delta.shape)
diff --git a/usadelndsoc/plotting.py b/usadelndsoc/plotting.py
index 4708ad5cea2f76cfdc68b938a359ed2d8e8a5f79..a09e8934dc8392c6f0dd2798711e8461d880a3bf 100644
--- a/usadelndsoc/plotting.py
+++ b/usadelndsoc/plotting.py
@@ -26,8 +26,8 @@ def plot_Delta_J_2d(x, y, Delta, xc, yc, Ix, Iy, ax=None):
m = pcolormesh_complex(x, y, Delta, ax=ax)
lw = np.hypot(Ix, Iy)
- lw = 6 * (lw / max(1.0, lw.max())) ** 0.5
- ax.streamplot(xc, yc, Ix.T, Iy.T, color=(0.5, 0.5, 0.5), linewidth=lw.T)
+ lw = 4 * (lw / max(1e-6, lw.max())) ** 1.0
+ ax.streamplot(xc, yc, Ix.T, Iy.T, color=(0.5, 0.5, 0.5), linewidth=lw.T, density=2)
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")
diff --git a/usadelndsoc/solver.py b/usadelndsoc/solver.py
index 6e9af73f4adc0b06117dc710fa7e684ed86546e7..f8ca00968b33b2a73e389a0af1e9a9518cad7228 100644
--- a/usadelndsoc/solver.py
+++ b/usadelndsoc/solver.py
@@ -152,6 +152,7 @@ class Solver:
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_)
shape = _core_property("shape")
Lx = _core_property("Lx")
@@ -170,11 +171,20 @@ class Solver:
Phi = _array_property(
"_Phi",
- doc="""Current saddle-point solution in Riccati parametrization. Array of complex, shape (nx, ny, 2, 2, 2).""",
+ doc=r"""Current saddle-point solution in Riccati parametrization. Array of complex, shape (nx, ny, 2, 2, 2).""",
)
J_tot = _array_property(
"_J_tot",
- doc="""Matsubara-summed matrix current between cells. Array of complex, shape (nx, ny, 4(direction), 4, 4)""",
+ doc=r"""
+ Matsubara-summed matrix current between cells, :math:`J = -i \pi T\sum_{\omega_n} \mathcal{J}(\omega_n)`.
+ Array of complex, shape (nx, ny, 4(direction), 4, 4).
+ """,
+ )
+ S_tot = _array_property(
+ "_S_tot",
+ doc=r"""Matsubara-summed density, :math:`S = \pi T \sum_{\omega_n} [Q(\omega_n) - \tau_3]`.
+ Array of complex, shape (nx, ny, 4, 4).
+ """,
)
def reset(self):
@@ -661,10 +671,10 @@ class Solver:
return Result(self, omega=omega0, Phi=Phi)
def self_consistency(self, T, T_c0, continuation=False, **kw):
- """
+ r"""
Solve self-consistency equations for Delta.
- Assumes simple singlet weak coupling.
+ Assumes s-wave singlet weak coupling.
Supports additional constraint variables and equations.
@@ -695,6 +705,8 @@ class Solver:
Order parameter matrix
J_tot
Matsubara-summed matrix current
+ S_tot
+ Matsubara-summed matrix :math:`Q-\tau_3`
cx
Final constraint variable values
succes : bool
@@ -708,8 +720,8 @@ class Solver:
mask = self.mask == MASK_NONE
self.Delta[mask] = 2 * T_c0[mask, None, None] * np.eye(2)
res = _self_consistency(self, T, T_c0, **kw)
- self.Delta[...], self.J_tot[...], cx, success = res
- return self.Delta, self.J_tot, cx, success
+ self.Delta[...], self.J_tot[...], self.S_tot[...], cx, success = res
+ return self.Delta, self.J_tot, self.S_tot, cx, success
def _bcast_left(x, shape):
@@ -791,6 +803,18 @@ class Result:
with np.errstate(invalid="ignore"):
return 2 * s * np.linalg.solve(I + gt @ g, gt)
+ @property
+ def Q(self):
+ """
+ Full Green function / Q-matrix. Shape ([nw,] nx, ny, 4, 4).
+ """
+ r = np.zeros(self.shape + (4, 4), dtype=complex)
+ r[..., :2, :2] = self.G
+ r[..., 2:, 2:] = self.Gc
+ r[..., :2, 2:] = self.F
+ r[..., 2:, :2] = self.Fc
+ return r
+
def _get_J(self, omega, Phi):
old_omega = self._core.omega
self._core.omega = omega
@@ -1002,6 +1026,7 @@ def _self_consistency(
success = True
J = None
+ S = None
solver_kw["tol"] = 0.1 * tol
@@ -1011,7 +1036,7 @@ def _self_consistency(
# Solve linearized problem
def ev(z):
- nonlocal J
+ nonlocal J, S
for k, v in cache:
if (z == k).all():
@@ -1026,7 +1051,7 @@ def _self_consistency(
# Constraint function updates terminal phases
solver.Delta[...] = singlet_m(Delta1)
- res, J, m = _self_consistent_Delta_f(
+ res, J, S, m = _self_consistent_Delta_f(
solver, Delta1, T, T_c0, workers=workers, **solver_kw
)
res = res[m].ravel().view(float)
@@ -1129,7 +1154,7 @@ def _self_consistency(
if success:
_log.info(f"Self-consistent iteration converged with {count} ops")
- return Delta, J, cx, success
+ return Delta, J, S, cx, success
@_with_workers
@@ -1148,6 +1173,7 @@ def _self_consistent_Delta_f(
rtot = np.zeros(tuple(solver.shape), dtype=complex)
Jtot = np.zeros(tuple(solver.shape) + (4, 4, 4), dtype=complex)
+ Stot = np.zeros(tuple(solver.shape) + (4, 4), dtype=complex)
mask1 = solver.mask == MASK_NONE
mask2 = (T_c0 > 0) & (T < T_c0)
@@ -1170,13 +1196,14 @@ def _self_consistent_Delta_f(
else:
work = lambda: (_mp_one((w, a, solver, solver_kw)),)
- for rtotx, Jtotx in work():
+ for rtotx, Jtotx, Stotx in work():
rtot[mask] += rtotx[mask]
Jtot += Jtotx
+ Stot += Stotx
rtot[mask] -= np.log(T_c0[mask] / T) * singlet(solver.Delta[mask])
- return rtot, Jtot, mask
+ return rtot, Jtot, Stot, mask
_prev_solver = None
@@ -1198,22 +1225,26 @@ def _mp_one(args, solver=None, solver_kw=None):
else:
_prev_solver = solver
+ tau_3 = np.diag([1, 1, -1, -1])
+
old_level = _log_solve.level
try:
_log_solve.setLevel(_log.getEffectiveLevel() + 10)
rtot = 0
Jtot = 0
+ Stot = 0
for wx, ax in zip(w, a):
res = solver.solve(omega=wx, **solver_kw)
r = singlet(solver.Delta) / abs(wx) - singlet(res.F)
rtot += (2 * np.pi * ax) * r
Jtot += (-2j * np.pi * ax) * res.J
+ Stot += (2 * np.pi * ax) * (res.Q - tau_3)
finally:
_log_solve.setLevel(old_level)
- return rtot, Jtot
+ return rtot, Jtot, Stot
class CPRState:
@@ -1296,7 +1327,11 @@ def cpr(
"""
Compute current-phase relation with selfconsistency.
- Adjust superconducting phases in cells indicated by `phase_mask`.
+ Adjust superconducting phases in cells indicated by `phase_mask`,
+ and trace the current-phase relation using pseudoarclength continuation.
+ This recovers also the unstable branches of the CPR.
+
+ Assumes s-wave singlet weak coupling.
Parameters
----------
@@ -1310,13 +1345,23 @@ def cpr(
Cells where to tune superconducting phase
max_points : int
ds : float
- Step size
+ Arclength step size
phi0 : float
filename : str, optional
File to resume or save progress. Computation is resumed if
the file exists and parameter values in it match.
**selfcons_kw
Parameters passed on to `Solver.self_consistency`.
+
+ Returns
+ -------
+ phis : array, shape (nphi,)
+ Phase differences calculation was done at.
+ Deltas : array, shape (nphi, nx, ny, 2, 2)
+ Self-consistent order parameter matrix, corresponding to each phase difference.
+ Js : array, shape (nphi, nx, ny, 4(direction), 4, 4)
+ Matsubara summed matrix current, corresponding to each phase difference.
+
"""
nx, ny = solver.shape