v0.7.0 core: Newton solitons and crystals, principal quadratures, thermal baths, marginal-mode guard

Author: Tanvir-Mahmud-Mahim

src/sqzcomb/__init__.py

@@ -22,17 +22,22 @@
                        logarithmic_negativity, ppt_symplectic_eigenvalue,
                        two_mode_reduction)
 from .gaussian import (covariance_xxpp, drift_from_qutip,
-                       intracavity_covariance, symplectic_eigenvalues)
+                       intracavity_covariance, principal_quadratures,
+                       symplectic_eigenvalues)
+from .soliton import continuation, newton_state, soliton_seed
+from .thermal import thermal_occupation
 
-__version__ = "0.6.0"
+__version__ = "0.7.0"
 __all__ = [
     "lle_evolve", "homogeneous_steady_states",
     "fluctuation_matrix", "single_mode_parametric",
     "photonic_molecule", "output_variance_ports", "molecule_threshold",
     "molecule_fluctuation_matrix",
+    "newton_state", "soliton_seed", "continuation",
     "output_quadrature_variance", "squeezing_db",
     "intracavity_covariance", "covariance_xxpp",
-    "symplectic_eigenvalues", "drift_from_qutip",
+    "symplectic_eigenvalues", "principal_quadratures",
+    "thermal_occupation", "drift_from_qutip",
     "two_mode_reduction", "ppt_symplectic_eigenvalue",
     "logarithmic_negativity", "duan_epr_sum", "entanglement_report",
     "detected_variance", "detected_squeezing_db", "dark_from_clearance_db",

src/sqzcomb/gaussian.py

@@ -19,8 +19,9 @@
 * Doubled basis z = (a_1..a_n, a*_1..a*_n); dz/dt = M z + inputs, with
   mode j carrying amplitude decay gamma_j (total input coupling
   sqrt(2 gamma_j), as everywhere in this package).
-* Vacuum input: <z_in(t) z_in(t')^dag> = N delta(t - t'), N = [[I, 0],
-  [0, 0]]; hence the diffusion matrix D = [[2 Gamma, 0], [0, 0]].
+* Bath input: <z_in(t) z_in(t')^dag> = N delta(t - t') with N =
+  [[(n_th + 1) I, 0], [0, n_th I]] (vacuum: n_th = 0); hence the
+  diffusion matrix D = [[2 Gamma (n_th + 1), 0], [0, 2 Gamma n_th]].
 * Steady covariance V = <z z^dag> solves M V + V M^dag + D = 0 (solved
   here with numpy alone via the Kronecker-vectorized linear system).
 * Quadratures x = (a + a^dag)/sqrt(2), p = -i (a - a^dag)/sqrt(2);
@@ -39,19 +40,27 @@
 import numpy as np
 
 
-def _diffusion(gammas):
+def _diffusion(gammas, n_th=0.0):
     gammas = np.asarray(gammas, dtype=float)
     n = gammas.size
+    nb = np.broadcast_to(np.asarray(n_th, dtype=float), (n,))
+    if np.any(nb < 0.0):
+        raise ValueError("thermal occupations must be non-negative")
     D = np.zeros((2 * n, 2 * n), dtype=complex)
-    D[:n, :n] = 2.0 * np.diag(gammas)
+    D[:n, :n] = 2.0 * np.diag(gammas * (nb + 1.0))
+    D[n:, n:] = 2.0 * np.diag(gammas * nb)
     return D
 
 
-def intracavity_covariance(M, gammas):
+def intracavity_covariance(M, gammas, n_th=0.0):
     """Steady-state complex covariance V = <z z^dag> of the doubled basis.
 
-    Solves M V + V M^dag + D = 0 with the vacuum-input diffusion matrix
-    D = diag(2 gamma, 0). Requires a stable M (all drift eigenvalues in
+    Solves M V + V M^dag + D = 0 with the bath diffusion matrix
+    D = diag(2 gamma (n_th + 1), 2 gamma n_th) -- vacuum baths for the
+    default n_th = 0, a Bose occupation per bath otherwise (scalar or
+    one value per mode; see `thermal_occupation` for the physical
+    number). A passive mode then holds exactly <a^dag a> = n_th, which
+    the tests assert. Requires a stable M (all drift eigenvalues in
     the open left half-plane); raises ValueError otherwise, for the same
     reason the spectra do: the linearized state does not exist above
     threshold.
@@ -67,7 +76,7 @@ def intracavity_covariance(M, gammas):
     gammas = np.asarray(gammas, dtype=float)
     if gammas.shape != (n,):
         raise ValueError("gammas must have one decay rate per mode")
-    D = _diffusion(gammas)
+    D = _diffusion(gammas, n_th)
     ident = np.eye(m2, dtype=complex)
     # row-major vec: vec(M V) = (M kron I) vec(V); vec(V M^dag) =
     # (I kron conj(M)) vec(V)
@@ -200,3 +209,45 @@ def fock(occ):
     B = -1j * G
     M = np.block([[A, B], [np.conj(B), np.conj(A)]])
     return M, gammas
+
+
+def principal_quadratures(sigma, hbar=2.0):
+    """Supermode decomposition of a multimode covariance matrix: the
+    principal quadratures and their variances.
+
+    For any real unit vector u, the generalized quadrature u . r
+    (r = (x_1..x_n, p_1..p_n)) has variance u^T sigma u, so the
+    eigendecomposition of sigma answers, exactly and completely, the
+    question "what is the most squeezed collective quadrature this
+    state contains, and along which mode combination does it lie" --
+    the smallest eigenvalue is the deepest squeezing any generalized
+    quadrature attains, its eigenvector the supermode that carries it.
+    That statement is linear algebra, not approximation, and the tests
+    pin it: the two-mode squeezed vacuum yields the exact pairs
+    (hbar/2) e^{-2r} and (hbar/2) e^{+2r} with the EPR combinations
+    (x_1 -/+ x_2)/sqrt(2), (p_1 +/- p_2)/sqrt(2) as supermodes, vacuum
+    yields hbar/2 in every direction, and on the photonic molecule the
+    principal variance is verified to lower-bound every tested
+    twin-beam quadrature.
+
+    Returns (variances, vectors): eigenvalues ascending, vectors[:, i]
+    the unit xxpp vector of principal quadrature i.
+
+    This is the orthogonal decomposition of the noise ellipsoid --
+    deliberately distinct from the *symplectic* (Williamson)
+    decomposition `symplectic_eigenvalues`, which measures mixedness:
+    a pure squeezed state has all symplectic eigenvalues at hbar/2
+    while its principal variances split as e^{-/+ 2r}. Both views are
+    exported because they answer different questions.
+    """
+    sigma = np.asarray(sigma, dtype=float)
+    m2 = sigma.shape[0]
+    if sigma.shape != (m2, m2) or m2 % 2:
+        raise ValueError("sigma must be a square (2n, 2n) xxpp matrix")
+    if np.abs(sigma - sigma.T).max() > 1e-9 * max(1.0, np.abs(sigma).max()):
+        raise ValueError("sigma must be symmetric")
+    w, v = np.linalg.eigh(0.5 * (sigma + sigma.T))
+    if w[0] < -1e-12 * max(1.0, abs(w[-1])):
+        raise ValueError("sigma is not positive semidefinite; not a "
+                         "covariance matrix")
+    return w, v

src/sqzcomb/soliton.py

@@ -0,0 +1,173 @@
+"""Localized steady states of the LLE: solitons, soliton crystals,
+and parameter continuation.
+
+`newton_state` solves the stationary Lugiato-Lefever equation
+
+    0 = -(1 + i alpha) psi + i |psi|^2 psi + i D psi + F
+
+by Newton's method in Fourier space -- and its Jacobian is the same
+linearization the package's `fluctuation_matrix` builds around a
+steady state, evaluated at the current iterate: the operator that
+damps a perturbation of a steady state and the operator Newton inverts
+to find that steady state are one and the same.  (The solver keeps the
+grid's wrapped convolution entries that the physically mode-truncated
+`fluctuation_matrix` sets to zero -- the exact Jacobian of the
+*discrete* equations, needed for machine-precision convergence.  The
+wrapped entries only transfer amplitude across the band edge, so on
+the resolved soliton's own modes the two operators act identically:
+the test suite asserts that both annihilate the translation mode to
+1e-9.)
+
+Exact facts the test suite asserts, rather than states:
+
+* the converged residual is below the requested tolerance (the solver
+  refuses to return otherwise), and the converged state, handed to the
+  independent split-step evolver `lle_evolve`, stays put;
+* from a flat seed, Newton lands on the exact root of the homogeneous
+  cubic;
+* an N-fold soliton crystal on the 2 pi ring equals, to machine
+  precision, the single soliton of the same equation with dispersion
+  d_m scaled by N^m, resampled onto the crystal grid -- the exact
+  rescaling theta -> N theta of the periodic domain;
+* the linearization around any localized steady state carries an
+  *exact* zero eigenvalue whose eigenvector is the translation mode
+  d psi / d theta (Goldstone mode of the broken translation symmetry);
+  the drift matrix is therefore marginally, not asymptotically, stable,
+  and the spectra machinery refuses it unless told that the marginal
+  direction is understood (`allow_marginal`, see `spectra`).
+
+The sech seed uses the standard soliton asymptotics of the anomalous
+LLE (amplitude ~ sqrt(2 alpha), width ~ sqrt(|d2| / alpha), pump phase
+cos phi_0 = sqrt(8 alpha) / (pi F); T. Herr et al., Nature Photonics 8,
+145 (2014)); the *converged* answers never depend on the seed quality,
+only Newton's basin does.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+from .lle import _linear_symbol, homogeneous_steady_states
+
+__all__ = ["newton_state", "soliton_seed", "continuation"]
+
+
+def _residual_k(psi, F, alpha, dispersion):
+    """Stationary-LLE residual in Fourier space (FFT convention)."""
+    n = psi.size
+    k = np.fft.fftfreq(n, d=1.0 / n)
+    L = _linear_symbol(alpha, dispersion, k)
+    Rk = L * np.fft.fft(psi) + np.fft.fft(1j * np.abs(psi) ** 2 * psi)
+    Rk[0] += F * n
+    return Rk
+
+
+def _jacobian(psi, alpha, dispersion):
+    """Exact Jacobian of the discrete (pseudospectral) residual, in the
+    doubled Fourier basis.
+
+    Identical to `fluctuation_matrix` except that the products
+    2 |psi|^2 delta and psi^2 conj(delta) are the *grid's* circular
+    convolutions (mode transfers taken modulo n): the exact Jacobian
+    of the discrete pseudospectral residual, which Newton needs to
+    converge to machine precision.  The entries where the two differ
+    are exactly the wrapped transfers across the band edge; on
+    perturbations resolved inside the band the two operators act
+    identically, which the test suite asserts on the soliton's
+    translation mode.
+    """
+    psi = np.asarray(psi, dtype=complex)
+    n = psi.size
+    kgrid = np.fft.fftfreq(n, d=1.0 / n).astype(int)
+    f_abs2 = np.fft.fft(np.abs(psi) ** 2) / n
+    f_sq = np.fft.fft(psi ** 2) / n
+    L = _linear_symbol(alpha, dispersion, kgrid.astype(float))
+    # FFT layout: the component of mode value v sits at slot v mod n
+    A = 2j * f_abs2[(kgrid[:, None] - kgrid[None, :]) % n]
+    B = 1j * f_sq[(kgrid[:, None] + kgrid[None, :]) % n]
+    A[np.arange(n), np.arange(n)] += L
+    return np.block([[A, B], [np.conj(B), np.conj(A)]])
+
+
+def newton_state(psi0, F, alpha, dispersion=(0.0,), tol=1e-12,
+                 maxiter=60):
+    """Newton solution of the stationary LLE from the seed psi0.
+
+    Returns (psi, info) with info = {"residual", "iterations"}; the
+    residual is the max-norm of the stationary equation on the theta
+    grid.  Raises RuntimeError if the tolerance is not reached -- a
+    state is never returned unverified.  The Jacobian at each iterate
+    is the exact discrete form of `fluctuation_matrix(psi, alpha,
+    dispersion)` (see `_jacobian`).
+    """
+    psi = np.asarray(psi0, dtype=complex).copy()
+    n = psi.size
+    res = np.inf
+    for it in range(int(maxiter)):
+        Rk = _residual_k(psi, F, alpha, dispersion)
+        res = float(np.abs(np.fft.ifft(Rk)).max())
+        if res < tol:
+            return psi, {"residual": res, "iterations": it}
+        J = _jacobian(psi, alpha, dispersion)
+        rhs = np.concatenate([Rk / n, np.conj(Rk) / n])
+        # minimal-norm least-squares step: at a localized state the
+        # Jacobian is exactly singular along the translation (Goldstone)
+        # mode, and the minimal-norm solution is the Newton step that
+        # does not slide along the soliton position
+        delta = np.linalg.lstsq(J, -rhs, rcond=1e-9)[0]
+        psi = psi + np.fft.ifft(delta[:n] * n)
+    raise RuntimeError(
+        f"Newton did not reach tol = {tol:g} in {maxiter} iterations "
+        f"(residual {res:.3e}); refusing to return an unconverged "
+        "state")
+
+
+def soliton_seed(n, F, alpha, dispersion, n_pulses=1):
+    """Sech ansatz seed for `newton_state`: the lower flat state plus
+    ``n_pulses`` equally spaced sech pulses with the standard
+    asymptotic amplitude sqrt(2 alpha), width sqrt(|d2| / (2 alpha))
+    and pump phase cos phi_0 = sqrt(8 alpha) / (pi F) (Herr et al.,
+    Nature Photonics 8, 145 (2014)).  Anomalous dispersion (d2 < 0 in
+    this package's sign convention) required.
+    """
+    d2 = float(dispersion[0])
+    if d2 >= 0.0:
+        raise ValueError("bright-soliton seeds need anomalous "
+                         "dispersion (d2 < 0 in this convention)")
+    roots = homogeneous_steady_states(F, alpha)
+    rho_low = float(roots[0])
+    psi_low = F / (1.0 + 1j * (alpha - rho_low))
+    arg = np.sqrt(8.0 * alpha) / (np.pi * F)
+    phi0 = np.arccos(min(arg, 1.0))
+    theta = 2.0 * np.pi * np.arange(n) / n
+    width = np.sqrt(abs(d2) / (2.0 * alpha))
+    psi = np.full(n, psi_low, dtype=complex)
+    for p in range(int(n_pulses)):
+        center = 2.0 * np.pi * (p + 0.5) / n_pulses
+        d = np.angle(np.exp(1j * (theta - center)))   # wrapped distance
+        psi += np.sqrt(2.0 * alpha) * np.exp(1j * phi0) / np.cosh(d / width)
+    return psi
+
+
+def continuation(psi0, F, alphas, dispersion=(0.0,), tol=1e-12,
+                 maxiter=60):
+    """Sweep the detuning: Newton-solve at each alpha in ``alphas``,
+    seeding each step with the previous converged state.
+
+    Returns (states, infos) as lists.  Every step is verified to the
+    tolerance or the sweep raises (naming the failing alpha), so a
+    returned branch contains no unconverged member.  Sweeping F at
+    fixed alpha works by transposing the roles: call in a loop, this
+    function stays deliberately simple.
+    """
+    psi = np.asarray(psi0, dtype=complex).copy()
+    states, infos = [], []
+    for a in alphas:
+        try:
+            psi, info = newton_state(psi, F, float(a), dispersion,
+                                     tol=tol, maxiter=maxiter)
+        except RuntimeError as exc:
+            raise RuntimeError(
+                f"continuation failed at alpha = {a}: {exc}") from exc
+        states.append(psi.copy())
+        infos.append(info)
+    return states, infos

src/sqzcomb/spectra.py

@@ -17,38 +17,72 @@
 
 import numpy as np
 
-from .linearize import is_stable
 
+def _bath_covariance(m2, n_bar):
+    """<z_in z_in^dagger> of a thermal bath in the doubled ordering:
+    <a a^dagger> = n_bar + 1, <a^dagger a> = n_bar (vacuum: n_bar = 0)."""
+    half = m2 // 2
+    N = np.zeros((m2, m2), dtype=complex)
+    N[:half, :half] = (float(n_bar) + 1.0) * np.eye(half)
+    N[half:, half:] = float(n_bar) * np.eye(half)
+    return N
+
+
+def _check_spectra_stability(M, allow_marginal, tol=1e-6):
+    lam = float(np.max(np.linalg.eigvals(np.asarray(M)).real))
+    if lam < -tol:
+        return
+    if allow_marginal and lam < tol:
+        return
+    if lam < tol:
+        raise ValueError(
+            "drift matrix is marginally stable (an eigenvalue's real "
+            "part is numerically zero). Around a localized steady "
+            "state this is the exact translation (Goldstone) mode of "
+            "the soliton; spectra at omega != 0 remain finite, so pass "
+            "allow_marginal=True if that marginal direction is "
+            "understood. A genuinely positive growth rate is still "
+            "refused.")
+    raise ValueError("drift matrix is unstable (above threshold); "
+                     "linearized spectra are meaningless there")
 
-def _output_covariance(M, eta, omega):
+
+def _output_covariance(M, eta, omega, n_th_port=0.0, n_th_loss=0.0):
     m2 = M.shape[0]
     ident = np.eye(m2, dtype=complex)
     G = np.linalg.inv(-1j * omega * ident - M)
     T_ex = 2.0 * eta * G - ident
     T_0 = 2.0 * np.sqrt(eta * (1.0 - eta)) * G
-    # vacuum covariance <z_in z_in^dagger> in the (a, a*) doubled ordering:
-    # <a a^dagger> = 1, all other blocks zero
-    half = m2 // 2
-    N = np.zeros((m2, m2), dtype=complex)
-    N[:half, :half] = np.eye(half)
-    S = T_ex @ N @ T_ex.conj().T + T_0 @ N @ T_0.conj().T
+    S = T_ex @ _bath_covariance(m2, n_th_port) @ T_ex.conj().T \
+        + T_0 @ _bath_covariance(m2, n_th_loss) @ T_0.conj().T
     return S
 
 
 def output_quadrature_variance(M, eta, omega, mode_index, n_modes,
-                               phi=0.0, mode_index_b=None):
+                               phi=0.0, mode_index_b=None,
+                               n_th_port=0.0, n_th_loss=0.0,
+                               allow_marginal=False):
     """Symmetrized variance of an output quadrature at frequency omega.
 
     mode_index : index of the mode (within the retained mode list) whose
         quadrature is detected. If mode_index_b is given, the joint
         two-mode quadrature (a + b)/sqrt(2) rotated by phi is used, the
         natural variable for twin-beam squeezing.
+    n_th_port, n_th_loss : Bose occupations of the extraction-port
+        input and of the intrinsic-loss bath (default vacuum, 0; see
+        `thermal_occupation` for the physical number). A passive cavity
+        with both baths at n_bar emits exactly (2 n_bar + 1)/2 at every
+        frequency, coupling and phase -- asserted in the tests.
+    allow_marginal : accept a drift matrix whose largest eigenvalue
+        real part is numerically zero (the soliton's exact translation
+        Goldstone mode); spectra at omega != 0 stay finite. Genuinely
+        unstable matrices are refused regardless.
     Vacuum level is 0.5. Requires a stable M.
     """
-    if not is_stable(M):
-        raise ValueError("drift matrix is unstable (above threshold); "
-                         "linearized spectra are meaningless there")
-    S = _output_covariance(M, eta, omega)
+    if float(n_th_port) < 0.0 or float(n_th_loss) < 0.0:
+        raise ValueError("thermal occupations must be non-negative")
+    _check_spectra_stability(M, allow_marginal)
+    S = _output_covariance(M, eta, omega, n_th_port, n_th_loss)
     u = np.zeros(2 * n_modes, dtype=complex)
     if mode_index_b is None:
         u[mode_index] = np.exp(-1j * phi) / np.sqrt(2.0)