Sodium D2 LGS

Simulate a sodium LGS beacon

Imports

from pylgs.all import *

System setup

Define an LGS system on the D_2 line with pump and repump beams, with all parameters fixed except for the light intensities and the magnetic-field strength:

lgs = LGSSystem(
    'NaD2_Repump', 
    {
        'EllipticityDegrees1': 45.,
        'PolarizationAngleDegrees1': 0,
        'DetuningHz1': 1.0832e9,
        'LaserWidthHz1': 10.0e6,
        'EllipticityDegrees2': 45.,
        'PolarizationAngleDegrees2': 0,
        'DetuningHz2': -6.268e8,
        'LaserWidthHz2': 10.0e6,
        'MagneticZenithDegrees': 45.,
        'MagneticAzimuthDegrees': 45.,
        'SDampingCollisionRatePerS': 4081.63,
        'BeamTransitRatePerS': 131.944,
        'VccRatePerS': 28571.,
        'TemperatureK': 185.,
        'RecoilParameter': 1
    }
)

The model contains equations for 374 density-matrix elements:

lgs.A_ind.range
Density matrix(374) 
<xarray.DataArray 'XarrayVectorSpace' (Density matrix: 374)> Size: 3kB
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
Coordinates:
  * Density matrix  (Density matrix) <U68 102kB 'ρ<sub>Re, (3S<sub>1/2</sub>,...

Define sample values for the variable parameters:

params = {'IntensitySI1': 5., 'IntensitySI2': 46., 'BFieldG': 0.5}

Steady state of a cw LGS

Evenly spaced velocity groups

Build a model for the steady state of a cw LGS with 200 evenly spaced velocity groups:

model = lgs.stationary_model(vg=200)

Find the steady-state return flux:

model.integrated_flux(params).array
<xarray.DataArray (Integral: 1, Transition: 1)> Size: 8B
array([[13624.43746615]])
Coordinates:
  * Integral    (Integral) int64 8B 1
  * Transition  (Transition) <U33 132B '3P<sub>3/2</sub>→3S<sub>1/2</sub>'
Attributes:
    long_name:  Flux
    units:      photon/atom/s/sr

Find the density matrix for each of the 200 velocity groups:

sol = model.solve(params)
sol
∅ ⛒ Atomic velocity(200) ⛒ Density matrix(374) 
<xarray.DataArray (Atomic velocity: 200, Density matrix: 374)> Size: 598kB
array([[ 4.39802225e-07,  1.36780608e-07, -1.36682759e-07, ...,
         4.11230968e-11,  2.14971110e-13,  3.90896036e-14],
       [ 5.25601399e-07,  1.63464515e-07, -1.63347172e-07, ...,
         4.96315102e-11,  2.62023741e-13,  4.73042111e-14],
       [ 6.26990263e-07,  1.94995681e-07, -1.94855332e-07, ...,
         5.97982711e-11,  3.18862485e-13,  5.71544167e-14],
       ...,
       [ 6.25137543e-07,  1.95191270e-07, -1.95024443e-07, ...,
        -6.26727641e-11,  3.50062625e-13,  6.07509715e-14],
       [ 5.24064828e-07,  1.63626936e-07, -1.63487640e-07, ...,
        -5.19924049e-11,  2.87389234e-13,  5.02281842e-14],
       [ 4.38543488e-07,  1.36920286e-07, -1.36804175e-07, ...,
        -4.30592416e-11,  2.35563258e-13,  4.14633341e-14]],
      shape=(200, 374))
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 2kB -2.985 -2.955 ... 2.955 2.985
  * Density matrix   (Density matrix) <U68 102kB 'ρ<sub>Re, (3S<sub>1/2</sub>...

The result is expressed as an XarrayVectorArray object representing 200 vectors each of size 374. The XarrayVectorArray is a wrapper for an xarray DataArray, which can be accessed via the XarrayVectorArray.array attribute:

sol.array
<xarray.DataArray (Atomic velocity: 200, Density matrix: 374)> Size: 598kB
array([[ 4.39802225e-07,  1.36780608e-07, -1.36682759e-07, ...,
         4.11230968e-11,  2.14971110e-13,  3.90896036e-14],
       [ 5.25601399e-07,  1.63464515e-07, -1.63347172e-07, ...,
         4.96315102e-11,  2.62023741e-13,  4.73042111e-14],
       [ 6.26990263e-07,  1.94995681e-07, -1.94855332e-07, ...,
         5.97982711e-11,  3.18862485e-13,  5.71544167e-14],
       ...,
       [ 6.25137543e-07,  1.95191270e-07, -1.95024443e-07, ...,
        -6.26727641e-11,  3.50062625e-13,  6.07509715e-14],
       [ 5.24064828e-07,  1.63626936e-07, -1.63487640e-07, ...,
        -5.19924049e-11,  2.87389234e-13,  5.02281842e-14],
       [ 4.38543488e-07,  1.36920286e-07, -1.36804175e-07, ...,
        -4.30592416e-11,  2.35563258e-13,  4.14633341e-14]],
      shape=(200, 374))
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 2kB -2.985 -2.955 ... 2.955 2.985
  * Density matrix   (Density matrix) <U68 102kB 'ρ<sub>Re, (3S<sub>1/2</sub>...

Use the density-matrix solution to find the velocity distribution of flux as an XarrayVectorArray:

flux = model.flux_distribution(sol)
flux
∅ ⛒ Atomic velocity(200) ⛒ Transition(1) 
<xarray.DataArray (Atomic velocity: 200, Transition: 1)> Size: 2kB
array([[1.45147242e-03],
       [1.75827346e-03],
       [2.12654245e-03],
       [2.56792042e-03],
       [3.09606914e-03],
       [3.72704816e-03],
       [4.47967265e-03],
       [5.37597113e-03],
       [6.44168446e-03],
       [7.70684136e-03],
       [9.20643170e-03],
       [1.09810853e-02],
       [1.30780162e-02],
       [1.55519210e-02],
       [1.84660850e-02],
       [2.18936189e-02],
       [2.59188032e-02],
       [3.06387153e-02],
       [3.61649436e-02],
       [4.26255713e-02],
...
       [4.28284385e-02],
       [3.63339722e-02],
       [3.07793129e-02],
       [2.60359423e-02],
       [2.19907904e-02],
       [1.85466171e-02],
       [1.56185621e-02],
       [1.31330824e-02],
       [1.10265241e-02],
       [9.24387305e-03],
       [7.73760305e-03],
       [6.46695876e-03],
       [5.39670670e-03],
       [4.49665980e-03],
       [3.74094398e-03],
       [3.10742552e-03],
       [2.57717897e-03],
       [2.13408281e-03],
       [1.76440606e-03],
       [1.45647896e-03]])
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 2kB -2.985 -2.955 ... 2.955 2.985
  * Transition       (Transition) <U33 132B '3P<sub>3/2</sub>→3S<sub>1/2</sub>'
Attributes:
    long_name:  Flux
    units:      photon/atom/s/sr

XarrayVectorArray objects have built-in visualization routines. Plot the velocity distribution of return flux:

flux.visualize(xaxis_range=(-.5, .5), markers=True)

Plot the level population distribution:

model.level_population_distribution(sol).visualize()

Adaptive refinement of velocity groups

Build a steady-state model with velocity groups chosen to meet an error tolerance:

model, sol = lgs.adaptive_stationary_model(params, rtol=.02)

Find the total flux for this model:

model.total_flux(params).item()
13416.964517051674

Find the density matrix. In this case the adaptive refinement has produced 114 velocity groups:

sol = model.solve(params)
sol
∅ ⛒ Atomic velocity(53) ⛒ Density matrix(374) 
<xarray.DataArray (Atomic velocity: 53, Density matrix: 374)> Size: 159kB
array([[ 1.07885598e-05,  3.31085510e-06, -3.30868542e-06, ...,
         1.05793441e-09,  5.78622221e-12,  1.01169323e-12],
       [ 5.32761395e-05,  1.63488503e-05, -1.63377131e-05, ...,
         5.82227891e-09,  3.54933923e-11,  5.79328138e-12],
       [ 2.20288305e-04,  6.75972718e-05, -6.75486641e-05, ...,
         2.71860270e-08,  1.87180066e-10,  2.86116139e-11],
       ...,
       [ 2.19440999e-04,  6.76882852e-05, -6.76273924e-05, ...,
        -2.88993891e-08,  2.11351096e-10,  3.13976216e-11],
       [ 5.30971002e-05,  1.63685018e-05, -1.63547795e-05, ...,
        -6.14461243e-09,  3.95063446e-11,  6.25572371e-12],
       [ 1.07563362e-05,  3.31444482e-06, -3.31181135e-06, ...,
        -1.11016371e-09,  6.36805509e-12,  1.07875489e-12]],
      shape=(53, 374))
Coordinates:
  * Atomic velocity  (Atomic velocity) float64 424B -2.85 -2.55 ... 2.55 2.85
  * Density matrix   (Density matrix) <U68 102kB 'ρ<sub>Re, (3S<sub>1/2</sub>...

The flux distribution shows that this is a more detailed description of the system than the model with a greater number of evenly spaced velocity groups:

model.flux_distribution(sol).visualize(xaxis_range=(-.5, .5), markers=True)

Intensity dependence

Find the steady-state density matrix for 10 different values of repump intensity, with the total pump + repump intensity held constant:

sol = model.solve([
    {'BFieldG': .5},
    {
        'IntensitySI1': np.linspace(0, 25, 10), 
        'IntensitySI2': 50 - np.linspace(0, 25, 10)
    }
])
sol
IntensitySI1(10) ⛒ Atomic velocity(53) ⛒ Density matrix(374) 
<xarray.DataArray (IntensitySI1: 10, Atomic velocity: 53, Density matrix: 374)> Size: 2MB
array([[[ 6.91869921e-06,  1.85660353e-06, -1.85542136e-06, ...,
         -4.31560752e-19, -2.42743834e-21,  1.00579811e-18],
        [ 3.41640958e-05,  9.16840580e-06, -9.16227450e-06, ...,
         -1.26254786e-17, -7.74196638e-20,  6.30558084e-18],
        [ 1.41252139e-04,  3.79119525e-05, -3.78848282e-05, ...,
         -4.60780331e-17, -3.19632873e-19,  3.39669637e-17],
        ...,
        [ 1.40733857e-04,  3.79780524e-05, -3.79433122e-05, ...,
          3.17725806e-16, -2.32148657e-18,  3.85240063e-17],
        [ 3.40545912e-05,  9.18255516e-06, -9.17482708e-06, ...,
         -7.51400881e-18,  4.86643960e-20,  7.13977227e-18],
        [ 6.89900639e-06,  1.85917605e-06, -1.85770812e-06, ...,
          5.22170306e-20, -2.40337675e-22,  1.12145701e-18]],

       [[ 9.74992893e-06,  2.93087677e-06, -2.92883547e-06, ...,
          8.76577272e-10,  4.79429808e-12,  6.21223940e-13],
        [ 4.81461956e-05,  1.44725263e-05, -1.44620587e-05, ...,
          4.82430750e-09,  2.94095088e-11,  3.56022430e-12],
        [ 1.99070910e-04,  5.98392744e-05, -5.97936447e-05, ...,
          2.25269839e-08,  1.55101131e-10,  1.75980139e-11],
...
         -3.99102815e-08,  2.91875629e-10,  9.12249120e-11],
        [ 7.45793442e-05,  2.40872613e-05, -2.40790705e-05, ...,
         -8.48715660e-09,  5.45672603e-11,  1.81700774e-11],
        [ 1.51055158e-05,  4.87756752e-06, -4.87604167e-06, ...,
         -1.53355360e-09,  8.79663587e-12,  3.13206937e-12]],

       [[ 1.56167934e-05,  5.04354884e-06, -5.04260514e-06, ...,
          1.46848535e-09,  8.03170427e-12,  3.11251592e-12],
        [ 7.71218630e-05,  2.49051511e-05, -2.49001093e-05, ...,
          8.07883294e-09,  4.92498631e-11,  1.78405754e-11],
        [ 3.18907056e-04,  1.02976505e-04, -1.02953349e-04, ...,
          3.77041581e-08,  2.59600308e-10,  8.81869223e-11],
        ...,
        [ 3.18208080e-04,  1.03055867e-04, -1.03022832e-04, ...,
         -4.01790575e-08,  2.93841371e-10,  9.70125115e-11],
        [ 7.69741554e-05,  2.49223063e-05, -2.49151843e-05, ...,
         -8.54467210e-09,  5.49370684e-11,  1.93064452e-11],
        [ 1.55902304e-05,  5.04669106e-06, -5.04537392e-06, ...,
         -1.54398882e-09,  8.85649618e-12,  3.32510498e-12]]],
      shape=(10, 53, 374))
Coordinates:
  * IntensitySI1     (IntensitySI1) float64 80B 0.0 2.778 5.556 ... 22.22 25.0
  * Atomic velocity  (Atomic velocity) float64 424B -2.85 -2.55 ... 2.55 2.85
  * Density matrix   (Density matrix) <U68 102kB 'ρ<sub>Re, (3S<sub>1/2</sub>...

Plot the total flux to show the optimal value of repump intensity for these conditions:

model.integrated_flux(sol).visualize()

Create an animation of the flux distribution as the repump intensity is varied:

model.flux_distribution(sol).visualize(xaxis_range=(-.5, .5))

Animate the level population distribution as a function of repump intensity:

model.level_population_distribution(sol).visualize()

Plot return flux as pump and repump intensities are varied together:

intensity = 10**np.linspace(1, 2, 5) 
sol = model.solve([
    {'BFieldG': .5},
    {
        'IntensitySI1': 0.25 * intensity, 
        'IntensitySI2': intensity
    }
])
sol
IntensitySI1(5) ⛒ Atomic velocity(53) ⛒ Density matrix(374) 
<xarray.DataArray (IntensitySI1: 5, Atomic velocity: 53, Density matrix: 374)> Size: 793kB
array([[[ 1.26767368e-05,  3.80361958e-06, -3.80348690e-06, ...,
          6.56480432e-10,  3.59050175e-12,  4.31815157e-13],
        [ 6.25925724e-05,  1.87805354e-05, -1.87797672e-05, ...,
          3.61269651e-09,  2.20233117e-11,  2.48215957e-12],
        [ 2.58769554e-04,  7.76417811e-05, -7.76379207e-05, ...,
          1.68674832e-08,  1.16134255e-10,  1.23062155e-11],
        ...,
        [ 2.58559617e-04,  7.76696127e-05, -7.76627208e-05, ...,
         -1.79225872e-08,  1.31075297e-10,  1.35240189e-11],
        [ 6.25481810e-05,  1.87865008e-05, -1.87850950e-05, ...,
         -3.81132525e-09,  2.45048562e-11,  2.68442346e-12],
        [ 1.26687505e-05,  3.80470572e-06, -3.80445867e-06, ...,
         -6.88677517e-10,  3.95038187e-12,  4.61147998e-13]],

       [[ 1.27177366e-05,  3.93995887e-06, -3.93953914e-06, ...,
          8.45735826e-10,  4.62562799e-12,  7.53356011e-13],
        [ 6.27973126e-05,  1.94540766e-05, -1.94518044e-05, ...,
          4.65417659e-09,  2.83724550e-11,  4.32110109e-12],
        [ 2.59628351e-04,  8.04285755e-05, -8.04179718e-05, ...,
          2.17300088e-08,  1.49614347e-10,  2.13765757e-11],
...
         -4.19345628e-08,  3.06680218e-10,  7.68785388e-11],
        [ 6.06271301e-05,  1.91444129e-05, -1.91231452e-05, ...,
         -8.91585110e-09,  5.73235698e-11,  1.53376809e-11],
        [ 1.22816666e-05,  3.87627337e-06, -3.87217500e-06, ...,
         -1.61079923e-09,  9.23973401e-12,  2.64828430e-12]],

       [[ 1.20015785e-05,  3.79449448e-06, -3.78626290e-06, ...,
          2.15886209e-09,  1.18075639e-11,  4.60530854e-12],
        [ 5.92869948e-05,  1.87413154e-05, -1.86996192e-05, ...,
          1.18796015e-08,  7.24194759e-11,  2.63778243e-11],
        [ 2.45256152e-04,  7.75149680e-05, -7.73362063e-05, ...,
          5.54599155e-08,  3.81849299e-10,  1.30301317e-10],
        ...,
        [ 2.43092318e-04,  7.76384816e-05, -7.74285779e-05, ...,
         -5.91161588e-08,  4.32335492e-10,  1.43314085e-10],
        [ 5.88298834e-05,  1.87686758e-05, -1.87204344e-05, ...,
         -1.25667552e-08,  8.07968508e-11,  2.85370847e-11],
        [ 1.19193863e-05,  3.79959331e-06, -3.79018826e-06, ...,
         -2.27011763e-09,  1.30216832e-11,  4.91824662e-12]]],
      shape=(5, 53, 374))
Coordinates:
  * IntensitySI1     (IntensitySI1) float64 40B 2.5 4.446 7.906 14.06 25.0
  * Atomic velocity  (Atomic velocity) float64 424B -2.85 -2.55 ... 2.55 2.85
  * Density matrix   (Density matrix) <U68 102kB 'ρ<sub>Re, (3S<sub>1/2</sub>...
model.integrated_flux(sol).visualize()

Detuning dependence

lgs = LGSSystem(
    'NaD2_Repump', 
    {
        'EllipticityDegrees1': 45.,
        'PolarizationAngleDegrees1': 0,
        'IntensitySI1': 5.,
        'DetuningHz1': 1.0832e9,
        'LaserWidthHz1': 10.0e6,
        'EllipticityDegrees2': 45.,
        'PolarizationAngleDegrees2': 0,
        'IntensitySI2': 46.,
        # 'DetuningHz2': -6.268e8,
        'LaserWidthHz2': 10.0e6,
        'BFieldG': 0.5,
        'MagneticZenithDegrees': 45.,
        'MagneticAzimuthDegrees': 45.,
        'SDampingCollisionRatePerS': 4081.63,
        'BeamTransitRatePerS': 131.944,
        'VccRatePerS': 28571.,
        'TemperatureK': 185.,
        'RecoilParameter': 1
    }
)
results = []
deltas = np.linspace(-2e9, 1e9, 100)
for delta in deltas:
    model, sol = lgs.adaptive_stationary_model({'DetuningHz2': delta}, rtol=.005)
    results.append(model.integrated_flux(sol).array)
results = xr.concat(results, dim='delta').assign_coords(delta=deltas).squeeze()
plotly_dataarray(results, markers=True)

Transient dynamics

Transient dynamics from a laser pulse (Adams solver)

Build a dynamical model (InstationaryModel) using the velocity groups from the adaptively refined steady-state model:

lgs = LGSSystem(
    'NaD2_Repump', 
    {
        'EllipticityDegrees1': 45.,
        'PolarizationAngleDegrees1': 0,
        'DetuningHz1': 1.0832e9,
        'LaserWidthHz1': 10.0e6,
        'EllipticityDegrees2': 45.,
        'PolarizationAngleDegrees2': 0,
        'DetuningHz2': -6.268e8,
        'LaserWidthHz2': 10.0e6,
        'MagneticZenithDegrees': 45.,
        'MagneticAzimuthDegrees': 45.,
        'SDampingCollisionRatePerS': 4081.63,
        'BeamTransitRatePerS': 131.944,
        'VccRatePerS': 28571.,
        'TemperatureK': 185.,
        'RecoilParameter': 1,
        'BFieldG': .5
    }
)
model = lgs.instationary_model(
    vg=100, 
    T=1e-8, 
    num_values=100,
    time_stepper=AdamsTimeStepper()
)
params = {
    'IntensitySI1': "5.*sin(1e6*t)**2", 
    'IntensitySI2': "46.*sin(1e6*t)**2",
}

Solve for the density-matrix evolution. The solution contains the density matrix for each velocity group at 100 time values:

sol = model.solve(params)

Plot the total flux as a function of time:

model.integrated_flux(sol).visualize()

Animate the flux distribution as a function of time:

model.flux_distribution(sol).visualize()

Animate the level population distribution:

model.level_population_distribution(sol).visualize(xaxis_range=(-1, 1))
model.total_population(sol).visualize()

Transient dynamics from a laser pulse

Build a dynamical model (InstationaryModel) using the velocity groups from the adaptively refined steady-state model:

model = lgs.instationary_model(
    vg=model.data['velocity_groups'], 
    T=1e-5, 
    num_values=100
)

Set tolerances for the BDF solver:

pymor.basic.set_defaults({
    'pylgs.pymor.timestepping.cvode_solver_options.cvode_bdf_atol': 1e-6,
    'pylgs.pymor.timestepping.cvode_solver_options.cvode_bdf_rtol': 1e-4,
})

Solve for the density-matrix evolution. The solution contains the density matrix for each velocity group at 100 time values:

sol = model.solve(params)
sol
Time(57) ⛒ Atomic velocity(100) ⛒ Density matrix(374) 
<xarray.DataArray (Time: 57, Atomic velocity: 100, Density matrix: 374)> Size: 17MB
array([[[ 6.27745202e-07,  0.00000000e+00,  0.00000000e+00, ...,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
        [ 8.93121976e-07,  0.00000000e+00,  0.00000000e+00, ...,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
        [ 1.26157506e-06,  0.00000000e+00,  0.00000000e+00, ...,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
        ...,
        [ 1.26157506e-06,  0.00000000e+00,  0.00000000e+00, ...,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
        [ 8.93121976e-07,  0.00000000e+00,  0.00000000e+00, ...,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
        [ 6.27745202e-07,  0.00000000e+00,  0.00000000e+00, ...,
          0.00000000e+00,  0.00000000e+00,  0.00000000e+00]],

       [[ 6.27745197e-07,  6.64372848e-19,  6.37767192e-19, ...,
          2.24369430e-12,  2.76625963e-14,  0.00000000e+00],
        [ 8.93121969e-07,  9.84131179e-19,  9.44981078e-19, ...,
          3.25650041e-12,  4.09584164e-14,  0.00000000e+00],
        [ 1.26157505e-06,  1.44853373e-18,  1.39130928e-18, ...,
          4.69450367e-12,  6.02588426e-14,  0.00000000e+00],
...
         -1.73853390e-10,  1.02500680e-12,  6.44068195e-14],
        [ 9.00328894e-07,  2.54194536e-09, -1.21382873e-08, ...,
         -1.20478847e-10,  6.95323830e-13,  4.42363750e-14],
        [ 6.32807969e-07,  1.78677591e-09, -8.52665523e-09, ...,
         -8.29291448e-11,  4.68716078e-13,  3.01950954e-14]],

       [[ 6.35112431e-07,  2.02660953e-09, -8.24916037e-09, ...,
          1.04831075e-10,  5.64363891e-13,  4.67700957e-14],
        [ 9.03611931e-07,  2.88509730e-09, -1.17422919e-08, ...,
          1.52153533e-10,  8.35629924e-13,  6.85238154e-14],
        [ 1.27640258e-06,  4.07779555e-09, -1.65946701e-08, ...,
          2.19343763e-10,  1.22940878e-12,  9.97668830e-14],
        ...,
        [ 1.27639115e-06,  3.93735934e-09, -1.66262198e-08, ...,
         -2.30089689e-10,  1.35278523e-12,  1.07248167e-13],
        [ 9.03604161e-07,  2.78777468e-09, -1.17638635e-08, ...,
         -1.59450177e-10,  9.17676074e-13,  7.34989180e-14],
        [ 6.35108311e-07,  1.95968385e-09, -8.26410513e-09, ...,
         -1.09754311e-10,  6.18603432e-13,  5.00590650e-14]]],
      shape=(57, 100, 374))
Coordinates:
  * Time             (Time) float64 456B 0.0 1.193e-08 ... 9.855e-06 1.001e-05
  * Atomic velocity  (Atomic velocity) float64 800B -2.97 -2.91 ... 2.91 2.97
  * Density matrix   (Density matrix) <U68 102kB 'ρ<sub>Re, (3S<sub>1/2</sub>...

Plot the total flux as a function of time:

model.total_flux(sol).visualize()

Animate the flux distribution as a function of time:

model.flux_distribution(sol).visualize()

Animate the level population distribution:

model.level_population_distribution(sol).visualize(xaxis_range=(-1, 1))
model.total_population(sol).visualize()

Transient dynamics with modulated laser intensity

Solve the InstationaryModel with sinusoidally modulated light intensity:

sol = model.solve({
    'IntensitySI1': "5.*sin(1e6*t)**2", 
    'IntensitySI2': "46.*sin(1e6*t)**2",
})

The total flux as a function of time:

model.total_flux(sol).visualize()

The animated flux distribution:

model.flux_distribution(sol).visualize(xaxis_range=(-.25, .25))

The animated level population distribution:

model.level_population_distribution(sol).visualize(xaxis_range=(-.25, .25))

Periodic state of a modulated system

The StationaryFloquetModel solves for the periodic state of a modulated system after the transient dynamics have died out.

Build a stationary Floquet model for six Fourier harmonics of a system with sinusoidally modulated intensity:

model = lgs.stationary_floquet_model(
    {
        'IntensitySI1': "5.*sin(nu*t)**2", 
        'IntensitySI2': "46.*sin(nu*t)**2"
    },
    time=np.linspace(0, 1e-5, 100),
    n_vector=6,
    vg=model.data['velocity_groups']
)

Define values for the modulation frequency and magnetic field:

mu = {'nu': 1e6}

Solve for the periodic state. The solution contains 13 Fourier coefficients for each density-matrix element in each velocity group:

sol = model.solve(mu)
sol
∅ ⛒ Atomic velocity(100) ⛒ Fourier coefficient(13) ⛒ Density matrix(374) 
<xarray.DataArray (Atomic velocity: 100, Fourier coefficient: 13,
                   Density matrix: 374)> Size: 8MB
array([[[-3.33367347e-12-2.47873923e-12j,
         -5.53416897e-12+2.76553604e-12j,
          6.95526475e-13+1.65771275e-13j, ...,
          1.91710895e-15-8.44773128e-16j,
          1.03835666e-17-3.64382081e-18j,
          1.47643825e-15-1.28216448e-16j],
        [ 0.00000000e+00+0.00000000e+00j,
          0.00000000e+00+0.00000000e+00j,
          0.00000000e+00+0.00000000e+00j, ...,
          0.00000000e+00+0.00000000e+00j,
          0.00000000e+00+0.00000000e+00j,
          0.00000000e+00+0.00000000e+00j],
        [-7.54791468e-11+9.56668441e-11j,
          1.58595587e-10+4.15846078e-11j,
          3.17070272e-11+1.37616560e-10j, ...,
         -4.11980479e-12-4.30059377e-15j,
         -2.16407808e-14-1.41428212e-15j,
         -8.37021058e-16-3.02715067e-16j],
        ...,
        [-7.54764350e-11-9.56822983e-11j,
...
         -8.78795776e-16-3.03163727e-16j],
        ...,
        [-8.57213246e-11-8.81627867e-11j,
          1.52781103e-10-5.54094578e-11j,
          2.40109399e-11-1.42535875e-10j, ...,
          4.30144959e-12-6.33183866e-15j,
         -2.36515217e-14+1.55683335e-15j,
         -8.78787370e-16+3.03149094e-16j],
        [ 0.00000000e+00+0.00000000e+00j,
          0.00000000e+00+0.00000000e+00j,
          0.00000000e+00+0.00000000e+00j, ...,
          0.00000000e+00+0.00000000e+00j,
          0.00000000e+00+0.00000000e+00j,
          0.00000000e+00+0.00000000e+00j],
        [ 3.69219016e-12+2.17471410e-12j,
         -5.94433688e-12+7.04489108e-12j,
          1.85848150e-12+9.20116130e-12j, ...,
         -4.33573739e-15-1.25641127e-15j,
          2.43749224e-17+5.04146035e-18j,
          1.57807894e-15+1.33303376e-16j]]], shape=(100, 13, 374))
Coordinates:
  * Atomic velocity      (Atomic velocity) float64 800B -2.97 -2.91 ... 2.97
  * Fourier coefficient  (Fourier coefficient) int64 104B -6 -5 -4 -3 ... 4 5 6
  * Density matrix       (Density matrix) <U68 102kB 'ρ<sub>Re, (3S<sub>1/2</...

The real and imaginary parts of the Fourier coefficients of the time-dependent return flux:

model.total_flux(sol).visualize()

The total return flux reconstructed as a function of time:

model.reconstructed_total_flux(sol, mu).real.visualize()

The level population distribution as a function of time:

model.level_population_distribution(model.reconstruct(sol, mu)).real.visualize()