Introducing Qiskit Dynamics, a New Qiskit Module for Simulating Quantum Systems

By Abby Mitchell, Qiskit Developer Advocate, and Daniel Puzzuoli, Software Developer

In October of 2021, the Qiskit team released a first version of the Qiskit Dynamics module. The team developed this module to enable continuous-time simulation of quantum systems. In this blog, we’ll walk you through the basic need-to-know information to help you get started using this shiny new package!

We’ll cover:

What is Qiskit Dynamics?

What’s in the first release?

• The Signal Class
• The Solver Interface
• Optimize your Calculations using JAX

What’s Coming Next?

What is Qiskit Dynamics?

We describe the physics of quantum systems using differential equations, the most well-known of which is the Schrödinger equation. A key component of understanding these systems is building models of them in terms of these differential equations, which we can then solve to gain insight into how the state of the system changes over time, also known as its evolution. Through studying these models, and building successively better models, we gain insight into the real-world behavior of quantum systems, such as quantum processors. We can use these insights to design better gates, improve the overall performance of quantum devices, or study other relevant problems in quantum physics.

Qiskit Dynamics aims to provide tools for building and numerically solving models of quantum systems, including tools for automatically performing common mathematical transformations on these models.

A key design decision in Qiskit Dynamics’ development was that it be modular and extensible with respect to the numerical methods used. This gives users the flexibility to adapt or switch between different numerical methods to suit their research needs, and will enable easy integration of new methods in the future. Furthermore, we have ensured that at every stage, a user can choose to use the JAX Python library for numerical computing to represent and solve models, making Qiskit Dynamics code automatically differentiable, just-in-time compilable, and executable on GPUs.

What’s in the first release?

This first release aims to establish the core functionality of the package, and lay the groundwork for future enhancements. The key features enable users to construct, solve, and transform models of quantum systems, and offers users the flexibility to choose different numerical methods depending on their needs.

Simulating quantum systems using Qiskit Dynamics generally consists of 4 parts:

1. Preparing the inputs for a model, i.e. creating objects representing the different terms in the model. Currently both Hamiltonian and Lindblad (dissipative) terms are supported.
2. Constructing the model using the prepared inputs.
3. Applying any transformations and approximations (e.g. Rotating Frames, Rotating Wave Approximation) to that model.
4. Solving the model over a specified time interval.

For example, let’s say we want to use Qiskit Dynamics to solve the Schrödinger equation for a system:

To do so we need to construct the components of the Hamiltonian, represented by H(t) in the above equation.

Qiskit Dynamics breaks a Hamiltonian down into the following standard linear decomposition:

Hence, to specify a Hamiltonian, we need to build the above components.

Signal Class

To prepare the signals for our Hamiltonian we can use Qiskit Dynamics’ built in Signal Class.

A Signal object represents a complex mixed signal (i.e. a combination of envelope and carrier waveforms) that can be represented by a time dependent function of the form:

Where f(t) represents the envelope and ν represents the carrier frequency (i.e. the frequency of the carrier waveform).

To construct a Signal object we first define the envelope (in this case a Gaussian) with the specifications of our choosing:

Then we construct the signal object by calling the Signal class with our envelope and a carrier frequency of our choosing:

from qiskit_dynamics import Signal

w = 1. # carrier frequency
gauss_signal = Signal(envelope=gaussian_envelope, carrier_freq=w) 

If you would like to visualize the signal you can do so using the .draw() function, specifying the initial time (t0), final time (tf) and the number of points to sample in that time frame (n).

gauss_signal.draw(t0=0, tf=T, n=100)

Note: you can also set the function argument to 'signal', 'envelope', or 'complex_value' depending on which waveform you want to visualize.

Enabling developers to construct Signal objects with arbitrary envelopes was an important design decision made in order to facilitate a level of customization necessary for more advanced simulations, e.g. including arbitrary signal processing.

Static Hamiltonian and Hamiltonian Operators

To construct the objects we need for the static Hamiltonian and Hamiltonian operators, we can use the Operator class from Qiskit (if you've used Qiskit's quantum_info module this should look familiar to you):

To refresh your knowledge on the Qiskit Operator Class you can check out the documentation here.

Solver Class

Once we’ve constructed the inputs from step 1, then we can combine the remaining steps into one by using the Solver Class, which is a higher level interface that automatically combines various lower level components, like model construction, transformations, and solving, as well as automatic handling of various quantum_info state types when solving. While Qiskit Dynamics does offer an interface layer for building model objects and working with them directly, for most users we recommend using the Solver class as best practice.

Constructing the Solver class with the Hamiltonian model data above simply requires passing the components into the constructor:

Once the solver has been initialized, we can finally solve for the time evolution of our model by calling the solve function.

Note: As we supplied only Hamiltonian information, calls to solve will integrate the Schrödinger equation, whereas if we had supplied dissipators as well, it would integrate the Lindblad equation.

For example, here we solve over the course of a 𝜋 -pulse, starting with an initial point y0 representing the ground state of the Hamiltonian, and saving snapshots of the state over 500 points between 0 and T.

You can access the quantum state values for any point in your solution:

To visualize the results, we plot the population in the excited state over the course of the time evolution of the pulse.

import matplotlib.pyplot as plt
excited_state_pop = [np.abs(y.data[0])**2 for y in sol.y]
plt.plot(times, excited_state_pop)

By zooming in on the above graph we can observe small ripples in the solution:

plt.plot(times[230:260], excited_state_pop[230:260])

Using the Solver class with transformations

The example we gave above is the simplest implementation of the Solver class, however you can also apply model transformations using the class as well.

For example, this time we will construct the same Hamiltonian solver as previously, but with the addition of a Rotating Wave Approximation (RWA). This approximation is a common technique used to simplify Hamiltonians by setting high frequency terms (i.e. any terms in your equation representing fast oscillating waves) to zero.

This time, when we initialize our solver, we can add a RWA by setting the rotating_frame argument to be our static_hamiltonian. Then we set a value for the rwa_cutoff_freq parameter, meaning any terms with frequencies above this value will be set to 0.

rwa_hamiltonian_solver = Solver(static_hamiltonian=drift, 
                                hamiltonian_operators=operators, 
                                hamiltonian_signals=signals, 
                                rotating_frame=drift,
                                rwa_cutoff_freq=2. * w)

Next, we solve the new problem in the same way as before, by calling the .solve() function:

y0 = Statevector([0., 1.])
times = np.linspace(0., T, 500)
rwa_sol = rwa_hamiltonian_solver.solve(t_span=[0., T], y0=y0, t_eval=times)
#plot the results:
rwa_pops = [np.abs(y.data[0])**2 for y in rwa_sol.y]
plt.plot(times, rwa_pops)

As you can see from the differences in the result plots, the effect of the RWA is to approximate the evolution with one that is more smooth. This produces a trade-off between solving time and accuracy: the smoother evolution is simpler and the solver therefore takes less time to come to a solution, although the ultimate result is slightly less accurate. Depending on the application, this discrepancy between the un-approximated and RWA approach may or may not be within acceptable tolerances.

For more details on the Rotating Wave Approximation, see the documentation here

Optimize Calculations using JAX

Aside from the core functionality, one of the most exciting features prioritized in Qiskit Dynamic’s first release is the capability for users to switch from the default backend numerical methods (that use numpy and scipy) to a higher performance JAX implementation internally.

Switching to the JAX backend enables Qiskit Dynamics code to be:

• Just-in-time compiled for faster execution,
• Automatically differentiated (a key component for optimization applications like pulse engineering),
• and executed on GPUs, taking advantage of accelerated linear algebraic computations.

We won’t go into detail on the specifics of using Qiskit Dynamics with JAX but you can check out the JAX user guide for more info, as well as the tutorial on gradient-based pulse optimization.

What’s coming next?

This first release has mainly focused on establishing core functionality for the Qiskit Dynamics module, laying the groundwork for more exciting things to come. The roadmap is still being set for the future of the package, but we couldn’t resist giving you a little teaser of some of the most exciting ideas and plans we have! Keep an eye out for these upcoming themes in future releases:

• Pulse Simulation — expanded integration of Qiskit Pulse, for simulation of experiments from Qiskit Experiments (previously Ignis).
• Model building tools — instead of needing to fully specify operators and matrices for your model, future releases of Qiskit Dynamics aim to include tools to help users describe models more efficiently.
• Matrix Product Operator (MPO) simulation — integration of a new simulator for quantum dynamics, based on matrix product operators, developed in collaboration with Grégoire Misguich, Institut de Physique Théorique, CEA — Saclay (France).
• Perturbation theory computation — new functions for computing Dyson or Magnus series expansions, and solvers based on these.

If you would like to find out more about what is happening with Qiskit Dynamics, and even get involved in it’s development you can check out the following resources:

• Qiskit Dynamics Tutorials page — see example applications of Qiskit Dynamics for simulating Rabi oscillations, simulating a spin chain with noise, simulating a Qiskit Pulse schedule, and gradient optimization of a pulse sequence.
• Qiskit Dynamics GitHub repo — here you can find issues that are being worked on, submit bug reports or feature requests, and contribute to the package itself
• #qiskit-dev slack channel — got a question about the development/roadmap of Qiskit Dynamics, you can ask here!
• Qiskit Community Feedback repo — use this repo to find out about past and future Demo Days (for all qiskit modules, including Dynamics), and get involved in the discussion forum if you’d like to start a topic.