Dear Qiskit team,

I guess there is an issue with the recently updated documentation of the Gaussian (https://qiskit.org/documentation/stubs/qiskit.pulse.library.Gaussian.html) and Drag (https://qiskit.org/documentation/stubs/qiskit.pulse.library.Drag.html) pulse. In the following, I will comment on both documentations. Please let me know if my observations are correct.

1) Gaussian:

With the auxiliary function

f’(x)=exp(-1/2 * (x - duration/2)^2/sigma^2) (identical to the current documentation),

the 'LiftedGaussian‘ should follow

f(x)=amp * (f'(x + 1/2) - f’(-1))/(1 - f’(-1)).

The difference to the current documentation is the '+1/2‘ shift, which is due to the function 'midpoint_sample‘ of qiskit.pulse.library.samplers.strategies.py. 'midpoint_sample' is needed for the function 'gaussian‘ (and 'drag') of qiskit.pulse.library.discrete.py and shifts the timesteps 'times‘ from np.arange(0, duration) to np.arange(1/2, duration + 1/2).

2) Drag:

In the definition of the function 'drag‘ in qiskit.pulse.library.continuous.py the Drag pulse is defined as 'gauss + 1j * beta * gauss_deriv‘. In the function 'gaussian_deriv‘ of the same file 'gauss_deriv‘ is defined as '- x / sigma * gauss‘. The 'LiftedDrag‘ is therefore proportional to the 'LiftedGaussian‘ defined in the function '_fix_gaussian_width‘ of the continuous.py file. The correct formula of the 'LiftedDrag‘ is

d(x)=f(x) * (1 + 1j * ( -(x + 1/2 - duration/2)/sigma^2) = amp * (f'(x + 1/2) - f’(-1))/(1 - f’(-1)) * (1 + 1j * ( -(x + 1/2 - duration/2)/sigma^2).

3) Verfification:

I numerically checked the accuracy of these observations by comparing them to the waveform samples of both Gaussian and Drag using the following code. The result is a sum of absolute deviations (g_dev and d_dev) of the order 1e-14.

'''
from qiskit import pulse
import numpy as np

d, a, s, b = 320, 0.9, 80, 30

g_data = pulse.Gaussian(duration=d, amp=a, sigma=s).get_waveform().samples
d_data = pulse.Drag(duration=d, amp=a, sigma=s, beta=b).get_waveform().samples

g_func = lambda x, duration, amp, sigma: amp * (np.exp(-1/2 * (x + 1/2 - duration/2)2/sigma2) - np.exp(-1/2 * (-1 - duration/2)2/sigma2))/(1 - np.exp(-1/2 * (-1 - duration/2)2/sigma2))

d_func = lambda x, duration, amp, sigma, beta: g_func(x, duration, amp, sigma) * (1 + 1j * beta * ( - (x + 1/2 - duration/2)/sigma**2))

times = np.arange(0, d)

g_dev = sum([abs(g_func(t, d, a, s) - g_data[t]) for t in times])
d_dev = sum([abs(np.real(d_func(t, d, a, s, b)) - np.real(d_data[t])) for t in times]) + 1j * sum([abs(np.imag(d_func(t, d, a, s, b)) - np.imag(d_data[t])) for t in times])

print(g_dev)
print(d_dev)
'‘'

Thanks for your help,
Tobias