ramseyinterferometry2d.cpp

Go to the documentation of this file.

/**
 * @file ramseyinterferometry2d.cpp
 * @brief Implements the ramseyinterferometry2d physics/numerical component used by QPE experiments.
 * @ingroup qpe_quantum
 *
 * Source: @c core/quantum/ramseyinterferometry2d.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 "ramseyinterferometry2d.h"
 
#include "core/quantum/boundstates1d.h"
#include "core/quantum/timedependentpotential2d.h"
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <limits>
#include <vector>
 
namespace qpe {
namespace {
 
const double Pi = 3.1415926535897932384626433832795;
const double EnsembleSpanSigma = 3.0;
 
struct State2
{
    std::complex<double> a;
    std::complex<double> b;
};
 
struct Operator2
{
    std::complex<double> aa;
    std::complex<double> ab;
    std::complex<double> ba;
    std::complex<double> bb;
};
 
struct PreparedMember
{
    double energyA = 0.0;
    double energyB = 0.0;
    double coupling = 0.0;
    double weight = 1.0;
    Operator2 firstPulse;
    Operator2 secondPulse;
    Operator2 echoPulse;
};
 
/**
 * @brief Performs finite.
 * @param v The v value used by the operation.
 */
bool finite(double v) noexcept
{
    return std::isfinite(v);
}
 
/**
 * @brief Performs valid pair.
 * @param spectrum The spectrum value used by the operation.
 * @param a The a value used by the operation.
 * @param b The b value used by the operation.
 */
bool validPair(const BoundStateSpectrum1D &spectrum, int a, int b) noexcept
{
    return spectrum.valid && a >= 0 && b >= 0 && a != b &&
           a < static_cast<int>(spectrum.states.size()) &&
           b < static_cast<int>(spectrum.states.size());
}
 
/**
 * @brief Performs propagate constant.
 * @param state The state value used by the operation.
 * @param energyA The energy a value used by the operation.
 * @param energyB The energy b value used by the operation.
 * @param coupling The coupling value used by the operation.
 * @param driveValue The drive value value used by the operation.
 * @param dt Time increment in solver units.
 * @param hbar The hbar value used by the operation.
 */
State2 propagateConstant(const State2 &state,
                         double energyA,
                         double energyB,
                         double coupling,
                         double driveValue,
                         double dt,
                         double hbar) noexcept
{
    const double mean = 0.5 * (energyA + energyB);
    const double hz = 0.5 * (energyA - energyB);
    const double hx = coupling * driveValue;
    const double radius = std::sqrt(hx * hx + hz * hz);
    const std::complex<double> global =
        std::exp(std::complex<double>(0.0, -mean * dt / hbar));
    if (!(radius > 1.0e-15)) {
        State2 out;
        out.a = global * state.a;
        out.b = global * state.b;
        return out;
    }
    const double angle = radius * dt / hbar;
    const double c = std::cos(angle);
    const double s = std::sin(angle) / radius;
    const std::complex<double> minusI(0.0, -1.0);
    State2 out;
    out.a = global * (c * state.a + minusI * s * (hz * state.a + hx * state.b));
    out.b = global * (c * state.b + minusI * s * (hx * state.a - hz * state.b));
    return out;
}
 
/**
 * @brief Applies operator.
 * @param op The op value used by the operation.
 * @param state The state value used by the operation.
 */
State2 applyOperator(const Operator2 &op, const State2 &state) noexcept
{
    State2 out;
    out.a = op.aa * state.a + op.ab * state.b;
    out.b = op.ba * state.a + op.bb * state.b;
    return out;
}
 
/**
 * @brief Performs pulse operator.
 * @param energyA The energy a value used by the operation.
 * @param energyB The energy b value used by the operation.
 * @param coupling The coupling value used by the operation.
 * @param sequence The sequence value used by the operation.
 * @param duration The duration value used by the operation.
 * @param phaseDegrees The phase degrees value used by the operation.
 * @param hbar The hbar value used by the operation.
 * @param maximumStep The maximum step value used by the operation.
 */
Operator2 pulseOperator(double energyA,
                        double energyB,
                        double coupling,
                        const RamseyPulseSequence2DConfig &sequence,
                        double duration,
                        double phaseDegrees,
                        double hbar,
                        double maximumStep) noexcept
{
    Operator2 op;
    State2 columnA;
    columnA.a = std::complex<double>(1.0, 0.0);
    columnA.b = std::complex<double>(0.0, 0.0);
    State2 columnB;
    columnB.a = std::complex<double>(0.0, 0.0);
    columnB.b = std::complex<double>(1.0, 0.0);
 
    if (!(duration > 0.0)) {
        op.aa = columnA.a;
        op.ba = columnA.b;
        op.ab = columnB.a;
        op.bb = columnB.b;
        return op;
    }
 
    const int steps = std::max(80, static_cast<int>(std::ceil(duration / maximumStep)));
    const double dt = duration / static_cast<double>(steps);
    const double phase0 = phaseDegrees * Pi / 180.0;
    for (int i = 0; i < steps; ++i) {
        const double localTime = (static_cast<double>(i) + 0.5) * dt;
        const double phase = sequence.angularFrequency * localTime + phase0;
        const double drive = sequence.meanValue + sequence.amplitude * std::sin(phase);
        columnA = propagateConstant(columnA, energyA, energyB, coupling, drive, dt, hbar);
        columnB = propagateConstant(columnB, energyA, energyB, coupling, drive, dt, hbar);
    }
 
    op.aa = columnA.a;
    op.ba = columnA.b;
    op.ab = columnB.a;
    op.bb = columnB.b;
    return op;
}
 
/**
 * @brief Prepares member.
 * @param spectrum The spectrum value used by the operation.
 * @param sequence The sequence 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 hbar The hbar value used by the operation.
 * @param transitionOffset The transition offset value used by the operation.
 * @param weight The weight value used by the operation.
 * @param maximumPulseStep The maximum pulse step value used by the operation.
 */
PreparedMember prepareMember(const BoundStateSpectrum1D &spectrum,
                             const RamseyPulseSequence2DConfig &sequence,
                             int stateA,
                             int stateB,
                             double hbar,
                             double transitionOffset,
                             double weight,
                             double maximumPulseStep) noexcept
{
    PreparedMember member;
    const double baseA = spectrum.states[static_cast<std::size_t>(stateA)].energy;
    const double baseB = spectrum.states[static_cast<std::size_t>(stateB)].energy;
    // Keep the mean energy fixed and shift only the transition frequency.
    const double halfShift = 0.5 * hbar * transitionOffset;
    member.energyA = baseA - halfShift;
    member.energyB = baseB + halfShift;
    member.coupling = DrivenTwoLevel2D::matrixElementX(spectrum, stateA, stateB);
    member.weight = weight;
    member.firstPulse = pulseOperator(member.energyA, member.energyB, member.coupling,
                                      sequence, sequence.pulseDuration,
                                      sequence.firstPulsePhaseDegrees,
                                      hbar, maximumPulseStep);
    member.secondPulse = pulseOperator(member.energyA, member.energyB, member.coupling,
                                       sequence, sequence.pulseDuration,
                                       sequence.secondPulsePhaseDegrees,
                                       hbar, maximumPulseStep);
    if (sequence.echoEnabled) {
        member.echoPulse = pulseOperator(member.energyA, member.energyB, member.coupling,
                                         sequence, sequence.echoPulseDuration,
                                         sequence.echoPulsePhaseDegrees,
                                         hbar, maximumPulseStep);
    }
    return member;
}
 
/**
 * @brief Performs predict prepared.
 * @param member The member value used by the operation.
 * @param sequence The sequence value used by the operation.
 * @param freeEvolutionDuration The free evolution duration value used by the operation.
 * @param hbar The hbar value used by the operation.
 * @param populationA The population a value used by the operation.
 * @param populationB The population b value used by the operation.
 */
void predictPrepared(const PreparedMember &member,
                     const RamseyPulseSequence2DConfig &sequence,
                     double freeEvolutionDuration,
                     double hbar,
                     double &populationA,
                     double &populationB) noexcept
{
    State2 state;
    state.a = std::complex<double>(1.0, 0.0);
    state.b = std::complex<double>(0.0, 0.0);
    state = applyOperator(member.firstPulse, state);
 
    if (sequence.echoEnabled) {
        const double halfFree = 0.5 * freeEvolutionDuration;
        if (halfFree > 0.0) {
            state = propagateConstant(state, member.energyA, member.energyB,
                                      member.coupling, sequence.meanValue,
                                      halfFree, hbar);
        }
        state = applyOperator(member.echoPulse, state);
        if (halfFree > 0.0) {
            state = propagateConstant(state, member.energyA, member.energyB,
                                      member.coupling, sequence.meanValue,
                                      halfFree, hbar);
        }
    } else if (freeEvolutionDuration > 0.0) {
        state = propagateConstant(state, member.energyA, member.energyB,
                                  member.coupling, sequence.meanValue,
                                  freeEvolutionDuration, hbar);
    }
 
    state = applyOperator(member.secondPulse, state);
    populationA = std::norm(state.a);
    populationB = std::norm(state.b);
}
 
/**
 * @brief Performs sequence valid for reference.
 * @param sequence The sequence value used by the operation.
 */
bool sequenceValidForReference(const RamseyPulseSequence2DConfig &sequence) noexcept
{
    return sequence.enabled &&
           sequence.parameter == TimePotentialParameter2D::LinearSlopeX &&
           finite(sequence.meanValue) && finite(sequence.amplitude) &&
           finite(sequence.angularFrequency) && sequence.angularFrequency > 0.0 &&
           finite(sequence.pulseDuration) && sequence.pulseDuration > 0.0 &&
           finite(sequence.freeEvolutionDuration) && sequence.freeEvolutionDuration >= 0.0 &&
           finite(sequence.echoPulseDuration) &&
           (!sequence.echoEnabled || sequence.echoPulseDuration > 0.0) &&
           finite(sequence.firstPulsePhaseDegrees) &&
           finite(sequence.secondPulsePhaseDegrees) &&
           finite(sequence.echoPulsePhaseDegrees);
}
 
std::vector<PreparedMember> prepareEnsemble(const BoundStateSpectrum1D &spectrum,
                                            const RamseyPulseSequence2DConfig &sequence,
                                            const RamseyEnsemble2DConfig &ensemble,
                                            int stateA,
                                            int stateB,
                                            double hbar)
{
    std::vector<PreparedMember> members;
    if (!validPair(spectrum, stateA, stateB) || !sequenceValidForReference(sequence) ||
        !RamseyInterferometry2D::validateEnsembleConfig(ensemble) || !(hbar > 0.0)) {
        return members;
    }
 
    const int count = ensemble.detuningStdDev > 0.0 ? ensemble.memberCount : 1;
    members.reserve(static_cast<std::size_t>(count));
    double totalWeight = 0.0;
    for (int i = 0; i < count; ++i) {
        double z = 0.0;
        if (count > 1) {
            z = -EnsembleSpanSigma + 2.0 * EnsembleSpanSigma *
                static_cast<double>(i) / static_cast<double>(count - 1);
        }
        const double weight = std::exp(-0.5 * z * z);
        const double offset = z * ensemble.detuningStdDev;
        // Step 28 is diagnostic only.  A 0.01 pulse step is sufficiently fine
        // for the ensemble overlay and makes interactive scans inexpensive.
        members.push_back(prepareMember(spectrum, sequence, stateA, stateB, hbar,
                                        offset, weight, 0.01));
        totalWeight += weight;
    }
    if (!(totalWeight > 0.0)) {
        members.clear();
        return members;
    }
    for (PreparedMember &member : members) {
        member.weight /= totalWeight;
    }
    return members;
}
 
/**
 * @brief Performs average prepared.
 * @param members The members value used by the operation.
 * @param sequence The sequence value used by the operation.
 * @param freeEvolutionDuration The free evolution duration value used by the operation.
 * @param hbar The hbar value used by the operation.
 * @param meanA The mean a value used by the operation.
 * @param meanB The mean b value used by the operation.
 * @param minimumB The minimum b value used by the operation.
 * @param maximumB The maximum b value used by the operation.
 */
void averagePrepared(const std::vector<PreparedMember> &members,
                     const RamseyPulseSequence2DConfig &sequence,
                     double freeEvolutionDuration,
                     double hbar,
                     double &meanA,
                     double &meanB,
                     double *minimumB,
                     double *maximumB) noexcept
{
    meanA = 0.0;
    meanB = 0.0;
    double minB = std::numeric_limits<double>::infinity();
    double maxB = -std::numeric_limits<double>::infinity();
    for (const PreparedMember &member : members) {
        double a = 0.0;
        double b = 0.0;
        predictPrepared(member, sequence, freeEvolutionDuration, hbar, a, b);
        meanA += member.weight * a;
        meanB += member.weight * b;
        minB = std::min(minB, b);
        maxB = std::max(maxB, b);
    }
    if (minimumB) *minimumB = std::isfinite(minB) ? minB : 0.0;
    if (maximumB) *maximumB = std::isfinite(maxB) ? maxB : 0.0;
}
 
} // namespace
 
RamseyTheory2D RamseyInterferometry2D::analyze(
    const BoundStateSpectrum1D &spectrum,
    const RamseyPulseSequence2DConfig &sequence,
    int stateA,
    int stateB,
    double hbar) noexcept
{
    RamseyTheory2D result;
    if (!validPair(spectrum, stateA, stateB) || !(hbar > 0.0) ||
        !sequenceValidForReference(sequence)) {
        return result;
    }
 
    result.stateA = stateA;
    result.stateB = stateB;
    result.energyA = spectrum.states[static_cast<std::size_t>(stateA)].energy;
    result.energyB = spectrum.states[static_cast<std::size_t>(stateB)].energy;
    result.transitionAngularFrequency = std::abs(result.energyB - result.energyA) / hbar;
    result.driveAngularFrequency = sequence.angularFrequency;
    result.detuning = sequence.angularFrequency - result.transitionAngularFrequency;
    result.couplingMatrixElementX = DrivenTwoLevel2D::matrixElementX(spectrum, stateA, stateB);
    result.resonantRabiAngularFrequency =
        std::abs(sequence.amplitude * result.couplingMatrixElementX) / hbar;
    if (result.resonantRabiAngularFrequency > 0.0) {
        result.idealPiOver2Duration = Pi / (2.0 * result.resonantRabiAngularFrequency);
    }
    result.firstPulseEnd = sequence.startTime + sequence.pulseDuration;
    result.echoEnabled = sequence.echoEnabled;
    if (sequence.echoEnabled) {
        result.echoPulseStart = result.firstPulseEnd + 0.5 * sequence.freeEvolutionDuration;
        result.echoPulseEnd = result.echoPulseStart + sequence.echoPulseDuration;
        result.secondPulseStart = result.echoPulseEnd + 0.5 * sequence.freeEvolutionDuration;
    } else {
        result.echoPulseStart = 0.0;
        result.echoPulseEnd = 0.0;
        result.secondPulseStart = result.firstPulseEnd + sequence.freeEvolutionDuration;
    }
    result.secondPulseEnd = result.secondPulseStart + sequence.pulseDuration;
    result.freePhaseRadians = result.transitionAngularFrequency * sequence.freeEvolutionDuration;
 
    const PreparedMember nominal = prepareMember(spectrum, sequence, stateA, stateB,
                                                 hbar, 0.0, 1.0, 0.0025);
    predictPrepared(nominal, sequence, sequence.freeEvolutionDuration, hbar,
                    result.predictedFinalPopulationA,
                    result.predictedFinalPopulationB);
    result.valid = true;
    return result;
}
 
DrivenTwoLevelProjection2D RamseyInterferometry2D::project(
    const ComplexField2D &field,
    const BoundStateSpectrum1D &spectrum,
    int stateA,
    int stateB)
{
    return DrivenTwoLevel2D::project(field, spectrum, stateA, stateB);
}
 
std::vector<RamseyFringeSample2D> RamseyInterferometry2D::scanFreeEvolution(
    const BoundStateSpectrum1D &spectrum,
    const RamseyPulseSequence2DConfig &sequence,
    double minimumFreeTime,
    double maximumFreeTime,
    int sampleCount,
    int stateA,
    int stateB,
    double hbar)
{
    std::vector<RamseyFringeSample2D> result;
    if (!validPair(spectrum, stateA, stateB) || !sequenceValidForReference(sequence) ||
        !finite(minimumFreeTime) || !finite(maximumFreeTime) ||
        minimumFreeTime < 0.0 || !(maximumFreeTime > minimumFreeTime) ||
        sampleCount < 2 || !(hbar > 0.0)) {
        return result;
    }
 
    const double omega12 = std::abs(
        spectrum.states[static_cast<std::size_t>(stateB)].energy -
        spectrum.states[static_cast<std::size_t>(stateA)].energy) / hbar;
    const PreparedMember nominal = prepareMember(spectrum, sequence, stateA, stateB,
                                                 hbar, 0.0, 1.0, 0.0025);
    result.reserve(static_cast<std::size_t>(sampleCount));
    for (int i = 0; i < sampleCount; ++i) {
        const double u = static_cast<double>(i) / static_cast<double>(sampleCount - 1);
        RamseyFringeSample2D sample;
        sample.freeEvolutionDuration =
            minimumFreeTime + (maximumFreeTime - minimumFreeTime) * u;
        sample.phaseRadians = omega12 * sample.freeEvolutionDuration;
        predictPrepared(nominal, sequence, sample.freeEvolutionDuration, hbar,
                        sample.populationA, sample.populationB);
        result.push_back(sample);
    }
    return result;
}
 
bool RamseyInterferometry2D::validateEnsembleConfig(
    const RamseyEnsemble2DConfig &config) noexcept
{
    return finite(config.detuningStdDev) && config.detuningStdDev >= 0.0 &&
           config.detuningStdDev <= 10.0 && config.memberCount >= 3 &&
           config.memberCount <= 201;
}
 
RamseyEnsembleResult2D RamseyInterferometry2D::analyzeStaticDetuningEnsemble(
    const BoundStateSpectrum1D &spectrum,
    const RamseyPulseSequence2DConfig &sequence,
    const RamseyEnsemble2DConfig &ensemble,
    int stateA,
    int stateB,
    double hbar)
{
    RamseyEnsembleResult2D result;
    if (!ensemble.enabled || !validateEnsembleConfig(ensemble) ||
        !validPair(spectrum, stateA, stateB) || !sequenceValidForReference(sequence) ||
        !(hbar > 0.0)) {
        return result;
    }
 
    const std::vector<PreparedMember> members =
        prepareEnsemble(spectrum, sequence, ensemble, stateA, stateB, hbar);
    if (members.empty()) return result;
 
    averagePrepared(members, sequence, sequence.freeEvolutionDuration, hbar,
                    result.meanPopulationA, result.meanPopulationB,
                    &result.minimumPopulationB, &result.maximumPopulationB);
    result.valid = true;
    result.memberCount = ensemble.detuningStdDev > 0.0
        ? static_cast<int>(members.size()) : 1;
    result.detuningStdDev = ensemble.detuningStdDev;
    return result;
}
 
std::vector<RamseyEnsembleFringeSample2D>
RamseyInterferometry2D::scanFreeEvolutionEnsemble(
    const BoundStateSpectrum1D &spectrum,
    const RamseyPulseSequence2DConfig &sequence,
    const RamseyEnsemble2DConfig &ensemble,
    double minimumFreeTime,
    double maximumFreeTime,
    int sampleCount,
    int stateA,
    int stateB,
    double hbar)
{
    std::vector<RamseyEnsembleFringeSample2D> result;
    if (!ensemble.enabled || !validateEnsembleConfig(ensemble) ||
        !validPair(spectrum, stateA, stateB) || !sequenceValidForReference(sequence) ||
        !finite(minimumFreeTime) || !finite(maximumFreeTime) ||
        minimumFreeTime < 0.0 || !(maximumFreeTime > minimumFreeTime) ||
        sampleCount < 2 || !(hbar > 0.0)) {
        return result;
    }
 
    const std::vector<PreparedMember> members =
        prepareEnsemble(spectrum, sequence, ensemble, stateA, stateB, hbar);
    if (members.empty()) return result;
 
    result.reserve(static_cast<std::size_t>(sampleCount));
    for (int i = 0; i < sampleCount; ++i) {
        const double u = static_cast<double>(i) / static_cast<double>(sampleCount - 1);
        RamseyEnsembleFringeSample2D sample;
        sample.freeEvolutionDuration =
            minimumFreeTime + (maximumFreeTime - minimumFreeTime) * u;
        averagePrepared(members, sequence, sample.freeEvolutionDuration, hbar,
                        sample.meanPopulationA, sample.meanPopulationB,
                        nullptr, nullptr);
        result.push_back(sample);
    }
    return result;
}
 
} // namespace qpe