landauzener2d.cpp

Go to the documentation of this file.

/**
 * @file landauzener2d.cpp
 * @brief Implements the landauzener2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/landauzener2d.cpp. This file is part of Quantum Phenomena Explorer (QPE),
 * maintained by Officina Turini. See LICENSE_IT.txt / LICENSE.txt for the
 * QPE freeware source-available terms and THIRD_PARTY_NOTICES.txt for
 * third-party components.
 */
 
#include "landauzener2d.h"
 
#include "core/math/complexfield2d.h"
#include "core/quantum/boundstates1d.h"
#include "core/quantum/coherenttunneling2d.h"
#include "core/quantum/timedependentpotential2d.h"
 
#include <algorithm>
#include <cmath>
 
namespace qpe {
namespace {
 
const double Pi = 3.1415926535897932384626433832795;
 
/**
 * @brief Performs valid pair.
 * @param spectrum The spectrum value used by the operation.
 * @param stateA The state a value used by the operation.
 * @param stateB The state b value used by the operation.
 */
bool validPair(const BoundStateSpectrum1D &spectrum, int stateA, int stateB) noexcept
{
    return spectrum.valid && stateA >= 0 && stateB >= 0 && stateA != stateB &&
           stateA < static_cast<int>(spectrum.states.size()) &&
           stateB < static_cast<int>(spectrum.states.size()) &&
           spectrum.states[static_cast<std::size_t>(stateA)].amplitude.size() == spectrum.x.size() &&
           spectrum.states[static_cast<std::size_t>(stateB)].amplitude.size() == spectrum.x.size();
}
 
/**
 * @brief Performs matrix element x.
 * @param spectrum The spectrum value used by the operation.
 * @param stateA The state a value used by the operation.
 * @param stateB The state b value used by the operation.
 */
double matrixElementX(const BoundStateSpectrum1D &spectrum,
                      int stateA,
                      int stateB) noexcept
{
    if (!validPair(spectrum, stateA, stateB) || !(spectrum.dx > 0.0)) {
        return 0.0;
    }
    const BoundState1D &a = spectrum.states[static_cast<std::size_t>(stateA)];
    const BoundState1D &b = spectrum.states[static_cast<std::size_t>(stateB)];
    double sum = 0.0;
    for (std::size_t i = 0u; i < spectrum.x.size(); ++i) {
        sum += static_cast<double>(a.amplitude[i]) * spectrum.x[i] *
               static_cast<double>(b.amplitude[i]) * spectrum.dx;
    }
    return sum;
}
 
/**
 * @brief Performs project stationary pair.
 * @param field Field data to read or modify.
 * @param spectrum The spectrum value used by the operation.
 * @param stateA The state a value used by the operation.
 * @param stateB The state b value used by the operation.
 * @param ca The ca value used by the operation.
 * @param cb The cb value used by the operation.
 */
bool projectStationaryPair(const ComplexField2D &field,
                           const BoundStateSpectrum1D &spectrum,
                           int stateA,
                           int stateB,
                           std::complex<double> &ca,
                           std::complex<double> &cb) noexcept
{
    if (!validPair(spectrum, stateA, stateB) || field.empty() ||
        field.width() != spectrum.x.size()) {
        return false;
    }
    const double yMeasure = static_cast<double>(field.height()) * static_cast<double>(field.dy());
    if (!(yMeasure > 0.0)) {
        return false;
    }
    const BoundState1D &a = spectrum.states[static_cast<std::size_t>(stateA)];
    const BoundState1D &b = spectrum.states[static_cast<std::size_t>(stateB)];
    const double yScale = 1.0 / std::sqrt(yMeasure);
    const double cellArea = static_cast<double>(field.dx()) * static_cast<double>(field.dy());
    ca = std::complex<double>(0.0, 0.0);
    cb = std::complex<double>(0.0, 0.0);
    for (std::size_t y = 0u; y < field.height(); ++y) {
        for (std::size_t x = 0u; x < field.width(); ++x) {
            const ComplexField2D::Value value = field(x, y);
            const std::complex<double> psi(static_cast<double>(value.real()),
                                           static_cast<double>(value.imag()));
            ca += static_cast<double>(a.amplitude[x]) * yScale * psi * cellArea;
            cb += static_cast<double>(b.amplitude[x]) * yScale * psi * cellArea;
        }
    }
    return true;
}
 
} // namespace
 
LandauZenerTheory2D LandauZener2D::analyze(
    const BoundStateSpectrum1D &zeroBiasSpectrum,
    const LinearPotentialSweep2DConfig &sweep,
    int stateA,
    int stateB,
    double hbar) noexcept
{
    LandauZenerTheory2D result;
    if (!validPair(zeroBiasSpectrum, stateA, stateB) || !(hbar > 0.0) ||
        !sweep.enabled || sweep.parameter != TimePotentialParameter2D::LinearSlopeX ||
        !std::isfinite(sweep.startValue) || !std::isfinite(sweep.endValue) ||
        !std::isfinite(sweep.startTime) || !std::isfinite(sweep.duration) ||
        !(sweep.duration > 0.0) || sweep.startTime < 0.0 ||
        sweep.startValue == sweep.endValue) {
        return result;
    }
 
    result.stateA = stateA;
    result.stateB = stateB;
    result.energyA = zeroBiasSpectrum.states[static_cast<std::size_t>(stateA)].energy;
    result.energyB = zeroBiasSpectrum.states[static_cast<std::size_t>(stateB)].energy;
    result.minimumGap = std::abs(result.energyB - result.energyA);
    result.couplingMatrixElementX = matrixElementX(zeroBiasSpectrum, stateA, stateB);
    result.parameterSweepRate = (sweep.endValue - sweep.startValue) / sweep.duration;
    result.diabaticDetuningRate =
        2.0 * std::abs(result.couplingMatrixElementX * result.parameterSweepRate);
    if (!(result.minimumGap > 1.0e-12) ||
        !(std::abs(result.couplingMatrixElementX) > 1.0e-12) ||
        !(result.diabaticDetuningRate > 1.0e-15)) {
        return result;
    }
 
    const double u = -sweep.startValue / (sweep.endValue - sweep.startValue);
    if (!(u >= 0.0 && u <= 1.0)) {
        return result; // this recipe does not cross zero bias
    }
    result.crossingTime = sweep.startTime + u * sweep.duration;
 
    const double exponent = Pi * result.minimumGap * result.minimumGap /
                            (2.0 * hbar * result.diabaticDetuningRate);
    result.diabaticProbability = std::exp(-std::max(0.0, exponent));
    result.diabaticProbability = std::max(0.0, std::min(1.0, result.diabaticProbability));
    result.adiabaticFollowingProbability = 1.0 - result.diabaticProbability;
    result.leftPhaseRadians = CoherentTunneling2D::preferredLeftPhase(
        zeroBiasSpectrum, stateA, stateB, nullptr);
    result.leftSign = std::cos(result.leftPhaseRadians) >= 0.0 ? 1.0 : -1.0;
    result.valid = true;
    return result;
}
 
LandauZenerProjection2D LandauZener2D::project(
    const ComplexField2D &field,
    const BoundStateSpectrum1D &zeroBiasSpectrum,
    const LandauZenerTheory2D &theory,
    double currentBias)
{
    LandauZenerProjection2D result;
    if (!theory.valid || !std::isfinite(currentBias)) {
        return result;
    }
 
    std::complex<double> ca;
    std::complex<double> cb;
    if (!projectStationaryPair(field, zeroBiasSpectrum,
                               theory.stateA, theory.stateB, ca, cb)) {
        return result;
    }
    result.amplitudeA = ca;
    result.amplitudeB = cb;
    result.capturedPopulation = std::norm(ca) + std::norm(cb);
 
    const double invSqrt2 = 1.0 / std::sqrt(2.0);
    const std::complex<double> cLeft =
        (ca + theory.leftSign * cb) * invSqrt2;
    const std::complex<double> cRight =
        (ca - theory.leftSign * cb) * invSqrt2;
    result.leftPopulation = std::norm(cLeft);
    result.rightPopulation = std::norm(cRight);
 
    // Instantaneous effective Hamiltonian in the zero-bias stationary basis:
    // [ E1, F*x12 ; F*x12, E2 ].
    const double a = theory.energyA;
    const double d = theory.energyB;
    const double g = currentBias * theory.couplingMatrixElementX;
    const double mean = 0.5 * (a + d);
    const double halfDifference = 0.5 * (a - d);
    const double radius = std::sqrt(halfDifference * halfDifference + g * g);
    result.lowerEnergy = mean - radius;
    result.upperEnergy = mean + radius;
 
    double v1 = 1.0;
    double v2 = 0.0;
    if (std::abs(g) > 1.0e-14) {
        v1 = g;
        v2 = result.lowerEnergy - a;
        const double n = std::sqrt(v1 * v1 + v2 * v2);
        if (n > 0.0) {
            v1 /= n;
            v2 /= n;
        }
    } else if (d < a) {
        v1 = 0.0;
        v2 = 1.0;
    }
 
    const std::complex<double> cLower = v1 * ca + v2 * cb;
    const std::complex<double> cUpper = -v2 * ca + v1 * cb;
    result.lowerAdiabaticPopulation = std::norm(cLower);
    result.upperAdiabaticPopulation = std::norm(cUpper);
    result.valid = true;
    return result;
}
 
} // namespace qpe