quantumzitterbewegung2d.cpp

Go to the documentation of this file.

/**
 * @file quantumzitterbewegung2d.cpp
 * @brief Implements the quantumzitterbewegung2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/quantumzitterbewegung2d.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 "quantumzitterbewegung2d.h"
 
#include "core/math/fft2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
 
namespace qpe {
namespace {
constexpr double kTwoPi = 6.283185307179586476925286766559005768;
 
/**
 * @brief Performs finite.
 * @param x X coordinate/value.
 */
bool finite(double x) noexcept { return std::isfinite(x); }
 
/**
 * @brief Performs coordinate.
 * @param i The i value used by the operation.
 * @param n The n value used by the operation.
 * @param d The d value used by the operation.
 */
double coordinate(std::size_t i, std::size_t n, double d) noexcept
{
    return (static_cast<double>(i) - 0.5 * static_cast<double>(n)) * d;
}
 
/**
 * @brief Performs fft k.
 * @param i The i value used by the operation.
 * @param n The n value used by the operation.
 * @param d The d value used by the operation.
 */
double fftK(std::size_t i, std::size_t n, double d) noexcept
{
    const long long m = (i <= n / 2u) ? static_cast<long long>(i)
                                      : static_cast<long long>(i) - static_cast<long long>(n);
    return kTwoPi * static_cast<double>(m) / (static_cast<double>(n) * d);
}
 
/**
 * @brief Performs mean position.
 * @param s The s value used by the operation.
 * @param mx The mx value used by the operation.
 * @param my The my value used by the operation.
 */
void meanPosition(const DiracSpinor2D &s, double &mx, double &my) noexcept
{
    mx = my = 0.0;
    if (s.empty()) return;
    const double cell = static_cast<double>(s.dx()) * static_cast<double>(s.dy());
    double n = 0.0;
    for (std::size_t y = 0; y < s.height(); ++y) {
        const double yy = coordinate(y, s.height(), s.dy());
        for (std::size_t x = 0; x < s.width(); ++x) {
            const double xx = coordinate(x, s.width(), s.dx());
            const double rho = std::norm(s.upper()(x,y)) + std::norm(s.lower()(x,y));
            const double w = rho * cell;
            n += w;
            mx += xx * w;
            my += yy * w;
        }
    }
    if (n > 0.0) { mx /= n; my /= n; }
}
}
 
bool QuantumZitterbewegung2D::validateConfig(const QuantumZitterbewegung2DConfig &c,
                                              std::string *reason) noexcept
{
    const bool ok = finite(c.lightSpeed) && finite(c.restMass) && finite(c.centerX) &&
                    finite(c.centerY) && finite(c.sigma) && finite(c.momentumX) &&
                    finite(c.momentumY) && finite(c.analysisEndTime) &&
                    c.lightSpeed >= 0.05 && c.lightSpeed <= 20.0 &&
                    c.restMass >= 0.02 && c.restMass <= 20.0 &&
                    c.sigma >= 0.25 && c.sigma <= 20.0 &&
                    std::abs(c.momentumX) <= 20.0 && std::abs(c.momentumY) <= 20.0 &&
                    std::abs(c.centerX) <= 1000.0 && std::abs(c.centerY) <= 1000.0 &&
                    (c.preparationMode == 0 || c.preparationMode == 1) &&
                    c.analysisEndTime >= 0.2 && c.analysisEndTime <= 200.0;
    if (!ok && reason) *reason = "QuantumZitterbewegung2D: invalid configuration";
    return ok;
}
 
bool QuantumZitterbewegung2D::prepareState(DiracSpinor2D &s,
                                            const QuantumZitterbewegung2DConfig &c,
                                            double hbar,
                                            std::string *reason)
{
    if (s.empty() || !FFT2D::supportsSize(s.width(), s.height()) ||
        !(hbar > 0.0) || !finite(hbar) || !validateConfig(c, reason)) {
        if (reason && reason->empty()) *reason = "QuantumZitterbewegung2D: invalid spinor or hbar";
        return false;
    }
 
    const double inv4Sigma2 = 1.0 / (4.0 * c.sigma * c.sigma);
    const double kx0 = c.momentumX / hbar;
    const double ky0 = c.momentumY / hbar;
    for (std::size_t y = 0; y < s.height(); ++y) {
        const double yy = coordinate(y, s.height(), s.dy());
        for (std::size_t x = 0; x < s.width(); ++x) {
            const double xx = coordinate(x, s.width(), s.dx());
            const double dx = xx - c.centerX;
            const double dy = yy - c.centerY;
            const double amp = std::exp(-(dx*dx + dy*dy) * inv4Sigma2);
            const double ph = kx0 * dx + ky0 * dy;
            const std::complex<double> a = amp * std::complex<double>(std::cos(ph), std::sin(ph));
            s.upper()(x,y) = ComplexField2D::Value(static_cast<float>(a.real()), static_cast<float>(a.imag()));
            s.lower()(x,y) = ComplexField2D::Value(0.0f, 0.0f);
        }
    }
 
    if (c.preparationMode == 1) {
        // Project each Fourier component onto the positive-energy eigenspace:
        // P+ = (I + H/E)/2.  Starting from (a,0) gives the two entries below.
        FFT2D::forward(s.upper());
        FFT2D::forward(s.lower());
        const double cc = c.lightSpeed;
        const double rest = c.restMass * cc * cc;
        for (std::size_t y = 0; y < s.height(); ++y) {
            const double py = hbar * fftK(y, s.height(), s.dy());
            for (std::size_t x = 0; x < s.width(); ++x) {
                const double px = hbar * fftK(x, s.width(), s.dx());
                const double ax = cc * px;
                const double ay = cc * py;
                const double e = std::sqrt(ax*ax + ay*ay + rest*rest);
                const std::complex<double> a(s.upper()(x,y).real(), s.upper()(x,y).imag());
                const std::complex<double> up = 0.5 * (1.0 + rest / e) * a;
                const std::complex<double> lo = 0.5 * std::complex<double>(ax, ay) / e * a;
                s.upper()(x,y) = ComplexField2D::Value(static_cast<float>(up.real()), static_cast<float>(up.imag()));
                s.lower()(x,y) = ComplexField2D::Value(static_cast<float>(lo.real()), static_cast<float>(lo.imag()));
            }
        }
        FFT2D::inverse(s.upper());
        FFT2D::inverse(s.lower());
    }
 
    if (!s.normalize()) {
        if (reason) *reason = "QuantumZitterbewegung2D: normalization failed";
        return false;
    }
    return true;
}
 
DiracEnergyFractions2D QuantumZitterbewegung2D::energyFractions(
    const DiracSpinor2D &s, const QuantumZitterbewegung2DConfig &c, double hbar) noexcept
{
    DiracEnergyFractions2D out;
    if (s.empty() || !(hbar > 0.0) || !finite(hbar)) return out;
    try {
        ComplexField2D u = s.upper();
        ComplexField2D v = s.lower();
        FFT2D::forward(u);
        FFT2D::forward(v);
        const double cc = c.lightSpeed;
        const double rest = c.restMass * cc * cc;
        double pos = 0.0, neg = 0.0;
        for (std::size_t y = 0; y < s.height(); ++y) {
            const double py = hbar * fftK(y, s.height(), s.dy());
            for (std::size_t x = 0; x < s.width(); ++x) {
                const double px = hbar * fftK(x, s.width(), s.dx());
                const double ax = cc * px;
                const double ay = cc * py;
                const double e = std::sqrt(ax*ax + ay*ay + rest*rest);
                const std::complex<double> a(u(x,y).real(), u(x,y).imag());
                const std::complex<double> b(v(x,y).real(), v(x,y).imag());
                const double rho = std::norm(a) + std::norm(b);
                const std::complex<double> h01(ax, -ay);
                const std::complex<double> h10(ax, ay);
                const std::complex<double> ha = rest * a + h01 * b;
                const std::complex<double> hb = h10 * a - rest * b;
                const double hExp = std::real(std::conj(a)*ha + std::conj(b)*hb);
                pos += 0.5 * (rho + hExp/e);
                neg += 0.5 * (rho - hExp/e);
            }
        }
        const double total = pos + neg;
        if (total > 0.0) { out.positive = pos/total; out.negative = neg/total; }
    } catch (...) {}
    return out;
}
 
double QuantumZitterbewegung2D::centralEnergy(const QuantumZitterbewegung2DConfig &c) noexcept
{
    const double rest = c.restMass * c.lightSpeed * c.lightSpeed;
    const double cp2 = c.lightSpeed*c.lightSpeed *
                       (c.momentumX*c.momentumX + c.momentumY*c.momentumY);
    return std::sqrt(cp2 + rest*rest);
}
 
double QuantumZitterbewegung2D::predictedZitterOmega(const QuantumZitterbewegung2DConfig &c,
                                                      double hbar) noexcept
{
    return hbar > 0.0 ? 2.0 * centralEnergy(c) / hbar : 0.0;
}
 
double QuantumZitterbewegung2D::predictedZitterPeriod(const QuantumZitterbewegung2DConfig &c,
                                                       double hbar) noexcept
{
    const double w = predictedZitterOmega(c, hbar);
    return w > 0.0 ? kTwoPi / w : 0.0;
}
 
double QuantumZitterbewegung2D::predictedDriftX(const QuantumZitterbewegung2DConfig &c) noexcept
{
    const double e = centralEnergy(c);
    if (!(e > 0.0)) return 0.0;
    if (c.preparationMode == 1)
        return c.lightSpeed*c.lightSpeed*c.momentumX/e;
    return c.restMass*std::pow(c.lightSpeed,4.0)*c.momentumX/(e*e);
}
 
double QuantumZitterbewegung2D::predictedDriftY(const QuantumZitterbewegung2DConfig &c) noexcept
{
    const double e = centralEnergy(c);
    if (!(e > 0.0)) return 0.0;
    if (c.preparationMode == 1)
        return c.lightSpeed*c.lightSpeed*c.momentumY/e;
    return c.restMass*std::pow(c.lightSpeed,4.0)*c.momentumY/(e*e);
}
 
QuantumZitterbewegungSample2D QuantumZitterbewegung2D::measure(
    const DiracSpinor2D &s, const QuantumZitterbewegung2DConfig &c, double hbar,
    double time, double x0, double y0) noexcept
{
    QuantumZitterbewegungSample2D out;
    out.time = time;
    if (s.empty()) return out;
    meanPosition(s, out.meanX, out.meanY);
    const double cell = static_cast<double>(s.dx()) * static_cast<double>(s.dy());
    double vx = 0.0, vy = 0.0, n = 0.0;
    for (std::size_t i = 0; i < s.upper().size(); ++i) {
        const std::complex<double> a(s.upper().data()[i].real(), s.upper().data()[i].imag());
        const std::complex<double> b(s.lower().data()[i].real(), s.lower().data()[i].imag());
        const double w = (std::norm(a)+std::norm(b))*cell;
        n += w;
        vx += c.lightSpeed * 2.0 * std::real(std::conj(a)*b) * cell;
        vy += c.lightSpeed * 2.0 * std::imag(std::conj(a)*b) * cell;
    }
    if (n > 0.0) { vx /= n; vy /= n; }
    out.meanVx = vx;
    out.meanVy = vy;
    out.residualX = out.meanX - (x0 + predictedDriftX(c)*time);
    out.residualY = out.meanY - (y0 + predictedDriftY(c)*time);
    out.upperPopulation = s.upperPopulation();
    out.lowerPopulation = s.lowerPopulation();
    out.norm = s.norm2();
    (void)hbar;
    return out;
}
 
void QuantumZitterbewegungHistory2D::clear()
{
    configured_ = false;
    config_ = QuantumZitterbewegung2DConfig();
    hbar_ = 1.0;
    initialMeanX_ = initialMeanY_ = 0.0;
    maximumAbsResidualX_ = 0.0;
    fittedDriftX_ = fittedDriftY_ = 0.0;
    energyFractions_ = DiracEnergyFractions2D();
    samples_.clear();
}
 
bool QuantumZitterbewegungHistory2D::configure(const DiracSpinor2D &s,
                                                const QuantumZitterbewegung2DConfig &c,
                                                double hbar)
{
    clear();
    if (s.empty() || !QuantumZitterbewegung2D::validateConfig(c, nullptr) ||
        !(hbar > 0.0) || !finite(hbar)) return false;
    config_ = c;
    hbar_ = hbar;
    meanPosition(s, initialMeanX_, initialMeanY_);
    energyFractions_ = QuantumZitterbewegung2D::energyFractions(s, c, hbar);
    configured_ = true;
    return true;
}
 
void QuantumZitterbewegungHistory2D::append(double time, const DiracSpinor2D &s)
{
    if (!configured_ || s.empty()) return;
    samples_.push_back(QuantumZitterbewegung2D::measure(
        s, config_, hbar_, time, initialMeanX_, initialMeanY_));
    trim();
 
    // Detrend against the best linear trajectory actually carried by the
    // finite-width packet.  Using only the central-momentum group velocity
    // creates a false residual for the positive-energy control because its
    // momentum distribution has finite width.  A least-squares secular line
    // leaves the genuinely oscillatory Zitterbewegung component while making
    // the positive-energy control essentially flat.
    const std::size_t n = samples_.size();
    if (n < 2u) {
        fittedDriftX_ = fittedDriftY_ = 0.0;
        maximumAbsResidualX_ = 0.0;
        if (!samples_.empty()) { samples_.back().residualX = 0.0; samples_.back().residualY = 0.0; }
        return;
    }
 
    double st = 0.0, stt = 0.0, sx = 0.0, sy = 0.0, stx = 0.0, sty = 0.0;
    for (const auto &a : samples_) {
        st += a.time; stt += a.time * a.time;
        sx += a.meanX; sy += a.meanY;
        stx += a.time * a.meanX; sty += a.time * a.meanY;
    }
    const double dn = static_cast<double>(n);
    const double denominator = dn * stt - st * st;
    if (!(std::abs(denominator) > 1.0e-18)) return;
    fittedDriftX_ = (dn * stx - st * sx) / denominator;
    fittedDriftY_ = (dn * sty - st * sy) / denominator;
    const double interceptX = (sx - fittedDriftX_ * st) / dn;
    const double interceptY = (sy - fittedDriftY_ * st) / dn;
    maximumAbsResidualX_ = 0.0;
    for (auto &a : samples_) {
        a.residualX = a.meanX - (interceptX + fittedDriftX_ * a.time);
        a.residualY = a.meanY - (interceptY + fittedDriftY_ * a.time);
        maximumAbsResidualX_ = std::max(maximumAbsResidualX_, std::abs(a.residualX));
    }
}
 
void QuantumZitterbewegungHistory2D::trim()
{
    if (maximumSamples_ == 0u) { samples_.clear(); return; }
    if (samples_.size() > maximumSamples_)
        samples_.erase(samples_.begin(), samples_.begin() + static_cast<std::ptrdiff_t>(samples_.size()-maximumSamples_));
}
 
} // namespace qpe