Establishing and achieving the ultimate limit of optical resolution is long-standing problems in physics with profound and widespread impacts on a number of disciplines, from engineering to medical practice. Continuous advances in physics have the disruptive capability of redefining such limit pushing the boundaries of image resolution.

New imaging systems are commonly benchmarked against the well-known Rayleigh resolution criterion. Due to the wavelike nature of light, the image of a point-source is a spot characterized by the point-spread function (PSF) of the optical system. The Rayleigh criterion establishes that two point-sources are difficult to resolve if their transverse separation is smaller than the width of the PSF, which in turn is determined by the Rayleigh length ${\mathit{x}}_{\mathit{R}}=\lambda D/R$, where $\lambda$ is the wavelength, $R$ the radius of the pupil of the optical system, and $D$ the distance to the object [1]. Notable methodologies that bypass the Rayleigh resolution criterion include fluorescence microscopy, which exploits controlled activation and deactivation of neighboring emitters [2], and quantum optics, which exploits the fact that states with exactly $N$ photons have an effective $N$-times smaller wavelength [3,4].

Inspired by quantum metrology, whose general goal is to exploit quantum physics to boost measurement devices and methodologies, a modern formulation of the Rayleigh criterion has been recently proposed [5] inspired by the seminal work of Tsang et al. [6]. The key observation is that the classic Rayleigh resolution criterion applies only to the standard imaging approach where the intensity of the focused light is measured pixel-by-pixel on the image plane. However, optics offers much more than direct intensity detection, and a structured measurement may be able to extract the information carried by the phase of the light field. In a typical structured measurement, light is first preprocessed, e.g., in a multiport interferometer, and then measured [7,8]. This accounts to sorting the light field into a suitably defined set of optical modes, followed by mode-wise photodetection. Spatial-mode demultiplexing (SPADE) is a particularly efficient way of measuring the focused optical field. It has been shown to be optimal for a number of problems [5,6,9–14], including the textbook problem of resolving two neighboring pointlike sources. In SPADE, the field is demultiplexed in its transverse degrees of freedom. Oftentimes the transverse field is decomposed into the Hermite–Gaussian (HG) modes. Experimentally, SPADE can be implemented using a number of methodologies, including spatial light modulators [15–20], image inversion [21], photonic lantern [22], and multiplane light conversion [23–26].

Theoretical analyses show that SPADE is in principle the overall optimal measurement for a number of estimation and discrimination problems, with a notable improvement in performance with respect to direct detection (DD) in the single-photon regime and for bright sources [18,20,23–25]. Practical implementations are affected by cross talk [23,25,27–29], which represents the most important limitation of this methodology.

The natural field of application of SPADE as a super-resolution detection methodology is within astronomical imaging. Proposed use cases include the estimation of the angular separation between stars [6], or between a star and its exoplanet [30], and the detection of the presence of an exoplanet in the vicinity of a star [11,14]. All these applications refer to spatial degrees of freedom with the goal of extracting information about the spatial distribution of the source intensity in the transverse plane. Going beyond spatial degrees of freedom, SPADE may also be exploited to enhance spectroscopy. In a recent work, some of us have proposed to use SPADE to boost exoplanet spectroscopy [31], with a potential use to detect biosignature on distant planets. In this work, we develop on this idea and present, for the first time to the best of our knowledge, an experimental demonstration of super-resolved spectroscopy. The key enabler is the correlation between spatial and frequency degrees of freedom, e.g., when the spectrum of the exoplanet is slightly different from the spectrum of the star due to absorption lines. SPADE is useful, as it helps in sorting which photons are coming from the exoplanet and distinguishing them from those coming directly from the star. These photons are difficult to distinguish using direct detection, since the star is typically much brighter than the planet, and the two may have small angular separation with respect to the PSF of the optical system.

This paper develops as follows. In Section 2, we introduce the physical model and our assumptions. In Section 3, we present an analysis of exoplanet spectroscopy developed in terms of Fisher information and quantum Fisher information. This allows us to show how SPADE allows super-resolution spectroscopy. Section 4 introduces our experimental setup and discusses how it can be used to simulate the observation of a planetary system. The experiment, including calibration and measurement outcomes, is discussed in Section 5. Conclusions are presented in Section 6.

Consider a pointlike source located in position $x$ in the object plane that emits incoherent, monochromatic light at wavelength $\lambda$. Through an optical system, a focused image is produced in the image plane. In the far-field and paraxial approximation, the field in the image plane is characterized by the PSF of the optical system, denoted as $T$ [32]. The width of the PSF is determined by the Rayleigh length, and its shape depends on the pupil of the optical system. As long as the source is much larger than the wavelength, it is sufficient to use a scalar theory for the light field. Therefore, the scalar field in position $y$ in the image field is proportional to $T(y/M-x)$, where $M$ is the magnification factor.

In linear optics, the classical scalar theory can be directly lifted into the quantum theory, and the PSF gives the rule of transformation of the field operators in the object plane into those in the image plane [33–35]. Here, we are interested in the single-photon regime, where at most one photon is detected in the image plane per detection window. It is therefore sufficient to consider a single photon emitted by the source and how this photon arrives on the image plane. A single-photon state emitted at point $x$ in the object plane is described by the state $${|1}_{x,\lambda }⟩=a_{x,\lambda }^{†}|0⟩,$$(1)where $\{a_{x,\lambda },a_{x,\lambda }^{†}\}$ are the canonical operators that annihilate and create a photon of wavelength $\lambda$ at position $x$ in the object plane. When such a state passes through a diffraction-limited optical system, characterized by the PSF $T$, the single-photon state in the image plane reads $$|T_{x,\lambda }⟩=\int \mathrm{d}yT(y-x)b_{y,\lambda }^{†}|0⟩,$$(2)where $\{b_{y,\lambda },b_{y,\lambda }^{†}\}$ are the canonical operators that annihilate and create a photon of wavelength $\lambda$ at position $y$ in the image plane. Here, to simplify the notation, we have put $M=1$.

From a single pointlike source, we now move to the case of an extended incoherent source. Consider that a photon of wavelength $\lambda$ is emitted from position $x$ with probability $p(x,\lambda )$. Then, the state of the single-photon in the image plane is represented by the density matrix $$\rho =\sum _{x,\lambda }p(x,\lambda )|T_{x,\lambda }⟩⟨T_{x,\lambda }|.$$(3)

Our goal is to model the observation of light from an exoplanetary system, where the light coming from the star (primary source) is also scattered by the exoplanet (secondary source). Due to absorption through a planet’s atmosphere, the two sources may have different spectra. Modeling stars and planets as pointlike sources, the state of a single photon collected in the image plane is $$\rho =(1-ϵ)\sum _{\lambda }{f}_s(\lambda )|T_{x_s,\lambda }⟩⟨T_{x_s,\lambda }|+ϵ\sum _{\lambda }{f}_p(\lambda )|T_{x_p,\lambda }⟩⟨T_{x_p,\lambda }|.$$(4)Here, $ϵ$ and $1-ϵ$ are the relative intensities of planet and star, $|T_{x_{\star },\lambda }⟩$ represents the state of a single photon at wavelength $\lambda$, emitted either from the star ($\star =s$) or planet ($\star =p$), and $f_{\star }(\lambda )$ is the corresponding spectrum. The objective of exoplanet spectroscopy is to obtain information about the spectrum $f_p(\lambda )$ of the planet. This task comes with at least two challenges: (i) the star is much brighter than the planet, and therefore most of the photons detected come from the star; and (ii) if the transverse separation between star and planet is below the Rayleigh length, then the PSFs $T_{x_p,\lambda }$ and $T_{x_s,\lambda }$ will overlap substantially, making it difficult to distinguish which photons are coming from the planet.

To make our presentation more concrete, we assume a Gaussian PSF. Formally, this would follow from a pupil function with a Gaussian profile. However, a Gaussian PSF is also a good approximation for other PSFs with enough regularity [24] in the subdiffraction regime where the object is much smaller than the width of the PSF: this is exactly the regime we are interested in. A Gaussian PSF reads $$T(y-x)=\mathcal{N}{e}^{-\frac{{(y-x)}^2}{4{\sigma }^2}},$$(5)where the factor $\mathcal{N}$ is to normalize the PSF to 1, i.e., $\int \mathrm{d}y|T(y-x{)|}^2=1$. The width of the Gaussian profile can be identified with the Rayleigh length, $\sigma \equiv {\mathit{x}}_{\mathit{R}}$. It follows that in general the PSF depends on the wavelength through its width. Here, we assume a regime where the observed spectral range is narrow enough that the width of the PSF may be assumed constant.

Direct detection of the state in Eq. (4) amounts to address the field in each pixel in the image plane and to apply a spectral analyzer pixel by pixel. The outcome of such a measurement strategy is described by the probability of detecting a photon at position $y$ with wavelength $\lambda$, $$p_{\mathrm{DD}}(y,\lambda )=(1-ϵ)f_s(\lambda )|⟨1_{y,\lambda }|T_{x_s,\lambda }⟩{|}^2+ϵf_p(\lambda )|⟨1_{y,\lambda }|T_{x_p,\lambda }⟩{|}^2$$(6)$$=(1-ϵ)f_s(\lambda ){\mathcal{N}}^2{e}^{-{(y-x_s)}^2/2{\sigma }^2}+ϵf_p(\lambda ){\mathcal{N}}^2{e}^{-{(y-x_p)}^2/2{\sigma }^2}.$$(7)

We will compare direct detection with a methodology to obtain spectral measurements based on HG SPADE. In a structured measurement of the optical field, we first sort the transverse field in the image plane along basis elements ${\{|{\mathrm{\Psi }}_u⟩\}}_{u-\mathrm{0,1,2,}\dots }$. Then, we apply a spectral analyzer to each mode. The outcomes of this measurement are described by the probability of measuring a photon of wavelength $\lambda$ in mode $u$: $$p_{\mathrm{SPADE}}(u,\lambda )=(1-ϵ)f_s(\lambda ){|⟨{\mathrm{\Psi }}_u|T_{x_s,\lambda }⟩|}^2+ϵf_p(\lambda ){|⟨{\mathrm{\Psi }}_u|T_{x_p,\lambda }⟩|}^2.$$(8)HG SPADE is obtained when the basis elements are HG modes, i.e., $$|{\mathrm{\Psi }}_u⟩=\int \mathrm{d}y{\mathrm{\Psi }}_u(y)b_{y,\lambda }^{†}|0⟩,$$(9)with $${\mathrm{\Psi }}_u(y)={\mathrm{HG}}_u(y)=\frac{\mathcal{N}}{\sqrt{2^{u}u!}}{e}^{-\frac{y^2}{4{\sigma }^2}}{H}_u\left(\frac{y}{\sqrt{2}\sigma }\right),$$(10)and $H_u$ denotes the Hermite polynomials. Note that the width $\sigma$ is chosen to match that of the Gaussian PSF.

To study the problem of spectral analysis with subdiffraction resolution and compare different detection strategy, we may use the Fisher information as a figure of merit. Spectral analysis can be understood as a problem of multiparameter estimation, whose goal is to estimate the spectrum of the exoplanet, $f_p(\lambda )$, for a range of discrete values of the wavelength ${\lambda }_1$, ${\lambda }_2$, $\dots$, ${\lambda }_N$.

For given $n$ photons detected, the uncertainty in the estimation of the spectrum due to statistical fluctuations is quantified by the covariance matrix $\mathrm{\Sigma }$, with matrix elements $${\mathrm{\Sigma }}_{ij}=\mathbb{E}[f_p({\lambda }_i)f_p({\lambda }_j)]-\mathbb{E}[f_p({\lambda }_i)]\mathbb{E}[f_p({\lambda }_j)],$$(11)where $\mathbb{E}$ denotes the expectation value. According to the Cramér–Rao bound, for a given measurement, the covariance matrix is lower-bounded by the Fisher information matrix $$\mathrm{\Sigma }\ge \frac{1}{n}{F}^{-1},$$(12)where $F^{-1}$ denotes the inverse matrix of the Fisher information matrix. For a given measurement described by the probability mass distribution $p(x)$, we have $$F_{ij}=\sum _{x}p(x)\frac{\partial \log p(x)}{\partial f_p({\lambda }_i)}\frac{\partial \log p(x)}{\partial f_p({\lambda }_j)}.$$(13)

In direct detection, we have $x\equiv (y,\lambda )$. From Eq. (7), we obtain $$\frac{\partial \log p_{\mathrm{DD}}(y,\lambda )}{\partial f_p({\lambda }_i)}={\delta }_{\lambda ,{\lambda }_i}ϵ{\mathcal{N}}^2{e}^{-{(y-x_p)}^2/2{\sigma }^2},$$(14)where ${\delta }_{\lambda ,{\lambda }_i}$ is the Kronecker delta function. Due to the presence of the ${\delta }_{\lambda ,{\lambda }_i}$, the Fisher information matrix is diagonal, with elements $$F_{jj}^{\mathrm{DD}}=\sum _y{p}_{\mathrm{DD}}(y,{\lambda }_j){\left(\frac{1}{p_{\mathrm{DD}}(y,{\lambda }_j)}\frac{\partial p_{\mathrm{DD}}(y,{\lambda }_j)}{\partial f_p({\lambda }_j)}\right)}^2.$$(15)For a Gaussian PSF, this reads $$F_{jj}^{\mathrm{DD}}=\sum _y{p}_{\mathrm{DD}}(y,{\lambda }_j){\left(\frac{ϵ{\mathcal{N}}^2{e}^{-{(y-x_p)}^2/2{\sigma }^2}}{p_{\mathrm{DD}}(y,{\lambda }_j)}\right)}^2$$(16)$$=\sum _y{p}_{\mathrm{DD}}(y,{\lambda }_j){\left(\frac{ϵe^{-{(y-x_p)}^2/2{\sigma }^2}}{(1-ϵ)f_s({\lambda }_j)e^{-{(y-x_s)}^2/2{\sigma }^2}+ϵf_p({\lambda }_j)e^{-{(y-x_p)}^2/2{\sigma }^2}}\right)}^2$$(17)$$=\sum _y\frac{{ϵ}^2{e}^{-{(y-x_p)}^2/{\sigma }^2}}{(1-ϵ)f_s({\lambda }_j)e^{-{(y-x_s)}^2/2{\sigma }^2}+ϵf_p({\lambda }_j)e^{-{(y-x_p)}^2/2{\sigma }^2}}.$$(18)

In the following, to make the notation lighter, we put $F^{\mathrm{DD}}({\lambda }_j):=F_{jj}^{\mathrm{DD}}$.

In Eq. (17), the Fisher information is written as the average of the squared signal-to-noise ratio (SNR). The maximum value of the SNR is when $y=x_p$, yielding $$F^{\mathrm{DD}}(\lambda )\le {\left(\frac{ϵ}{(1-ϵ)f_s(\lambda )e^{-{(x_p-x_s)}^2/2{\sigma }^2}+ϵf_p(\lambda )}\right)}^2,$$(19)which follows from the fact that the average is always smaller than the maximum. This lower bound shows that the Fisher information from direct detection is of order ${ϵ}^2$ unless the sources are well separated, i.e., $$e^{-{(x_p-x_s)}^2/2{\sigma }^2}\ll \frac{ϵ}{1-ϵ}\frac{f_p(\lambda )}{f_s(\lambda )}.$$(20)If the sources overlap substantially, then $F^{\mathrm{DD}}\sim {ϵ}^2$. Otherwise, if they are well separated, then $F^{\mathrm{DD}}\sim ϵ$. In fact, from Eq. (18), we obtain $$F^{\mathrm{DD}}(\lambda )\ge \frac{{ϵ}^2}{(1-ϵ)f_s(\lambda )e^{-{(x_p-x_s)}^2/2{\sigma }^2}+ϵf_p(\lambda )}\simeq \frac{ϵ}{f_p(\lambda )}.$$(21)Here, we used the fact that the sum is always larger than one of the addends.

The form of the state in Eq. (4) implies that the Fisher information matrix is diagonal for any measurement. In fact, the argument of Section 3.A applies to any measurement. For HG SPADE, we have [to make the notation lighter, we use $p(u,\lambda )$ instead of $p_{\mathrm{SPADE}}(u,\lambda )$] $$F^{\mathrm{HG}}(\lambda )=\sum _{u}p(u,\lambda ){\left(\frac{1}{p(u,\lambda )}\frac{\partial p(u,\lambda )}{\partial f_p(\lambda )}\right)}^2$$(22)$$=\sum _{u}p(u,\lambda ){\left(\frac{ϵ|⟨{\mathrm{\Psi }}_u|T_{x_p,\lambda }⟩{|}^2}{(1-ϵ)f_s(\lambda ){|⟨{\mathrm{\Psi }}_u|T_{x_s,\lambda }⟩|}^2+ϵf_p(\lambda ){|⟨{\mathrm{\Psi }}_u|T_{x_p,\lambda }⟩|}^2}\right)}^2.$$(23)

If the optical system is carefully aligned toward the star, i.e., the PSF matching the lower HG modes, most of the photons from the star will couple into the fundamental mode ${\mathrm{\Psi }}_0$. Also, in principle, the photons from the star would not couple to mode ${\mathrm{\Psi }}_1$, as the latter is an odd function with respect to the position of the star. Finally, higher modes may be neglected in first approximation, as they are highly suppressed in the sub-Rayleigh regime [36]. Overall, by considering only the lower HG modes, we have $$\begin{gathered} F^{\mathrm{HG}}(\lambda )\simeq p(0,\lambda ){\left(\frac{ϵ|⟨{\mathrm{\Psi }}_0|T_{x_p,\lambda }⟩{|}^2}{(1-ϵ)f_s(\lambda ){|⟨{\mathrm{\Psi }}_0|T_{x_s,\lambda }⟩|}^2+ϵf_p(\lambda ){|⟨{\mathrm{\Psi }}_0|T_{x_p,\lambda }⟩|}^2}\right)}^2 \\ +p(1,\lambda ){\left(\frac{ϵ|⟨{\mathrm{\Psi }}_1|T_{x_p,\lambda }⟩{|}^2}{(1-ϵ)f_s(\lambda ){|⟨{\mathrm{\Psi }}_1|T_{x_s,\lambda }⟩|}^2+ϵf_p(\lambda ){|⟨{\mathrm{\Psi }}_1|T_{x_p,\lambda }⟩|}^2}\right)}^2 \end{gathered}$$(24)$$\simeq \frac{{ϵ}^2|⟨{\mathrm{\Psi }}_0|T_{x_p,\lambda }⟩{|}^4}{(1-ϵ)f_s(\lambda ){|⟨{\mathrm{\Psi }}_0|T_{x_s,\lambda }⟩|}^2+ϵf_p(\lambda ){|⟨{\mathrm{\Psi }}_0|T_{x_p,\lambda }⟩|}^2}+\frac{ϵ|⟨{\mathrm{\Psi }}_1|T_{x_p,\lambda }⟩{|}^2}{f_p(\lambda )}$$(25)$$\ge \frac{ϵ|⟨{\mathrm{\Psi }}_1|T_{x_p,\lambda }⟩{|}^2}{f_p(\lambda )}.$$(26)

Note that the first term in Eq. (25), which comes from photon detection in the fundamental mode ${\mathrm{\Psi }}_0$, is proportional to ${ϵ}^2$, whereas the second term, which comes from detection in mode ${\mathrm{\Psi }}_1$, is proportional to $ϵ$. The latter dominates in the regime where the planet is much dimmer than the star.

Comparing Eq. (26) with Eq. (19), we see that HG SPADE offers in principle a quadratic improvement in the scaling of the Fisher information when the transverse separation between star and planet is below the Rayleigh length. In the limit of an ultraweak signal coming from the planet, with direct detection the Fisher information vanishes with ${ϵ}^2$, whereas the scaling is with $ϵ$ in the case of HG SPADE.

By inspection of Eqs. (24)–(26), we note that the $O(ϵ)$ scaling of the Fisher information crucially follows from the fact that there is a vanishing probability that a photon emitted from the star is detected in the mode ${\mathrm{\Psi }}_1$, i.e., $⟨{\mathrm{\Psi }}_1|T_{x_s,\lambda }⟩=0$ in Eq. (24). In practice, this condition cannot be matched exactly due to misalignment and cross talk [23,27–29]. In both cases, we will have a nonzero probability that a photon emitted from the star is detected in the mode ${\mathrm{\Psi }}_1$. This implies that the term $⟨{\mathrm{\Psi }}_1|T_{x_s,\lambda }⟩$ does not vanish, and the Fisher information will scale as ${ϵ}^2$.

In other words, the better performance of HG SPADE, when compared with DD, comes from the fact that, ideally, the HG SPADE detection probability $p(1,\lambda )$ in Eq. (8) does not depend on $f_s(\lambda )$. This implies that photon detection in mode $|{\mathrm{\Psi }}_1⟩$ gives information only about the planet spectrum without background noise from the star. This ideal condition is only approximately reached due to misalignment and cross talk. In our experimental setup, the two spectra are nearly mutually orthogonal (i.e., they are nearly nonoverlapping). In this case, when cross talk is negligible, we have that $f_s(\lambda )$ and $p_{\mathrm{SPADE}}(1,\lambda )$ are also mutually orthogonal, i.e., ${\sum }_{\lambda }{f}_s(\lambda )p(1,\lambda )=0$. Experimentally, we are not able to access $f_s(\lambda )$, but we can estimate $p(0,\lambda )$, which is proportional to $f_s(\lambda )$ up to a correction term of order $ϵ$. Therefore, it makes sense to consider the scalar product between the probability vectors $p(0,\lambda )$ and $p(1,\lambda )$. We define $$\mathcal{S}=\sum _{\lambda }p(0,\lambda )p(1,\lambda ).$$(27)For nonoverlapping spectra and negligible cross talk, this quantity is close to zero, $\mathcal{S}\simeq 0$; otherwise, it increases with increasing cross talk in principle up to $\mathcal{S}=1$ when cross talk is so intense that the signal is equally spread on both modes. In conclusion, $\mathcal{S}$ quantifies the distance of the experimental setup to the ideal case where HG SPADE allows for optimal spectral estimation.

In general, the Fisher information depends on the chosen measurement strategy. The global bound on the covariance matrix for all possible measurements allowed by the principles of quantum mechanics is given by the quantum Cramér–Rao bound $$\mathrm{\Sigma }\ge \frac{1}{n}{Q}^{-1},$$(28)where $Q$ denotes the quantum Fisher information matrix. Due to the form of the state in Eq. (4), the quantum Fisher information is necessarily diagonal. This follows from the fact that photonic states at different wavelengths are mutually orthogonal (this is an extension to the quantum domain of the argument in Sections 3.A). The elements of $Q$ corresponding to the estimation of $f_p(\lambda )$ can be computed from the results of Ref. [26] (which in turn are derived from Ref. [30]) after a simple change of variables. We obtain $$Q(\lambda )=\frac{ϵ}{1-ϵ}\frac{1-w_{\lambda }}{f_p(\lambda )},$$(29)where $w_{\lambda }=⟨T_{x_p,\lambda }|T_{x_s,\lambda }⟩$ [30]. This expression shows the same $ϵ/f_p(\lambda )$ scaling as Eq. (26), hence indicating the optimality of the linear scaling of the Fisher information with $ϵ$.

In a proof-of-principle experiment, we simulate the use of HG SPADE to achieve super-resolved spectroscopy of a planetary system. The experimental setup is shown in Fig. 1. We use a fiber-coupled light emitting diode (LED) at telecom wavelength (Item 1 in Fig. 1) to simulate two incoherent pointlike emitters. The LED is attenuated using a 50 dB fiber attenuator (Item 3 in Fig. 1) and, by means of a fiber beam splitter (FBS) (Item 4 in Fig. 1), it is split into two fiber-coupled beams: Beam S simulates the star and Beam P simulates the exoplanet. Beam S is filtered using a fiber-coupled filter (Item 5 in Fig. 1) to induce a spectral shape different from that of Beam P (described below), as shown in Fig. 2. Then, it is free-space collimated using a fiber collimator (Item S in Fig. 1) and a two-lens system (items L and L in Fig. 1) to obtain a beam waist $w_0\simeq 300\mathrm{\mu m}$. Beam S, through two steering mirrors (items M and M in Fig. 1), impinges on a cube beam splitter (BS) (Item 9 in Fig. 1) to be combined with Beam P and, after crossing a fixed-film polarizer (Item 7’ in Fig. 1), is coupled with the free-space input port of a 300 μm waist HG demultiplexer (PROTEUS-C model from Cailabs) (Item 10 in Fig. 1).

The other output of the FBS (Beam P) is collimated using a fiber collimator (Item P in Fig. 1) and a two-lens system (items L and L in Fig. 1) to obtain a beam waist $w_0\simeq 300\mathrm{\mu m}$. Then, it crosses a free-space spectral filter (Item 6 in Fig. 1) generating the spectral shape of the simulated exoplanet (shown in Fig. 2). Beam P crosses a polarizer mounted on a motorized rotation stage (Item 7 in Fig. 1) and two steering mirrors (Items M and M in Fig. 1), one of which is mounted on a micrometer translation stage (Item 8 in Fig. 1). Then, it is reflected from the cube BS (Item 9 in Fig. 1), where it is combined with Beam S and crosses the fixed polarizer (Item 7’ in Fig. 1) to overcome the dependency of detector efficiency on photon polarization when the stage is rotated. Finally, Beam P is coupled with the HG demultiplexer. In this way, by translating the translation stage (Item 8 in Fig. 1), it is possible to move Beam P then changing the separation $d_a$ between the two beams. Moreover, by rotating the rotation stage, Beam P intensity is tuned and the intensity ratio $ϵ=N_P/N_S$ changes, where $N_S$ and $N_P$ are the number of photons impinging on demultiplexer emitted by S and P sources, respectively. The position and the intensity of Beam S remain fixed during the experiment. Once the two beams enter the demultiplexer, they are decomposed in the lowest-order HG modes: ${\mathrm{HG}}_{00}$; ${\mathrm{HG}}_{10}$; ${\mathrm{HG}}_{01}$; ${\mathrm{HG}}_{11}$; ${\mathrm{HG}}_{20}$; and ${\mathrm{HG}}_{02}$. The modes are converted in ${\mathrm{HG}}_{00}$ mode again and coupled with six single-mode fibers. Finally, the modes ${\mathrm{HG}}_{00}$, ${\mathrm{HG}}_{01}$, and ${\mathrm{HG}}_{10}$ are coupled to He-cooled superconducting nanowires single-photon detectors (NSPDs) (Item 11 in Fig. 1) through single-mode fibers equipped with polarization paddles to rotate the photon polarization and maximize the NSPDs efficiency. The spectra are acquired by scanning a motorized tunable filter (MTF) (Item 2 in Fig. 1). When photons impinge on the nanowires, the detectors generate electric pulses that, if exceeding the setup threshold, are recorded and counted by a commercial time tagger (Item 12 in Fig. 1) in the set temporal window of 50 ms. The overall detectors quantum efficiency (system efficiency) is about 80%, the dark count is lower than 20 Hz, and the reset time is about 100 ns. When the beams completely overlap (the separation $d_a$ between simulated point like emitters reduces to zero), the overall S+P beam symmetry is circular and the power leaked into first-order modes (${\mathrm{HG}}_{01}$ and ${\mathrm{HG}}_{10}$) is minimum. Such a minimum value, named “cross talk,” is due to demultiplexer manufacturing imperfections.

In general, the cross talk between ${\mathrm{HG}}_{00}$ and ${\mathrm{HG}}_{\mathit{nm}}$ is defined by the ratio $P_{\mathit{nm}}/P_{00}$, where $P_{\mathit{nm}}$ is the power on the ${\mathrm{HG}}_{\mathit{nm}}$ output fiber when only ${\mathrm{HG}}_{00}$ is injected with a $P_{00}$ power in the input. The cross talk is the main limiting factor in these kinds of experiments and is due to some signal ending up in high-order modes, even if the incoming light is fully matched, in terms of waist and direction, with the fundamental mode (if the crosstalk were negligible only ${\mathrm{HG}}_{00}$ should be excited for a fully matched incoming radiation). The measured cross talk $\chi$ of the first order modes in our case is $\chi =P_{01}/P_{00}+P_{10}/P_{00}\simeq 0.0035$.

The aim of the experimental work is to understand how HG SPADE can be used to perform spectroscopy of a P source (representing the exoplanet) when a strong S source (representing the star) is at a separation shorter than Rayleigh distance. For this purpose, a motorized tunable filter (MTF) is inserted at the output of the LED (it is equivalent to put it at the input of the demultiplexer). It selects a given wavelength window ($\mathrm{FWHM}=1\mathrm{nm}$) and is scanned between 1520 and 1569 nm. In this way, the system is able to acquire the spectrum of ${\mathrm{HG}}_{00}$, ${\mathrm{HG}}_{10}$, and ${\mathrm{HG}}_{01}$.

In the experiment, we operate in the single-photon regime and align the optical system to the brighter source. This is in contrast to what is done in other works, where the alignment is on the centroid or median point [6,11,15,16,25,37]. In our case, as the S source is much brighter than the P source, the procedure essentially coincides with alignment with the centroid; see Ref. [26] for more details. We want to stress that, in the case of a real observation, the only way to proceed is by aligning the optical system to the centroid, as the median point is unknown and the sources may have an arbitrary shape.

To align the demultiplexer, we superimposed Beam P with Beam S (by setting the translation stage at 0 μm) and maximize the signal in the ${\mathrm{HG}}_{00}$ mode. We collected the counts $C_0(\lambda )$ and $C_1(\lambda )$ in fundamental (${\mathrm{HG}}_{00}$) and first-order modes (${\mathrm{HG}}_1\equiv {\mathrm{HG}}_{01}+{\mathrm{HG}}_{10}$) as a function of wavelength $\lambda$ and for different values of the (dimensionless) source separation $d_a=d/w_0$ and intensity ratio $ϵ=N_P/N_S$. The counts are then normalized, obtaining $$C_0^N(\lambda )=\frac{C_0(\lambda )}{\sqrt{{\sum }_{{\lambda }^{\prime }}{C}_0{({\lambda }^{\prime })}^2}},$$(30)$$C_1^N(\lambda )=\frac{C_1(\lambda )}{\sqrt{{\sum }_{{\lambda }^{\prime }}{C}_1{({\lambda }^{\prime })}^2}},$$(31)which can be interpreted as a pair of vectors with 50 elements (one for each wavelength). In the experiment, the total number of photons $N_S$ impinging on demultiplexer is fixed as $$N_S={\int }_{1520\mathrm{nm}}^{1569\mathrm{nm}}{N}_S(\lambda )\mathrm{d}\lambda \simeq \mathrm{148,000}$$(32)during the scanning time $t_s\simeq 2.5\mathrm{s}$ (50 ms time windows for each of the 50 wavelengths).

Figure 3 shows the counts $C_0^N(\lambda )$ (upper panel) and $C_1^N(\lambda )$ (lower panel) as functions of wavelength for five values of separation $d_a$. This figure shows that spectral information from the exoplanet is enhanced, in comparison with that from the star, when observing the mode ${\mathrm{HG}}_1$. The spectra $C_0^N(\lambda )$ are nearly independent on the separation, unlike $C_1^N(\lambda )$. This corresponds to the fact that the background from the star is mostly coupled into the fundamental mode ${\mathrm{HG}}_{00}$. In particular, the left parts of the spectra (where the emission of the exoplanet is centered) are strongly dependent on separation $d_a$ as highlighted in the inset where the portion of the spectrum between 1538 and 1547 nm is shown. The spectra appear quite rough because of filter minimum step (1 nm); in principle, however, it would be possible to obtain higher spectral resolution by using a finer step or a grating dispersing light on an array of single-photon detectors.

Similarly, Fig. 4 shows the counts $C_0^N(\lambda )$ (upper panel) and $C_1^N(\lambda )$ (lower panel) as functions of wavelength for five values of intensity ratio $ϵ$. As above, the spectra $C_0^N(\lambda )$ are nearly independent on $ϵ$, whereas $C_1^N(\lambda )$ are not. Also, the left parts of the spectra are strongly dependent on $ϵ$, as highlighted in the inset. Note that changes of the star component in ${\mathrm{HG}}_1$ spectra (Figs. 3 and 4) for different separation and intensity ratio are a mere effect of the normalization. Being the whole spectrum normalized, the ratio between star peak and planet peak changes.

The acquired single-photon spectroscopy data are noisy if compared with ordinary spectroscopy. Being the relative uncertainty in each point of the order of $\frac{1}{\sqrt{N}}$ (where $N$ is the collected photon in a given point), the spectra show large fluctuation when the number of collected photons is low. To quantify the distinguishability of the photons coming from the simulated exoplanet in terms of their spectral properties, we may consider the scalar product between the vectors $C_0^N(\lambda )$ and $C_1^N(\lambda )$: $${\mathcal{S}}_P=\sum _{\lambda }{C}_0^N(\lambda )C_1^N(\lambda ).$$(33)Up to a different choice of normalization, this quantity is equivalent to the scalar product introduced in Eq. (27). Since in our setup the spectra of planet and star are nearly orthogonal (they have nearly nonoverlapping support in the spectral domain), this scalar product is close to zero when the system is able to distinguish the photons coming from the exoplanet from those coming directly from the star. Otherwise, as discussed in Section 3, cross talk may increase the value of ${\mathcal{S}}_P$ up to the point where the photons become hardly distinguishable and ${\mathcal{S}}_P\simeq 1$.

Figures 5 and 6 show the quantity ${\mathcal{S}}_P(d_a,ϵ)$ as a function of $d_a$ at fixed $ϵ$ and ${\mathcal{S}}_P(d_a,ϵ)$ as a function of $ϵ$ at fixed $d_a$, with ${\mathcal{S}}_P(d_a,ϵ=0.021)$ and ${\mathcal{S}}_P(d_a=\mathrm{0.33,}ϵ)$, respectively. To compare with direct detection, we have computed the quantity ${\mathcal{S}}_P$ that would be obtained by implementing the spectral analysis through direct detection. To this goal, we have simulated two direct detection systems centered with an exoplanet and a star. Using the same number of photons of the HG SPADE experiment ($N_S\simeq \mathrm{148,000}$), the direct detectors aligned with the exoplanet (at coordinate $y_1=0,y_2=w_0{d}_a$ on image plane) and the star (at coordinate $y_1=0,y_2=0$ on image plane) acquired ${\mathrm{dd}}_p$ and ${\mathrm{dd}}_s$ photons, respectively. The two systems collect the photons over two spatial regions with radius $w_0$ centered to the exoplanet and star and spectrally disperse them with the same resolution used for HG SPADE $\simeq 1\mathrm{nm}$. We compute $${\mathrm{dd}}_p(\lambda )=I_1{f}_s(\lambda )N_A(1-ϵ)+I_2{f}_p(\lambda )N_Aϵ,$$(34)$${\mathrm{dd}}_s(\lambda )=I_3{f}_s(\lambda )N_A(1-ϵ)+I_4{f}_p(\lambda )N_Aϵ,$$(35)where $$I_1={\int }_{-w_0}^{w_0}\mathrm{d}{y}_1{\int }_{w_0{d}_a-w_0}^{w_0{d}_a+w_0}{G}_s(y_1,y_2)\mathrm{d}{y}_2,$$(36)$$I_2={\int }_{-w_0}^{w_0}\mathrm{d}{y}_1{\int }_{w_0{d}_a-w_0}^{w_0{d}_a+w_0}{G}_p(y_1,y_2)\mathrm{d}{y}_2,$$(37)$$I_3={\int }_{-w_0}^{w_0}\mathrm{d}{y}_1{\int }_{w_0}^{+w_0}{G}_s(y_1,y_2)\mathrm{d}{y}_2,$$(38)$$I_4={\int }_{-w_0}^{w_0}\mathrm{d}{y}_1{\int }_{w_0}^{+w_0}{G}_p(y_1,y_2)\mathrm{d}{y}_2,$$(39)and $$G_s(y_1,y_2)=\frac{1}{2\pi w_0^2}{e}^{-\frac{-y_1^2+y_2^2}{2w_0^2}},$$(40)$$G_p(y_1,y_2)=\frac{1}{2\pi w_0^2}{e}^{-\frac{-y_1^2+{(y_2-d_a{w}_0)}^2}{2w_0^2}}.$$(41)

In analogy with Eqs. (30) and (31), we define the normalized vectors $${\mathrm{DD}}_p(\lambda ):=\frac{{\mathrm{dd}}_p(\lambda )}{\sqrt{{\sum }_{{\lambda }^{\prime }}{\mathrm{dd}}_p{({\lambda }^{\prime })}^2}},$$(42)$${\mathrm{DD}}_s(\lambda ):=\frac{{\mathrm{dd}}_s(\lambda )}{\sqrt{{\sum }_{{\lambda }^{\prime }}{\mathrm{dd}}_s{({\lambda }^{\prime })}^2}}.$$(43)

Considering the same number of photons, the uncertainty is of the order of a few percent and, as shown in the Figs. 5 and 6, the scalar product ${\mathcal{S}}_P^{\mathrm{DD}}={\sum }_{{\lambda }_i}{\mathrm{DD}}_s({\lambda }_i){\mathrm{DD}}_p({\lambda }_i)$ for direct detection spectroscopy remains close to 1 in the region of interest and is much larger than scalar product ${\mathcal{S}}_P$ in the HG SPADE case. This corresponds to the fact that, in the sub-Rayleigh regime ($d_a<1$), the HG SPADE is more capable than direct detection to distinguish the photons coming from the secondary source, hence allowing better estimation of their spectral properties.

In this work, we demonstrate a scheme for spectroscopy with spatial super-resolution. The scheme is based on spatial-mode demultiplexing (SPADE) of the optical field in its transverse components. We make use of the demultiplexer PROTEUS-C model from Cailabs, which decomposes the transverse field along the HG components. HG SPADE has been widely studied in the past few years, following the seminal work of Tsang et al. [6], as a means to achieve sub-Rayleigh estimation and discrimination of quantum states. Potential applications have been proposed to enhance astronomic observations, including the observation of exoplanets [11,14,30]. SPADE-based super-resolution imaging has been implemented in several ways, including spatial light modulators [15–20], image inversion [21], photonic lantern [22], and multiplane light conversion [23–26]; so far, none of these have ever been applied to spectroscopy. Here, we demonstrate the first application of HG SPADE to achieve super-resolved spectroscopy.

One of the challenges addressed by exoplanet science is to determine the atmospheric makeup of exoplanets. In particular, one searches for biomarkers such as oxygen or methane, whose presence is witnessed by absorption lines in the visible and near-infrared spectrum [38,39]. Gaining information about the spectrum of the exoplanet is difficult because of the angular vicinity to the star. Most of the photons collected come directly from the star, and it is difficult to discriminate the relatively few photons that contain information about the planetary atmosphere. In our proof-of-principle experiment, we simulate the observation and spectral analysis of a system composed of a primary source (the star) and a secondary source (the planet).

In principle, HG SPADE allows us to completely decouple the photons coming from the secondary source. Due to symmetry, if the demultiplexer is aligned toward the primary source, then only the photons from the secondary source are collected into the excited mode ${\mathrm{HG}}_1$. In practice, the scheme is limited by cross talk due to experimental imperfections in the demultiplexer as well as residual misalignment. Here, we introduce a quantity, which can be directly measured experimentally, that quantifies the capability of the system to distinguish the photons emitted by the secondary source and estimate their spectral properties. The crucial parameters are the relative intensity of the planet compared with the star and their angular separation. We show that there exists a regime, where the relative intensity and angular separation are not too small, where HG SPADE is capable of extracting substantial information about the spectrum of the secondary source, in particular much more than direct detection. The investigated relative brightness $ϵ$ range is between 0.007 and 0.095 because of low SNR for lower $ϵ$. The low SNR is due to several reasons, some of which could be easily improved in the next experiments, and others can also be improved but with greater effort. In our case, we use an LED to simulate the sources, whose intensity, polarization, and spectral shape may slightly change over the observation time, which in turn affects the experiment but not the performance of the developed spectroscopic technique. In fact, if the same system were used to observe a star–planet system, which is more stable than LED, the SNR degradation due to it would not be present. Another cause of SNR degradation is cross talk, which in our case is about 0.0035. Even if this value may appear large, it must be noted that the spatial-mode demultiplexing based on multiplane light conversion is a young technology with ample room for improvement. It is important to emphasize that, to detect an exoplanet, e.g., ${10}^{-8}$ less intense than a star, is not necessary to have a cross talk of ${10}^{-8}$ since: (i) if the sources are stable and the number of collected photons is large, a cross talk subtraction can be performed with great accuracy, and so the planet signal in ${\mathrm{HG}}_1$ mode must be comparable with the fluctuation of $\chi N_S$ and not with $\chi N_S$; and (ii) the planet motion can be used to track the spectral shape change of the ${\mathrm{HG}}_{01}/{\mathrm{HG}}_{10}$ modes compared with the ${\mathrm{HG}}_{00}$ mode, allowing an even stronger common mode noise rejection. These results demonstrate the potential usefulness of SPADE in exoplanet spectroscopy and pave the way to experiments beyond proof of principle, toward in-field demonstration of super-resolution spectroscopy from spatial-mode demultiplexing.

Acknowledgment. We thank Cailabs, 38 Boulevard Albert 1er, 35200 Rennes, France.

Author Contributions.L. S. A. designed the experimental setup and realized the apparatus; C. L. developed the theory and mathematical modeling, assisted with the experimental design, and supervised the experiment; L. S. A. and F. S. performed data acquisition and analysis; and D. P. developed the instruments control software. All authors discussed the results and contributed to the final manuscript.