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