Add muon-catalyzed fusion physics

[stognini]

Physical processes

Necessary for a proof-of-principle

• Muonic atom capture (G4/NK)
  • Formation of $d_\mu$ and $t_\mu$
• Muonic molecule formation (G4/NK)
  • Formation of $(dd)_\mu$ , $(dt)_\mu$ , or $(tt)_\mu$
• MuCF (NK)
  • Fusion of of $(dd)_\mu$ , $(dt)_\mu$ , and $(tt)_\mu$
• Muon decay (G4)
• Muon basic EM physics
• Neutron, alpha, proton, $^3\text{He}$, and $^4\text{H}$ transport

Good to have

• Muonic atom decay (G4)
• Isotopic transfer (NK)
  • Muon tends to be transferred to higher Z atoms
• Muon atom stripping (NK)
• Alpha deexcitation cascade
• Spin flip
• Muon, alpha, neutron (G4)
  • Energy-loss
  • Ionization
  • Scattering

[stognini]

Physics cycle

• Physical interaction times

flowchart LR
mu[Muon] --"10^{-12} - 10^{-13} s"--> muatom["Atom formation"]
muatom --"10^{-8} - 10^{-10} s"--> molecule["Molecule formation"]
molecule --"10^{-12} s"--> Fusion

• Cycle times (atom formation to fusion) heavily depend on material temperature [PhysRevLett.51.1757]

  • Ref. above states the $(dt)_\mu$ molecular formation $\lambda_{dt\mu}$ being 3 to 7 $\times 10^{-8}$ s, but theoretical limit is at $10^{-10}$ s [10.1016/0920-3796(89)90023-9]

• The muon can decay at any time

• Fusion channels:

  • $(dd)_\mu$
    • $\longrightarrow ^3\text{He} + \mu + n + 3.27 \ \text{MeV}$
    • $\longrightarrow (^3\text{He})_\mu + n + 3.27 \ \text{MeV}$
    • $\longrightarrow t + \mu + p + 4.03 \ \text{MeV}$
    • $\longrightarrow (t)_\mu + p + 4.03 \ \text{MeV}$
  • $(dt)_\mu$
    • $\longrightarrow \alpha + \mu + n + 17.6 \ \text{MeV}$
    • $\longrightarrow (\alpha)_\mu + n + 17.6 \ \text{MeV}$
  • $(tt)_\mu$
    • $\longrightarrow \alpha + \mu + 2n + 11.33 \ \text{MeV}$
    • $\longrightarrow (\alpha)_\mu + 2n + 11.33 \ \text{MeV}$
    • There is also the resonant channel $^5\text{He} + \mu + n \longrightarrow t + \mu + 2n + 8.7 \ \text{MeV}$
      • Probably below second order effects in terms of relevance

• Experimental data is broken down into grid tables with temperature (x) vs. rate (y). This applies to mean cycle times, mean atom transfers, and mean atom spin flip.

• Cycle time depends on spin as well. Theoretically, there are many possible states (see table), but only a few are reactive, i.e. have the correct rotational and nuclear-spin symmetry to tunnel through the Coulomb barrier. All these project onto $J = 1$.

  • States not containing $J = 1$ have fusion rates orders of magnitude lower
  • $(dt)_\mu$ is the dominant channel by far (> 80%).

  Molecule	Nuclei	$I_N$	$F = I_N \pm 1/2$	Reactive states (F)
  $(dd)_\mu$	1, 1	0, 1, 2	1/2, 3/2, 5/2	1/2, 3/2
  $(dt)_\mu$	1, 1/2	1/2, 3/2	0, 1, 2	0, 1
  $(tt)_\mu$	1/2, 1/2	0, 1	1/2, 3/2	1/2

Only the reactive states are implemented

Implementation

Geant4

• In Geant4, the MuonCatalyzedFusion::AtRestDotIt sequence is

---
title: DD
---
flowchart LR
a["calculate cycle time"] --> b["sample sticking"] --> d["sample secondaries"]
d --> he3["He-3"]
d --> h["H-3"]
he3 --> he3s["(He-3)mu + n"]
he3 --> he3ns["He-3 + mu + n"]
h -. not implemented .-> hs["(H-3)mu + p"]
h --> hns["H-3 + mu + p"]

---
title: DT
---
flowchart LR
a["calculate cycle time"] --> b["sample sticking"] --> d["sample secondaries"]
d --> e
d --> f
e["(alpha)mu + n"]
f["alpha + mu + n"]

---
title: TT
---
flowchart LR
a["calculate cycle time"] --> b["sample sticking"] --> d["sample secondaries"]
d --> e
d --> f
e["(alpha)mu + n + n"]
f["alpha + mu + n + n"]

Celeritas

• In Celeritas, one large Executor does all the work
• At rest process:
  • Form muonic atom
  • May execute atom spin flip or atom transfer
  • Calculate mean cycle time (time it takes from atom formation to fusion)
  • Form muonic molecule and select its spin
  • Confirm if fusion happens or the if the muon decays before that
  • Select channel and call appropriate Interactor

flowchart LR

%% Main flowchart
atom[muonic atom]
molecule[muonic molecule]
spinflip[spin flip]
transfer[isotopic transfer]
cycle[cycle time]

atom --> molecule --> cycle --> dd & dt & tt
atom --> spinflip & transfer --> molecule
molecule -.-> decay

subgraph second[Second order effects]
    spinflip
    transfer
end

subgraph mucf[muCF]
    direction LR

    %% DD fusion
    dd[DD]
    he3[$$^3$$He]
    h[H]
    ddhe3s["(He-3)mu + n"]
    ddhe3ns["He-3 + mu + n"]
    ddhs["(H-3)mu + p"]
    ddhns["H-3 + mu + p"]
    dd --> he3 & h
    he3 --> ddhe3s & ddhe3ns
    h -.not
        implemented.-> ddhs
    h --> ddhns

    %% DT fusion
    dt[DT]
    dts["(alpha)mu + n"]
    dtns["alpha + mu + n"]
    dt --> dts & dtns

    %% TT fusion
    tt[TT]
    tts["(alpha)mu + n + n"]
    ttns["alpha + mu + n + n"]
    tt --> tts & ttns
end

ddhe3s & ddhs & dts & tts --> stripping["Muon stripping"]

• Each Interactor (one for each: dd, dt, tt) takes as input the selected fusion channel and it simply returns an Interaction::from_absorption() and samples the outgoing secondaries.

Model/Executor class diagram

classDiagram
class MuonCatalyzedFusionExecutor {
    At-rest muon-catalyzed fusion executor
    Interaction operator()(celeritas::CoreTrackView const& track)
}

MuonCatalyzedFusionExecutor <-- DTMixMuonicAtomSelector
MuonCatalyzedFusionExecutor <-- DTMixMuonicAtomSpinSelector
MuonCatalyzedFusionExecutor <-- DTMixMuonicMoleculeSelector
MuonCatalyzedFusionExecutor <-- DTMixMuonicMoleculeSpinSelector
MuonCatalyzedFusionExecutor <-- DDMuonCatalyzedFusionInteractor
MuonCatalyzedFusionExecutor <-- DTMuonCatalyzedFusionInteractor
MuonCatalyzedFusionExecutor <-- TTMuonCatalyzedFusionInteractor

DDMuonCatalyzedFusionInteractor <-- DDChannelSelector
DTMuonCatalyzedFusionInteractor <-- DTChannelSelector
TTMuonCatalyzedFusionInteractor <-- TTChannelSelector

DTMixMuonCatalyzedFusionModel <-- MuonCatalyzedFusionExecutor
DTMixMuonCatalyzedFusionModel <-- DTMixMuonCatalyzedFusionData

DTMixMuonCatalyzedFusionData <-- LiquidHydrogenDensityCalculator
DTMixMuonCatalyzedFusionData <-- EquilibrateHydrogenDensityCalculator
DTMixMuonCatalyzedFusionData <-- MeanCycleTimeCalculator


%%---------------------------------------------------------------------------%%
%% Purpose-specific helper classes
namespace detail {
    class DTMixMuonicAtomSelector {
        Muonic atom formation
        DTMixMuonicAtomSelector(...)
        Result operator(Engine& rng)
    }
    class DTMixMuonicAtomSpinSelector {
        Atom spin selector
        DTMixMuonicAtomSpinSelector(...)
        Result operator(Engine& rng)
    }
    class DTMixMuonicMoleculeSelector {
        Muonic molecule formation
        DTMixMuonicMoleculeSelector(...)
        Result operator(Engine& rng)
    }
    class DTMixMuonicMoleculeSpinSelector {
        Molecule spin selector
        DTMixMuonicMoleculeSpinSelector(...)
        Result operator(Engine& rng)
    }
    class DDChannelSelector {
        Select channel for dd fusion
        DDChannelSelector(...)
        Result operator(Engine& rng)
    }
    class DTChannelSelector {
        Select channel for dt fusion
        DTChannelSelector(...)
        Result operator(Engine& rng)
    }
    class TTChannelSelector {
        Select channel for tt fusion
        TTChannelSelector(...)
        Result operator(Engine& rng)
    }
}

%%---------------------------------------------------------------------------%%
%% Return Interaction result
namespace Interaction {
    class DDMuonCatalyzedFusionInteractor {
        Interactor for DD muCF
        DDMuonCatalyzedFusionInteractor(...)
        Result operator(Engine& rng)
    }
    class DTMuonCatalyzedFusionInteractor {
        Interactor for DT muCF
        DTMuonCatalyzedFusionInteractor(...)
        Result operator(Engine& rng)
    }
    class TTMuonCatalyzedFusionInteractor {
        Interactor for TT muCF
        DTMuonCatalyzedFusionInteractor(...)
        Result operator(Engine& rng)
    }
}

%%---------------------------------------------------------------------------%%
%% Process / Model
namespace Model {
    class DTMixMuonCatalyzedFusionModel {
        DT muCF model
        DTMixMuonCatalyzedFusionModel()
        void step(CoreParams const&, CoreState&)
        DTMixMuonCatalyzedFusionData data_
    }
    class DTMixMuonCatalyzedFusionData {
        Static scalars
        Static grid tables
        Material-dependent data
    }
    class EquilibrateHydrogenDensityCalculator {
        Ensure hydrogen mixture is at equilibrium
        EquilibrateHydrogenDensityCalculator(...)
        Result operator()
    }
    class LiquidHydrogenDensityCalculator {
        Calculate dt fractions in units of LHD
        LiquidHydrogenDensityCalculator(...)
        Result operator()
    }
    class MeanCycleTimeCalculator {
        Calculate fusion cycle time
        MeanCycleTimeCalculator(...)
        Result operator()
    }
}

Model/Executor host/device data

flowchart LR
subgraph memspace["DTMixMuonCatalyzedFusionData"]
    subgraph static scalars
        pids[Particle Ids
             - Elementary
             - Nuclei
             - Muonic atoms]
    end
    subgraph per-material scalars
        dtfractions["DT LHD fractions"]
        dtequilibrium["DT equilibrium fractions"]
        cycle["Cycle time per molecule"]
        spinflip["Spin flip per combination"]
        transfer["Atom transfer per combination"]
    end
    subgraph static tables
        energy["Muon energy CDF"]
    end
end

MuCF in the stepping loop

flowchart LR
rest[At rest]
executor[MuCF Executor]
mucf[MuCF Interactor]
decay[decay Interactor]

muon --EM physics--> rest --> executor --> mucf
muon -.-> decay
executor -.-> decay