Gaussian State¶
- class GaussianState(d: int, connector: BaseConnector, config: Config | None = None)¶
Class to represent a Gaussian state.
- Parameters:
d (int) – The number of modes.
- validate() None¶
Validates the state.
- Raises:
InvalidState – Raised if the underlying Gaussian state is invalid, which could mean - ill-shaped mean and covariance; - non-symmetric covariance matrix; - the covariance matrix doesn’t fulfill the Robertson-Schrödinger uncertainty relations.
- property xxpp_mean_vector: ndarray¶
The state’s mean in the xxpp-ordered basis.
The expectation value of the quadrature operators in xxpp basis, i.e. \(\operatorname{Tr} \left ( \rho Y \right )\), where \(Y = (x_1, \dots, x_d, p_1, \dots, p_d)^T\).
- Returns:
A \(2d\)-vector where \(d\) is the number of modes.
- Return type:
- property xxpp_covariance_matrix: ndarray¶
The xxpp-ordered covariance matrix of the state.
The xxpp-ordered covariance matrix \(\sigma_{xp}\) is defined by
\[\sigma_{xp, ij} := \langle Y_i Y_j + Y_j Y_i \rangle_\rho - 2 \langle Y_i \rangle_\rho \langle Y_j \rangle_\rho,\]where
\[Y = (x_1, \dots, x_d, p_1, \dots, p_d)^T,\]and \(\rho\) is the density operator of the currently represented state.
- Returns:
The \(2d \times 2d\) xp-ordered covariance matrix.
- Return type:
- property xxpp_correlation_matrix: ndarray¶
The state’s correlation matrix in the xxpp basis.
Let \(M_{(xxpp)}\) be the correlation matrix in the xxpp basis. Then
\[M_{ij (xxpp)} = \langle Y_i Y_j + Y_j Y_i \rangle_\rho,\]where \(M_{ij (xxpp)}\) denotes the matrix element in the \(i\)-th row and \(j\)-th column,
\[Y = (x_1, \dots, x_d, p_1, \dots, p_d)^T,\]and \(\rho\) is the density operator of the currently represented state.
- Returns:
The \(2d \times 2d\) correlation matrix in the xxpp-basis.
- Return type:
- property xxpp_representation: Tuple[ndarray, ndarray]¶
The state’s mean and correlation matrix ordered in the xxpp basis.
- Returns:
- Return type:
- property xpxp_mean_vector: ndarray¶
Returns the xpxp-ordered mean of the state.
- Returns:
A \(2d\)-vector. The expectation value of the quadrature operators in xxpp-ordering, i.e. \(\operatorname{Tr} \rho R\), where \(R = (x_1, p_1, \dots, x_d, p_d)^T\).
- Return type:
- property xpxp_covariance_matrix: ndarray¶
The xpxp-ordered covariance matrix of the state.
The xpxp-ordered covariance matrix \(\sigma\) is defined by
\[\sigma_{ij} := \langle R_i R_j + R_j R_i \rangle_\rho - 2 \langle R_i \rangle_\rho \langle R_j \rangle_\rho,\]where
\[R = (x_1, p_1, \dots, x_d, p_d)^T,\]and \(\rho\) is the density operator of the currently represented state.
- Returns:
The \(2d \times 2d\) xpxp-ordered covariance matrix in xpxp-ordered basis.
- Return type:
- property xpxp_correlation_matrix: ndarray¶
The xpxp-ordered correlation matrix of the state.
Let \(M\) be the correlation matrix in the xpxp basis. Then
\[M_{ij} = \langle R_i R_j + R_j R_i \rangle_\rho,\]where \(M_{ij}\) denotes the matrix element in the \(i\)-th row and \(j\)-th column,
\[R = (x_1, p_1, \dots, x_d, p_d)^T,\]and \(\rho\) is the density operator of the currently represented state.
- Returns:
The \(2d \times 2d\) xpxp-ordered correlation matrix.
- Return type:
- property xpxp_representation: Tuple[ndarray, ndarray]¶
The state’s mean and correlation matrix ordered by the xpxp basis.
- Returns:
mean(),corr().- Return type:
- property complex_displacement: ndarray¶
The complex displacement of the state.
The complex displacement is defined by
\[\mu_{c} := \begin{bmatrix} \langle a_1 \rangle_{\rho}, \dots, \langle a_d \rangle_{\rho}, \langle a^\dagger_1 \rangle_{\rho}, \dots, \langle a^\dagger_d \rangle_{\rho} \end{bmatrix}.\]Equivalently, one can write
\[\mu_{c} := W \mu_{xxpp},\]where \(\mu_{xxpp}\) is the xxpp-ordered mean vector
xxpp_mean_vectorand\[\begin{split}W = \frac{1}{\sqrt{2}} \begin{bmatrix} I_{d} & i I_{d} \\ I_{d} & -i I_{d} \end{bmatrix}.\end{split}\]- Returns:
The complex displacement.
- Return type:
- property complex_covariance: ndarray¶
The complex covariance of the state.
The complex covariance is defined by
\[\sigma_{c, ij} = \langle \xi_i \xi^\dagger_j + \xi_j \xi^\dagger_i \rangle_{\rho} - 2 \langle \xi_i \rangle_{\rho} \langle \xi^\dagger_j \rangle_{\rho},\]where
\[\xi = \begin{bmatrix} a_1, \dots a_d, a^\dagger_1, \dots, a^\dagger_d \end{bmatrix}.\]Equivalently, one can write
\[\sigma_{c} = \frac{1}{\hbar} W \sigma_{xxpp} W^{\dagger},\]where \(\sigma_{xxpp}\) is the xxpp-ordered covariance matrix
xxpp_covand\[\begin{split}W = \frac{1}{\sqrt{2}} \begin{bmatrix} I_{d} & i I_{d} \\ I_{d} & -i I_{d} \end{bmatrix}.\end{split}\]- Returns:
The complex covariance.
- Return type:
- rotated(phi: float) GaussianState¶
Returns a copy of the current Gaussian state, rotated by an angle phi.
Let \(\phi \in [ 0, 2 \pi )\). Let us define the following:
\[x_{i, \phi} = \cos\phi~x_i + \sin\phi~p_i,\]which is a generalized quadrature operator. One could rotate the state by a simple complex phase, which can be shown by using the transformation rules between the ladder operators and quadrature operators, i.e.
\[\begin{split}x_i &= \sqrt{\frac{\hbar}{2}} (a_i + a_i^\dagger) \\ p_i &= -i \sqrt{\frac{\hbar}{2}} (a_i - a_i^\dagger),\end{split}\]which we could rewrite \(x_{i, \phi}\) to
\[x_{i, \phi} = \sqrt{\frac{\hbar}{2}} \left( a_i \exp(-i \phi) + a_i^\dagger \exp(i \phi) \right),\]meaning that e.g. the annihilation operators \(a_i\) are transformed just by multiplying it with a phase factor \(\exp(-i \phi)\) i.e.
\[\begin{split}(\langle a_i \rangle_{\rho} =: )~m_i &\mapsto \exp(-i \phi) m_i \\ (\langle a^\dagger_i a_j \rangle_{\rho} =: )~C_{ij} &\mapsto C_{ij} \\ (\langle a_i a_j \rangle_{\rho} =: )~G_{ij} &\mapsto \exp(-i 2 \phi) G_{ij}.\end{split}\]- Parameters:
phi (float) – The angle to rotate the state with.
- Returns:
The rotated GaussianState instance.
- Return type:
- reduced(modes: Tuple[int, ...]) GaussianState¶
Returns a copy of the current state, reduced to the given modes.
This method essentially preserves the modes specified from the representation of the Gaussian state, but cuts out the other modes.
- Parameters:
- Returns:
The reduced GaussianState instance.
- Return type:
- xpxp_reduced_rotated_mean_and_covariance(modes: Tuple[int, ...], phi: float) Tuple[ndarray, ndarray]¶
The mean and covariance on a rotated and reduced state.
Let the index set \(\vec{i}\) correspond to modes, and the angle \(\phi\) correspond to phi. The current
GaussianStateinstance is reduced to modes and rotated by phi in a new instance, and let that state be denoted by \(\rho_{\vec{i}, \phi}\).Then the xpxp-ordered mean and covariance can be calculated by
\[\begin{split}\mu_{\vec{i}, \phi} &:= \langle R_{\vec{i}} \rangle_{\rho_{\vec{i}, \phi}}, \\ \sigma_{\vec{i}, \phi} &:= \langle R_{\vec{i}} R_{\vec{i}}^T \rangle_{\rho_{\vec{i}, \phi}} - \mu_{\vec{i}, \phi} \mu_{\vec{i}, \phi}^T,\end{split}\]where
\[R = (x_1, p_1, \dots, x_d, p_d)^T,\]and \(R_{\vec{i}}\) is just the same vector, reduced to a subsystem specified by \(\vec{i}\).
- Parameters:
- Returns:
xpxp-ordered mean and covariance of the reduced and rotated version of the current
GaussianState.- Return type:
- mean_photon_number(modes: Tuple[int, ...] | None = None) float¶
This method returns the mean photon number of the given modes. The mean photon number \(\bar{n} = \langle \hat{n} \rangle\) can be calculated in terms of the ladder operators by the following expression
\[\sum_{i=1}^{d} \langle a_{i}^\dagger a_{i} \rangle = \operatorname{Tr}(\rho \hat{n}),\]where \(a\), \(a ^\dagger\) are the annihilation and the creation operators respectively, \(\rho\) is the density operator of the currently represented state and \(d\) is the number of modes. For a general displaced squeezed gaussian state, the mean photon number is
\[\langle \hat{n} \rangle = \operatorname{Tr}(\langle a^\dagger a \rangle) + \mu_{c}^ \dagger \cdot \mu_{c},\]where \(\mu_{c}\) is the
complex_displacement.Note
This method can also be used to return the summation of the mean photon number for multiple modes if the mode parameter contains more than one integer e.g \((0,1,...)\).
- variance_photon_number(modes: Tuple[int, ...] | None = None) float¶
This method calculates the variance of the photon number operator as follows:
\[\operatorname{Var}(\hat{n}) = \mathbb{E}((\hat{n} - \bar{n})^2),\]where \(\bar{n}\) is the expectation value of the photon number, given by
mean_photon_number().See, e.g., [Means and covariances of photon numbers in multimode Gaussian states](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.99.023817) for details on how to calculate this for Gaussian states.
- fidelity(state: GaussianState) float¶
Calculates the state fidelity between two quantum states.
The state fidelity \(F\) between two density matrices \(\rho_1, \rho_2\) is given by:
\[\operatorname{F}(\rho_1, \rho_2) = \left [ \operatorname{Tr} \sqrt{\sqrt{\rho_1} \rho_2 \sqrt{\rho_1} } \right ]\]A Gaussian state can be represented by its covariance matrix and displacement vector.
Let \(\mu_1, \mu_2\) be the displacement vectors and \(\sigma_1, \sigma_2\) be the covariance matrices of the \(\rho_1, \rho_2\) Gaussian states, respectively. Define \(\hat{\sigma} = \frac{\sigma_1 + \sigma_2}{2}\) and \(\Delta\mu = \mu_2 - \mu_1\). The fidelity can be written as
\[\operatorname{F}(\rho_1, \rho_2) = \operatorname{F_0}(\sigma_1, \sigma_2) \exp \left( -\frac{1}{4} \Delta\mu^T (\hat{\sigma})^{-1} \Delta\mu \right),\]where \(F_0\) can be written as
\[\operatorname{F_0} = \frac{ \prod_{i=1}^d \left [w_i + \sqrt{w_i^2 - 1} \right] }{ \sqrt{\det \hat{\sigma}} },\]where \(w_i \geq 1\) and \(\pm w_i\) are the eigenvalues of the matrix
\[W := - \frac{i}{2} \Omega^T \hat{\sigma}^{-1} \left( I - \sigma_2 \Omega \sigma_1 \Omega \right)\]and
\[\begin{split}\Omega = \begin{bmatrix} 0 & 1 \\-1 & 0 \end{bmatrix} \otimes I_{2d \times 2d}.\end{split}\]References
Note
In this notation \(\sqrt{\hbar} \mu_i\) is equivalent to
xpxp_mean_vectorand \(\hbar \sigma_i\) is equivalent toxpxp_covariance_matrix.- Parameters:
state – Another
GaussianStateinstance.- Returns:
The calculated fidelity.
- Return type:
- quadratic_polynomial_expectation(A: ndarray, b: ndarray, c: float = 0.0, phi: float = 0.0) float¶
The expectation value of the specified quadratic polynomial.
A quadratic polynomial can be written as
\[f(R) = R^T A R + R \cdot b + c,\]where \(R = (x_1, p_1, \dots, x_d, p_d)^T\) is the vector of the quadrature operators where \(d\) is the number of modes, \(A \in \mathbb{R}^{2d \times 2d}\) is a symmetric matrix, \(b \in \mathbb{R}^{2d}\), and \(c\in\mathbb{R}\).
This method returns the expectation value \(E[f(R)]\) using the following equation:
\[\operatorname{E}[f(R)] = \operatorname{Tr}[ A\sigma ] + \mu^T A \mu + \mu^T b + c,\]where \(\sigma\) is the covariance matrix, \(\mu = E[R]\) is the mean of the quadrature operators and
\[\operatorname{E}[\cdot] = \operatorname{Tr}[\cdot \rho],\]where \(\rho\) is the density matrix of the state.
- Parameters:
A (numpy.ndarray) – A \(2d \times 2d\) real symmetric matrix corresponding to the quadratic coefficients, where \(d\) is the number of modes.
b (numpy.ndarray) – A one-dimensional \(2d\)-length real-valued vector that corresponds to the first order terms of the quadratic polynomial.
c (float) – The constant term in the quadratic polynomial. Defaults to 0.
phi (float) – Rotation angle, by which the state is rotated. Defaults to 0.
- Returns:
The expectation value of the quadratic polynomial.
- Return type:
- wigner_function(positions: List[List[float]], momentums: List[List[float]], modes: Tuple[int, ...] | None = None) ndarray¶
This method calculates the Wigner function values at the specified position and momentum vectors, according to the following equation:
\[W(r) = \frac{1}{\pi^d \sqrt{\mathrm{det} \sigma}} \exp \big ( - (r - \mu)^T \sigma^{-1} (r - \mu) \big ).\]- Parameters:
- Returns:
The Wigner function values in the shape of a grid specified by the input.
- Return type:
- plot_wigner(positions: List[List[float]], momentums: List[List[float]], mode: int | None = None) None¶
Plots the Wigner function in phase space for a single mode using the given positions and momentums.
- Parameters:
Note
Only a single mode is supported.
- Returns:
This method generates a plot and does not return a value.
- Return type:
None
- get_particle_detection_probability(occupation_number: ndarray) float¶
Returns the particle number detection probability using the occupation number specified as a parameter.
- get_threshold_detection_probability(occupation_number: ndarray | Tuple[int, ...]) float¶
Calculates the probability of threshold particle number detection events.
The occupation number needs to be specified as consecutive 0s and 1s, where 0 corresponds to no particle detection and 1 to detecting some (>0) particles.
- Parameters:
occupation_number (Union[np.ndarray, Tuple[int, ...]]) – Array of 0s and 1s corresponding to particle detection events per mode.
- Raises:
PiquassoException – When an invalid input ‘occupation_number’ is specified.
- Returns:
The threshold particle number detection probability.
- Return type:
- property fock_probabilities: ndarray¶
Returns the particle detection probabilities.
Note
The ordering of the Fock basis is increasing with particle numbers, and in each particle number conserving subspace, lexicographic ordering is used.
- Returns:
The particle detection probabilities.
- Return type:
- get_marginal_fock_probabilities(modes: Tuple[int, ...]) ndarray¶
Return the particle number probabilities on a given subsystem.
- get_xp_string_moment(string: List[int]) complex¶
Moment corresponding to a product of quadrature operators (XP-string).
The indices of string correspond to the xxpp-ordering, i.e., \(0, \dots, d-1\) denotes \(\hat{x}_1, \dots \hat{x}_d\) and \(d, \dots, 2d-1\) denotes \(\hat{p}_0, \dots, \hat{p}_d\).
We use the following formula:
\[\begin{split}\left\langle Y_{\alpha_1} Y_{\alpha_2} \cdots Y_{\alpha_n} \right\rangle = \sum_{\substack{\text{all pairings}\\ \text{and singletons}}} \Biggl( \prod_{\text{single }r} \mu_{\alpha_r} \Biggr) \Biggl( \prod_{\text{pairing }(i,j)} V_{\alpha_i, \alpha_j} \Biggr),\end{split}\]where \(Y_{\alpha_i}\) are quadrature operators in xxpp-ordering, \(\mu\) is the mean vector from
xxpp_mean_vector(), \(\Sigma\) is the covariance matrix fromxxpp_covariance_matrix(), and\[\begin{split}V &= \Sigma + \frac{i\hbar}{2} \Omega, \\ \Omega &= \begin{bmatrix} 0 & I \\ -I & 0 \end{bmatrix},\end{split}\]where \(\Omega\) encodes the canonical commutation relations.
The reasoning for this formula goes as follows: suppose that the characteristic function giving the Weyl ordering of the quadrature operators is \(\chi_W(\mathbf{x}, \mathbf{p})\). Then, writing the Weyl operator in terms of \(\hat{x}_i\) and \(\hat{p}_i\) and using the Baker-Campbell-Hausdorff formula, we get that for normal ordering (\(\hat{x}_i\) always in the front) we need to use \(\chi_W(\mathbf{x}, \mathbf{p}) e^{i \mathbf{x}\cdot\mathbf{p}}\), which induces a shift in the covariance matrix \(\Sigma_{\text{normal}} := \Sigma + \frac{i}{2}X\), where
\[\begin{split}X = \begin{bmatrix} 0 & I \\ I & 0 \end{bmatrix}.\end{split}\]By differentiating the characteristic function correspondng to the normal order, one can get the normal-ordered moments. Hence, inserting \(\Sigma_{\text{normal}}\) as the covariance matrix in Isserlis’ formula, one would get the value for the normal ordered moment. However, the terms form \((\Sigma_{\text{normal}})_{i+d, i}\) do not appear (which is also equal to \((\Sigma_{\text{normal}})_{i, i+d}\) due to symmetry). Therefore, we can just encode the canonical commutation relations by applying the following transformation:
\[(\Sigma_{\text{normal}})_{i+d, i} \mapsto = (\Sigma_{\text{normal}})_{i, i+d} - i \hbar.\]This yields a formula which works for normal ordered moments, but also incorporates the commutation relations, therefore it must work for any ordering.
Furthermore, by differentiating the normal-ordered characteristic function, we get—as an extension to the centered Isserlis’ theorem—the following formula:
\[\begin{split}\left\langle Y_{\alpha_1} Y_{\alpha_2} \cdots Y_{\alpha_n} \right\rangle = \sum_{\substack{\text{all pairings}\\ \text{and singletons}}} \Biggl( \prod_{\text{single }r} \mu_{\alpha_r} \Biggr) \Biggl( \prod_{\text{pairing }(i,j)} (\Sigma_\text{normal})_{\alpha_i, \alpha_j} \Biggr),\end{split}\]and here, one can just replace \(\Sigma_\text{normal}\) with \(V\) by the previous argument to obtain the formula for all orderings.
- Sources:
- purify() GaussianState¶
Returns a purification corresponding to the current state on a double-sized system.
Note, that the purification results in a state where the mean vector on the auxiliary system is identical with the mean vector in the original system.
- Returns:
The purified Gaussian state.
- Return type:
- get_parity_operator_expectation_value() float¶
Expectation value of the parity operator.
The parity operator \(\hat{\Pi}\) is defined as
\[\hat{\Pi} = \exp \left( i \pi \sum_{i=1}^{d} a_i^\dagger a_i \right),\]where \(a_i\), \(a_i^\dagger\) are the annihilation and creation operators respectively, and \(d\) is the number of modes.
The expectation value of the parity operator can be calculated for Gaussian states by the following formula:
\[\langle \hat{\Pi} \rangle = \frac{ \exp \left( - \mu^\dagger \sigma^{-1} \mu \right) }{ \sqrt{\det \sigma} },\]where \(\mu\) is the mean vector
complex_displacementand \(\sigma\) is the covariance matrixcomplex_covariance.- Returns:
The expectation value of the parity operator.
- Return type: