import numpy as np

Spectral Analysis

We can study the response of a quantised system to an input field using spectral analysis. For linear systems, analytic response functions are known. We can perform spectral analysis by transforming into the frequency domain.

Here we’ll look at how to use the spectral module to obtain the absorption and dispersion profiles by transforming into the frequency domain using a discrete Fourier transform.

First we’ll define a typical two-level problem. We’ll need to use a high resolution in the time domain if we want to consider high frequencies above and below resonance. So here, note that we use a lot of time steps. We wouldn’t usually want this! It will take a long time to save and result in a large savefile.

mb_solve_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": 1000,
  "z_min": -0.2,
  "z_max": 1.2,
  "z_steps": 50,
  "interaction_strengths": [
    1.0
  ],
  "savefile": "spectral-analysis"
}
"""

from maxwellbloch import mb_solve

mbs = mb_solve.MBSolve().from_json_str(mb_solve_json)

Omegas_zt, states_zt = mbs.mbsolve(recalc=False)

/home/docs/checkouts/readthedocs.org/user_builds/maxwellbloch/envs/v0.11.0/lib/python3.11/site-packages/maxwellbloch/mb_solve.py:344: UserWarning: Savefile was built with maxwellbloch==0.10.0, current version is 0.11.0.
  self.load_results()

Results in the Time Domain

First we’ll look at the results in the time domain as usual.

from maxwellbloch.plot import fields

fields.field_spacetime(mbs).show()

Data type cannot be displayed: application/vnd.plotly.v1+json

As we can see, it’s the typical absorption picture we’ve seen before.

Results in the Frequency Domain

The envelope of a short Gaussian pulse on resonance is equivalent to a superposition of waves across a wide range of frequencies. We can analyse the spectral response of the medium to the field by transforming the simulated Rabi frequency in the time domain, \(\Omega(z, t)\), to the frequency domain, \(\Omega(z, \omega)\).

We determine the absoprtion profile by measuring the intensity of the pulse that arrives at the back of the medium relative to the intensity put in at the front. This we can relate to the imaginary part of the optical susceptibility \(\chi_I\) via

\[\frac{k}{2} \chi_I(\omega) z = -\log \frac{| \Omega (z, \omega) |}{| \Omega(z = 0, \omega) |}\]

where \(k\) is the wavenumber.

We determine the dispersion profile by measuring the phase shift of the pulse that arrives at the back of the medium relative to the pulse that arrived at the front. This we can relate to the real part of optical susceptibility \(\chi_R\) via

\[\frac{k}{2} \chi_R(\omega) z = \phi(z, \omega) - \phi(z=0, \omega)\]

where \(\phi\) is the complex phase defined via

\[\Omega(z, \omega) = |\Omega|\mathrm{e}^{\mathrm{i}\phi}.\]

See section 2.5 of my thesis for an explanation of how the susceptibility is derived.

Methods to provide the dispersion and absorption of a given field are provided in the spectral module.

from maxwellbloch import spectral, utility
from maxwellbloch.plot import spectra
import plotly.graph_objects as go

freq_list = spectral.freq_list(mbs)

fig = spectra.spectrum(mbs, show_dispersion=True, freq_range=3.0)

# Overlay the normalised input pulse
rabi_freq_abs_0 = np.abs(spectral.rabi_freq(mbs, 0))[0]
mask = np.abs(freq_list) <= 3.0
fig.add_trace(go.Scatter(
    x=freq_list[mask],
    y=(rabi_freq_abs_0 / np.max(rabi_freq_abs_0))[mask],
    mode="lines",
    line=dict(dash="dash"),
    name="pulse (normalised)",
))
fig.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

The input pulse profile (normalised) is shown in grey. In blue we see the familiar Lorentzian absorption profile of homogeneous broadening due to the spontaneous decay. The full-width at half maximum (FWHM) of the peak is the decay rate, \(1 \Gamma\). The associated dispersion profile is shown in orange, describing the phase of the input field relative to the frequency of the oscillator as it passes over resonance from lagging to leading.

Comparing with Analytic Results

For a weak input field like this one, the field is in the regime of linear optics and it is possible to use the weak probe approximation to derive an analytic lineshape

\[\frac{k}{2} \chi({\omega}) = -Ng \frac{1}{\mathrm{i}\frac{\Gamma}{2} + \omega}.\]

For details on how this is derived, see equation (2.61) of my thesis.

These analytic absorption and dispersion profiles for weak fields are also included in the spectral module, so you can check if you are really in the linear regime.

interaction_strength = mbs.interaction_strengths[0]
decay_rate = mbs.atom.decays[0]["rate"]
absorption_linear_known = spectral.absorption_two_linear_known(
    freq_list, interaction_strength, decay_rate
)
dispersion_linear_known = spectral.dispersion_two_linear_known(
    freq_list, interaction_strength, decay_rate
)

fig = spectra.spectrum(mbs, show_dispersion=True, freq_range=3.0)

mask = np.abs(freq_list) <= 3.0
fig.add_trace(go.Scatter(
    x=freq_list[mask],
    y=absorption_linear_known[mask],
    mode="lines",
    line=dict(dash="dash", color="white"),
    name="analytic absorption",
))
fig.add_trace(go.Scatter(
    x=freq_list[mask],
    y=dispersion_linear_known[mask],
    mode="lines",
    line=dict(dash="dash", color="white"),
    yaxis="y2",
    name="analytic dispersion",
))

hm, r1, r2 = utility.half_max_roots(freq_list, spectral.absorption(mbs, field_idx=0))
fig.add_shape(type="line", x0=r1, x1=r2, y0=hm, y1=hm, line=dict(dash="dot"))
fig.add_annotation(x=r2, y=hm, text=f"FWHM: {r2 - r1:.2f}", showarrow=False, xshift=40)
fig.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

fig = go.Figure(
    layout=go.Layout(
        title="Absorption residuals",
        xaxis=dict(title="frequency (γ)", range=[-3.0, 3.0]),
        yaxis=dict(title="residual (a.u.)", range=[-3e-2, 3e-2]),
        width=700,
        height=200,
        template="maxwellbloch",
    )
)
fig.add_trace(go.Scatter(
    x=freq_list,
    y=spectral.absorption(mbs, 0, -1) - absorption_linear_known,
    mode="lines",
    name="absorption",
))
fig.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

The known analytic lineshapes are dashed white over the simulated lineshapes. For absorption, we see a maximum absolute residual of 3% and a simualated full width at half maximum (FWHM) of \(1 \Gamma\), which is the decay rate. This tells us we are in the linear regime with negligible additional broadening. At higher field intensities we will see power broadening. Inclusion of effects like inhomogeneous broadening due to the thermal motion and collisions will also lead to broadened lineshapes.