DysonSolver

class DysonSolver(operators, rotating_frame, dt, carrier_freqs, chebyshev_orders, expansion_order=None, expansion_labels=None, integration_method=None, include_imag=None, **kwargs)

Bases: _PerturbativeSolver

Solver for linear matrix differential equations based on the Dyson series.

This class implements the Dyson-series based solver presented in [[1]], which is a variant of the Dysolve algorithm originally introduced in [[2]]. This solver applies to linear matrix differential equations with generators of the form:

\[G(t) = G_0 + \sum_{j=1}^s \textnormal{Re}[f_j(t) e^{i 2 \pi \nu_j t}]G_j,\]

and solves the LMDE in the rotating frame of \(G_0\), which is assumed to be anti-Hermitian. I.e. it solves the LMDE with generator:

\[\tilde{G}(t) = \sum_{j=1}^s \textnormal{Re}[f_j(t) e^{i 2 \pi \nu_j t}]\tilde{G}_j(t),\]

with \(\tilde{G}_i(t) = e^{-t G_0} G_i e^{tG_0}\). The solver is fixed-step, with step size \(\Delta t\) being defined at instantiation, and solves over each step by computing a truncated Dyson series. See the Time-dependent perturbation theory and multi-variable series expansions review for a description of the Dyson series.

At instantiation, the following parameters, which define the structure of the Dyson series used, are fixed:

• The step size \(\Delta t\),

• The operator structure \(G_0\), \(G_i\),

• The reference frequencies \(\nu_j\),

• Approximation schemes for the envelopes \(f_j\) over each time step (see below), and

• The Dyson series terms to keep in the truncation.

A ‘compilation’ or ‘pre-computation’ step computing the truncated expansion occurs at instantiation. Once instantiated, the LMDE can be solved repeatedly for different lists of envelopes \(f_1(t), \dots, f_s(t)\) by calling the solve() method with the initial time t0 and number of time-steps n_steps of size \(\Delta t\). The list of envelopes are specified as Signal objects, whose carrier frequencies are automatically shifted to the reference frequencies \(\nu_j\), with the frequency difference, and any phase, being absorbed into the envelope.

More explicitly, the process of solving over an interval \([t_0, t_0 + \Delta t]\) is as follows. After shifting the carrier frequencies, the resulting envelopes are approximated using a discrete Chebyshev transform, whose orders for each signal is given by chebyshev_orders. That is, for \(t \in [t_0, t_0 + \Delta t]\), each envelope is approximated as:

\[f_j(t) \approx \sum_{m=0}^{d_j} f_{j,m}T_m(t-t_0)\]

where \(T_m(\cdot)\) are the Chebyshev polynomials over the interval, \(f_{j,m}\) are the approximation coefficients attained via Discrete Chebyshev Transform, and \(d_j\) is the order of approximation used for the given envelope. Using:

\[\begin{split}\textnormal{Re}[f_{j,m}T_m(t-t_0)e^{i2 \pi \nu_j t}] = \textnormal{Re}[f_{j,m}e^{i 2 \pi \nu_j t_0}] \cos(2 \pi \nu_j (t-t_0))T_m(t-t_0) \\ + \textnormal{Im}[f_{j,m}e^{i 2 \pi \nu_j t_0}] \sin(-2 \pi \nu_j (t-t_0))T_m(t-t_0)\end{split}\]

The generator is approximately decomposed as

\[\begin{split}\tilde{G}(t) \approx \sum_{j=1}^s \sum_{m=0}^{d_j} \textnormal{Re}[f_{j,m}e^{i 2 \pi \nu_j t_0}] \cos(2 \pi \nu_j (t-t_0))T_m(t-t_0) \tilde{G}_j \\ + \sum_{j=1}^s \sum_{m=0}^{d_j} \textnormal{Im}[f_{j,m}e^{i 2 \pi \nu_j t_0}] \sin(- 2 \pi \nu_j (t-t_0))T_m(t-t_0) \tilde{G}_j\end{split}\]

The multivariable Dyson series is then computed relative to the above decomposition, with the variables being the \(\textnormal{Re}[f_{j,m}e^{i 2 \pi \nu_j t_0}]\) and \(\textnormal{Im}[f_{j,m}e^{i 2 \pi \nu_j t_0}]\), and the operators being \(\cos(2 \pi \nu_j (t-t_0))T_m(t-t_0) G_j\) and \(\sin(- 2 \pi \nu_j (t-t_0))T_m(t-t_0) G_j\). As shown in [[1], [2]], the multivariable Dyson series for intervals of length \(\Delta t\) with different starting times are related via a simple frame change, and as such these need only be computed once, and this makes up the ‘pre-computation’ step of this object.

The optional argument include_imag can be used to control, on a signal by signal basis, whether or not the imaginary terms

\[\textnormal{Im}[f_{j,m}e^{i 2 \pi \nu_j t_0}] \sin(- 2 \pi \nu_j (t-t_0))T_m(t-t_0) \tilde{G}_j\]

are included in the scheme. In generality they are required, but in special cases they are not necessary, such as when \(\nu_j = 0\), or if \(f_j(t)\), including phase, is purely real. By default all such terms are included.

References

[1] (1,2)

D. Puzzuoli, S. F. Lin, M. Malekakhlagh, E. Pritchett, B. Rosand, and C. J. Wood. Algorithms for perturbative analysis and simulation of quantum dynamics. Journal of Computational Physics, 2023. doi:10.1016/j.jcp.2023.112262.

[2] (1,2)

R. Shillito, J. A. Gross, A. Di Paolo, É. Genois, and A. Blais. Fast and differentiable simulation of driven quantum systems. Physical Review Research, 3(3):033266, 2021. Publisher: American Physical Society. doi:10.1103/PhysRevResearch.3.033266.

Initialize.

Parameters:

• operators (List[Operator]) – List of constant operators specifying the operators with signal coefficients.

• rotating_frame (Union[ndarray, number, int, float, complex, Tracer, Array, Array, spmatrix, BCOO, list, Operator, RotatingFrame, None]) – Rotating frame to setup the solver in. Must be Hermitian or anti-Hermitian.

• dt (float) – Fixed step size to compile to.

• carrier_freqs (Union[ndarray, number, int, float, complex, Tracer, Array, Array, spmatrix, BCOO, list]) – Carrier frequencies of the signals in the generator decomposition.

• chebyshev_orders (List[int]) – Approximation degrees for each signal over the interval [0, dt].

• expansion_order (Optional[int]) – Order of perturbation terms to compute up to. Specifying this argument results in computation of all terms up to the given order. Can be used in conjunction with expansion_terms.

• expansion_labels (Optional[List[Multiset]]) – Specific perturbation terms to compute. If both expansion_order and expansion_terms are specified, then all terms up to expansion_order are computed, along with the additional terms specified in expansion_terms. Labels are specified either as Multiset or as valid arguments to the Multiset constructor. This function further requires that Multisets consist only of non-negative integers.

• integration_method (Optional[str]) – ODE solver method to use when computing perturbation terms.

• include_imag (Optional[List[bool]]) – List of bools determining whether to keep imaginary components in the signal approximation. Defaults to True for all signals.

• kwargs – Additional arguments to pass to the solver when computing perturbation terms.

Methods

solve(t0, n_steps, y0, signals, jax_control_flow=None)

Solve given an initial time, number of steps, signals, and initial state.

Note that this method can be used to solve a list of simulations at once, by specifying one or more of the arguments t0, n_steps, y0, or signals as a list of valid inputs. For this mode of operation, all of these arguments must be either lists of the same length, or a single valid input, which will be used repeatedly.

Parameters:

• t0 (Union[float, List[float]]) – Initial time.

• n_steps (Union[int, List[int]]) – Number of time steps to solve for.

• y0 (Union[ndarray, number, int, float, complex, Tracer, Array, Array, spmatrix, BCOO, list, List[Union[ndarray, number, int, float, complex, Tracer, Array, Array, spmatrix, BCOO, list]]]) – Initial state at time t0.

• signals (Union[List[Signal], List[List[Signal]]]) – List of signals.

• jax_control_flow (Optional[bool]) – Whether or not to use JAX control flow during solver runs. If unspecified it will be automatically determined based on the input state type, whether the call is taking place within a JAX transformation, or if the underlying expansion polynomial consists of JAX arrays.

Returns:

Results object, or list of results objects.

Return type:

OdeResult

Raises:

QiskitError – If improperly formatted arguments.

Attributes

model

Model object storing expansion details.