Proposal for a continuous wave laser with linewidth well below the standard quantum limit

Due to their high coherence, lasers are ubiquitous tools in science. We show that by engineering the coupling between the gain medium and the laser cavity as well as the laser cavity and the output port, it is possible to eliminate most of the noise due to photons entering as well as leaving the laser cavity. Hence, it is possible to reduce the laser linewidth by a factor equal to the number of photons in the laser cavity below the standard quantum limit. We design and theoretically analyze a superconducting circuit that uses Josephson junctions, capacitors and inductors to implement a microwave laser, including the low-noise couplers that allow the design to surpass the standard quantum limit. Our proposal relies on the elements of superconducting quantum information, and thus is an example of how quantum engineering techniques can inspire us to re-imagine the limits of conventional quantum systems.

In our main text, we used an incoherent drive to pump the transmon qubit from its ground state to the excited state. In this note we provide a specific design recipe for how to build such an incoherent pump.
Our goal is to build an effective three-level-atom (see Fig. 1a), where the transition between the ground state |g and the second excited state |f is coherently driven while the second excited state experiences a fast decay to the first excited state |e . If the decay process is sufficiently fast, as the population of the atom is driven to the state |f , it quickly relaxes to |e and we achieve population inversion on the two lasing levels |g and |e . However, the single-photon transition between |g and |f for a transmon qubit is forbidden by selection rule. Therefore, we propose coupling a SNAIL qubit to the transmon qubit to form a composite system (see Fig. 1b). The key feature of the SNAIL qubit is that it has third order nonlinearity that makes the |g → |f transition allowed.
The level structure of the two qubit system is shown in Fig. 1c. We use |g t , |e t and |f t to represent the ground, first, and second excited states of the transmon qubit. For the SNAIL qubit we use |0 s , |1 s and |2 s , etc., to represent the ground, first excited, second excited, etc. states.
Coupling of the SNAIL and transmon qubits results in the hybridization of their states. Consequently, third order non-linearity of the SNAIL qubit can be used to drive the |0 s , g t → |1 s , e t transition as the state |1 s , e t is hybridized with the state |2 s , g t . This type of transitions are labeled by blue arrows in Fig. 1c. If the SNAIL qubit is also coupled to an output port, such that the relaxation of the SNAIL qubit [see Fig. 1c, black dashed arrows] is fast compared to the pump process (and also the transmon-cavity coupling, see Fig. 1c red arrows), then the two qubit system can form an effective three-level atom, in which |0 s , g t plays the role of the ground state (|g in Fig. 1a), |1 s , e t the role of the second excited state (|f in Fig. 1a), and |0 s , e t the role of the first excited state (|e in Fig. 1(a)).
The Hamiltonian of a SNAIL qubit coupled to a transmon qubit is H T = ω tt †t + k tt †t †tt whereŝ (t) is the photon annihilation operator for the SNAIL (transmon) qubit. We have also truncated the nonlinear parts of the Hamiltonians of both qubits to the lowest non-trivial order and dropped the rapidly rotating terms likê s †ŝ †ŝ † andŝŝŝ in SNAIL qubit Hamiltonian. Finally, we make the assumption that the SNAIL and transmon qubits are strongly detuned as compared with the strength of the linear coupling, i.e., 2∆ = |ω s − ω t | g 2 , and hence the modes of the two qubits are only weakly hybridized. In the dressed basis, with respect to g 2 coupling, the third order nonlinearity of the bare SNAIL mode results in a third order nonlinear coupling between the dressed SNAIL and transmon modes H 3 =ŝ †ŝ t + h.c., whereŝ andt are the dressed SNAIL and transmon operators. Applying classical drive to the bare SNAIL mode induces the two-photon pump process (see blue arrows in Fig. 1) via H 3 .
To demonstrate that the proposed two-qubit system functions as a three level atom, we numerically integrate the the master equation, Eq. (3). As the higher levels of the composite systems are weakly populated, we truncated the Hilbert space of the SNAIL qubit to allow a maximum of 6 Table I for values of parameters used to construct this figure. that the coherent drive applied to the SNAIL qubit, together with the photon loss from the SNAIL qubit, induce an effectively population transfer from |0 s , g t to |0 s , e t . We observe population inversion after ∼ 3µs. At long times the system reaches a steady state with a significant population inversion (with roughly 90% occupancy of the state |0 s , e t and 10% of the state |0 s , g t ). The residual population of the state |0 s , g t is caused by the decay of the |0 s , e t via its hybridization with the state |1 s , g t . The ripples on the population curves are caused by the classical drive on the SNAIL qubit which causes the SNAIL mode to have a fast oscillating component at the frequency (ω s + ω t ). In Fig. 2b, we plot the population of the ground, first and second excited states of the transmon qubit after tracing over the SNAIL qubit degrees of freedom. From this plot we observe that the transmon qubit is effectively being pumped from the ground state |g t to the first excited state |e t . From the two plots in Fig. 2, we observe that the higher excited states of the composite systems (e.g., |0 s , f t and |1 s , g t ) have very little population, especially the second excited state of the transmon qubit |f t , which justifies the truncation of the composite system Hilbert space in our numerical calculation. In this note we derive the effective photon loss operator for the laser cavity induced by a linear inductive coupling between the laser cavity (LC resonator) and the transmission line. We extend this description to the ABOCC in Sec. Supplementary Note 4.
The quantization of the LC resonator and the transmission line is discussed in Ref. [1,2]. The canonical position and momentum of the LC resonator are the node superconducting phase ϕ c and charge Q c . Using these coordinates, the quadratic Hamiltonian of the LC resonator can be quantized, similar to the Harmonic oscillator, viâ where Z c is the characteristic impedance of the LC resonator, Z c = L c /C c , and the Hamiltonian of the LC resonator, in second quantized form, is H c = ω câ †â , where the frequency of the LC resonator is ω c = 1/ √ L c C c . The voltage on the LC resonator is V c =Φ c and the current flow in the LC resonator is I c =Q c . We can express the voltage and current operators using the raising and lowering operators viâ Here we consider a single-mode transmission line that couples to the LC resonator via a linear inductor. The generalized flux along the transmission line is and the charge density along the transmission line q(x) can be quantized. After the quantization of the transmission line fields, φ(x) and q(x) areq where v p = 1/ √ LC is the wave speed along the transmission line and the dispersion relation of the mode with momentum k is ω 2 k = v 2 p k 2 , Z tl = L/C is the characteristic impedance of the transmission line, l is the total length of the transmission line.
We further assume that the linear inductive coupling element (L c-tl ) that couples the LC resonator and the transmission line, connects to the x = 0 point of the transmission line. The coupling Hamiltonian is, where the dimensionless parametersφ c andφ tl (k) are defined as Here we only consider the LC resonator and the transmission line parts of the Josephson micromaser to understand the dissipation induced on LC resonator by the transmission line. The Hamiltonian of the system under rotating wave approximation is where κ k = φ 2 0 L c-tlφ cφtl (k). We further treat the transmission line as a vacuum bath, and apply the Born-Markov approximation to simplify the dynamics of the cavity field as in Ref. [3,Chap. 8]. The dynamics of the LC resonator mode can be described by the master equation where ρ is the density operator for the LC resonator, and the decay rate The phase noise of the laser cavity field causes the phase of the laser light to fluctuate, which gives a finite linewidth to the laser. Before proceeding with a detailed analysis, we begin by summarizing the key points. The master equation Eq. (15) can be thought of as an eigenvalue problem Eq. (19). The spectrum of eigenvalues has one zero eigenvalue λ 0 = 0, which corresponds to the steady state solution of the laser, and a number of negative eigenvalues λ i =0 < 0 which correspond to the decaying modes. As we show in this note, the spectral representation of the two-time correlation function G(t + τ, t) that describes the decay of coherence consists of a linear combination of decaying exponentials with the decay time set by these negative eigenvalues, see Eq. (22). To obtain the laser line shape we go to Fourier space. In the Fourier representation, G(ω) consists of a linear combination of Lorentzians with width set by the negative eigenvalues. Moreover, almost all of the weight (∼ 97.32 % for the ABOCC laser example in this note) is carried by the Lorentzian with the largest nonzero eigenvalue (that is narrowest of the Lorentzians). It is precisely this Lorentzian that forms the central peak of the laser line shape, while the remaining Lorentzians contribute to slightly broadening the "pedestal" at the base of the central peak, see Fig. 3.
For a conventional laser system phase noise results in the two-time correlation function decaying as [3] The power spectrum of the laser is given by the Fourier transform of the two-time correlation function G(t + τ, t), which is a Lorentzian with full-width at half minimum 2D, where D is the linewidth of the laser field.
For a system that couples to a Markovian bath, which can be described by the master equation whereL is the super-operator acting on the system density operator (ρ), the time-evolution of the density operator can be formally solved by where V (t, t 0 ) is a time-evolution super-operator that acts on the system degrees of freedom and ρ(t 0 ) is the initial state of the laser system [4]. Also notice that the two-time correlator can be obtained via where U is the time-evolution operator for the system and the bath, R is the density operator for both the system and the bath and Tr S , Tr B , Tr S+B are trace over system, bath, system and bath degrees of freedom, respectively. The term Tr B U (τ )aR(t)U † (τ ) can be thought of as a time-evolution of the "state" aR(t) by a time period τ , which is {V (τ, 0)[aρ(t)]} according to Eq. (16). For the laser system, we are only interested in the laser linewidth when the system is stable, i.e., Tr S a † (τ )aρ s , where ρ s is the steady state of the system, then the two-time correlation function can be written as Similar to the time-evolution of a closed quantum system, where we study the eigenstates of the Hamiltonian operator to understand the system dynamics, we can also find the eigen-spectrum of super-operatorL (so called damping basis) to study the time-evolution of the open quantum systems whose dynamics is described by the master equation Eq. (15) [5][6][7]. The right-eigenstates of the super-operatorL can be defined aŝ where λ i is the corresponding eigenvalue. Notice that for any physical systems, the steady states are invariant under the time evolution, which means they are the nullspace of the super-operators. These states should be valid quantum states, i.e., they should have unit trace. Further, because the master equation should preserve the trace of the density operator, all the eigenstates with nonzero eigenvalues should have zero trace. Suppose aρ s in Eq. (18) can be expanded using the the right eigenstates as where c's are the corresponding expansion coefficients and c 0û0 corresponding to the steady state solution, the twotime correlation function can be calculated as, where we take the fact that the density operator for steady state of a typical laser system is purely diagonal in the Fock basis, which kills the first term. Within all the right eigenstates of the super-operatorL which have significant contribution to the decay of the two-time correlation function (c i Tr S [a †û(i) ] term is large), the largest eigenvalue (L is non-positive) contribute to the long-time performance, i.e., it controls the linewidth of the laser system. To decompose the operator using the eigenstates of the super-operatorL, notice that the super-operator is, in general, not Hermitian, the left eigenstates are needed to define the orthogonality relation. Following Ref. [5], the dual operator of the super-operator is defined as for any operator A and state B. Then the left-eigenstates of the super-operatorL iš with orthonormal relation Tr(v (i)û(j) ) = δ i,j .
To numerically solve the eigenstates of the super-operatorL, we notice in typical laser systems, the super-operator contains terms in the form ofL whereP (i) andQ (i) operate on the system Hilbert space. After expanding the operatorsP (i) andQ (i) using a suitable basis of the system Hilbert space and acting on the system density operator, in the matrix form we get where ρ (i) ≡P (i) ρQ (i) . We can redefine the density operator as a vector, and hence the super-operator becomes a matrix, whose eigenvectors can be easily solved. To explicitly show the matrix representation of the master equation, we re-group the indices as and the super-operator becomesL The left and right eigenvectors of the matrix representation of the super-operatorL can be re-mapped back to the original index convention to get back the matrix representation in the system Hilbert space. The vector product of left and right eigenvectors are equivalent to the product definition of the operators, i.e., the trace of the two operator product. When we consider the laser system, especially we consider the cavity field, the Hilbert space is infinite dimensional. We can truncate the Hilbert space of the cavity field, but, in practice, the dimension is still very large and hence it is not practical to find all the eigenvectors of the super-operator. However, we notice that the right eigenvector with the lowest nonzero |λ (i) | is very close to aρ s . The other eigenstates, with larger |λ (i) |, decay faster than this eigenstate and do not contribute significantly to the linewidth of the laser.
In Fig. 3a, we compare the numerical calculated normalized two-time correlation function g (τ ) = G(τ )/ n using direct time-evolution of the master equations versus the function d 1 e λ (1) τ where d 1 is a coefficient that we fit and λ (1) is the largest nonzero eigenvalue of the super-operatorL. In both the conventional laser system and ABOCC laser system, the decay of the normalized two-time correlation function from the numerical time-evolution method (solid lines) match well with the super-operator eigenvalues (dashed lines). For the ABOCC laser system, we also notice that at the small time-scale, the decay of the two-time correlation function slightly differs from the single eigenvalue fit (see the inset of Fig. 3a). This is caused by the finite overlapping to the other eigenstates which a responsible for the rapid initial decay. We fit the numerical time-evolved g (1) (τ ) function using first 20 eigenvalues (decay rates) of the super-operator of the ABOCC laser (blue dashed line), which gives a good fit to the direct time-evolution curve. From the fitting, we extract the overlap constant for the first nonzero eigenvalue of the super-operator is 0.9732, which means the main contribution of the long-time decay is given by the eigenstates with the largest nonzero eigenvalue. In 51 Supplementary Figure 3. Comparison between the exact numerical time evolution versus the eigenspectrum method. In (a), we show the comparison between the numerical calculated two-time correlation function g (1) (τ ) versus the first eigenvalue of the laser systems. The numerical calculation is shown as solid lines (blue for conventional laser and orange for ABOCC laser), while the dashed lines are the exponential decay with the decay rate calculated by the eigen-spectrum of the super-operator. In the inset of (a), we zoom in on the ABOCC laser g (1) for short time period. The blue dashed line is the fitted decay of g (1) using the first 51 nonzero eigenvalues of the ABOCC laser system. In (b), we show the Fourier transform of the g (1) function. The fitted g (1) decay (orange line) matches with the numerical time evolution result (blue dots) well. The result using only the first nonzero eigenvalue is shown as the green line. Parameters: for conventional laser, Γp = 30, Γc = 0.1 and g ∼ 1.013; for ABOCC laser: m0 = 43 (Zc = 150 Ω), Γ1 = 1.0, Γp = 3.58 and g = 0.8747. Fig. 3b, we further examine the contribution to the linewidth in frequency domain. We Fourier transform the g (1) (τ ) data from both time-evolved method (blue dots) and the eigen-spectrum method (orange and green curves). We notice that the Fourier transform of fitted g (1) (τ ) using first 51 eigenvalues of the super-operator (the orange line) gives good agreement with the time-evolve method (blue dots). If we just extract the lineshape from the first eigenvalue ( Fig. 3b green curve), the lineshape is slightly different from the time-evolution result (see inset). The difference is due to the finite overlap with the eigenstates that have a larger decay constant. However, the faster dacaying states only contribute to the short-time dynamics (see Fig. 3a inset), which gives a slightly larger background on the Lorentzian lineshape and slightly increase the linewidth. However, the central peak of the spectrum is dominated by the first nonzero eigenvalue. In Fig. 3b, we observe that the Lorentzian based on the first nonzero eigenvalue (green curve) and the fitted time-evolution using the first 51 eigenvalues, and the numerically calculated lineshape all match well.
Compared to direct time-evolution method of calculating the laser dynamics, especially the steady state photon distribution and the laser linewidth or two-time correlation function, the eigen-spectrum method is significantly more computationally efficient. To further improve computational efficiency, we truncate the density matrices to include only the main diagonal and the first few minor diagonals. We verify that this truncation does not affect the part of the eigen-spectrum that we are interested in by increasing the number of minor diagonals until convergence to the exact solution is reached for all values of n ≤ 1000 for the conventional and SGBO laser and m 0 ≤ 1000 for the ABOCC laser. In our main text, the ABOCC coupling circuit between the cavity and the transmission line is shown in Fig. 1. The Hamiltonian for the ABOCC coupling is where E J:c-tl is the coupling Josephson junction energy, L c-tl is the coupling linear inductance, E J:c is the π-junction Josephson energy, and the node phases operators areφ c andφ tl as labeled in Fig. 4. The π-junctions are used to correct the dispersion given by the nonlinear coupling between the cavity and the transmission line. Here, we will ignore these two π junctions at the beginning of the discussion. Further, we define a dimensionless parameter r = With the second quantization of the transmission line field and the LC resonator mode (see Appendix Supplementary Note 2), the phase across the ABOCC circuit iŝ The coupling Hamiltonian (without the π junctions) becomes where we have used the fact that the resonator operators and the transmission line field operators commute. We use the Baker-Campbell-Hausdorff formula to transform the exponential of the operators We will use this result in order to reach a normal ordering in which the LC resonator operators are to the right of the transmission line operators. The expansion of Eq. (31) to third order in the transmission line field operators yields where h c and h tl are the dimensionless Hamiltonian acting solely on the cavity field and the transmission line, respectively. The cavity-transmission line coupling is expanded in orders of the transmission line field operators and the first, second and third order terms are labeled as h 1 , h 2 and h 3 .
where, C tl and C c are two constants. Notice that h sys and h tl , induced by the coupling circuit, contributes new nonlinearities to the cavity and the transmission line. These nonlinearities, especially the nonlinearity of the cavity field, will degrade the laser performance by shifting the cavity frequency as the number of photons in the cavity increases. To compensate for these dispersive effects on the cavity, we include a π junction as shown in Fig. 4, in which the Josephson energies satisfy This π junctions cancel out the non-linear contributions of h tl thus reducing dephasing. In the following discussion we focus on the remaining terms, which are the cavity-transmission line coupling terms h 1 , h 2 and h 3 . Expansion of Eq. (32) yields where the summation k is from −∞ to +∞, the summations with prime k,q and k,q,p omit the terms in which any of the summation indices are equal.
We further assume that the coupling strength between the LC resonator and the transmission line, which is controlled by E J:c-tl , is small compared to the LC resonator mode frequency and the transmission line dynamics, which allows us to apply the rotating-wave approximation (RWA). Applying the rotating wave approximation to h 1 of Eq. (36a), and restoring dimensions of energy, we obtain Defining the nonlinear operator for the cavity field and applying the Born-Markov approximation to trace over the transmission line degrees of freedom, the first order coupling contributes to a resonator dissipator term where the decay rate is Before we move to the higher order nonlinear terms in the ABOCC Hamiltonian, we calculate the constant parameter C TL , which is given by Eq. (34). In the quantum optics regime, we assume that the transmission line has a large length, where we can assume l → ∞. In this limit, we can approximate the summation of k by the integral of k as Further, in quantum optics regime, especially for the system that Born-Markov approximation applies, the cavity frequency is the dominant frequency to the coupling bandwidth θ and the system-bath coupling strength. Here we explicitly assume that the cutoff frequencies for the system-bath coupling are ω L = ω c − θ/2 and ω H = ω c + θ/2, while the corresponding cutoff wave-vectors are k L and k H 1 . The constant C TL can be calculated as where ω H and ω L are the high and low frequencies of the bandwidth. Next, we adopt the assumption that θ/ω c 1, so the ratio ω H /ω L can be expanded in the order of θ/ω c , and C tl , If we choose the characteristic impedance of the transmission line as Z tl = 50 Ω. the lowest order approximation of the parameter C tl is 0.9961. The second order term in the expansion of the coupling Hamiltonian is given by Eq. (36b). Under the rotating wave approximation (and restoring dimensions), the second order term is where the summation does not contain the terms that have the same indices. Suppose the density operator of the cavity is given by ρ(t) and the transmission line is assumed to be a vacuum bath. the master equation given by the nonlinear coupling Hamiltonian H 2 in Eq. (44) is where R is the density operator for the system (cavity) and the bath (transmission line), which can be approximated by R(t) ∼ ρ(t) ⊗ ρ B and ρ B = |vac vac| where |vac is the vacuum state of the bath (transmission line). The partial trace of the bath DOFs is where |n k1 , n k2 , ... = |n k1 |n k2 ..., |n ki is the n ki -photon Fock state of the mode k = k i . Note that because in H 2 , all the terms are aligned in the normal order, the first term in Eq. (45) is zero. We will focus on the second term of Eq. (45). After expansion of the commutation relation, the Eq. (45) iṡ Note that the Hamiltonian H 2 should be considered as the interaction picture Hamiltonian, where the transformation Next, we will work term by term in Eq. (47) to get master equation for the cavity field. We start from the term As the transmission line is assumed to be a vacuum bath, and the coupling Hamiltonian H 2 is in normal order, the partial trace will kill all the terms that contain lowering operators for the bath DOFs. The Hamiltonian terms that survive in T 3 partial trace are where the cavity nonlinear operatorB 2 is defined aŝ and the Eq. (49) is where we apply the bath state orthogonality relations and remove all the zero terms. Further, the term For the first term, which is involved in the time integral of the first line of Eq. (52) (noted as T 3,1 ), after applying the Born-Markov approximation, and define k c = ω c /v p and take Eq. (41), the term T 3,1 is Similar to the definition of the nonlienar operator of the cavity field in first order coupling, we can redefine the nonlinear operator for second order asB and then the associated rate isΓ Compared with the rate associated with the first order coupling Hamiltonian, Γ 1 [see Eq. (40)], the rate associated with this term is For a realistic setup, where we assume the transmission is 1 m long, the cavity frequency is 7.5 GHz, and the speed of microwave along the transmission line is speed of light, and the characterestic impedance of the transmission line is 50 Ω, the quantization parameterφ tl (k c ) ∼ 0.0124. So this process is much slower than the first order coupling, which is controlled by the small parameterφ 2 tl (k c ), which is equivalent to the small parameter The second line of the Eq. (52) (noted as T 3,2 ) is Note the second summation term is similar to the calculation in Eq. (54), and it is 2Γ 2,1B2 ρ(t)B † 2 . The first summation term With Born-Markov approximation, we replaceφ T k andφ T q by the central frequency mode k c = ω k /v p , and because of the fast oscillation term e i(ω k +ωq−2ωc)t , only the modes that satisfies ω k + ω q = 2ω c will have large contribution, we can approximate the integral of two modes frequencies by θ dω whereω = (ω k + ω q )/2, and θ is the coupling bandwidth. Then the integral in Eq. (59) is Note that this term is also slow compared to the first order dynamics. Similarly, to consistently compare with the first order rate in Eq. (40), we redefine the cavity operatorB as Eq. (55) and the rate associated rate Γ 2,2 as and then the ratio for the rates where if the transmission line impedance is 50 Ω, the term Z tl 2φ 2 0 ∼ 0.1560 and in the quantum optics system assumption, θ/ω c 1. So this second order coupling dynamics is also slower than the first order coupling dynamics, and is controlled by the small parameter θ/ω c .
Similarly, we can perform the same procedure for the other three terms and obtain the master equation induced by Finally, I want to note that the above derivation is valid when the length of the transmission line is large. This is consistent with the Born-Markov approximation. We assume in a coupling bandwidth θ ω c , the number of modes in this bandwidth is still much greater than the system DOFs, so the transmission line must be considered to be long, in which the integer such that we can find an other integer n θ which satisfies |n θ − n c | 1 and n θ /n c 1. In the regime where l → ∞, the rate of the quantum process given by H 2 nonlinear system-bath coupling is given by Γ 2,2 term [Eq. (61)], and is controlled by small parameter θ/ω c which does not depend on the length of the transmission line.
Similar to the discussion for the second order term in the ABOCC coupling circuit, after we apply the rotating-wave approximation, the third order Hamiltonian is given by Following the same argument as in the discussion of the second order terms, the only Hamiltonian term that contributes to the system dynamics when the bath is in vacuum state is the first term in Eq. (65a). We can define a system nonlinear operatorB In Eq. (65a), there are three terms, the first term, kφ 3 tl (k) 6 (b † k ) 3 term, will give a Lindblad term in master equation of cavity field as D[B 3 ]ρ(t) with rate Γ 3,1 . This process is further suppressed by the small parameter 1/n c [see Eq. (64)] as The second term, k,q (b † k ) 2 b † q term, will give a Lindblad term D[B 3 ]ρ(t) with rate Γ 3,2 , where θ is the coupling bandwidth. The third term k,q,p b † k b † q b † p givens the same Lindblad term with rate Γ 3,3 , In the limit where l → ∞, the third term is dominant, but is still further suppressed by θ/ω c , even compared with the second order dynamics. In the main text, we ignore the second and third order terms in the system-bath coupling Hamiltonian and only focus on the first order terms. Further, in the main text, we normalize the nonlinear cavity operatorB 1 by introducing a normalization constant N , to make the flat region of function n|B 1 |n + 1 to be unity. The normalized nonlinear cavity operator and associated decay rate (see Eq. (38) and Eq. (12) in main text) iŝ where the dimensionless parameter and the corresponding decay rate (see Eq. (40) and Eq. (11) in Main text) is When building the ABOCC laser system to achieve the narrow linewidth as we proposed in our paper, the fabrication and tuning errors may appear in the system. In this note, we consider several possible perturbations on the parameters in our ABOCC laser system, and explore how these perturbations can change the linewidth of our ABOCC laser from the perfect situation.
Notice that in real experiments, there are possible fabrication errors, e.g., the critical currents of the Josephson junctions, the linear inductance and the π-junctions inside the two ABOCCs may not match the perfect designed values; the transmon qubit frequency may mismatch the cavity frequency, and the transmon qubit impedance may also be shifted from the value we chose, which affect the parameters of our ABOCC laser. For example, the critical currents of the Josephson junctions inside the ABOCCs control the atom-cavity coupling strength g and the cavity loss rate Γ 1 , while combined with the linear inductance of the ABOCCs, they also control the shape of the ABOCC operators and change the slope of the first off-diagonal elements of the ABOCC operators. Since these experimental errors can have different effects on the laser parameters simultaneously and the authenticated simulation of the laser system on device level is beyond the scope of the current paper, we focus on the perturbation on the ABOCC laser parameters, and examine the robustness of the ultra-narrow linewidth which suppresses the standard limit against the perturbation on the laser parameters.

A. Perturbation on the coupling strengths and pump strength
In the fabrication of the ABOCC laser, the control on the critical current of the Josephson junctions inside the ABOCCs may not be accurate enough. As we mentioned above, the critical currents of the Josephson junctions control the coupling strength between the segments that the ABOCCs connect. In this subsection, we consider the perturbation on the atom-cavity coupling strength g and the cavity loss rate Γ 1 . Further, as we use Γ 1 as a frequency unit in our discussion, we effectively perturb the ratio g/Γ 1 . From Fig. 4b in main text, we notice that there is a valley that the laser can still have ultra-narrow linewidth as we tune the parameter g/Γ 1 . To reach this valley, for different coupling strength g, the incoherent pump strength Γ p needs to be fine tuned. In Fig. 5, we consider shifting the laser atom-cavity coupling g from the optimum value g c by ±0.05. We show how the laser linewidth (log(D/D ST ), shown as blue lines in Fig. 5a-c) and mean photon number ( n , shown as orange lines in Fig. 5a-c) vary as the incoherent pump strength Γ p is tuned. Notice that for the optimum coupling strength g/g c = 1, as we increase the pump strength to Γ p /Γ p,c ∼ 1.3, the incoherent pump strength is strong compared to the coupling strength, which causes the self-quenching on the single-atom laser system, as pointed out in Ref. [8] for single-atom laser. In this regime, as we increase the pump strength, the fast incoherent drive destroys the coherence between the ground state and the excited state of the atom, which reduces the photon number inside the cavity. Similar performance can be found in the case where g/g c = 1.05. Around the optimum pump strength Γ p,c for both g/g c = 1.0 and 1.05, the optimum linewidth is achieved by fine-tuning the incoherent pump strength Γ p to make the photon distribution living inside the flatten region of the ABOCC operator. This can be seen from the photon distributions shown in Fig. 5d and Fig. 5e green lines. As we decrease the atom-cavity coupling strength, the threshold for the incoherent pump induced self-quenching also decreases, which makes the achievable maximum mean photon number inside the cavity lower than the plateau in the ABOCC operators. Then the ABOCC laser cannot achieve the large suppression on the laser linewidth (see Fig. 5c and e). With a given cavity-atom coupling strength g, as long as the coupling strength is strong enough such that the maximum mean photon number can reach the designed m 0 value of the ABOCC operator, as shown in Fig. 5a and Fig. 5b, the optimum linewidth can achieved by tuning the pump strength. There are ±5% to ±10% error allowed to the pump strength (Γ p ) around the optimum pumping strength Γ p,c , where we can still get a decent suppression (> 10 times) beyond the standard limit.
Because the incoherent pump strength is not a parameter determined by the device, we can tune the pump strength experimentally, we further explore in detail on the robustness of the ultra-narrow linewidth as we perturb the cavityatom coupling strength g but leave Γ p as a tunable parameter. In Fig. 6, we focus on two ABOCC lasers, whose ABOCC operators have central photon number for the flat plateau m 0 = 1000 (see Fig. 6a) and m 0 = 500 (see Fig. 6b). As we perturb the coupling strength g, we notice that the optimum cavity-atom coupling strength g c is very close to the edge of the strong self-quenching regime. For the ABOCC laser with m 0 = 1000, the coupling strength that unable to get cavity photon distributed on the ABOCC operator plateau is g/g c = 0.99 while for m 0 = 500, the coupling strength is g/g c = 0.995 (see the red vertical lines in Fig. 6). However, as we go beyond this point, there is a large parametric space for coupling strength g to vary, which means the ABOCC laser can still get a decent suppression of the laser linewidth beyond the standard limit.

B. Perturbation on the inductance ratio in ABOCC operators
In the main text, we focus on the ABOCC laser system whose two ABOCC operators are identical and fine-tuned. Compare the ABOCC operator between the atom and the cavity (Â c,m0 ) given in Eq. (9) in main text and the operator between the cavity and the transmission lineB c,m0 given in Eq. (70) (and also Eq. (12) in main text), to make identical ABOCC operatorsÂ andB, we can fine-tune the impedance of the transmon qubit Z a , the impedance of the cavity Z c and the transmission line impedance Z tl and the inductance ratio r a-c ≡ L J:a-c /L a-c and r c-tl ≡ L J:c-tl /L c-tl . In our proposal, we manually set to ensure the two ABOCC operators to be the same. We further set Z tl = 50 Ω, which gives Z a = 47.71 Ω.
The impedance ratios r a-c and r c-tl , controls the slope of the flat region of the ABOCC operators, while the cavity impedance Z c controls the center of the flat plateau of the ABOCC operators, m 0 (see Fig. 4a in main text). To get a perfect ABOCC operator, r a-c = r c-tl = 0.418. As long as we ensure Eq. (74), m 0 for two ABOCC operators are always equal. In this subsection, we at first discuss the robustness of the ultra-narrow ABOCC laser linewidth against on one of the ABOCC operator r values. We then set the cavity-atom ABOCC operator to be the perfect (with r = 0.418), while we tune the loss ABOCC operator by tuning the r c-tl . In Fig. 7a, for each r c-tl , we optimize the laser linewidth by tuning the atom-cavity coupling strength g and the incoherent pump strength Γ p . We notice that the ultra-narrow laser linewidth does indeed require fine tuning of the ABOCC operator r c-tl . However, it can still be tolerant to small perturbation. Further, as we increase the r c-tl from the optimum value (0.418, shown by the red vertical line), the linewidth ratio D/D ST decreases slightly beyond the linewidth achieved by two fine-tuned balanced ABOCC operators. Notice that the best ratio can be achieved at a slightly larger r c-tl , i.e., the slope on the cavity loss ABOCC operator is slightly larger (see Fig. 7c orange line). The larger slope on cavity loss enhance the Boson amplification rate slightly around m 0 , which shift the steady state of the laser photon distribution. The photon distribution is shift towards the center of the plateau of the perfect ABOCC operator (the ABOCC operator for the cavity-atom coupling). This effect slightly decreases the noise added by the ABOCC operator on pump side, which enhance the linewidth performance of the ABOCC laser. On the other hand, in terms of the fabrication error, the ABOCC laser can maintain a decent suppression on the linewidth beyond the ST limit with ∼ 10% (0.42 to 0.48) on the r c-tl value.

Supplementary
Next, we consider the case that the possible fabrication errors that cause the Eq. (74) not to be held, which makes the two ABOCC operators (ABOCC between the atom and the cavity, ABOCC between the cavity and the transmission line) to have different m 0 values. Here we assume the atom-cavity ABOCC operator always have m 0 = 1000 and r a-c = r c-tl = 0.418, but manually tune the cavity-transmission line ABOCC operator m 0 . As g/Γ 1 is determined by the fabricated circuit, but Γ p can be tuned by the pump strength to drive the transmon qubit, we manually set the g/Γ 1 to be the optimum ratio for the ABOCC laser with two balanced ABOCC operators (m 0 = 1000, r a-c = r c-tl = 0.418, then g/Γ 1 ∼ 1.43). In Fig. 8, we show the ABOCC laser linewidth as we tune the cavity-transmission line ABOCC m 0 and the incoherent pump strength Γ p . We notice that as we decrease the m 0 away from the balanced value m 0 = 1000, the optimum laser linewidth slightly improves. The best linewidth ratio is achieved at Γ p /Γ 1 = 3.59, m 0 = 925 with log(D/D ST ) ∼ −5.02 (see Fig. 8b). As we further decrease the m 0 value, the optimum linewidth starts to decrease. The optimum pump strength Γ p also shifts towards Γ p /Γ 1 ∼ 4. For the ABOCC operators with m 0 ≤ 888, as we tune the pump strength away from Γ p /Γ 1 = 4, the laser system starts to be tuned into a "meta-stable" regime, where there are multiple stable states solved from the eigen-spectrum of the super-operator of the ABOCC laser systems. This regime is labeled as "Meta-stable" and filled in white in Fig. 8a and b. As we increase m 0 from 1000, the best linewidth of the ABOCC laser by tuning pump strength can still reach log(D) < −20, which gives a decent suppression beyond the standard limit. In summary, by tuning the pump strength, the ABOCC laser can maintain the narrow linewidth with decent perturbation on ABOCC operator m 0 value.

C. Perturb the atom frequency
In this subsection, we explore the perturbation on the atomic frequency to detune it from the cavity frequency. For a typical laser system, when the cavity frequency is shifted away from the atomic frequency, the laser will be degraded. Similarly, for the ABOCC laser system, as we shift the frequency of the atom away from the cavity frequency, the laser linewidth increases, while the photon distribution does not change significantly (see Fig. 9). As we increase the detuning ∆, eventually the ABOCC laser will exceed the standard limit (∆/Γ 1 ∼ 0.0178. In terms of fabrication of ABOCC laser, the cavity frequency needs to be fine-tuned to match the atomic frequency to achieve the ultra-narrow linewidth. This can be done by using a magnetic field tunable SQUID inside the transmon qubit, instead of the single Josephson junction, to make the frequency of the qubit to be tunable. By tuning the magnetic flux, the frequency of the transmon qubit can be manually tuned to match the cavity frequency with a relative high precision.

Supplementary Note 6. PROPERTIES OF THE ABOCC LOWERING OPERATOR
In this note, we consider the properties of the ABOCC lowering operator eigenstates and demonstrate that the eigenstates which correspond to the ultranarrow linewidth laser are highly squeezed. Here, we focus on the ABOCC operatorÂ c,m0 (see Eq. (9) in main text) that is associated with atom-cavity coupling. The operatorB c,m0 (see Eq. (12) in main text), which is associated with cavity-transmission line coupling, can be made to matchÂ c,m0 by tuning the properties of the output ABOCC coupler.
Our first goal is to obtain the eigenstates |ψ α ofÂ c,m0 . The eigen-equation iŝ where α is the eigenvalue. We express where |n 's are the photon number states and ψ α (n)'s are the coefficients. Since the operatorÂ c,m0 always lowers the photon number by 1, we can write it asÂ where the coefficients f (n)'s are determined by m 0 (see Fig. 4a in main text). The coefficients appearing in the eigenstates ψ α (n) obey the recurrence relation Using this recursion relation, along with the normalization condition n |ψ(n)| 2 = 1, the eigenstate coefficients can be determined if the eigenvalue α is given. Now, consider a laser tuned so that the ABOCC operator for both the cavity-atom coupling and the cavity loss term are identical with m 0 = 43 (Z c = 150 Ω). We use the cavity loss constant Γ 1 as a frequency unit and minimize the laser linewidth to obtain Γ p /Γ 1 = 3.58 and g/Γ 1 = 1.468. The corresponding photon distribution is shown in Fig. 10a (blue line). Next, we find the eigenstate of the ABOCC operator that best matches this photon distribution by tuning α in Eq. (75) and using the cost function where p(n) is the photon distribution in the laser state. In this case, we find that the optimum eigenvalue is |α| = 1.0225, where the photon distribution given by the eigenstate is shown in Fig. 10a (orange line). Without loss of generality (as we show later) we focus on the case in which the eigenvalue α is a positive real number. The corresponding Wigner distribution for the eigenstate of the ABOCC operator is plotted in Fig. 10b. For comparison, we show the Wigner distribution of a coherent state with the same phase and mean photon number in Fig. 10c. Notice that the Wigner distribution of the eigenstate is squeezed along the phase direction (and expanded along the photon number direction) as compared to a coherent state. Our next goal is to quantify the amount of squeezing. Here we consider two approaches: • Fluctuations of position and momentum [q = 1 2 (â +â † ),p = 1 2i (â −â † )]. • Fluctuations of phase and number.
The second approach suffers from the well known problem that a Hermitian phase operator of the optical field is hard to define [9][10][11]. The standard solution, which we adopt here, is to use the sine and cosine of the phase, which can be defined as the following Hermitian operatorŝ The resulting operators obey the uncertainty relation [9][10][11] For a state with a Wigner distribution that lies along the positive real axis we can estimate ∆φ ∼ ∆S/ Ĉ . Naively, the uncertainty relation between the phase and photon number of the light field implies that ∆n∆φ ∼ ∆n∆S/| Ĉ | ≥ 1/2. However, this notion does not quite hold, even for states with a Wigner distribution that lies along the positive real axis, because sinφ/ cosφ is only approximatelyφ. For example, coherent state |β saturates the uncertainty bound Eq. 81 only as β → ∞.
In Table II, we compare the fluctuations in the narrowest linewidth eigenstate ofÂ c,m0 [column labeled Eigen-NLW] and the coherent state with the same number of photons [column labeled Coh-NLW]. Comparing the momentum and position fluctuations we observe that Coh-NLW is not squeezed while Eigen-NLW state is indeed squeezed in the momentum direction. Similarly, comparing fluctuations ofŜ andn we observe that the Eigen-NLW state is squeezed as compared to the Coh-NLW state. Finally, we compute the uncertainties. We observe that the Coh-NLW state is a minimum uncertainty state that saturates the ∆p∆q uncertainty and almost saturates Eq. (81) uncertainty. On the other hand, Eigen-NLW state is not a minimum uncertainty state, which has low uncertainty but does not saturate either measure of uncertainty.
To explore how much squeezing can be obtained in the laser system, we slightly decrease the atom-cavity coupling to the value g/Γ 1 = 1.424, which corresponds to the case in which the photon number distribution of the laser state has the flattest-top (see Fig. 10d (blue line)). Again, we use the cost function in Eq. (79) to find the eigenstate of the ABOCC lowering operator that best matches the photon number distribution from the master equation and find the corresponding eigenvalue |α| ∼ 1.0059. The photon number distribution of this eigenstate is shown in Fig. 10d (orange line). The corresponding Wigner function is plotted in Fig. 10e and the coherent state with the same photon number is shown in Fig. 10f. We observe that this eigenstate, which we label "Eigen-FT" for short, is even stronger squeezed than the Eigen-NLW state. This notion is confirmed in Table II.
Having confirmed that the eigenstates of the ABOCC lowering operator are indeed squeezed, we investigate how the squeezing depends on the number of photons in the laser cavity. In Fig. 11, we plot the mean photon number n and the noise inp direction of the eigenstates of the ABOCC operator with two different m 0 values (which is controlled by the cavity effective impedance Z c ), as a function of the eigenvalue |α|. In Fig. 11a, we focus on the case Supplementary Table II. The properties of the eigenstates of the ABOCC operators and the corresponding coherent states with the same mean photon number. State "Eigen-NLW" is for eigenstate of the ABOCC operator that best-matches the narrowest linewidth laser state (see Fig. 10a). State "Eigen-FT" is for the eigenstate of the ABOCC opreator that best-matches the flattest-top laser state (see Fig. 10d). State "Coh-NLW" and "Coh-FT" are the corresponding coherent state with the same mean photon number.   In Fig. 12 and Fig. 13, we show a few more states with different ABOCC operators. The target photon number m 0 of the ABOCC operators can be controlled by tuning the cavity effective impedance Z c (also see Fig. 4a in main text). In Fig. 12 and Fig. 13, we consider three ABOCC operators: the ABOCC operators with cavity impedance Z c = 100 Ω (m 0 = 64, Fig. 12a and d, Fig. 13a and d), with cavity impedance Z c = 50 Ω (m 0 = 128, Fig. 12b and e, Fig. 13b and e) and with cavity impedance Z c = 40 Ω (m 0 = 159, Fig. 12c and f, Fig. 13c and f). In each case we find the eigenstate of the corresponding ABOCC operator that best matches the photon number distribution of the narrowest linewidth laser state. The photon distribution of the laser states and the ABOCC eigenstates are shown in Fig. 12a-c and the Wigner distribution of the corresponding ABOCC eigenstates are shown in Fig. 12d-f. As the mean photon number increases with decreasing cavity impedance, the Wigner distribution width drops along thep direction, while the width along theq direction increases. This indicates that photon number noise increases while the phase noise drops. Similarly, we also calculate the ABOCC eigenstates that best match the flat-top laser states in Fig. 13. We again focus on these three cavity impedance Z c = 100 Ω, 50 Ω and 40 Ω. In this case, we observe the same trend: the Wigner distributions of the ABOCC lowering operator eigenstates become wider in theq direction and narrower along thep direction as we increase m 0 .
To systematically investigate the noise and the squeezing along the phase direction of the ABOCC eigenstates, we sweep the eigenstates of the ABOCC operators with different m 0 , and calculate the number-phase uncertainty (∆n∆S/| Ĉ |) and the momentum-position uncertainty relation ∆p∆q in Fig. 14a and Fig. 14b. For each point, we use the ABOCC operator with the desired m 0 value to construct the ABOCC laser system. We solve the laser state that gives the best linewidth performance (labeled as "NL") and the flat-top photon distribution (labeled as "FT"). Then we solve the ABOCC operator eigenstates and fit the eigenvalues by minimizing the cost function in Eq. (79). With the best matched eigenstates of the ABOCC operator, we calculate the phase noise and the photon number noise, then the photon-phase uncertainty relation (in Fig. 14a) and the momentum-position uncertainty in (in   14b). Notice that the dashed red lines in both plots are the quantum limit (0.5). The eigenstates of the ABOCC lowering operators under investigation are not minimum-uncertainty states. However, the uncertainty relation of both NL case and the FT cases change only slightly within the range shown in Fig. 14a and b. Next, we consider the photon number noise for both the NL and the FT cases. Notice that the photon noise (∆n) increases as we increase m 0 (see Fig. 14c). We fit the scaling parameter and the dependence on the photon number is ∆n ∼ m 0.761 0 (see red dashed line in Fig. 14c). This scaling performance can be analyzed using a hypothetical ABOCC operator with the first diagonal elements in Fock basis satisfies f (n) = 1 + (n − m 0 ) 3 and minimize the laser performance by also tuning the small parameter . 2 Further, we can estimate the scaling parameter and get the relation ∆n ∼ n 0.820 . As the photon number-phase uncertainty relation only change slightly as the m 0 value is changed, the phase noise ∆S/| C | ∼ m −0.761 0 ∼ n −0.820 . Notice that for a coherent state, the photon number noise is ∆n = n 0.5 , so the corresponding phase noise should be ∼ n −0.5 . The eigenstates of our proposed ABOCC lowering operator will have a larger photon number noise, but squeezed along the phase direction. As we increase the m 0 value, the phase noise is further squeezed.
As a final sanity check, we verify that the squeezing direction of the eigenstates of the ABOCC lowering operator is always along the phase direction (rather than thep direction). To verify this, we construct an eigenstate with an eigenvalue α → αe iθ . We expect that the resulting Wigner distribution is rotated around the origin ψ(n) → ψ(n)e inθ . In Fig. 15, we plot the Wigner distribution of the eigenstates of ABOCC lowering operator with Z c = 100 Ω (m 0 = 64) with θ = π/4. We observe that the Wigner distribution indeed rotated by π/4, as has the direction of squeezing, as compared to Fig. 12d. ] of the ABOCC operator eigenstates that best-matched to the narrowest-linewidth laser states (labeled as "NL") and the flat-top photon distribution (labeled as "FT"), as we increase the m0 of the ABOCC operator. The dashed lines show the minimum uncertainty relation. In (c), we calculate the photon number noise in the corresponding eigenstates of the ABOCC lowering operators as we increase the m0 value. The dashed line shows the linear fit in the log-log scale, which gives a scaling parameter 0.761.