Squeezed Coherent States in Double Optical Resonance

: In this work, we consider a “ Λ -type” three-level system where the ﬁrst transition is driven by a radiation ﬁeld initially prepared in a squeezed coherent state, while the second one by a weak probe ﬁeld. If the squeezed ﬁeld is sufﬁciently strong to cause Stark splitting of the states it connects, such a splitting can be monitored through the population of the probe state, a scheme also known as “double optical resonance”. Our results deviate from the well-studied case of coherent driving indicating that the splitting proﬁle shows great sensitivity to the value of the squeezing parameter, as well as its phase difference from the complex displacement parameter. The theory is cast in terms of the resolvent operator where both the atom and the radiation ﬁeld are treated quantum mechanically, while the effects of squeezing are obtained by appropriate averaging over the photon number distribution of the squeezed coherent state.


Introduction
One of the most intriguing accomplishments of quantum optics over the last 35 years is the ability to generate quantum states of light that do not have any classical analog and tailor their statistical properties. Among these states of non-classical radiation, there is a class of states that have received considerable attention both from a theoretical, as well as from an experimental viewpoint due to their numerous applications [1] in modern quantum technologies, i.e., the "squeezed" states of radiation [2]. A quantum state is referred to as "squeezed" if the variance of one of the quadrature amplitudes of the state is modified in such a way that it becomes smaller than the respective variance of a vacuum or a coherent state. Such states are referred in the bibliography under the terms "Squeezed Vacuum State" (SQVS) and "Squeezed Coherent State" (SQCS), respectively. Due to Heisenberg's uncertainty principle, squeezing always results in the increase of the conjugated quadrature variance above the variance of a vacuum, in the case of the former, and the variance of a coherent state, in the case of the latter [3].
The first experimental observation of squeezed light was reported back in 1985 in the pioneering work of Slusher et al. [4], who achieved four-wave-mixing in an atomic vapor of sodium atoms. Since then, significant advances in the generation and detection of squeezed light [5] gave the green light for the implementation of squeezing as a resource in a series of different applications in basic research and technology. Among these applications, the ones that stand out are related to the reduction of quantum noise in optical communications [6], the detection of sub-shot-noise phase shifts [7,8], the ability to achieve maximum sensitivity in interferometry at lower laser powers [9], which was implemented in the detection scheme of gravitational waves [10], the storage of quantum memory [11,12], necessary for quantum information tasks, and most recently, the ability of quantum-enhanced microscopy and effective bio-imaging without the danger of damaging the cell sample [13,14].
The underlying physical mechanism behind the last application is thoroughly intriguing since it is based on the effective yield enhancement that one can achieve by inducing non-linear processes with squeezed radiation. As has been known since the 1960s [15,16], any non-linear light-matter interaction depends on the quantum statistical properties of the radiation that are embodied in its correlation functions. This realization naturally led to a series of studies [17][18][19][20][21][22][23][24][25][26][27][28][29][30] on non-linear phenomena induced by fields that may exhibit photon bunching or superbunching properties, resulting in values of the intensity correlation functions higher than those of a coherent state. The simplest example illustrating this dependence is the transition from a bound state to a continuum via the absorption of N photons, i.e., N-photon ionization. The derivation of the transition probability per unit time for N-photon ionization when all real bound atomic intermediate states are assumed sufficiently far from resonance [31][32][33][34] indicates that the rate of the process is proportional to an effective N-photon matrix element multiplied by the Nth-order intensity correlation function [15,16]. In view of the dependence of the correlation functions on the stochastic properties of the radiation, the rate of such a process can be affected dramatically by the intensity fluctuations of the source [35]. As a prototype example, one can take an 11-photon ionization process induced by thermal radiation whose Nth order intensity correlation function is given by N! × I N , where I is the average intensity, i.e., N! times larger than the respective correlation function of the coherent field. For N = 11, the N! factor arising from the strong intensity fluctuations of the chaotic field gives us an enormous enhancement of around seven orders of magnitude, which in fact has been observed experimentally in the past for a process of such a high order [36,37]. The situation becomes even more interesting if one considers squeezed radiation, which under certain values of the parameters may exhibit superbunching properties resulting in correlation functions even larger than those of the chaotic field [38]. As an example, the Nth order correlation function of a bright (high intensity) SQVS is (2N − 1)!! times larger than the respective correlation function of the coherent field [39], an enhancement factor even larger than the N! of the chaotic field. This enhancement factor has been measured experimentally in a recent work by Spasibko et al. [40] by generating optical harmonics of the order of 2-4 from a bright SQVS.
The pronounced interplay between the non-linear character of an atomic transition [41] and the photon correlation effects of the radiation that induces it has also been investigated in the context of the dynamics of an atomic transition between two bound states that are connected via a strong field with various quantum photon statistics. In this scenario, the dynamics can be monitored either through the calculation of the spontaneous emission spectrum of the upper state (fluorescence) or through the calculation of the population of a third state weakly connected to one of the two bound states, serving as a probe. The second scheme is often found in the bibliography under the term "Double Optical Resonance" (DOR) [42,43]. If the driving field connecting the two bound states is sufficiently strong, then each of the two states is split into doublets, whose energy separation depends on the coupling strength with the field. This effect is widely known as Stark or Autler-Townes splitting [44]. The case of strong driving by a chaotic field of arbitrary bandwidth has been investigated in great detail [45][46][47][48] and predicts a variety of interesting discrepancies from the well-studied case of coherent driving, which to the best of our knowledge still remain to be explored experimentally. The problems of resonance fluorescence into a squeezed vacuum reservoir [49][50][51][52][53][54][55][56][57], as well as the weak-field absorption of a two-level atom by various quantum states of radiation in the narrow bandwidth limit [58] have also been investigated in depth.
In a recent work [59], we studied the single-and two-photon Autler-Townes splitting profile in realistic atomic transitions driven by bunched and superbunched radiation and obtained a counterintuitive behavior resulting from the complex interplay between the form of the photon probability distributions of the squeezed vacuum and chaotic fields and the compounded non-linearity of the system, arising both from the strong driving and the order of the process. Our results indicated that much has to be still discovered concerning the strong interaction between atoms and non-classical radiation. The main hindrance in the experimental investigation of such systems so far has been the difficulty of producing stable, high intensity squeezed sources with controllable stochastic properties in the lab. However, in view of recent works [60] providing novel methods to overcome experimental difficulties associated with the generation of squeezed coherent radiation, the study of the strong driving of a bound-bound transition by a radiation field initially prepared in an SQCS is timely. Since squeezed light is not amenable to simulation in terms of classical stochastic processes, we adopt a fully quantum mechanical treatment in terms of the resolvent operator, involving averaging over the photon number distribution of the SQCS, a method valid in the "zero bandwidth" limit, i.e., in the limit where the bandwidth of the source is sufficiently smaller than the natural decay of the excited state.

Materials and Methods
We begin by considering an atom initially resting in its ground state |g in the presence of a quantized radiation field in a two-mode Fock state with n photons in the first mode with frequency ω 1 and m photons in the second mode with frequency ω 2 . The atom can absorb one photon from the first mode and move to the excited state |a , while the latter is coupled to another excited state |b via the emission of one of the m photons. The Hamiltonian of the system consists of three parts, namely: the atomic Hamiltonian H A = ∑ jh ω j |j j| (j = g, a, b), the Hamiltonian of the radiation field H R = ∑ 2 i=1h ω i a † i a i , and the interaction Hamiltonian under the rotating wave approximation, − , wherē hω g ,hω a , andhω b are the energies of the atomic states |g , |a , and |b , respectively, g 1 and g 2 represent the coupling strengths between those states in units of frequency, while a i and a † i are photon annihilation and creation operators, respectively. The atomic operators σ The eigenstates of the unperturbed Hamiltonian H 0 ≡ H A + H R of the compound system "atom + radiation" are |I = |g |n 1 |m 2 , |A = |a |n − 1 1 |m 2 , and |B = |b |n − 1 1 |m + 1 2 , with energieshω I =hω g + nhω 1 + mhω 2 ,hω A =hω a + (n − 1)hω 1 + mhω 2 , andhω B =hω b + (n − 1)hω 1 + (m + 1)hω 2 , respectively. The detunings from resonance of the two transitions are defined as To account for the spontaneous decay of the excited state |a , we make the substitution ω A →ω A = ω A − i 2 γ, in view of which, our system now becomes open. Note that this method introduces a decay rate without accounting for the repopulation of the ground state. It is however a good approximation as long as the first transition is sufficiently strong so that the Rabi frequency Ω 1 of the first transition is sufficiently larger than the decay rate γ. In such a case, the interaction between the uncoupled states of the compound system |I and |A causes a splitting into doublets energetically separated byhΩ 1 (Stark splitting). Such a splitting can be monitored through the calculation of the probe state population (state |b ) as a function of ∆ 2 . Our ultimate goal is the study of the splitting profile in the case where the |g ↔ |a transition is driven by a strong field prepared in a squeezed coherent state. The switch over from the initial Fock state to a squeezed coherent state can be realized through appropriate averaging of the probe state population over the photon statistics distribution of the squeezed coherent field, as will be described in detail below. A schematic presentation of the system we study is depicted in Figure 1.
Our problem is formulated in terms of the resolvent operator, which is the Laplace transform of the time evolution operator U(t) = exp(−iHt). Taking the Laplace transform with s being the usual Laplace variable and making the change of variables s = −iz, we obtain the relation G(z) = (z − H) −1 , where G(z) is the resolvent operator [29,61,62]. Using the relation H = H 0 + V, it is straightforward to show that G(z) obeys the equation z−H 0 is the unperturbed resolvent operator. In view of this equation, the matrix elements of the resolvent operator in the compound system basis satisfy the equations: Field Field Δ 1 Figure 1. Schematic presentation of the system of study. A strong radiation field prepared initially in a squeezed coherent state drives the |g ↔ |a transition, resulting in a Stark splitting of both states. This splitting is monitored through the calculation of the population of state |b as a function of ∆ 2 , which is assumed to be weakly coupled to state |a through a probe field.
Solving for G BI , one obtains: The matrix elements of the time evolution operator in the compound system basis, i.e., U ij (t), i, j = I, A, B, are related to the respective matrix elements of the resolvent operator through the inverse transform [29,61,62]: where x + = x + iη, with η → 0 + . In order to calculate the integral of Equation (5), one should first calculate the roots of the third order polynomial appearing in the denominator of Equation (4). If we denote these three roots by z 1 , z 2 , and z 3 , the resulting expression is: The population of the probe state |b at times t > 0 is given by: The population of |b depends non-linearly on the photon numbers n and m through the expressions of the compound system energieshω I ,hω A , andhω B , as well as the matrix elements V AI and V BA that reflect the Rabi frequencies of the |g ↔ |a and |a ↔ |b transitions, respectively, via the relations 2V AI =hΩ 1 =hg 1 √ n and 2V BA =hΩ 2 = hg 2 √ m + 1.
We focus on the dependence of P b (t) on n and t, by adopting the notation P b (n, t). The dependence on m is not of particular importance as long as the probe field coupling strength g 2 is chosen to be sufficiently smaller than all the other rates that appear in the problem at hand. To capture the case where the field driving the |g ↔ |a is initially prepared in a state other than Fock, one can average the desired quantity (probe state population in our problem) over the corresponding photon number distribution of the initial field state. This is a standard method in quantum optics [63]. However, one should note that it is strictly valid in the "zero bandwidth" approximation, i.e., in the limit where the field bandwidth is much smaller than the natural decay γ. We are particularly interested in the case where the field driving the first transition is initially prepared in an SQCS.
An SQCS is defined as the state resulting upon acting on the vacuum state with the squeezing operatorŜ(ζ) followed by the displacement operatorD(a), i.e.: |SQCS ≡D(a)Ŝ(ζ) |0 (8) The squeezing and displacement operators are given by the relationŝ the complex squeezing parameter and a = |a|e iφ is the complex displacement parameter. The photon number distribution of an SQCS is given by [64]: where H n (x) is the n th -order Hermite polynomial. By substituting a = |a|e iφ in Equation (9) and adopting the definition Φ ≡ φ − θ/2, one finally gets: In Figure 2a, we plot the photon number distribution of an SQCS with a squeezing parameter r = 0.5 and compare it to the distribution of a coherent field (r = 0). As becomes evident, depending on the phase difference Φ, the distribution deviates from its Poissonian form, by exhibiting super-Poissonian (red line) or sub-Poissonian statistics (teal line). The situation becomes even more interesting for larger values of r. In Figure 2b, we plot the photon number distribution of an SQCS with a squeezing parameter r = 2 for two values of Φ and observe an obvious non-classical behavior with vivid oscillations in both cases. These oscillation were interpreted in the past by Schleich and Wheeler [65,66] as resulting from the interference of error contours in phase space. In the Φ = 0 case (red line), the distribution is peaked at n = 0, and the oscillations are much faster compared to the Φ = π/2 case (teal line), where the distribution is peaked at a slightly higher photon number than the peak of the Poissonian distribution and the oscillations appear only for photon numbers larger than the position of the peak. It should be noted that the photon number distributions of Figure 2 acquire discrete values for each value of the integer n, but are depicted as continuous since the photon number scale over which they are plotted is large. To avoid any misconception, there is another form of Equation (10) often found in the bibliography where the distribution is essentially the same, but it exhibits sub-Poissonian statistics for Φ = 0 and super-Poissonian statistics for Φ = π/2. This is a reflection of whether the SQCS is defined as |SQCS ≡D(a)Ŝ(ζ) |0 or |SQCS ≡Ŝ(ζ)D(a) |0 , which are not equal since the operatorsD(a) andŜ(ζ) do not commute. In any case, one can go from one definition to the other through the relation: Returning back to our problem, we average the population of the probe state |b over the photon number distribution of Equation (10), according to: In what follows, we will be concerned about the behavior of P SQCS b (t) as a function of ∆ 2 for different parameters of the SQCS, or put otherwise, the effects of various schemes of squeezing on the resulting Stark splitting profile imprinted on the population of the probe state.

Results and Discussion
Among other parameters, Equation (12) depends also on time. Therefore, in order to study the behavior of the probe state population as a function of ∆ 2 , one should first make a choice of the interaction time t = T during which the driven system is exposed to the radiation. In Figure 3a, we plot P SQCS b (t) as a function of time for various detunings ∆ 2 and r = 0. Under this choice of parameters, we notice that the population of the probe state reaches its steady-state value at about t = 25γ −1 . Numerical investigations of the behavior of P SQCS b (t) as a function of time for different combinations of the parameters ∆ 2 , |a| 2 , r, and Φ revealed that the choice of t = 30γ −1 ≡ T always guarantees that the system is well within its steady-state regime. Therefore, we adopt this choice of time for the calculations throughout the paper. In Figure 3b, we plot P SQCS b (T) as a function of ∆ 2 for various |a| 2 and r = 0 parameters corresponding to a coherent state. Note that, as expected, for r = 0, the distribution does not depend on the choice of Φ. We confirm that as |a| 2 increases and the mean Rabi frequencyΩ 1 = g 1 |a| becomes larger than the spontaneous decay rate γ, the single peak structure (black line) splits into two peaks forming the well-known Autler-Townes doublet structure [44]. The distance between the two peaks is equal to the mean Rabi frequencyΩ 1 , which is proportional to the square root of the mean photon number of the coherent staten = |a| 2 . As long asΩ 1 is smaller than γ, the increase of |a| 2 will result in power broadening of the profile [67] until it splits into the doublet.
T=30γ -1 Before continuing to the case of strong driving by an SQCS, we should note that the exact choice of the Rabi frequency of the probe transition is not of direct relevance to our problem, as long as it is sufficiently weaker than all of the other rates that appear in the derivation. If this condition is satisfied, then the only effect of changing the Rabi frequency of the probe will be a linear change in the total population of the probe state (y-axis), but not in the ratio between individual peaks or widths. Since our aim is the study of the form of the splitting profile, the exact values of the probe state population are of no importance. For the same reason, we are not concerned about the possible natural decay of the probe state |b to some other unobserved state.
In Figure 4a, we plot the population of the probe state as a function of ∆ 2 for various degrees of squeezing and Φ = 0. The splitting profile shows great sensitivity to the value of the squeezing parameter r, becoming significantly broader as the latter increases. The physics behind this broadening of the profile is related to the strong amplitude fluctuations of the field that are translated into fluctuations in the value of the Rabi frequency of the transition that cause partial smearing of the doublet structure after taking the average over the SQCS distribution. However, as seen in Figure 4b, under the same choice of parameters, but with Φ = π/2, the profile only slightly deviates from the profile exhibited under coherent state driving (r = 0). An effective way of interpreting this sensitivity on the phase difference Φ is through the form of the SQCS photon number variance. In particular, the photon number variance of the SQCS is given by the relation [68]: (∆n) 2 = |a| 2 e 2r cos 2 Φ + e −2r sin 2 Φ + 2(sinh r cosh r) 2 (13) For Φ = 0, Equation (13) reduces to (∆n) 2 = |a| 2 e 2r + 2(sinh r cosh r) 2 . In this case, it is straightforward to show that the increase of the squeezing parameter will always lead to the increase of the variance of the SQCS (inset of Figure 4a) and therefore to the broadening of the total profile after averaging over the SQCS distribution. However, the picture is drastically different if Φ = π/2, where Equation (13) takes the form (∆n) 2 = |a| 2 e −2r + 2(sinh r cosh r) 2 . It is easy to check that in this case, the variance is not increasing monotonically as a function of r, but it acquires a minimum value at a position that depends on the choice of |a| 2 . In the considered case |a| 2 = 60, the minimum is positioned at r = 0.91. The variance remains smaller than its value at r = 0 (corresponds to a coherent state) up to the position r = 1.53 and becomes larger than that thereafter (inset of Figure 4b). Furthermore, the value of the variance in the vicinity of the minimum is not much smaller than its r = 0 value. This behavior explains why for the squeezing parameters considered in Figure 4b, the SQCS splitting profile does not exhibit significant deviations from the corresponding coherent profile. It also explains why the peaks of the teal and orange lines corresponding to the values r = 1 and r = 1.5, respectively, appear larger than the peaks of the black line (coherent field), due to the sub-Poissonian form of the SQCS photon number distribution.
As r increases, for Φ = 0, the peaks of the splitting profile appear at slightly smaller detunings than the coherent splitting, contrary to the Φ = π/2 case where the peaks tend towards higher detunings. This can be attributed to the behavior of the mode of the SQCS distribution (most probable value) as a function of r in each case. In particular, for Φ = 0, the mode of the distribution tends towards the zero photon number, while for Φ = π/2, it is positioned at photon numbers slightly higher than |a| 2 , depending on the value of r. This tendency is also evident in Figure 2. We should note that the peaks of the resulting profile are not separated by the Rabi frequency corresponding to the mode of the photon number distribution, nor the Rabi frequency corresponding to its mean photon number, given by the relationn = |a| 2 + (sinh r) 2 , for every Φ. As argued in previous work, the exact shape of the resulting profile stems from a complex interplay between the properties of the probe state population function associated with the order of the process and the structure of the photon probability distributions of the driving field [59]. The behavior of the mode as a function of r and different values of Φ can however give us good evidence for the expected behavior of the position of the peaks. This tendency of the peaks towards slightly higher (absolute) detunings for Φ = π/2 as r is increased, i.e., the increasing in the frequency distance between the two peaks, becomes even more evident by inspecting Figure 5a. In this figure, we examine the resulting profile for very large values of the squeezing parameter corresponding to well beyond the state-of-the-art squeezing. At this point, we should mention that squeezing is most usually measured in the bibliography in terms of the squeezing factor in dB units. The connection between the squeezing factor R and the squeezing parameter r used in our work is given by the relation: Therefore, the very high squeezing factor of 20dB reported back in 2016 [69] corresponds approximately to a squeezing parameter r = 2.3. As becomes evident in Figure 5a, under extreme squeezing, the resulting profile is distorted in a rather unusual way. The peaks of the profile are positioned at even higher detunings from those of Figure 4b, and the total profile is broadened mainly towards one direction. This effect can be interpreted through the particular form of the SQCS photon number distribution for Φ = π/2. As seen in Figure 2b, the distribution exhibits a large peak at a photon number slightly higher than |a| 2 , followed by an oscillatory behavior for larger photon numbers. As r increases, the distribution maintains its form qualitatively, but the probability of the mode (most probable value) of the distribution tends towards larger values, while the ratio of the probabilities between subsequent peaks of the distribution is decreased. This indicates that as r is increased, the weight of the distribution is transferred from its mode towards the peaks of higher photon numbers, leading also to the increase of the variance of the distribution. Therefore, the splitting profile after averaging over the SQCS distribution exhibits smaller peaks, but it is more broadened towards higher (absolute) detunings.
Another interesting scenario occurs for |a| 2 = 0, but a non-zero squeezing parameter. This choice of parameters corresponds to what is widely known as a squeezed vacuum state (SQVS). The SQVS is defined as |SQVS ≡Ŝ(ζ) |0 , and its photon number distribution is given by: wheren = (sinh r) 2 is the mean photon number of the distribution. As seen from Equation (15), the SQVS distribution exhibits oscillatory behavior with zero values for odd photon numbers. The mode of the distribution is sharply peaked at n = 0, and its variance is given by the relation (∆n) 2 = 2n(n + 1). In Figure 5b, we plot the population of the probe state when the strong field is initially prepared in an SQVS. In agreement with previous work [59], we notice that the positions of the peaks are rather insensitive to the increasing of the squeezing parameter (translated to intensity in [59]), while the widths of the peaks are significantly increased, as one would expect by inspecting the SQVS variance for increasing r. From a physical viewpoint, this sharp increase of the width is due to the superbunching effect inherent in squeezed vacuum sources, translated to strong intensity fluctuations that smear out the profile. Note that the exact shape of the profile also depends on the choice of the coupling strength g 1 , which is however chosen to be the same through all of our calculations (g 1 = 0.7γ), providing us the ability for straightforward comparison between the profiles of different figures. Lastly, in Figure 6a,b, we consider the case where the squeezed field is detuned from resonance with the |g ↔ |a transition. Similar to the case of coherent driving, the profile exhibits an asymmetry with the larger peak appearing at positive detunings for positive ∆ 1 , and vice versa. This asymmetry can be interpreted in terms of the dressed states of the strong field transition [42]. The resulting doublets are expressed as linear combinations of the ground and the first excited state of the atom, with equal weights in the case of exact resonance (∆ 1 = 0). However, as ∆ 1 is increased, one can show that the upper dressed state of the doublet contains more of the state |a and less of the ground state, while the opposite is true for the lower dressed state. In this case, the probe state is connected through a larger dipole moment to the upper state of the doublet than the lower, resulting in an asymmetry in the splitting profile. The asymmetry is examined for r = 1.5 and various values of the phase difference Φ. As Φ is increased, the width of the two peaks is decreased, and the peaks move towards slightly higher detunings (their frequency distance is increased). The width of each peak also depends on ∆ 1 . Increasing ∆ 1 results in the decrease of the width of the large peak and the increase of the width of the smaller one. Finally, we report that the asymmetry ratio (peak-to-peak) does not remain constant, but it increases with the increase of Φ. Probe state population  Figure 6. (a) Population of the probe state |b as a function of ∆ 2 for various angles Φ. The squeezed field is detuned from resonance with the |g ↔ |a transition by ∆ 1 = γ. (b) Population of the probe state |b as a function of ∆ 2 for various angles Φ. The squeezed field is detuned from resonance with the |g ↔ |a transition by ∆ 1 = 3γ. In both panels, |a| 2 = 60, T = 30γ −1 , g 1 = 0.7γ, and r = 1.5. Black lines: Φ = 0, teal lines: Φ = π/6, orange lines: Φ = π/4, and purple lines: Φ = π/2.

Concluding Remarks
In this work, we have examined the Autler-Townes splitting profile resulting from the strong driving of two bound atomic states by an SQCS. The splitting was probed through a third state assumed to be weakly coupled to the upper bound state in a three-level "Λ-type" system. Since squeezed radiation is not amenable to modeling in terms of classical stochastic processes, we adopted a fully quantum mechanical treatment in terms of the resolvent operator and averaged the population of the probe state over the photon number distribution of the SQCS in the zero-bandwidth limit, i.e., in the limit where the bandwidth of the squeezed field is sufficiently smaller than the natural decay of the excited state. Our results indicate that the resulting splitting profile is greatly affected by the parameters of the SQCS, i.e., the squeezing parameter r and its phase difference Φ from the complex displacement parameter. In particular, we showed that for Φ = 0, the profile becomes significantly broader with the increase of r, while under the same values of parameters, but with Φ = π/2, increasing r up to two results in a profile resembling the one acquired with coherent driving. This degree of squeezing is well within today's squeezing capabilities, while the intensities necessary for the observation of single-photon Stark splitting in atomic systems are of the order of 1 W/cm 2 [59]. The case of extreme squeezing was also investigated, revealing unusual distortions of the splitting profile, as well as the cases of SQVS driving and off-resonant strong driving resulting in an asymmetric profile. To the best of our knowledge, the investigation of the strong driving of a bound-bound atomic transition by an SQCS remains still an open and interesting experimental problem. At the same time, recent progress in quantum engineering using ion traps opens up the possibility to efficiently prepare squeezed states of ionic motion that do not posses any classical analog [70][71][72]. Although such systems are not directly analogous to our system of study, they should be highlighted before closing this work, due to their numerous applications in spectroscopy, in quantum metrology, and in quantum information technology [73].
Author Contributions: Both authors contributed in the conceptualization, methodology and preparation of the manuscript. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.