Deterministic three-photon down-conversion by a passive ultrastrong cavity-QED system

In ultra- and deep-strong cavity quantum electrodynamics (QED) systems, many intriguing phenomena that do not conserve the excitation number are expected to occur. In this study, we theoretically analyze the optical response of an ultrastrong cavity-QED system in which an atom is coupled to the fundamental and third harmonic modes of a cavity, and report the possibility of deterministic three-photon down-conversion of itinerant photons upon reflection at the cavity. In the conventional parametric down-conversion, a strong input field is needed because of the smallness of the transition matrix elements of the higher order processes. However, if we use an atom-cavity system in an unprecedentedly strong-coupling region, even a weak field in the linear-response regime is sufficient to cause this rare event involving the fourth order transitions.


I. INTRODUCTION
The history of cavity quantum electrodynamics (QED) parallels with the enhancement of the atom-cavity coupling g. From the observations of the suppressed or enhanced atomic decay rate in the weak coupling regime (g < κ, γ, where κ and γ are the loss rates of the cavity and the atom, respectively) [1][2][3][4][5], more than a decade was required to reach the strong coupling regime (g > κ, γ), where the Rabi splitting or oscillation becomes observable [6][7][8][9][10][11]. Recently, the ultrastrong coupling regime (g ω a /10, ω c /10, where ω a and ω c are the resonance frequencies of the atom and cavity, respectively) and even the deep-strong coupling regime (g ω a , ω c ) have been realized in various physical platforms such as polaritons, superconducting qubits, and molecules [12][13][14][15][16][17][18]. In the ultra-and deep-strong coupling regimes, various novel phenomena originating in the counter-rotating terms of the atom-cavity coupling are expected to become observable.
One of such phenomena is the deterministic nonlinear optics [19][20][21]. The generation efficiency of entangled photons, which is typically of the order of 10 −6 for spontaneous parametric down-conversion with a bulk nonlinear crystal [22][23][24], might be drastically improved by this scheme. However, the proposed deterministic nonlinear-optical processes are for intracavity photons. In order to apply this scheme for itinerant photons, deterministic capturing of propagating photons into a cavity is indispensable. This is in principle possible but requires a precise dynamic control of the external cavity loss rate in accordance with the incoming photon shape [25,26]. In contrast, deterministic downconversion of itinerant photons is possible and has been demonstrated in a passive waveguide QED setup [27][28][29].  1: Schematic of the considered setup. A qubit (resonance frequency ωq) interacts with the first (ω1) and third (ω3) cavity modes with coupling constants g1 and g3, respectively. Dissipation channels are as follows: qubit radiative decay (rate γ), external and internal losses of the first cavity mode (κ1e, κ1i) and those for the third one (κ3e, κ3i). A weak monochromatic field (frequency ωin ∼ ω3) is input through the waveguide.
The physical origin of the deterministic conversion is the destructive interference between incoming field and radiation from the atom. Accordingly, the efficiency is sensitive to the external loss rates of the relevant transitions of the atom.
In this study, we investigate an ultrastrong cavity QED system in which an atom is placed at the center of the cavity and is thus coupled to the fundamental and third-harmonic modes of the cavity (Fig. 1) [30]. We show that, if the atom-cavity and cavity-waveguide couplings are adequately chosen, an input photon (resonant to the third mode) is down-converted nearly deterministically to three daughter photons (resonant to the fundamental mode) upon reflection at the cavity . The drastic enhancement of the conversion efficiency in comparison with the prior demonstrations of triplet-photon generation [31][32][33][34][35][36][37][38] originates in the waveguide QED effect, in other words, the engineered dissipation rates of the optical system. Such deterministic conversion is possible even in the conventional strong-coupling cavity QED. However, considering the required intrinsic loss rates of the cavity modes, this phenomenon is characteristic to the ultrastrong cavity QED.
The rest of this paper is organized as follows. We present the theoretical model of the atom-cavity-waveguide coupled system in Sec. II. We observe the internal dynamics of the atom-cavity system and evaluate the effective coupling between the two levels relevant to three-photon down-conversion in Sec. III. We develop the input-output formalism applicable to the ultrastrong coupling regime in Sec. IV, and apply to the investigated setup in Sec. V. We numerically show that the deterministic three-photon down-conversion is possible in Sec. VI and clarify the required conditions. Section VII is devoted to a summary.

II. SYSTEM
In this study, we investigate an optical response of a cavity QED system schematically illustrated in Fig. 1. A two-level system (qubit) is placed at the center of a cavity and interacts with the first-and third-harmonic cavity modes. These cavity modes are coupled to an external waveguide field through the right mirror (capacitor, in circuit QED implementation), and a monochromatic field close to the resonance of the third cavity mode is applied through the waveguide.
We consider five dissipation channels of this system: the external/internal decay of the first/third cavity mode and the longitudinal decay of the qubit. We label these dissipation channels as 1e, 1i, 3e, 3i and q, respectively, and denote their rates as κ 1e , κ 1i , κ 3e , κ 3i , and γ, respectively. As we observe later (Fig. 6), the investigated phenomenon is robust against the qubit dissipation, since the qubit excited state is used only virtually. We therefore neglect the qubit pure dephasing for simplicity, which plays essentially the same role as the longitudinal decay in the present phenomenon. The Hamiltonians describing the 1e and 3e channels are given bŷ whereĉ k (d k ) is the annihilation operator of a waveguide mode with wavenumber k coupled to the first (third) cavity mode and ω k is its frequency. The commutators forĉ k andd k are given by . Although these modes can be treated as a single continuum in principle, we may safely treat them as independent continua because of large separation between the relevant frequencies. The dispersion relation in the waveguide is linear, ω k = k, where the velocity of waveguide photons is set to unity for simplicity. A real quantity ξ (je) k (j = 1, 3) represents the cavity-waveguide coupling. By naively applying the Fermi golden rule, the coupling and the decay rate are related by The other loss channels are modeled similarly. For example, the longitudinal relaxation of the qubit is modeled bŷ  The relation to the qubit decay rate is γ = 2πξ In this section, we investigate the properties of the qubit-cavity system, neglecting dissipation for the moment. We denote the state vector of the system by |qmn , where q(= g, e) specifies the qubit state and m, n (= 0, 1, · · · ) specify the photon numbers of the first and third cavity modes, respectively. In particular, we focus on the coupling between |g01 and |g30 , which is essential for the three-photon down conversion. Figure 2 shows the transition paths between |g01 and |g30 . We observe that |g01 and |g30 are coupled through the fourth order process in the qubit-cavity coupling g, and that transitions non-conserving the excitation number (dotted arrows in Fig. 2) are indispensable for their coupling.
For a strong coupling between these two states, degeneracy of these two states is required. The eigenenergies of |g01 and |g30 measured from |g00 , which we denote by ε g01 and ε g30 , are renormalized by the dispersive qubitcavity coupling and depend on g 1, 3 and ω q . In Fig. 3(a), by numerically diagonalizingĤ s (taking into account up to the 6 (2) photon states forâ 1 (â 3 ) mode), ε g01 and ε g30 are plotted as functions of the qubit frequency ω q . We observe an anticrossing between them, and from this plot we can identify the following quantities. (i) The optimal qubit frequency ω opt q , at which the level spacing is minimized. (ii) The effective coupling g eff between |g01 and |g30 , which is half of the level spacing at ω opt q . (iii) The optimal input photon frequency ω opt in , which is the average of the two transition frequencies at ω opt q . In Figs. 3(b) and (c), assuming 3ω 1 = ω 3 = 2π × 9 GHz and g 1 = g 3 (= g), we plot ω opt q , ω opt in and g eff as functions of g. We can confirm in Fig. 3(c) that g eff is proportional to g 4 , which is consistent with the fact that |g01 and |g30 are coupled through the fourth order process. We also observe in Fig. 3(b) that, in the weak-coupling limit of g → 0, ω opt in → ω 3 = 3ω 1 = 2π × 9 GHz and ω opt q → (3ω 2 3 − ω 2 1 )/2 = 2π × 10.82 GHz. This optimal qubit frequency in the weak-coupling limit is derived as follows. Within the second-order perturbation in g 1 and g 3 , ε g01 and ε g30 , respectively. The degeneracy condition, ε g01 = ε g30 , reduces to ω q = (3ω 2 3 − ω 2 1 )/2 for the present case of ω 3 = 3ω 1 and g 1 = g 3 . In Fig. 3(d), coherent time evolution of the system starting from |g01 is shown. We observe the vacuum Rabi oscillation between |g01 and |g30 . The oscillation is not pure sinusoidal, since other states than |g01 and |g30 (such as |e00 ) are involved in forming the eigenstates. The Rabi oscillation period measured in Fig. 3(d) is T = 2.247 µs. This is compatible with g eff = 2π × 222.5 kHz evaluated from Fig. 3(a), because g eff T /π ≈ 1.

IV. INPUT-OUTPUT FORMALISM FOR ULTRASTRONG CAVITY QED
In the following part of this paper, we analyze the response of the qubit-cavity system to a waveguide field. For this purpose, the input-output formalism is useful. However, although highly useful in the weak-and (usual) strongcoupling regimes of cavity QED, the conventional formalism based on the white reservoir approximation has several difficulties in treating the ultrastrong cavity QED [39]. The input-output formalism applicable to the ultrastrong cavity QED has been discussed first in linear systems [40] and later extended to nonlinear systems [41,42]. However, assuming weak dissipation, the counter-rotating terms in the system-environment coupling have been neglected in prior works. In this section, in order to extend the formalism applicable to highly dissipative cases, we develop an input-output formalism incorporating the counter-rotating terms in the system-environment coupling.

A. Heisenberg equations
In this section, for simplicity, we consider a case in which the cavity QED system is coupled only to the 1e dissipation channel. Namely, the overall Hamiltonian is given byĤ =Ĥ s +Ĥ 1e . We also omit the superscript (1e) throughout this section (for example, ξ The Heisenberg equation of an arbitrary system operatorŜ is given by . Using Eq. (6), this becomes a delay differential equation, whereΓ(t) is a noise operator defined byΓ In Eq. (7) and hereafter, we omit explicit time-dependence for the operators at time t in the Heisenberg picture. Note that the noise operatorΓ(t) is not in the Heisenberg picture.
In order to convert a delay differential equation (7) into an instantaneous one, we employ a free-evolution approximation for the time evolution of the system operator during the delay time τ [43][44][45]. We denote the eigenstates and eigenenergies ofĤ s by |j and ε j (j = 0, 1, · · · ) in the energy-increasing order, and define the transition operator bŷ s ij = |i j|. We then approximateX(t − τ ) aŝ where x ij = i|X|j and ε ji = ε j − ε i . Then, Eq. (7) is rewritten as and h ji = t 0 dτ ∞ 0 dk ξ 2 k e i(εji −k)τ . Strictly speaking, this quantity depends on t. However, since ∞ 0 dk ξ 2 k e −ikτ is nonzero only for short τ (∼ inversed bandwidth of ξ 2 k ), we can safely treat h ji as a t-independent quantity, where θ is the Heaviside step function and P represents the principal value. This is a self energy correction to the transition frequency: the real part corresponds to half of the decay rate and the imaginary part corresponds the Lamb shift [46]. Note that the real part vanishes for negative ε ji , reflecting the prohibited decay from a lower level to an upper level.

B. Waveguide-field operator
For a semi-infinite waveguide (Fig. 1), the waveguide eigenmodes are the standing wave extending in the r > 0 region. Assuming an open boundary condition at r = 0, the eigenmode function f k with wavenumber k is given by f k (r) = 2/π cos(kr)θ(r), which is orthonormalized as . Accordingly, the waveguide mode operatorĉ r in the real-space representation is defined in the r > 0 region byĉ r (t) = 2/π ∞ 0 dk cos(kr)ĉ k (t), which satisfies the commutation relation [ĉ r ,ĉ † r ′ ] = δ(r − r ′ ). The incoming field propagates in the r > 0 region into the negative direction (k < 0).
However, it is convenient to treat the incoming field as if it propagates in the r < 0 region into the positive direction (k > 0). We therefore introduce the real-space representation by where r runs over the full one-dimensional space (−∞ < r < ∞) [47]. From Eq. (6), we have This equation represents the waveguide-field operator in terms of the input-field operator and the system operator, and enables, for example, evaluation of the output field amplitude and flux.

C. Input-output relation
We can further simplify Eq. (14) under some approximations. Introducing τ ′ = τ − r and employing the freeevolution approximation, (14) is rewritten as Note that the integrand in the right-hand side of Eq. (16) is not singular at k = ε. We can approximately evaluate f (ε, r, t) as follows. Since the main contribution of this integral comes from the k ≈ ε region, we set ξ k ≈ ξ ε θ(ε) and remove the lower limit of k integral as In Appendix A, we observe fairly good agreement between f (ε, r, t) and f app (ε, r, t), assuming a concrete form of ξ k [Eq. (28)]. Thus we have Note that the summation over i and j is conditioned by i < j in Eq. (18), which is due to the following reasons. (i) For i > j, ε ji is negative and accordingly θ(ε ji ) = 0. (ii) For i = j, x ij = 0 due to the parity selection rule. Defining the input and output field operators byĉ in (t) = c −0 (t) = c −t (0) andĉ out (t) = c +0 (t), Eq. (18) is rewritten into a more familiar form,ĉ V. OPTICAL RESPONSE THEORY

A. Initial state vector
In this study, instead of a single photon pulse, we apply a weak classical monochromatic field close to the resonance of the third cavity mode [E in (t) = E in e −iωint with ω in ∼ ω 3 ] to the cavity (Fig. 1). Assuming that, at the initial moment (t = 0), the overall system is in the vacuum state except the applied field, the initial state vector is written as where |vac is the overall vacuum state. Note that this is an eigenstate of the noise operators: ωin E in (t)|ψ i andΓ (j) (t)|ψ i = 0 for j = 1e, 1i, 3i, and q.

B. Density matrix elements
The Heisenberg equation for a system operator is given by Eq. (10) with the dissipators and the noise operators corresponding to the five decay channels (1e, 1i, 3e, 3i and q). The equation of motion for s ij (t) = ψ i |ŝ ij (t)|ψ i , which is identical to the density matrix element ρ ji (t) = j|ρ(t)|i in the Schrödinger picture, is then given by where the coefficients η (1,2,3) ijmn are given by where h ij . In this study, we apply a continuous field and observe the stationary response of the system. Therefore, we numerically determine the stationary solution of these simultaneous equations by perturbation with respect to E in (t). Further details on analysis are presented in Appendix B.

C. Photon flux
Since we treat a stationary input/output field, we quantify the amount of photons by the photon flux, namely, the rate of incoming/outgoing photons per unit time. The input photon flux is evaluated by F in = ψ i |d † in (t)d in (t)|ψ i . From Eq. (20), this quantity reduces to In the output port, the fluxes of down-converted and unconverted photons are respectively evaluated by F 1 out = ψ i |ĉ † out (t)ĉ out (t)|ψ i and F 3 out = ψ i |d † out (t)d out (t)|ψ i . From Eq. (19) and its counterpart ford out , these quantities are given by

VI. NUMERICAL RESULTS
In this section, we present the numerical results on the optical response, fixing the bare cavity frequencies at 3ω 1 = ω 3 = 2π × 9 GHz. For reduction of parameters, we restrict ourselves to the case of g 1 = g 3 , κ 1e = κ 3e and κ 1i = κ 3i and denote them by g, κ e and κ i , respectively. Furthermore, regarding the coupling for the 1e decay channel for example, we assume the following form, where k x is the cutoff wavenumber. Note that this coupling satisfies Eq. (4). We fix k x at 2π × 20 GHz and confirmed that numerical results are mostly insensitive to k x . The other system-environment couplings are defined similarly and with the same cutoff wavenumber. From Eq. (12), h (1e) ji is analytically given by Regarding the input field power, we assume the weak-field limit in Secs. VI A-VI C, and discuss the input power dependence in Sec. VI D.

A. Optimal condition for κe
First, we search the optimal value of κ e assuming no intrinsic losses (κ i = γ = 0). Figures 4(a) and (b) show the dependence of the down-converted flux F 1 out on κ e and ω in , fixing g and ω q . It is observed that F 1 out has two peaks for small κ e . This is due to the Rabi splitting of |g30 and |g01 , and the frequency difference of the two peaks agrees with 2g eff in Fig. 3(c). F 1 out is maximized for a larger κ e , at which the two peaks become spectrally indistinguishable. Therefore, the optimal condition for the external loss rate of the cavity is given by κ opt e ∼ g eff . We confirm in Appendix C that this condition is identical to the impedance-matching condition of a linear optical system composed of oscillators and waveguides. Actually, we can confirm in Figs. 4(a) and (b) that κ opt e is 2π × 255 kHz (35.4 MHz) for g = 2π×0.3 GHz (1.0 GHz). This is almost identical to g eff = 2π×223 kHz (29.2 MHz) in Fig. 3(c). Figures 4(c) and (d) are the cross section of Figs. 4(a) and (b) at κ opt e . It is observed that the deterministic downconversion (F 1 out ≈ 3F in and F 3 out ≈ 0) is attained regardless of the value of g, when the input photon frequency ω in is optimally chosen. Furthermore, reflecting the absence of intrinsic loss channels, we can also confirm the energy conservation, F 1 out /3 + F 3 out ≈ F in , for any input photon frequency. However, by carefully examining the numerical results, this conservation law is slightly broken at the order of 10 −5 [10 −3 ] in Fig. 4(c) [Fig. 4(d)]. We attribute the main reason for this slight discrepancy to the free-evolution approximation [Eq. (9)], whose validity is gradually lost for larger dissipation rates.

B. Qubit detuning
Here, assuming again the absence of intrinsic losses, we observe the effects of the qubit detuning from its optimal value. Figure 5 shows the dependence of the down-converted flux F 1 out on ω q and ω in , fixing κ e at its optimal value [255 kHz in (a) and 35.4 MHz in (b)]. It is observed that, as the qubit-cavity coupling g increases, the deterministic down-conversion becomes more robust against the qubit detuning. This is because of the increase of the optimal qubit linewidth for larger g. The allowed qubit detuning (the full width in ω q at the half maximum of the cross sectional plot at the optimal ω in ) is about 20 MHz (240 MHz) for g = 2π × 0.3 GHz (1.0 GHz). Here, we investigate the effects of intrinsic losses of the qubit and cavity. Figure 6 shows the dependence of the down-converted flux F 1 out on γ and κ i , fixing the other parameters at their optimal values. When g = 2π×0.3 GHz, the deterministic down-conversion is highly vulnerable to the intrinsic losses. The condition for achieving 50% conversion (F 1 out > 1.5) is κ i 2π×92.9 kHz (intrinsic quality factor Q i 3.23 × 10 4 for the first cavity mode) and γ 2π×5.55 MHz (lifetime T 1 28.7 ns). These conditions are drastically relaxed for g = 2π×1.0 GHz: κ i 2π×12.7 MHz (Q i 236) and γ 2π×70.5 MHz (T 1 2.26 ns). We observe that the condition for the cavity is tighter than that for the qubit. This is because, in the present phenomenon, the qubit excited state is used only virtually to realize the effective coupling between |g01 and |g30 states.

D. Dependence on input photon rate
In the previous subsections, we discussed the down-conversion efficiency assuming a low input photon rate, in other words, the linear-response limit. Here, we observe the conversion efficiency for a higher input photon rate. In Fig. 7, we plot the dependence of the conversion efficiency on the input photon rate for various detuning of the input field. We observe that the efficiency decreases gradually for higher input photon rate. This is due to saturation of the atom-cavity system, which originates from the nonlinearity of the qubit. The star symbols in Fig. 7 represent the onset of saturation, which is given by where ∆ω is the detuning of the input field frequency from its optimal value. This is derived as follows. When one applies a monochromatic field E(t) = E in e −iωint to an empty one-sided cavity with an external decay rate κ e , the mean intracavity photon number n is proportional to the drive photon rate F in = |E in | 2 and is given by n = κ e F in /|κ e /2 + i∆ω| 2 . In the present system, the saturation effect due to nonlinearity would appear when the cavity is populated substantially. If we set this criterion at n ∼ 0.1 for example, the onset of saturation is estimated by Eq. (30). This explains the fact that the onset of saturation occurs at a higher input photon rate for a larger detuning.

VII. SUMMARY
We theoretically proved the possibility of the deterministic three-photon down-conversion of itinerant photons using a passive ultrastrong cavity QED system, in which an atom is coupled to the fundamental and third-harmonic cavity modes. For this purpose, we developed an input-output formalism applicable to highly dissipative cavity QED systems. The conditions for the deterministic conversion are as follows: (i) the frequencies of the qubit and the cavity modes are adequately chosen so that the two relevant levels (|g30 and |g01 ) are coupled effectively, and (ii) the  Fig. 4(a). (b) Results under the optimal condition in Fig. 4(b). cavity loss rates are adequately chosen so that they are comparable to the effective coupling. Such down-conversion is characteristic to the ultrastrong coupling regime of cavity QED, considering the upper limit of the intrinsic loss rates of the cavity. In this Appendix, we present the method to determine the stationary solution of Eq. (21) perturbatively. As the stationary solution, we employ the following form, where s The stationary solution is readily obtained by the replacement of d/dt → −iω in in Eqs. (C5) and (C6). The transmission coefficient, T = b out,2 / b in,1 , is then given by b out,2 b in,1 = i √ κ 1 κ 2 g eff κ 1 κ 2 /4 + g 2 eff . (C9) The impedance-matching condition, |T | = 1, reduces to √ κ 1 κ 2 = 2g eff . This is in agreement with the optimal condition of the external cavity decay rate, κ opt e ∼ g eff , derived in Sec. VI A.