density_matrix_evolution_pwc

Graph.density_matrix_evolution_pwc(initial_density_matrix, hamiltonian, lindblad_terms, sample_times=None, error_tolerance=1e-06, *, name=None)

Calculate the state evolution of an open system described by the GKS–Lindblad master equation.

The controls that you provide to this function have to be in piecewise-constant format. If your controls are smooth sampleable tensor-valued functions (STFs), you have to discretize them with discretize_stf before passing them to this function. You may need to increase the number of segments that you choose for the discretization depending on the sizes of oscillations in the smooth controls.

By default, this function computes an approximate piecewise-constant solution for the consideration of efficiency, with the accuracy controlled by the parameter error_tolerance. If your system is small, you can set the error_tolerance to None to obtain an exact solution.

Note that when using the exact method, both hamiltonian and lindblad_terms are converted to the dense representation, regardless of their original formats. This means that the computation can be slow and memory intensive when applied to large systems.

When using the approximate method, the sparse representation is used internally if hamiltonian is a SparsePwc, otherwise the dense representation is used.

Parameters:

initial_density_matrix (np.ndarray or Tensor) – A 2D array of the shape (D, D) representing the initial density matrix of the system, \(\rho_{\rm s}\). You can also pass a batch of density matrices and the input data shape must be (B, D, D) where B is the batch dimension.
hamiltonian (Pwc or SparsePwc) – A piecewise-constant function representing the effective system Hamiltonian, \(H_{\rm s}(t)\), for the entire evolution duration.
lindblad_terms (list[tuple[float, np.ndarray or Tensor or scipy.sparse.spmatrix]]) – A list of pairs, \((\gamma_j, L_j)\), representing the positive decay rate \(\gamma_j\) and the Lindblad operator \(L_j\) for each coupling channel \(j\). You must provide at least one Lindblad term.
sample_times (np.ndarray or None, optional) – A 1D array like object of length \(T\) specifying the times \(\{t_i\}\) at which this function calculates system states. Must be ordered and contain at least one element. Note that increasing the density of sample times does not affect the computation precision of this function. If omitted only the evolved density matrix at the final time of the system Hamiltonian is returned.
error_tolerance (float or None, optional) – Defaults to 1e-6. This option enables an approximate method to solve the master equation, meaning the 2-norm of the difference between the propagated state and the exact solution at the final time (and at each sample time if passed) is within the error tolerance. Note that, if set, this value must be smaller than 1e-2 (inclusive). However, setting it to a too small value (for example below 1e-12) might result in slower computation, but would not further improve the precision, since the dominating error in that case is due to floating point error. You can also explicitly set this option to None to find the exact piecewise-constant solution. Note that using the exact solution can be very computationally demanding in calculations involving a large Hilbert space or a large number of segments.
name (str or None, optional) – The name of the node.

Returns:

The time-evolved density matrix, with shape (D, D) or (T, D, D), depending on whether you provided sample times. If you provide a batch of initial states, the shape is (B, T, D, D) or (B, D, D).

Return type:

Tensor

See also

Graph.jump_trajectory_evolution_pwc: Trajectory-based state evolution of an open quantum system.
Graph.steady_state: Compute the steady state of open quantum system.

Notes

Under the Markovian approximation, the dynamics of an open quantum system can be described by the GKS–Lindblad master equation [1] [2]

\[\frac{{\rm d}\rho_{\rm s}(t)}{{\rm d}t} = -i [H_{\rm s}(t), \rho_{\rm s}(t)] + \sum_j \gamma_j {\mathcal D}[L_j] \rho_{\rm s}(t) ,\]

where \({\mathcal D}\) is a superoperator describing the decoherent process in the system evolution and defined as

\[{\mathcal D}[X]\rho := X \rho X^\dagger - \frac{1}{2}\left( X^\dagger X \rho + \rho X^\dagger X \right)\]

for any system operator \(X\).

This function uses sparse matrix multiplication when the Hamiltonian is passed as a SparsePwc and the Lindblad operators as sparse matrices. This leads to more efficient calculations when they involve large operators that are relatively sparse (contain mostly zeros). In this case, the initial density matrix is still a densely represented array or tensor.

References

Examples

Simulate a trivial decay process for a single qubit described by the following master equation \(\dot{\rho} = -i[\sigma_z / 2, \, \rho] + \mathcal{D}[\sigma_-]\rho\).

>>> duration = 20
>>> initial_density_matrix = np.array([[0, 0], [0, 1]])
>>> hamiltonian = graph.constant_pwc_operator(
...     duration=duration, operator=graph.pauli_matrix("Z") / 2
... )
>>> lindblad_terms = [(1, graph.pauli_matrix("M"))]
>>> graph.density_matrix_evolution_pwc(
...     initial_density_matrix, hamiltonian, lindblad_terms, name="decay"
... )
<Tensor: name="decay", operation_name="density_matrix_evolution_pwc", shape=(2, 2)>
>>> result = bo.execute_graph(graph=graph, output_node_names="decay")
>>> result["output"]["decay"]["value"]
array([[9.99999998e-01+0.j, 0.00000000e+00+0.j],
       [0.00000000e+00+0.j, 2.06115362e-09+0.j]])

See more examples in the How to simulate open system dynamics and How to simulate large open system dynamics user guides.