quspin.tools.measurements.diag_ensemble

quspin.tools.measurements.diag_ensemble(N, system_state, E2, V2, density=True, alpha=1.0, rho_d=False, Obs=False, delta_t_Obs=False, delta_q_Obs=False, Sd_Renyi=False, Srdm_Renyi=False, Srdm_args={})[source]

Calculates expectation values in the Diagonal ensemble of the initial state.

Equivalently, these are also the infinite-time expectation values after a sudden quench from a Hamiltonian \(H_1\) to a Hamiltonian \(H_2\). Let us label the two eigenbases by

\[V_1=\{|n_1\rangle: H_1|n_1\rangle=E_1|n_1\rangle\} \qquad V_2=\{|n_2\rangle: H_2|n_2\rangle=E_2|n_2\rangle\}\]

See eg. arXiv:1509.06411 for the physical definition of Diagonal Ensemble.

Note: All expectation values depend statistically on the symmetry block used via the available number of states, due to the generic system-size dependence!

Parameters:

Nint

System size/dimension (e.g. number of sites).

system_state{array_like,dict}

State of the quantum system; can be either of:

• numpy.ndarray: pure state, shape = (Ns,) or (,Ns).

• numpy.ndarray: density matrix (DM), shape = (Ns,Ns).

• dict: mixed DM as dictionary {“V1”:V1, “E1”:E1, “f”:f, “f_args”:f_args, “V1_state”:int, “f_norm”:`False}` to define a diagonal DM in the basis \(V_1\) of the Hamiltonian \(H_1\). The meaning of the keys (keys CANNOT be chosen arbitrarily) is as flollows:

  • numpy.ndarray: V1 (required) contains eigenbasis of \(H_1\) in the columns.

  • numpy.ndarray: E1 (required) eigenenergies of \(H_1\).

  • function ‘f’ (optional) is a function which represents the distribution of the spectrum

    used to define the mixed DM of the initial state (see example).

    Default is a thermal distribution with inverse temperature beta: f = lambda E,beta: numpy.exp(-beta*(E - E[0]) ).

  • list(float): f_args (required) list of arguments for function f.

    If f is not defined, by default we have \(f(E)=\exp(-\beta(E - E_\mathrm{GS}))\), and f_args=[beta] specifies the inverse temeprature.

  • list(int): V1_state (optional) is a list of integers to specify arbitrary states of V1

    whose pure expectations are also returned.

  • bool: f_norm (optional). If set to False the mixed DM built from f is NOT normalised

    and the norm is returned under the key f_norm.

    Use this option if you need to average your results over multiple symmetry blocks, which require a separate normalisations.

  If this option is specified, then all Diagonal Ensemble quantities are averaged over the energy distribution \(f(E_1,f\_args)\):

  \[\overline{\mathcal{M}_d} = \frac{1}{Z_f}\sum_{n_1} f(E_{n_1},f\_args)\mathcal{M}^{n_1}_d, \qquad \mathcal{M}^{\psi}_d = \langle\mathcal{O}\rangle_d^\psi,\ \delta_q\mathcal{O}^\psi_d,\ \delta_t\mathcal{O}^\psi_d,\ S_d^\psi,\ S_\mathrm{rdm}^\psi\]

V2numpy.ndarray

Contains the basis of the Hamiltonian \(H_2\) in the columns.

E2numpy.ndarray

Contains the eigenenergies corresponding to the eigenstates in V2.

This variable is only used to check for degeneracies, in which case the function is NOT expected to produce correct resultsand raises an error.

rho_dbool, optional

When set to True, returns the Diagonal ensemble DM. Default is False.

Adds the key “rho_d” to output.

For example, if system_state is the pure state \(|\psi\rangle\):

\[\rho_d^\psi = \sum_{n_2} \left|\langle\psi|n_2\rangle\right|^2\left|n_2\rangle\langle n_2\right| = \sum_{n_2} \left(\rho_d^\psi\right)_{n_2n_2}\left|n_2\rangle\langle n_2\right|\]

Obs:obj:, optional

Hermitian matrix of the same shape as V2, to calculate the Diagonal ensemble expectation value of.

Adds the key “Obs” to output. Can be of type numpy.ndarray or an instance of the hamiltonian class.

For example, if system_state is the pure state \(|\psi\rangle\):

\[\langle\mathcal{O}\rangle_d^\psi = \lim_{T\to\infty}\frac{1}{T}\int_0^T\mathrm{d}t \frac{1}{N}\langle\psi\left|\mathcal{O}(t)\right|\psi\rangle = \frac{1}{N}\sum_{n_2}\left(\rho_d^\psi\right)_{n_2n_2} \langle n_2\left|\mathcal{O}\right|n_2\rangle\]

delta_q_Obsbool, optional

QUANTUM fluctuations of the expectation of Obs at infinite times. Requires Obs. Calculates temporal fluctuations delta_t_Obs for along the way (see above).

Adds keys “delta_q_Obs” and “delta_t_Obs” to output.

For example, if system_state is the pure state \(|\psi\rangle\):

\[\delta_q\mathcal{O}^\psi_d = \frac{1}{N}\sqrt{\lim_{T\to\infty}\frac{1}{T}\int_0^T\mathrm{d}t \langle\psi\left| \mathcal{O}(t)\right| \psi\rangle^2 - \langle\mathcal{O}\rangle_d^2}= \frac{1}{N}\sqrt{ \sum_{n_2\neq m_2} \left(\rho_d^\psi\right)_{n_2n_2} [\mathcal{O}]^2_{n_2m_2} \left(\rho_d^\psi\right)_{m_2m_2} }\]

delta_t_Obsbool, optional

TEMPORAL fluctuations around infinite-time expectation of Obs. Requires Obs.

Adds the key “delta_t_Obs” to output.

For example, if system_state is the pure state \(|\psi\rangle\):

\[\delta_t\mathcal{O}^\psi_d = \frac{1}{N}\sqrt{ \lim_{T\to\infty}\frac{1}{T}\int_0^T\mathrm{d}t \langle\psi\left|[\mathcal{O}(t)]^2\right|\psi\rangle - \langle\psi\left|\mathcal{O}(t)\right|\psi\rangle^2} = \frac{1}{N}\sqrt{\langle\mathcal{O}^2\rangle_d - \langle\mathcal{O}\rangle_d^2 - \left(\delta_q\mathcal{O}^\psi_d\right)^2 }\]

alphafloat, optional

Renyi \(alpha\) parameter. Default is alpha = 1.0.

Sd_Renyibool, optional

Computes the DIAGONAL Renyi entropy in the basis of \(H_2\). The default Renyi parameter is alpha=1.0 (see below). Adds the key “Sd_Renyi” to output. For example, if system_state is the pure state \(|\psi\rangle\):

\[S_d^\psi = \frac{1}{1-\alpha}\log\mathrm{tr}\left(\rho_d^\psi\right)^\alpha\]

Srdm_Renyibool, optional

Computes ENTANGLEMENT Renyi entropy of a subsystem (see Srdm_args for subsystem definition).

Requires passing the (otherwise optional) argument Srdm_args (see below).

The default Renyi parameter is alpha=1.0 (see below).

Adds the key “Srdm_Renyi” to output.

For example, if system_state is the pure state \(|\psi\rangle\) (see also notation in documentation of ent_entropy):

\[S_\mathrm{rdm}^\psi = \frac{1}{1-\alpha}\log \mathrm{tr}_{A} \left( \mathrm{tr}_{A^c} \rho_d^\psi \right)^\alpha\]

Srdm_argsdict, semi-optional

Dictionary which contains all arguments required for the computation of the entanglement Renyi entropy. Required when Srdm_Renyi = True. The following keys are allowed/supported:

• “basis”: obj(basis), required

  Basis used to build system_state in. Must be an instance of the basis class.

• “chain_subsys”list, optional

  Lattice sites to specify the chain subsystem of interest. Default is:

  – [0,1,…,N/2-1,N/2] for spin_basis_1d, fermion_basis_1d, boson_basis_1d.

  – [0,1,…,N-1,N] for photon_basis.

densitybool, optional

If set to True, all observables are normalised by the system size N, except for the Srdm_Renyi which is normalised by the subsystem size, i.e. by the length of chain_subsys. Default is ‘True’.

Returns:

dict

The following keys of the output are possible, depending on the choice of flags:

• “rho_d”: density matrix of Diagonal Ensemble.

• “Obs…”: infinite-time expectation of observable Obs.

• “delta_t_Obs…”: infinite-time temporal fluctuations of Obs.

• “delta_q_Obs…”: infinite-time quantum fluctuations of Obs.

• “Sd…” (“Sd_Renyi…” for \(\alpha\neq1.0\)): Renyi diagonal entropy of density matrix of rho_d with parameter alpha.

• “Srdm…” (“Srdm_Renyi…” for \(\alpha\neq1.0\)): Renyi entanglement entropy of reduced DM of`rho_d` (rho_d is a mixed DM itself) with parameter alpha.

Replace “…” above by ‘pure’, ‘thermal’ or ‘mixed’ depending on input parameters.

Examples

We prepare a quantum system in an eigenstate \(\psi_1\) of the Hamiltonian \(H_1=\sum_j hS^x_j + g S^z_j\). At time \(t=0\) we quench to the Hamiltonian \(H_2=\sum_j JS^z_{j+1}S^z_j+ hS^x_j + g S^z_j\), and evolve the initial state \(\psi_1\) with it. We compute the infinite-time (i.e. Diagonal Ensemble) expectation value of the Hamiltonian \(H_1\), and it’s infinite-time temporal fluctuations \(\delta_t\mathcal{O}^\psi_d\) (see above for the definition).

 1from quspin.basis import spin_basis_1d  # Hilbert space spin basis
 2from quspin.tools.measurements import diag_ensemble
 3import numpy as np  # generic math functions
 4
 5#
 6L = 12  # syste size
 7# coupling strenghts
 8J = 1.0  # spin-spin coupling
 9h = 0.8945  # x-field strength
10g = 0.945  # z-field strength
11# create site-coupling lists
12J_zz = [[J, i, (i + 1) % L] for i in range(L)]  # PBC
13x_field = [[h, i] for i in range(L)]
14z_field = [[g, i] for i in range(L)]
15# create static and dynamic lists
16static_1 = [["x", x_field], ["z", z_field]]
17static_2 = [["zz", J_zz], ["x", x_field], ["z", z_field]]
18dynamic = []
19# create spin-1/2 basis
20basis = spin_basis_1d(L, kblock=0, pblock=1)
21# set up Hamiltonian
22H1 = hamiltonian(static_1, dynamic, basis=basis, dtype=np.float64)
23H2 = hamiltonian(static_2, dynamic, basis=basis, dtype=np.float64)
24# compute eigensystems of H1 and H2
25E1, V1 = H1.eigh()
26psi1 = V1[:, 14]  # pick any state as initial state
27E2, V2 = H2.eigh()
28#
29# calculate long-time (diagonal ensemble) expectations of H1 and its temporal fluctuations
30Diag_Ens = diag_ensemble(L, psi1, E2, V2, Obs=H1, delta_t_Obs=True)
31print(Diag_Ens["Obs_pure"], Diag_Ens["delta_t_Obs_pure"])