Quantum-Enhanced Imaging Model Based on Squeezed States

1. Introduction

Optical imaging technology plays a crucial role in humanity’s exploration of the unknown world. However, the resolution of traditional imaging systems is limited by the Rayleigh diffraction limit. Based on quantum Fisher information theory, Tsang et al. were the first to investigate the fundamental limits of imaging resolution and the method for achieving two-point super-resolution imaging [1]. Subsequently, research on super-resolution imaging grounded in quantum information theory has been extensively developed, evolving from initial one-dimensional two-point imaging to the estimation of multiple points and three-dimensional imaging [2,3,4]. This research has further extended to practical imaging scenarios involving factors such as detection noise and mode crosstalk [5,6,7]. Currently, super-resolution quantum imaging has become a focal point of research. Nevertheless, existing studies remain primarily concentrated in the field of imaging with incoherent thermal light sources.

Squeezed states serve as a crucial quantum resource for enhancing the precision of quantum measurements [8], exemplified by their application in gravitational-wave laser interferometers [9]. With the continuous advancement of research, the application of squeezed states in quantum imaging has emerged as a new and prominent research direction. This is evidenced by studies in areas such as Raman microscopy imaging [10,11,12] and stimulated Brillouin scattering imaging [13]. These studies demonstrate that quantum imaging based on squeezed states can effectively overcome the shot noise limit and improve imaging quality.

Can squeezed states enable quantum imaging that simultaneously overcomes both the shot noise limit and the Rayleigh diffraction limit? Mikhail I. Kolobov et al. theoretically explored quantum imaging based on spatially multimode squeezed vacuum states [14]. Their results indicate that squeezed states can achieve super-resolution quantum imaging surpassing the shot noise limit. However, they did not provide the ultimate quantum theoretical limit for such imaging. Giacomo Sorelli et al. [15] investigated the super-resolution problem for point sources using squeezed states based on quantum Fisher information theory and derived a quantum limit for resolution. Nevertheless, their system model is relatively simple and cannot be readily extended to the study of arbitrary images. Therefore, how to establish a comprehensive imaging theoretical model based on squeezed states and conduct research on super-resolution imaging utilizing squeezed states remains an open and pressing issue to be addressed.

To address the aforementioned issues related to super-resolution imaging based on squeezed states, this paper constructs a corresponding imaging theory model and presents numerical calculation results. By employing discretization and quantum Fisher information theory, the problem of two-point resolution is investigated. The correctness of the established model is verified through comparison with results from existing literature. Building upon this theoretical model, the paper further explores the problem of spatial multi-parameter estimation in imaging for light fields with arbitrary spatial distributions, demonstrating the advantages of squeezed coherent states in quantum imaging. This quantum imaging model and its computational methodology can be extended to more complex scenarios, such as imaging schemes based on multimode squeezed states, and even to incoherent imaging systems [16]. It thus lays a theoretical foundation for future experimental research on super-resolution imaging utilizing squeezed states.

2. Theoretical Methods

2.1. Description of the Imaging System

We consider a simple imaging system as shown in Figure 1. The multimode optical field at the object plane, represented by ${\widehat{E}}_{\mathrm{Source}}\left(x\right)={\widehat{\mathbf{a}}}^{\left(u\right)\mathrm{T}}·\mathbf{u}\left(x\right)$, is composed of a finite number of modes $\left\{{\mathbf{u}}_j\left(x\right)\right\}$. This field corresponds to the optical modes excited at the object plane (typically resulting from the reflection or transmission of an illumination source), where ${\widehat{\mathbf{a}}}^{\left(u\right)}={[{\widehat{a}}_0^{\left(u\right)},{\widehat{a}}_1^{\left(u\right)},\dots ]}^{\mathrm{T}}$ denotes the annihilation operators of the optical field, and $\mathbf{u}\left(x\right)={[u_0\left(x\right),u_1\left(x\right),\dots ]}^{\mathrm{T}}$ represents the corresponding spatial modes of the field. It carries the desired complex amplitude distribution information, which is the target to be reconstructed from the measurement results. In addition to ${\widehat{E}}_{\mathrm{Source}}\left(x\right)$, the vacuum field $\widehat{\mathbf{b}}$ is also present at the object plane. The complete optical field at the object plane, ${\widehat{E}}_{\mathrm{Object}}\left(x\right)$ is the superposition of the two. In our imaging model, ${\widehat{\mathbf{E}}}_{\mathrm{Object}}\left(x\right)$ is expanded using a set of modes $\left\{{\varphi }_j\left(\mathbf{x}\right)\right\}$ as ${\widehat{E}}_{\mathrm{Object}}\left(\mathbf{x}\right)={\widehat{\mathbf{a}}}^{\left(\varphi \right)\mathrm{T}}·\varphi \left(x\right)$, where the relationship between ${\widehat{\mathbf{a}}}^{\left(\varphi \right)}$ and ${\widehat{\mathbf{a}}}^{\left(u\right)}$ is given by

$${\widehat{\mathbf{a}}}^{\left(\varphi \right)}=X_1{\widehat{\mathbf{a}}}^{\left(u\right)}+Y_1\widehat{\mathbf{b}}$$

(1)

$$X_1=\left[\begin{array}{ccc}\langle u_0\left(x\right)\mid {\varphi }_0\left(x\right)\rangle & \langle u_1\left(x\right)\mid {\varphi }_0\left(x\right)\rangle & \dots \\ \langle u_0\left(x\right)\mid {\varphi }_1\left(x\right)\rangle & \langle u_1\left(x\right)\mid {\varphi }_1\left(x\right)\rangle & \\ ⋮& & \ddots \end{array}\right]$$

(2)

The expression $Y_1\widehat{\mathbf{b}}$ refers to Appendix A. For scalar, quasi-monochromatic classical imaging problems, the propagation of the optical field’s complex amplitude from the object plane to the image plane is regarded as the action of an integral operator with the amplitude point spread function $g(x,y)$ as its kernel on $E_{\mathrm{Source}}\left(x\right)$, yielding the optical field at the image plane.

$$E_{\mathrm{Image}}\left(y\right)=\underset{\mathbb{S}}{\int \int }g(x,y)E_{\mathrm{Source}}\left(x\right)dx$$

(3)

Where $\mathbb{S}$ represents the distribution area of the optical field. The quantum version of this form was provided by Jeffrey H. Shapiro [17]. Through singular value decomposition, the amplitude point spread function can be expressed as

$$g(x,y)=\sum _j\sqrt{{\eta }_j}{\Phi }_j\left(x\right){\Psi }_j\left(y\right)$$

(4)

Then, the diffraction process can be described as a transformation from the vector of annihilation operators ${\widehat{\mathbf{a}}}^{(\Phi )}$ corresponding to the modes $\left\{{\Phi }_j\left(x\right)\right\}$ to the vector of annihilation operators ${\widehat{\mathbf{a}}}^{(\Psi )}$ corresponding to the modes $\left\{{\Psi }_j\left(y\right)\right\}$, given by ${\widehat{\mathbf{a}}}^{(\Psi )}=X_d{\widehat{\mathbf{a}}}^{(\Phi )}+Y_d{\widehat{\mathbf{b}}}^{\prime }$, where

$$X_d=\left[\begin{array}{ccc}\sqrt{{\eta }_0}& & \\ & \sqrt{{\eta }_1}& \\ & & \ddots \end{array}\right]$$

(5)

Where ${\widehat{\mathbf{b}}}^{\prime }$ represents the vacuum state introduced due to diffraction losses. In practical computations, the sets $\left\{{\Phi }_j\left(x\right)\right\}$ and $\left\{{\Psi }_j\left(y\right)\right\}$ cannot be obtained directly. To perform the singular value decomposition, we first need to select a known set of basis functions for the object plane and image plane, $\left\{{\varphi }_k\left(x\right)\right\}$ and $\left\{{\psi }_j\left(y\right)\right\}$, and express $g(x,y)$ in matrix form. Based on the calculation results, we then reselect $\left\{{\varphi }_k\left(x\right)\right\}$ and $\left\{{\psi }_j\left(y\right)\right\}$ into the desired form. A suitable choice is the rectangular function basis. For our one-dimensional problem, the rectangular function basis $\left\{R_b(x;x_k,\Delta )\right\}$ is defined as

$$R_b(x;x_k,\Delta )=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\sqrt{\Delta }}},\hfill & x_k-{\displaystyle \frac{\Delta }{2}}\le x<x_k+{\displaystyle \frac{\Delta }{2}}\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(6)

Here, $\Delta$ is the width of the rectangle. As $\Delta \to 0$, we have

$$\langle u_j\left(x\right)\mid R_b(x;x_k,\Delta )\rangle \approx \sqrt{\Delta }{u}_j\left(x_k\right)$$

(7)

$$\langle g(x,y)\mid R_b(x;x_k,\Delta )R_b(y;y_j,\Delta )\rangle \approx \Delta k(x_k,y_j)$$

(8)

Its core advantage lies in enabling efficient inner product operations to avoid numerical integration, thereby improving computational efficiency. We first expand $g(x,y)$ using $\left\{{\varphi }_k\left(x\right)\right\}=\left\{R_b(x;x_k,\Delta )\right\}$ and $\left\{{\psi }_j\left(y\right)\right\}=\left\{R_b(y;y_j,\Delta )\right\}$ to obtain the matrix

$$g=\left[\begin{array}{ccc}\Delta g(x_0,y_0)& \Delta g(x_1,y_0)& \dots \\ \Delta g(x_0,y_1)& \Delta g(x_1,y_1)& \\ ⋮& & \ddots \end{array}\right]$$

(9)

Performing a truncated singular value decomposition on it, we obtain $g=U\Lambda V^{†}$ and discard the smaller singular values. The basis functions $\left\{{\Psi }_j\left(y\right)\right\}$ corresponding to the right singular vectors are approximated as ${\Psi }_j\left(y\right)={\sum }_k{U}_{kj}^{*}{R}_b(y;y_k,\Delta )$, where $U_{kj}^{*}$ represents the matrix elements of U. At this point, the corresponding annihilation operator vector ${\widehat{\mathbf{a}}}^{(\Psi )}$ has a lower dimensionality but contains the vast majority of the information from the image plane optical field. We use $\left\{{\varphi }_j\left(x\right)\right\}=\left\{R_b(x;x_j,\Delta )\right\}$ as the basis functions for expanding the object plane field to reduce the computational complexity of projection operations. The corresponding annihilation operator vector is ${\widehat{\mathbf{a}}}^{\left(\varphi \right)}$. The image plane basis functions are expanded using $\left\{{\psi }_j\left(y\right)\right\}=\left\{{\Psi }_j\left(y\right)\right\}$ to reduce the scale of the output information. The diffraction process is then expressed as the transformation from the annihilation operator vector ${\widehat{\mathbf{a}}}^{\left(\varphi \right)}$ corresponding to $\left\{{\varphi }_k\left(x\right)\right\}$ to ${\widehat{\mathbf{a}}}^{(\Psi )}$, given by ${\widehat{\mathbf{a}}}^{(\Psi )}=X_2{\widehat{\mathbf{a}}}^{\left(\varphi \right)}+Y_2{\widehat{\mathbf{b}}}^{\prime }$, where

$$X_2=\Lambda V^{†}$$

(10)

Here, $Y_2{\widehat{\mathbf{b}}}^{\prime }$ represents the vacuum state introduced due to losses (see Appendix A). The optical field at the image plane is obtained as ${\widehat{E}}_{\mathrm{Object}}\left(x\right)={\widehat{\mathbf{a}}}^{(\Psi )\mathrm{T}}·\Psi \left(x\right)$. The overall transformation from the source field ${\widehat{E}}_{\mathrm{Source}}\left(x\right)$ to the image plane field ${\widehat{E}}_{\mathrm{Object}}\left(x\right)$ is described as follows.

$${\widehat{\mathbf{a}}}^{(\Psi )}=X_2(X_1{\widehat{\mathbf{a}}}^{\left(u\right)}+Y_1\widehat{\mathbf{b}})+Y_2{\widehat{\mathbf{b}}}^{\prime }$$

(11)

2.2. Gaussian States in Linear Systems

In quantum optics, the quadrature amplitude operator ${\widehat{q}}_j$ and quadrature phase operator ${\widehat{p}}_j$ of an optical field are defined through the annihilation operator as ${\widehat{q}}_j={\displaystyle \frac{{\widehat{a}}_j+{\widehat{a}}_j^{†}}{\sqrt{2}}}$ and ${\widehat{p}}_j={\displaystyle \frac{{\widehat{a}}_j-{\widehat{a}}_j^{†}}{\sqrt{2}\mathrm{i}}}$. The operator vector for a multimode Gaussian state is defined as $\widehat{\mathbf{d}}={[{\widehat{q}}_0,{\widehat{p}}_0,{\widehat{q}}_1,{\widehat{p}}_1,{\widehat{p}}_1,\dots ]}^{\mathrm{T}}$. An important property of Gaussian states is that all statistical moments of the probability distributions for measurements of their quadrature amplitudes and phases can be completely described by their expectation value $\langle \widehat{\mathbf{d}}\rangle$ and covariance matrix $\mathbb{V}=\langle \{\widehat{\mathbf{d}}-\langle \widehat{\mathbf{d}}\rangle ,{\widehat{\mathbf{d}}}^{\mathrm{T}}-{\langle \widehat{\mathbf{d}}\rangle }^{\mathrm{T}}\}\rangle$. Therefore, these contain all the statistical information of the Gaussian state [18]. For a transformation

$${\widehat{\mathbf{a}}}_{\mathrm{out}}=X{\widehat{\mathbf{a}}}_{\mathrm{in}}+Y\widehat{\mathbf{b}}$$

(12)

where X represents the transformation acting on the non-vacuum input ${\widehat{\mathbf{a}}}_{\mathrm{in}}$, and Y is the transformation acting on the vacuum input. The following relationships hold:

$$\begin{array}{cc}\hfill {\langle \widehat{\mathbf{a}}\rangle }_{\mathrm{out}}& =X{\langle \widehat{\mathbf{a}}\rangle }_{\mathrm{in}}\hfill \\ \hfill {\langle {\widehat{\mathbf{a}}}^{†}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{out}}& =\overline{X}{\langle {\widehat{\mathbf{a}}}^{†}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{in}}{X}^{\mathrm{T}}\hfill \\ \hfill {\langle \widehat{\mathbf{a}}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{out}}& =X{\langle \widehat{\mathbf{a}}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{in}}{X}^{\mathrm{T}}\hfill \end{array}$$

(13)

where ⊗ denotes the Kronecker product. Both $\langle \widehat{\mathbf{d}}\rangle$ and $\mathbb{V}$ can be expressed in terms of $\langle \widehat{\mathbf{a}}\rangle$, $\langle {\widehat{\mathbf{a}}}^{†}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle$, and $\langle \widehat{\mathbf{a}}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle$:

$$\begin{array}{c}\begin{array}{cc}\hfill \langle \widehat{\mathbf{d}}\rangle & ={\displaystyle \frac{\langle \widehat{\mathbf{a}}\rangle +{\langle \widehat{\mathbf{a}}\rangle }^{*}}{\sqrt{2}}}\otimes \left[\begin{array}{c}1\\ 0\end{array}\right]+{\displaystyle \frac{\langle \widehat{\mathbf{a}}\rangle -{\langle \widehat{\mathbf{a}}\rangle }^{*}}{\sqrt{2}\mathrm{i}}}\otimes \left[\begin{array}{c}0\\ 1\end{array}\right]\hfill \end{array}\\ \begin{array}{cc}\hfill \mathbb{V}& =2\mathrm{Re}\left[\langle {\widehat{\mathbf{a}}}^{†}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle \right]\otimes \left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]+2\mathrm{Im}\left[\langle {\widehat{\mathbf{a}}}^{†}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle \right]\otimes \left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\hfill \\ \hfill & +2\mathrm{Re}\left[\langle \widehat{\mathbf{a}}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle \right]\otimes \left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]+2\mathrm{Im}\left[\langle \widehat{\mathbf{a}}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle \right]\otimes \left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\hfill \\ \hfill & -4\langle \widehat{\mathbf{d}}\rangle \otimes {\langle \widehat{\mathbf{d}}\rangle }^{\mathrm{T}}+I\hfill \end{array}\end{array}$$

(14)

where I denotes the identity matrix. Moreover, Equation (14) is invertible; its inverse operation is provided in Appendix B. Therefore, for a multi-mode Gaussian state, the calculation proceeds as follows: first, transform the input state’s ${\langle \widehat{\mathbf{d}}\rangle }_{\mathrm{in}}$ and ${\mathbb{V}}_{\mathrm{in}}$ into ${\langle \widehat{\mathbf{a}}\rangle }_{\mathrm{in}}$, ${\langle {\widehat{\mathbf{a}}}^{†}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{in}}$ and ${\langle \widehat{\mathbf{a}}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{in}}$ using the inverse of Equation (14). Next, compute the output state’s ${\langle \widehat{\mathbf{a}}\rangle }_{\mathrm{out}}$, ${\langle {\widehat{\mathbf{a}}}^{†}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{out}}$ and ${\langle \widehat{\mathbf{a}}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{out}}$ via Equation (13). Finally, convert these results into ${\langle \widehat{\mathbf{d}}\rangle }_{\mathrm{out}}$ and ${\mathbb{V}}_{\mathrm{out}}$ using Equation (14). In this process, the calculations require only the classical transformation X from Equation (13), with no $Y\widehat{\mathbf{b}}$ term involved, thereby simplifying the analysis. For more complex systems, the overall transformation can be constructed through matrix multiplication or direct sums of multiple X matrices, thereby optimizing the computational cost.

2.3. Quantum Fisher Information for Gaussian States

The quantum Fisher information (QFI), denoted as F, is a physical quantity that quantifies the achievable precision for parameter estimation given a specific quantum state and measurement problem. It yields the quantum Cramér-Rao bound (QCRB) [18]$\mathbb{V}\left(\mathbf{t}\right)\ge {\displaystyle \frac{1}{NF\left(\mathbf{t}\right)}}$, (where N is the number of repetitions of the experiment) which represents the lower bound on the variance of an unbiased estimator for a given probe state. Consider a prepared probe state described by the density matrix ${\widehat{\rho }}_0$. After interacting with the target of interest, the probe state evolves into $\widehat{\rho }$, which now encodes information about the parameter t. The symmetric logarithmic derivative (SLD) operator ${\widehat{L}}_t$, required for computing the QFI, can then be defined and solved via the Lyapunov equation:

$${\displaystyle \frac{\partial \widehat{\rho }}{\partial t}}={\displaystyle \frac{{\widehat{L}}_t\widehat{\rho }+\widehat{\rho }{\widehat{L}}_t}{2}}$$

(15)

This leads to the QFI $F=\mathrm{Tr}\left[\widehat{\rho }{\widehat{L}}_t^2\right]$. For multi-parameter estimation problems [19], the quantum Fisher information is generalized to the quantum Fisher information matrix (QFIM), with the corresponding QCRB expressed as a matrix inequality. Specifically, let the $\mathcal{N}$-dimensional parameter vector to be estimated be $\mathbf{t}=(t_0,t_1,\dots ,t_{\mu },t_{\nu },\dots ,t_{\mathcal{N}-1})$. Then the covariance matrix of $\mathbf{t}$ satisfies

$$\mathbb{V}\left(\mathbf{t}\right)\ge {\displaystyle \frac{1}{N}}{F}^{-1}\left(\mathbf{t}\right)$$

(16)

Generally, calculating the QFIM is quite complex. For scenarios where information is encoded into Gaussian states, Rosanna Nichols et al. [19] proposed an efficient and precise computational method that determines the QFIM based on the expectation values and covariance matrix. The matrix elements $F_{\mu \nu }\left(\mathbf{t}\right)$ of the QFIM $F\left(\mathbf{t}\right)$ are given by:

$$F_{\mu \nu }={\displaystyle \frac{1}{2}}\sum _{j,k=0}^{\mathcal{N}-1}\sum _{l=0}^3{\displaystyle \frac{Tr\left[S^{{\mathrm{T}}^{-1}}{M}_l^{jk}{S}^{-1}\frac{\partial \mathbb{V}}{\partial t_{\mu }}\right]Tr\left[S^{{\mathrm{T}}^{-1}}{M}_l^{jk}{S}^{-1}\frac{\partial \mathbb{V}}{\partial t_{\nu }}\right]}{{\lambda }_j{\lambda }_k-{(-1)}^k}}+2{\displaystyle \frac{\partial {\langle \widehat{\mathbf{d}}\rangle }^{\mathrm{T}}}{\partial t_{\mu }}}{\mathbb{V}}^{-1}{\displaystyle \frac{\partial \langle \widehat{\mathbf{d}}\rangle }{\partial t_{\nu }}}$$

(17)

where S and ${\lambda }_i$ are derived from the Williamson decomposition of the covariance matrix $\mathbb{V}$.

$$\mathbb{V}=SD{S}^{\mathrm{T}}$$

(18)

where $D=diag({\lambda }_0,{\lambda }_0,{\lambda }_1,{\lambda }_1,\dots )$, S belongs to the real symplectic group $Sp(2\mathcal{N},\mathbb{R})=\{S\mid S\Omega S^{\mathrm{T}}=\Omega \}$, and $\Omega ={\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]}^{\oplus \mathcal{N}}$ is the symplectic structure matrix. $M_l^{jk}$ is a block matrix composed of $2\times 2$ blocks, where all blocks are zero matrices except the block indexed by $jk$. The block labeled $jk$, denoted as $m_l^{jk}$, takes different forms depending on l: $m_0^{jk}={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]$, $m_1^{jk}={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]$, $m_2^{jk}={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$, $m_3^{jk}={\displaystyle \frac{1}{\sqrt{2}}}\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]$.

For high-dimensional covariance matrices, conventional methods cannot directly yield a valid Williamson decomposition and instead only provide constraints on the solution [18,20]. Here, we refer to the work of Martin Houdede et al. [21], which presents a robust and accurate numerical method to achieve this decomposition. The details are provided in Appendix C.

Therefore, as illustrated in the Figure 2, the method for analyzing the imaging process is as follows:

1.

Determine the amplitude point spread function $g(x,y)$ based on the actual imaging system.

2.

Discretize $g(x,y)$ using $\left\{{\varphi }_j\left(x\right)\right\}=\left\{R_b(x;x_j,\Delta )\right\}$ and $\left\{{\psi }_j\left(y\right)\right\}=\left\{R_b(y;y_j,\Delta )\right\}$.

3.

Obtain $X_2=\Lambda V^{†}$ via truncated singular value decomposition.

4.

Determine $X_1$ based on $\left\{u_j\left(x\right)\right\}$ and $\left\{{\varphi }_j\left(x\right)\right\}$.

The algorithm for calculating the QFIM $F\left({\mathbf{t}}_0\right)$ is:

1.

Input the parameter-dependent functions: ${\langle \widehat{\mathbf{d}}\rangle }_{\mathrm{in}}\left(\mathbf{t}\right)$, ${\mathbb{V}}_{\mathrm{in}}\left(\mathbf{t}\right)$, and the transformation matrix $X\left(\mathbf{t}\right)$, along with the parameter value ${\mathbf{t}}_0$.

2.

Calculate ${\langle \widehat{\mathbf{a}}\rangle }_{\mathrm{in}}\left(\mathbf{t}\right)$, ${\langle {\widehat{\mathbf{a}}}^{†}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{in}}\left(\mathbf{t}\right)$, and ${\langle \widehat{\mathbf{a}}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{in}}\left(\mathbf{t}\right)$ using the inverse of Equation (14).

3.

Compute ${\langle \widehat{\mathbf{a}}\rangle }_{\mathrm{out}}\left(\mathbf{t}\right)$, ${\langle {\widehat{\mathbf{a}}}^{†}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{out}}\left(\mathbf{t}\right)$, and ${\langle \widehat{\mathbf{a}}\otimes {\widehat{\mathbf{a}}}^{\mathrm{T}}\rangle }_{\mathrm{out}}\left(\mathbf{t}\right)$ via Equation (13).

4.

Determine ${\langle \widehat{\mathbf{d}}\rangle }_{\mathrm{out}}\left(\mathbf{t}\right)$ and ${\mathbb{V}}_{\mathrm{out}}\left(\mathbf{t}\right)$ using Equation (14).

5.

Compute the partial derivatives ${\displaystyle \frac{\partial {\langle \widehat{\mathbf{d}}\rangle }_{\mathrm{out}}\left(\mathbf{t}\right)}{\partial t_j}}$ and ${\displaystyle \frac{\partial {\mathbb{V}}_{\mathrm{out}}\left(\mathbf{t}\right)}{\partial t_j}}$.

6.

Assign the value ${\mathbf{t}}_0$ to the parameter $\mathbf{t}$ and perform the Williamson decomposition on ${\mathbb{V}}_{\mathrm{out}}\left({\mathbf{t}}_0\right)$ to obtain $S^{-1}$, $S^{-\mathrm{T}}$, and the symplectic eigenvalues $\left\{{\lambda }_j\right\}$ using the method described in Appendix C.

7.

Calculate and return the QFIM $F\left(\mathbf{t}\right)$ based on Equation (17).

Based on this framework, it is possible to automatically and rapidly compute the QFIM corresponding to any finite-dimensional parameter vector. Furthermore, by increasing the numerical precision and modeling accuracy, any desired finite computational precision can be achieved. In terms of describing interactions, compared to methods such as Gaussian channels [19] or symplectic transformations [18], this approach is not compatible with nonlinear processes like squeezing transformations. However, it offers advantages including convenient analysis, the ability to quickly construct composite maps via matrix multiplication and direct sums, efficient use of computational memory, and adaptability to cases where the numbers of input and output modes are unequal. Regarding the computation of the QFIM, this method integrates two efficient algorithms to achieve fast calculation of the matrix for multi-mode Gaussian states.

3. Application

3.1. Two-Point Resolution

Giacomo Sorelli et al. [15] discussed a two-dimensional imaging problem—the resolution of two point sources. The scenario involves two point-like emitters in the object plane, each in a Gaussian state. After passing through an imaging system described by a Gaussian point spread function, the light reaches the image plane, where the task is to estimate the separation distance between the two points using the imaging apparatus. Let the amplitude point spread function from object plane coordinates $(x,x^{\prime })$ to image plane coordinates $(y,y^{\prime })$ be

$$G(x,y,x^{\prime },y^{\prime })=\sqrt{\kappa }{\displaystyle \frac{1}{\sigma \sqrt{\pi }}}{e}^{-{\displaystyle \frac{{(x-y)}^2+{(x^{\prime }-y^{\prime })}^2}{2{\sigma }^2}}}$$

(19)

The object plane contains two point-like squeezed light sources located at $({\displaystyle \frac{s}{2}},0)$ and $(-{\displaystyle \frac{s}{2}},0)$, respectively. The objective is to estimate their separation distance s. Since only the two points on the x-axis are considered, the problem can be decoupled as $G(x,y,x^{\prime },y^{\prime })=\sqrt{\kappa }g(x,y)g(x^{\prime },y^{\prime })$, where

$$g(x,y)={\displaystyle \frac{1}{\sqrt{\sigma \sqrt{\pi }}}}{e}^{-{\displaystyle \frac{{(x-y)}^2}{2{\sigma }^2}}}$$

(20)

is the normalized amplitude point spread function. With reference to the method of Giacomo Sorelli et al. [15], we calculated the QFI for the separation distance s between two point sources $\delta \left(x-{\displaystyle \frac{s}{2}}\right)$ and $\delta \left(x+{\displaystyle \frac{s}{2}}\right)$ (Note: A point source is described by a Dirac $\delta$ function, which is not square-integrable. Its handling is detailed in Appendix D). The calculation was performed for two thermal states, each with a mean photon number $\overline{n}=100$.

$$\langle \widehat{\mathbf{d}}\rangle =\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right],\mathbb{V}=(2\overline{n}+1)\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(21)

coherent states $\overline{n}=100$

$$\langle \widehat{\mathbf{d}}\rangle =\sqrt{2\overline{n}}\left[\begin{array}{c}1\\ 0\\ -1\\ 0\end{array}\right],\mathbb{V}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(22)

squeezed vacuum states $\overline{n}={\overline{n}}_s={sin}^{2}r=3$

$$\langle \widehat{\mathbf{d}}\rangle =\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right],\mathbb{V}=\left[\begin{array}{cccc}{e}^{-2r}& 0& 0& 0\\ 0& e^{2r}& 0& 0\\ 0& 0& e^{2r}& 0\\ 0& 0& 0& e^{-2r}\end{array}\right]$$

(23)

squeezed coherent states $\overline{n}={\overline{n}}_c+{\overline{n}}_s$, where ${\overline{n}}_c=100$, ${\overline{n}}_s={sin}^{2}r=3$

$$\langle \widehat{\mathbf{d}}\rangle =\sqrt{2{\overline{n}}_c}\left[\begin{array}{c}1\\ 0\\ -1\\ 0\end{array}\right],\mathbb{V}=\left[\begin{array}{cccc}{e}^{-2r}& 0& 0& 0\\ 0& e^{2r}& 0& 0\\ 0& 0& e^{-2r}& 0\\ 0& 0& 0& e^{2r}\end{array}\right]$$

(24)

For the numerical calculations, we set the parameters as follows: The amplitude point spread function is given by $\sqrt{\kappa }g(x,y)=\sqrt{\kappa }{\displaystyle \frac{1}{\sqrt{\sigma \sqrt{\pi }}}}{e}^{-{\displaystyle \frac{{(x-y)}^2}{2{\sigma }^2}}}$ with $\sigma =1$ and the transmission coefficient $\kappa =0.35$. Correspondingly, the basis function width is set to $\Delta =0.01$. The basis functions at the object plane cover a range of width $L_{\mathrm{obj}}=5$, and those at the image plane cover a range of width $L_{\mathrm{obj}}=L_{\mathrm{obj}}+6\sigma =11$. After performing the singular value decomposition, the largest singular value is $\sqrt{{\eta }_0}=0.997$. The smallest singular value retained is the 20th one, with $\sqrt{{\eta }_{19}}=1.08\times {10}^{-10}$. In the computation process, the partial derivatives in Equation (17) are approximated using the central difference formula, ${\displaystyle \frac{df\left(x\right)}{dx}}\approx {\displaystyle \frac{f(x+\delta )-f(x-\delta )}{2\delta }}$, with a step size of $\delta =2\Delta =0.02$.

Under this parameter configuration, the QFI for the two-point separation ${\displaystyle \frac{s}{\sigma }}$ is calculated and plotted as a connected curve, with the results shown in the upper panel of Figure 3. For each calculated data point, the QFI per effective photon ${\displaystyle \frac{F{\sigma }^2}{\eta \overline{n}}}$, obtained using the method from reference [15] (denoted as $F_1$) and that from our proposed method (denoted as $F_2$) are compared. The relative deviation ${\displaystyle \frac{F_1-F_2}{F_1}}$ is illustrated in the lower panel of Figure 3.

Under this parameter configuration, the discrepancy between the two methods is $\left|{\displaystyle \frac{F_1-F_2}{F_1}}\right|<3\%$, thereby validating the correctness of our approach using existing results. For the two-point resolution problem, the estimation precision for a coherent state exhibits different curves depending on the relative phase between the two light sources. For estimating small separations ${\displaystyle \frac{s}{\sigma }}$, the coherent state achieves optimal precision when the relative phase is $\pi$, and it reaches a higher precision limit compared to the case when ${\displaystyle \frac{s}{\sigma }}\gg 1$, demonstrating a super-resolution phenomenon. This result is close to the upper bound of QFI given by Cosmo Lupo et al. [5]. For the squeezed vacuum state, it shows no advantage over the coherent state in the region of small${\displaystyle \frac{s}{\sigma }}$, which is of primary interest. Building upon the work in reference [15], we consider coupling the optimal coherent state with a squeezed vacuum state to form a squeezed coherent state. By suppressing its displacement fluctuations, the precision can be further enhanced, surpassing the shot noise level.

3.2. Spatial Multi-Parameter Estimation of Optical Fields

To extend imaging from the two-point problem to broader quantum imaging applications, we now investigate the estimation of the unknown spatial distribution of an optical field. An optical field with a specific profile, after passing through an imaging system, forms a distorted spot on the image plane. Our task is to estimate the spatial distribution of the optical field at the object plane based on the field measured at the image plane.

Consider an optical field with an unknown spatial profile, represented by $\widehat{E}\left(x\right)={\widehat{\mathbf{a}}}_0^{\left(u\right)}{u}_0\left(x\right)$. The function $u_0\left(x\right)$ is expanded using a set of basis functions $\left\{h_j\left(x\right)\right\}$

$$\widehat{E}\left(x\right)={\widehat{\mathbf{a}}}_0^{\left(u\right)}\left[c_0{h}_0\left(x\right)+c_1{h}_1\left(x\right)+\dots \right]$$

(25)

where the expansion coefficients $[c_0,c_1,\dots ]$ are our estimation targets. Together with the known basis functions $\left\{h_j\left(x\right)\right\}$, they determine the characteristics of $u_0\left(x\right)$. We consider three schemes for mode expansion using $\left\{h_j\left(x\right)\right\}$ and compare them. The first is the Hermite-Gaussian mode, commonly used in quantum optics and laser physics: $h_j\left(x\right)={\displaystyle \frac{1}{\sqrt{2^{j}j!4\sqrt{2}\sqrt{\pi }}}}{e}^{-{\displaystyle \frac{x^2}{2\times {\left(4\sqrt{2}\right)}^2}}}{H}_j\left({\displaystyle \frac{x}{4\sqrt{2}}}\right)$. The second is a typical basis function on a closed interval, the Legendre polynomials: $h_j\left(x\right)=\sqrt{{\displaystyle \frac{2j+1}{L_{\mathrm{obj}}}}}{L}_j\left({\displaystyle \frac{2x}{L_{\mathrm{obj}}}}\right)$. The third is obtained from the singular value decomposition of the amplitude point spread function $\sqrt{\kappa }g(x,y)=\sqrt{\kappa }{\displaystyle \frac{1}{\sqrt{\sigma \sqrt{\pi }}}}{e}^{-{\displaystyle \frac{{(x-y)}^2}{2{\sigma }^2}}}$, giving the right singular vectors: $h_j\left(x\right)={\sum }_{k=0}^{k_{\max }}{V}_{kj}^{†}{R}_b(x;x_k,\Delta )$.

We choose a representative and intuitive hypothetical field distribution for $u_0\left(x\right)$: $E_{\mathrm{test}}\left(x\right)=C_{\mathrm{test}}\left(x+{\displaystyle \frac{L_{\mathrm{obj}}}{2}}\right){\left(1-{\displaystyle \frac{4x^2}{L_{\mathrm{obj}}^2}}\right)}^2$, where $C_{\mathrm{test}}$ is a normalization constant. Its profile is shown in Figure 4 (a). We consider reconstructing $E_{\mathrm{test}}\left(x\right)$ using the first 6 modes of the three aforementioned schemes. All three can achieve reconstruction with high accuracy; therefore, neglecting the influence of higher-order modes, we calculate the QFIM per average photon number ${\displaystyle \frac{F}{\overline{n}}}$ for the parameter vector $[c_0,c_1,c_2,c_3,c_4,c_5]$.

Our parameter settings are as follows: the basis function width is $\Delta =0.01$, $L_{\mathrm{obj}}=25$, and correspondingly $L_{\mathrm{img}}=L_{\mathrm{obj}}+6\sigma =49$. For the amplitude point spread function, $\sigma =4$ and the transmission coefficient is $\kappa =0.074$ (chosen to avoid singular values greater than 1, which would cause model failure). After singular value decomposition, the largest singular value is $\sqrt{{\eta }_0}=0.947$. The smallest singular value retained is the 15th one, with $\sqrt{{\eta }_{14}}=6.70\times {10}^{-6}$. The step size for the central difference approximation is $\delta =2\Delta =0.02$.

Under the three basis function schemes, the QFIM per average photon number ${\displaystyle \frac{F}{\overline{n}}}$, yielded by the coherent state is shown in Figure 4, respectively.

Furthermore, we calculated the diagonal elements of the QCRB $F^{-1}$ for the three basis function schemes and compared them with the results ${\left.F_{jj}^{-1}\right|}_{\mathrm{SVD}}$ obtained from the right singular vectors. This comparison reveals that the QCRB ${\left.F_{jj}^{-1}\right|}_{\mathrm{SVD}}$ provided by the right singular vector basis expansion is optimal for every parameter’s precision lower bound. This conclusion serves as a verification, using quantum estimation theory, of a classical optics finding [22]. Additionally, after testing with various field profiles, we found that these computational results are independent of the specific values of the parameter vector $[c_0,c_1,c_2,c_3,c_4,c_5]$ and are general in nature.

Next, we consider further enhancing the estimation precision for the expansion coefficients within the framework of the right singular vector basis. We note that a bright squeezed coherent state $|\alpha ,r\rangle$, where the displacement direction aligns with the squeezing direction, can effectively improve estimation accuracy. The average photon number of a squeezed coherent state is given by $\overline{n}={\left|\alpha \right|}^2+{sinh}^{2}r$. Considering practical experimental constraints, the coherent-part mean photon number ${\overline{n}}_c={\left|\alpha \right|}^2$ is relatively easy to prepare. However, the squeezed-part mean photon number ${\overline{n}}_s={sinh}^{2}r$ is typically limited to a smaller mean photon number due to technical challenges. Here, we focus on the scenario where $\alpha \gg r$.

Unlike the classical case, when the bright squeezed coherent state $|\alpha ,r\rangle$ is used as the probe state,

$$\langle \widehat{\mathbf{d}}\rangle =\left[\begin{array}{c}\sqrt{2}\alpha \\ 0\end{array}\right],\mathbb{V}=\left[\begin{array}{cc}{e}^{-2r}& 0\\ 0& e^{2r}\end{array}\right]$$

(26)

the calculated QFIM becomes non-diagonal. This is because, as a non-classical optical field, the squeezed state exhibits quantum correlations between photons in different modes $\left\{h_j\left(x\right)\right\}$. However, the diagonal elements of $F^{-1}$ are significantly reduced, and the matrix remains approximately diagonal. Its off-diagonal elements are of the same order of magnitude or one order smaller than the smallest diagonal element, $F_{00}^{-1}$, yet are numerically significantly smaller than $F_{00}^{-1}$. Thus, overall, the squeezed coherent state provides an improvement in precision.

Under the premise that $\overline{n}\gg 1$, we fixed the squeezed-part mean photon number ${\overline{n}}_s$ to a fixed value and incrementally increased the total photon number $\overline{n}$ to calculate the QFIM. The observed trend in its values demonstrates an improvement characterized by a specific coefficient relative to the standard quantum limit. Keeping the total photon number fixed at $\overline{n}={10}^4$ while varying ${\overline{n}}_s$, we found that the degree of improvement offered by the squeezed coherent state is jointly influenced by the values of the parameter vector $[c_0,c_1,c_2,c_3,c_4,c_5]$ and the singular values $\sqrt{{\eta }_j}$. We investigated the enhancement effect provided by the squeezed-part mean photon number under two different optical field energy distributions.

From Figure 5, it can be observed that the performance gain offered by the squeezed state over the coherent state is more pronounced for larger singular values $\sqrt{{\eta }_j}$, and larger coefficients $c_j$. However, as the mean photon number in the squeezed-part ${\overline{n}}_s$ increases, the marginal improvement in estimation precision achieved by further squeezing gradually diminishes.

In general, for coherent imaging, one must account for not only the intensity distribution but also the phase distribution of the optical field, both of which are crucial. Therefore, we consider a more general complex-valued scenario represented by the parameter set $[c_0{e}^{i{\theta }_0},c_1{e}^{i{\theta }_1},c_2{e}^{i{\theta }_2},c_3{e}^{i{\theta }_3},c_4{e}^{i{\theta }_4},c_5{e}^{i{\theta }_5}]$. Using the same imaging system and parameters, with the right singular vectors as the basis set for calculation, we obtained the following result:

For the coherent state, the QCRB $F^{-1}$ remains diagonal. This indicates that when using classical resources, there exists no statistical correlation among the parameters $[c_0,c_1,c_2,c_3,c_4,c_5,{\theta }_0,{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5]$.

For the squeezed coherent state, as shown in Figure 6, $F^{-1}$ is a block-diagonal matrix. Quantum entanglement induces statistical correlations among the estimates for different amplitude distribution parameters $[c_0,c_1,c_2,c_3,c_4,c_5]$ and among the estimates for different phase parameters $[{\theta }_0,{\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5]$. However, no statistical correlation exists between the amplitude part and the phase part. Therefore, estimating the amplitude distribution and the phase distribution independently is optimal. However, due to the presence of correlations within each part, performing separate measurements on individual modes cannot extract complete quantum information. Consequently, the optimal measurement scheme that approaches the QCRB requires a specific form of joint measurement.

Furthermore, we also examined the achievable estimation precision for the phase distribution of the optical field. Limited by the number-phase uncertainty relation, the estimation precision for the amplitude distribution and the phase distribution cannot be enhanced simultaneously. Employing a phase-squeezed state can likewise improve the estimation of the phase distribution. The degree of this improvement is influenced by the singular values $\sqrt{{\eta }_j}$ and the energy distribution of the optical field, exhibiting a trend similar to that observed for amplitude estimation.

Overall, for estimating the spatial distribution of $u_0\left(x\right)$ in the unknown optical field $\widehat{E}\left(x\right)={\widehat{\mathbf{a}}}_0^{\left(u\right)}{u}_0\left(x\right)$, we transform the problem into a multi-parameter estimation task via the basis function expansion method. Among various expansion schemes, we confirm that using the right singular vectors as the basis functions achieves the optimal estimation precision.

Compared to the coherent state, the squeezed coherent state can further enhance the estimation precision for either the amplitude or the phase of the basis function expansion coefficients. The degree of this improvement increases with larger singular values $\sqrt{{\eta }_j}$, larger expansion coefficients $c_j$, and a higher mean photon number in the squeezed-part ${\overline{n}}_s$.

Due to the quantum correlations among different modes, when utilizing a squeezed coherent state as the resource, estimating the amplitude or phase distribution parameters requires a joint measurement scheme to approach the QCRB.

4. Conclusions

This paper establishes an imaging theoretical model based on Gaussian states and a corresponding method for calculating the quantum Fisher information of an imaging system, investigating quantum imaging problems utilizing squeezed states. First, we investigated the two-point imaging problem. Consistent with existing literature, our study confirms that a squeezed vacuum state does not enhance coherent imaging quality. We further explored two-point imaging with a squeezed coherent state, and the results demonstrate its superior performance compared to both thermal and coherent states.

Subsequently, we extended the research to the imaging of general patterns, specifically addressing the multi-parameter estimation problem for the spatial distribution features of an image. The squeezed coherent state exhibits the potential to surpass the shot noise limit in multi-parameter estimation, with the optical field’s energy distribution and system losses identified as key factors influencing the degree of squeezing-enhanced performance. However, due to potential quantum correlations among different parameters, separate measurements targeting individual modes are insufficient for extracting complete quantum information. Therefore, realizing the quantum Cramér-Rao bound by devising optimal experimental measurement schemes, and further proceeding to a comprehensive evaluation of the quantum enhancement effect, remain important open questions for future exploration.

The quantum imaging model and computational methodology presented here can be extended to more complex scenarios, such as imaging schemes based on multimode squeezed states, as well as to incoherent imaging systems [16]. Furthermore, it paves the way for exploring applications in practical imaging contexts, including Raman microscopy imaging [10,11,12] and stimulated Brillouin scattering imaging [13].