🧰 Importing Required Libraries

🧰 Importing Required Libraries#

In this notebook, we’ll use standard Python packages for numerical calculations and plotting:

  • NumPy — numerical arrays and linear algebra

  • Matplotlib — plotting and visualization

  • Matplotlib.cm — access to built-in colormaps

Run the cell below to import all required modules and set default plot styles.

# === Import Required Libraries ===

import numpy as np
from numpy import linalg as la
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib import rcParams

# --- Plot Styling ---
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 12

print("✅ Libraries imported and plot settings configured.")

✅ Libraries imported and plot settings configured.

🌟 Overview#

In this notebook, we will explore vibrational strong coupling between a cavity and the HF molecule, where the cavity is tuned to the fundamental transition of HF within the harmonic oscillator approximation.

We will implement and compare two forms of the light–matter Hamiltonian:

  • the minimal coupling (or p·A) Hamiltonian, and

  • the Pauli–Fierz (d·E) Hamiltonian.

These correspond to Eqs. (2.50) and (2.51) of Chapter 2 in your course materials.

⚙️ The Coupled Hamiltonians#

The minimal coupling Hamiltonian is written as:

\[ \mathbf{H}_{p \cdot A} = \mathbf{H}_{\mathrm{m}} \otimes \mathbf{I}_{\mathrm{P}} + \mathbf{I}_{\mathrm{M}} \otimes \mathbf{H}_{\mathrm{c}} + \mathbf{I}_{\mathrm{M}} \otimes \mathbf{H}_{\mathrm{dia}} + \mathbf{H}_{\mathrm{int}} \]

The Pauli–Fierz Hamiltonian is written as:

\[ \mathbf{H}_{\mathrm{PF}} = \mathbf{H}_{\mathrm{m}} \otimes \mathbf{I}_{\mathrm{P}} + \mathbf{I}_{\mathrm{M}} \otimes \mathbf{H}_{\mathrm{c}} + \mathbf{H}_{\mathrm{dse}} \otimes \mathbf{I}_{\mathrm{P}} + \mathbf{H}_{\mathrm{int}} \]

🧠 Notes on the Components#

  • The bare matter Hamiltonian (\(\mathbf{H}_{\mathrm{m}}\)) is defined in Eq. (2.14).

  • The bare cavity Hamiltonian (\(\mathbf{H}_{\mathrm{c}}\)) is given in Eq. (2.18).

  • For the p·A form: the interaction and diamagnetic terms are given in Eq. (2.25).

  • For the d·E (Pauli–Fierz) form: the interaction and dipole self-energy (DSE) terms appear in Eq. (2.43).

Each Hamiltonian can be expressed in terms of creation and annihilation (ladder) operators for the matter and photonic subsystems.

🧮 Building the Ladder Operators#

A key first step in constructing the Hamiltonian matrices is to build matrix representations of the ladder operators for each subsystem.
To do this, we’ll define a helper function:

build_ladder_operator_matrix(dimension)


# 🧩 Step: Build matrix representations of ladder (creation/annihilation) operators

def build_ladder_operator_matrices(dimension):
    """
    Construct the matrix representations of the bosonic raising (creation) and lowering
    (annihilation) operators. These operators will be used for both:
      • the matter subsystem (within the harmonic oscillator approximation), and
      • the photonic subsystem.

    Parameters
    ----------
    dimension : int
        Dimension of the truncated Hilbert space (number of energy levels).

    Returns
    -------
    raise_mat : np.ndarray
        The raising (creation) operator matrix, shape (dimension, dimension)
    lower_mat : np.ndarray
        The lowering (annihilation) operator matrix, shape (dimension, dimension)

    Example
    -------
    >>> a_dag, a = build_ladder_operator_matrices(3)
    >>> a
    array([[0.        , 1.        , 0.        ],
           [0.        , 0.        , 1.41421356],
           [0.        , 0.        , 0.        ]])

    >>> a_dag
    array([[0.        , 0.        , 0.        ],
           [1.        , 0.        , 0.        ],
           [0.        , 1.41421356, 0.        ]])
    """

    # ✏️ TODO 1: Initialize the lowering operator as a dim × dim zero matrix
    # lower_mat = np.zeros((dimension, dimension), dtype=complex)

    # ✏️ TODO 2: Populate elements such that lower_mat[i, i+1] = sqrt(i+1)
    # 💡 Hint: Use a for-loop over range(dimension-1)
    # for i in range(dimension - 1):
    #     lower_mat[i, i+1] = np.sqrt(i + 1)

    # ✏️ TODO 3: Define the raising operator as the Hermitian adjoint of lowering
    # raise_mat = lower_mat.conj().T

    # ✏️ TODO 4: Return both matrices
    # return raise_mat, lower_mat

# 🧩 Step: Build and display ladder operator matrices for matter and photonic subsystems

# --- Matter subsystem: Hilbert space truncated to 4 states ---
a_dag, a = build_ladder_operator_matrices(4)

print("🧱 Matter Subsystem (4-level truncation)\n")
print("Raising operator (a†):")
print(np.round(a_dag, 3))
print("\nLowering operator (a):")
print(np.round(a, 3))
print("\n" + "="*60 + "\n")

# --- Photonic subsystem: Hilbert space truncated to 2 states ---
b_dag, b = build_ladder_operator_matrices(2)

print("💡 Photonic Subsystem (2-level truncation)\n")
print("Raising operator (b†):")
print(np.round(b_dag, 3))
print("\nLowering operator (b):")
print(np.round(b, 3))

⚛️ Details of the Matter Subsystem within the Harmonic Approximation#

We will model the matter subsystem as the vibrational structure of the hydrogen fluoride (HF) molecule within the harmonic approximation.
This section defines the physical parameters that enter both the bare matter Hamiltonian and the light–matter coupling terms in the two Hamiltonian forms.

🧮 Reduced Mass of HF#

In the coupled Hamiltonians:

  • \(z\) — effective charge of the HF molecule

  • \(\mu\) — reduced mass of HF, given by

\[ \mu = \frac{m_F \, m_H}{m_F + m_H} \]

The harmonic potential for HF is then:

\[ V(x) = \tfrac{1}{2} k x^2 = \tfrac{1}{2} \mu \omega^2 x^2 \]

⚙️ Unit System#

We will work entirely in atomic units (a.u.), so the following quantities must all be expressed in atomic units:

  • Masses of H and F

  • Harmonic force constant \(k\)

  • Cavity frequency \(\omega_{\mathrm{cav}}\)

  • Effective charge associated with the dipole operator \(z\)

  • Vector potential amplitude \(A\)

We will vary the vector potential amplitude, but keep all other parameters fixed as follows:

Quantity

Symbol

Value (a.u.)

Mass of F

\(m_F\)

34616.6811

Mass of H

\(m_H\)

1837.4731

Force constant of HF*

\(k\)

0.6377

Equilibrium bond length*

\(r_{\mathrm{eq}}\)

1.7310

Effective charge*

\(z\)

–0.4688

Cavity frequency*

\(\omega_{\mathrm{cav}}\)

0.01911

*Values obtained from ab initio quantum chemistry calculations.

🧠 Notes on Derived Quantities#

  • Force constant (\(k\)): obtained from numerical derivatives of CCSD(T)/cc-pVTZ energies around the optimized equilibrium geometry.

  • Effective charge (\(z\)): derived from the RHF/cc-pVTZ dipole moment at the optimized geometry, divided by the equilibrium bond length.

  • Cavity frequency (\(\omega_{\mathrm{cav}}\)): tuned to the fundamental vibrational frequency of HF, computed from the harmonic model:

\[ \omega_{\mathrm{cav}} = \omega_{\mathrm{f}} = \sqrt{\frac{k}{\mu}} \]

where \(k\) is the force constant and \(\mu\) is the reduced mass.

🧾 Next Step#

In the following code block, we will:

  1. Assign or compute these parameters in atomic units, and

  2. Plot the harmonic potential \(V(x)\) and the associated vibrational energy levels for HF.

# 🧮 Step: Define physical parameters for the HF molecule (in atomic units)

# Force constant of HF (in atomic units)
# Derived from finite-difference CCSD(T)/cc-pVTZ calculations about r_eq = 1.73106 a.u. (≈ 0.916 Å)
# ✏️ TODO: Define k_au here
# k_au = ...

# Atomic masses (in atomic units)
# ✏️ TODO: Define mF_au and mH_au here
# mF_au = ...
# mH_au = ...

# Reduced mass of HF
# μ = (m_F * m_H) / (m_F + m_H)
# ✏️ TODO: Compute mu_au here
# mu_au = ...

# Fundamental vibrational frequency of HF (in a.u.)
# ω = sqrt(k / μ)
# ✏️ TODO: Compute omega_au here
# omega_au = ...

# Set the cavity frequency equal to the fundamental frequency
# ✏️ TODO: Define omega_cav here
# omega_cav = ...

# Effective charge of HF (in atomic units)
# ✏️ TODO: Define z_au here
# z_au = ...

# Print computed cavity frequency
print(f"✅ Cavity frequency (atomic units): {omega_cav:.6f}")

# Reduced Planck's constant (arbitrary units for consistency)
hbar = 1.0

🎵 Visualizing the Harmonic Potential and Vibrational Energy Levels#

In this step, we’ll define a few helper functions to compute and plot the harmonic potential of HF, along with its lowest few energy eigenstates.

You’ll implement:

  • potential(omega, mu, x) — computes \(V(x) = \tfrac{1}{2}\mu \omega^2 x^2\)

  • energy(omega, n) — computes \(E_n = \hbar \omega (n + \tfrac{1}{2})\)

  • turning_points(omega, mu, E_n) — finds the classical turning points where \(V(x) = E_n\)

These will be used to:

  • plot the potential energy curve,

  • overlay the first few vibrational energy levels, and

  • visualize how quantization arises in the harmonic oscillator model.

# 🧩 Step: Define functions for potential energy and energy levels of the harmonic oscillator

def potential(omega, mu, x):
    """
    Compute the harmonic oscillator potential:
        V(x) = 1/2 * mu * omega^2 * x^2

    Parameters
    ----------
    omega : float
        Fundamental frequency (in a.u.)
    mu : float
        Reduced mass (in a.u.)
    x : float or np.ndarray
        Displacement coordinate(s)

    Returns
    -------
    V_x : float or np.ndarray
        Potential energy at displacement(s) x
    """
    # ✏️ TODO: Compute V_x = 1/2 * mu * omega**2 * x**2
    # V_x = ...

    # ✏️ TODO: Return V_x
    # return V_x


def energy(omega, n):
    """
    Compute the energy eigenvalue of the n-th harmonic oscillator level:
        E_n = hbar * omega * (n + 1/2)

    Parameters
    ----------
    omega : float
        Fundamental frequency (in a.u.)
    n : int
        Quantum number (n = 0, 1, 2, ...)

    Returns
    -------
    E_n : float
        Energy of level n (in a.u.)
    """
    hbar = 1.0  # in atomic units
    # ✏️ TODO: Compute E_n
    # E_n = ...

    # ✏️ TODO: Return E_n
    # return E_n


def turning_points(omega, mu, E_n):
    """
    Compute the classical turning points for a given energy E_n.
    Turning points satisfy:
        E_n = 1/2 * mu * omega^2 * x_0^2
        => x_0 = sqrt(2 * E_n / (mu * omega^2))

    Parameters
    ----------
    omega : float
        Fundamental frequency (in a.u.)
    mu : float
        Reduced mass (in a.u.)
    E_n : float
        Energy eigenvalue (in a.u.)

    Returns
    -------
    x_min, x_max : floats
        Negative and positive turning points
    """
    x_max = np.sqrt(2 * E_n / (mu * omega**2))
    return -x_max, x_max


# --- Plotting the potential and energy levels ---

# Define range for x (Bohr radii)
x = np.linspace(-0.5, 0.5, 1000)

# Compute potential energy curve
V = potential(omega_au, mu_au, x)

plt.figure(figsize=(7, 5))
plt.plot(x, V, label="Potential Energy", color="black")

# Overlay first few eigenstates
n_max = 3
for n in range(n_max):
    E_n = energy(omega_au, n)
    x_min, x_max = turning_points(omega_au, mu_au, E_n)
    plt.hlines(E_n, x_min, x_max, color="red", linestyle="--",
               label="Energy Levels" if n == 0 else None)

# --- Labels, legend, and style ---
plt.title("Harmonic Oscillator Potential with Energy Eigenstates", fontsize=13, weight="bold")
plt.xlabel("Displacement (Bohr radii)")
plt.ylabel("Energy (Hartrees)")
plt.ylim(0, energy(omega_au, n_max))
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

💭 Reflection: Understanding the Harmonic Potential#

Take a moment to interpret your plot of the harmonic oscillator potential and energy eigenlevels.

Questions to Consider:

  1. What do you notice about the spacing between the energy levels?

    • Is it uniform or changing with \(n\)?

    • How does this reflect the assumptions of the harmonic approximation?

  2. How do the turning points (where \(V(x) = E_n\)) change with \(n\)?

    • What does this say about the amplitude of vibration as the energy increases?

  3. If the potential were anharmonic (as in a Morse potential),

    • how would the level spacing and turning points differ?

💡 Hint:
In the harmonic model, each level is separated by a constant energy spacing of \(\hbar \omega\),
but in real molecular potentials, levels get closer together as \(n\) increases.

⚙️ Building the Bare Matter Hamiltonian#

With the fundamental molecular parameters now defined, we can build matrix representations of each term in the two Hamiltonians under consideration.

Each function we define will:

  • take arguments specifying the dimension of the matter and/or photonic subspaces, and

  • include all relevant parameters (\(\omega_f\), \(\omega_c\), \(E_0\), etc.)
    needed to construct the corresponding Hamiltonian matrix.

These functions will call on the helper function
build_ladder_operator_matrices() to generate the matrix representations of the ladder operators.

🧱 The Matter Hamiltonian#

The bare matter Hamiltonian is given by:

\[ \mathbf{H}_{\mathrm{m}} = \hbar \omega_f \left( \mathbf{A}^\dagger \mathbf{A} + \tfrac{1}{2}\mathbf{I} \right) \]

where:

  • \(\mathbf{A}^\dagger\) and \(\mathbf{A}\) are the creation and annihilation operator matrices in the matter Hilbert space, and

  • \(\mathbf{I}\) is the identity matrix in that subspace.

💡 Note:
In this simple harmonic oscillator model, \(\omega_f\) corresponds to the fundamental vibrational frequency of the HF molecule, computed as \(\omega_f = \sqrt{k / \mu}\).

# 🧩 Step: Build the bare matter Hamiltonian matrix

def compute_matter_energy_matrix(matter_dim, k, mu):
    """ 
    Build the matrix representation of the matter (vibrational) Hamiltonian.

    Parameters
    ----------
    matter_dim : int
        Dimension of the matter Hilbert space
    k : float
        Force constant (in atomic units)
    mu : float
        Reduced mass (in atomic units)

    Returns
    -------
    matter_Hamiltonian : np.ndarray
        Matrix representation of the bare matter Hamiltonian

    Example
    -------
    >>> compute_matter_energy_matrix(2, 1, 1)
    array([[0.5, 0.0],
           [0.0, 1.5]])
    """

    # ✏️ TODO 1: Build the ladder operators (a†, a) for the matter subspace
    # a_dag, a = build_ladder_operator_matrices(matter_dim)

    # ✏️ TODO 2: Compute the fundamental frequency ω_f = sqrt(k / μ)
    # omega_f = ...

    # ✏️ TODO 3: Build the matter Hamiltonian
    # H_m = omega_f * (a_dag @ a + 0.5 * np.eye(matter_dim))

    # ✏️ TODO 4: Return the Hamiltonian matrix
    # return H_m

💡 Building the Bare Cavity Hamiltonian#

Next, we’ll construct the bare cavity (photonic) Hamiltonian.
This term represents the energy stored in the cavity photon mode and is mathematically analogous to the harmonic oscillator Hamiltonian for the matter subsystem.

The Hamiltonian is given by:

\[ \mathbf{H}_{\mathrm{c}} = \hbar \omega_{\mathrm{c}} \left( \mathbf{B}^\dagger \mathbf{B} + \tfrac{1}{2}\mathbf{I} \right) \]

where:

  • \(\mathbf{B}^\dagger\) and \(\mathbf{B}\) are the creation and annihilation operator matrices for the photonic Hilbert space, and

  • \(\mathbf{I}\) is the identity matrix in that subspace.

Just as before, we’ll use the helper function build_ladder_operator_matrices() to generate these operators.

💡 Note:
\(\omega_{\mathrm{c}}\) is the cavity frequency, which we have already set equal to the fundamental vibrational frequency \(\omega_{\mathrm{f}}\) in the previous section.

# 🧩 Step: Build the bare cavity Hamiltonian matrix

def compute_cavity_energy_matrix(cavity_dim, omega_cav):
    """ 
    Build the matrix representation of the cavity (photonic) Hamiltonian.

    Parameters
    ----------
    cavity_dim : int
        Dimension of the cavity (photon) Hilbert space
    omega_cav : float
        Cavity frequency (in atomic units)

    Returns
    -------
    cavity_Hamiltonian : np.ndarray
        Matrix representation of the bare cavity Hamiltonian

    Example
    -------
    >>> compute_cavity_energy_matrix(2, 1)
    array([[0.5, 0.0],
           [0.0, 1.5]])
    """

    # ✏️ TODO 1: Build the ladder operators (b†, b) for the cavity subspace
    # b_dag, b = build_ladder_operator_matrices(cavity_dim)

    # ✏️ TODO 2: Construct the cavity Hamiltonian
    # H_c = omega_cav * (b_dag @ b + 0.5 * np.eye(cavity_dim))

    # ✏️ TODO 3: Return the Hamiltonian matrix
    # return H_c

🧲 Building the Diamagnetic Term for the p·A Hamiltonian#

The diamagnetic term appears only in the minimal coupling (p·A) Hamiltonian.
It accounts for the quadratic interaction between the electromagnetic vector potential and the charge distribution of the molecule.

The diamagnetic Hamiltonian is given by:

\[ \mathbf{H}_{\mathrm{dia}} = \frac{z^2}{2 \mu} \mathbf{A}_0^2 \left( \mathbf{B}^\dagger + \mathbf{B} \right)^2 \]

where:

  • \(z\) is the effective charge of the matter subsystem,

  • \(\mu\) is the reduced mass of the matter subsystem, and

  • \(\mathbf{A}_0\) is the magnitude of the vector potential for the cavity mode.

🧠 Physical Meaning#

This term ensures gauge invariance and prevents the total energy from becoming unphysically negative as the light–matter coupling strength increases.
It effectively contributes a photon-number-dependent energy shift.

💡 Note:
The ladder operators \(\mathbf{B}^\dagger\) and \(\mathbf{B}\) act in the photonic Hilbert space, and the matrix form of \((\mathbf{B}^\dagger + \mathbf{B})^2\) will depend on the dimension you choose for that space.

# 🧩 Step: Build the diamagnetic term for the p·A Hamiltonian

def compute_diamagnetic_matrix_p_dot_A(cavity_dim, z_charge, A0, mu):
    """ 
    Compute the matrix representation of the diamagnetic term in the *p·A* Hamiltonian.

    Parameters
    ----------
    cavity_dim : int
        Dimension of the photon Hilbert space
    z_charge : float
        Effective charge of the matter subsystem (in a.u.)
    A0 : float
        Magnitude of the cavity vector potential (in a.u.)
    mu : float
        Reduced mass of the matter subsystem (in a.u.)

    Returns
    -------
    diamag_Hamiltonian : np.ndarray
        Matrix representation of the diamagnetic term

    Example
    -------
    >>> compute_diamagnetic_matrix_p_dot_A(3, 1, 1, 1)
    array([[0.5      , 0.        , 0.70710678],
           [0.        , 1.5      , 0.        ],
           [0.70710678, 0.        , 1.        ]])
    """

    # ✏️ TODO 1: Build the ladder operators (b†, b) for the cavity subspace
    # b_dag, b = build_ladder_operator_matrices(cavity_dim)

    # ✏️ TODO 2: Compute prefactor  (z^2 / (2 * mu)) * A0^2
    # fac = ...

    # ✏️ TODO 3: Compute (b† + b)
    # b_sum = ...

    # ✏️ TODO 4: Compute (b† + b)^2  →  b_sum @ b_sum
    # b_sum_squared = ...

    # ✏️ TODO 5: Compute the diamagnetic Hamiltonian
    # H_dia = fac * b_sum_squared

    # ✏️ TODO 6: Return the result
    # return H_dia

⚡ Building the Interaction Term for the p·A Hamiltonian#

The interaction Hamiltonian represents the coupling between the molecular momentum operator and the quantized vector potential of the electromagnetic field.
It appears only in the minimal coupling (p·A) Hamiltonian and has the form:

\[ \mathbf{H}_{\mathrm{int}} = -i\,\frac{z}{\mu} \sqrt{\frac{\mu \hbar \omega_f}{2}}\, \mathbf{A}_0 \left( \mathbf{A}^\dagger - \mathbf{A} \right) \otimes \left( \mathbf{B}^\dagger + \mathbf{B} \right) \]

where:

  • \(z\) is the effective charge of the matter subsystem,

  • \(\mu\) is the reduced mass,

  • \(\omega_f = \sqrt{k / \mu}\) is the fundamental vibrational frequency, and

  • \(\mathbf{A}_0\) is the amplitude of the vector potential.

The tensor product \(\otimes\) ensures that this operator acts simultaneously in the matter and photon Hilbert spaces.

🧠 Physical Meaning#

This term describes the exchange of quanta (energy) between the light and matter subsystems —
it’s the mathematical origin of Rabi splitting and polariton formation in strongly coupled systems.

💡 Reminder:
This interaction term, together with the diamagnetic term \(\mathbf{H}_{\mathrm{dia}}\), ensures that the full p·A Hamiltonian remains gauge invariant.

# 🧩 Step: Build the interaction term for the p·A Hamiltonian

def compute_interaction_matrix_p_dot_A(matter_dim, cavity_dim, mu, k, z_charge, A0):
    """ 
    Build the matrix representation of the interaction term in the *p·A* Hamiltonian.

    Parameters
    ----------
    matter_dim : int
        Dimension of the matter Hilbert space
    cavity_dim : int
        Dimension of the photon Hilbert space
    mu : float
        Reduced mass of the matter subsystem (in a.u.)
    k : float
        Force constant of the matter subsystem (in a.u.)
    z_charge : float
        Effective charge of the matter subsystem (in a.u.)
    A0 : float
        Magnitude of the cavity vector potential (in a.u.)

    Returns
    -------
    interaction_Hamiltonian : np.ndarray
        Matrix representation of the interaction term

    Example
    -------
    >>> compute_interaction_matrix_p_dot_A(2, 2, 1, 1, 1, 0.1)
    array([[0.+0.j        , 0.+0.j        , 0.+0.j        , 0.+0.07071068j],
           [0.+0.j        , 0.+0.j        , 0.+0.07071068j, 0.+0.j        ],
           [0.+0.j        , 0.-0.07071068j, 0.+0.j        , 0.+0.j        ],
           [0.-0.07071068j, 0.+0.j        , 0.+0.j        , 0.+0.j        ]])
    """

    # ✏️ TODO 1: Build the ladder operators (a†, a) for the matter subspace
    # a_dag, a = build_ladder_operator_matrices(matter_dim)

    # ✏️ TODO 2: Build the ladder operators (b†, b) for the cavity subspace
    # b_dag, b = build_ladder_operator_matrices(cavity_dim)

    # ✏️ TODO 3: Compute (a† - a)
    # a_diff = ...

    # ✏️ TODO 4: Compute (b† + b)
    # b_sum = ...

    # ✏️ TODO 5: Compute the fundamental frequency ω_f = sqrt(k / μ)
    # omega_f = ...

    # ✏️ TODO 6: Compute prefactor  fac = -1j * (z / μ) * sqrt(μ * hbar * ω_f / 2) * A0
    # fac = ...

    # ✏️ TODO 7: Build the Kronecker product for the tensor interaction term
    # H_int = fac * np.kron(a_diff, b_sum)

    # ✏️ TODO 8: Return the result
    # return H_int

⚡ Building the Bilinear Interaction Term for the Pauli–Fierz (d·E) Hamiltonian#

The Pauli–Fierz (PF) Hamiltonian describes light–matter coupling in the electric dipole gauge.
Here, the interaction between the molecular dipole and the quantized electric field is expressed as:

\[ \mathbf{H}_{\mathrm{int}} = -\,z \frac{\hbar}{2} \sqrt{\frac{\omega_{\mathrm{cav}}}{\mu \omega_f}}\, \mathbf{E}_0 \left( \mathbf{A}^\dagger + \mathbf{A} \right) \otimes \left( \mathbf{B}^\dagger + \mathbf{B} \right) \]

where:

  • \(z\) — effective charge of the matter subsystem

  • \(\mu\) — reduced mass of the molecule

  • \(\omega_f = \sqrt{k / \mu}\) — fundamental vibrational frequency of the matter subsystem

  • \(\omega_{\mathrm{cav}}\) — cavity frequency

  • \(\mathbf{E}_0\) — amplitude of the quantized cavity electric field

🧠 Physical Interpretation#

This bilinear term describes dipole–field coupling:
the molecular dipole moment operator \((\mathbf{A}^\dagger + \mathbf{A})\) interacts directly with
the electric field operator \((\mathbf{B}^\dagger + \mathbf{B})\).

  • The strength of coupling depends on \(\sqrt{\omega_{\mathrm{cav}} / (\mu \omega_f)}\).

  • The prefactor ensures the correct scaling between field amplitude and molecular motion.

💡 Reminder:
Unlike the p·A Hamiltonian, this term does not include a diamagnetic contribution,
but the Pauli–Fierz form includes a dipole self-energy (DSE) term, which you’ll add later to maintain physical consistency.

# 🧩 Step: Build the bilinear light–matter coupling term for the Pauli–Fierz (d·E) Hamiltonian

def compute_interaction_matrix_PF(matter_dim, cavity_dim, omega_cav, mu, k, z_charge, A0):
    """
    Build the matrix representation of the bilinear interaction term in the Pauli–Fierz (d·E) Hamiltonian.

    Parameters
    ----------
    matter_dim : int
        Dimension of the matter Hilbert space
    cavity_dim : int
        Dimension of the photon Hilbert space
    omega_cav : float
        Cavity frequency (in a.u.)
    mu : float
        Reduced mass of the matter subsystem (in a.u.)
    k : float
        Force constant of the matter subsystem (in a.u.)
    z_charge : float
        Effective charge of the matter subsystem (in a.u.)
    A0 : float
        Magnitude of the vector potential (in a.u.)

    Returns
    -------
    interaction_Hamiltonian : np.ndarray
        Matrix representation of the PF bilinear coupling term
    """

    # ✏️ TODO 1: Compute the matter fundamental frequency ω_f = sqrt(k / μ)
    # omega_f = ...

    # ✏️ TODO 2: Compute the electric field amplitude from the vector potential
    # E0 = np.sqrt(2 * omega_cav) * A0

    # ✏️ TODO 3: Compute the prefactor:
    # fac = -z_charge * 0.5 * np.sqrt(omega_cav / (mu * omega_f)) * E0

    # ✏️ TODO 4: Build the ladder operators for the matter subsystem
    # a_dag, a = build_ladder_operator_matrices(matter_dim)

    # ✏️ TODO 5: Build the ladder operators for the photonic subsystem
    # b_dag, b = build_ladder_operator_matrices(cavity_dim)

    # ✏️ TODO 6: Compute (a† + a)
    # a_sum = ...

    # ✏️ TODO 7: Compute (b† + b)
    # b_sum = ...

    # ✏️ TODO 8: Compute the full tensor-product interaction term
    # H_int = fac * np.kron(a_sum, b_sum)

    # ✏️ TODO 9: Return the interaction Hamiltonian
    # return H_int

💫 Building the Dipole Self-Energy (DSE) Term for the Pauli–Fierz Hamiltonian#

The dipole self-energy (DSE) term appears only in the Pauli–Fierz (PF) form of the Hamiltonian.
It represents the energy stored in the polarization field created by the molecule’s own dipole moment.
This term is essential for ensuring that the PF Hamiltonian remains stable and gauge invariant.

The DSE Hamiltonian is written as:

\[ \mathbf{H}_{\mathrm{DSE}} = z^2 \frac{\hbar}{4 \mu \omega_f}\, \mathbf{E}_0^2 \left( \mathbf{A}^\dagger + \mathbf{A} \right)^2 \]

where:

  • \(z\) — effective charge of the matter subsystem

  • \(\mu\) — reduced mass of the matter subsystem

  • \(\omega_f = \sqrt{k / \mu}\) — fundamental vibrational frequency

  • \(\mathbf{E}_0\) — amplitude of the cavity electric field

🧠 Physical Meaning#

This term adds a quadratic dipole contribution to the molecular potential.
It ensures the Pauli–Fierz Hamiltonian produces physically meaningful (non-diverging) energy levels at strong coupling.
Without it, the model would overcount the attractive light–matter interaction.

💡 Hint:
The DSE operator acts only in the matter Hilbert space, unlike the bilinear coupling term which couples both the matter and photonic subspaces.

# 🧩 Step: Build the dipole self-energy term for the Pauli–Fierz Hamiltonian

def compute_dipole_self_energy_PF(matter_dim, omega_cav, mu, k, z_charge, A0):
    """ 
    Compute the matrix representation of the dipole self-energy operator in the matter basis.

    Parameters
    ----------
    matter_dim : int
        Dimension of the truncated matter Hilbert space
    omega_cav : float
        Cavity frequency (in a.u.)
    mu : float
        Reduced mass of the matter subsystem (in a.u.)
    k : float
        Force constant of the matter subsystem (in a.u.)
    z_charge : float
        Effective charge of the matter subsystem (in a.u.)
    A0 : float
        Magnitude of the cavity vector potential (in a.u.)

    Returns
    -------
    dipole_self_energy_matrix : np.ndarray
        Matrix representation of the dipole self-energy operator

    Example
    -------
    >>> compute_dipole_self_energy_PF(2, 1, 1, 1, 1, 0.1)
    array([[0.005, 0.   ],
           [0.   , 0.005]])
    """

    # ✏️ TODO 1: Compute the fundamental frequency ω_f = sqrt(k / μ)
    # omega_f = ...

    # ✏️ TODO 2: Compute the electric field amplitude from the vector potential
    # E0 = np.sqrt(2 * omega_cav) * A0

    # ✏️ TODO 3: Compute prefactor  z^2 * ħ / (4 * μ * ω_f) * E0^2
    # fac = ...

    # ✏️ TODO 4: Build matter ladder operators (a†, a)
    # a_dag, a = build_ladder_operator_matrices(matter_dim)

    # ✏️ TODO 5: Compute (a† + a)
    # a_sum = ...

    # ✏️ TODO 6: Square the sum:  (a† + a)^2
    # a_sum_sq = ...

    # ✏️ TODO 7: Compute the dipole self-energy matrix
    # H_dse = fac * a_sum_sq

    # ✏️ TODO 8: Return the matrix
    # return H_dse

🧩 Constructing and Comparing the Full Hamiltonians#

Now that we have defined all of the individual components, we can assemble the complete Hamiltonians for both the minimal coupling (p·A) and the Pauli–Fierz (d·E) forms.

We will:

  1. Construct each Hamiltonian in the tensor product basis of the matter and photon subspaces.

  2. Compute and compare their eigenvalues to confirm that they agree (as they must, when constructed correctly).

  3. Visualize how polariton states split away from the first excited vibrational energy of the bare HF molecule.

⚙️ Adjustable Parameters#

You can change the basis dimensions to explore convergence behavior:

  • matter_dim — size of the uncoupled matter basis

  • photon_dim — size of the uncoupled photon basis

Try experimenting with different values of A0_au, the vector potential magnitude, to see how the light–matter coupling strength affects the energy level splitting.

The following parameters were already set earlier in the notebook:

  • Force constant k_au

  • Reduced mass mu_au

  • Effective charge z_au

💡 Tip:
For testing, start with small Hilbert space dimensions (e.g., 2×2).
For higher precision, you can gradually increase them — but note that the matrix size grows quickly!

# 🧮 Step: Construct and compare the full Hamiltonians (p·A vs. Pauli–Fierz)

# --- Define Hilbert space dimensions ---
photon_dim = 2  # Try increasing up to 5–10 for convergence testing
matter_dim = 2  # Can increase to explore larger bases

# --- Identity operators for tensor product spaces ---
I_matter = np.eye(matter_dim)
I_photon = np.eye(photon_dim)

# --- Frequencies and vector potential ---
omega_p_au = np.sqrt(k_au / mu_au)  # photon frequency (matched to HF fundamental)
A0_au = 2.0                          # vector potential magnitude (in a.u.)

# --- Matter subsystem ---
H_matter = compute_matter_energy_matrix(matter_dim, k_au, mu_au)
H_matter_tp = np.kron(H_matter, I_photon)

# --- Photon subsystem ---
H_photon = compute_cavity_energy_matrix(photon_dim, omega_p_au)
H_photon_tp = np.kron(I_matter, H_photon)

# --- Diamagnetic term (p·A only) ---
H_diamag = compute_diamagnetic_matrix_p_dot_A(photon_dim, z_au, A0_au, mu_au)
H_diamag_tp = np.kron(I_matter, H_diamag)

# --- Interaction terms ---
H_interaction_pda = compute_interaction_matrix_p_dot_A(matter_dim, photon_dim, mu_au, k_au, z_au, A0_au)
H_interaction_PF = compute_interaction_matrix_PF(matter_dim, photon_dim, omega_p_au, mu_au, k_au, z_au, A0_au)

# --- Dipole self-energy (PF only) ---
H_dse_PF = compute_dipole_self_energy_PF(matter_dim, omega_p_au, mu_au, k_au, z_au, A0_au)
H_dse_PF_tp = np.kron(H_dse_PF, I_photon)

# --- Build total Hamiltonians ---
H_pda_tp = H_matter_tp + H_photon_tp + H_diamag_tp + H_interaction_pda
H_PF_tp  = H_matter_tp + H_photon_tp + H_dse_PF_tp + H_interaction_PF

# --- Compute eigenvalues and eigenvectors ---
vals_pda, vecs_pda = la.eigh(H_pda_tp)
vals_pf,  vecs_pf  = la.eigh(H_PF_tp)

# --- Display results ---
print("📈 Eigenvalues of the Pauli–Fierz (d·E) Hamiltonian:")
print(np.round(vals_pf, 8))

print("\n📈 Eigenvalues of the minimal coupling (p·A) Hamiltonian:")
print(np.round(vals_pda, 8))

# --- Optional: verify against expected minimal-basis values ---
expected_eigenvalues = np.array([0.00968507, 0.02673352, 0.03112254, 0.04817098])

if matter_dim == 2 and photon_dim == 2:
    print("\n✅ Minimal basis selected: comparing eigenvalues to expected results...")
    print(f"Pauli–Fierz matches expected: {np.allclose(vals_pf, expected_eigenvalues)}")
    print(f"p·A matches expected:          {np.allclose(vals_pda, expected_eigenvalues)}")

🌈 Visualizing Polaritonic Splitting and the Rabi Effect#

In this final step, we’ll replot the harmonic potential of HF along with the eigenvalues of the coupled system.

We’ll overlay:

  • The uncoupled (bare) harmonic oscillator energy levels in red, and

  • The coupled (Pauli–Fierz / p·A) polaritonic energy levels in purple.

🧠 Interpretation#

You should observe that:

  • The ground state of the coupled system lies slightly above the uncoupled ground state (the vacuum Lamb shift).

  • The first excited manifold splits into two states — one below and one above the bare first excited vibrational level.
    These are the upper and lower polaritons, and their energy separation is known as the Rabi splitting.

This splitting is the hallmark of vibrational strong coupling, where light and molecular vibrations hybridize to form new, mixed states of light and matter.

# 🪞 Step: Plot the harmonic potential and polaritonic energy levels

plt.figure(figsize=(7, 5))
plt.plot(x, V, label="Potential Energy (HF)", color="black")

# --- Overlay bare and coupled eigenvalues ---
n_max = 3  # number of energy levels to plot
for n in range(n_max):
    E_n = energy(omega_au, n)       # bare vibrational energy
    pol_En = vals_pda[n]            # coupled system energy (Pauli–Fierz / p·A)
    
    # Turning points for each level
    x_min, x_max = turning_points(omega_au, mu_au, E_n)
    pol_x_min, pol_x_max = turning_points(omega_au, mu_au, pol_En)
    
    # Plot horizontal lines for energy levels
    plt.hlines(E_n, x_min, x_max, color="red", linestyle="--",
               label="Uncoupled Levels" if n == 0 else None)
    plt.hlines(pol_En, pol_x_min, pol_x_max, color="purple", linestyle="--",
               label="Coupled Levels" if n == 0 else None)

# --- Aesthetics ---
plt.title("Polaritonic Splitting in Vibrational Strong Coupling", fontsize=13, weight="bold")
plt.xlabel("Displacement (Bohr radii)")
plt.ylabel("Energy (Hartrees)")
plt.ylim(0, energy(omega_au, n_max))
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

Printing matrices#

It can be instructive to inspect the forms of each submatrix for the p dot A and Pauli-Fierz Hamiltonian. The code block below repeats the builds for these matrices in a minimal basis for the matter and photon subsystems and prints each matrix with a heading.

# define dimensions of each subspace
photon_dim = 2 # can make this as small as 2 or as large as you dare... 100 is pretty large and will be slow
matter_dim = 2 # can make this as small as 2 or as large as you dare... 100 is pretty large and will be slow

# compute identity matrix on the matter subspace
I_matter = np.eye(matter_dim)

# compute identity matrix on the photon subspace
I_photon = np.eye(photon_dim)



# set k to be the force constant of HF in atomic units from the top of this notebook
k_val = k_au

# set mu to be the reduced mass of HF in atomic units from the top of this notebook
mu_val = mu_au

# set z to be the charge of HF in atomic units from the top of this notebook
z_val = z_au

# set the photon frequency to be equal to the fundamental frequency of HF
omega_p_val = np.sqrt( k_au / mu_au)

# assign a magnitude of the vector potential
A0_val = 2.0

# compute matter Hamiltonian in the matter basis
H_matter = compute_matter_energy_matrix(matter_dim, k_val, mu_val)

print("Bare matter Hamiltonian in the matter basis")
print(H_matter)

# compute matter Hamiltonian in the tensor product basis
H_matter_tp = np.kron(H_matter, I_photon)

print("Matter Hamiltonian in the tensor product basis")
print(H_matter_tp)

# compute photon Hamiltonian in the photon basis
H_photon = compute_cavity_energy_matrix(photon_dim, omega_cav)

print("Bare photon Hamiltonian in the photon basis")
print(H_photon)

# compute photon Hamiltonian in the tensor product basis
H_photon_tp = np.kron(I_matter, H_photon)

print("Photon Hamiltonian in the tensor product basis")
print(H_photon_tp)

# compute diamagnetic Hamiltonian in the photon basis
H_diamag = compute_diamagnetic_matrix_p_dot_A(photon_dim, z_val, A0_val, mu_val)

print("Diamagnetic Hamiltonian in the photon basis")
print(H_diamag)

# compute diamagnetic Hamiltonian in the tensor product basis
H_diamag_tp = np.kron(I_matter, H_diamag)

print("Diamagnetic Hamiltonian in the tensor product basis")
print(H_diamag_tp)

# compute the interaction Hamiltonian for p dot A
H_interaction_pda = compute_interaction_matrix_p_dot_A(matter_dim, photon_dim, mu_val, k_val, z_val, A0_val)

print("p dot A Interaction Hamiltonian in the tensor product basis")
print(H_interaction_pda)

# compute the total p dot A Hamiltonian in the tensor product basis
H_pda_tp = H_matter_tp + H_photon_tp + H_diamag_tp + H_interaction_pda

print("Total p dot A Hamiltonian in the tensor product basis")
print(H_pda_tp)


# compute the interaction Hamiltonian for PF
H_interaction_PF = compute_interaction_matrix_PF(matter_dim, photon_dim, omega_cav, mu_val, k_val, z_val, A0_val)

print("PF Interaction Hamiltonian in the tensor product basis")
print(H_interaction_PF)

# compute the dipole self energy for PF in the matter basis
H_dipole_self_energy_PF = compute_dipole_self_energy_PF(matter_dim, omega_cav, mu_val, k_val, z_val, A0_val)

print("Dipole self energy in the matter basis")
print(H_dipole_self_energy_PF)

# compute the dipole self energy for PF in the tensor product basis
H_dipole_self_energy_PF_tp = np.kron(H_dipole_self_energy_PF, I_photon)

print("Dipole self energy in the tensor product basis")
print(H_dipole_self_energy_PF_tp)

# compute the total Pauli-Fierz Hamiltonian
H_PF_tp = H_matter_tp + H_photon_tp + H_dipole_self_energy_PF_tp + H_interaction_PF

print("Total Pauli-Fierz Hamiltonian in the tensor product basis")
print(H_PF_tp)

Extensions to Anharmonic Oscillators#

The Harmonic approximation neglects important features of the vibrational structure of molecules. In particular, the Harmonic approximation does not permit dissociation of chemical bonds, does not adequately capture the fact that higher vibrational states are associated with average and most probable bond lengths that are stretched from equilibrium, does predict distinct vibrational transition energies seen in hot bands, and does not predict that overtone transitions can be spectroscopically excited.

While beyond the scope of this primer, we wish to point out that the code above can be adapted to model potentials that capture anhoarmonic effects and that can recover the features described above. There are a couple of trade offs that we wish to point out if one develops these models for general anharmonic oscillator potential. The first trade off is that the simple rules of the ladder operators that enabled simple construction of the bare matter Hamiltonian, the dipole self energy, and the coupling matrices will no longer be suitable for anharmonic oscillators, and so one will have to adapt the construction of these operators accordingly. The second tradeoff is that the eigenspectra of the p dot A and Pauli-Fierz Hamiltonians will typically not agree with each other at the same level of basis truncation. They will converge to the same spectra in the complete basis limit, but depending on how strong the coupling is (i.e. the magnitude of \({\bf A}_0\), agreement between the two spectra may require a very large number of basis states for the matter and photonic subsystems.

We direct the interested reader to python code treating the p dot A, d dot E, and Pauli-Fierz Hamiltonians using the Morse potential for the vibrational structure here with an associated notebook implementing this approach here

Modeling Electronic Strong Coupling#

While we have focused specifically on vibrational strong coupling in this notebook and in Chapter 2 of this primer, the formalism presented can also be applied to electronic strong coupling. We direct the interested reader to an example where the magnesium hydride cation is coupled through an electronic transition to a cavity mode at ultraviolet frequencies here. Just as in the Morse oscillator case, one will sacrifice simple analytical terms for the bilinear coupling and dipole self energy terms, but further, one will sacrifice analytical expressions even for the matter energy eigenstates. In the linked notebook, we obtain energies and transition dipole moments from the time-dependent density functional theory formalism for the magnesium hydride subsystem, and use those quantities to parameterize our coupled model. We also note that in this notebook, we effectively implement the Jaynes-Cummings model, which neglects the dipole self energy (and certain bilinear coupling elements known as the counter-rotating terms). We encourage the interested reader to work through this example, and consider adding dipole self energy and counter rotating terms to evaluate their impact for different coupling strengths!

Contents