import matplotlib.pyplot as plt

%matplotlib inline
import numpy as np

from maxwellbloch import plot

Getting Started with the Two-Level Atom

The first model we’ll consider is a medium of two-level atoms, coupled by a weak pulse.

____ |1>
 |
 | Ω
 |
____ |0>

We’ll need to define a field, and the two-level atom for the field to address.

Defining the Field and Atom

from maxwellbloch import field, ob_atom

First we’ll define the field. The quantity we care about is the input Rabi frequency as a function of time, \(\Omega(t)\). Let’s say our input pulse has a max \(\Omega_0 = 2\pi \cdot 0.001 \mathrm{~ MHz}\), and is a Gaussian with a full-width at half-maximum (FWHM) of \(1 \mathrm{~ \mu s}\).

the_field = field.Field(
    coupled_levels=[[0, 1]],
    rabi_freq=1.0e-3,  # [2π MHz]
    rabi_freq_t_func="gaussian",
    rabi_freq_t_args={
        "ampl": 1.0,
        "centre": 0.0,  # [μs]
        "fwhm": 1.0,
    },  # [μs]
)

The first property coupled_levels is a list of pairs to be coupled by the field. In this case we only have two levels, indexed 0 (the ground state) and 1 (the excited state). The properties rabi_freq, rabi_freq_t_func and rabi_freq_t_args together define the field profile entering the medium. What does the profile look like in this case?

tlist = np.linspace(-2, 10, 100)
fig = plot.field_profile(the_field, tlist)
fig.show(renderer='notebook_connected')

It’s a Gaussian-shaped pulse, centred at \(t=0\). Many more pre-defined profiles are defined in the t_funcs module, see [REF t_funcs doc] for details.

So that’s the field. Now we need to define an atom to address. Here we just need to state that num_states=2. The atom also takes a list of fields, in this case we’ve only got one, the_field, so we’ll add that.

the_atom = ob_atom.OBAtom(
    num_states=2,
)
the_atom.fields.append(the_field)
the_atom.build_operators()

Solving the Maxwell-Bloch Equations

Now we have everything required to solve the problem. We just need to define the space \(z\) and time \(t\) dimesions we want to solve the problem over. Let’s say that our medium is \(1 \mathrm{~ cm}\) long and has constant number density, for example a vapour cell.

from maxwellbloch import mb_solve

mbs = mb_solve.MBSolve(
    t_min=-2.0,  # [μs]
    t_max=10.0,  # [μs]
    t_steps=100,
    z_min=-0.5,  # [cm]
    z_max=1.5,  # [cm]
    z_steps=10,
    interaction_strengths=[0.1],  # [2π MHz /cm]
)
mbs.atom = the_atom

Here we’ve set z_min=-0.5 and z_max=1.5 so that we can see how the pulse travels before and after it is in the medium, i.e. as it travels in vacuum.

Interaction Strength

The interaction strength (\(Ng\)) for each field describes the strength of the coupling of the field to the medium. You can think of it as the amount of absorption per unit length of medium, so in this system it has units of \(\mathrm{2\pi~MHz~/cm}\).

Now we can solve the system.

Omegas_zt, states_zt = mbs.mbsolve()

100%|██████████| 10/10 [00:00<00:00, 100.47z/s, Ω_max=0.00524]

The solver indicates progress as it goes, in 10% increments by default.

Plotting and Analysing Results

The solver returns two objects. Omegas_zt contains the complex value of each field, at each \(z\)-point and each \(t\)-point. In this case we have only one field, \(50+1\) \(z\)-points (for 50 z_steps) and \(100+1\) \(t\)-points (for 100 t_steps).

Omegas_zt.shape

(1, 11, 101)

states_zt contains the density matrix representation of the states of the atoms at each \(z\)-point and each \(t\)-point. The density matrix \(\rho\) is of size \(n \times n\) where \(n\) is the number of states in the atom, in this case \(n=2\).

states_zt.shape

(11, 101, 2, 2)

Both the Rabi frequency \(\Omega\) and the density matrix \(\rho\) are complex-valued, so the elements of both Omegas_zt and states_zt are of np.complex datatype.

fig = plot.field_spacetime(mbs, show_z_bounds=(0.0, 1.0))
fig.show(renderer='notebook_connected')

The plot above shows \(|\Omega(z, t)|\) over the full simulation domain. Time \(t\) is on the \(x\)-axis and propagation distance \(z\) on the \(y\)-axis, so the field enters at the bottom. The dotted grey lines mark \(z=0\) and \(z=1\), the boundaries of the medium. The horizontal slice at \(z=0\) represents the input field pulse.

This result is presented in the speed-of-light reference frame: a pulse travelling at \(c\) moves vertically in this plot. See [REF reference frame] for how to shift to the lab frame.

Analysis

What happens in this simulation? We can see that the pulse is attenuated by the medium, in that the maximum of the peak leaving the medium at \(z=1~\mathrm{cm}\) is lower than the maximum peak entering at \(z=0~\mathrm{cm}\).

We also see that the pulse is slightly fast, such that the peak arrives at the back of the medium, \(z=1~\mathrm{cm}\), at a time \(t < 0~\mathrm{\mu s}\) in this speed-of-light reference frame. This is typical of propagation through a medium with a normal dispersion profile. Spectral Analysis. There’s a little bit of ringing after (to the right) of the pulse, seen in the lighter region.

Defining the System with JSON

We went through defining Field, Atom and MBSolve objects separately so I could introduce these concepts, but in practice it is rare that we will need to interact with these directly. The whole problem can be defined via a JSON file or string, and this is the approach used throughout the rest of the documentation.

Spontaneous Decay

To include spontaneous decay in the Lindblad master equation we have to add collapse operators. We add a list of decays into the atom, each with a list of channels describing the decays (ordered lower then upper) and a decay rate in the same units as the Rabi frequency (in this case \(2\pi \mathrm{~ MHz}\)). Here the decay rate is \(2\pi \cdot 1 \mathrm{~ MHz}\). We also increase interaction_strength to \(1.0\) so that absorption is clearly visible in the analysis plots below.

mbs_json = """
{
  "atom": {
    "decays": [
      {
        "channels": [[0, 1]],
        "rate": 1.0
      }
    ],
    "fields": [
      {
        "coupled_levels": [[0, 1]],
        "rabi_freq": 1.0e-3,
        "rabi_freq_t_args": {
          "ampl": 1.0,
          "centre": 0.0,
          "fwhm": 1.0
        },
        "rabi_freq_t_func": "gaussian"
      }
    ],
    "num_states": 2
  },
  "t_min": -2.0,
  "t_max": 10.0,
  "t_steps": 100,
  "z_min": -0.5,
  "z_max": 1.5,
  "z_steps": 20,
  "interaction_strengths": [
    1.0
  ]
}
"""

mbs = mb_solve.MBSolve().from_json_str(mbs_json)
Omegas_zt, states_zt = mbs.mbsolve()

100%|██████████| 20/20 [00:00<00:00, 86.61z/s, Ω_max=0.00289]

fig = plot.field_spacetime(mbs, show_z_bounds=(0.0, 1.0))
fig.show(renderer='notebook_connected')

# Animated spatial profile — use the slider or Play to watch the pulse propagate
fig = plot.field_z_profile_anim(mbs)
fig.show(renderer='notebook_connected')

# Pulse area vs z — should decrease as the probe is absorbed
fig = plot.pulse_area(mbs)
fig.show(renderer='notebook_connected')

# Input (z=z_min) vs output (z=z_max) pulse shape
fig = plot.field_envelope(mbs, z_indices=[0, -1])
fig.show(renderer='notebook_connected')

Using a Natural Unit System

It is possible to use any reciprocal angular momentum and time units (e.g. \(2\pi \mathrm{~ MHz}\) and \(\mathrm{~ \mu s}\) in the above example). The length unit is also arbitrary (e.g. \(\mathrm{cm}\)) as long as the interaction strength \(Ng\) is given in the same units (e.g. \(\mathrm{2\pi~MHz~/cm}\)).

For a two-level system we have a single natural linewidth, and so it is convenient to introduce a natural unit system, with frequencies in units of the natural linewidth \(\Gamma\), times in units of the reciprocal spontaneous decay lifetime \(\tau = 1/\Gamma\) and distances in units of the length of the medium \(L\).

By introducing this natural unit system we are able to reduce the number of parameters involved in the mathematical problem. For example, it becomes clear that increasing the length of the medium ten times is equivalent to raising the number density by the same scale, or by choosing a system with a suitably higher dipole moment.

For most examples in this documentation from now on, we will use this natural unit system.

Shifting to the Fixed Frame of Reference

We solve the problem in a frame of reference that moves with the speed of light across the medium. This means we can keep time in the system separate from distance, which makes the coupled equations easier to solve. We can also see all the important details of the solution in this frame of reference, as in the plots above. But if we want to look at how a field actually propagates in time, we need to shift back to a fixed (or laboratory) frame of reference. To do that we need to connect the time and space dimensions via the speed of light in the system. The utilities for this shift are in the fixed module.

from maxwellbloch import fixed

We define the speed of light based on the units in the system. For example, in the units above, \(c \approx 3 \times 10^4 ~ \mathrm{cm / \mu s }\). If I put that speed in to the simulation results above, you will see no difference in the output. It will be too fast, and appear that the pulse progresses through the medium instantanously. Instead we can put in a slow speed-of-light to get an idea of what the shift does, say \(10^{-5} c\). Something like looking at the pulse in slow motion.

speed_of_light = 0.3  # [cm /μs] THIS IS 10^5 SLOWER THAN C!

tlist_fixed_frame = fixed.t_list(mbs, speed_of_light)
field_fixed_frame = fixed.rabi_freq(mbs, 0, speed_of_light, part="abs")

fig = plot.field_spacetime(mbs, speed_of_light=speed_of_light, show_z_bounds=(0.0, 1.0))
fig.show(renderer='notebook_connected')

So we see that in this simulation with a slow speed-of-light, the pulse would arrive at the back of the medium around \(3.33 \mathrm{~ \mu s}\) after it hit the front of the medium. The shift to the fixed reference frame has the effect of skewing the results rightward as we move through the medium, representing progression through time as well as space.