0
# Phase Space Functions
1

2
Phase space representations and quasi-probability distributions for continuous variable quantum systems.
3

4
## Capabilities
5

6
### Wigner Function
7

8
Calculate and manipulate Wigner functions for quantum states.
9

10
```python { .api }
11
def wigner(rho: Qobj, xvec: ArrayLike, yvec: ArrayLike, g: float = np.sqrt(2), 
12
           sparse: bool = False, method: str = 'clenshaw') -> np.ndarray:
13
    """
14
    Calculate Wigner function on phase space grid.
15
    
16
    Parameters:
17
    - rho: Density matrix or state vector
18
    - xvec: Array of x-coordinates (position-like)
19
    - yvec: Array of y-coordinates (momentum-like)
20
    - g: Scaling factor (default √2)
21
    - sparse: Use sparse matrices
22
    - method: 'clenshaw', 'iterative', or 'laguerre'
23
    
24
    Returns:
25
    - ndarray: Wigner function W(x,y) on grid
26
    """
27

28
def wigner_transform(psi: Qobj, op: Qobj, g: float = np.sqrt(2)) -> Qobj:
29
    """
30
    Calculate Wigner transform of operator with respect to state.
31
    
32
    Parameters:
33
    - psi: Reference state
34
    - op: Operator to transform
35
    - g: Scaling factor
36
    
37
    Returns:
38
    - Qobj: Wigner-transformed operator
39
    """
40
```
41

42
### Q-Function (Husimi Distribution)
43

44
Calculate Q-function and related phase space distributions.
45

46
```python { .api }
47
def qfunc(rho: Qobj, xvec: ArrayLike, yvec: ArrayLike, g: float = np.sqrt(2), 
48
          sparse: bool = False) -> np.ndarray:
49
    """
50
    Calculate Q-function (Husimi distribution) on phase space grid.
51
    
52
    Parameters:
53
    - rho: Density matrix or state vector
54
    - xvec: Array of x-coordinates
55
    - yvec: Array of y-coordinates  
56
    - g: Scaling factor (default √2)
57
    - sparse: Use sparse matrices
58
    
59
    Returns:
60
    - ndarray: Q-function Q(α) = ⟨α|ρ|α⟩/π on grid
61
    """
62

63
class QFunc:
64
    """
65
    Q-function calculation class for efficient repeated calculations.
66
    """
67
    def __init__(self, rho: Qobj = None, xvec: ArrayLike = None, yvec: ArrayLike = None):
68
        """
69
        Initialize Q-function calculator.
70
        
71
        Parameters:
72
        - rho: Density matrix
73
        - xvec: X-coordinate array
74
        - yvec: Y-coordinate array
75
        """
76
    
77
    def qfunc(self, rho: Qobj = None) -> np.ndarray:
78
        """
79
        Calculate Q-function for given or stored state.
80
        
81
        Parameters:
82
        - rho: Density matrix (uses stored if None)
83
        
84
        Returns:
85
        - ndarray: Q-function values
86
        """
87
```
88

89
### Spin Phase Space Functions
90

91
Phase space representations for spin systems.
92

93
```python { .api }
94
def spin_wigner(rho: Qobj, theta: ArrayLike, phi: ArrayLike) -> np.ndarray:
95
    """
96
    Calculate spin Wigner function on sphere.
97
    
98
    Parameters:
99
    - rho: Spin density matrix
100
    - theta: Polar angle array (0 to π)
101
    - phi: Azimuthal angle array (0 to 2π)
102
    
103
    Returns:
104
    - ndarray: Spin Wigner function on sphere
105
    """
106

107
def spin_q_function(rho: Qobj, theta: ArrayLike, phi: ArrayLike) -> np.ndarray:
108
    """
109
    Calculate spin Q-function on sphere.
110
    
111
    Parameters:
112
    - rho: Spin density matrix
113
    - theta: Polar angle array
114
    - phi: Azimuthal angle array
115
    
116
    Returns:
117
    - ndarray: Spin Q-function on sphere
118
    """
119
```
120

121
### Usage Examples
122

123
```python
124
import qutip as qt
125
import numpy as np
126
import matplotlib.pyplot as plt
127

128
# Phase space grid
129
x_max = 4
130
N_points = 100
131
xvec = np.linspace(-x_max, x_max, N_points)
132
yvec = np.linspace(-x_max, x_max, N_points)
133

134
# Various quantum states for Wigner function analysis
135

136
# 1. Vacuum state
137
N_levels = 20
138
vacuum = qt.basis(N_levels, 0)
139
rho_vacuum = vacuum * vacuum.dag()
140

141
W_vacuum = qt.wigner(rho_vacuum, xvec, yvec)
142
print(f"Vacuum Wigner function:")
143
print(f"  Min value: {W_vacuum.min():.3f}")
144
print(f"  Max value: {W_vacuum.max():.3f}")
145
print(f"  Should be Gaussian with max ≈ 0.318")
146

147
# 2. Coherent state
148
alpha = 1.5 + 1.0j
149
coherent = qt.coherent(N_levels, alpha)
150
rho_coherent = coherent * coherent.dag()
151

152
W_coherent = qt.wigner(rho_coherent, xvec, yvec)
153
Q_coherent = qt.qfunc(rho_coherent, xvec, yvec)
154

155
print(f"Coherent state α = {alpha}:")
156
print(f"  Wigner min/max: [{W_coherent.min():.3f}, {W_coherent.max():.3f}]")
157
print(f"  Q-function min/max: [{Q_coherent.min():.3f}, {Q_coherent.max():.3f}]")
158

159
# Find peak positions
160
X, Y = np.meshgrid(xvec, yvec)
161
max_idx = np.unravel_index(np.argmax(W_coherent), W_coherent.shape)
162
x_peak, y_peak = xvec[max_idx[1]], yvec[max_idx[0]]
163
print(f"  Peak at (x,p) ≈ ({x_peak:.2f}, {y_peak:.2f})")
164
print(f"  Expected: ({np.sqrt(2)*alpha.real:.2f}, {np.sqrt(2)*alpha.imag:.2f})")
165

166
# 3. Squeezed state
167
squeeze_param = 0.8
168
squeezed = qt.squeeze(N_levels, squeeze_param) * vacuum
169
rho_squeezed = squeezed * squeezed.dag()
170

171
W_squeezed = qt.wigner(rho_squeezed, xvec, yvec)
172
print(f"Squeezed vacuum (r={squeeze_param}):")
173
print(f"  Wigner min/max: [{W_squeezed.min():.3f}, {W_squeezed.max():.3f}]")
174

175
# 4. Fock state (non-classical)
176
n_photons = 2
177
fock = qt.basis(N_levels, n_photons)
178
rho_fock = fock * fock.dag()
179

180
W_fock = qt.wigner(rho_fock, xvec, yvec)
181
Q_fock = qt.qfunc(rho_fock, xvec, yvec)
182

183
print(f"Fock state |{n_photons}⟩:")
184
print(f"  Wigner min/max: [{W_fock.min():.3f}, {W_fock.max():.3f}]")  
185
print(f"  Has negative values: {W_fock.min() < -1e-10}")  # Non-classical
186
print(f"  Q-function is always ≥ 0: {Q_fock.min() >= -1e-10}")
187

188
# 5. Cat state (superposition)
189
cat_amp = 1.5
190
cat_state = (qt.coherent(N_levels, cat_amp) + qt.coherent(N_levels, -cat_amp)).unit()
191
rho_cat = cat_state * cat_state.dag()
192

193
W_cat = qt.wigner(rho_cat, xvec, yvec)
194
print(f"Schrödinger cat state (α = ±{cat_amp}):")
195
print(f"  Wigner min/max: [{W_cat.min():.3f}, {W_cat.max():.3f}]")
196
print(f"  Exhibits interference fringes with negative values")
197

198
# 6. Thermal state
199
n_thermal = 1.5
200
rho_thermal = qt.thermal_dm(N_levels, n_thermal)
201

202
W_thermal = qt.wigner(rho_thermal, xvec, yvec)
203
Q_thermal = qt.qfunc(rho_thermal, xvec, yvec)
204

205
print(f"Thermal state (⟨n⟩ = {n_thermal}):")
206
print(f"  Wigner min/max: [{W_thermal.min():.3f}, {W_thermal.max():.3f}]")
207
print(f"  Q-function max: {Q_thermal.max():.3f}")
208

209
# Comparison of calculation methods
210
# Different methods for Wigner function calculation
211
W_clenshaw = qt.wigner(rho_coherent, xvec, yvec, method='clenshaw')
212
W_iterative = qt.wigner(rho_coherent, xvec, yvec, method='iterative')
213

214
method_diff = np.abs(W_clenshaw - W_iterative).max()
215
print(f"Clenshaw vs iterative method difference: {method_diff:.2e}")
216

217
# Spin Wigner functions
218
# Spin-1/2 system
219
j = 0.5
220
theta = np.linspace(0, np.pi, 50)
221
phi = np.linspace(0, 2*np.pi, 100)
222

223
# Spin-up state
224
spin_up = qt.spin_state(j, j)  # |j,j⟩ = |1/2,1/2⟩
225
rho_spin_up = spin_up * spin_up.dag()
226

227
W_spin_up = qt.spin_wigner(rho_spin_up, theta, phi)
228
Q_spin_up = qt.spin_q_function(rho_spin_up, theta, phi)
229

230
print(f"Spin-up Wigner function:")
231
print(f"  Min/max: [{W_spin_up.min():.3f}, {W_spin_up.max():.3f}]")
232
print(f"  Peak should be at north pole (θ=0)")
233

234
# Spin coherent state
235
theta_coherent, phi_coherent = np.pi/3, np.pi/4  # Pointing direction
236
spin_coherent = qt.spin_coherent(j, theta_coherent, phi_coherent)
237
rho_spin_coherent = spin_coherent * spin_coherent.dag()
238

239
W_spin_coherent = qt.spin_wigner(rho_spin_coherent, theta, phi)
240
print(f"Spin coherent state:")
241
print(f"  Direction: θ={theta_coherent:.2f}, φ={phi_coherent:.2f}")
242
print(f"  Wigner min/max: [{W_spin_coherent.min():.3f}, {W_spin_coherent.max():.3f}]")
243

244
# Spin-1 system (larger)
245
j1 = 1.0
246
spin_states_j1 = [qt.spin_state(j1, m) for m in [1, 0, -1]]  # |1,1⟩, |1,0⟩, |1,-1⟩
247

248
# Superposition state
249
superpos_j1 = (spin_states_j1[0] + spin_states_j1[2]).unit()  # (|1,1⟩ + |1,-1⟩)/√2
250
rho_superpos_j1 = superpos_j1 * superpos_j1.dag()
251

252
W_superpos_j1 = qt.spin_wigner(rho_superpos_j1, theta, phi)
253
print(f"Spin-1 superposition state:")
254
print(f"  Wigner min/max: [{W_superpos_j1.min():.3f}, {W_superpos_j1.max():.3f}]")
255

256
# Phase space integral properties
257
# Wigner function should integrate to 1
258
X, Y = np.meshgrid(xvec, yvec)
259
dx = xvec[1] - xvec[0]
260
dy = yvec[1] - yvec[0]
261

262
integral_W_vacuum = np.sum(W_vacuum) * dx * dy / np.pi
263
integral_W_coherent = np.sum(W_coherent) * dx * dy / np.pi
264
integral_Q_coherent = np.sum(Q_coherent) * dx * dy / np.pi
265

266
print(f"Phase space integrals (should be 1):")
267
print(f"  Vacuum Wigner: {integral_W_vacuum:.3f}")
268
print(f"  Coherent Wigner: {integral_W_coherent:.3f}")
269
print(f"  Coherent Q-function: {integral_Q_coherent:.3f}")
270

271
# Q-function class usage
272
qfunc_calc = qt.QFunc(rho_coherent, xvec, yvec)
273
Q_from_class = qfunc_calc.qfunc()
274

275
# Verify consistency
276
consistency_check = np.abs(Q_from_class - Q_coherent).max()
277
print(f"Q-function class consistency: {consistency_check:.2e}")
278

279
# Time evolution in phase space
280
# Harmonic oscillator evolution
281
H = qt.num(N_levels)  # Number operator
282
times = np.linspace(0, 2*np.pi, 5)  # One period
283

284
print(f"Time evolution of coherent state:")
285
for i, t in enumerate(times):
286
    # Evolve state
287
    U = (-1j * H * t).expm()
288
    psi_t = U * coherent
289
    rho_t = psi_t * psi_t.dag()
290
    
291
    # Calculate Wigner function
292
    W_t = qt.wigner(rho_t, xvec, yvec)
293
    
294
    # Find peak
295
    max_idx = np.unravel_index(np.argmax(W_t), W_t.shape)
296
    x_peak, y_peak = xvec[max_idx[1]], yvec[max_idx[0]]
297
    
298
    print(f"  t = {t:.2f}: peak at ({x_peak:.2f}, {y_peak:.2f})")
299

300
# The coherent state should rotate in phase space with frequency ω=1
301
# Expected position: α(t) = α * exp(-iωt) = α * exp(-it)
302
alpha_expected = alpha * np.exp(-1j * times[-1])
303
x_expected = np.sqrt(2) * alpha_expected.real
304
y_expected = np.sqrt(2) * alpha_expected.imag
305
print(f"  Expected final position: ({x_expected:.2f}, {y_expected:.2f})")
306

307
# Visualization example (pseudo-code, actual plotting would need matplotlib)
308
def plot_phase_space_comparison(states_dict, xvec, yvec):
309
    """
310
    Example function showing how to visualize phase space functions.
311
    """
312
    print("Phase space visualization:")
313
    for name, rho in states_dict.items():
314
        W = qt.wigner(rho, xvec, yvec)
315
        Q = qt.qfunc(rho, xvec, yvec)
316
        
317
        print(f"  {name}:")
318
        print(f"    Wigner: min={W.min():.3f}, max={W.max():.3f}")
319
        print(f"    Q-func: min={Q.min():.3f}, max={Q.max():.3f}")
320
        print(f"    Non-classical (W<0): {W.min() < -1e-10}")
321

322
# Example usage
323
states_to_plot = {
324
    'Vacuum': rho_vacuum,
325
    'Coherent': rho_coherent, 
326
    'Squeezed': rho_squeezed,
327
    'Fock n=2': rho_fock,
328
    'Cat state': rho_cat,
329
    'Thermal': rho_thermal
330
}
331

332
plot_phase_space_comparison(states_to_plot, xvec, yvec)
333
```
334

335
## Types
336

337
```python { .api }
338
# Phase space functions return numpy ndarrays containing
339
# the calculated quasi-probability distributions on the specified grids
340
```