Create Physical Review Letters-Style Plots in Python

Last Updated on December 7, 2024 by Max

Scientific publications demand high-quality figures to effectively communicate research findings, ensuring clarity, precision, and aesthetic appeal. Figures often serve as focal points for readers, encapsulating complex data in accessible visual formats.

Publications like Physical Review Letters (PRL) emphasize visually consistent, publication-grade plots to uphold their reputation for precision and professionalism. High-quality figures not only enhance understanding but also increase the likelihood of acceptance in reputable journals.

Creating such figures requires adherence to stringent design standards, including proper labeling, scalable fonts, and consistent color schemes, which enhance interpretability.

Python, with libraries like Matplotlib and Seaborn, offers versatile tools to generate professional-grade plots. This article guides you in efficiently crafting Physical Review Letters-style plots using Python.

Let us begin by exploring the figure guidelines for PRL publications.

General Recommendations:

  1. Clarity and Quality:
    • Figures must be sharp and clear to ensure readability. All elements, including text, symbols, and lines, should be distinctly visible.
    • Use high-resolution images for all figures. For raster graphics (e.g., photographs, heatmaps), maintain a resolution of at least 300 DPI.
  2. Accepted Formats:
    • Submit figures in standard formats such as EPS, PDF, PNG, or TIFF.
    • Use EPS or PDF for vector graphics and PNG or TIFF for images that require rasterization.
  3. Dimensions and Layout:
    • Figures should be designed to fit either a single-column width (8.6 cm) or double-column width (17.8 cm).
    • Maintain uniform font sizes and proportional element spacing within and across all figures.
  4. Font and Labels:
    • Employ clean, professional fonts such as Arial or Helvetica.
    • Font sizes should be legible and typically between 6–10 points to match the text in the manuscript.
    • Label all axes and key features clearly, ensuring consistency in symbols and notation.
  5. Color Usage:
    • Figures must be interpretable in both color and grayscale.
    • Use different line styles, markers, or patterns to differentiate data instead of relying on color alone.

Additional Specifics:

  1. Line Styles and Resolution:
    • Use a minimum line thickness of 0.5 points to ensure visibility.
    • Rasterized figures should have resolutions between 300 and 600 DPI to prevent pixelation.
  2. Composite Figures:
    • For figures with multiple panels, label each panel with (a), (b), (c), etc., in the upper-left corner.
    • Ensure that all subpanels are evenly spaced and aligned.
  3. Axes and Scale Information:
    • Axes should include evenly spaced tick marks and appropriate labels.
    • Include scale bars for images or micrographs where applicable.
  4. Legends:
    • Legends should be placed where they do not obstruct the figure’s content. Concise descriptions within the legend are preferred.

Density Plot

A density plot is a visual representation of the spatial distribution of a scalar field, often used to describe probabilities, intensities, or densities in physical systems. In quantum mechanics, density plots are commonly employed to visualize the squared magnitude of wave functions (\(|\psi(x, y)|^2\)), which represents the probability density of particles in a given region.

In scientific research, density plots are crucial for analyzing quantum wave functions, electromagnetic fields, and fluid dynamics, offering insight into spatial variations and phase structures.

Density plot of a 2D Gaussian
Figure: Density plot of a 2D Gaussian wave function.

Python Code:
The code below generates a density plot for a 2D quantum wave function, modeled as a Gaussian envelope with phase modulation (\( \psi(x, y) = e^{-\frac{x^2 + y^2}{2 \sigma^2}} e^{i(k_x x + k_y y + \theta)} \)). The squared magnitude of the wave function (\(|\psi(x, y)|^2\)) is computed and plotted using a contour fill map. The plot is styled to meet Physical Review Letters (PRL) standards, including LaTeX labels, single-column dimensions, and a high-quality colorbar for clarity.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define the computational grid
x = np.linspace(-5, 5, 500)
y = np.linspace(-5, 5, 500)
x, y = np.meshgrid(x, y)

# Define the wave function (example: 2D Gaussian with phase modulation)
kx, ky = 2, 3  # Wave vector components
sigma = 1.0  # Width of the Gaussian
theta = np.pi / 4  # Phase angle

wave_function = np.exp(-((x ** 2 + y ** 2) / (2 * sigma ** 2))) * np.exp(1j * (kx * x + ky * y + theta))
density = np.abs(wave_function) ** 2  # Density is the squared magnitude of the wave function

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width (3.375 inches for PRL)

# Plot the density
c = ax.contourf(x, y, density, levels=100, cmap='viridis', extend='both')

# Add a colorbar
cbar = fig.colorbar(c, ax=ax, orientation='vertical', pad=0.04)
cbar.set_label(r'Density', fontsize=12)
cbar.ax.tick_params(labelsize=10, direction='in', length=4, width=0.5)

# Set axis labels and ticks
ax.set_xlabel(r'$x \ \mathrm{(units)}$', fontsize=12)
ax.set_ylabel(r'$y \ \mathrm{(units)}$', fontsize=12)
ax.tick_params(axis='both', which='major', labelsize=10, direction='in', length=4, width=0.5)

# Adjust axis limits and aspect ratio
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
ax.set_aspect('equal', adjustable='box')

# Format the ticks to align with PRL style
ax.xaxis.set_major_locator(MaxNLocator(integer=True))
ax.yaxis.set_major_locator(MaxNLocator(integer=True))
# Add gridlines for better visualization

ax.grid(visible=True, which='major', color='gray', linestyle='--', linewidth=0.5, alpha=0.5)

# Reduce whitespace around the figure for publication
plt.tight_layout(pad=0.5)

# Show the plot
plt.show()

                

Contour Plots

A contour plot is a two-dimensional graphical representation of a three-dimensional scalar field, where constant values of the scalar field are represented by contour lines.

These plots are invaluable in visualizing potential energy surfaces, wave functions, or other physical properties in quantum mechanics. Contour plots highlight critical points, symmetry, and gradients in the data, making them an essential tool for analyzing quantum systems.

Figure: Contour plot of a 2D harmonic oscillator potential.

Python Code:
The 2D harmonic oscillator potential is given by the equation,
\[V(x, y) = \frac{1}{2}m\omega^2(x^2 + y^2),\]
where \(m\) is the mass of the particle and \(\omega\) is the angular frequency. This function represents a paraboloid potential well, and the contour lines correspond to regions of constant potential energy.

The following Python code generates a contour plot of the 2D harmonic oscillator potential. The plot adheres to Physical Review Letters (PRL) standards with proper axis labels, scalable fonts, and consistent formatting.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define the computational grid
x = np.linspace(-5, 5, 500)
y = np.linspace(-5, 5, 500)
x, y = np.meshgrid(x, y)

# Define the 2D harmonic oscillator potential
m = 1.0  # Mass of the particle
omega = 1.0  # Angular frequency
V = 0.5 * m * omega**2 * (x**2 + y**2)  # Potential energy

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm width)

# Plot the contour with sharp contrast
contour = ax.contour(x, y, V, levels=10, colors='black', linewidths=1.0)

# Label contours for better scientific clarity
ax.clabel(contour, inline=True, fontsize=8, fmt=r'$%.1f$')

# Set axis labels and ticks
ax.set_xlabel(r'$x \ \mathrm{(units)}$', fontsize=10)
ax.set_ylabel(r'$y \ \mathrm{(units)}$', fontsize=10)
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.8)

# Set aspect ratio and limits
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
ax.set_aspect('equal', adjustable='box')

# Minimize whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_style_plot.png", dpi=300)

# Show the plot
plt.show()
                

Vector Field Plots

A vector field plot is a graphical representation of a vector function, where vectors are depicted at various points in space to indicate their magnitude and direction. These plots are extensively used in physics to visualize electric and magnetic fields, velocity fields, or force distributions.

In quantum mechanics and matter-wave optics, vector fields play a crucial role in understanding gradient forces, phase gradients, and probability flux.

Figure: Probability current density for a two-dimensional Gaussian wavefunction with phase modulation.

Python Code:
In matter-wave optics, vector field plots can be used to visualize the probability current density, a fundamental quantity representing the flow of probability for quantum particles.

For a given wavefunction \(\psi(x, y)\), the probability current density \(\mathbf{J}\) is defined as:

\[\mathbf{J} = \frac{\hbar}{m} \mathrm{Im} \left( \psi^* \nabla \psi \right),\]

where \(\hbar\) is the reduced Planck’s constant, \(m\) is the mass of the particle, and \(\nabla\) is the gradient operator.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define the computational grid
x = np.linspace(-2, 2, 10)
y = np.linspace(-2, 2, 10)
x, y = np.meshgrid(x, y)

# Define the wavefunction (example: 2D Gaussian with phase modulation)
kx, ky = 2, 3  # Wave vector components
sigma = 1.0  # Width of the Gaussian
theta = np.pi / 4  # Phase angle

psi = np.exp(-((x ** 2 + y ** 2) / (2 * sigma ** 2))) * np.exp(1j * (kx * x + ky * y + theta))
psi_conj = np.conj(psi)  # Complex conjugate of the wavefunction

# Compute the gradient of the wavefunction
psi_x = np.gradient(psi, axis=1)
psi_y = np.gradient(psi, axis=0)

# Compute the probability current density
hbar = 1.0  # Reduced Planck's constant
m = 1.0  # Mass of the particle
Jx = (hbar / m) * np.imag(psi_conj * psi_x)
Jy = (hbar / m) * np.imag(psi_conj * psi_y)

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot the vector field
quiver = ax.quiver(
    x, y, Jx, Jy,
    color='black', angles='xy', scale_units='xy', scale=1.5,
    linewidth=0.5
)

# Add axis labels with proper LaTeX formatting
ax.set_xlabel(r'$x \ \mathrm{(units)}$', fontsize=10)
ax.set_ylabel(r'$y \ \mathrm{(units)}$', fontsize=10)

# Set ticks with inward direction and scientific style
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)

# Set axis limits and ensure equal aspect ratio
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)
ax.set_aspect('equal', adjustable='box')

# Reduce whitespace around the figure
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_style_probability_current.png", dpi=300)

# Show the plot
plt.show()
                

Phase Space Plots

A phase space plot visually represents the states of a dynamical system by showing the relationship between position and momentum (or other related variables).

In quantum mechanics and matter-wave optics, these plots are valuable tools for studying wave packets, quantum oscillators, and coherent states. They provide a link between classical and quantum descriptions, helping to illustrate concepts like coherence, entanglement, and quantum uncertainty in a clear and intuitive way.

Figure: Wigner distribution for a Gaussian wave packet.

Python Code:
The following Python code generates a phase space plot of the Wigner distribution for a 1D Gaussian wave packet.

The plot represents the Wigner distribution \( W(x, p) \), a quasi-probability function in quantum mechanics that describes the state of a Gaussian wave packet in phase space. It is mathematically defined as \( W(x, p) = \frac{1}{\pi \hbar} \exp\left(-\frac{(x – x_0)^2}{2 \sigma^2}\right) \exp\left(-\frac{2 \sigma^2 (p – p_0)^2}{\hbar^2}\right) \), where \( x \) and \( p \) are position and momentum, \( x_0 \) and \( p_0 \) are their respective expectation values, \( \sigma \) represents the width of the Gaussian in position space, and \( \hbar \) is the reduced Planck’s constant. The position term \( \exp\left(-\frac{(x – x_0)^2}{2 \sigma^2}\right) \) defines the spatial spread, while the momentum term \( \exp\left(-\frac{2 \sigma^2 (p – p_0)^2}{\hbar^2}\right) \) governs the momentum spread, which is inversely proportional to \( \sigma \) and depends on \( \hbar \). The prefactor \( \frac{1}{\pi \hbar} \) ensures the normalization of the Wigner distribution, making its integral over all \( x \) and \( p \) equal to one. Together, these components describe the quantum state in phase space, highlighting the wave packet’s behavior in both position and momentum domains.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define the computational grid
x = np.linspace(-5, 5, 300)
p = np.linspace(-5, 5, 300)
x, p = np.meshgrid(x, p)

# Define the parameters of the Gaussian wave packet
x0 = 0.0  # Position expectation value
p0 = 1.0  # Momentum expectation value
sigma = 1.0  # Width of the Gaussian
hbar = 1.0  # Reduced Planck's constant

# Compute the Wigner distribution
wigner = (1 / (np.pi * hbar)) * np.exp(-((x - x0)**2 / (2 * sigma**2))) * np.exp(-(2 * sigma**2 * (p - p0)**2) / hbar**2)

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.6))  # Single-column width for PRL (8.6 cm)

# Plot the Wigner distribution
c = ax.contourf(x, p, wigner, levels=50, cmap='viridis')

# Add a colorbar
cbar = fig.colorbar(c, ax=ax, orientation='vertical', pad=0.08)
cbar.set_label(r'Wigner Distribution $W(x, p)$', fontsize=10)
cbar.ax.tick_params(labelsize=10, direction='in', length=3, width=0.5)

# Set axis labels and ticks
ax.set_xlabel(r'$x \ \mathrm{(units)}$', fontsize=12)
ax.set_ylabel(r'$p \ \mathrm{(units)}$', fontsize=12)
ax.tick_params(axis='both', which='major', labelsize=12, direction='out', length=4, width=0.5)

# Adjust axis limits and aspect ratio
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
ax.set_aspect('equal', adjustable='box')

# Reduce whitespace around the figure for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_Wigner_Distribution.png", dpi=300)

# Show the plot
plt.show()
                

Logarithmic Plots

A logarithmic plot is a graphical representation in which one or both axes are scaled logarithmically, enabling the visualization of data spanning several orders of magnitude. Such plots are essential for revealing underlying trends and scaling behaviors that may remain hidden on linear axes.

In quantum mechanics and matter-wave optics, logarithmic plots find extensive use in analyzing tunneling probabilities, energy spectra, and intensity variations, offering critical insights into the behavior of quantum systems.

Python Code:
The tunneling probability \(T(E)\) for a particle of energy \(E\) encountering a potential barrier of height \(V_0\) and width \(a\) is given by:

\[T(E) = \exp\left(-2a \sqrt{\frac{2m}{\hbar^2}(V_0 – E)}\right),\]

for \(E < V_0\). Here, \(m\) is the mass of the particle, and \(\hbar\) is the reduced Planck’s constant.

The following Python code generates a logarithmic plot of \(T(E)\) as a function of \(E\), illustrating the exponential decrease in tunneling probability with increasing barrier height.

Figure: Tunneling probability versus energy for a finite potential barrier.
Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define constants
m = 1.0  # Mass of the particle
hbar = 1.0  # Reduced Planck's constant
V0 = 10.0  # Barrier height
a = 1.0  # Barrier width

# Define energy range
E = np.linspace(0.01, V0 - 0.01, 500)  # Energies below V0

# Compute tunneling probability
T = np.exp(-2 * a * np.sqrt((2 * m / hbar**2) * (V0 - E)))

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot the tunneling probability
ax.plot(E, T, color='black', linewidth=1.0)

# Set axis labels and scale
ax.set_xlabel(r'Energy $E \ \mathrm{(units)}$', fontsize=10)
ax.set_ylabel(r'Tunneling Probability $T(E)$', fontsize=10)
ax.set_xscale('linear')  # Linear scale for energy
ax.set_yscale('log')  # Logarithmic scale for tunneling probability

# Customize tick parameters
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)
ax.tick_params(axis='both', which='minor', labelsize=8, direction='in', length=2, width=0.5)

# Add gridlines for better readability
ax.grid(visible=True, which='major', linestyle='--', linewidth=0.5, alpha=0.7)

# Adjust axis limits and layout
ax.set_xlim(0, V0)
ax.set_ylim(1e-4, 1)
ax.set_aspect('auto', adjustable='box')

# Reduce whitespace and save for publication
plt.tight_layout(pad=0.3)
plt.savefig("PRL_Tunneling_Probability.png", dpi=300)

# Show the plot
plt.show()
                

Power Spectral Density (PSD) Plot

A Power Spectral Density (PSD) plot illustrates the distribution of power or intensity of a signal across its frequency components.

In quantum mechanics and matter-wave optics, PSD plots are essential tools for analyzing the frequency spectra of time-dependent wave functions, interference phenomena, and experimental noise. These plots provide critical insights into periodicities, coherence properties, and the energy distribution within quantum systems, which help to understand their dynamical behavior.

Figure: Power Spectral Density (PSD) of a time-dependent oscillatory signal.

Python Code:
The Gaussian wave packet in a harmonic potential is described by:
\[\psi(x, t) = \left(\frac{1}{\sqrt{2\pi\sigma^2(t)}}\right)^{1/2} \exp\left(-\frac{(x – x_0)^2}{4\sigma^2(t)} + i\left(kx – \frac{\omega t}{2}\right)\right),\]
where \(\sigma(t)\) is the time-dependent width, \(x_0\) is the initial position, \(k\) is the wave vector, and \(\omega\) is the angular frequency of the harmonic potential.

The following Python code computes and plots the PSD of the wave packet’s time evolution.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Parameters for the Gaussian wave packet
x0 = 0.0  # Initial position
k = 1.0  # Wave vector
omega = 30.0 * 2 * np.pi  # Angular frequency
sigma0 = 1.0  # Initial width
t_max = 10.0  # Maximum time
dt = 0.01  # Time step

# Time array
time = np.arange(0, t_max, dt)

# Simulate the wave packet's time evolution
amplitude = np.cos(omega * time)  # Example oscillatory component

# Compute the Fourier transform for PSD
fft_result = fft(amplitude)
freqs = fftfreq(len(time), dt)
psd = np.abs(fft_result) ** 2

# Retain only positive frequencies
positive_freqs = freqs[freqs >= 0]
positive_psd = psd[freqs >= 0]

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot the PSD
ax.plot(positive_freqs, positive_psd, color='black', linewidth=1.0)

# Set axis labels and scale
ax.set_xlabel(r'Frequency $f \ \mathrm{(Hz)}$', fontsize=10)
ax.set_ylabel(r'Power Spectral Density $|\psi(f)|^2$', fontsize=10)
ax.set_xscale('linear')
ax.set_yscale('log')

# Customize tick parameters
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)
ax.tick_params(axis='both', which='minor', labelsize=8, direction='in', length=2, width=0.5)

# Add gridlines for better readability
ax.grid(visible=True, which='major', linestyle='--', linewidth=0.5, alpha=0.7)

# Reduce whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_PSD.png", dpi=300)

# Show the plot
plt.show()
                

Heatmaps

A heatmap is a graphical representation of data in a matrix format, where individual values are encoded using a color scale.

In quantum mechanics and matter-wave optics, heatmaps are widely employed to visualize spatial or temporal variations in physical quantities such as intensity, probability densities, or energy distributions. These visualizations provide critical insights into phenomena like interference patterns, quantum wavefunctions, and the spatial structure of energy states, making them indispensable for analyzing complex quantum systems.

Figure: Interference pattern of matter waves through double slits.

Python Code:
The interference pattern for a double-slit experiment is described by the probability density:

\[P(x, y) \propto \left| e^{ik(x – d/2)^2/2y} + e^{ik(x + d/2)^2/2y} \right|^2,\]

where \(d\) is the slit separation, \(k = 2\pi / \lambda\) is the wave number, and \(y\) is the distance from the slits to the screen. The pattern results from the constructive and destructive interference of the two wavefunctions.

The following Python code generates a heatmap for the interference pattern.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define constants
wavelength = 1.0  # Wavelength of the matter wave
k = 2 * np.pi / wavelength  # Wave number
d = 1.0  # Separation between slits
y = 10.0  # Distance from the slits to the screen

# Define the computational grid
x = np.linspace(-5, 5, 500)  # Screen width
z = np.linspace(0, 10, 500)  # Distance from slits
x, z = np.meshgrid(x, z)

# Compute the interference pattern
slit1 = np.exp(1j * k * ((x - d / 2)**2) / (2 * z))
slit2 = np.exp(1j * k * ((x + d / 2)**2) / (2 * z))
interference = np.abs(slit1 + slit2)**2  # Probability density

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot the heatmap
c = ax.imshow(interference, extent=(-5, 5, 0, 10), origin='lower', aspect='auto', cmap='viridis')

# Add a colorbar
cbar = fig.colorbar(c, ax=ax, orientation='vertical', pad=0.04)
cbar.set_label(r'Probability Density $P(x, z)$', fontsize=10)
cbar.ax.tick_params(labelsize=8, direction='in', length=4, width=0.5)

# Set axis labels and ticks
ax.set_xlabel(r'$x \ \mathrm{(units)}$', fontsize=10)
ax.set_ylabel(r'$z \ \mathrm{(units)}$', fontsize=10)
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)
ax.tick_params(axis='both', which='minor', labelsize=8, direction='in', length=2, width=0.5)

# Adjust layout and reduce whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_Interference_Pattern.png", dpi=300)

# Show the plot
plt.show()
                

Parametric Plots

A parametric plot represents the relationship between two or more physical quantities as a function of a parameter. In quantum mechanics and matter-wave optics, such plots are indispensable for visualizing phase space trajectories, quantum state evolution, and correlations between conjugate variables like position and momentum. These plots provide critical insights into the dynamical behavior and underlying symmetries of quantum systems.

Figure : Phase space trajectory of a harmonic oscillator.

Python Code:
The trajectory of a particle in a harmonic potential is described by,
\[x(t) = A \cos(\omega t + \phi), \quad p(t) = -m \omega A \sin(\omega t + \phi),\]
where \(A\) is the amplitude, \(\omega\) is the angular frequency, \(\phi\) is the phase, and \(m\) is the mass of the particle. The parametric plot in the \(x\)-\(p\) plane forms a closed elliptical trajectory, representing the phase space dynamics of the particle.

The following Python code generates a parametric plot for the trajectory.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define constants
A = 1.0  # Amplitude of oscillation
omega = 2.0 * np.pi  # Angular frequency
phi = 0.0  # Initial phase
m = 1.0  # Mass of the particle

# Define the time array
t = np.linspace(0, 2 * np.pi, 500)

# Compute position and momentum as functions of time
x = A * np.cos(omega * t + phi)
p = -m * omega * A * np.sin(omega * t + phi)

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot the parametric trajectory
ax.plot(x, p, color='black', linewidth=1.0)

# Set axis labels and ticks
ax.set_xlabel(r'Position $x \ \mathrm{(units)}$', fontsize=10)
ax.set_ylabel(r'Momentum $p \ \mathrm{(units)}$', fontsize=10)
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)
ax.tick_params(axis='both', which='minor', labelsize=8, direction='in', length=2, width=0.5)

# Set axis limits and aspect ratio
ax.set_xlim(-1.2, 1.2)
ax.set_ylim(-1.2 * m * omega, 1.2 * m * omega)
#ax.set_aspect('equal', adjustable='box')

# Add gridlines for better visualization
ax.grid(visible=True, which='major', linestyle='--', linewidth=0.5, alpha=0.5)

# Reduce whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_Phase_Space_Trajectory.png", dpi=300)

# Show the plot
plt.show()

                

Eigenvalue Distribution Plot

An eigenvalue distribution plot represents the spectrum of eigenvalues of an operator or matrix, offering valuable insights into the properties of a system, such as energy levels, stability, and quantum transitions.

In quantum mechanics, such plots are widely used to analyze the energy spectra of quantum systems, characterize waveguide modes, and investigate the stability of dynamical systems. These visualizations provide a fundamental understanding of the discrete or continuous nature of eigenvalue distributions and their physical significance.

Figure: Discrete energy levels of a quantum harmonic oscillator.

Python Code:
The quantum harmonic oscillator has eigenvalues given by,
\[E_n = \hbar \omega \left(n + \frac{1}{2}\right), \quad n = 0, 1, 2, \dots\]
where \(\omega\) is the angular frequency and \(n\) is the quantum number. The eigenvalue distribution highlights the evenly spaced energy levels characteristic of the harmonic oscillator.

The following Python code generates an eigenvalue distribution plot for the first \(n = 10\) energy levels of a quantum harmonic oscillator.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define constants
hbar = 1.0  # Reduced Planck's constant
omega = 1.0  # Angular frequency
n_levels = 10  # Number of energy levels

# Compute eigenvalues
n = np.arange(n_levels)
eigenvalues = hbar * omega * (n + 0.5)

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL

# Plot the eigenvalues
ax.stem(n, eigenvalues, basefmt=" ", linefmt="black", markerfmt="ko", use_line_collection=True)

# Set axis labels and ticks
ax.set_xlabel(r'Quantum Number $n$', fontsize=10)
ax.set_ylabel(r'Energy Levels $E_n \ \mathrm{(units)}$', fontsize=10)
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)
ax.tick_params(axis='both', which='minor', labelsize=8, direction='in', length=2, width=0.5)

# Add gridlines for better visualization
ax.grid(visible=True, which='major', linestyle='--', linewidth=0.5, alpha=0.5)

# Set axis limits
ax.set_xlim(-0.5, n_levels - 0.5)
ax.set_ylim(0, eigenvalues[-1] + hbar * omega)

# Reduce whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_Energy_Levels.png", dpi=300)

# Show the plot
plt.show()
                

Phase-Amplitude Plots

A phase-amplitude plot shows the relationship between a wavefunction’s amplitude and phase, helping to analyze coherence, phase stability, and oscillatory behavior in quantum systems.

These plots are valuable for visualizing wavefunction dynamics, including oscillations and phase variations, providing insights into the properties of quantum states.

Figure : Phase-amplitude relationship showing a decaying amplitude with increasing phase.

Python Code:
The wave function of a coherent state in a harmonic oscillator can be described by,
\[\psi(x, t) = \exp\left(-\frac{x^2}{2\sigma^2}\right) \exp\left(i(kx – \omega t)\right),\]
where \(\sigma\) is the width of the Gaussian wave packet, \(k\) is the wave vector, and \(\omega\) is the angular frequency. The amplitude is related to the width \(\sigma\), and the phase varies linearly with time.

The following Python code generates a phase-amplitude plot for a coherent state.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define constants
omega = 2.0 * np.pi  # Angular frequency
sigma = 1.0  # Width of the wave packet
k = 1.0  # Wave vector
t = np.linspace(0, 2 * np.pi, 500)  # Time array

# Define amplitude and phase
amplitude = np.exp(-t / 2)  # Example decaying amplitude
phase = omega * t  # Linearly increasing phase

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot amplitude vs phase
ax.plot(phase, amplitude, color='black', linewidth=1.0)

# Set axis labels and ticks
ax.set_xlabel(r'Phase $\phi \ \mathrm{(radians)}$', fontsize=10)
ax.set_ylabel(r'Amplitude $|\psi(t)|$', fontsize=10)
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)
ax.tick_params(axis='both', which='minor', labelsize=8, direction='in', length=2, width=0.5)

# Add gridlines for better visualization
ax.grid(visible=True, which='major', linestyle='--', linewidth=0.5, alpha=0.5)

# Set limits for better visualization
ax.set_xlim(0, 2 * np.pi * omega)
ax.set_ylim(0, 1.2)

# Reduce whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_Phase_Amplitude.png", dpi=300)

# Show the plot
plt.show()
                

Histogram Plots

A histogram plot represents the distribution of data by grouping values into bins and displaying their frequencies.

In quantum mechanics, histogram plots are widely used to analyze measurement outcomes, visualize probability distributions, and study the statistical properties of eigenvalues or energy levels. These plots provide a clear and intuitive representation of data, offering insights into the underlying quantum phenomena and statistical behaviors of physical systems.

Figure: Histogram of sampled positions with theoretical Gaussian fit.

Python Code:
The probability density for the position measurement of a Gaussian wavefunction is given by,
\[P(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{x^2}{2\sigma^2}\right),\]
where \(\sigma\) is the standard deviation of the Gaussian wavefunction. Sampling from this probability density generates position measurements, which can be visualized as a histogram.

The following Python code generates a histogram plot for position measurements sampled from a Gaussian wave function.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define constants
sigma = 1.0  # Standard deviation of the Gaussian wavefunction
num_samples = 10000  # Number of position measurements

# Generate position measurements sampled from a Gaussian distribution
positions = np.random.normal(loc=0, scale=sigma, size=num_samples)

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot the histogram
n, bins, patches = ax.hist(
    positions, bins=50, density=True, color='#4A90E2', alpha=0.7, edgecolor='black'
)

# Plot the theoretical Gaussian probability density function for comparison
x = np.linspace(-4 * sigma, 4 * sigma, 1000)
P_x = (1 / (np.sqrt(2 * np.pi * sigma**2))) * np.exp(-x**2 / (2 * sigma**2))
ax.plot(x, P_x, color='darkred', linestyle='-', linewidth=1.5, label='Theoretical PDF')

# Set axis labels and legend
ax.set_xlabel(r'Position $x \ \mathrm{(units)}$', fontsize=10)
ax.set_ylabel(r'Probability Density $P(x)$', fontsize=10)
ax.legend(fontsize=8, loc='upper left')

# Customize tick parameters
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)
ax.tick_params(axis='both', which='minor', labelsize=8, direction='in', length=2, width=0.5)

# Add gridlines for better visualization
ax.grid(visible=True, which='major', linestyle='--', linewidth=0.5, alpha=0.5)

# Set limits for better visualization
ax.set_xlim(-4 * sigma, 4 * sigma)
ax.set_ylim(0, max(n) * 1.2)

# Reduce whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_Histogram_Plot.png", dpi=300)

# Show the plot
plt.show()
                

Lyapunov Exponent Plot

A Lyapunov exponent plot quantifies the sensitivity of a dynamical system to initial conditions, providing a measure of its stability and predictability. Positive Lyapunov exponents indicate chaotic behavior, while negative or zero values correspond to regular dynamics.

In quantum mechanics, Lyapunov exponents are particularly useful in studying the stability of quantum systems, semiclassical dynamics, and signatures of quantum chaos. These plots serve as a critical tool for analyzing the interplay between classical and quantum stability in complex systems.

Figure: Lyapunov exponent versus driving frequency for driven oscillator.

Python Code:
The time evolution of a particle in a driven harmonic oscillator is described by,
\[\ddot{x} + \omega^2 x = A \sin(\Omega t),\]
where \(\omega\) is the natural frequency, \(A\) is the amplitude of the driving force, and \(\Omega\) is the driving frequency. Small differences in initial conditions can grow exponentially in chaotic regimes, quantified by the Lyapunov exponent,
\[\lambda = \lim_{t \to \infty} \frac{1}{t} \ln \frac{|\Delta x(t)|}{|\Delta x(0)|}.\]
The Lyapunov exponent is plotted as a function of the driving frequency \(\Omega\).

The following Python code computes and plots the Lyapunov exponent for a range of driving frequencies.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Parameters
omega = 1.0  # Natural frequency
A = 1.0  # Amplitude of driving force
t_span = (0, 100)  # Time span
dt = 0.01  # Time step
initial_conditions = [1.0, 0.0]  # Initial position and velocity

# Driving frequencies to evaluate Lyapunov exponents
driving_frequencies = np.linspace(0.5, 2.0, 50)
lyapunov_exponents = []

# Function to compute the system's dynamics
def driven_harmonic_oscillator(t, y, omega, A, Omega):
    x, v = y
    dxdt = v
    dvdt = -omega**2 * x + A * np.sin(Omega * t)
    return [dxdt, dvdt]

# Calculate Lyapunov exponent for each driving frequency
for Omega in driving_frequencies:
    # Solve the differential equation
    solution = solve_ivp(
        driven_harmonic_oscillator,
        t_span,
        initial_conditions,
        args=(omega, A, Omega),
        t_eval=np.arange(t_span[0], t_span[1], dt),
        rtol=1e-6,
        atol=1e-9
    )
    x = solution.y[0]
    delta_x = np.abs(x[1:] - x[:-1])  # Approximate differences
    delta_x[delta_x == 0] = 1e-10  # Avoid log(0)
    lyapunov_exponent = np.mean(np.log(delta_x)) / dt
    lyapunov_exponents.append(lyapunov_exponent)

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot the Lyapunov exponents
ax.plot(driving_frequencies, lyapunov_exponents, color='#1f77b4', linewidth=1.5, label='Lyapunov Exponent')

# Set axis labels
ax.set_xlabel(r'Driving Frequency $\Omega \ \mathrm{(units)}$', fontsize=10)
ax.set_ylabel(r'Lyapunov Exponent $\lambda$', fontsize=10)

# Customize tick parameters
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)

# Add gridlines for better visualization
ax.grid(visible=True, which='major', linestyle='--', linewidth=0.5, alpha=0.7)

# Add legend for clarity
ax.legend(fontsize=8, loc='upper right', frameon=False)

# Reduce whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_Lyapunov_Exponents.png", dpi=300)

# Show the plot
plt.show()
                

Fidelity Plot

A fidelity plot quantifies the similarity between two quantum states, typically as a function of time, system parameters, or external perturbations.

Fidelity serves as a critical measure in quantum mechanics and quantum information theory, providing insights into the coherence, stability, and robustness of quantum states and operations. Such plots are essential for evaluating quantum protocols, characterizing decoherence, and assessing the performance of quantum systems.

Figure: Fidelity versus time for perturbed quantum state evolution.

Python Code:
In quantum mechanics, fidelity is often used to evaluate the effect of perturbations on quantum systems. For instance, a wave packet evolving in a harmonic oscillator potential can experience deviations due to external perturbations or decoherence. Fidelity quantifies how well the perturbed state matches the reference (unperturbed) state.

The fidelity between two quantum states \(|\psi_1\rangle\) and \(|\psi_2\rangle\) is defined as,
\[F(t) = |\langle \psi_1(t) | \psi_2(t) \rangle|^2,\]
where \(F(t)\) ranges from 0 (orthogonal states) to 1 (identical states).

Consider a Gaussian wave packet in a harmonic potential where the reference state evolves under an unperturbed Hamiltonian, and the perturbed state evolves under a Hamiltonian with an additional sinusoidal driving force. The fidelity between these two states changes over time, reflecting the effect of the perturbation.

The following Python code generates a fidelity plot for this scenario.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define constants
omega = 2.0 * np.pi  # Angular frequency of the harmonic oscillator
epsilon = 0.05  # Perturbation strength
t = np.linspace(0, 10, 500)  # Time array

# Define reference and perturbed states
def reference_state(t):
    return np.exp(-1j * omega * t)  # Simplified reference evolution

def perturbed_state(t, epsilon):
    return np.exp(-1j * omega * t) * np.exp(-1j * epsilon * np.sin(omega * t))  # Perturbed evolution

# Compute fidelity
fidelity = np.abs(np.conj(reference_state(t)) * perturbed_state(t, epsilon))**2

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot fidelity
ax.plot(t, fidelity, color='#1f77b4', linewidth=1.5, label=r'Fidelity $F(t)$')

# Set axis labels
ax.set_xlabel(r'Time $t \ \mathrm{(units)}$', fontsize=10)
ax.set_ylabel(r'Fidelity $F(t)$', fontsize=10)

# Customize tick parameters
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)
ax.tick_params(axis='both', which='minor', labelsize=8, direction='in', length=2, width=0.5)

# Add gridlines for better visualization
ax.grid(visible=True, which='major', linestyle='--', linewidth=0.5, alpha=0.7)

# Add a legend
ax.legend(fontsize=8, loc='lower left', frameon=False)

# Set limits for better visualization
ax.set_xlim(0, t[-1])
ax.set_ylim(0, 1.05)

# Reduce whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_Fidelity_Plot.png", dpi=300)

# Show the plot
plt.show()
                

State Probability Distribution Plot

A state probability distribution plot depicts the probabilities of a quantum system being in specific eigenstates or configurations as a function of time or system parameters. Such plots are fundamental in quantum mechanics for analyzing state transitions, energy level populations, and the impact of external perturbations. They provide critical insights into the dynamics and stability of quantum systems, enabling a deeper understanding of coherence, decoherence, and quantum evolution.

Figure: Time evolution of state probabilities in quantum system.

Python Code:
In quantum mechanics, the state probability distribution is often used to analyze the distribution of a system’s wavefunction over a set of eigenstates of an observable, such as energy or position. For example, a quantum harmonic oscillator initialized in a superposition of eigenstates can have a dynamically evolving probability distribution over time.
The probability of finding the system in the \(n\)-th eigenstate is given by,

\[P(n) = |\langle \phi_n | \psi(t) \rangle|^2,\]
where \(|\phi_n\rangle\) is the \(n\)-th eigenstate of the system and \(|\psi(t)\rangle\) is the time-evolving wavefunction.

Consider a quantum harmonic oscillator initialized in a superposition of the ground state and the first excited state. The time-evolved wavefunction is,
\[|\psi(t)\rangle = \frac{1}{\sqrt{2}} \left( |\phi_0\rangle + e^{-iE_1t/\hbar} |\phi_1\rangle \right),\]
where \(|\phi_0\rangle\) and \(|\phi_1\rangle\) are the ground and first excited states, respectively, and \(E_1\) is the energy of the first excited state. The probability distribution evolves periodically in time.

The following Python code generates the state probability distribution plot for this scenario.

Wave Packet Simulation

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc

# Use LaTeX fonts for PRL-style visualization
rc('font', **{'family': 'serif', 'serif': ['Times']})
rc('text', usetex=True)

# Define constants
hbar = 1.0  # Reduced Planck's constant
omega = 1.0  # Angular frequency of the harmonic oscillator
t = np.linspace(0, 10, 500)  # Time array

# Compute state probabilities
P0 = 0.5  # Constant probability for the ground state
P1 = 0.5 * (1 + np.cos(omega * t))  # Time-dependent probability for the first excited state

# Create the figure and axis
fig, ax = plt.subplots(figsize=(3.375, 2.9))  # Single-column width for PRL (8.6 cm)

# Plot the probabilities
ax.plot(t, P0 * np.ones_like(t), label=r'$P(0)$', color='#1f77b4', linewidth=1.5)  # Blue for ground state
ax.plot(t, P1, label=r'$P(1)$', color='#ff7f0e', linewidth=1.5)  # Orange for first excited state

# Set axis labels
ax.set_xlabel(r'Time $t \ \mathrm{(units)}$', fontsize=10)
ax.set_ylabel(r'Probability $P(n)$', fontsize=10)

# Add a legend
ax.legend(fontsize=8, loc='upper right', frameon=False)

# Customize tick parameters
ax.tick_params(axis='both', which='major', labelsize=8, direction='in', length=4, width=0.5)

# Add gridlines for better visualization
ax.grid(visible=True, which='major', linestyle='--', linewidth=0.5, alpha=0.7)

# Set limits for better visualization
ax.set_xlim(0, t[-1])
ax.set_ylim(0, 1.05)

# Reduce whitespace for publication
plt.tight_layout(pad=0.3)

# Save the plot as a high-resolution image for publication
plt.savefig("PRL_State_Probabilities.png", dpi=300)

# Show the plot
plt.show()
                

This article provides a guide to creating publication-grade figures following Physical Review Letters (PRL) standards, focusing on clarity, precision, and adherence to formatting guidelines.

Using Python libraries like Matplotlib, it demonstrates generating high-quality plots, including density, contour, vector field, phase space, and eigenvalue distributions, with examples drawn from quantum mechanics.

Key techniques include LaTeX integration, scalable fonts, and high-resolution outputs. The article highlights the importance of visually consistent and scientifically rigorous figures in communicating complex data effectively.

These methods aim to help researchers produce impactful visuals that meet the exacting standards of high-impact scientific journals like PRL.

  1. Physical Review Journals Information for Authors, American Physical Society, https://journals.aps.org/revtex (accessed November 30, 2024).
  2. American Physical Society (APS) Submission Guidelines, https://journals.aps.org/authors (accessed November 30, 2024).
  3. Springer Submission Guidelines, Springer Nature, Tips for Preparing Figures for Publication, available at https://www.springernature.com/gp/authors/campaigns/submission-guidelines (accessed November 30, 2024).

Leave a Comment

Your email address will not be published. Required fields are marked *

Scroll to Top