# Source code for oqupy.system

# Copyright 2022 The TEMPO Collaboration
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
"""
Module on physical information of the system.
"""

from typing import Callable, List, Optional, Text, Tuple
from copy import copy
from functools import lru_cache

import numpy as np
from numpy import ndarray

from scipy.linalg import expm
from scipy import integrate

from oqupy.base_api import BaseAPIClass
from oqupy.config import NpDtype
import oqupy.operators as opr

[docs]class BaseSystem(BaseAPIClass):
"""Base class for systems. """
def __init__(
self,
dimension: int,
name: Optional[Text] = None,
description: Optional[Text] = None) -> None:
"""Create a BaseSystem object."""
self._dimension = dimension
super().__init__(name, description)

@property
def dimension(self) -> ndarray:
"""Hilbert space dimension of the system. """
return self._dimension

[docs]class System(BaseSystem):
r"""
Represents a system (without any coupling to a non-Markovian bath).
It is possible to include Lindblad terms in the master equation.
The equations of motion for a system density matrix (without any coupling
to a non-Markovian bath) is then:

.. math::

\frac{d}{dt}\rho(t) = &-i [\hat{H}, \rho(t)] \\
&+ \sum_n^N \gamma_n \left(
\hat{A}_n \rho(t) \hat{A}_n^\dagger
- \frac{1}{2} \hat{A}_n^\dagger \hat{A}_n \rho(t)
- \frac{1}{2} \rho(t) \hat{A}_n^\dagger \hat{A}_n \right)

with hamiltionian :math:\hat{H}, the rates gammas :math:\gamma_n and
linblad_operators :math:\hat{A}_n.

Parameters
----------
hamiltonian: ndarray
System-only Hamiltonian :math:\hat{H}.
gammas: List(float)
The rates :math:\gamma_n.
The Lindblad operators :math:\hat{A}_n.
name: str
An optional name for the system.
description: str
An optional description of the system.
"""

def __init__(
self,
hamiltonian: ndarray,
gammas: Optional[List[float]] = None,
name: Optional[Text] = None,
description: Optional[Text] = None) -> None:
"""Create a System object. """

# input check for Hamiltonian.
self._hamiltonian = _check_hamiltonian(hamiltonian)
tmp_dimension = self._hamiltonian.shape[0]
# input check gammas and lindblad_operators
super().__init__(tmp_dimension, name, description)

[docs]    @lru_cache(4)
def liouvillian(self) -> ndarray:
r"""
Returns the Liouvillian super-operator :math:\mathcal{L} with

.. math::

\mathcal{L}\rho = -i [\hat{H}, \rho]
+ \sum_n^N \gamma_n \left(
\hat{A}_n \rho \hat{A}_n^\dagger
- \frac{1}{2} \hat{A}_n^\dagger \hat{A}_n \rho
- \frac{1}{2} \rho \hat{A}_n^\dagger \hat{A}_n
\right) .

Returns
-------
liouvillian : ndarray
Liouvillian :math:\mathcal{L}.
"""
return _liouvillian(self._hamiltonian,
self._gammas,

[docs]    def get_propagators(self, dt, start_time, subdiv_limit, epsrel):
"""Prepare propagator functions for the system. """
first_step = expm(self.liouvillian()*dt/2.0)
second_step = expm(self.liouvillian()*dt/2.0)
def propagators(step: int):
"""Create the system propagators (first and second half) for
the time step step  """
return first_step, second_step
return propagators

@property
def hamiltonian(self) -> ndarray:
"""The system Hamiltonian."""
return copy(self._hamiltonian)

@property
def gammas(self) -> List[float]:
"""List of gammas."""
return copy(self._gammas)

@property

[docs]class TimeDependentSystem(BaseSystem):
r"""
Represents an explicitly time dependent system (without any coupling to a
non-Markovian bath). It is possible to include (also explicitly
time dependent) Lindblad terms in the master equation.
The equations of motion for a system density matrix (without any coupling
to a non-Markovian bath) is then:

.. math::

\frac{d}{dt}\rho(t) = &-i [\hat{H}(t), \rho(t)] \\
&+ \sum_n^N \gamma_n(t) \left(
\hat{A}_n(t) \rho(t) \hat{A}_n(t)^\dagger
- \frac{1}{2} \hat{A}_n^\dagger(t) \hat{A}_n(t) \rho(t)
- \frac{1}{2} \rho(t) \hat{A}_n^\dagger(t) \hat{A}_n(t) \right)

with the time dependent hamiltionian :math:\hat{H}(t), the  time
dependent rates gammas :math:\gamma_n(t) and the time dependent
linblad_operators :math:\hat{A}_n(t).

Parameters
----------
hamiltonian: callable
System-only Hamiltonian :math:\hat{H}(t).
gammas: List(callable)
The rates :math:\gamma_n(t).
The Lindblad operators :math:\hat{A}_n(t).
name: str
An optional name for the system.
description: str
An optional description of the system.
"""
def __init__(
self,
hamiltonian: Callable[[float], ndarray],
gammas: \
Optional[List[Callable[[float], float]]] = None,
Optional[List[Callable[[float], ndarray]]] = None,
name: Optional[Text] = None,
description: Optional[Text] = None) -> None:
"""Create a TimeDependentSystem object."""
# input check for Hamiltonian.
self._hamiltonian = _check_tdependent_hamiltonian(hamiltonian)
tmp_dimension = self._hamiltonian(1.0).shape[0]
# input check gammas and lindblad_operators
gammas,

super().__init__(tmp_dimension, name, description)

[docs]    def liouvillian(self, t: float) -> ndarray:
r"""
Returns the Liouvillian super-operator :math:\mathcal{L}(t) with

.. math::

\mathcal{L}(t)\rho = -i [\hat{H}(t), \rho]
+ \sum_n^N \gamma_n \left(
\hat{A}_n(t) \rho \hat{A}_n^\dagger(t)
- \frac{1}{2} \hat{A}_n^\dagger(t) \hat{A}_n(t) \rho
- \frac{1}{2} \rho \hat{A}_n^\dagger(t) \hat{A}_n(t)
\right),

with time :math:t.

Parameters
----------
t: float (default = None)
time :math:t.

Returns
-------
liouvillian : ndarray
Liouvillian :math:\mathcal{L}(t) at time :math:t.
"""
hamiltonian = self._hamiltonian(t)
gammas = [gamma(t) for gamma in self._gammas]

[docs]    def get_propagators(self, dt, start_time, subdiv_limit, epsrel):
"""Prepare propagator functions for the system according to
subdiv_limit. """
if subdiv_limit is None:
# Sample Liouvillian at dt/4, 3dt/4 to make propagators for first-
# and second-half timesteps
def propagators(step: int):
"""Create the system propagators (first and second half) for
the time step step  """
t = start_time + step * dt
first_step = expm(self.liouvillian(t+dt/4.0)*dt/2.0)
second_step = expm(self.liouvillian(t+dt*3.0/4.0)*dt/2.0)
return first_step, second_step
else:
# Integrate Liouvillian to make propagators for first- and
# second-half timesteps
def propagators(step: int):
"""Create the system propagators (first and second half) for
the time step step  """
t = start_time + step * dt
a=t,
b=t+dt/2.0,
epsrel=epsrel,
limit=subdiv_limit)[0])
a=t+dt/2.0,
b=t+dt,
epsrel=epsrel,
limit=subdiv_limit)[0])
return first_step, second_step
return propagators

@property
def hamiltonian(self) -> Callable[[float], ndarray]:
"""The system Hamiltonian. """
return copy(self._hamiltonian)

@property
def gammas(self) -> List[Callable[[float], float]]:
"""List of gammas. """
return copy(self._gammas)

@property

[docs]class TimeDependentSystemWithField(BaseSystem):
r"""
Represents a system which depends on time and an auxiliary field
(complex scalar). Forms one component of a MeanFieldSystem.

It is possible to include time (but not field) dependent Lindblad
terms in the master equation. The equations of motion for the system
density matrix (without any coupling to a non-Markovian bath) is
then:

.. math::

\frac{d}{dt}\rho(t) = &-i [\hat{H}(t, \langle a \rangle), \rho(t)] \\
&+ \sum_n^N \gamma_n(t) \left(
\hat{A}_n(t) \rho(t) \hat{A}_n(t)^\dagger
- \frac{1}{2} \hat{A}_n^\dagger(t) \hat{A}_n(t) \rho(t)
- \frac{1}{2} \rho(t) \hat{A}_n^\dagger(t) \hat{A}_n(t) \right)

with the  hamiltionian :math:\hat{H}(t, \langle a \rangle)
depending on both time :math:t and field :math:\langle
a \rangle, the  time dependent rates gammas
:math:\gamma_n(t) and the time dependent linblad_operators
:math:\hat{A}_n(t).

Parameters
----------
hamiltonian: callable
System-only Hamiltonian :math:\hat{H}(t, \langle a \rangle)
where :math:\langle a \rangle is the field at time :math:t.
gammas: list(callable)
The rates :math:\gamma_n(t).
The Lindblad operators :math:\hat{A}_n(t).
name: str
An optional name for the system.
description: str
An optional description of the system.
"""
def __init__(
self,
hamiltonian: Callable[[float, complex], ndarray],
gammas: \
Optional[List[Callable[[float], float]]] = None,
Optional[List[Callable[[float], ndarray]]] = None,
name: Optional[Text] = None,
description: Optional[Text] = None) -> None:
"""Create a TimeDependentSystemWithField object."""

# input check for Hamiltonian
self._hamiltonian = _check_tfielddependent_hamiltonian(hamiltonian)
tmp_dimension = self._hamiltonian(1.0, 1.0+1.0j).shape[0]

# input check gammas and lindblad_operators
gammas,

super().__init__(tmp_dimension, name, description)

def _linearised_hamiltonian(self, t0: float, t: float,
field: complex,
field_derivative: complex) -> complex:
r"""
Return value of the system Hamiltonian at time t using a linearisation
of the field coupled to the subsystem from its value at time t0.
"""
return self._hamiltonian(t,
self._linearised_field(t0, t, field, field_derivative))

@staticmethod
def _linearised_field(t0: float, t: float,
field: complex,
field_derivative: complex):
r"""
Return the value of the field at time (t-t0) given the value at t0
in a linear approximation using the value of the time derivative at
t0.
"""
return field + field_derivative * (t-t0)

[docs]    def liouvillian(self,
t0: float,
t: float,
field: complex,
field_derivative: complex) -> ndarray:
r"""
Returns the Liouvillian super-operator
:math:\mathcal{L}(t, \langle a \rangle) such that

.. math::

\mathcal{L}(t, \langle a \rangle)\rho
= -i [\hat{H}(t, \langle a \rangle), \rho]
+ \sum_n^N \gamma_n \left(
\hat{A}_n(t) \rho \hat{A}_n^\dagger(t)
- \frac{1}{2} \hat{A}_n^\dagger(t) \hat{A}_n(t) \rho
- \frac{1}{2} \rho \hat{A}_n^\dagger(t) \hat{A}_n(t)
\right),

with time :math:t.

Parameters
----------
t0: float
Start time of the current step.
t: float
Current time :math:t.
field: complex
Field value at time :math:t obtained from the
linearisation of the field at :math:t using the field
equation of motion.
field_derivative: complex
Value of the time derivative of the field at time t0

Returns
-------
liouvillian : ndarray
Liouvillian :math:\mathcal{L}(t, \langle a \rangle) at time
:math:t using a linearisation of the field \langle a \rangle
from its value at t0 to time t.
"""
try:
t0 = float(t0)
except Exception as e:
raise TypeError("Argument t0 must be float") from e
try:
t = float(t)
except Exception as e:
raise TypeError("Argument t must be float") from e
assert t >= t0, "Argument t must equal or exceed t0"
try:
field = complex(field)
except Exception as e:
raise TypeError("Argument field must be complex") from e
try:
field_derivative = complex(field_derivative)
except Exception as e:
raise TypeError("Argument field_derivative must be complex") from e
hamiltonian = self._linearised_hamiltonian(t0, t, field,
field_derivative)
gammas = [gamma(t) for gamma in self._gammas]

[docs]    def get_propagators(self, dt, start_time, subdiv_limit, epsrel):
"""Prepare propagator functions for the system according to
subdiv_limit. """
if subdiv_limit is None:
# Sample Liouvillian at dt/4, 3dt/4 to make propagators for first-
# and second-half timesteps
def propagators(step: int, field: complex,
field_derivative: complex):
t = start_time + step * dt
first_step = expm(self.liouvillian(t, t+dt/4.0,
field, field_derivative)*dt/2.0)
second_step = expm(self.liouvillian(t, t+dt*3.0/4.0,
field, field_derivative)*dt/2.0)
return first_step, second_step
else:
# Integrate Liouvillian to make propagators for first- and
# second-half timesteps
def propagators(step: int, field: complex,
field_derivative: complex):
t = start_time + step * dt
liouvillian = lambda tau: self.liouvillian(t, tau,
field, field_derivative)
a=t,
b=t+dt/2.0,
epsrel=epsrel,
limit=subdiv_limit)[0])
a=t+dt/2.0,
b=t+dt,
epsrel=epsrel,
limit=subdiv_limit)[0])
return first_step, second_step
return propagators

@property
def hamiltonian(self) -> Callable[[float, complex], ndarray]:
"""The system Hamiltonian. """
return copy(self._hamiltonian)

@property
def gammas(self) -> List[Callable[[float], float]]:
"""List of gammas. """
return copy(self._gammas)

@property

[docs]class MeanFieldSystem(BaseAPIClass):
r"""Represents a collection of time dependent systems interacting
with a common field. The systems are encoded as
TimeDependentSystemWithField objects, and the field as a complex
scalar :math:\langle a \rangle that evolves according to a
specified equation of motion :math:\partial_t\langle a \rangle.

Parameters
----------
system_list: List[TimeDependentSystemWithField],
List of TimeDependentSystemWithField objects interacting with
a common field :math:\langle a \rangle.
field_eom: callable
Field equation of motion :math:\partial_t
\langle a \rangle(t, [\rho], \langle a \rangle)
where :math:[\rho] is a list of square matrices for the state
of each system in system_list at time :math:t and
:math:\langle a \rangle the field at time :math:t.
name: str
An optional name for the mean-field system.
description: str
An optional description of the mean-field system.
"""

def __init__(self,
system_list: List[TimeDependentSystemWithField],
field_eom: Callable[[float, List[ndarray], complex], complex],
name: Optional[Text] = None,
description: Optional[Text] = None) -> None:

super().__init__(name, description)
tmp_system_list = _check_mean_field_system_list(system_list)
self._system_list = tmp_system_list

# input check for field equation of motion
tmp_dimension_list = [system.hamiltonian(1.0, 1.0+1.0j).shape[0]
for system in self.system_list]
tmp_field_eom = _check_mean_field_system_eom(tmp_dimension_list,
field_eom)
self._field_eom = tmp_field_eom

@property
def system_list(self) -> List[TimeDependentSystemWithField]:
"""The list of systems interacting with a common field. """
return self._system_list

@property
def field_eom(self) -> Callable[[float, List[ndarray], complex], complex]:
"""The field equation of motion. """
return copy(self._field_eom)

[docs]class SystemChain(BaseAPIClass):
"""
Represents a 1D chain of systems with nearest neighbor interactions.

Parameters
----------
hilbert_space_dimensions: List[int]
Hilbert space dimension for each chain site.
name: str
An optional name for the system chain.
description: str
An optional description of the system chain.
"""
def __init__(
self,
hilbert_space_dimensions: List[int],
name: Optional[Text] = None,
description: Optional[Text] = None) -> None:
"""Create a SystemChain object. """
tmp_hs_dims = np.array(hilbert_space_dimensions, int)
assert len(tmp_hs_dims.shape) == 1
assert len(hilbert_space_dimensions) >= 1
assert np.all(tmp_hs_dims > 0)
self._hs_dims = tmp_hs_dims

self._site_liouvillians = []
for hs_dim in self._hs_dims:
self._site_liouvillians.append(
np.zeros((hs_dim**2, hs_dim**2), dtype=NpDtype))

self._nn_liouvillians = []
for hs_dim_l, hs_dim_r in zip(self._hs_dims[:-1], self._hs_dims[1:]):
self._nn_liouvillians.append(
np.zeros((hs_dim_l**2 * hs_dim_r**2, hs_dim_l**2 * hs_dim_r**2),
dtype=NpDtype))

super().__init__(name, description)

def __len__(self):
"""Chain length. """
return len(self._hs_dims)

@property
def hs_dims(self):
"""Hilbert space dimension for each chain site. """
return self._hs_dims

@property
def site_liouvillians(self):
"""The single site Liouvillians. """
return self._site_liouvillians

@property
def nn_liouvillians(self):
"""The nearest neighbor Liouvillians. """
return self._nn_liouvillians

self,
site: int,
hamiltonian: ndarray) -> None:
r"""
Add a hamiltonian term to a single site Liouvillian

.. math::

\mathcal{L} \rho_n = -i [\hat{H}, \rho_n]

with site :math:n and hamiltonian :math:\hat{H}.

Parameters
----------
site: int
Index of the site.
hamiltonian: ndarray
Hamiltonian acting on the single site.
"""
assert isinstance(site, int)
assert site >= 0
assert site < len(self)
op = np.array(hamiltonian, dtype=NpDtype)
assert len(op.shape) == 2
assert op.shape[0] == op.shape[1]
assert self._hs_dims[site] == op.shape[0]

self._site_liouvillians[site] += (0.0-1.0j) * opr.commutator(op)

self,
site: int,
liouvillian: ndarray) -> None:
"""

Parameters
----------
site: int
Index of the site.
liouvillian: ndarray
Liouvillian acting on the single site.
"""
self._site_liouvillians[site] += np.array(liouvillian, dtype=NpDtype)

self,
site: int,
gamma: Optional[float] = 1.0) -> None:
r"""

.. math::

\mathcal{L} \rho_n = \gamma \left(
\hat{A} \rho_n \hat{A}^\dagger
- \frac{1}{2} \hat{A}^\dagger \hat{A} \rho_n
- \frac{1}{2} \rho_n \hat{A}^\dagger \hat{A} \right)

with site :math:n, lindblad_operator :math:\hat{A},
and gamma :math:\gamma.

Parameters
----------
site: int
Index of the site.
Lindblad dissipator acting on the single site.
gamma: float
Optional multiplicative factor :math:\gamma.
"""
op_dagger = op.conjugate().T
self._site_liouvillians[site] += \
gamma * (opr.left_right_super(op, op_dagger) \
- 0.5 * opr.acommutator(np.dot(op_dagger, op)))

self,
site: int,
hamiltonian_l: ndarray,
hamiltonian_r: ndarray) -> None:
r"""
Add a hamiltonian term to the Liouvillian of two neighboring sites:

.. math::

\mathcal{L} \rho_{n,n+1} =
-i [\hat{H}_l \otimes \hat{H}_r, \rho_{n,n+1}]

with site :math:n, hamiltonian_l :math:\hat{H}_l and
hamiltonian_r :math:\hat{H}_r.

Parameters
----------
site: int
Index of the left site :math:n.
hamiltonian_l: ndarray
Hamiltonian acting on the left site :math:n.
hamiltonian_r: ndarray
Hamiltonian acting on the right site :math:n+1.
"""
assert isinstance(site, int)
assert site >= 0
assert site < len(self) - 1
op_l = np.array(hamiltonian_l, dtype=NpDtype)
op_r = np.array(hamiltonian_r, dtype=NpDtype)
assert len(op_l.shape) == 2
assert len(op_r.shape) == 2
assert op_l.shape[0] == op_l.shape[1]
assert op_r.shape[0] == op_r.shape[1]
assert self._hs_dims[site] == op_l.shape[0]
assert self._hs_dims[site+1] == op_r.shape[0]

self._nn_liouvillians[site] += (0.0-1.0j) \
* opr.cross_commutator(op_l, op_r)

self,
site: int,
liouvillian_l_r: ndarray) -> None:
"""
Add Liouvillian of for the two neighboring sites site and site +1.

Parameters
----------
site: int
Index of the left site :math:n.
liouvillian_l_r: ndarray
Liouvillian acting on sites :math:n and :math:n+1.
"""
self._nn_liouvillians[site] += np.array(liouvillian_l_r, dtype=NpDtype)

self,
site: int,
gamma: Optional[float] = 1.0) -> None:
r"""

.. math::

\mathcal{L} \rho_{n,n+1} = \gamma \left(
\hat{A} \rho_{n,n+1} \hat{A}^\dagger
- \frac{1}{2} \hat{A}^\dagger \hat{A} \rho_{n,n+1}
- \frac{1}{2} \rho_{n,n+1} \hat{A}^\dagger \hat{A} \right)

where :math:\hat{A}=\hat{A}_l\otimes\hat{A}_r, with site :math:n,
lindblad_operator_l :math:\hat{A}_l,
lindblad_operator_r :math:\hat{A}_r, and gamma :math:\gamma.

Parameters
----------
site: int
Index of the left site :math:n.
Lindblad dissipator acting on the left site :math:n.
Lindblad dissipator acting on the right site :math:n+1.
gamma: float
Optional multiplicative factor :math:\gamma.
"""
assert isinstance(site, int)
assert site >= 0
assert site < len(self) - 1
assert len(op_l.shape) == 2
assert len(op_r.shape) == 2
assert op_l.shape[0] == op_l.shape[1]
assert op_r.shape[0] == op_r.shape[1]
assert self._hs_dims[site] == op_l.shape[0]
assert self._hs_dims[site+1] == op_r.shape[0]

cross_lr = opr.cross_left_right_super(
operator_1_l=op_l,
operator_1_r=op_l.T.conjugate(),
operator_2_l=op_r,
operator_2_r=op_r.T.conjugate())
cross_acomm = opr.cross_acommutator(
operator_1=op_l.T.conjugate() @ op_l,
operator_2=op_r.T.conjugate() @ op_r)

self._nn_liouvillians[site] += \
gamma * (cross_lr - 0.5 * cross_acomm)

[docs]    def get_nn_full_liouvillians(self) -> List[ndarray]:
"""
Return the list of nearest neighbor Liouvillians
(incorporating single site terms).
"""
assert len(self) >= 2, \
"To return a full set of nearest neighbor liouvillians, " \
+ "the chain has to be at least two sites long."

nn_full_liouvillians = []
for i in range(len(self)-1):
factor_l = 1 if i == 0 else 0.5
factor_r = 1 if i == len(self)-2 else 0.5

liouv_l = self._site_liouvillians[i]
id_l = np.identity(self._hs_dims[i]**2)
liouv_r = self._site_liouvillians[i+1]
id_r = np.identity(self._hs_dims[i+1]**2)
liouv_nn = self._nn_liouvillians[i]

nn_full_liouvillian = \
factor_l * np.kron(liouv_l, id_r) \
+ factor_r * np.kron(id_l, liouv_r) \
+ liouv_nn

nn_full_liouvillians.append(nn_full_liouvillian)

return nn_full_liouvillians

def _check_hamiltonian(hamiltonian) -> ndarray:
"""Input checking for a single Hamiltonian. """
try:
tmp_hamiltonian = np.array(hamiltonian, dtype=NpDtype)
tmp_hamiltonian.setflags(write=False)
except Exception as e:
raise AssertionError("Coupling operator must be numpy array") from e
assert len(tmp_hamiltonian.shape) == 2, \
"Coupling operator is not a matrix."
assert tmp_hamiltonian.shape[0] == \
tmp_hamiltonian.shape[1], \
"Coupling operator is not a square matrix."
return tmp_hamiltonian

def _check_tdependent_hamiltonian(hamiltonian) -> Callable[[float],
ndarray]:
"""Input checking for a time-dependent Hamiltonian. """
try:
tmp_hamiltonian = np.vectorize(hamiltonian)
_check_hamiltonian(tmp_hamiltonian(1.0))
except Exception as e:
raise AssertionError(
"Time dependent Hamiltonian must be vectorizable callable.") \
from e
return tmp_hamiltonian

def _check_tfielddependent_hamiltonian(hamiltonian) -> Callable[[float,
complex], ndarray]:
try:
tmp_hamiltonian = np.vectorize(hamiltonian)
_check_hamiltonian(tmp_hamiltonian(1.0, 1.0+1.0j))
except Exception as e:
raise AssertionError(
"Time and field dependent Hamiltonian must be vectorizable "\
"callable.") from e
return tmp_hamiltonian

def _check_dissipator_lists(gammas, lindblad_operators) -> Tuple[List, List]:
"""Check gammas and lindblad operators are lists of equal length."""
if gammas is None:
gammas = []
assert isinstance(gammas, list), \
"Argument gammas must be a list)]."
"Argument lindblad_operators must be a list."
"Lists gammas and lindblad_operators must have the same length."

List[float], List[ndarray]]:
"""Input check for time-independent gammas and lindblad_operators"""
# firstly check both are lists of the same length
try:
tmp_gammas = []
for gamma in gammas:
tmp_gammas.append(float(gamma))
except Exception as e:
raise AssertionError("All elements of gammas must be floats.") \
from e
try:
except Exception as e:
raise AssertionError(
"All elements of lindblad_operators must be numpy arrays.") \
from e

gammas,
List[Callable[[float], ndarray]]]:
"""Input check for time-dependent gammas and lindblad_operators"""
# firstly check both are lists of the same length
gammas,
try:
tmp_gammas = []
for gamma in gammas:
float(gamma(1.0))
tmp_gamma = np.vectorize(gamma)
tmp_gammas.append(tmp_gamma)
except Exception as e:
raise AssertionError(
"All elements of gammas must be vectorizable " \
+ "callables returning floats.") from e
try:
except Exception as e:
raise AssertionError(
"All elements of lindblad_operators must be vectorizable " \
+ "callables returning numpy arrays.") from e

def _check_mean_field_system_list(system_list):
assert isinstance(system_list, list), "Parameter system_list must "\
"be a list of TimeDependentSystemWithField objects."
for obj in system_list:
assert isinstance(obj, TimeDependentSystemWithField), "Each "\
"element of system_list must be a "\
"TimeDependentSystemWithField object."
return system_list

def _check_mean_field_system_eom(dim_list, field_eom):
"""Input check a field equation of motion for a mean-field-system"""
test_matrix_list = [_create_density_matrix(dim) for dim in dim_list]
test_field = 1.0+1.0j
test_time = 1.0
try:
value = field_eom(test_time, test_matrix_list, test_field)
complex(value)
except Exception as e:
raise AssertionError("Field equation of motion must "\
"take a time, a list of matrices with shapes\n "\
+ str([f"({dim}, {dim})" for dim in dim_list]) \
+ " and return a complex scalar.") from e
return field_eom

"""
liouvillian = -1j * opr.commutator(hamiltonian)
for gamma, op in zip(gammas, lindblad_operators):
op_dagger = op.conjugate().T
liouvillian += gamma * (opr.left_right_super(op, op_dagger) \
- 0.5 * opr.acommutator(np.dot(op_dagger, op)))
return liouvillian

def _create_density_matrix(dim, seed=1):
r"""Create a repeatable (dim,dim) matrix that represents
a valid density matrix :math:rho"""
rng = np.random.default_rng(seed)
a = rng.random((dim, dim)) + 1j*rng.random((dim,dim))
b = np.matmul(a, a.conj().T)
rho = b / b.trace()
return rho