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:

HQobj, 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.

bathBath 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_depthint

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

state0Qobj 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.

tlistlist

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

e_opsQobj / 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).

argsdict, optional

Change the args of the RHS for the evolution.

optionsdict, 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:

HQobj, 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.

bathBath 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_parityBool

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_depthint

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

optionsdict, optional

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

Attributes:

adosHierarchyADOs

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

rhsQobjEvo

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:

state0Qobj 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.

tlistlist

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

argsdict, optional {None}

Change the args of the RHS for the evolution.

e_opsQobj / 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:

state0Qobj

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.

t0double

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_mklbool, 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_refineint

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_matchingbool

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_stateQobj

The steady state density matrix of the system.

steady_adosHierarchyADOsState

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

property sys_dims

Dimensions of the space that the system use, excluding any environment:

qutip.basis(sovler.dims) will create a state with proper dimensions for this solver.

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_sysQobj or QobjEvo or list

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

coup_opQobj

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

coup_strengthfloat

Coupling strength. Referred to as lam in DrudeLorentzEnvironment.

temperaturefloat

Bath temperature. Referred to as T in DrudeLorentzEnvironment.

N_cutint

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

N_expint

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_freqfloat

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

bnd_cut_approxbool

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.

optionsdict, optional

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

combinebool, 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:

exponentslist of BathExponent

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

max_depthint

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

Attributes:

exponentslist of BathExponent

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

max_depthint

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

dimslist of int

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

vklist of complex

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

cklist 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:

labeltuple

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:

levelint

The hierarchy depth to return ADOs from.

tagslist of object or None

Filter parameter that matches the .tag attribute of exponents.

dimslist of int

Filter parameter that matches the .dim attribute of exponents.

typeslist 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:

labeltuple

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:

labeltuple

The ADO label to add an excitation to.

kint

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:

labeltuple

The ADO label to remove the excitation from.

kint

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:

rhoQobj

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

adosHierarchyADOs

The description of the hierarchy.

ado_statenumpy.array

The full state of the hierarchy.

Attributes:

rhoQobj

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:

idxint 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 “+”).

dimint 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.

QQobj

The coupling operator for this excitation mode.

vkcomplex

The frequency of the exponent of the excitation term.

ckcomplex

The coefficient of the excitation term.

ck2optional, 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_offsetoptional, 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.

tagoptional, str, tuple or any other object

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

Attributes:

fermionicbool

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.

types

alias of ExponentType

class Bath(exponents)

Represents a list of bath expansion exponents.

Parameters:

exponentslist 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:

QQobj

The coupling operator for the bath.

ck_reallist 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_reallist 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_imaglist 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_imaglist 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.

combinebool, default True

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

tagoptional, 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:

envExponentialBosonicEnvironment

The bath.

QQobj

The coupling operator for the bath.

dimoptional, 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:

QQobj

Operator describing the coupling between system and bath.

lamfloat

Coupling strength.

gammafloat

Bath spectral density cutoff frequency.

Tfloat

Bath temperature.

Nkint

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

combinebool, default True

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

tagoptional, 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.

terminatorQobj

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:

QQobj

Operator describing the coupling between system and bath.

lamfloat

Coupling strength.

gammafloat

Bath spectral density cutoff frequency.

Tfloat

Bath temperature.

Nkint

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

combinebool, default True

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

tagoptional, 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.

terminatorQobj

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:

QQobj

Operator describing the coupling between system and bath.

lamfloat

Coupling strength.

gammafloat

Bath spectral density cutoff frequency.

w0float

Bath spectral density resonance frequency.

Tfloat

Bath temperature.

Nkint

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

combinebool, default True

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

tagoptional, 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:

QQobj

The coupling operator for the bath.

ck_pluslist of complex

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

vk_pluslist 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_minuslist of complex

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

vk_minuslist of complex

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

tagoptional, 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:

envExponentialFermionicEnvironment

The bath.

QQobj

The coupling operator for the bath.

dimoptional, 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

This 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:

QQobj

Operator describing the coupling between system and bath.

gammafloat

The coupling strength between the system and the bath.

wfloat

The width of the environment.

mufloat

The chemical potential of the bath.

Tfloat

Bath temperature.

Nkint

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

tagoptional, 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:

QQobj

Operator describing the coupling between system and bath.

gammafloat

The coupling strength between the system and the bath.

wfloat

The width of the environment.

mufloat

The chemical potential of the bath.

Tfloat

Bath temperature.

Nkint

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

tagoptional, 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.