Source code for qutip.continuous_variables

# This file is part of QuTiP: Quantum Toolbox in Python.
#
#    Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
#    All rights reserved.
#
#    Redistribution and use in source and binary forms, with or without
#    modification, are permitted provided that the following conditions are
#    met:
#
#    1. Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#
#    2. Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#
#    3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
#       of its contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
#    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
#    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
#    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
#    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
#    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
#    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
#    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
#    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
#    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###############################################################################

"""
This module contains a collection functions for calculating continuous variable
quantities from fock-basis representation of the state of multi-mode fields.
"""

__all__ = ['correlation_matrix', 'covariance_matrix',
           'correlation_matrix_field', 'correlation_matrix_quadrature',
           'wigner_covariance_matrix', 'logarithmic_negativity']

from qutip.expect import expect
import numpy as np


[docs]def correlation_matrix(basis, rho=None): """ Given a basis set of operators :math:`\\{a\\}_n`, calculate the correlation matrix: .. math:: C_{mn} = \\langle a_m a_n \\rangle Parameters ---------- basis : list of :class:`qutip.qobj.Qobj` List of operators that defines the basis for the correlation matrix. rho : :class:`qutip.qobj.Qobj` Density matrix for which to calculate the correlation matrix. If `rho` is `None`, then a matrix of correlation matrix operators is returned instead of expectation values of those operators. Returns ------- corr_mat: *array* A 2-dimensional *array* of correlation values or operators. """ if rho is None: # return array of operators return np.array([[op1 * op2 for op1 in basis] for op2 in basis], dtype=object) else: # return array of expectation values return np.array([[expect(op1 * op2, rho) for op1 in basis] for op2 in basis], dtype=object)
[docs]def covariance_matrix(basis, rho, symmetrized=True): """ Given a basis set of operators :math:`\{a\}_n`, calculate the covariance matrix: .. math:: V_{mn} = \\frac{1}{2}\\langle a_m a_n + a_n a_m \\rangle - \\langle a_m \\rangle \\langle a_n\\rangle or, if of the optional argument `symmetrized=False`, .. math:: V_{mn} = \\langle a_m a_n\\rangle - \\langle a_m \\rangle \\langle a_n\\rangle Parameters ---------- basis : list of :class:`qutip.qobj.Qobj` List of operators that defines the basis for the covariance matrix. rho : :class:`qutip.qobj.Qobj` Density matrix for which to calculate the covariance matrix. symmetrized : *bool* Flag indicating whether the symmetrized (default) or non-symmetrized correlation matrix is to be calculated. Returns ------- corr_mat: *array* A 2-dimensional *array* of covariance values. """ if symmetrized: return np.array([[0.5 * expect(op1 * op2 + op2 * op1, rho) - expect(op1, rho) * expect(op2, rho) for op1 in basis] for op2 in basis], dtype=object) else: return np.array([[expect(op1 * op2, rho) - expect(op1, rho) * expect(op2, rho) for op1 in basis] for op2 in basis], dtype=object)
[docs]def correlation_matrix_field(a1, a2, rho=None): """ Calculate the correlation matrix for given field operators :math:`a_1` and :math:`a_2`. If a density matrix is given the expectation values are calculated, otherwise a matrix with operators is returned. Parameters ---------- a1 : :class:`qutip.qobj.Qobj` Field operator for mode 1. a2 : :class:`qutip.qobj.Qobj` Field operator for mode 2. rho : :class:`qutip.qobj.Qobj` Density matrix for which to calculate the covariance matrix. Returns ------- cov_mat: *array* of complex numbers or :class:`qutip.qobj.Qobj` A 2-dimensional *array* of covariance values, or, if rho=0, a matrix of operators. """ basis = [a1, a1.dag(), a2, a2.dag()] return correlation_matrix(basis, rho)
[docs]def correlation_matrix_quadrature(a1, a2, rho=None): """ Calculate the quadrature correlation matrix with given field operators :math:`a_1` and :math:`a_2`. If a density matrix is given the expectation values are calculated, otherwise a matrix with operators is returned. Parameters ---------- a1 : :class:`qutip.qobj.Qobj` Field operator for mode 1. a2 : :class:`qutip.qobj.Qobj` Field operator for mode 2. rho : :class:`qutip.qobj.Qobj` Density matrix for which to calculate the covariance matrix. Returns ------- corr_mat: *array* of complex numbers or :class:`qutip.qobj.Qobj` A 2-dimensional *array* of covariance values for the field quadratures, or, if rho=0, a matrix of operators. """ x1 = (a1 + a1.dag()) / np.sqrt(2) p1 = -1j * (a1 - a1.dag()) / np.sqrt(2) x2 = (a2 + a2.dag()) / np.sqrt(2) p2 = -1j * (a2 - a2.dag()) / np.sqrt(2) basis = [x1, p1, x2, p2] return correlation_matrix(basis, rho)
[docs]def wigner_covariance_matrix(a1=None, a2=None, R=None, rho=None): """ Calculate the Wigner covariance matrix :math:`V_{ij} = \\frac{1}{2}(R_{ij} + R_{ji})`, given the quadrature correlation matrix :math:`R_{ij} = \\langle R_{i} R_{j}\\rangle - \\langle R_{i}\\rangle \\langle R_{j}\\rangle`, where :math:`R = (q_1, p_1, q_2, p_2)^T` is the vector with quadrature operators for the two modes. Alternatively, if `R = None`, and if annihilation operators `a1` and `a2` for the two modes are supplied instead, the quadrature correlation matrix is constructed from the annihilation operators before then the covariance matrix is calculated. Parameters ---------- a1 : :class:`qutip.qobj.Qobj` Field operator for mode 1. a2 : :class:`qutip.qobj.Qobj` Field operator for mode 2. R : *array* The quadrature correlation matrix. rho : :class:`qutip.qobj.Qobj` Density matrix for which to calculate the covariance matrix. Returns ------- cov_mat: *array* A 2-dimensional *array* of covariance values. """ if R is not None: if rho is None: return np.array([[0.5 * np.real(R[i, j] + R[j, i]) for i in range(4)] for j in range(4)], dtype=object) else: return np.array([[0.5 * np.real(expect(R[i, j] + R[j, i], rho)) for i in range(4)] for j in range(4)], dtype=object) elif a1 is not None and a2 is not None: if rho is not None: x1 = (a1 + a1.dag()) / np.sqrt(2) p1 = -1j * (a1 - a1.dag()) / np.sqrt(2) x2 = (a2 + a2.dag()) / np.sqrt(2) p2 = -1j * (a2 - a2.dag()) / np.sqrt(2) return covariance_matrix([x1, p1, x2, p2], rho) else: raise ValueError("Must give rho if using field operators " + "(a1 and a2)") else: raise ValueError("Must give either field operators (a1 and a2) " + "or a precomputed correlation matrix (R)")
[docs]def logarithmic_negativity(V): """ Calculate the logarithmic negativity given the symmetrized covariance matrix, see :func:`qutip.continous_variables.covariance_matrix`. Note that the two-mode field state that is described by `V` must be Gaussian for this function to applicable. Parameters ---------- V : *2d array* The covariance matrix. Returns ------- N: *float*, the logarithmic negativity for the two-mode Gaussian state that is described by the the Wigner covariance matrix V. """ A = V[0:2, 0:2] B = V[2:4, 2:4] C = V[0:2, 2:4] sigma = np.linalg.det(A) + np.linalg.det(B) - 2 * np.linalg.det(C) nu_ = sigma / 2 - np.sqrt(sigma ** 2 - 4 * np.linalg.det(V)) / 2 if nu_ < 0.0: return 0.0 nu = np.sqrt(nu_) lognu = -np.log(2 * nu) logneg = max(0, lognu) return logneg