dIdV3

dIdV3

Calculating differential conductence (dI/dV) spectra in 3rd order.

This module provides routines to compute differential conductance spectra using third-order perturbation theory.

Functions:

Name	Description
calc_IETS_3rd	Calculate conductance in 3rd order.
calc_IETS_3rdR	Calculate third-order spectra from corrected rates and populations.

Dependencies

• Uses spinfinity.SpinSys.SpinSys class to represent the spin system and its properties.
• Uses spinfinity.SpinHamiltonian for eigenvalue and eigenvector calculations.
• Uses functions from spinfinity.Bath for population and rate calculations.

References

• Ternes, M., Spin excitations and correlations in scanning tunneling spectroscopy. New J. Phys. 17, 063016 (2015). https://dx.doi.org/10.1088/1367-2630/17/6/063016

calc_IETS_3rd(spinsys, V_array=None, norm=False, energy_step=None)

Calculate conductance in 3rd order.

Parameters:

Name	Type	Description	Default
spinsys	SpinSys	Spin system. J contains signed dimensionless substrate couplings per spin; J_tip is the signed dimensionless tip coupling of the readout spin. T is in K; w and lt are in meV. The Hamiltonian is rebuilt.	required
V_array	array_like	Bias voltages in mV. Default is 200 points from -30 to 30 mV.	None
norm	bool	Divide by the first conductance value if True. Default is False.	False
energy_step	float	Logarithm convolution grid spacing in meV. Default is one quarter of the smaller nonzero scale among lt and kB*T.	None

Returns:

Name	Type	Description
V_array	ndarray	Bias voltages in mV.
conductance	ndarray	Differential conductance in S, or dimensionless when norm=True.

Raises:

Type	Description
ValueError	If the ground-state elastic normalization or the requested first-point normalization is zero.

Notes

Uses thermal populations and factorized broadened-step times logarithm contributions, as in Ternes (2015), Eqs. (24)-(25). J=J_tip=0 exactly recovers calc_IETS_2nd. Eigenstates and the existing second-order thermal caches are refreshed.

Source code in spinfinity\dIdV3.py

def calc_IETS_3rd(spinsys: SpinSys, V_array=None, norm=False, energy_step=None):
    """Calculate conductance in 3rd order.

    Parameters
    ----------
    spinsys : SpinSys
        Spin system. J contains signed dimensionless substrate couplings
        per spin; J_tip is the signed dimensionless tip coupling of the
        readout spin. T is in K; w and lt are in meV.
        The Hamiltonian is rebuilt.
    V_array : array_like, optional
        Bias voltages in mV. Default is 200 points from -30 to 30 mV.
    norm : bool, optional
        Divide by the first conductance value if True. Default is False.
    energy_step : float, optional
        Logarithm convolution grid spacing in meV. Default is one quarter
        of the smaller nonzero scale among lt and kB*T.

    Returns
    -------
    V_array : ndarray
        Bias voltages in mV.
    conductance : ndarray
        Differential conductance in S, or dimensionless when norm=True.

    Raises
    ------
    ValueError
        If the ground-state elastic normalization or the requested
        first-point normalization is zero.

    Notes
    -----
    Uses thermal populations and factorized broadened-step times logarithm
    contributions, as in Ternes (2015), Eqs. (24)-(25).
    J=J_tip=0 exactly recovers calc_IETS_2nd.
    Eigenstates and the existing second-order thermal caches are refreshed.
    """
    V, thermal_energy, _, grid, logarithm, paths, tip_paths = _prepare(
        spinsys, V_array, energy_step
    )

    rate_0 = bath.calc_TunnelingMatrixElements(spinsys, "st")[0, 0, 0]
    if rate_0 <= 0:
        raise ValueError("The ground-state elastic normalization is zero.")

    _, conductance = dIdV.calc_IETS_2nd(spinsys, V_array=V, norm=False)
    derivatives = _factorized_conductance(
        spinsys.E_All,
        V,
        paths,
        thermal_energy,
        grid,
        logarithm,
    )

    # The tip-intermediate calculation uses reversed reservoir roles.
    derivatives += _factorized_conductance(
        spinsys.E_All,
        -V,
        tip_paths,
        thermal_energy,
        grid,
        logarithm,
    )[::-1]

    thermal = bath.calc_ThermalPopulations(spinsys)
    correction = np.einsum("i,cvij->v", thermal, derivatives)
    conductance += (1 - spinsys.b0) * spinsys.G_st / rate_0 * correction

    if norm:
        if conductance[0] == 0:
            raise ValueError("Cannot normalize by a zero first conductance value.")
        conductance = conductance / conductance[0]

    return V, conductance

calc_IETS_3rdR(spinsys, V_array=None, AllowPumping=True, norm=False, energy_step=None)

Calculate third-order spectra from corrected rates and populations.

Parameters:

Name	Type	Description	Default
spinsys	SpinSys	Spin system. J contains signed dimensionless substrate couplings per spin; J_tip is the signed dimensionless tip coupling of the readout spin. T is in K; w and lt are in meV. The Hamiltonian is rebuilt.	required
V_array	array_like	Strictly increasing bias voltages in mV, at least two points. Default is 200 points from -30 to 30 mV.	None
AllowPumping	bool	Solve the corrected stationary rate equation if True (default); otherwise use thermal populations at every bias.	True
norm	bool	Divide by the first conductance value if True. Default is False.	False
energy_step	float	Internal energy spacing in meV; default as in calc_IETS_3rd.	None

Returns:

Name	Type	Description
V_array	ndarray	Bias voltages in mV.
conductance	ndarray	Differential conductance in S, or dimensionless when norm=True.
populations	ndarray	Eigenstate populations, shape (dimensionOfMatrix, len(V_array)).

Raises:

Type	Description
ValueError	If the voltage grid is not strictly increasing with at least two points, a normalization is zero, enabled G_tt has nonpositive readout G_ss, total rates become negative, or the steady state is not unique.

Source code in spinfinity\dIdV3.py

def calc_IETS_3rdR(
    spinsys: SpinSys,
    V_array=None,
    AllowPumping=True,
    norm=False,
    energy_step=None,
):
    """Calculate third-order spectra from corrected rates and populations.

    Parameters
    ----------
    spinsys : SpinSys
        Spin system. J contains signed dimensionless substrate couplings
        per spin; J_tip is the signed dimensionless tip coupling of the
        readout spin. T is in K; w and lt are in meV.
        The Hamiltonian is rebuilt.
    V_array : array_like, optional
        Strictly increasing bias voltages in mV, at least two points.
        Default is 200 points from -30 to 30 mV.
    AllowPumping : bool, optional
        Solve the corrected stationary rate equation if True (default);
        otherwise use thermal populations at every bias.
    norm : bool, optional
        Divide by the first conductance value if True. Default is False.
    energy_step : float, optional
        Internal energy spacing in meV; default as in calc_IETS_3rd.

    Returns
    -------
    V_array : ndarray
        Bias voltages in mV.
    conductance : ndarray
        Differential conductance in S, or dimensionless when norm=True.
    populations : ndarray
        Eigenstate populations, shape (dimensionOfMatrix, len(V_array)).

    Raises
    ------
    ValueError
        If the voltage grid is not strictly increasing with at least two
        points, a normalization is zero, enabled G_tt has nonpositive
        readout G_ss, total rates become negative, or the steady state is not unique.
    """
    V, thermal_energy, step, grid, logarithm, paths, tip_paths = _prepare(
        spinsys, V_array, energy_step
    )
    if len(V) < 2 or np.any(np.diff(V) <= 0):
        raise ValueError(
            "V_array must be strictly increasing with at least two points."
        )

    energies = spinsys.E_All
    dimension = len(energies)
    transport = _phase_integrals(
        energies, V, paths, thermal_energy, grid, logarithm, step
    )
    transport += _phase_integrals(
        energies, -V, tip_paths, thermal_energy, grid, logarithm, step
    )[::-1]

    matrix_st = bath.calc_TunnelingMatrixElements(spinsys, "st")[0]
    matrix_ss = bath.calc_TunnelingMatrixElements(spinsys, "ss")
    matrix_tt = bath.calc_TunnelingMatrixElements(spinsys, "tt")[0]

    rate_0 = matrix_st[0, 0]
    if rate_0 <= 0:
        raise ValueError("The ground-state elastic normalization is zero.")

    factor = 1e-3 / (const.e * rate_0)
    effective = (1 - spinsys.b0) * spinsys.G_st

    # Local substrate corrections, summed incoherently over scattering sites.
    substrate = np.zeros((dimension, dimension))
    for site in range(spinsys.NSpins):
        local_paths = _pathways(
            spinsys,
            site,
            spinsys.SamplePolarization,
            spinsys.SamplePolarization,
            spinsys.J[site],
            potential=False,
        )
        correction = _phase_integrals(
            energies,
            np.array([0.0]),
            local_paths,
            thermal_energy,
            grid,
            logarithm,
            step,
        )[0, 0]
        substrate += spinsys.G_ss[site] * correction

    tip_paths = _pathways(
        spinsys,
        spinsys.ReadoutSpin,
        spinsys.TipPolarization,
        spinsys.TipPolarization,
        spinsys.J_tip,
        potential=False,
    )
    tip_correction = _phase_integrals(
        energies,
        np.array([0.0]),
        tip_paths,
        thermal_energy,
        grid,
        logarithm,
        step,
    )[0, 0]

    tip_scale = 0.0
    if spinsys.G_tt:
        if spinsys.G_ss[spinsys.ReadoutSpin] <= 0:
            raise ValueError("G_tt requires positive G_ss for the readout spin.")
        tip_scale = effective**2 / spinsys.G_ss[spinsys.ReadoutSpin] * spinsys.G_tt

    old_voltage = spinsys.V_DC
    populations = np.zeros((dimension, len(V)))
    current = np.zeros(len(V))
    thermal = bath.calc_ThermalPopulations(spinsys)

    try:
        bath_integral = bath.calc_RateIntegrals(spinsys, "ss")
        fixed_ss = factor * (
            np.einsum("k,kij->ij", spinsys.G_ss, matrix_ss) * bath_integral + substrate
        )
        fixed_tt = factor * tip_scale * (matrix_tt * bath_integral + tip_correction)

        for index, voltage in enumerate(V):
            spinsys.V_DC = voltage
            rates = np.array(
                [
                    factor
                    * effective
                    * (
                        matrix_st.T * bath.calc_RateIntegrals(spinsys, "ts")
                        + transport[0, index]
                    ),
                    factor
                    * effective
                    * (
                        matrix_st * bath.calc_RateIntegrals(spinsys, "st")
                        + transport[1, index]
                    ),
                    fixed_ss,
                    fixed_tt,
                ]
            )

            if np.any(rates < 0):
                raise ValueError(
                    "Negative total rates: reduce |J|/|J_tip| or increase T."
                )

            if AllowPumping:
                total = np.sum(rates, axis=0)
                np.fill_diagonal(total, 0)
                generator = total.T.copy()
                np.fill_diagonal(generator, -np.sum(total, axis=1))
                steady = null_space(generator)

                if steady.shape[1] != 1:
                    raise ValueError("The stationary population is not unique.")

                population = steady[:, 0] / np.sum(steady[:, 0])
            else:
                population = thermal

            populations[:, index] = population
            current[index] = (
                const.e * np.sum(population[:, None] * (rates[0] - rates[1]))
                + spinsys.b0 * spinsys.G_st * voltage * 1e-3
            )
    finally:
        spinsys.V_DC = old_voltage

    conductance = np.gradient(current, V * 1e-3)
    if norm:
        if conductance[0] == 0:
            raise ValueError("Cannot normalize by a zero first conductance value.")
        conductance = conductance / conductance[0]

    return V, conductance, populations