차세대 스핀트로닉스 SOT-MRAM 스핀궤도토크 자화반전 & 1,000Hz 0-Alloc C# 란다우-리프시츠-길버트(LLGS) 동역학 시뮬레이터

using System;
using System.Runtime.CompilerServices;

namespace ZeliDesk.Stem.Spintronics;

public readonly ref struct SpinSwitchingState
{
    public readonly double MagnetizationZ;
    public readonly double SwitchingTimePicoseconds;
    public readonly double CriticalCurrentDensityRatio;
    public readonly double ThermalStabilityDelta;
    public readonly bool IsSwitchedSuccess;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public SpinSwitchingState(double mz, double timePs, double jRatio, double delta, bool switched)
    {
        MagnetizationZ = mz;
        SwitchingTimePicoseconds = timePs;
        CriticalCurrentDensityRatio = jRatio;
        ThermalStabilityDelta = delta;
        IsSwitchedSuccess = switched;
    }
}

public ref struct SotMramSimulator
{
    private const double Gamma = 1.76e11; // rad/(s*T)
    private const double GilbertAlpha = 0.02; // 감쇠 상수
    private const double Hbar = 1.0545718e-34;
    private const double ElectronCharge = 1.602176634e-19;
    private const double Mu0 = 1.256637e-6;

    // 자성 자유층 물성 (CoFeB)
    private const double Ms = 1.0e6; // A/m 포화자화
    private const double Hк = 0.35; // 0.35 T 수직이방성 자계
    private const double ThicknessF = 1.2e-9; // 1.2 nm 자유층 두께
    private const double SpinHallTheta = 0.35; // 텅스텐(W) beta-상 스핀홀 각도

    private double _mx, _my, _mz;

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public SotMramSimulator(double initialMz = 0.99)
    {
        _mz = initialMz;
        double trans = Math.Sqrt(Math.Max(0.0001, 1.0 - _mz * _mz));
        _mx = trans;
        _my = 0.0;
    }

    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    public SpinSwitchingState SimulatePulse(double appliedCurrentDensityMaCm2, double pulseWidthPs = 500.0)
    {
        double currentDensityA_m2 = appliedCurrentDensityMaCm2 * 1e4; // A/m^2

        // 1. 댐핑 유사(Damping-Like) 스핀 토크 크기 산출
        double tauDl = (Hbar * SpinHallTheta * currentDensityA_m2) / (2.0 * ElectronCharge * Mu0 * Ms * ThicknessF);

        // 2. 임계 스위칭 전류 밀도 산출 (단위: MA/cm^2)
        double jCritA_m2 = (2.0 * ElectronCharge * Mu0 * Ms * ThicknessF * GilbertAlpha * Hк) / (Hbar * SpinHallTheta);
        double jCritMaCm2 = jCritA_m2 * 1e-4;
        double jRatio = appliedCurrentDensityMaCm2 / Math.Max(1.0, jCritMaCm2);

        // 3. 피코초 단위 LLGS 적분 (dt = 0.2 ps, 최대 2,500 스텝)
        double dtSec = 0.2e-12;
        int maxSteps = (int)(pulseWidthPs / 0.2);
        double switchedTimePs = 0.0;
        bool switched = false;

        // 스핀 분극 방향 sigma = +y (면내 직교 주입)
        double sigX = 0.0, sigY = 1.0, sigZ = 0.0;

        for (int step = 0; step < maxSteps; step++)
        {
            // 유효 자기장 (수직 이방성 + 면내 외부 보조 자계 0.02T)
            double heffX = 0.02; // 무자장 반전을 위한 미세 면내 자계
            double heffY = 0.0;
            double heffZ = Hк * _mz;

            // 라모어 세차 토크 T_prec = -gamma * (m x Heff)
            double tPrecX = -Gamma * (_my * heffZ - _mz * heffY);
            double tPrecY = -Gamma * (_mz * heffX - _mx * heffZ);
            double tPrecZ = -Gamma * (_mx * heffY - _my * heffX);

            // SOT 댐핑 토크 T_sot = -gamma * tau_dl * (m x (m x sigma))
            // m x sigma
            double mCrossSigX = _my * sigZ - _mz * sigY;
            double mCrossSigY = _mz * sigX - _mx * sigZ;
            double mCrossSigZ = _mx * sigY - _my * sigX;
            // m x (m x sigma)
            double mCross2X = _my * mCrossSigZ - _mz * mCrossSigY;
            double mCross2Y = _mz * mCrossSigX - _mx * mCrossSigZ;
            double mCross2Z = _mx * mCrossSigY - _my * mCrossSigX;

            double tSotX = -Gamma * tauDl * mCross2X;
            double tSotY = -Gamma * tauDl * mCross2Y;
            double tSotZ = -Gamma * tauDl * mCross2Z;

            // 길버트 댐핑 근사 반영
            double dmX = (tPrecX + tSotX + GilbertAlpha * (_my * tPrecZ - _mz * tPrecY)) * dtSec;
            double dmY = (tPrecY + tSotY + GilbertAlpha * (_mz * tPrecX - _mx * tPrecZ)) * dtSec;
            double dmZ = (tPrecZ + tSotZ + GilbertAlpha * (_mx * tPrecY - _my * tPrecX)) * dtSec;

            _mx += dmX;
            _my += dmY;
            _mz += dmZ;

            // 정규화 |m| = 1
            double norm = Math.Sqrt(_mx * _mx + _my * _my + _mz * _mz);
            _mx /= norm;
            _my /= norm;
            _mz /= norm;

            if (!switched && _mz < 0.0)
            {
                switched = true;
                switchedTimePs = step * 0.2;
            }
        }

        // 4. 열적 안정성 지수 Delta = (Ku * V) / (kB * T)
        double volume = Math.PI * 25e-9 * 25e-9 * ThicknessF;
        double ku = 0.5 * Mu0 * Ms * Hк;
        double deltaThermal = (ku * volume) / (1.38e-23 * 300.0);

        return new SpinSwitchingState(_mz, switchedTimePs, jRatio, deltaThermal, switched && _mz < -0.7);
    }
}