# Hierarchical Equations of Motion

## HEOM Solvers

This module provides solvers for system-bath evolution using the HEOM (hierarchy equations of motion).

See https://en.wikipedia.org/wiki/Hierarchical_equations_of_motion for a very basic introduction to the technique.

The implementation is derived from the BoFiN library (see https://github.com/tehruhn/bofin) which was itself derived from an earlier implementation in QuTiP itself.

For backwards compatibility with QuTiP 4.6 and below, a new version of HSolverDL (the Drude-Lorentz specific HEOM solver) is provided. It is implemented on top of the new HEOMSolver but should largely be a drop-in replacement for the old HSolverDL.

heomsolve(

  H,
  bath,
  max_depth,
  state0,
  tlist,
  *,
  e_ops=None,
  args=None,
  options=None,

)

  Hierarchical Equations of Motion (HEOM) solver that supports multiple baths.

  Each bath must be either bosonic or fermionic, but bosonic and fermionic baths can be mixed.

  If you need to run many evolutions of the same system and bath, consider using `HEOMSolver` directly to avoid having to continually reconstruct the equation hierarchy for every evolution.

  Parameters:
    H `Qobj`, `QobjEvo`

      Possibly time-dependent system Liouvillian or Hamiltonian as a Qobj or QobjEvo. List of [`Qobj`, `Coefficient`], or callable that can be made into `QobjEvo` are also accepted.

    bath Bath specification or list of Bath specifications

      A bath containing the exponents of the expansion of the bath correlation funcion and their associated coefficients and coupling operators, or a list of baths.

      Each bath can be specified as either an object of type `Bath`, `BosonicBath`, `FermionicBath`, or their subtypes, or as a tuple `(env, Q)`, where `env` is an `ExponentialBosonicEnvironment` or an `ExponentialFermionicEnvironment` and `Q` the system coupling operator.

    max_depth int

      The maximum depth of the heirarchy (i.e. the maximum number of bath exponent “excitations” to retain).

    state0 `Qobj` or `HierarchyADOsState` or array-like

      If `rho0` is a `Qobj` the it is the initial state of the system (i.e. a `Qobj` density matrix).

      If it is a `HierarchyADOsState` or array-like, then `rho0` gives the initial state of all ADOs.

      Usually the state of the ADOs would be determine from a previous call to `.run(...)` with the solver results option `store_ados` set to True. For example, `result = solver.run(...)` could be followed by `solver.run(result.ado_states[-1], tlist)`.

      If a numpy array-like is passed its shape must be `(number_of_ados, n, n)` where `(n, n)` is the system shape (i.e. shape of the system density matrix) and the ADOs must be in the same order as in `.ados.labels`.

    tlist list

      An ordered list of times at which to return the value of the state.

    e_ops Qobj / QobjEvo / callable / list / dict / None, optional

      A list or dictionary of operators as `Qobj`, `QobjEvo` and/or callable functions (they can be mixed) or a single operator or callable function. For an operator `op`, the result will be computed using `(state * op).tr()` and the state at each time `t`. For callable functions, `f`, the result is computed using `f(t, ado_state)`. The values are stored in the `expect` and `e_data` attributes of the result (see the return section below).

    args dict, optional

      Change the `args` of the RHS for the evolution.

    options dict, optional

      Generic solver options.

      - store_final_state : bool
        Whether or not to store the final state of the evolution in the result class.
      - store_states : bool, None
        Whether or not to store the state vectors or density matrices. On None the states will be saved if no expectation operators are given.
      - store_ados : bool
        Whether or not to store the HEOM ADOs.
      - normalize_output : bool
        Normalize output state to hide ODE numerical errors. Only normalize the state if the initial state is already normalized.
      - progress_bar : str {‘text’, ‘enhanced’, ‘tqdm’, ‘’}
        How to present the solver progress. ‘tqdm’ uses the python module of the same name and raise an error if not installed. Empty string or False will disable the bar.
      - progress_kwargs : dict
        kwargs to pass to the progress_bar. Qutip’s bars use chunk_size.
      - state_data_type: str {‘dense’, ‘CSR’, ‘Dia’, }
        Name of the data type of the state used during the ODE evolution. Use an empty string to keep the input state type. Many integrator can only work with Dense.
      - method : str [“adams”, “bdf”, “lsoda”, “dop853”, “vern9”, etc.]
        Which differential equation integration method to use.
      - atol, rtol : float
        Absolute and relative tolerance of the ODE integrator.
      - nsteps : int
        Maximum number of (internally defined) steps allowed in one `tlist` step.
      - max_step : float,
        Maximum lenght of one internal step. When using pulses, it should be less than half the width of the thinnest pulse.
  Returns:
    `HEOMResult`

      The results of the simulation run, with the following important attributes:

      - `times`: the times `t` (i.e. the `tlist`).

      - `states`: the system state at each time `t` (only available if `e_ops` was `None` or if the solver option `store_states` was set to `True`).

      - `ado_states`: the full ADO state at each time (only available if the results option `ado_return` was set to `True`). Each element is an instance of `HierarchyADOsState`. The state of a particular ADO may be extracted from `result.ado_states[i]` by calling `extract`.

      - `expect`: a list containing the values of each `e_ops` at time `t`.

      - `e_data`: a dictionary containing the values of each `e_ops` at tme `t`. The keys are those given by `e_ops` if it was a dict, otherwise they are the indexes of the supplied `e_ops`.

      See `HEOMResult` and `Result` for the complete list of attributes.

class HEOMSolver(H, bath, max_depth, *, odd_parity=False, options=None)

  HEOM solver that supports multiple baths.

  Each bath must be either bosonic or fermionic, but bosonic and fermionic baths can be mixed.

  Parameters:
    H `Qobj`, `QobjEvo`

      Possibly time-dependent system Liouvillian or Hamiltonian as a Qobj or QobjEvo. list of [`Qobj`, `Coefficient`] or callable that can be made into `QobjEvo` are also accepted.

    bath Bath specification or list of Bath specifications

      A bath containing the exponents of the expansion of the bath correlation funcion and their associated coefficients and coupling operators, or a list of baths.

      Each bath can be specified as either an object of type `Bath`, `BosonicBath`, `FermionicBath`, or their subtypes, or as a tuple `(env, Q)`, where `env` is an `ExponentialBosonicEnvironment` or an `ExponentialFermionicEnvironment` and `Q` the system coupling operator.

    odd_parity Bool

      For fermionic baths only. Default is “False”. “Parity” refers to the parity of the initial system state used with the HEOM. An example of an odd parity state is one made from applying an odd number of fermionic creation operators to a physical density operator. Physical systems have even parity, but allowing the generalization to odd-parity states allows one to calculate useful physical quantities like the system power spectrum or density of states. The form of the HEOM differs depending on the parity of the initial system state, so if this option is set to “True”, a different RHS is constructed, which can then be used with a system state of odd parity.

    max_depth int

      The maximum depth of the hierarchy (i.e. the maximum number of bath exponent “excitations” to retain).

    options dict, optional

      Generic solver options. If set to None the default options will be used. Keyword only. Default: None.

  Attributes:
    ados `HierarchyADOs`

      The description of the hierarchy constructed from the given bath and maximum depth.

    rhs `QobjEvo`

      The right-hand side (RHS) of the hierarchy evolution ODE. Internally the system and bath coupling operators are converted to `qutip.data.CSR` instances during construction of the RHS, so the operators in the `rhs` will all be sparse.

  property options

    Options for HEOMSolver:

    store_final_state: bool, default: False

      Whether or not to store the final state of the evolution in the result class.

    store_states: bool, default: None

      Whether or not to store the state vectors or density matrices. On None the states will be saved if no expectation operators are given.

    normalize_output: bool, default: False

      Normalize output state to hide ODE numerical errors.

    progress_bar: str {‘text’, ‘enhanced’, ‘tqdm’, ‘’}, default: “text”

      How to present the solver progress. ‘tqdm’ uses the python module of the same name and raise an error if not installed. Empty string or False will disable the bar.

    progress_kwargs: dict, default: {“chunk_size”: 10}

      Arguments to pass to the progress_bar. Qutip’s bars use `chunk_size`.

    method: str, default: “adams”

      Which ordinary differential equation integration method to use.

    state_data_type: str, default: “dense”

      Name of the data type of the state used during the ODE evolution. Use an empty string to keep the input state type. Many integrators support only work with Dense.

    store_adosbool, default: False

      Whether or not to store the HEOM ADOs. Only relevant when using the HEOM solver.

  run(state0, tlist, *, args=None, e_ops=None)

    Solve for the time evolution of the system.

    Parameters:
      state0 `Qobj` or `HierarchyADOsState` or array-like

        If `rho0` is a `Qobj` the it is the initial state of the system (i.e. a `Qobj` density matrix).

        If it is a `HierarchyADOsState` or array-like, then `rho0` gives the initial state of all ADOs.

        Usually the state of the ADOs would be determine from a previous call to `.run(...)` with the solver results option `store_ados` set to True. For example, `result = solver.run(...)` could be followed by `solver.run(result.ado_states[-1], tlist)`.

        If a numpy array-like is passed its shape must be `(number_of_ados, n, n)` where `(n, n)` is the system shape (i.e. shape of the system density matrix) and the ADOs must be in the same order as in `.ados.labels`.

      tlist list

        An ordered list of times at which to return the value of the state.

      args dict, optional {None}

        Change the `args` of the RHS for the evolution.

      e_ops Qobj / QobjEvo / callable / list / dict / None, optional

        A list or dictionary of operators as `Qobj`, `QobjEvo` and/or callable functions (they can be mixed) or a single operator or callable function. For an operator `op`, the result will be computed using `(state * op).tr()` and the state at each time `t`. For callable functions, `f`, the result is computed using `f(t, ado_state)`. The values are stored in the `expect` and `e_data` attributes of the result (see the return section below).

    Returns:
      `HEOMResult`

        The results of the simulation run, with the following important attributes:

        - `times`: the times `t` (i.e. the `tlist`).

        - `states`: the system state at each time `t` (only available if `e_ops` was `None` or if the solver option `store_states` was set to `True`).

        - `ado_states`: the full ADO state at each time (only available if the results option `ado_return` was set to `True`). Each element is an instance of `HierarchyADOsState`. The state of a particular ADO may be extracted from `result.ado_states[i]` by calling `extract`.

        - `expect`: a list containing the values of each `e_ops` at time `t`.

        - `e_data`: a dictionary containing the values of each `e_ops` at tme `t`. The keys are those given by `e_ops` if it was a dict, otherwise they are the indexes of the supplied `e_ops`.

        See `HEOMResult` and `Result` for the complete list of attributes.

  start(state0, t0)

    Set the initial state and time for a step evolution.

    Parameters:
      state0 `Qobj`

        Initial state of the evolution. This may provide either just the initial density matrix of the system, or the full set of ADOs for the hierarchy. See the documentation for `rho0` in the `.run(...)` method for details.

      t0 double

        Initial time of the evolution.

  steady_state(
    use_mkl=True,
    mkl_max_iter_refine=100,
    mkl_weighted_matching=False,
  )

    Compute the steady state of the system.

    Parameters:
      use_mkl bool, default=False

        Whether to use mkl or not. If mkl is not installed or if this is false, use the scipy splu solver instead.

      mkl_max_iter_refine int

        Specifies the the maximum number of iterative refinement steps that the MKL PARDISO solver performs.

        For a complete description, see iparm(7) in https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-0/pardiso-iparm-parameter.html

      mkl_weighted_matching bool

        MKL PARDISO can use a maximum weighted matching algorithm to permute large elements close the diagonal. This strategy adds an additional level of reliability to the factorization methods.

        For a complete description, see iparm(12) in https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2023-0/pardiso-iparm-parameter.html

    Returns:
      steady_state Qobj

        The steady state density matrix of the system.

      steady_ados `HierarchyADOsState`

        The steady state of the full ADO hierarchy. A particular ADO may be extracted from the full state by calling `extract`.

class HSolverDL(

  H_sys,
  coup_op,
  coup_strength,
  temperature,
  N_cut,
  N_exp,
  cut_freq,
  *,
  bnd_cut_approx=False,
  options=None,
  combine=True,

)

  A helper class for creating an `HEOMSolver` that is backwards compatible with the `HSolverDL` provided in `qutip.nonmarkov.heom` in QuTiP 4.6 and below.

  See `HEOMSolver` and `DrudeLorentzBath` for more descriptions of the underlying solver and bath construction.

  Note

  Unlike the version of `HSolverDL` in QuTiP 4.6, this solver supports supplying a time-dependent or Liouvillian `H_sys`.

  Note

  For compatibility with `HSolverDL` in QuTiP 4.6 and below, the parameter `N_exp` specifying the number of exponents to keep in the expansion of the bath correlation function is one more than the equivalent `Nk` used in the `DrudeLorentzBath`. I.e., `Nk = N_exp - 1`. The `Nk` parameter in the `DrudeLorentzBath` does not count the zeroeth exponent in order to better match common usage in the literature.

  Note

  The `stats` and `renorm` arguments accepted in QuTiP 4.6 and below are no longer supported.

  Parameters:
    H_sys Qobj or QobjEvo or list

      The system Hamiltonian or Liouvillian. See `HEOMSolver` for a complete description.

    coup_op Qobj

      Operator describing the coupling between system and bath. See parameter `Q` in `BosonicBath` for a complete description.

    coup_strength float

      Coupling strength. Referred to as `lam` in `DrudeLorentzEnvironment`.

    temperature float

      Bath temperature. Referred to as `T` in `DrudeLorentzEnvironment`.

    N_cut int

      The maximum depth of the hierarchy. See `max_depth` in `HEOMSolver` for a full description.

    N_exp int

      Number of exponential terms used to approximate the bath correlation functions. The equivalent `Nk` in `DrudeLorentzBath` is one less than `N_exp` (see note above).

    cut_freq float

      Bath spectral density cutoff frequency. Referred to as `gamma` in `DrudeLorentzEnvironment`.

    bnd_cut_approx bool

      Use boundary cut off approximation. If true, the Matsubara terminator is added to the system Liouvillian (and H_sys is promoted to a Liouvillian if it was a Hamiltonian). Keyword only. Default: False.

    options dict, optional

      Generic solver options. If set to None the default options will be used. Keyword only. Default: None.

    combine bool, default: True

      Whether to combine exponents with the same frequency (and coupling operator). See `ExponentialBosonicEnvironment.combine` for details. Keyword only. Default: True.

class HierarchyADOs(exponents, max_depth)

  A description of ADOs (auxilliary density operators) with the hierarchical equations of motion.

  The list of ADOs is constructed from a list of bath exponents (corresponding to one or more baths). Each ADO is referred to by a label that lists the number of “excitations” of each bath exponent. The level of a label within the hierarchy is the sum of the “excitations” within the label.

  For example the label `(0, 0, ..., 0)` represents the density matrix of the system being solved and is the only 0th level label.

  The labels with a single 1, i.e. `(1, 0, ..., 0)`, `(0, 1, 0, ... 0)`, etc. are the 1st level labels.

  The second level labels all have either two 1s or a single 2, and so on for the third and higher levels of the hierarchy.

  Parameters:
    exponents list of `BathExponent`

      The exponents of the correlation function describing the bath or baths.

    max_depth int

      The maximum depth of the hierarchy (i.e. the maximum sum of “excitations” in the hierarchy ADO labels or maximum ADO level).

  Attributes:
    exponents list of `BathExponent`

      The exponents of the correlation function describing the bath or baths.

    max_depth int

      The maximum depth of the hierarchy (i.e. the maximum sum of “excitations” in the hierarchy ADO labels).

    dims list of int

      The dimensions of each exponent within the bath(s).

    vk list of complex

      The frequency of each exponent within the bath(s).

    ck list of complex

      The coefficient of each exponent within the bath(s).

    ck2: list of complex

      For exponents of type “RI”, the coefficient of the exponent within the imaginary expansion. For other exponent types, the entry is None.

    sigma_bar_k_offset: list of int

      For exponents of type “+” or “-” the offset within the list of modes of the corresponding “-” or “+” exponent. For other exponent types, the entry is None.

    labels: list of tuples

      A list of the ADO labels within the hierarchy.

  exps(label)

    Converts an ADO label into a tuple of exponents, with one exponent for each “excitation” within the label.

    The number of exponents returned is always equal to the level of the label within the hierarchy (i.e. the sum of the indices within the label).

    Parameters:
      label tuple

        The ADO label to convert to a list of exponents.

    Returns:
      tuple of `BathExponent`

        A tuple of BathExponents.

    Examples

    `ados.exps((1, 0, 0))` would return `[ados.exponents[0]]`

    `ados.exps((2, 0, 0))` would return `[ados.exponents[0], ados.exponents[0]]`.

    `ados.exps((1, 2, 1))` would return `[ados.exponents[0], ados.exponents[1], ados.exponents[1], ados.exponents[2]]`.

  filter(level=None, tags=None, dims=None, types=None)

    Return a list of ADO labels for ADOs whose “excitations” match the given patterns.

    Each of the filter parameters (tags, dims, types) may be either unspecified (None) or a list. Unspecified parameters are excluded from the filtering.

    All specified filter parameters must be lists of the same length. Each position in the lists describes a particular excitation and any exponent that matches the filters may supply that excitation. The level of all labels returned is thus equal to the length of the filter parameter lists.

    Within a filter parameter list, items that are None represent wildcards and match any value of that exponent attribute

    Parameters:
      level int

        The hierarchy depth to return ADOs from.

      tags list of object or None

        Filter parameter that matches the `.tag` attribute of exponents.

      dims list of int

        Filter parameter that matches the `.dim` attribute of exponents.

      types list of BathExponent types or list of str

        Filter parameter that matches the `.type` attribute of exponents. Types may be supplied by name (e.g. “R”, “I”, “+”) instead of by the actual type (e.g. `CFExponent.types.R`).

    Returns:
      list of tuple

        The ADO label for each ADO whose exponent excitations (i.e. label) match the given filters or level.

  idx(label)

    Return the index of the ADO label within the list of labels, i.e. within `self.labels`.

    Parameters:
      label tuple

        The label to look up.

    Returns:
      int

        The index of the label within the list of ADO labels.

    Notes

    This implementation of the `.idx(...)` method is just for reference and documentation. To avoid the cost of a Python function call, it is replaced with `self._label_idx.__getitem__` when the instance is created.

  next(label, k)

    Return the ADO label with one more excitation in the k’th exponent dimension or `None` if adding the excitation would exceed the dimension or maximum depth of the hierarchy.

    Parameters:
      label tuple

        The ADO label to add an excitation to.

      k int

        The exponent to add the excitation to.

    Returns:
      tuple or None

        The next label.

  prev(label, k)

    Return the ADO label with one fewer excitation in the k’th exponent dimension or `None` if the label has no exciations in the k’th exponent.

    Parameters:
      label tuple

        The ADO label to remove the excitation from.

      k int

        The exponent to remove the excitation from.

    Returns:
      tuple or None

        The previous label.

class HierarchyADOsState(rho, ados, ado_state)

  Provides convenient access to the full hierarchy ADO state at a particular point in time, `t`.

  Parameters:
    rho `Qobj`

      The current state of the system (i.e. the 0th component of the hierarchy).

    ados `HierarchyADOs`

      The description of the hierarchy.

    ado_state numpy.array

      The full state of the hierarchy.

  Attributes:
    rho Qobj

      The system state.

  Notes

  In addition, all of the attributes of the hierarchy description, i.e. `HierarchyADOs`, are provided directly on this class for convenience. E.g. one can access `.labels`, or `.exponents` or call `.idx(label)` directly.

  See `HierarchyADOs` for a full list of the available attributes and methods.

  extract(idx_or_label)

    Extract a Qobj representing the specified ADO from a full representation of the ADO states.

    Parameters:
      idx int or label

        The index of the ADO to extract. If an ADO label, e.g. `(0, 1, 0, ...)` is supplied instead, then the ADO is extracted by label instead.

    Returns:
      Qobj

        A `Qobj` representing the state of the specified ADO.

class HEOMResult(

  e_ops: dict[Any, Qobj | QobjEvo | Callable[[float, Qobj], Any]],
  options: ResultOptions,
  *,
  solver: str = None,
  stats: dict[str, Any] = None,
  **kw,

)

## Baths

class BathExponent(type, dim, Q, ck, vk, ck2=None, sigma_bar_k_offset=None, tag=None)

  Represents a single exponent (naively, an excitation mode) within the decomposition of the correlation functions of a bath.

  This class extends the `CFExponent` of the environment-module for use with the HEOM solver. In addition to the exponent itself, the BathExponent keeps track of the corresponding system coupling operator `Q`, as well as the parameters `dim` and `sigma_bar_k_offset`.

  Parameters:
    type{“R”, “I”, “RI”, “+”, “-”} or `CFExponent.ExponentType`

      The type of bath exponent.

      “R” and “I” are bosonic bath exponents that appear in the real and imaginary parts of the correlation expansion.

      “RI” is combined bosonic bath exponent that appears in both the real and imaginary parts of the correlation expansion. The combined exponent has a single `vk`. The `ck` is the coefficient in the real expansion and `ck2` is the coefficient in the imaginary expansion.

      “+” and “-” are fermionic bath exponents. These fermionic bath exponents must specify `sigma_bar_k_offset` which specifies the amount to add to `k` (the exponent index within the bath of this exponent) to determine the `k` of the corresponding exponent with the opposite sign (i.e. “-” or “+”).

    dim int or None

      The dimension (i.e. maximum number of excitations for this exponent). Usually `2` for fermionic exponents or `None` (i.e. unlimited) for bosonic exponents.

    Q Qobj

      The coupling operator for this excitation mode.

    vk complex

      The frequency of the exponent of the excitation term.

    ck complex

      The coefficient of the excitation term.

    ck2 optional, complex

      For exponents of type “RI” this is the coefficient of the term in the imaginary expansion (and `ck` is the coefficient in the real expansion).

    sigma_bar_k_offset optional, int

      For exponents of type “+” this gives the offset (within the list of exponents within the bath) of the corresponding “-” bath exponent. For exponents of type “-” it gives the offset of the corresponding “+” exponent.

    tag optional, str, tuple or any other object

      A label for the exponent (often the name of the bath). It defaults to None.

  Attributes:
    fermionic bool

      True if the type of the exponent is a Fermionic type (i.e. either “+” or “-”) and False otherwise.

    All of the parameters are also available as attributes.
  rescale(alpha)

    Rescale the coefficient of the exponent by a factor of alpha.

  types

    alias of `ExponentType`

class Bath(exponents)

  Represents a list of bath expansion exponents.

  Parameters:
    exponents list of `BathExponent`

      The exponents of the correlation function describing the bath.

class BosonicBath(Q, ck_real, vk_real, ck_imag, vk_imag, combine=True, tag=None)

  A helper class for constructing a bosonic bath from the expansion coefficients and frequencies for the real and imaginary parts of the bath correlation function.

  If the correlation functions `C(t)` is split into real and imaginary parts:

      C(t) = C_real(t) + i * C_imag(t)

  then:

      C_real(t) = sum(ck_real * exp(- vk_real * t))
      C_imag(t) = sum(ck_imag * exp(- vk_imag * t))

  Defines the coefficients `ck` and the frequencies `vk`.

  Note that the `ck` and `vk` may be complex, even through `C_real(t)` and `C_imag(t)` (i.e. the sum) is real.

  Parameters:
    Q Qobj

      The coupling operator for the bath.

    ck_real list of complex

      The coefficients of the expansion terms for the real part of the correlation function. The corresponding frequencies are passed as vk_real.

    vk_real list of complex

      The frequencies (exponents) of the expansion terms for the real part of the correlation function. The corresponding ceofficients are passed as ck_real.

    ck_imag list of complex

      The coefficients of the expansion terms in the imaginary part of the correlation function. The corresponding frequencies are passed as vk_imag.

    vk_imag list of complex

      The frequencies (exponents) of the expansion terms for the imaginary part of the correlation function. The corresponding ceofficients are passed as ck_imag.

    combine bool, default True

      Whether to combine exponents with the same frequency (and coupling operator). See `ExponentialBosonicEnvironment.combine` for details.

    tag optional, str, tuple or any other object

      A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from.

  Notes

  This class is part of the “bath” API, which is now mirrored by the newer “environment” API. The bath classes are kept in QuTiP for reasons of backwards compatibility and convenience. This class is an extended version of the `ExponentialBosonicEnvironment`, adding the parameter `Q`, which is not included in the newer “environment” API.

  classmethod from_environment(env, Q, dim=None)

    Converts from the “environment” API to the “bath” API. A BosonicBath combines the information from an ExponentialBosonicEnvironment and a coupling operator.

    Parameters:
      env `ExponentialBosonicEnvironment`

        The bath.

      Q Qobj

        The coupling operator for the bath.

      dim optional, int or `None` (default `None`)

        The maximum number of excitations for each exponent. Usually `None` (i.e. unlimited).

class DrudeLorentzBath(Q, lam, gamma, T, Nk, combine=True, tag=None)

  A helper class for constructing a Drude-Lorentz bosonic bath from the bath parameters (see parameters below).

  Parameters:
    Q Qobj

      Operator describing the coupling between system and bath.

    lam float

      Coupling strength.

    gamma float

      Bath spectral density cutoff frequency.

    T float

      Bath temperature.

    Nk int

      Number of exponential terms used to approximate the bath correlation functions.

    combine bool, default True

      Whether to combine exponents with the same frequency (and coupling operator). See `ExponentialBosonicEnvironment.combine` for details.

    tag optional, str, tuple or any other object

      A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from.

  Notes

  This class is part of the “bath” API, which is now mirrored by the newer “environment” API. The bath classes are kept in QuTiP for reasons of backwards compatibility and convenience. Creating a DrudeLorentzBath is equivalent to creating a `DrudeLorentzEnvironment`, performing a `Matsubara` approximation, and finally bundling the result together with the coupling operator `Q` for convenient use with the HEOM solver.

  terminator()

    Return the Matsubara terminator for the bath and the calculated approximation discrepancy.

    Returns:
      delta: float

        The approximation discrepancy. That is, the difference between the true correlation function of the Drude-Lorentz bath and the sum of the `Nk` exponential terms is approximately `2 * delta * dirac(t)`, where `dirac(t)` denotes the Dirac delta function.

      terminator Qobj

        The Padé terminator – i.e. a liouvillian term representing the contribution to the system-bath dynamics of all exponential expansion terms beyond `Nk`. It should be used by adding it to the system liouvillian (i.e. `liouvillian(H_sys)`).

class DrudeLorentzPadeBath(Q, lam, gamma, T, Nk, combine=True, tag=None)

  A helper class for constructing a Padé expansion for a Drude-Lorentz bosonic bath from the bath parameters (see parameters below).

  A Padé approximant is a sum-over-poles expansion (see https://en.wikipedia.org/wiki/Pad%C3%A9_approximant).

  The application of the Padé method to spectrum decompoisitions is described in “Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems” [1].

  The implementation here follows the approach in the paper.

  [1] J. Chem. Phys. 134, 244106 (2011); https://doi.org/10.1063/1.3602466

  This is an alternative to the `DrudeLorentzBath` which constructs a simpler exponential expansion.

  Parameters:
    Q Qobj

      Operator describing the coupling between system and bath.

    lam float

      Coupling strength.

    gamma float

      Bath spectral density cutoff frequency.

    T float

      Bath temperature.

    Nk int

      Number of Padé exponentials terms used to approximate the bath correlation functions.

    combine bool, default True

      Whether to combine exponents with the same frequency (and coupling operator). See `ExponentialBosonicEnvironment.combine` for details.

    tag optional, str, tuple or any other object

      A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from.

  Notes

  This class is part of the “bath” API, which is now mirrored by the newer “environment” API. The bath classes are kept in QuTiP for reasons of backwards compatibility and convenience. Creating a DrudeLorentzPadeBath is equivalent to creating a `DrudeLorentzEnvironment`, performing a `Pade` approximation, and finally bundling the result together with the coupling operator `Q` for convenient use with the HEOM solver.

  terminator()

    Return the Padé terminator for the bath and the calculated approximation discrepancy.

    Returns:
      delta: float

        The approximation discrepancy. That is, the difference between the true correlation function of the Drude-Lorentz bath and the sum of the `Nk` exponential terms is approximately `2 * delta * dirac(t)`, where `dirac(t)` denotes the Dirac delta function.

      terminator Qobj

        The Padé terminator – i.e. a liouvillian term representing the contribution to the system-bath dynamics of all exponential expansion terms beyond `Nk`. It should be used by adding it to the system liouvillian (i.e. `liouvillian(H_sys)`).

class UnderDampedBath(Q, lam, gamma, w0, T, Nk, combine=True, tag=None)

  A helper class for constructing an under-damped bosonic bath from the bath parameters (see parameters below).

  Parameters:
    Q Qobj

      Operator describing the coupling between system and bath.

    lam float

      Coupling strength.

    gamma float

      Bath spectral density cutoff frequency.

    w0 float

      Bath spectral density resonance frequency.

    T float

      Bath temperature.

    Nk int

      Number of exponential terms used to approximate the bath correlation functions.

    combine bool, default True

      Whether to combine exponents with the same frequency (and coupling operator). See `ExponentialBosonicEnvironment.combine` for details.

    tag optional, str, tuple or any other object

      A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from.

  Notes

  This class is part of the “bath” API, which is now mirrored by the newer “environment” API. The bath classes are kept in QuTiP for reasons of backwards compatibility and convenience. Creating an UnderDampedBath is equivalent to creating an `UnderDampedEnvironment`, performing a `Matsubara` approximation, and finally bundling the result together with the coupling operator `Q` for convenient use with the HEOM solver.

class FermionicBath(Q, ck_plus, vk_plus, ck_minus, vk_minus, tag=None)

  A helper class for constructing a fermionic bath from the expansion coefficients and frequencies for the `+` and `-` modes of the bath correlation function.

  There must be the same number of `+` and `-` modes and their coefficients must be specified in the same order so that `ck_plus[i], vk_plus[i]` are the plus coefficient and frequency corresponding to the minus mode `ck_minus[i], vk_minus[i]`.

  In the fermionic case the order in which excitations are created or destroyed is important, resulting in two different correlation functions labelled `C_plus(t)` and `C_plus(t)`:

      C_plus(t) = sum(ck_plus * exp(- vk_plus * t))
      C_minus(t) = sum(ck_minus * exp(- vk_minus * t))

  where the expansions above define the coeffiients `ck` and the frequencies `vk`.

  Parameters:
    Q Qobj

      The coupling operator for the bath.

    ck_plus list of complex

      The coefficients of the expansion terms for the `+` part of the correlation function. The corresponding frequencies are passed as vk_plus.

    vk_plus list of complex

      The frequencies (exponents) of the expansion terms for the `+` part of the correlation function. The corresponding ceofficients are passed as ck_plus.

    ck_minus list of complex

      The coefficients of the expansion terms for the `-` part of the correlation function. The corresponding frequencies are passed as vk_minus.

    vk_minus list of complex

      The frequencies (exponents) of the expansion terms for the `-` part of the correlation function. The corresponding ceofficients are passed as ck_minus.

    tag optional, str, tuple or any other object

      A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from.

  Notes

  This class is part of the “bath” API, which is now mirrored by the newer “environment” API. The bath classes are kept in QuTiP for reasons of backwards compatibility and convenience. This class is an extended version of the `ExponentialFermionicEnvironment`, adding the parameter `Q`, which is not included in the newer “environment” API.

  classmethod from_environment(env, Q, dim=2)

    Converts from the “environment” API to the “bath” API. A FermionicBath combines the information from an ExponentialFermionicEnvironment and a coupling operator.

    Note that the “environment” API does not require fermionic exponents to come in pairs. This method will add additional exponents with zero coefficient in order to achieve the same amount of `+` and `-` exponents.

    Parameters:
      env `ExponentialFermionicEnvironment`

        The bath.

      Q Qobj

        The coupling operator for the bath.

      dim optional, int or `None` (default `2`)

        The maximum number of excitations for each exponent. Usually `2`.

class LorentzianBath(Q, gamma, w, mu, T, Nk, tag=None)

  A helper class for constructing a Lorentzian fermionic bath from the bath parameters (see parameters below).

  Note

  The Matsubara expansion used in this bath converges very slowly and `Nk > 20` may be required to get good convergence. The Padé expansion used by `LorentzianPadeBath` converges much more quickly.

  Parameters:
    Q Qobj

      Operator describing the coupling between system and bath.

    gamma float

      The coupling strength between the system and the bath.

    w float

      The width of the environment.

    mu float

      The chemical potential of the bath.

    T float

      Bath temperature.

    Nk int

      Number of exponential terms used to approximate the bath correlation functions.

    tag optional, str, tuple or any other object

      A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from.

  Notes

  This class is part of the “bath” API, which is now mirrored by the newer “environment” API. The bath classes are kept in QuTiP for reasons of backwards compatibility and convenience. Creating a LorentzianBath is equivalent to creating a `LorentzianEnvironment`, performing a `Matsubara` approximation, and finally bundling the result together with the coupling operator `Q` for convenient use with the HEOM solver.

class LorentzianPadeBath(Q, gamma, w, mu, T, Nk, tag=None)

  A helper class for constructing a Padé expansion for Lorentzian fermionic bath from the bath parameters (see parameters below).

  A Padé approximant is a sum-over-poles expansion ( see https://en.wikipedia.org/wiki/Pad%C3%A9_approximant).

  The application of the Padé method to spectrum decompoisitions is described in “Padé spectrum decompositions of quantum distribution functions and optimal hierarchical equations of motion construction for quantum open systems” [1].

  The implementation here follows the approach in the paper.

  [1] J. Chem. Phys. 134, 244106 (2011); https://doi.org/10.1063/1.3602466

  This is an alternative to the `LorentzianBath` which constructs a simpler exponential expansion that converges much more slowly in this particular case.

  Parameters:
    Q Qobj

      Operator describing the coupling between system and bath.

    gamma float

      The coupling strength between the system and the bath.

    w float

      The width of the environment.

    mu float

      The chemical potential of the bath.

    T float

      Bath temperature.

    Nk int

      Number of exponential terms used to approximate the bath correlation functions.

    tag optional, str, tuple or any other object

      A label for the bath exponents (for example, the name of the bath). It defaults to None but can be set to help identify which bath an exponent is from.

  Notes

  This class is part of the “bath” API, which is now mirrored by the newer “environment” API. The bath classes are kept in QuTiP for reasons of backwards compatibility and convenience. Creating a LorentzianPadeBath is equivalent to creating a `LorentzianEnvironment`, performing a `Pade` approximation, and finally bundling the result together with the coupling operator `Q` for convenient use with the HEOM solver.